M Sanchez
email: msanchez@toulouse.inra.fr

D Allouche
email: allouche@toulouse.inra.fr

S De Givry
email: degivry@toulouse.inra.fr

T Schiex
email: tschiex@toulouse.inra.fr

Poupées russes et décomposition arborescente

Les problèmes d'optimisation dans les modèles graphiques ont été étudiés dans différents formalismes de l'intelligence artificielle tels que les CSP pondérés (weighted CSP ou WCSP), le maximum de satisfiabilité (MaxSAT) ou dans les réseaux probabilistes (réseaux Bayésiens et champs markoviens). En identifiant des sous-problèmes conditionnellement indépendants,qui sont résolus de façon indépendante et dont l'optimum est mémorisé, il est possible de rendre les algorithmes de type Séparation-Évaluation (Branch and Bound) asymptotiquement plus efficaces. Mais la localité des bornes induites par la décomposition affaiblit les effets réels de ce résultat asymptotique car de nombreux sous-problèmes, ne participant pas à la construction d'une solution optimale, sont résolus à l'optimum inutilement.

En s'inspirant de l'algorithme des poupées russes (RDS pour Russian Doll Search), une approche possible pour surmonter cette faiblesse consiste à (récursivement) résoudre une relaxation de chacun des sous-problèmes pour obtenir des bornes renforcées. L'algorithme ainsi défini généralise à la fois l'algorithme RDS mais aussi les algorithmes de types Branch and Bound exploitant une décomposition arborescente tels que BTD ou AND-OR Branch and Bound. Nous étudions son efficacité sur différents problèmes et fermons une instance très dure de problème d'affectation de fréquences ouverte depuis plus de 10 ans.

Introduction

Le traitement de modèles graphiques est un problème central de l'intelligence artificielle. L'optimisation d'une combinaison de fonctions de coût local est en particulier étudiée dans les réseaux de contraintes valués [START_REF] Schiex | Valued constraint satisfaction problems : hard and easy problems[END_REF] . Elle permets de capturer des problèmes tels que le problème MaxSAT pondéré, CSP pondérés, recherche d'une explication de probabilité maximum (MPE pour Maximum Probability Explanation) dans les réseaux Bayésiens et les champs markoviens. Ses applications directes sont nombreuses notamment en affectation de ressources [START_REF] Verfaillie | Russian Doll Search for Solving Constraint Optimization Problems[END_REF][START_REF] Cabon | Radio link frequency assignment[END_REF] et en bioinformatique [START_REF] Sanchez | Mendelian error detection in complex pedigrees using weighted constraint satisfaction techniques[END_REF].

Ces dernières années, afin de résoudre des problèmes de satisfaction, d'optimisation ou de comptage, différents algorithmes de recherche arborescente exploitant le couplage d'une décomposition arborescente du problème avec une propagation des informations dures du réseau ont été proposés. Cette classe d'algorithme couvre l'algorithme de Recursive Conditioning [START_REF] Darwiche | Recursive Conditioning[END_REF], l'algorithme Backtrack bounded by Tree Decomposition (BTD) [START_REF] Terrioux | Hybrid backtracking bounded by tree-decomposition of constraint networks[END_REF] et l'algorithme AND-OR graph search [START_REF] Marinescu | Memory intensive branch-and-bound search for graphical models[END_REF], tous héritiers de l'algorithme de Pseudo-Tree Search [START_REF] Freuder | Taking advantage of stable sets of variables in constraint satisfaction problems[END_REF]. En comparaison avec des algorithmes de recherche arborescente traditionnels, ces algorithmes offrent des complexités asymptotiques améliorées, seulement exponentielles dans la largeur d'arbre du graphe du problème. Ce gain s'obtient toutefois au prix d'une restriction dans l'ordre d'affectation des variables (voir [START_REF] Jégou | Dynamic management of heuristics for solving structured CSPs[END_REF]).

Dans le contexte de l'optimisation (sans perte de généralité, nous considérons des problèmes de minimisation), l'exploitation d'une décomposition du problème a pour effet secondaire la perte de bornes globales. Dans les algorithmes Branch and Bound en profondeur d'abord traditionnels (DFBB pour Depth First Branch and , Bound), l'élagage se produit dès que le coût de la meilleure solution connue (le majorant) rencontre un minorant spécifique, habituellement calculé à partir d'une relaxation du problème global. Lorsqu'une décomposition arborescente est exploitée et qu'une affectation des variables connectant un sous-problème au reste du problème rend ce sous-problème indépendant, une solution optimale du sous-problème (conditionné par cette affectation) est recherchée. Ce calcul de solution optimale est coûteux, et souvent inutile car la solution optimale du sous-problème obtenue ne fait pas nécessairement partie de la solution optimale globale (du fait du reste du problème). Pour limiter ce comportement de "trashing" désagréable, il est possible d'informer l'algorithme en charge de la recherche d'une solution optimale qu'il doit absolument produire une solution suffisamment bonne pour pouvoir s'intégrer dans une solution globale ou sinon s'arrêter dès qu'il prouve que c'est infaisable. Ce majorant initial est obtenu en soustrayant du majorant global une combinaison des minorants disponibles sur le reste du problème.

Dans cet article, nous avons tenté d'améliorer l'algorithme décrit dans [START_REF] De Givry | Exploiting Tree Decomposition and Soft Local Consistency in Weighted CSP[END_REF] qui introduit des minorants obtenus par un filtrage cohérence locale pondérée [START_REF] Cooper | Arc consistency for soft constraints[END_REF] dans l'algorithme BTD. Dans l'objectif d'une meilleure prise en compte de la contribution du problème global dans la définition des majorants locaux, nous remplaçons le moteur de type "Branch and Bound" utilisé dans BTD par une algorithme de type "Poupées Russes" [START_REF] Verfaillie | Russian Doll Search for Solving Constraint Optimization Problems[END_REF]. L'algorithme RDS est étendu afin de pouvoir exploiter une décomposition arborescente et un niveau de cohérence local plus élevé. Dans l'algorithme RDS-BTD ainsi défini, des relaxations des sous-problèmes conditionnellement indépendants identifiés par la décomposition arborescente sont résolues de façon récursive. Les coûts des solutions optimales de ces relaxations définissent de puissants minorants et fournissent donc des majorants locaux améliorés, permettant ainsi de réduire de façon drastique le comportement de "trashing" observé.

Dans la suite de l'article, nous introduisons d'abord le cadre des CSP pondérés, l'algorithme BTD et l'algorithme RDS. Nous présentons ensuite notre algorithme générique qui peut se spécialiser en BTD ou en RDS. Enfin, l'algorithme RDS-BTD est évalué sur différentes instances de problèmes réels issues du domaine de l'affectation de fréquences et de la sélection de "Tag SNP" (en bioinformatique).

Weighted Constraint Satisfaction Problem

Un CSP pondéré (Weighted CSP, WCSP) est un quadruplet (X, D, W, m). X et D sont des ensembles de n variables et domaines, comme dans les CSP classiques. Le domaine de la variable i est noté D i , avec une taille de domaine maximum notée d. Étant donné un sous-ensemble de variables S ⊆ X, on note ℓ(S) l'ensemble des n-uplets (tuples) définis sur S. W est un ensemble de fonctions de coût. Chaque fonction de coût (ou contrainte molle) w S de W est définie sur un ensemble de variables S appelé sa portée. Cette dernière est supposée différent pour chaque fonction de coût. Une fonction de coût w S assigne un coût à chaque affectation des variables de S i.e. : Pour les fonctions de coût binaires et unaires, nous utilisons des notations simplifiées : une fonction de coût binaire entre les variables i et j est notée w ij . Une fonction de coût unaire sur la variable i est notée w i . Nous faisons l'hypothèse qu'il existe pour chaque variable, une fonction de coût unaire w i ainsi qu'une fonction de coût d'arité nulle notée w ∅ correspondant à un coût constant payé par toutes affectations.

Le coût d'une affectation complète t ∈ ℓ(X) dans un problème P = (X, D, W, m) est égal à Val P (t) = w S ∈W w S (t[S]) où t[S] représente la projection habituelle d'un tuple sur l'ensemble de variables S. La minimisation de Val P (t) défini un problème d'optimisation avec un problème de décision associé qui est NP-complet.

L'établissement d'une cohérence locale donnée sur un problème P consiste à transformer P = (X, D, W, m) en un problème P ′ = (X, D, W ′ , m) qui est équivalent à P (Val P = Val P ′ ) et qui satisfait la propriété de cohérence locale considérée. Ce processus peut augmenter w ∅ et fournir ainsi un minorant amélioré du coût optimum. Il s'appuie sur l'application de transformations préservant l'équivalence (EPT) qui déplace les coûts entre des portées différentes [START_REF] Schiex | Arc consistency for soft constraints[END_REF][START_REF] Cooper | Arc consistency for soft constraints[END_REF][START_REF] De Givry | Existential arc consistency : Getting closer to full arc consistency in weighted CSPs[END_REF][START_REF] Cooper | Optimal soft arc consistency[END_REF]. Un algorithme de recherche peut exploiter la propriété précédente pour autant qu'un ordre des variables adapté soit utilisé : les variables d'un cluster C e quelconque doivent être instanciées avant les variables restantes dans ses clusters fils. Dans ce cas, pour tout cluster C f ∈ Fils(C e ), une fois le séparateur S f affecté, le sous-problème P f conditionné par l'affectation courante A f de S f (notée P f /A f ) peut être résolue à l'optimalité indépendemment du reste du problème. Le coût optimal de P f /A f peut alors être enregistré avec pour conséquence qu'il ne sera jamais plus résolu pour cette affectation du séparateur S f . Ainsi, des algorithmes tels que BTD ou AND-OR graph search sont capables d'avoir une complexité temporelle exponentielle seulement dans la taille du plus grand cluster.

