DYNAMIQUE ET MICROPHYSIQUE DES SYSTEMES CONVECTIFS DE L’AFRIQUE DE L’OUEST

Complément sur la reconstruction par régularisation quadratique (RA)

Dans ce chapitre, on s’intéresse à la reconstruction d’images hyperspectrales par une régularisation quadratique de type Thikonov [TA77] préservant les contours, comme présenté section 1.4 et étudiée en détail dans l’article du chapitre 3. La solution est donc définie par : ob = arg min o ||m − Ho||2 Γ−1 + µx||Dxo||2 + µy||Dyo||2 + µλ||Dλo||2 . (4.1) La solution de ce problème de minimisation d’un critère quadratique possède une expression analytique : (HtΓ −1H + µxDt xDx + µyDt yDy + µλDt λDλ)ob = HtΓ −1m. (4.2) Nous présentons § 4.1 une implémentation efficace du calcul de cette solution par un algorithme de type gradients conjugués, l’algorithme CGNR, dont chaque itération a un coût de calcul faible si l’on prend en compte les caractéristiques de l’instrument. Nous étudions ensuite § 4.2 la vitesse de convergence de l’algorithme pour le dispositif utilisé et nous proposons des solutions permettant d’accélérer la convergence de l’algorithme sans pour autant augmenter son coût de calcul. Nous proposons en particulier d’exploiter certaines configurations du DMD permettent d’améliorer le conditionnement des matrices entrant en jeu dans cet algorithme et ainsi accélérer la convergence de l’algorithme, en particulier dans le cas d’un modèle de bruit blanc : les configurations orthogonales normalisées.

Implémentation de la méthode RA par gradients conjugués (CGNR)

Implémentation de la méthode RA par gradients conjugués (CGNR)

Bien que la reconstruction du cube HS par inversion directe de la matrice dans l’équation (4.2) soit envisageable, elle n’est cependant pas possible dans notre cas à cause des dimensions de la matrice à inverser. En effet, la taille de celle-ci pour un cube de R×C ×W est de RCW ×RCW, valeur supérieure à 1013 éléments pour un cube de dimension 200 × 200 × 100, ce qui requiert plusieurs gigaoctets d’espace mémoire de stockage, et ce pour des cubes de faibles dimensions. L’inversion directe nécessite alors une quantité de mémoire et un temps de calcul prohibitifs. Il est immédiat d’écrire cette solution comme solution d’une équation normale et exploiter des algorithmes itératifs spécifiques à ce type d’équations [Saa03]. Plus particulièrement, nous proposons de mettre en œuvre un algorithme itératif de type gradients conjugués pour la résolution des équations normales (en anglais Conjugate Gradient for Normal Residual ou CGNR) (voir par exemple [Saa03, Chap. 8]). Cet algorithme est particulièrement adapté à notre problème car on peut exploiter les caractéristiques du dispositif imageur pour une implémentation efficace. 

 Équation normale et algorithme CGNR

 On notera dans la suite L la matrice (HtΓ −1H + µxDt xDx + µyDt yDy + µλDt λDλ) et b le vecteur HtΓ −1m. Sous cette notation la solution de l’équation (4.2) est solution du système linéaire : Lob = b (4.3) Il est alors immédiat de noter que l’équation (4.3) est une équation normale. En effet, la matrice L peut s’écrire : L = HtΓ −1H + µxDt xDx + µyDt yDy + µλDt λDλ = AtA (4.4) avec A la matrice définie par blocs : A = [HtΓ − 1 2 √ µxDt x √µyDt y √ µλDt λ ] t (4.5) Par ailleurs, on peut exprimer le vecteur b en fonction de la matrice A : b = HtΓ −1m = Ateb (4.6) avec le vecteur eb construit par bloc : eb = [Γ − 1 2m 0 0 0]. Ainsi, la solution de l’équation (4.2) recherchée est solution de l’équation normale : AtAob = At˜b (4.7) L’algorithme CGNR (voir Algorithme 1) est un algorithme itératif de type gradient conjugué spécialisé pour la résolution des équations normales telles que l’équation (4.7). D’un point de vue calculatoire, cet algorithme nécessite principalement le calcul de deux produits matriciels faisant intervenir la matrice A, à savoir A · et At · où · dénote un vecteur quelconque de dimension appropriée. Les autres opérations correspondent à des calculs de normes et d’additions vectorielles. Il est important de noter que pour mettre en œuvre cet algorithme, il n’est pas nécessaire de stocker en mémoire la matrice A mais qu’il suffit de pouvoir calculer les produits A · et At ·. On peut donc pour ces calculs, bénéficier des caractéristiques de notre modèle instrumental (matrice H) et de la régularisation quadratique choisie (matrices Dx, Dy et Dλ) pour effectuer ces calculs sans stocker en mémoire ces matrices, ce qui réduit considérablement le besoin en mémoire vive. Algorithme 1 Algorithme de gradients conjugués pour les équations normales (CGNR) x0 : solution initiale q0 = AT (b − Ax0) : p0 = q0 for j = 0, 1 . . . jusqu’à convergence do zj = Apj αj = kqjk 2 2 kzjk 2 2 xj+1 = xj + αjpj qj+1 = qj − αjAT zj βj = kqj+1k 2 2 kqjk 2 2 pj+1 = qj+1 + βjpj end for 

