Démixage en présence de pixels rares

Le démixage hyperspectral, un problème d’optimisation

 Dans le cadre du démixage d’images hyperspectrales on cherche à estimer deux matrices A ∈ R P ×k et S ∈ R k×L telles que : Y = AS + N (4.1) S ≥ 0, A ≥ 0, A1K = 1P Ce problème d’estimation peut-être réécrit sous la forme d’un problème d’optimisation dont l’objectif est de minimiser une fonction de coût. Une fonction de coût souvent utilisée est la norme de Frobenius, ou encore l’erreur quadratique de reconstruction (EQR) entre l’observation Y et sa factorisation AS. J(A, S) = EQR(A, S) = kY − ASk 2 (4.2) Lorsque l’on se place dans l’hypothèse qui est la notre d’un bruit gaussien, minimiser l’erreur quadratique de reconstruction entre Y et AS est optimal [88]. 41 Démixage en présence de pixels rares .Nous cherchons alors à minimiser la fonction de coût (4.2) en fonction de A et S, sous contrainte de positivité et de somme à l’unité. Le problème d’optimisation contraint devient donc : minimiserA,S J(A, S) (4.3) sous les contraintes S ≥ 0, A ≥ 0, A1K = 1P Cependant la minimisation de J(A, S) sous les seules contraintes de nonnégativité des matrices A, et S et de somme à l’unité n’a pas de solution unique [44]. En effet, il existe des solutions dégénérées à ce problème. Par exemple si le couple A ∈ R P ×k + et S ∈ R k×L + est solution de l’équation (4.3), alors toute matrice carrée inversible P ∈ R k×k + telle que AP respecte la contrainte de somme à l’unité, permet d’obtenir une infinité de couples solution tels que : Y = AS + N = AP|{z} A0 P −1S | {z } S0 +N (4.4)

Contrainte de somme à l’unité relaxée

D’un point de vue théorique la somme à l’unité est une contrainte élégante et naturelle, puisqu’elle revient à dire que les abondances correspondent aux proportions de chaque matériau présent dans le pixel, et qu’en conséquent la somme de ces proportions doit être égale à 1. Cependant, en pratique un matériau avec une réflexivité qui varie en fonction de l’angle d’incidence aura un unique endmember dont l’intensité dépendra de son orientation. De même, la mesure du spectre d’un matériau situé dans une zone ombragée différera du spectre mesuré pour le même matériau illuminé directement. Ce phénomène appelé variabilité d’endmember [28, 89] fait que l’hypothèse de somme à l’unité n’est pas toujours exacte. De ce fait on préférera appliquer cette contrainte de manière moins stricte en favorisant plutôt les solutions dont le résultat est proche de cette contrainte à l’aide d’un terme de régularisation : J2(A, S) = kY − ASk 2 + α 2 kA1K − 1P k 2 (4.5) Le paramètre α permet d’imposer de façon plus ou moins stricte la contrainte. Une valeur de α élevée indique une contrainte imposée strictement alors qu’une valeur nulle indique que la contrainte n’est pas imposée. L’équation (4.5) se dé42 Démixage en présence de pixels rares veloppe alors comme suit : J2(A, S) = kY − ASk 2 + α 2 kA1K − 1P k 2 = tr  YTY − 2S TATY + S TATAS + α 2 tr  1 T P1P − 21 T KAT1P + 1 T KATA1K  = tr  (YTY + 1 T P1P) − 2(YST + α 21P1 T K)AT + (SST + α 21K1 T K)ATA  = tr  [Y|α1P] T [Y|α1P] − 2[Y|α1P] T [S|α1K]AT + [S|α1K] TATA[S|α1K]  = tr  YT c Yc − 2S T c ATYc + S T c ATASc  = kYc − ASck 2 où tr(·) est l’opérateur trace et : Yc = [Y|α1P] et Sc = [S|α1K] (4.6) sont les matrices Y et S auxquelles ont a ajouté une colonne de valeur constante α. On remarque que la fonction de coût J2 donnée par l’équation (4.5) est identique à la fonction de coût J de l’équation (4.2) dans laquelle on aurait remplacé les matrices Y et S par Yc et Sc. En conséquent, les méthodes qui permettent de résoudre le problème d’optimisation donné en (4.3) sans contrainte de somme à l’unité, permettent également de résoudre le problème avec une contrainte de somme à l’unité relaxée. 