Dans [START_REF] De Givry | Exploiting Tree Decomposition and Soft Local Consistency in Weighted CSP[END_REF], l'algorithme BTD exploite de minorants produits par le maintien d'une propriété de cohérence locale pondérée. Pour tout cluster C e , une fonction de coût constante (d'arité nulle) spécifique au cluster (notée w e ∅ ) est utilisée. Elle peut être incrémentée par le mécanisme de filtrage par cohérence locale et ainsi fournir pour chaque cluster un minorant du coût d'une solution optimale. Au delà d'un meilleur élagage, ceci permet également d'éviter de résoudre systématiquement les sous-problèmes à l'optimalité en transmettant une exigence de coût maximum (un majorant) pour chaque sous-problème P f . Pour P 1 , au départ, le coût de la meilleure solution connue définit le majorant initial cub 1 . Il faut ensuite améliorer ce coût. Si l'on considère un cluster arbitraire C f , fils de C e avec un majorant associé cub e , le majorant du sous-problème P f est obtenu est en soustrayant de cub e les minorants associés à tous les autres clusters fils de C e . Ainsi, des minorants plus forts sur les autres sous-problèmes induit un majorant plus fort pour le problème courant.

Du fait de ce majorant initial, l'optimum de P f n'est pas nécessairement calculé mais un minorant est toujours obtenu. Ce minorant et son éventuelle optimalité sont enregistrés dans LB P f /A f et Opt P f /A f . Ils sont initialement fixés respectivement à 0 et faux.

Pour un élagage amélioré et de meilleurs majorants, BTD utilise le maximum du minorant w e ∅ (fourni par la cohérence locale) et du minorant enregistré (LB P f /A f ). Les minorants enregistrés ne peuvent être utilisés que lorsque le séparateur S f est totalement instancié par l'affectation courante A. Ceci permet récursivement de définir le minorant utilisé dans [START_REF] De Givry | Exploiting Tree Decomposition and Soft Local Consistency in Weighted CSP[END_REF] comme lb(P e /A) = w e ∅ ⊕ C f ∈Fils(Ce) max(lb(P f /A), LB P f /A f ).

Exemple 2 Dans l'exemple 1, les variables {1, 2, 3, 4} de C 1 sont affectées en premier, en utilisant par exemple l'heuristique d'ordonnancement dynamique min domain/max degree. L'algorithme BTD ainsi défini a une complexité temporelle qui, dans le pire des cas, est exponentielle dans la taille du plus grand cluster moins un, un nombre aussi appelé largeur d'arbre de la décomposition arborescente (C, T ). Sa complexité spatiale est, dans le pire des cas, exponentielle en s où s = max Ce∈C |S e |, taille du plus grand séparateur [START_REF] De Givry | Exploiting Tree Decomposition and Soft Local Consistency in Weighted CSP[END_REF].

Soit A = {(4, a), (1, a), (2, b), (3, b)} l'affectation courante. Le filtrage par la cohérence locale EDAC [10] sur P 1 /A produit w 1 ∅ = 2, w 2 ∅ = w 4 ∅ = 1, w 3 ∅ = w 5 ∅ = 0, indui- sant lb(P 1 /A) = Ce∈C w e ∅ =

Poupées russes et décomposition arborescentes

L'algorithme original des poupées russes (RDS) [START_REF] Verfaillie | Russian Doll Search for Solving Constraint Optimization Problems[END_REF], a été introduit en 1996, alors qu'aucune cohérence locale pour les CSP pondérés n'avait été définie. Son mécanisme de base a été élaboré pour construire des minorants forts à partir de minorants plus faibles en utilisant une recherche arborescente. L'application de RDS consiste à résoudre n sous-problèmes imbriqués d'un problème initial P avec n variables. Étant donné un ordre des variables statique, RDS commence par résoudre le sous-problème défini par la dernière variable seule. Ensuite, il inclut la variable précédente dans l'ordre et résout le sous-problème avec deux variables. Ce processus est répété jusqu'à la résolution compléte du problème. Chacun des sous-problèmes est résolu avec un algorithme de type DFBB utilisant un ordre statique des variables suivant l'ordre de la décomposition en sous-problèmes imbriqués. Le minorant amélioré utilisé dans l'algorithme DFBB est obtenu en combinant le minorant faible fourni par "Partial Forward Checking" (similaire au filtrage par cohé-rence de noeud [START_REF] Larrosa | Solving Weighted CSP by Maintaining Arc-consistency[END_REF]) avec l'optimum du sous-problème résolu à l'étape précédente par RDS.

Nous nous proposons d'exploiter le principe de RDS en s'appuyant sur une décomposition arborescente (RDS-BTD). La principale différence avec RDS réside dans le fait que l'ensemble des sous-problèmes résolus est défini par une décomposition arborescente enracinée (C, T ) : RDS-BTD résout |C| sous-problème ordonnés par un parcours en profondeur de T , démarrant aux feuilles jusqu'à la racine P rds L'avantage de l'utilisation de BTD réside aussi dans le fait que les minorants enregistrés peuvent aussi être exploités dans les itérations suivantes de RDS-1 En fait, du fait que le filtrage par cohérence locale peut déplacer des coûts entre les clusters, le minorant LB rds P f doit être ajusté. Ceci est rendu possible par l'utilisation de la structure de données ∆W introduite dans [START_REF] De Givry | Exploiting Tree Decomposition and Soft Local Consistency in Weighted CSP[END_REF] en soustrayant Le pseudo-code de l'algorithme RDS-BTD est présenté ci-aprés comme l'algorithme 1. L'algorithme BTD utilisé est le même que dans [START_REF] De Givry | Exploiting Tree Decomposition and Soft Local Consistency in Weighted CSP[END_REF] Les complexités temporelles et spatiales de RDS-BTD sont celles de BTD. Notez que sans mémorisation, en utilisant uniquement un filtrage par cohérence de noeud et une décomposition induite par un pseudo-arbre (avec un cluster pour chaque variable et un ordre statique des variables), RDS-BTD est équivalent à l'algorithme Pseudo-Tree RDS [START_REF] Larrosa | Pseudo-tree search with soft constraints[END_REF]. Si

L i∈S f max a∈D i ∆W f i (a)). 2 Une valeur a ∈ D i totalement supportée vérifie w i (a) = 0 et ∀w S ∈ W with i ∈ S, ∃t ∈ ℓ(S) with t[i] = a tel que w S (t) = 0.

Résultats expérimentaux

Nous avons codé DFBB, BTD, RDS-BTD, et RDS dans toulbar2 un solveur C++ de CSP pondérés 3 . A l'intérieur des clusters, l'heuristique dynamique de choix de variable min domain / max degree est utilisée (par BTD et RDS-BTD) et par DFBB, les exaequo sont classés par les coûts unaires maximum. Cette heuristique est modifiée par une heuristique de type "conflict back-jumping" comme cela est proposé dans [START_REF] Lecoutre | Last conflict based reasoning[END_REF]. Durant la recherche, un niveau de cohérence locale EDAC [START_REF] De Givry | Existential arc consistency : Getting closer to full arc consistency in weighted CSPs[END_REF] est maintenu à chaque étape. RDS s'appuie sur la cohérence de noeud [START_REF] Larrosa | Solving Weighted CSP by Maintaining Arc-consistency[END_REF] seulement. Les décompositions arborescentes sont construites par l'heuristique "Maximum Cardinality Search" (MCS), en utilisant le plus grand cluster comme racine (excepté pour le problème scen07). À partir de la décomposition arborescente fournie par MCS, nous avons dérivé des décompositions arborescentes alternatives définies par une taille de séparateur maximum s max comme cela est proposé dans [START_REF] Jégou | Dynamic management of heuristics for solving structured CSPs[END_REF] : en partant des feuilles de la décomposition, nous fusionnons un cluster avec son père si la taille du séparateur dépasse s max . Un ordre des variables compatible avec la décomposition arborescente enracinée exploitée est utilisé pour l'établissement de la cohérence locale DAC [START_REF] Cooper | Arc consistency for soft constraints[END_REF] ainsi que par RDS.

Les minorants enregistrés (et si disponibles, fournis par RDS) sont exploités par le filtrage par cohérence locale dés que le séparateur correspondant est totalement affecté. Si le minorant enregistré est optimal ou supérieur au minorant fourni par EDAC, le sous-problème (P e /A e ) correspondant est déconnecté du filtrage par cohérence locale et la différence positive entre les minorants est ajoutée au minorant du cluster père (w Père(Ce)

∅

), conduisant à de nouveaux effacements de valeurs par cohérence de noeud.

Toutes les méthodes s'appuient sur un schéma de branchement binaire. Si d > 10, le domaine ordonné est coupé en deux parties (autour de la valeur centrale) sinon, la variable est affectée à son support complet ou cette valeur est retirée du domaine. Dans les deux cas, la branche contenant la valeur support complet est explorée en premier, excepté pour RDS et les méthodes dérivées de BTD qui choisissent en priorité (si elle existe) la branche qui contient la valeur de la dernière solution trouvée. Sauf mention précise, aucun majorant initial n'est fourni. Les temps CPU présentés correspondent au temps nécessaire pour trouver une solution optimale et prouver son optimalité.

Affectation de fréquences

Parmi les différentes instances du CELAR [START_REF] Cabon | Radio link frequency assignment[END_REF] qui peuvent être décrites comme un CSP pondéré binaire, nous avons sélectionné deux instances difficiles : scen06 and scen07. L'instance scen07 a été réduite par différentes règles de prétraitement (filtrage des domaines par singleton EDAC, par les règles de dominance proposées par Koster [START_REF] De Givry | Singleton consistency and dominance for weighted CSP[END_REF], élimination des variables de degré faible) conduisant à une instance où n = 162, d = 44 et e = 764. La décomposition arborescente utilisée a un s max = 3 et une largeur d'arbre de 53. Un majorant initial (354008) a été fourni à toutes les méthodes. Ce majorant et la solution correspondante ont été trouvés par un algorithme de type DFBB utilisant un mécanisme de type "Limited Discrepancy Search". Cette solution est également utilisée par RDS-BTD comme choix de valeur initial (les performances de DFBB et BTD sont dégradées par ce choix). Le cluster racine a été choisi sur la base du "First Fail Principle" (cluster de coût optimal maximum). BTD et RDS-BTD trouvent l'optimum de 343592 et prouvent l'optimalité en respectivement 6 et 4.5 jours (amélioration de 26% pour RDS-BTD) sur un processeur à 2.6 GHz équipé de 32GB de mémoire. BTD a enregistré 90528 minorants (20% de l'espace des affectations de séparateurs) comparés aux 29453 (6.4%) mémorisés par RDS-BTD. DFBB n'a pas terminé en 50 jours. C'est la première résolution avec preuve d'optimalité de ce problème ouvert depuis plus de 10 ans.

Nous avons aussi résolu l'instance scen06 (100 variables avec une largeur d'arbre de 11) sans aucun prétraitement ni majorant initial. BTD et RDS-BTD ont pris 221 et 316 secondes respectivement pour trouver l'optimum et prouver l'optimalité. Ils ont respective-ment enregistré 5633 et 7145 minorants malgré le fait que la taille des séparateurs n'ait pas été restreinte (on observe un s max = 8). DFBB a pris 2588 secondes et RDS n'a pas terminé en 10 heures.

Sélection de Tag SNP

Ce problème apparaît en génétique et analyse du polymorphisme. Les SNP (ou Single Nucleotide Polymorphism, prononcé snip) sont des variations ponctuelles dans le génome d'individus d'une même espèce. Il se caractérise par une altération d'un seul nucléotide (A,T,C,or G) dans la séquence. Par exemple, un SNP peut changer la séquence D'ADN AAGGCTAA en ATGGCTAA. Pour qu'une variation soit considérée comme un SNP, elle doit apparaître dans au moins 1% de la population considérée. Dans les trois milliards de nucléotides que compte le génome humain, on dénombre plusieurs millions de SNP. Ils expliquent jusqu'à 90% de toutes les variations génétiques humaines, notamment une partie du risque héritable de maladies communes ainsi que la susceptibilité de réponses aux pathogènes, produits chimiques, médicaments vaccins et autres agents.