Étude du calcul de A 

 Le calcul de Ao pour un vecteur quelconque o se ramène aux calculs de HtΓ − 1 2 o, Dt xo et Dt yo. Chaque ligne des trois derniers termes consiste en une simple différence de valeurs entre deux éléments de o (différence entre pixels voisins spatialement ou spectralement), donc se calcule sans aucun produit, mais avec respectivement R(C − 1)W, (R − 1)CW et RC(W − 1) différences. Nous avons vu au Chapitre 1, § 1.4 et dans l’article du chapitre 3 que les propriétés du dispositif d’acquisition, en particulier l’indépendance de l’acquisition ligne par ligne et la colocalisation, permettaient d’obtenir un modèle direct simple pour le dispositif. Le calcul de H · correspond exactement au calcul de ce modèle direct. Le schéma de la figure 4.1 illustre la façon dont on peut implémenter ce modèle direct à l’aide d’un cube de filtrage F C (Filtering Cube), correspondant, pour chaque longueur d’onde, à une version décalée de la configuration du DMD : F C(r, c, w) = DMD(r, c−w+w0, w), où r, c, et w sont respectivement les indices de lignes, de colonne et de longueur d’onde (w0 étant une longueur d’onde arbitraire s’imageant au centre du DMD). Ainsi, pour un objet observé O(r, c, w) et une configuration donnée du DMD, on obtient les mesures M(r, c) = P w F C(r, c, w)O(r, c, w). Aussi, le calcul du modèle direct H · peut se faire pour chaque acquisition n ∈ {1, . . . , N} par τRCW additions, où τ est le taux de miroirs ouverts et R et C le nombre de lignes et de colonnes sur le détecteur. Objet O(r, c, w) Cube de filtrage Mesures M(r, c) F C(r, c, w) w r c c r P w r cw r cw • = Figure 4.1 – Illustration de l’effet de l’opérateur H· en une position (r, c) de l’image mesurée. – 81 – 4.1. Implémentation de la méthode RA par gradients conjugués (CGNR) 

 Étude du calcul de At 

Le calcul de Ht · est encore plus simple puisqu’il ne nécessite que des affectations mémoires. Pour comprendre cela, il suffit de remarquer que H et Ht jouent deux rôles opposés, car si H permet de passer d’un cube HS à N images acquises sur le CCD, Ht fait l’inverse, tel qu’illustré sur la figure 4.2 pour une acquisition, où l’image CCD est sur la droite et le cube HS obtenu est sur la gauche. On remarque alors que le contenu du cube HS obtenu par application de Ht à une image CCD, correspond simplement pour chaque position spatiale de l’image, à affecter la valeur du pixel (son intensité) à toutes les bandes spectrales du cube HS correspondant aux éléments non-nuls du filtre hyperspectral. Cette propriété permet de mettre en œuvre l’opération Ht · sans aucun calcul, mais avec de simples affectations en mémoire. Figure 4.2 – Illustration de l’effet de l’opérateur adjoint Ht · en une position (r, c) de l’image mesurée. On peut faire une analyse similaire pour le calcul des termes de régularisation via les produits impliquant les matrices Dx, Dy et Dλ et leurs transposées, qui se calculent uniquement via des différences de valeurs entre pixels voisins. 

Coût de calcul de l’algorithme CGNR 