Démixage par la NMF

 La factorisation en matrices non-négatives (NMF) est une méthode de démixage qui consiste à factoriser une matrice non-négative Y ∈ R P ×L + en deux matrices facteurs A ∈ R P ×k et S ∈ R k×L tels que AS ≈ Y pour un rang de factorisation k donné, sous contrainte que tous les coefficients de A et de S soient positifs ou nuls. La NMF vise à résoudre le problème d’optimisation suivant : minimiserA,S kY − ASk 2 (4.7) sous les contraintes S ≥ 0, A ≥ 0 43 Ce problème d’optimisation est identique à celui du démixage évoqué dans l’équation (4.3) lorsque la contrainte de somme à l’unité n’est pas appliquée strictement. Le problème ainsi posé est non-convexe, ce qui rend la recherche d’un minimum global difficile. Cependant si l’une des deux matrices facteur A ou S est connue, alors le problème devient convexe. Une méthode de résolution consiste alors à séparer le problème en deux sous-problèmes convexes. L’un consistant à estimer A en supposant S fixée et l’autre à estimer S en supposant A fixée. Ces deux sous problèmes s’expriment sous la forme : min S≥0 kY − ASk 2 (4.8) min A≥0 kY − ASk 2 (4.9) Il s’agit là de deux problèmes des moindres carrés. De nombreux algorithmes ont été proposés pour les résoudre. Leur adaptation à la NMF est basée sur la méthode des moindres carrés alternés (ALS) qui consiste à résoudre itérativement chacun des deux sous-problèmes. Elle est donnée dans l’algorithme 1. Algorithme 1 NMF – ALS /* Initialisation */ l ← 1 A(l) ≥ 0 S (l) ≥ 0 /* Jusqu’à convergence */ tant que EQR(A(l) , S (l) ) n’a pas convergé répéter A(l+1) ← minA≥0 EQR(A, S (l) ) S (l+1) ← minS≥0 EQR(A(l+1) , S) l ← l + 1 fin tant que Sortie : A(l) , S (l) 4.1.2.1 Problème des moindres carrés Le problème des moindres carrés consiste à résoudre le problème d’optimisation suivant : minimiserX kY − AXk 2 (4.10) Où X ∈ R k×L . La méthode NMF nécessite cependant que les matrices soient positives. Une façon de garantir ce critère est de projeter dans R k×L + le résultat 44 Démixage en présence de pixels rares obtenu par la méthode des moindres carrés en annulant les coefficients négatifs avec l’opérateur [ · ]+. Plusieurs méthodes utilisent cette approche : Algorithme multiplicatif : Lee et Seung furent parmis les premiers à populariser la méthode NMF . Ils proposent un algorithme de NMF basé sur un schéma d’algorithme des moindres carrés alternés qui utilise des règles multiplicatives pour mettre à jour les matrices A et S à chaque itération [63]. A ← h YST (SST ) −1 i + (4.11) S ← h (ATA) −1ATY i + (4.12) La démonstration permettant d’aboutir aux équations (4.11) et (4.12) est détaillée dans [90]. La convergece de cet algorithme n’est cependant pas garantie [91, 92]. Algorithme du gradient projeté : L’algorithme du gradient projeté, PG (pour Projected Gradient) [62], utilise un processus de mise à jour additif des matrices A et S. A partir de la fonction de régularisation EQR(A, S) = kY − ASk 2 , chaque itération de l’algorithme du gradient projeté peut être résumée ainsi : A ←  » A − µA ∂ ∂A EQR(A, S) # + (4.13) S ←  » S − µS ∂ ∂S EQR(A, S) # + (4.14) La démonstration permettant d’aboutir aux équations (4.13) et (4.14) est présentée dans [90]. µA et µS sont les pas du gradient. Ces pas peuvent être fixes ou variables en fonction des algorithmes dérivant du gradient projeté et auront une influence sur la précision ou la rapidité de convergence de l’algorithme. Enfin les gradients de la fonction de régularisation EQR(A, S) par rapport à A et S sont donnés dans [90] : ∂ ∂A EQR(A, S) = 2(AS − Y)S T (4.15) ∂ ∂S EQR(A, S) = 2AT (AS − Y) (4.16) D’autres méthodes de résolution du problème des moindres carrées ont également été proposées, en utilisant par exemple une approche quasi-newtonienne [93, 94]. Pour garantir les contraintes de positivité de la NMF, certaines méthodes utilisent le problème de moindres carrés contraints.