Le problème TagSNP est une forme de compression d'information avec perte consistant à sélectionner un sous-ensemble de SNP tel que les SNP sélectionnés (appelés tag SNP) capturent l'essentiel de l'information génétique. Le but est de capturer un ensemble de petite taille le plus informatif possible afin de rendre possible le criblage et l'analyse statistique d'une population de grande taille [START_REF] Hirschhorn | Genome-wide association studies for common diseases and complex traits[END_REF].

La mesure de corrélation r 2 entre une paire de SNP peut dans un premier temps être déterminée sur une petite population. Un tag SNP est considéré comme "représentatif" d'un autre SNP si les deux SNP sont suffisamment corrélés. Le problème de sélection de tag SNP le plus simple consiste à sélectionner un nombre minimum de tag SNP de façon à ce que tous les SNP soient représentés. Ceci est capturé par le fait que la mesure r 2 entre deux SNP est supérieure à un seuil θ (souvent fixé à θ = 0.8 [START_REF] Carlson | Selecting a maximally informative set of single-nucleotide polymorphisms for association analyses using linkage disequilibrium[END_REF]). Nous considérons donc un graphe dans lequel chaque sommet est un SNP et où les arêtes sont pondérées par la mesure r 2 entre les SNP des sommets. Les arêtes dont le seuil est plus petit que θ sont supprimées. Le graphe filtré ainsi obtenu peut avoir plusieurs composantes connexes. Le problème TagSNP se réduit alors à un problème de couverture d'ensemble (set covering, NP-dur) sur chacune des composantes. Cette formulation simple du problème a été étudiée avec de bons résultats dans [START_REF] Choi | Efficient Genome Wide Tagging by Reduction to SAT[END_REF] via le compilateur c2d .