La charge calculatoire de l’algorithme CGNR est résumée dans le tableau 4.1 qui récapitule des principales opérations exécutées à chaque itération, pour la reconstruction d’un cube HS de taille R×C ×W à partir de N acquisitions de taille R×C pour un taux d’ouverture des miroirs τ et une matrice de covariance Γ diagonale. Complément sur la reconstruction par régularisation quadratique (RA) Étapes Sous-Étapes Opérations + ou – × ou ÷ Affectation zj = Apj Γ − 1 2 Hpj τNWRC NRC Dxpj WRC Dypj WRC Dλpj WRC αj = kqjk 2 2 kzjk 2 2 (3W+N)RC (3W+N)RC xj+1 = xj + αjpj WRC WRC AT zj HtΓ − 1 2 zj NRC τNWRC Dt xzj WRC Dt yzj WRC Dt λ zj WRC qj+1 = qj − αjAT zj WRC WRC βj = kqj+1k 2 2 kqjk 2 2 WRC WRC pj+1 = qj + αjpj WRC WRC Table 4.1 – Nombres d’opérations requises pour une itération du CGNR La charge de calcul totale est donc répartie comme suit (on néglige ici N devant W, étant donné que dans le cadre de cette étude N  W) : — 7WRC produits et additions pour le calcul de αj et βj ainsi que la mise à jour de xj , qj et pj ; — 6WRC additions pour le calcul des termes de régularisation Dx·, Dy· et Dλ· directs et transposés. — 2NRC produits et τNWRC additions pour le calcul de Γ − 1 2 H· et transposé. On peut ainsi constater que, pour une valeur typique de τ = 1 N le coût de calcul du problème direct Γ − 1 2 H· et adjoint HtΓ − 1 2 · est de l’ordre de 2NRC produits et WRC additions et n’est donc pas le principal responsable du coût de calcul de l’algorithme. Une telle répartition de la charge de calcul, due aux propriétés de l’imageur hyperspectral considéré dont nous avons su profiter, est peu commune pour la résolution de problèmes inverses. 

