Accueil CV Thèse Théâtre Photos
1. Introduction

  • Introduction[RB1]
  • Les méthodes de reconstruction standard qui sont utilisées en routine ne passent pas par la résolution d’un système tel que nous l’avons posé au chapitre précédent. La reconstruction s’effectue par une rétroprojection filtrée.

    A coté de ces méthodes dites standards, les techniques de reconstruction basées sur une modélisation statistique des données connaissent un essor considérable, notamment avec l’emploi des champs de Markov qui permettent de glisser de l’information à priori lors de la reconstruction. En raison des nombreux paramètres spécifiques à chaque image reconstruite, elles restent mal commodes pour un emploi standard. En effet, pour améliorer la qualité de reconstruction, il est nécessaire d’ajuster ces paramètres au cas par cas.

    Une autre catégorie de reconstruction regroupe les méthodes dites algébriques qui elles, passent par la résolution du système p=Hc. Au sein de ce groupe, on distingue les méthodes en fonction du type de résolution : itérative ou directe. On va étudier plus en détail des méthodes relatives à cette dernière catégorie.

    Les méthodes directes ou matricielles de reconstruction d’image n’ont pas été utilisées jusqu’à présent d’une manière générale, du fait de la taille importante du système considéré, du comportement des matrices lors de l’inversion et du succès des résolution utilisant une rétroprojection filtrée. Les progrès dans les capacités de calculs et de mémoire des ordinateurs rendent une telle résolution envisageable.

    Dans le principe, le système que l’on cherche à résoudre se résume à :

    p = Hc

    et l’on souhaite trouver une matrice K telle que :

    c = Kp

    Des problèmes se posent assez rapidement puisque la matrice H est en général, rectangulaire et singulière d’une part, et que d’autre part notre système est bruité.

    Il existe un ensemble de méthodes puissantes pour résoudre les systèmes d’équations qui correspondent à des matrices H singulières où numériquement proches d’être singulières. Parmi celles ci, la décomposition en valeurs singulières (SVD) reste la plus générale, mais nous verrons que dans certains cas, le problème tel que nous l’avons posé peut être résolu en utilisant une factorisation QR Householder. En tous les cas, les méthodes de résolution directe impliquent bien souvent une factorisation de la matrice H.

  • Factorisation et Résolution
    1. Décomposition en Valeurs Singulières[RB14] [RB15][RB18]
      1. Définition
      2. La décomposition en valeurs singulières est basée sur le théorème suivant :

        Toute matrice H de dimension M*N pour laquelle le nombre de lignes M est supérieur ou égal au nombre de colonnes N peut s’écrire sous la forme du produit d’une matrice V de dimension M*N, d'une matrice diagonale dont les éléments sont positifs ou nuls et de la transposée d'une matrice U de dimension N*N

        où les matrices U et V sont orthonormées :

        Il est toujours possible de faire cette décomposition quel que soit le "degré de singularité" de la matrice et celle ci est unique.

         

      3. Observations
      4. La matrice H est rectangulaire de dimension M*N, elle est donc au plus de rang N. Supposons par la suite que H est de rang R<N.

        HTH est une matrice de dimension N*N dont les valeurs propres l i sont telles que :

        et dont les N vecteurs propres associés sont

        HHT est une matrice de dimension M*M dont les valeurs propres l i sont telles que :

        et dont les M vecteurs propres associés sont

        Remarquons au passage que les valeurs propres l i sont les mêmes pour HTH et HHT et que :

        De plus, on montre aisément que

        Les vecteurs propres de HTH , ie les vecteurs {ui} forment une base complète de décomposition pour tout vecteur de dimension N. Ainsi le vecteur c peut se décomposer sur cette base :

        Au chapitre précédent nous avons montré que l'information relative à l'action physique de notre système était contenue dans la matrice H. En appliquant cette transformation linéaire à notre vecteur c on obtient :

        Comme la matrice H est supposée être de rang R, les valeurs propres de la matrice sont nulles pour tout i>R, et p se résume à :

        Si on réécrit la décomposition de c dans la base des {ui} en la scindant en deux termes tels que :

        On s'aperçoit que le premier terme cIm spécifie les composantes de c qui peuvent être vues par le système. Le deuxième terme est transparent à H, ie il n'apporte pas de contribution à p.

        Si on applique H à l’un des vecteurs propres ui, on retrouve alors, à cette constante multiplicative l i près, le vecteur propre vi. Autrement dit, l'image par H d'un vecteur du sous-espace de dimension N engendré par les vecteurs {ui} est un vecteur du sous-espace de dimension R engendré par les R premiers vecteurs {vi}.

        Dans notre système p=Hc, si on considère H comme un mapping linéaire de l'espace W vers Y , et si de plus H est singulière, alors il existe un sous-espace W 2 de W dont l'image est zéro dans Y . Il s'agit du noyau de l'application H, noté Ker(H). D'autre part, il existe un sous espace Y 1 de Y qui peut être atteint par H. Ce sous espace est l'image de l'application H, noté Im(H)

        Qu'est-ce que cela vient faire avec la SVD ? La SVD construit explicitement une base orthonormale pour le noyau et l'image de H. Autrement dit, les colonnes de V qui correspondent à des indices i de valeurs propres l i¹ 0 définissent l'image de H et les colonnes de U qui correspondent à des indices i de valeurs propres l i=0 définissent le noyau de H.

      5. Résolution
      6. Dans notre système p=Hc, le jeu de mesures n’a aucune raison de se trouver dans l’espace Im(H). Il n’existe donc pas une unique solution à ce système. En revanche même si p n’appartient pas à Im(H) on peut séparer la décomposition de ce vecteur dans la base des {vi} en deux termes :

        l’un pIm compris dans Im(H), l’autre non.

        On cherche alors tel que pIm=HcIm .

        On obtient ainsi une solution particulière au système p=Hc. Ensuite toutes les solutions s’obtiennent par c = cim + ckercker est un élément de W 2. Pour trouver cIm on ne considère que le système réduit suivant :

        VR (respectivement UR) de dimension M*R (respectivement N*R) sont obtenues à partir de V (respectivement U) en ne conservant que les R premières colonnes et est une sous matrice diagonale de dimension R*R extraite de .

        et la solution est alors :

        Même si une solution unique n’existe pas, restons positifs : la SVD permet de construire une solution la plus proche possible de la solution exacte au sens des moindres carrés.

        Moindres Carrés, le mot est lâché...

      7. Solution au sens des Moindres Carrés

      Nous cherchons à minimiser (Chapitre précédent) l’écart entre le modèle de projection et les mesures effectuées. La solution c choisie sera la solution la plus proche de l’activité réelle au sens des moindres carrés pour ce modèle. C’est à dire c sera tel que soit minimum.

      Remarque :

      On envisage dans notre cas uniquement la norme euclidienne.

      Or

      et

      d’où

       

      Soit une colonne de H, les N colonnes induisent un espace Y 1 (=Im(H)).

      Minimiser

      revient à déterminer le vecteur de Y 1 qui approche le mieux p, C’est à dire la projection orthogonale de p sur Y 1. Or le vecteur c solution s’écrit :

      c=(HTH)-1HTp

      ou encore :

      Hc=H(HTH)-1HTp

      et H(HTH)-1HT est connue pour être une matrice de projection orthogonale.

      On cherche donc une solution de telle sorte que Hc soit une projection orthogonale sur Im(H). Celle-ci s’obtient en inversant HTH, matrice carrée de plus petite taille que H.

      La SVD, quant à elle, calcule explicitement une base de Y . En scindant en deux termes la décomposition de P dans cette base, on peut extraire directement la projection sur Y 1 en ne conservant que le premier terme.

      On désigne donc par Moindres Carrés une solution c qui minimise . En revanche, si on trouve pour H une factorisation telle que H=QRQ est orthonormée, alors :

      Il se peut alors que le système à résoudre soit plus simple. Voyons donc plus en détail une telle factorisation.

    2. Factorisation QR et Résolution[RB5] [RB16]
    3. Il existe un certain nombre de transformations orthogonales élémentaires sur lesquelles la factorisation QR va s’appuyer. Pour réaliser une factorisation, il est nécessaire d’introduire des changements de base qui annulent certains coefficients. Les transformations de Givens annulent certains coefficients par le biais de rotations, les transformations de Householder par le biais de réflexions. On va se limiter à l’emploi de ces dernières

      1. Vecteur et Transformation de Housholder
        1. Définition et Propriétés
        2. On appelle matrice élémentaire de Householder une matrice Ho de la forme Ho= Id - 2uuTu est un vecteur de dimension M de norme égale à 1. Ce vecteur est le vecteur de Householder.

          Cette matrice Ho est symétrique et orthogonale.

        3. Interprétation géométrique
        4. Ho représente la symétrie par rapport au plan P passant par l’origine et perpendiculaire à la direction u. Appliquons Ho à un vecteur a :

          Hoa = (Id - 2uuT)a = a - 2 uuTa = a - 2(uTa)u

          Géométriquement le point M’ symétrique du point M par rapport au plan P est défini par :

           

          Si OM = a, comme u unitaire est représenté par

          Donc :

           

        5. Intérêt Pratique
          1. Intérêt
          2. Tout l’intérêt d’une telle transformation réside dans la remarque suivante :

            Etant donnés deux vecteurs a et c de dimension M, non colinéaires, on peut déterminer un plan P ( Il est possible de trouver un vecteur u de dimension M , u unitaire ) pour que les directions des deux vecteurs soient symétriques par rapport à P ( tel que si Ho= Id - 2uuT alors Hoa = c )

            Corollaire :

            Si on se fixe un vecteur b, unitaire on peut trouver un u et un a tel que :

            Hoa = a b

          3. Formulation

    Si on pose v = 2l u et b =2l 2 alors :

    Dans la pratique on utilise cette transformation élémentaire pour b=e1. Dans ce cas, à une nuance près sur le signe de a qui est fonction du signe de a(1), c’est sous cette forme qu’est implémentée l’étape élémentaire de la factorisation QR Householder.

      1. Pour un vecteur a donné on cherche le vecteur de Householder v associé ( Procédure House ) tel que v(1) = 1.

      1. Ensuite pour un vecteur d connu, on peut appliquer la transformation de Householder Ho très simplement :

      • Calcul de
      • Calcul de Ho d = d - ( vTd )v

    Plus généralement on applique la transformation élémentaire de Householder basé sur le vecteur v à une matrice A dans laquelle on souhaite créer des zéros sur la première colonne. (Procedure RowHouse). Voir Annexe pour plus de détails concernant l’implémentation de ces procédures.

          1. Factorisation QR Householder
          2. Soit une matrice A de dimension M*N avec M>N

            Remarque :

            On ne prend pas H pour éviter certaines confusions avec la matrice de Householder.

            De même, le système p=Hc va être remplacé par b=Ax

            On suppose dans un premier temps que la matrice A est de rang N.

            1. Définition
            2. La factorisation QR Householder va trouver N matrices de Householder élémentaires H1,....,HN. Si on forme Q = H1H2....HN on obtient par QTA une matrice R, triangulaire supérieure et non singulière.

            3. Principe
            4. Décrivons la première étape . Soit a1 de dimension M la première colonne d’une matrice A0=A. Nous avons vu comment déterminer H1 = I - 1/b ( v1 v1T ) pour que H1a1=a 1e1. La matrice A1= H1A0 aura sa première colonne égale à a 1e1, ie toutes ses composantes sont nulles sauf la première. La ième colonne de A1 est égale à H1ai ai est la ième colonne de A0.

              Pour continuer, il suffit d’appliquer une transformation analogue à la sous- matrice formée des (M-1) dernières lignes et des (N-1) dernières colonnes de A1. Et ainsi de suite jusqu’à An qui est rectangulaire supérieure. On a ainsi la récurrence :

              A0= A

              Ai= HiAi -1

              Il est évident que pour l’implémentation il est inutile de stocker l’ensemble des matrices Hi, on a juste besoin de pouvoir calculer la transformation de Householder sur un vecteur et donc le vecteur vi est suffisant pour cela. Dans notre cas comme vi(i) =1, nous avons :

              avec :

              comme , Q est orthogonale puisque les Hi le sont.

            5. Propriétés

    Soient ai et qi les ièmes colonnes de A et Q, alors :

      • Comme R est triangulaire supérieure rij=0 pour j>i. on peut donc écrire

    On vérifie par récurrence que le vecteur qk est une combinaison linéaire des vecteurs a1,.....,ak .

      • Comme la matrice Q est orthogonale on a :

            1. Cas ou rang(A)<N
              1. Effet d’un rang < N
              2. Nous avons introduit au paragraphe 1.2.2.2.3, une récurrence. Cette dernière montrait qu’à l’itération i, on construisait le vecteur qi en considérant la colonne ai de notre matrice A et en calculant la différence:

                Qu’est ce que cela signifie géométriquement ?

                Après i-1 itérations, Les vecteurs q1,.....,qi-1 calculés(ie les i-1 premières colonnes de A) engendrent un espace Ei-1 . On regarde alors comment se situe le vecteur ai par rapport à cet espace.

                Dans le cas où la matrice A est de rang N, à chaque itération i le vecteur ai est linéairement indépendant des vecteurs précédemment calculés. Il n’appartient pas au sous-espace Ei-1 et donc le vecteur qi calculé à partir de cette différence apporte une nouvelle dimension.

                Dans le cas où le rang de la matrice A est inférieur à N, cela signifie que les colonnes de A ne sont plus linéairement indépendantes. Ainsi il existe forcément des indices w pour lesquels les vecteurs aw s’expriment comme combinaison linéaire des vecteurs a1,....,aw-1. En conséquence le vecteur aw appartient déjà à Ei-1 et n’apporte pas de dimension supplémentaire (qi serait nul).

                Pratiquement cela signifie, si k est le premier de ces indices w, qu’ayant effectué les (k-1) premières transformations de Householder, on n'a pas besoin de réduire la keme colonne qui est nulle. On la laisse et on poursuit en passant aux colonnes suivantes.

              3. Permutation des colonnes

    Dans notre récurrence du paragraphe 1.2.2.2.3, on construit progressivement une base {q}. Chaque élément de la base est construit en prenant successivement chacune des colonnes de A, et finalement si le rang de A est R, les R vecteurs {q} engendre un espace de dimension R. Mais il se peut que parcourir successivement les colonnes de A ne soit pas la meilleure façon de générer une base de cet espace de dimension R.

    On peut donc envisager de parcourir A de manière différente en permutant l’ordre des colonnes. Pour ce faire, on va choisir dans les vecteurs colonnes qui restent à orthogonaliser celui qui maximise di (voir figure ci dessus). Ainsi à l’itération i de notre factorisation QR , dans la matrice Ai restant à factoriser, on calcule la norme des différentes colonnes de Ai. On permute alors celle d’indice l qui correspond à la norme maximale avec la première colonne de Ai , on obtient alors AiPiPi représente la permutation à l’itération i. C'est sur cette nouvelle matrice que l’on applique la transformation élémentaire de Householder.

    Finalement on obtient la factorisation AP = QR P = P1....PR.

          1. Moindres Carrés et factorisation QR[RB2] [RB6]

    En quoi une aussi jolie décomposition peut-elle se révéler utile dans un problème des moindres carrés?...

    A l’issue de notre factorisation nous avons :

    où R est le rang de A, Q est une matrice orthogonale, R11 est triangulaire supérieure et non singulière et P représente l’ensemble des permutations.

    Partant de cela, on peut directement retranscrire le problème des moindres carrés :

    où :

    et :

    ainsi x est obtenu au sens des moindres carrés par :

    Si z est mis à zéro, on trouve xb :

    On n’obtient pas une solution de norme minimale sauf si R12 est nulle en revanche la solution obtenue xb est telle que :

    L’implémentation est maintenant quasi directe. Pour l’implémentation on se reporte à l’Annexe.

     

     

        1. Stabilité et régularisation
          1. Pour la factorisation SVD [RB19] [RB18]

    Nous avons vu que pour une résolution de notre système p=Hc, nous étions amenés à résoudre le système tronqué :

    Que se passe t-il si l’on tronque plus sévèrement avec Q<R et que l’on résout le système ainsi posé ?

    Pour cela, il est nécessaire de se rappeler que les valeurs singulières sont rangées par ordre décroissant dans notre matrice L . Or, du fait de la construction de H, on constate que plus l’indice de la valeur propre augmente, plus la sous image qui lui est associée contient de hautes fréquences. Cela revient à dire que notre troncature à Q va faire disparaître des hautes fréquences de notre image et donc des détails, mais que l’image globale sera conservée.

    Si maintenant, on regarde l’évolution de la distance de l’image reconstruite à l’image originale en fonction du nombre de valeurs singulières retenues, on obtient :

      • la courbe en trait plein pour des projections non bruitées
      • la courbe en pointillé si les projections sont bruitées.

    Pour une image non bruitée, il faut conserver toutes les valeurs singulières non nulles pour retrouver l’image originale. Les valeurs de haut indice permettent de retrouver les détails de l’image. En revanche, en présence de bruit, les valeurs singulières de plus haut indice perturbent la reconstruction. On ne retiendra donc que les Q premières valeurs, juste avant que la courbe en pointillé ne s’infléchisse, pour obtenir la distance à l’image originale minimale.

    Sur une image bruitée donc, lors de la décomposition SVD, de très petites valeurs singulières vont apparaître, qui, en l’absence de bruit auraient dû être exactement zéro. Si ces valeurs singulières sont conservées lors de l’inversion du système, cette erreur même petite va induire une forte erreur après inversion, ce qui aura un effet désastreux sur la représentation de l’image. En tronquant, on va donc limiter cet effet et par la même effectuer une certaine régularisation.

    D’une manière générale, la régularisation de notre système peut être réalisée en pondérant les différentes valeurs propres par une fonction L.

          1. Pour la factorisation QR

    D’une manière analogue, nous avons vu que la factorisation QR Householder avec permutation des colonnes revenait à construire à l’aide des colonnes de H une base et que cette permutation répondait à une certaine stratégie. De fait, on obtenait également un classement des différentes colonnes dans un ordre décroissant " d’importance ".

    Soit J un sous-ensemble d’indices contenus dans {1,....,N}. Appelons AJ la sous matrice extraite de A, en ne conservant que les colonnes aj pour j appartenant à J. Le vecteur, noté xJ , est obtenu, après inversion, en ne conservant que ces card(J) colonnes. On montre que le résidu qui traduit la différence entre la reconstruction et l’image réelle va diminuer à mesure que le nombre de colonnes J retenues augmente.

    Par conséquent, à l’issue de ces deux remarques il s’avère également possible lors de cette factorisation et en présence de bruit, d’effectuer une séléction des colonnes les plus significatives donc de régulariser notre solution.