En pratique, le nombre de solutions optimales peut être très important et les outils spécialisés du domaines tel que FESTA [START_REF] Qin | An efficient comprehensive search algorithm for tagsnp selection using linkage disequilibrium criteria[END_REF], (en s'appuyant sur deux algo-rithmes incomplets), optimisent, en plus du nombre de tagsnp selection, des critères secondaires : Entre deux tag SNP, une mesure r 2 faible est préférée afin de maximiser la dispersion des tag SNP. Entre un non tag SNP et un de ses représentants, une mesure r 2 élevée est préférée afin de maximiser la représentativité de la sélection de tag SNP effectuée.

Pour un graphe connexe donné G = (V, E), nous construisons un WCSP binaire avec des coûts entiers capturant le problème TagSNP avec les critères secondaires précédents. Pour chaque SNP i, deux variables i s et i r sont introduites. i s est une variable booléenne indiquant si le SNP est sélectionné comme tag SNP ou non. Le domaine de i r est formé par l'ensemble des voisins de i complété avec i lui-même. Il indique le tag SNP qui est chargé de représenter i. Pour un SNP i, des fonctions de coût dures (avec des coûts nuls ou infinis) établissent le fait que i s ⇒ (i r = i). Des fonctions de coût dures similaires établissent (i r = j) ⇒ j s avec les SNP j voisins dans G. Une fonction de coût unaire sur chaque variable i s produit un coût élémentaire U si la variable est à vrai. Le WCSP résultant capture déjà le problème de couverture d'ensembles pur défini par le problème TagSNP.

Pour prendre en compte le critère de représentativité, une fonction de coût unaire est associée avec chaque variable i r . Si i r = i, elle génère un coût non nul égal à ⌊100. Ce problème est similaire à un problème de couverture d'ensembles (set covering) mais avec des coûts additionnels binaires (quadratiques). Ces critères sont ignorés par [START_REF] Choi | Efficient Genome Wide Tagging by Reduction to SAT[END_REF]. Dans cet article, c2d donne une représentation compacte de l'ensemble des solutions du problème de couverture simple, mais le nombre de solutions optimales est si grand (typiquement supérieur à des milliards) que l'application des critères secondaires à des solutions générées par c2d serait trop coûteuse. Comme le mentionnent les auteurs dans leur conclusion, une compilation directe des critères secondaire dans le d-DNNF ne semble pas immédiate. Elle nécessiterait une reformulation du problème sous forme MaxSAT.

Les instances considérées dérivent de données obtenues sur chromosome humain numéro 1, aimable-Fig. 3 -Temps de résolution global et nombre de problèmes résolus en faisant varier la taille maximale du séparateur s max . ment fournies par Steve Qin [START_REF] Qin | An efficient comprehensive search algorithm for tagsnp selection using linkage disequilibrium criteria[END_REF]. Deux valeurs de seuil, θ = 0.8 et 0.5 ont été essayées. Pour θ = 0.8, une valeur usuelle dans le domaine, 43251 composantes connexes ont été identifiées parmi lesquelles nous avons sélectionné les 82 plus grandes. Ces problèmes, qui incluent de 33 à 464 SNP, définissent des WCSP avec des tailles de domaine entre 15 et 224. Ces instances sont relativement faciles. Leur résolution à l'optimalité sélectionne 359 tag SNPs en 2h37' au lieu des 487 SNP sélectionnés en 3' par FESTA (en mode greedy), soit 21% d'amélioration. En mode hybride, FESTA sélectionne 370 tag SNP en 39h17' soit un modeste gain de 3% dans notre cas, mais avec une amélioration d'un facteur 15 en temps. Afin d'obtenir des problèmes plus difficiles, nous avons abaissé le seuil θ à 0.5. Ceci défini 19, 750 composantes connexes parmi lesquelles 516non triviales sont résolues par FESTA par approches incomplètes. Nous avons sélectionnées les 25 plus grandes. Ces instances contiennent de 171 à 777 SNP et ont une densité d'arête qui varie entre 6% et 37%. Elles définissent des instances WCSP dont la taille des domaines varie entre 30 et 266 et contenant de 8000 à 250, 000 fonctions de coûts. La décomposabilité de ces problèmes, estimée par le ratio entre la largeur d'arbre d'une décomposition fournie par MCS (s max = +∞) et le nombre de variables, varie de 14% à 23%.

Tous les problèmes ont été traités avec un majorant initial fourmi par FESTA (mode greedy). L'éxperimentation a été réalisé sur machine équipée de CPU à 2.8 Ghz et 32 GB de RAM. Pour mettre en valeur l'importance de l'utilisation d'une taille de séparateurs bornée (via s max ), nous avons considéré des valeurs de s max allant de 0 (DFBB), 4, 8, 12, 24, 32 à +∞ pour les algorithmes BTD et RDS-BTD. Nous présentons le nombre de problèmes résolus dans une limite de temps de deux heures ainsi que le temps CPU global utilisé (une instance non résolue comptant pour deux heures). RDS présentant de mauvaises performances, les résultats correspondants ne sont pas présentés.

En utilisant s max = 4, notre implémentation améliore le taux de compression de FESTA-greedy par un rapport de 15% (en sélectionnant 2952 tag SNP au lieu de 3477 pour les 516-13 instances résolues). Notez que la différence en temps CPU entre BTD et RDS-BTD augmenterait si une plus grande limite de temps était utilisée. D'un point de vue pratique, le critère utilisé dans le problème TagSNP pourrait être raffinée encore en incluant des informations variées : informations d'annotation de séquences (e.g. préférant les tag SNP apparaissant dans les gènes), mesure entre triplets de SNP comme proposé dans [START_REF] Choi | Efficient Genome Wide Tagging by Reduction to SAT[END_REF] (SNP couverts par des paires de tag SNP). Les bonnes performances de RDS-BTD pourraient permettre de résoudre ces problèmes avec un seuil réaliste de θ = 0.8.

La Figure 4 représente, pour chaque instance (θ = 0.5) résolue au moins une fois lors des expérimentations précédentes, l'évolution individuelle du ratio largeur d'arbre de la décomposition sur taille du problème en fonction de s max . Le succès de la résolution (preuve d'optimalité obtenue dans le temps imparti) est matérialisé par un trait continu, à l'inverse les échecs sont représentés par des pointillés. L'unique résolution de l'instance "17034" (s max =4), est représentée par un carré gris.

Ainsi, on retrouve individuellement les 12 instances résolues pour s max = {4, 8}, mais il est intéréssant de noter que la nature des instances clôturées sont légèrement différentes. L'instance "17034" est résolue uniquement à s max =4, ce qui est compensé par l'instance "14359" , qui est résolue sur l'intervalle 8 et 12. Au total, 13 instances sont clôturées sur l'ensemble des 0 0,25 0,5 Pour la majorité des instances étudiées, reflétant la structure des problémes, les courbes montrent une zone de décroissance importante pour s max ∈ [0, 8], qui est sans doute à l'origine des bons résultats obtenus. La diminution de la largeur d'arbre est ensuite plus faible. Elle atteint progressivement (s max > 12), un plateau de pente faible, qui tends asymptotiquement vers le ratio associé à la décomposition complète du problème. Cette zone des courbes est souvent associée à un échec dans les expérimentations. Ainsi par exemple, l'instance 15757 est résolue pour des valeurs de s max de 4, 8, 12 et 24 avec respectivement des temps de 514, 7, 299 et 3810 secondes. Sur le plan de la décomposition, le ratio de largeur d'arbre passe de 0, 63 pour s max =4 à 0, 30 pour s max = 8 et s max = 12 pour enfin avoisiner 0, 2 à s max = 24. La résolution de l'instance est optimale à s max = 8. Le cas de cette instance n'est pas isolé. Comme le montre la Figure 3, la plupart des instances résolues présentent une valeur s max optimale pour la résolution. Cette situation correspond, probablement, à un bon compromis entre gains issus de la décomposition et pertes liées à une diminution dans la liberté de l'ordre de choix des variables.

Conclusion

L'exploitation pratique de décompositions arborescentes pour résoudre des problèmes d'optimisation combinatoire qui ont une structure n'est pas toujours immédiate. Les problèmes d'optimisation définissent de façon explicite un critère global qui doit être optimisé sur l'ensemble de l'instance. En résolvant des sous-problèmes conditionnellement indépendants, les algorithmes exploitant une décomposition arborescente des problèmes perdent la vision globale sur le critère en exploitant des bornes locales assez faibles. En fournissant des minorants forts pour chacun des sous-problèmes, l'approche des poupées russes permet d'injecter une information plus globale dans chaque résolution de sous-problème au travers d'un majorant initial renforcé, évitant ainsi nombre de résolution inutiles et menant à une efficacité accrue.

Au delà, nos expérimentations montrent que, même sur des problèmes ayant une structure forte, il est souvent profitable et parfois indispensable de restreindre la décomposition en bornant la taille des séparateurs. La théorie nous indique que la taille des séparateurs influe sur la complexité spatiale d'algorithmes tels que BTD ou RDS-BTD, mais même lorsque la taille des séparateurs n'est pas restreinte, aucun des instances considérées dans nos tests n'a pu épuiser la mémoire disponible. En pratique, l'amélioration en efficacité est d'abord explicable par la liberté supplémentaire dans l'ordonnancement des variables rendue possible par la fusion de clusters. Cette observation est cohérente avec les conclusions de [START_REF] Jégou | Dynamic management of heuristics for solving structured CSPs[END_REF]. On notera aussi que le nombre d'instanciations possibles d'un séparateur croît exponentiellement avec sa taille. Plus le séparateur est grand, et plus la probabilité qu'un minorant mémorisé serve à nouveau devient faible : la mémorisation est sans doute surtout utile dans les séparateurs de petite taille.

Concernant notre algorithme, la synergie entre RDS et BTD pourrait être améliorée en autorisant BTD à mémoriser des minorants associés à des instanciations partielles des séparateurs. Cela permettrait de réduire encore les recherches redondantes entre itérations successives de RDS-BTD et serait également profitable dans une stratégie d'approfondissement itératif (iterative deepening) comme cela a été proposé dans [START_REF] Cabon | Anytime Lower Bounds for Constraint Optimization Problems[END_REF].

Dans nos tests, le seuil s max a été appliqué de façon uniforme à toutes les instances. Il reste à étudier à quel point un ajustement de ce seuil par instance, en prenant en compte la taille des séparateurs existants dans le problème et les variations de largeur d'arbre qu'ils peuvent induire, pourrait améliorer les performances de l'algorithme.

  w S : ℓ(S) → [0, m]. L'ensemble des coûts possibles est [0, m] où m ∈ {1, . . . , +∞} représente un coût intolérable. Les coûts sont combinés par l'addition bornée ⊕, définie par a⊕b = min{m, a+b}, et comparés entre eux via ≥. Notez que le coût intolérable m peut être fini ou infini et qu'un coût b peut être soustrait d'un coût a plus large en utilisant l'opération ⊖ où a ⊖ b est égal à (a -b) si a = m . Sinon, il est égal à m.

Fig. 1 -Exemple 1

 11 Fig. 1 -Le graphe d'un WCSP et une décomposition arborescente.

4

 4 étant donné qu'aucun minorant n'a encore été enregistré. Ensuite, le sous-problème P 2 /{(4, a)} et P 4 /{(4, a)} sont résolus de façon indépendante et les solutions optimales correspondantes sont mémorisées dans LB P2/{(4,a)} = 1, LB P4/{(4,a)} = 2, Opt P2/{(4,a)} = Opt P4/{(4,a)} = vrai (aucun solution n'étant connue, aucun majorant n'est défini). Une première affectation complète de coût w 1 ∅ ⊕ LB P2/{(4,a)} ⊕ LB P4/{(4,a)} = 5 est obtenue.

Fig. 2 -Exemple 3 3 , P rds 2 , P rds 5 , P rds 4 et P rds 1 = 3 ,= 1 2 ⊕ LB P rds 4 = 2

 233254131242 Fig. 2 -Résolution récursive des sous-problèmes relaxés issus de la Figure 1 par RDS-BTD.

Algorithm 1 :/∅) ; 1 Soit

 11 RDS-BTD algorithmeFunction RDS-BTD(P rds e ) : [0, +∞] foreach C f ∈ Sons(Ce) do RDS-BTD(P rds f ) ; cube := cub1 -lb(P/∅) + lb(P rds e A l'affectation de Se à des valeurs totalement 2 supportées; LB rds Pe := BTD(P rds e , A, Ve, cube) ; 3 foreach C f descendant de Ce t.q. S f ∩ Se = ∅ do Fixer à faux les Opt P f /A enregistrés, pour tout A ∈ ℓ(S f ) ; return LB rds Pe ; de plus, nous restreignons l'algorithme à l'utilisation d'une décomposition arborescente spécifique (C, T ) telle que |C| = n, ∀e ∈ [1, n], C e = {1, . . . , e}, et ∀e ∈ [2, n], Père(C e ) = C e-1 , alors BTD-RDS se réduit à RDS.

  ⌋. Pour capturer la dispersion entre tag SNPs i et j, une fonction de coût binaire entre les booléens i s et j s est introduite. Elle produit un coût de ⌊100.

  ⌋ lorsque i s = j s = vrai. Le WCSP ainsi obtenu capture simultanément les critères de dispersion et de représentativité. Afin de garder à ces critères leur caractère secondaire, nous utilisons simplement une valeur de coût U (le coût de sélection d'un SNP) suffisamment grande. En pratique, supérieur à la somme des critère secondaire et primaire dans l'instance considérée.

  Ce∈C C e = X) et toutes les fonctions de coût (∀w S ∈ W, ∃C e ∈ C t.q. S ⊂ C e ). De plus, si une variable i apparaît dans deux clusters C e et C g , i doit aussi apparaître dans tout cluster C f sur l'unique chemin reliant C e et C g dans T . est noté Père(C e ) (resp. Fils(C e )). Le séparateur de C e est l'ensemble S e = C e ∩ Père(C e ). Les variables propres d'un cluster C e est l'ensemble V e = C e \ S e .

	3 Décompositions arborescentes et algo-
	rithme DFBB
	Une décomposition arborescente d'un WCSP
	connexe est définie par un arbre (C, T ). Chaque élé-
	ment C e de l'ensemble C = {C 1 , . . . , C k } des som-
	mets de l'arbre est appelé un "cluster". Il est formé
	d'un sous-ensemble de variables (C e ⊂ X). T est
	un ensemble d'arêtes définissant un graphe acyclique
	connexe (un arbre) sur C. L'ensemble des clusters
	C doit couvrir toutes les variables ( Pour un WCSP donné, nous considérons une décom-
	position arborescente enracinée (C, T ) définie par une
	racine notée C 1 . Dans T , le père (resp. les fils) d'un
	cluster C e

  comme le sous-problème défini par les variables propres de C e ainsi que par celles de tous ses clusters descendants dans T et par les fonctions de coût impliquant uniquement des variables propres de ces clusters. P rds e n'a aucun fonction de coût impliquant une variable de S e , le séparateur avec son cluster père, et ainsi son coût optimal est un minorant de P e conditionné par n'importe quel affectation de S e . Ce coût optimal sera donc noté LB rds BTD et non par DFBB. Ceci permet d'exploiter la décomposition et la mémorisation de minorants effectuée dans BTD, offrant ainsi une complexité asymptotique améliorée. Le minorant exploité peut alors s'appuyer sur les minorants LB rds

	Nous définissons P rds	1	= P 1 .
	déjà cal-culés sur des clusters déjà explorés, sur les mi-Pe
	norants fournis par la cohérence locale et sur les
	minorants enregistrés par BTD 1 . Le minorant cor-
	respondant à l'affectation courante A est main-
	tenant défini de façon récursive par lb(P e /A) = w e ∅ ⊕ C P f ),
	qui est évidemment plus fort.
	Dans BTD, la mémorisation n'est effectuée que
	lorsque les séparateurs sont complètement affectés alors que P rds e ne contient pas les variables du sé-
	parateur S e . Nous affectons donc S e avant de résoudre P rds en utilisant les valeurs totalement supportées e (full support) fournies par EDAC [10] 2 de chaque va-
	riable comme des valeurs temporaires utilisées unique-
	ment pour la mémorisation. Une approche alternative
	consisterait à mémoriser les minorants d'affectations
	partielles mais ceci nécessiterait un mécanisme de mé-
	morisation complexe, à considérer dans des travaux
	futurs.

e Pe Chaque sous-problème P rds e est résolu par f ∈Fils(Ce) max(lb(P f /A), LB P f /A f , LB rds

  si ce n'est pour le fait qu'il exploite le minorant RDS amélioré défini ci-dessus. Nous faisons l'hypothèse que la résolution de P rds e avec une affectation initiale A du séparateur S e et un majorant initial cub e obtenu par BTD est efféctuéé par l'appel à BTD(P rds

	RDS-BTD appelle BTD pour résoudre chaque sous-problème P rds e (ligne 3). Un majorant initial pour P rds e est déduit du majorant global et des minorants
	RDS disponibles (ligne 1). Comme expliqué ci-dessus,
	les variables de S e sont affectées à leurs valeurs to-
	talement supportées à la ligne 2 . L'appel initial est
	simplement RDS-BTD(P ). Il suppose que le problème P rds 1 = P est déjà localement cohérent et retourne
	son coût optimal. Notez que dès qu'une solution de P rds e le même coût optimal que lb(P rds e /∅)	ayant =
	f C f ∈Fils(Ce) LB P rds	est trouvée, la recherche
	s"arrête.

e

, A, V e , cub e ).

https ://mulcyber.toulouse.inra.fr/projects/toulbar2 version 0.8.

Acknowledgments Ce travail a été partiellement financé par l'Agence Nationale de la Recherche (projet STALDECOPT).

Références


Aurélie Favier
email: afavier@toulouse.inra.fr

Simon De Givry
email: degivry@toulouse.inra.fr

Philippe Jégou
email: philippe.jegou@univ-cezanne.fr

Comptage de solutions en exploitant la structure du graphe de contraintes

, qui exploitent la structure du graphe de contraintes dans le but de résoudre des CSPs efficacement. Nous proposons d'adapter BTD pour résoudre le problème #CSP. La méthode de comptage exacte résultante a dans le pire des cas une complexité temporelle exponentielle en un paramètre graphique spécifique appelé largeur d'arbre. Pour des problèmes ayant un graphe peu dense mais une grande largeur d'arbre, nous proposons une méthode itérative qui approche le nombre de solutions en résolvant une partition de l'ensemble des contraintes en une collection de sous-graphes partiels triangulés. Sa complexité temporelle est exponentielle en la taille de la plus grande clique (aussi appelé clique number ) du problème d'origine, taille qui peut être beaucoup plus petite que la largeur d'arbre. Des expérimentations sur des problèmes aléatoires CSP structurés et des benchmarks SAT montrent l'efficacité pratique de nos approches.

Introduction

Le formalisme des problèmes de satisfaction de contraintes (CSP) offre un cadre général pour représenter et résoudre de nombreux problèmes. Déterminer si une solution existe est un problème NP-complet. Un problème plus difficile encore consiste à compter le nombre de solutions. Ce problème noté #CSP est connu pour être #P-complet [START_REF] Valiant | The complexity of computing the permanent[END_REF]. Ce problème a de nombreuses applications en informatique, en particulier en IA, par exemple en raisonnement approché [START_REF] Roth | On the hardness of approximate reasonning[END_REF], en diagnostique [START_REF] Satish Kumar | A model counting characterization of diagnoses[END_REF], en révision de croyance [START_REF] Darwiche | On the tractable counting of theory models and its applications to truth maintenance and belief revision[END_REF], comme heuristique pour guider la recherche d'une solution dans les CSPs [START_REF] Kask | New look-ahead schemes for constraint satisfaction[END_REF], ainsi que dans d'autres domaines en dehors de l'informatique, tels que la physique statistique [START_REF] Burton | Nonuniqueness of measures of maximal entropy for subshifts of finite type[END_REF] ou en biochimie pour la prédiction de structures de protéines [START_REF] Mann | Decomposition during search for propagation-based constraint[END_REF].

Dans la littérature, deux principales approches ont été étudiées. D'un côté, les méthodes calculant le nombre exact de solutions et de l'autre côté des méthodes approchées. Pour les méthodes exactes, l'approche naturelle est d'étendre les méthodes systématiques telles que FC [START_REF] Haralick | Increasing tree search efficiency for constraint satisfaction problems[END_REF] ou MAC [START_REF] Sabin | Contradicting conventional wisdom in constraint satisfaction[END_REF] dans le but d'énumérer toutes les solutions. La complexité est bornée par O(m.d n ) avec n le nombre de variables, m le nombre de contraintes, et d la taille maximum des domaines. Il est évident qu'avec cette approche, plus il y a de solutions, plus cela prend du temps pour les énumérer.

Dans ce papier, nous sommes intéressés par des méthodes de recherche arborescente qui exploitent la structure des problèmes, offrant de meilleures bornes sur la complexité temporelle et spatiale. C'est le cas du compilateur de formules en logique propositionnelle en d-DNNF [START_REF] Darwiche | New advances in compiling cnf to decomposable negation normal form[END_REF] et de la recherche dans des graphes ET/OU [START_REF] Dechter | The impact of and/or search spaces on constraint satisfaction and counting[END_REF][START_REF] Dechter | And/or search spaces for graphical models[END_REF] qui font tous les deux du comptage de solutions.

Nous proposons d'adapter l'algorithme Backtracking with Tree-Decomposition (BTD) [START_REF] Jégou | Hybrid backtracking bounded by tree-decomposition of constraint networks[END_REF] à #CSP. BTD avait été initialement proposé pour résoudre des CSPs structurés. Notre modification sur BTD est similaire à celle effectuée dans le cadre de la recherche ET/OU [START_REF] Dechter | The impact of and/or search spaces on constraint satisfaction and counting[END_REF][START_REF] Dechter | And/or search spaces for graphical models[END_REF]. Cependant BTD utilise la notion d'arbre de décomposition de clusters au lieu d'un pseudo-arbre, ce qui conduit naturellement BTD à utiliser un ordre dynamique de choix de variables à l'intérieur des clusters, tandis que la recherche ET/OU utilise plutôt un ordre statique.

La plus part du travail réalisé sur le comptage a été fait sur #SAT, le problème de comptage de modèles associé à SAT [START_REF] Valiant | The complexity of computing the permanent[END_REF]. Les méthodes exactes pour #SAT étendent les résolveurs SAT systématiques, en ajoutant une analyse des composants [START_REF] Bayardo | Counting models using connected components[END_REF] et de la mémorisation [START_REF] Sang | Combining component caching and clause learning for effective model counting[END_REF] pour améliorer les performances.

Les approches qui réalisent une approximation proposent une estimation du nombre de solutions. Elles proposent des algorithmes en temps polynômial ou exponentiel qui doivent fournir des approximations de qualité raisonnable avec des garanties théoriques sur la qualité de l'approximation ou non. De même que pour les méthodes exactes, les principaux travaux sur les méthodes d'approximation ont été fait sur #SAT. Ces approches effectuent soit un échantillonnage dans l'espace de recherche d'origine [START_REF] Wei | A new approach to model counting[END_REF][START_REF] Gomes | From sampling to model counting[END_REF][START_REF] Gogate | Approximate counting by sampling the backtrack-free search space[END_REF][START_REF] Kroc | Leveraging belief propagation, backtrack search, and statistics for model counting[END_REF], soit dans l'espace de recherche ET/OU [START_REF] Gogate | Approximate solution sampling ( and counting) on and/or search spaces[END_REF]. Toutes ces méthodes, à l'exception de [START_REF] Wei | A new approach to model counting[END_REF] procurent un minorant probabiliste du nombre de solutions avec un intervalle de confiance très resserré obtenu par des affectations aléatoires des variables jusqu'à trouver des solutions. Un inconvénient possible de ces approches est qu'elles peuvent ne trouver aucune solution dans le temps de calcul imparti à cause du fait qu'elles rencontrent uniquement des affectations partielles incohérentes. Pour des problèmes complexes de grande taille, cela conduit à des minorants du nombre de solutions nuls. Pour tenter de remédier à ce problème, il faut alors réaliser un réglage délicat des paramètres de ces méthodes, par exemple en changeant le nombre d'échantillons. Une autre approche consiste à réduire l'espace de recherche en ajoutant par exemple des contraintes XOR [START_REF] Gomes | Model counting : A new strategy for obtaining good bounds[END_REF][START_REF] Gomes | Counting CSP solutions using generalized XOR constraints[END_REF]. Cependant, l'ajout de ces contraintes n'assure pas que le problème soit plus facile à résoudre.

Dans ce papier, nous proposons de relaxer le problème de comptage en effectuant un partitionnement des contraintes en une collection de sous-problèmes structurés triangulés. Chaque sous-problème est alors résolu en utilisant notre version modifiée de BTD. Cette tâche devrait être relativement facile si l'instance originale a un graphe peu dense1 . Finalement, une approximation du nombre de solutions du problème complet est obtenue en combinant les résultats obtenus pour chaque sousproblème. La méthode d'approximation résultante, appelée ApproxBTD par la suite, donne aussi un majorant trivial du nombre de solutions. Les résultats expérimentaux montrent qu'une telle approche est intéressante pour sa rapidité et la qualité de son approximation.

D'autres méthodes de comptage fondées sur des relaxations ont été étudiées dans la littérature, telle que l'élimination de variables approchée et la propagation itérative dans un graphe de jointure [START_REF] Kask | New look-ahead schemes for constraint satisfaction[END_REF], ou bien dans le contexte voisin de l'inférence Bayésienne, la propagation de croyance itérative et la méthode de suppression d'arêtes [START_REF] Choi | An edge deletion semantics for belief propagation and its practical impact on approximation quality[END_REF]

= (v 1 , . . . , v k ) sur d x1 × • • • × d x k . Une contrainte c est sa- tisfaite par A si c ⊆ Y, (v 1 , . . . , v k )[c] ∈ r c ,
elle est dite violée sinon. Nous écrirons l'affectation (v 1 , . . . , v k ) sous la forme plus explicite (x 1 ← v 1 , . . . , x k ← v k ). La structure d'un CSP peut être représentée par le graphe (X, C), appelé le graphe de contraintes, dont les sommets sont les variables de l'ensemble X et il y a une arête entre deux sommets s'il existe une constrainte entre les variables correspondantes.

Dans [START_REF] Jégou | Hybrid backtracking bounded by tree-decomposition of constraint networks[END_REF] une nouvelle méthode est proposée pour résoudre les CSP. Cette méthode verse, ApproxBTD commence directement avec un sous-problème triangulé pouvant être plus grand que le poly-arbre de la méthode précédente.

appelée BTD (pour Backtracking with Tree-Decomposition) est une méthode énumérative guidée par une décomposition arborescente du graphe de contraintes. Une décomposition arborescente [START_REF] Robertson | Graph minors II : Algorithmic aspects of treewidth[END_REF] d'un graphe G = (X, E) est une paire (C, T ) avec T = (I, F ) un arbre et C = {C i : i ∈ I} une famille de sous-ensemble de X, tels que chaque cluster C i est un noeud de T et vérifie : [START_REF] Arnborg | Complexity of finding embeddings in a k-tree[END_REF]

∪ i∈I C i = X, (2) Pour chaque arête {x, y} ∈ E, il existe i ∈ I avec {x, y} ⊆ C i , (3) Pour tout i, j, k ∈ I, si k est sur le chemin de i à j dans T , alors C i ∩C j ⊆ C k . La largeur d'arbre d'une décomposition arbo- rescente (C, T ) est égale à max i∈I |C i | -1.
La largeur d'arbre (tree-width) de G est la largeur minimale sur toutes les décompositions arborescentes de G. Notons que trouver une décomposition arborescente optimale est un problème NP-dur [START_REF] Arnborg | Complexity of finding embeddings in a k-tree[END_REF]. Cependant, nous pouvons facilement calculer un bonne décomposition arborescente en utilisant la notion de graphes triangulés. Une décomposition arborescente est calculée par triangulation (i.e ajout d'arêtes au) du graphe de contraintes pour qu'il devienne triangulé3 afin d'y rechercher les cliques maximales dans le graphe de contrainte triangulé. La figure 1(b) présente une décomposition arborescente possible pour le graphe de la figure 1(a). Nous avons

C 1 = {x 1 , x 2 , x 3 }, C 2 = {x 2 , x 3 , x 4 , x 5 }, C 3 = {x 4 , x 5 , x 6 } et C 4 = {x 3 , x 7 , x 8 },
et la largeur d'arbre est de 3. Dans la suite, à partir d'une décompostion arborecente, nous considérons un arbre enraciné (I, F ) où C 1 est la racine, nous noterons F ils(C i ) l'ensemble des clusters fils de C i et Desc(C j ) l'ensemble des variables qui appartiennent à C j ou à un descendant C k de C j dans l'arbre enraciné en C j . Par exemple,

Desc(C 2 ) = C 2 ∪ C 3 = {x 2 , x 3 , x 4 , x 5 , x 6 }.
La première étape de BTD consiste à calculer une décomposition arborescente du graphe de contraintes. Cette décomposition arborescente permet d'obtenir un ordre partiel sur les variables permettant à BTD d'exploiter les propriétés structurelles du graphe et de couper dans l'arbre de recherche. En effet, les variables sont affectées selon une recherche en profondeur d'abord à partir de la racine de la décomposition. Autrement dit, nous affectons les variables du cluster racine C 1 , puis celles de C 2 , celles de C 3 etc. . . Par exemple, x 1 , x 2 , . . . , x 8 est un ordre possible. De plus, la décomposition arborescente et l'ordre sur les variables permettent à BTD de diviser le problème P en plusieurs sous-problèmes. Etant donnés deux clusters C i et C j (avec C j un fils de C i ), le sous-problème enraciné en C j dépend de l'affectation courante A sur C i ∩ C j . On notera ce sous-problème P A,Ci/Cj . Son ensemble des variables est Desc(C j ). Le domaine de chaque variable appartenant à C i ∩ C j est réduit à sa valeur associée dans A. Concernant l'ensemble des contraintes, il contient les contraintes qui impliquent au moins une variable apparaissant exclusivement dans C j ou un de ses descendants. Considérons, par exemple, le CSP dont le graphe de contraintes est donné par la figure 1(a). Chaque domaine est {a, b, c, d} et chaque contrainte Dans la section suivante, nous définirons la notion de good structurel pour le comptage, qui est basée sur les mêmes principes que les goods et nogoods structurels.

c ij = {x i , x j } a une relation r cij telle que x i = x j . Soit A = (x 2 ← b, x 3 ← c), l'ensemble de variables de P A,C1/C2 est Desc(C 2 ), (avec d x2 = {b}, d x3 = {c} et d x4 = d x5 = d x6 = {a

Dénombrement de solutions avec BTD

La méthode BTD est une adaptation de BTD pour le comptage de solutions. Comme pour BTD, la première étape de BTD consiste à calculer une décomposition arborescente du graphe de contraintes. Cette décompostion arborescente induit un ordre partiel sur les variables pour exploiter les propriétés structurelles du graphe et permettre d'élaguer l'arbre de recherche. Tandis que BTD utilise la notion de good, pour le dénombrement de solutions nous utilisons la notion de #good. Le but est de sauver le nombre de solutions des sousproblèmes induit par la décomposition arborescente. En effet, un # good de C i par rapport à C j (avec C j un fils de C i ) est une paire (A, nb) où A est une affectation C i ∩C j et nb le nombre de solutions du sous-problème P A,Ci/Cj . Dans l'exemple de la Section 2, si nous considérons l'affectation A = (x 4 ← a, x 5 ← d) sur C 2 ∩ C 3 nous obtenons un #good (A, 2) car nous avons deux solutions pour le sous-problème P A,C2/C3 . BTD explore l'espace de recherche suivant l'ordre des variables induit par la décomposition arborescente. Ainsi, on commence par les En pratique, pour les problèmes ayant une grande largeur d'arbre, BTD explose en temps et en mémoire, voir à la Section 5. Dans ce cas, nous sommes intéressés par une méthode d'approximation.

Algorithme 1 : BTD(A, Ci, V C i ) : entier if V C i = ∅ then if F ils(Ci) = ∅ then return 1 else F ← F ils(Ci) N bSol ← 1 while F = ∅ et N bSol = 0 do Choisir Cj in F F ← F -{Cj } if (A[Ci ∩ Cj ], nb) n'est pas un #good dans P then nb ← BTD(A, Cj , V C j -(Ci ∩ Cj )) enregistre le #good (A[Ci ∩ Cj ], nb) de Ci/Cj dans P N bSol ← N bSol × nb ; return N bSol else Choisir x ∈ V C i N bSol ← 0 d ← dx while d = ∅ do Choisir v dans d d ← d-{x} if A ∪ {x ← v} ne viole aucune c ∈ C then 1 N bSol ← N bSol+BTD(A ∪ {x ← v}, Ci, V C i -
C i ) tel que A[C i -V Ci ] = B[C i -V Ci ]. V Ci est

Approximation avec BTD

S P ≈ k i=1 S Pi x∈Xi d x × x∈X d x
Notons que l'approximation retourne une réponse exacte si tous les sous-problèmes sont indépendants (∩X i = ∅) ou k = 1 (P est déjà triangulé) ou si il existe un sous-problème incohérent P i . De plus, nous pouvons fournir un majorant trivial du nombre de solutions dû au fait que chaque sous-problème P i est une relaxation de P (le même argument est utilisé dans [START_REF] Pesant | Counting solutions of CSPs : A structural approach[END_REF] pour construire un majorant).

S P ≤ min i∈[1,k] S Pi x∈Xi d x × x∈X d x
Notre méthode appelée ApproxBTD est décrite par l'algorithme 2. Appliqué sur un problème P avec un graphe de contraintes (X, C), Algorithme 2 : ApproxBTD(P) : integer Soit G ′ = (X ′ , C ′ ) un graphe de contraintes associé à

P ; i ← 0 ; while G ′ = ∅ do i ← i + 1 ;
Calculer un sous-graphe partiel triangulé (Xi, Ei) de G ′ ; Soit Pi le sous-problème associé à (Xi, Ei) ;

S P i ←BTD(∅, C ′ 1 , C ′ 1 ) avec C ′ 1 le cluster racine de la décomposition arborescente de Pi ; G ′ ← (X ′ , C ′ -Ei) avec X ′ l'ensemble des variables induites par C ′ -Ei ; k ← i ; return & Q k i=1 S P i Q x∈X i dx × Q x∈X dx ' ;
la méthode construit une partition {E 1 , ..., E k } de C telle que le graphe de contraintes (X i , E i ) est triangulé pour tout 1 i k. Les sousproblèmes associés aux (X i , E i ) sont résolus par BTD. La méthode retourne une approximation du nombre de solutions de P en utilisant la propriété 1.

Le nombre d'itérations de ApproxBTD est borné par n (nous avons au moins n -1 arêtes (un arbre) à chaque itération ou bien des sommets ont été effacés). Chaque sous-graphe triangulé et sa décomposition arborescente optimale associée peuvent être calculé en O(nm) opérations [START_REF] Dearing | Maximal chordal subgraphs[END_REF] 4 . De plus, nous garantissons que la largeur d'arbre w (plus un) de chaque sousgraphe triangulé produit est au plus égale à K, la taille de la plus grande clique (clique number ) de P. Soit w * la largeur d'arbre (minimale) de P, nous avons w + 1 ≤ K ≤ w * + 1. Finalement, la complexité temporelle de Ap-proxBTD est O(n 2 md K ) et sa complexité mémoire est O(nKd K-1 ).

Résultats expérimentaux

Nous avons effectué nos expérimentations sur des CSP aléatoires et sur des benchmarks SAT. Les expérimentations ont été réalisées sur un Pentium IV 3.2 GHz (resp. Xeon 2, 66 GHz) avec 1 Go (resp. 32 Go) pour les instances CSP (resp. SAT). Nous avons limité à une heure le temps autorisé pour résoudre chaque instance. Dans BTD à la ligne 1, nous utilisons Forward Checking au lieu de Backward Checking pour des raisons d'efficacité. A l'intérieur des clusters l'ordre dynamique min domaine / max degré est utilisé pour le choix de variables.

Instances aléatoires de CSP structurés

Nous étudions et comparons ApproxBTD avec BTD sur des instances structurées générées aléatoirement. Pour cela, nous considérons des instances aléatoires de k-arbres partiels avec les paramètres (n, d, r max , t, s max , nc, nr) où chaque graphe de contraintes est un arbre de nc cliques et nr est le pourcentage d'arêtes supprimées dans les cliques. Il y a n variables avec un domaine de taille d, la taille de la plus grande clique est r max , et la taille de la plus grande intersection est s max . Chaque contrainte a une dureté égale à t, le pourcentage de tuples interdits parmi les tuples possibles pour la contrainte.

Les résultats présentés dans les Figures 2 et 3 sont les résultats du nombre de solutions et du temps CPU en millisecondes pour les classes (50, 15, 8, t, 3, 10, 30%) avec la dureté t variant entre 40% et 65%.

Chaque point dans la Figure 2 représente l'évaluation trouvée par ApproxBTD par rapport au nombre de solutions trouvé par BTD pour chaque instance. Nous observons que plus il y a de solutions, meilleure est l'approximation.

Le temps requis pour ces approximations est illustré par la Figure 3. Lorsque la dureté est inférieur à 53%, ApproxBTD est plus rapide que BTD. Pour les plus grandes duretés BTD peut être plus rapide que ApproxBTD car il résoud seulement un problème contre une collection de sous-problèmes.

Benchmarks SAT

Les benchmarks SAT viennent de www. satlib.org. Nous avons sélectionné les instances satisfaisables académiques (Tours de nombre de solutions sur ces instances.

Discussion and conclusion

Dans cet article, nous avons proposé deux méthodes pour résoudre le problème de comptage du nombre de solutions d'un CSP. Ces méthodes exploitent une décomposition arborescente des CSPs. Nous avons présenté une méthode exacte qui est adaptée aux problèmes ayant une petite largeur d'arbre. Pour des problèmes avec une largeur d'arbre importante et un graphe de contraintes peu dense, nous avons présenté une nouvelle méthode d'approximation dont la qualité des résultats obtenus est proche des méthodes existantes mais qui sont obtenus bien plus rapidement et sans avoir à effectuer des réglages complexes de paramètres, à l'exception du choix de l'heuristique de décomposition arborescente.

Fig. 1 -

 1 Fig. 1 -(a) Un graphe de contraintes sur 8 variables. (b) Une décomposition arborescente de ce graphe de contraintes.

  , b, c, d}) et son ensemble de contraintes est {c 24 , c 25 , c 34 , c 35 , c 45 , c 46 , c 56 }. Nous définissons la notion de goods structurels. Un good structurel de C i par rapport à C j (avec C j un fils de C i ) est l'affectation courante A sur C i ∩ C j pouvant être étendue sur le sous-problème P A,Ci/Cj . Par exemple, si nous considérons l'affectation A = (x 2 ← b, x 3 ← c) sur C 1 ∩ C 2 nous obtenons le good (A) grâce à l'affectation sur Desc(C 2 ) défini par (x 2 ← b, x 3 ← c, x 4 ← a, x 5 ← d, x 6 ← b) satisfaisant toutes les contraintes de P A,C1/C2 . Dans un premier temps cette affectation (A) est explorée, puis lorsque son extension sur {x 4 , x 5 , x 6 } est valide, (A) est enregistrée comme un good. Ainsi, si durant la recherche, l'affectation A = (x 2 ← b, x 3 ← c) est à nouveau étudiée, la recherche n'explore pas le sous-problème P A,C1/C2 car nous avons déjà prouvé qu'il existe une solution compatible avec A. Au contraire, si aucune solution n'est trouvée pour une autre affectation A ′ sur P A ′ ,C1/C2 , un nogood structurel est alors enregistré. Ce nogood pourra être utilisé comme une nouvelle contrainte du problème.

  l'ensemble des variables non affectées de C i . Le premier appel est BTD(∅, C 1 , C 1 ) et il retourne le nombre de solutions. Notons S A,Ci/Cj le nombre de solutions d'un sous-problème P A,Ci/Cj . Ainsi, S A,Ci/Cj = B S B,Ci/Cj pour chaque affectation B sur Desc(C j ) tel que A[C i ∩ C j ] = B[C i ∩ C j ]. Par exemple, avec A = (x 1 ← a, x 2 ← b, x 3 ← c) comme affectation de C 1 alors S A,C1/C2 = S A∪(x4←a,x5←d),C1/C2 +S A∪(x4←d,x5←a),C1/C2 = 2 + 2 = 4.Soit Var(P A,Ci/Cj ) l'ensemble des variables de P A,Ci/Cj moins les variables de C i ∩ C j . Autrement dit, c'est l'ensemble des variables de P A,Ci/Cj non affectées. Si C i a k fils : C i1 , ..., C i k alors ∩ 1 j k Var(P A,Ci/Ci j ) = ∅. On note S A,Ci le nombre de solutions de Desc(C i ) avec l'affectation A sur C i . D'où, S A,Ci = 1 j k S A,Ci/Ci j . Par exemple, le nombre de solutions de P A,C1/C2 (avec l'affectation A = (x 1 ← a, x 2 ← b, x 3 ← c) est 4 sur (x 4 , x 5 , x 6 ) et 6 solutions pour P A,C1/C4 sur (x 7 , x 8 ). Donc, il y a 24 solutions sur (x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 , x 8 ) avec A. Notons que pour [C 1 ∩ C 4 ], ((x 3 ← c), 6) est un #good. Pour une autre affectation sur C 1 , e.g A ′ = (x 1 ← b, x 2 ← a, x 3 ← c), il n'est pas nécessaire de calculer les solutions sur C 4 car le #good ((x 3 ← c), 6) peut être utilisé pour A ′ . La complexité en espace de BTD est O(n.s.d s ) et celle en temps est O(n.m.d w+1 ) avec w + 1 la taille du plus grand cluster C k et s la taille de la plus grande intersection C i ∩ C j (C j un fils de C i ).