Vitesse de convergence du CGNR

 La vitesse de convergence d’une méthode itérative caractérise la diminution de la distance entre la solution itérative et la solution finale (c’est-à-dire la limite) en fonction du nombre d’itérations. Soit x (n) la solution obtenue à la n ème itération et soit x ∗ la solution optimale, caractériser la vitesse de convergence revient donc à majorer le plus précisément possible le terme kx (n) − x ∗kM où k.kM est une métrique quelconque équivalente à la distance quadratique. Le terme majorant doit dépendre de n et tendre nécessairement vers 0 lorsque n −→ ∞. On identifie donc cette vitesse de convergence vers 0 à celle de la convergence de la méthode itérative. Pour la famille des méthodes basées sur l’espace de Krylov, dont fait partie la méthode CGNR, il existe une majoration générale vérifiée par celles-ci, dont l’expression est [DCMR95, page 7] : kx (n) − x ∗ kM ≤ 2 p cond{L} − 1 p cond{L} + 1!n kr (0)k2 (4.8) – 83 – 4.2. Vitesse de convergence du CGNR où r (0) = b − Lx (0) est le résidu initial et cond{L} le conditionnement de la matrice L, c’està-dire le rapport entre sa plus grande et sa plus petite valeur propre. L’équation (4.8) montre l’impact du conditionnement de L sur la vitesse de convergence de l’algorithme et la nécessité de garder cette valeur le plus proche possible de 1. Un des moyens les plus utilisés pour assurer une telle condition sur cond{L}, est de modifier L en la multipliant par une matrice P choisie judicieusement de sorte que le conditionnement cond{PL} de la matrice PL soit plus petit que cond{L} et se rapproche donc de 1. Cette technique connue sous le nom de préconditionnement [BSvD99, DCMR95], permet de modifier l’équation (4.3) à résoudre en une équation équivalente, dont la résolution nécessite moins d’itérations et donc potentiellement moins de calcul. Néanmoins, le coût de calcul du pré-conditionneur s’ajoute au coût de calcul de l’algorithme. Remédier au problème de conditionnement passe donc par un compromis entre le coût calculatoire du pré-conditionneur et le gain en nombre d’itérations nécessaires à la convergence de l’algorithme. Dans le cadre de cette thèse, on peut jouer sur un autre paramètre pour améliorer le conditionnement de la matrice L. En effet, le dispositif imageur étant pilotable, par l’intermédiaire de sa matrice de micro-miroirs (DMD), il est possible de modifier la matrice H intervenant dans la matrice L et viser à améliorer le conditionnement de cette matrice sans avoir à faire appel à un préconditionneur. . . Il est donc important d’étudier ce conditionnement d’envisager des configurations du DMDpermettant de l’améliorer. La matrice HtΓ −1H correspond à la contribution de l’imageur que l’on peut modifier en pilotant le filtrage effectué par le DMD, afin d’améliorer le conditionnement de la matrice L et donc la vitesse de convergence de l’algorithme CGNR. L’idée est donc de déterminer les matrices HtΓ −1H qui soient réalisables par l’imageur hyperspectral tout en permettant de minimiser le conditionnement cond{L}. Pour contrôler le conditionnement cond{L} à l’aide du DMD, nous proposons dans un premier temps de le majorer par une fonction de cond{HtΓ −1H}. Ensuite, nous verrons comme certaines configurations particulières du DMD permettent d’exprimer aisément cond{HtΓ −1H} et de le réduire à 1 dans le cas d’un modèle de bruit blanc. Dans la suite, nous travaillerons uniquement avec des matrices symétriques à valeurs réelles, dont la décomposition en valeurs propres et vecteur propres est définie, avec des valeurs propres positives ou nulles (les matrices considérées s’écrivant sous la forme MtM elles sont définies nonnégatives). Pour une matrice symétrique M on notera eMn ses vecteurs propres de norme unité associés aux valeurs propres λMn . On notera par ailleurs, par eMmax et eM min les vecteurs propres associés aux valeurs propres maximale λMmax et minimale non nulle λM min. Ainsi, par définition, cond{M} = λMmax λM min . Afin d’alléger les notations, nous introduisons les matrices J et R telles que L = J + R avec J = HtΓ −1H et R = µxDt xDx + µyDt yDy + µλDt λDλ .

 Majoration de cond{L}

 Théorème 1 : D’après les notations précédentes, si l’on défini τ R,J max = λR max λJmax , on a : cond{L} ≤ (1 + τ R,J max) Prg(J) n=1 < e L min, e J n >2 cond{J} (4.9) Démonstration : Soit λ L ` une valeur propre quelconque de L et e L ` son vecteur propre associé de norme unité, on a : (e L ` ) tLe L ` = λ L ` (e L ` ) te L ` = λ L ` . (4.10) – 84 – Chapitre 4. Complément sur la reconstruction par régularisation quadratique (RA) En remplaçant L par J + R dans cette expression et en écrivant la décomposition en valeurs et vecteurs propres de R et L sous la forme J = rg X (J) n=1 λ J ne J n (e J n ) t et R = rg X (R) k=1 λ R k e R k (e R k ) t (4.11) on obtient : λ L ` = (e L ` ) t (J + R)e L ` = rg X (J) n=1 λ J n < e L ` , e J n > 2 + rg X (R) k=1 λ R k < e L ` , e R k > 2 . (4.12) Dans un premier temps on peut montrer que dans notre cas toutes les valeurs propres λ L n sont strictement positives, donc que la matrice L est inversible, hormis dans un cas extrêmement particulier mathématiquement. En effet, l’équation (4.12) nous dit que toute valeur propre s’écrit comme la somme de termes positifs ou nuls. Pour qu’une valeur propre soit nulle, il faut donc que tous les termes soient nuls. Les sommes se faisant uniquement sur les rg(J) et rg(R) valeurs propres non nulles des matrices J et R, cela signifie donc que tous les produits scalaires doivent être nuls. Or, par construction, si les paramètres de régularisation µx, µy et µλ sont non nuls la matrice R est de rang RCW − 1 et le seul vecteur propre e R RCW associé à la valeur propre nulle est le vecteur constant. En effet u tRu = 0 est nul uniquement pour les vecteurs u proportionnels à e R RCW puisque cette quantité correspond à la somme (pondérée) du carré des différences interpixels le long des lignes, des colonnes et des longueurs d’ondes. Aussi, les termes < e L ` , e R k >2 ne peuvent être tous nuls (pour k ∈ {1, . . . , RCW − 1}) que si e L ` = e R RCW . Pour que ce vecteur constant soit vecteur propre de la matrice L il faudrait, d’après la définition (4.5) de la matrice A que Ae R RCW = 0 et donc par construction, que He R RCW = 0. Autrement dit, que pour les N configurations du DMD des différentes acquisitions, les mesures soient insensibles à l’ajout d’une grandeur constante à l’objet observé. Cela ne sera donc bien évidemment jamais le cas en pratique.