Problème des moindres carrés contraint

 Le problème des moindres carrés contraint vise à minimiser la même fonction de coût que la version non-contrainte (4.10), en ajoutant une contrainte de positivité : minimiserX kY − AXk 2 (4.17) sous la contrainte X ≥ 0 L’utilisation d’une méthode de résolution des moindres carrés contraint pour résoudre les équations (4.8) et (4.9) permet d’obtenir à chaque itération la solution optimale. Plusieurs méthodes ont été développées pour résoudre ce problème [95, 96, 97]. Dans cette thèse nous avons principalement utilisé la méthode BPP (Block Principal Pivoting) . Remarque : L’équation (4.17) peut être décomposée colonne par colonne : min X≥0 kY − AXk 2 F = X j min x≥0 ky:j − Ax:jk 2 F (4.18) où y:j et x:j sont respectivement la je colonne de Y et X. Le problème des moindres carrés contraint revient donc à résoudre l’équation suivante : minimiserx ky − Axk 2 (4.19) sous la contrainte x ≥ 0 Algorithme Bloc Principal Pivoting : Résoudre le problème de moindre carrés contraint vecteur par vecteur n’est pas efficace d’un point de vu calculatoire, la résolution directe sous forme matricielle (4.17) est plus efficace. Cependant nous nous contenterons de présenter la méthode BPP pour résoudre le cas restreint de l’équation (4.19). Le lecteur 46 Démixage en présence de pixels rares intéressé par les améliorations que permet la résolution sous forme matricielle peut se référer à [64]. Les conditions d’optimalités de Karush-Kuhn-Tucker (KKT) pour l’équation (4.19) sont les suivantes : z = ATAx − AT y (4.20) z ≥ 0 (4.21) x ≥ 0 (4.22) xizi = 0 i = 1, …, k (4.23) Les conditions KKT sont des conditions nécessaires d’optimalité. Elles sont suffisantes lorsque la matrice A est de rang plein, ce qui est le cas pour la NMF. Toute solution qui satisfait les conditions KKT est donc une solution optimale de l’équation (4.19). La méthode BBP s’inspire des méthodes dites Active-Set [94], les variables sont divisées en deux groupes F et G complémentaires avec F ∪ G = {1, …, k} et F ∩ G = ∅. On note xF , zF et xG, zG les sous vecteurs correspondant aux indices dans F et G respectivement (c’est à dire, xF = {xi |i ∈ F} ). De même on note AF et AG les sous matrices de A dont seules les colonnes d’indice dans F (respectivement dans G) ont été conservées. Les vecteur xF et zG sont initialisés par un vecteur nul et seront mis à jour successivement par les formules : xF = min xF kAF xF − yk 2 (4.24) zG = AG(AF xF − y) (4.25) Si xF ≥ 0 et zG ≥ 0 alors le vecteur x = (xF , 0) est une solution optimale de l’équation (4.19). Dans le cas contraire, les groupes F et G vont être mis à jour pour échanger les variables qui ne respectent pas les conditions KKT (4.21) et (4.22). Les règles de mises à jour sont les suivantes : F = (F − H1) ∪ H2 (4.26) G = (G − H2) ∪ H1 (4.27) avec H1 = {i ∈ F|xi < 0} et H2 = {i ∈ G|zi < 0} Les vecteur xF et zG sont alors mis à jour avec les nouveaux groupes et le processus est répété jusqu’à ce que le critère d’arrêt xF ≥ 0 et zG ≥ 0 soit rempli.