Fig. 2 -

 2 Fig. 2 -Comparaison entre le nombre de solutions trouvé par BTD et l'estimation trouvée par ApproxBTD sur les CSP structurés aléatoires (les points d'abscisses -0.1 représentent les instances inconsistantes).

Fig. 3 -

 3 Fig. 3 -Temps CPU en millisecondes pour BTD et ApproxBTD sur les CSP structurés aléatoires.

2 . Ces approches ont le défaut de ne

  , . . . , x n } de n variables. Chaque variable x i prend ses valeurs dans le domaine fini d xi de D. Les variables sont soumises à un ensemble C de m contraintes. Chaque contrainte c est définie comme un ensemble {x c1 , . . . , x c k } de variables. Une relation r c (de R) est associée à chaque contrainte c telle que r c représente l'ensemble des tuples autorisés sur d xc 1 × • • • × d xc k . Notons que nous pouvons également définir les contraintes par des fonctions usuelles ou des prédicats, par exemple. Etant donné Y ⊆ X tel que Y = {x 1 , . . . , x k }, une affectation des variables de Y est un tuple A

	pas exploiter la structure locale (aussi appelée
	micro-structure) des instances comme c'est le
	cas pour BTD, grâce à son emploi d'une cohé-
	rence locale et d'un ordre de choix de variables
	dynamique au sein des clusters.
	Dans la section suivante, nous introduisons
	les notations et rappelons les principes de l'al-
	gorithme BTD. La section 3 décrit notre ver-
	sion modifiée de BTD. Dans la section 4, nous
	introduisons la méthode d'approximation Ap-
	proxBTD. Des résultats expérimentaux sur des
	CSPs aléatoires et des benchmarks SAT sont
	présentés dans la section 5. Enfin, nous don-
	nons une conclusion dans la section 6.
	2 Préliminaires
	Un problème de satisfaction de contraintes
	(CSP) est défini par un tuple (X, D, C, R). X
	est un ensemble {x 1

  Si c'est le cas, BTD multiplie le nombre de solutions enregistré avec le nombre de solutions de C i avec A comme affectation. Sinon, on étend A sur Desc(C i ) dans l'ordre pour compter le nombre nb d'extensions consistantes et on enregistre le #good (A[C i ∩C j ], nb). BTD calcule le nombre de solutions du sousproblème induit par le fils suivant de C i . Finalement lorsque chaque fils de C i a été examiné, BTD essaye de modifier l'affectation courante de C i . Le nombre de solutions de C i est la somme des solutions de chaque affectation.

	sur Desc(
	BTD est décrit par l'algorithme 1. Etant
	donnée une affectation A et un cluster C i ,
	BTD regarde le nombre d'extensions B de A

{x})

return N bSol variables du cluster racine C 1 . A l'intérieur du cluster C i , on affecte une valeur à une variable, on "backtrack" si une contrainte est violée. Ce schéma peut être amélioré avec le maintient de l'arc consistance. Lorsque toutes les variables de C i sont affectées, BTD calcule le nombre de solutions du sous-problème induit par le premier fils de C i , s'il en existe un. Plus généralement, considérons C j un fils de C i . Etant donnée une affectation courante A sur C i , BTD vérifie si l'affectation A[C i ∩ C j ] correspond à un #good.

  , E 1 ), . . . , (X k , E k ) tels que ∪X i = X, ∪E i = C et ∩E i = ∅ et tel que chaque sous-graphe (X i , E i ) soit triangulé. Ainsi, chaque (X i , E i ) peut être associé à un sous-problème structuré P i (avec l'ensemble de variables X i et l'ensemble de contraintes correspondant à E i ), qui peut être résolu efficacement en utilisant BTD. Ce partitionnement s'appuie sur le fait qu'il est facile de trouver une décomposition arborescente de largeur d'arbre minimum pour un graphe triangulé (théorème de Fulkerson et Gross, 1965). Supposons que S Pi est le nombre de solutions pour chaque sous-problème P

	cas de CSP binaires dans la présentation de
	la méthode. Nous pouvons définir une collec-
	tion de sous-problèmes en partitionnant l'en-
	semble des contraintes, c'est à dire dans le cas
	de CSP binaires, en partitionnant l'ensemble
	des arêtes du graphe de contraintes. Nous re-
	marquons que chaque graphe (X, C) peut être
	partitionné en k sous-graphes (X 1
	Nous considérons ici des CSPs qui ne sont
	pas nécessairement structurés. Pour les ré-
	soudre, nous proposons d'exploiter BTD en dé-
	finissant une nouvelle méthode d'approxima-
	tion appelée ApproxBTD.
	Sans perte de généralité, nous considérons le

i , 1 i k. Nous estimerons le nombre de solutions de P en exploitant la propriété suivante. Premièrement, nous dénotons P P (A)la probabilité de A est une solution de P . P P (A) = S P x∈X d x . Propriété 1. Soit un CSP donné P = (X, D, C, R) et une partition {P 1 , ..., P k } de P induite par une partition de C en k éléments.

Table 1

 1 résume nos résultats. Les colonnes sont le nom de l'instance, le nombre de Tab. 1 -Solution counting for SAT instances. Time in seconds. A "-" means the instance was not solved in less than 1 hour. variables booléennes, le nombre de clauses, la largeur d'arbre de la décomposition arborescente, le nombre exact de solutions, le temps CPU en secondes pour c2d et BTD, pour Ap-proxBTD : la largeur d'arbre maximale pour tous les sous-problèmes triangulés, le nombre de solutions approché et le temps, enfin pour SampleCount : le minorant du nombre de solutions et le temps. Nous remarquons que BTD peut résoudre des instances ayant une petite largeur d'arbre. c2d est généralement plus rapide (excepté pour hanoi5) mais il souffre également lorsqu'on a une grande largeur d'arbre (e.g. ais). Notre méthode d'approximation Ap-proxBTD exploite une partition du graphe de contraintes tels que les sous-problèmes résultants aient une petite largeur d'arbre (w ≤ 11) comme le montrent les résultats. En pratique la méthode permet d'obtenir des résultats relativement rapidement même si la largeur d'arbre originale est importante5 . La qualité de l'appoximation trouvée par ApproxBTD est relativement bonne et est comparable à Samplecount qui prend plus de temps et requiert un réglage de paramètres (nous utilisons les paramètres avec t = 7, s = 20, α = 1 les options par défaut). Le majorant calculé par ApproxBTD n'est pas mentionné car il est très supérieur au5 Une petite largeur d'arbre ne permet pas en pratique d'être rapide à chaque fois (voir e.g. les instances ssa )

	Instances	Vars Clauses	w Solutions	c2d	BTD		ApproxBTD	SampleCount
						Temps	Temps	w	Solutions	Temps	Solutions	Temps
	académiques
	ais6	61	581	41	24	0.04	0.06	5	≈ 1	0.09	≥ 6	0.184
	ais8	113	1520	77	44	0.51	2.34	7	≈ 1	0.52	≥ 8	41.31
	ais10	181	3151	116	296	17.14	390.32	9	≈ 1	2.69	≥ 38	384.8
	ais12	265	4269	181	1328	1162.75	-	11	≈ 1	2.69	≥ 0	17.6
	hanoi4	718	4934	46	1	3.41	1.72	6	≈ 1	1.57	≥ 0	5.26
	hanoi5	1931	14468	58	1	-	25.46	7	≈ 1	14.98	≥ 0	6.19
	industrielles
	ssa7552-038	1501	3575	25	2.84e40	0.15	0.72	7 ≈ 6.34e35	1.36	≥ 6.67e38	1763
	ssa7552-158	1363	3034	9	2.56e31	0.10	0.19	5 ≈ 2.59e27	0.79	≥ 3.57e29	175
	ssa7552-159	1363	3032	11	7.66e33	0.09	0.25	5 ≈ 1.09e30	0.84	≥ 1.62e32	237
	ssa7552-160	1391	3126	12	7.47e32	0.12	0.29	5 ≈ 2.88e33	0.99	≥ 2.02e31	1117
	2bitcomp 5	125	310	36	9.84e15	0.47	11.53	6 ≈ 2.23e15	0.02	≥ 1.80e15	1.42
	2bitmax 6	252	766	58	2.10e29	18.71	-	7 ≈ 1.98e28	0.1 ≥ 1.02e28	8.42

En fait, cela dépend de la largeur d'arbre des sousproblèmes, qui est bornée par la taille de la plus grande clique dans le problème complet. Dans le cas de graphes peu denses, nous nous attendons à ce que cette taille soit petite, ce qui est montré dans nos expérimentations. En pratique, notre méthode d'approximation ne sera efficace que si le CSP à résoudre ne contient pas de contraintes

globales.[START_REF] Bayardo | Counting models using connected components[END_REF] Cette méthode commence par résoudre un sousproblème structuré en poly-arbre, par la suite augmenté par la restauration d'arêtes supprimées, jusqu'à au final résoudre le problème complet. A l'in-

Un graphe est triangulé si chaque cycle de longueur supérieur à quatre a une corde, i.e une arête reliant deux sommets non consécutifs dans le cycle.

Notons que cette approche retourne un sousgraphe maximal au sens de l'inclusion en nombre de contraintes pour un CSP binaire seulement. Dans le cas d'un CSP non-binaire (également en SAT), nous ne garantissons pas la maximalité du sous-graphe et ajoutant dans le sous-problème toutes les contraintes / clauses dont la clique associée d'arêtes est totalement incluse dans le sous-graphe.

Remerciements

Ce travail est partiellement supporté par le projet ANR STAL-DEC-OPT.


Philippe Blanc

Lucien Wald

IMPORTANCE AND EFFECT OF CO-REGISTRATION QUALITY IN AN EXAMPLE OF « PIXEL TO PIXEL » FUSION PROCESS. (1,2)

By definition, as they require the images to merge to be superimposable, « pixel to pixel » fusion processes raise up the issue of the outcome of geometric co-registration quality on the resulting images. Indeed, local errors of geometric co-registration may introduce local errors in the merging process. Hence, we can wonder how can we define and reach the coregistration quality level that is required to minimise errors in « pixel to pixel » merging processes. The aim of this work is to understand, thanks to an example of such fusion process applied to images with different degrees of geometric distortions, the outcome of coregistration error on the fusion quality and, therefore, to measure the efficiency of some coregistration methods for the pre-processing of « pixel to pixel » merging application.

1.

INTRODUCTION

The paper deals with the issue of co-registration errors between the images to merge within a « pixel to pixel » fusion process. In other words, the aim of this work is, on the one hand, to underline and to quantify the outcome of relatively small geometric distortions on the fusion quality and, in on the other hand, to define (and reach) a co-registration quality level required for « pixel to pixel » fusion process. We have based our study on a particular fusion process: the ARSIS method. Indeed, this method is a « pixel to pixel » fusion process that is meant to increase the spatial resolution of an image (e.g. SPOT multispectral image XS) without any modification of its spectral contents by merging structures extracted from a higher resolution image of the same scene but in a different spectral band (e.g. SPOT Panchromatic image P). For more information about ARSIS, see [START_REF] Mangolini | Procédé et dispositif pour augmenter la résolution spatiale d'images à partir d'autres images de meilleure résolution spatiale, French patent n° 92-13691[END_REF](1993). It is important to note that the quality of such synthetic enhanced multispectral images can be efficiently assessed by the mean of the approach described in Wald et al. (1997) and in [START_REF] Mangolini | Evaluation de la qualité des images multispectrales à haute résolution spatiale dérivées de SPOT, Compte-rendus du colloque « Qualité de l'interprétation des images de télédétection pour la cartographie », Grignon, 1-3[END_REF]. This approach has been adopted here. More precisely, SPOT XS images are down-sampled to 40 m and then synthesised by the mean of the ARSIS method at 20 m. The synthesised image XS * are finally compared to the original image XS. In this paper, we deal with a SPOT image of the city of Barcelona, Spain. Multispectral (XS: resolution 20 m) and Panchromatic (P: resolution 10 m) images are available. The P image was acquired with an angle of 15° along the track, while the XS image was acquired at the nadir of the satellite (see Anonymous 1986). The P image has been geometrically corrected by the CNES to be superimposable to the XS one. A technique of sub-pixel correlation has been performed onto the XS image and the P image degraded at 20 m (P 20 ): for each pixel of XS1, the best correlation with P 20 pixels is looked for at sub-pixel level. Hence, this method provides a sub-pixel geometric deviation between each XS1 pixel with the corresponding P 20 pixel. After filtering in order to remove the bad correlation matches, a realistic and complex field of local geometric distortions was obtained. Table 1 displays the mean and the standard deviation of those residual distortions in row and column directions. These values are relatively small but it should be noted that the field of distortions exhibits structures which are correlated to the orographic features of the scene. Here, for our demonstration, we have used mostly XS1 and XS2 image. The P image was only used to produce a field of typical and realistic local distortions. In the following, XS1 is assumed to be the image with the highest resolution (20 m) while XS2 has been downsampled to 40 m to be the lowest resolution image noted XS2 40 .

∆x ∆y

Mean

More precisely, the study is carried out in two phases. The first one is meant to underline the influence of small geometric distortions between the images to merge with the ARSIS method. The second phase consists in the assessment of the benefits of a particular coregistration method as a pre-processing of the fusion.

INFLUENCE OF LOCAL GEOMETRIC DISTORTIONS.

The ARSIS method was performed on the XS1 image (resolution 20 m) and the downsampled XS2 40 image. It provided a synthesised image M1_XS2 * (resolution 20 m) which was compared, visually and quantitatively, to the original XS2 image.

Figure 2 shows the original XS2 image (a) and the synthesised M1_XS2 * image (b). One can note that there is no visual differences between the two images. Moreover, in the lines M1 of Tables 2 and3 are reported statistical criteria which quantify the differences between XS2 and M1_XS2 * following the work of Wald et al. (1997). Those values are very near to the ideal values, thus confirming the visual inspection and supporting the ARSIS method as an efficient one. This comparison between XS2 and M1_XS2* is meant to be a reference for the following demonstration.

To measure the influence of the local geometric distortions, the XS1 image was re-sampled into a new image S_XS1 image by injecting the realistic field of residual geometric distortions discussed in introduction. It means that S_XS1 and XS2 are close to be superimposable but small, realistic and spatially variable geometric distortions exist and are perfectly known at each point. The ARSIS method was then performed on S_XS1 and XS2 40 images to produce a second synthetic image M2_XS2 * to be compared to M1_XS2 * . This first phase of the study is illustrated in Figure 1. The quality of the M2_XS2* image can be firstly analysed and compared to XS2 and M1_XS2 * by a visual inspection (see Figure 2 c). This image clearly presents artefacts near discontinuities and structures that degrades the quality of the synthesised image. In other words, on account of the local geometric distortions, the ARSIS method injects into XS2 40 non-relevant high resolution structures from S_XS1 image because they are slightly shifted. The statistic criteria given in the lines M2 of the Tables 2 and 3 corroborate the visual inspection: according to the correlation coefficient, the standard deviation of the difference and the relative error probability, the small geometric distortion field degrades notably the quality of the results of the ARSIS method. Nevertheless, we have to note that the difference of variances, that is a measure, to some extent, of the quantity of information added (negative value) or lost (positive value) during the enhancement, is better for M2 than for M1. Indeed, in the M1 case, the ARSIS method slightly increases the amount of information (about 2.3 %). This increase of information is less important for M2 (about 0.4 %). This can be explained, in part, by the fact that, in this case, the ARSIS method injects high frequency structures extracted from S_XS1 that have been smoothed by the re-sampling. Table 3. Probability (in percent) for having in a pixel a relative error less than or equal to the thresholds noted in the first row for the different merging processes. The ideal value is 100 as early as the first threshold 0.001 %.

The comparisons between the lines M1 and M2 in Tables 2 and3 demonstrate the noticeable influence of small local geometric distortions. As a conclusion, small geometric residual distortions that are usually encountered in images made superimposable by standard methods have a noticeable influence on the quality of products resulting from fusion processes. Only the ARSIS method was applied here, but it is believed that the effects are similar for any « pixel to pixel » fusion process.

THE BENEFIT OF AN ANDVANCED CO-REGISTRATION METHOD

This section explores the potentialities of an advanced co-registration method to correct residual distortions and ,thus, providing better fusion products of better quality. For this work, we applied a co-registration method to the image S_XS1 to produce a re-sampled image R_XS 1 the most superimposable to XS 1 as possible. Once applied to this R_XS1 and XS2 40 images, the ARSIS method produces another synthetic image noted M3_XS2 * that will be compared to the original image XS2. This second phase of the study is illustrated in Figure 3. Figure 3. Illustration of the assessment of the benefit of the co-registration method as a pre-processing of the ARSIS method.

The applied co-registration method applied is a fully automatic registration method based on multiresolution analysis and local distortion model. This method is in fact an improved version of a co-registration method developed by [START_REF] Djamdji | A new method of automatic registration, based on a multiresolution decomposition of the images using the wavelet transform[END_REF]. The control points are detected at the lower resolution where the smoothed images present only few structures. The catching between the two sets of control points provides a first estimation of the deformation model. Then this estimation is refined as the resolution gets finer till the original resolution.

Validation tests have shown that this method provides a large number of good matched control points. This enables the use of local geometric distortion models. In other words, this method is well fitted to correct small and local geometric distortions.

In our case, we perfectly know the distortions between the two images S_XS1 and XS1. Therefore, we can accurately measure the co-registration error after the geometric correction. Table 4 gives some statistical criteria to describe it and can be compared to Table 1. First, we notice that the bias and the standard deviation resulting from this co-registration Compared to the Table 1, the gain in bias is spectacular. Nevertheless, the bias and the standard deviation are not totally relevant to describe local co-registration errors. Two other statistic criteria are then used:

∆x
• the difference of variances: this value is, to some extend, the difference of the quantity of information between the actual distortions and the estimated ones. • The correlation coefficient: it shows the similarity in shape between the actual and the estimated distortions. Those criteria are close to the ideal value. It means that the co-registration process has well corrected the local geometric distortions between S_XS1 and XS1. Let examine, now, the synthetic image M3_XS2 * . As far as visual analyses are concerned (Figure 4), this image does not present the same artefacts than in M2_XS2 * . Moreover, except for the difference of variances (for the same reason than in M2 case), the statistic criteria given in Tables 5 and6, show that the level of quality of the fusion is equivalent to the M1 case. In other words, the co-registration method has been almost efficient enough to correct the local geometric deviation for the fusion process.

CONCLUSION

We have studied the influence of geometric distortion in a very realistic, typical and controlled case. Despite the weakness of those small geometric distortions mean and standard deviation are, in this case, respectively less than 0.3 and 0.1 pixel, it was shown that it noticeable on the quality of products resulting from a "pixel to pixel" fusion process as far as visual and quantitative assessments are concerned.

The example shown here is based on SPOT images and ARSIS fusion method. It is important to note that those conclusions can be extended to other sensors and other fusion methods such as IHS method (see [START_REF] Harris | IHS transform for the integration of radar imagery with other remotely sensed data[END_REF] or the high frequency injection technique developed by Munechika et al. (1993). Moreover, we have pointed out the need for a better accuracy in co-registration prior to "pixel to pixel" fusion processes. An advanced method of automatic co-registration using multiresolution analysis and local geometric distortion models has been successfully tested.

It constitutes an answer to the need of high quality co-registration underlined previously.

REFERENCES:

Anonymous, 1986. Guide des utilisateurs de données SPOT, 3 volumes, Editors CNES and SPOT Image, Toulouse, France.

Figure 1 .

 1 Figure 1. Illustration of the assessment of the influence of geometric distortion errors in the ARSIS method.

  Figure 2. (a) Sub-region of the original XS2 image. (b) Same but for the synthesised M1_XS2 * image. (c) Same but for the synthesised M2_XS2 * image.

Figure 4 .

 4 Figure 4. As Figure 2, but for the synthetic image M3_XS2 * .

Table 4 .

 4 Means, standard deviations in pixel for the error of the co-registration method in the column and row directions (respectively ∆x and ∆y). Correlation and difference of variances are also reported (see text).

	∆y

Table 5 .

 5 As Table2, but for the ARSIS merging process M3.

		0.001	1	2	5	10	20
	M3	28	29	53	91	99	100

Table 6 .

 6 As Table3, but for the ARSIS merging process M3.