16 Jan 2019

Wiki

Python

Aide

edit SideBar

Search

Decomposition Matricielle


Dans ce qui suit, on présente une bibliothèque de calculs formels sous python, ainsi que différentes techniques de décompositions matricielles.

Ce TP s'inspire du livre Maths et Maple de Jean-Michel Ferrard : il condense le chapitre 5, l'applique à python, et propose des exercices en plus.

1.  sympy

1.1  La bibliothèque sympy

Présentation

Sympy est une bibliothèque de calcul formel, assez jeune mais en fort développement. A la différence de numpy, le calcul est ici exact (et non pas approché), et les variables peuvent être utilisées dans le calcul. On y perd une certaine rapidité, et moins de possibilités.

Il me semble intéressant de connaître les diverses possibilités qui vous sont offertes pour résoudre un problème, et de savoir passer d'une bibliothèque à une autre.

Si ces dernières sont propres à python, leurs équivalents existent pour d'autres langages, et certains logiciels ne se vouent qu'au calcul formel (tel est le cas de Maple, Mathematica, Octave, etc.) ou numérique (matlab, scilab, etc.). Si la syntaxe change, la philosophie reste la même.

Installation

Cas de la salle TP.

Normalement, la bibliothèque sympy est installée en salle machine. Si tel n'est pas le cas, essayez de l'installer dans un environnement virtuel de travail.

  • Si vous avez déjà votre environnement de travail, lancez-le.
  • Sinon, créer un environnement (à faire dans un terminal, en remplaçant monenv par le nom de votre choix) :
  guyeux@kerry:~$ virtualenv monenv

Cet environnement est à créer une seule fois, vous n'aurez plus à le refaire.

A chaque fois, par contre, il faudra l'activer :

  guyeux@kerry:~$ source monenv/bin/activate
en remplaçant, dans ce qui précède, monenv par le nom que vous avez choisi pour votre environnement virtuel de travail.

Vous êtes maintenant dans votre environnement virtuel, et vous pouvez installer les modules Python de votre choix, tel sympy, grâce à la commande easy_install...

  (mon_env)guyeux@kerry:~$ easy_install sympy

Si cela n'aboutit pas, appeler votre professeur à l'aide.

Chez vous.

On peut installer sympy ainsi sous une debian GNU/Linux, ou sous une Ubuntu :

  $ aptitude install python-sympy

dans un shell, en étant administrateur.

1.2  Création d'une matrice

Une fois python lancé, voici comment importer la bibliothèque sympy, et créer une première matrice :

  >>> from sympy import *
  >>> M=Matrix([ [1,2],[3,4] ])
  >>> M
  [1, 2]
  [3, 4]

Pour accéder à un élément, on procède comme dans numpy :

  >>> M[0,0]
  1

(Attention : l'élément supérieur gauche est en position (0,0)).

Dans ce qui suit, on va reprendre chaque point évoqué dans la section numpy, et montrer l'équivalent dans sympy. A charge pour vous de choisir la bibliothèque qui vous convient le mieux, selon vos goûts et vos besoins.

1.3  Matrices particulières

Venons-en à la création de matrices particulières. La matrice nulle se définit comme dans numpy :

  >>> zeros((3,4))
  [0, 0, 0, 0]
  [0, 0, 0, 0]
  [0, 0, 0, 0]

Les matrices

  • identité,
  • et constituée de uns,

sont ici forcément carrées :

  >>> eye(3)
  [1, 0, 0]
  [0, 1, 0]
  [0, 0, 1]

  >>> ones(3)
  [1, 1, 1]
  [1, 1, 1]
  [1, 1, 1]

Enfin, on peut créer une matrice aléatoire comme suit :

  >>> randMatrix(3,4,min=0, max=5)
  [3, 4, 4, 5]
  [3, 5, 0, 5]
  [1, 4, 5, 2]

Ici, la matrice générée est constituée d'entiers. Il faudra la diviser terme à terme pour obtenir des réels. Les arguments min et max sont optionnels, et fixés respectivement à 0 et 100 par défaut. Notons, pour finir, qu'il est possible de fournir la graine (argument seed) de génération des nombres pseudo-aléatoires.

1.4  Opérations matricielles

Considérons les deux matrices suivantes :

  >>> A = Matrix([[1,2,3],[4,5,6],[7,8,9]])
  >>> B = Matrix([[1,2,1],[2,1,2],[1,1,1]])

On a accès aux opérations matricielles élémentaires : somme et multiplication matricielle, multiplication par un scalaire...

  >>> A+B
  [2, 4,  4]
  [6, 6,  8]
  [8, 9, 10]
  >>> A*B
  [ 8,  7,  8]
  [20, 19, 20]
  [32, 31, 32]
  >>> 0.5*A
  [0.5,   1, 1.5]
  [2.0, 2.5, 3.0]
  [3.5, 4.0, 4.5]
  >>> A/2
  [1/2,   1, 3/2]
  [  2, 5/2,   3]
  [7/2,   4, 9/2]

On a encore la possibilité de transposer les matrices, ou d'en calculer leur trace :

  >>> A
  [1, 2, 3]
  [4, 5, 6]
  [7, 8, 9]
  >>> A.T
  [1, 4, 7]
  [2, 5, 8]
  [3, 6, 9]
  >>> A.trace()
  15

1.5  Les fonctions importantes

On peut directement calculer une inverse de matrice, un déterminant :

  >>> A
  [1, 2, 3]
  [4, 5, 6]
  [7, 8, 1]
  >>> A.det()
  24
  >>> A.inv()
  [-43/24, 11/12, -1/8]
  [ 19/12,  -5/6,  1/4]
  [  -1/8,   1/4, -1/8]

La différence entre calcul numérique et formel apparaît : on n'a pas ici une valeur approchée de $A^{-1}$, mais une valeur exacte. On peut aller un peu plus loin :

  >>> A*A.inv()
  [1, 0, 0]
  [0, 1, 0]
  [0, 0, 1]
  >>> A*A.inv()==eye(3)
  True

1.6  Systèmes linéaires

Voici comment procéder pour résoudre le système $\left\{ \begin{array}{rll} x + 4 y & = & 2 \\ -2 x + y & = & 14\end{array}\right.$

On introduit les variables $x$ et $y$, après avoir importé la bibliothèque sympy :

  >>> from sympy import *
  >>> x, y = symbols('xy')

On pose le système, puis on le résout :

  >>> system = Matrix(( (1, 4, 2), (-2, 1, 14)))
  >>> solve_linear_system(system, x, y)
  {x: -6, y: 2}

2.  Décomposition LU

2.1  Présentation

Soit $A$ une matrice carrée d'ordre $n$ et inversible. On appelle décomposition LU de $A$ toute égalité $A=LU$, où

  • $L$ est triangulaire inférieure avec des coefficients diagonaux égaux à 1,
  • $U$ est triangulaire supérieure.

Obtenir une décomposition LU d'une matrice $A$ est important lorsqu'on souhaite résoudre plusieurs fois à la suite des systèmes linéaires du type $Ax = b$. Il suffit alors en effet de résoudre deux systèmes triangulaires. En effet, le système devient

$LUx = b \Leftrightarrow \left\{\begin{array}{cc} Ly = b& (1),\\ Ux = y &(2). \end{array}\right.$

On résout le système (1) pour trouver le vecteur y, puis le système (2) pour trouver le vecteur x. La résolution est facilitée par la forme triangulaire des matrices.

Si donc on a a résoudre $Ax = b$, avec 10000 vecteurs $b$ différents, mais toujours la même matrice $A$, alors

  • on fait, une fois pour toute, la décomposition LU de $A$,
  • on résout les 20000 systèmes triangulaires, ce qui se fait en un temps "raisonnable".

La décomposition LU permet encore de calculer le déterminant de $A$, qui est égal au produit des éléments diagonaux de la matrice $U$ (si $A$ admet bien une telle décomposition).

2.2  Détermination avec sympy

Soit le bout de code suivant :

  >>> from sympy import *
  >>> M = Matrix([[1,4,3],[1,-2,1],[0,1,3]])
  >>> M
  [1,  4, 3]
  [1, -2, 1]
  [0,  1, 3]
  >>> (L,U,row) = M.LUdecomposition()
  >>> L
  [1,    0, 0]
  [1,    1, 0]
  [0, -1/6, 1]
  >>> U
  [1,  4,   3]
  [0, -6,  -2]
  [0,  0, 8/3]
  >>> L*U
  [1,  4, 3]
  [1, -2, 1]
  [0,  1, 3]

On voit que L et U ont bien la forme souhaitée, que leur produit donne bien M.

2.3  Existence et unicité

Résultat

On peut prouver que :

  • Cette décomposition, quand elle existe, est unique : s'il existe un couple de matrices $(L,U)$ vérifiant les propriétés ci-dessus, alors il n'y en a pas d'autres.
  • Cette décomposition existe si et seulement si les mineurs principaux de $A$ sont non nuls.

On rappelle, si besoin est, que le mineur principal $A_k$ d'ordre $k$ de $A$ est le déterminant obtenu à partir des $k$ premières lignes et des $k$ premières colonnes de $A$. Par exemple, pour

$A = \left(\begin{array}{ccc}1 & 4 & 3\\1 & -2 & 1\\0 & 1 & 3\end{array}\right)$

on trouve les mineurs principaux suivants :

$A_1 = |1|, A_2 = \left|\begin{array}{cc}1 & 4 \\1 & -2 \end{array}\right|, A_3 = det(A)$

Avec sympy

Sympy ne permet pas d'obtenir les mineurs principaux. Il ne qu'extraire les sous-matrices obtenues en enlevant une ligne et une colonne, ce qu'il appelle des mineurs (attention à la confusion des termes).

Cela dit, ces «mineurs» peuvent être utilisés pour récupérer les mineurs principaux :

  >>> M
  [1,  4, 3]
  [1, -2, 1]
  [0,  1, 3]
  >>> M.minorMatrix(2,2)
  [1,  4]
  [1, -2]
  >>> M.minorMatrix(2,2).minorMatrix(1,1)
  [1]

Il reste alors à calculer les déterminants de ce qui précède.

Travaux pratiques

Faire une fonction qui teste si une matrice donnée en argument possède une décomposition LU, ou pas.

On rappelle qu'il faut procéder de la façon suivante pour définir une fonction :

  >>> def f(x):
  ...     return 2*x+1
  ... 
  >>> 

2.4  Réalisation de la décomposition

Le but, maintenant, n'est plus d'utiliser des fonctions préexistantes, mais d'apprendre à faire soi-même une décomposition LU. On commence par utiliser la puissance d'une bibliothèque de calcul formel (sympy).

Utilisation du calcul formel

On montre ici comment faire une décomposition LU à la main, en utilisant à plein le calcul formel.

Après avoir importé la bibliothèque sympy, on crée une matrice

$ A = \left( \begin{array}{ccc} -85 & -55 & -37 \\ -35 & 97 & 50 \\ 79 & 56 & 49 \end{array}\right)$

comme suit :

  >>> from sympy import *
  >>> A = Matrix([[-85, -55, -37],[-35,97,50],[79,56,49]])

On vérifie que $A$ possède bien une décomposition LU :

  >>> A.det()
  -121529
  >>> A.minorMatrix(2,2).det()
  -10170
  >>> A.minorMatrix(2,2).minorMatrix(1,1).det()
  -85

Les mineurs principaux étant tous non nuls, $A$ possède bien une décomposition LU. Définissons formellement L et U, et créons leur produit :

  >>> L = eye(3)
  >>> for j in range(2):
  ...     for i in range(j+1,3):
  ...         L[i,j] = Symbol('L'+str(i)+str(j))
  ... 
  >>> L
  [  1,   0, 0]
  [L10,   1, 0]
  [L20, L21, 1]
  >>> U = zeros(3)
  >>> for i in range(3):
  ...     for j in range(i,3):
  ...         U[i,j] = Symbol('U'+str(i)+str(j))
  ... 
  >>> U
  [U00, U01, U02]
  [  0, U11, U12]
  [  0,   0, U22]
  >>> B=L*U
  >>> B
  [    U00,               U01,                     U02]
  [L10*U00,     U11 + L10*U01,           U12 + L10*U02]
  [L20*U00, L20*U01 + L21*U11, U22 + L20*U02 + L21*U12]

Ce produit B = LU devant être égal à A, on peut procéder par identification terme à terme des matrices A et B, ce qui donne un système qui peut être résolu :

  >>> X = solve(((A-B)[0,0],(A-B)[0,1],(A-B)[0,2],
  ...              (A-B)[1,0],(A-B)[1,1],(A-B)[1,2],
  ...              (A-B)[2,0],(A-B)[2,1],(A-B)[2,2]),
  ...              'L10','L20','L21',
  ...              'U00','U01','U02','U11','U12','U22')
  >>> X
  [(7/17, -79/85, 83/2034, -85, -55, -37, 2034/17, 1109/17, 121529/10170)]

Substituons les valeurs aux variables dans L et U :

  >>> X = X[0]
  >>> L=L.subs('L10',X[0])
  >>> L=L.subs('L20',X[1])
  >>> L=L.subs('L21',X[2])
  >>> U=U.subs('U00',X[3])
  >>> U=U.subs('U01',X[4])
  >>> U=U.subs('U02',X[5])
  >>> U=U.subs('U11',X[6])
  >>> U=U.subs('U12',X[7])
  >>> U=U.subs('U22',X[8])
  >>> L
  [     1,       0, 0]
  [  7/17,       1, 0]
  [-79/85, 83/2034, 1]
  >>> U
  [-85,     -55,          -37]
  [  0, 2034/17,      1109/17]
  [  0,       0, 121529/10170]

On a évidemment les formes attendues. On peut s'assurer que LU redonne bien A :

  >>> L*U
  [-85, -55, -37]
  [-35,  97,  50]
  [ 79,  56,  49]
  >>> A
  [-85, -55, -37]
  [-35,  97,  50]
  [ 79,  56,  49]

Travaux pratiques

  • Faire de même avec

$ A = \left( \begin{array}{ccccc} 8 & 5 & 1 & 7 & 9 \\ 9 & 2 & 1 & 7 & 1 \\ 6 & 1 & 4 & 4 & 8 \\ 7 & 5 & 4 & 9 & 6 \\ 9 & 6 & 7 & 9 & 6\end{array}\right)$

  • Faire une fonction qui retourne la décomposition LU d'une matrice de taille 4, si cette décomposition existe.

2.5  Utilisation des opérations élémentaires

Introduction

On sait qu'une matrice carrée A inversible est transformable en une matrice triangulaire supérieure à coefficients diagonaux non nuls par une succession d'opérations dites élémentaires sur les lignes.

Ces opérations peuvent être du type :

  • ajouter à la ligne i le produit de la ligne j par une coefficient $\beta$ : $L_i \leftarrow L_i + \beta L_j$,
  • multiplier la ligne $L_i$ par $\beta$ : $L_i \leftarrow \beta L_i$,
  • échanger deux lignes $L_i$ et $L_j$ : $L_i \leftrightarrow L_j$.

On trouve les mêmes opérations sur les colonnes.

Travaux pratiques

  1. Définir des fonctions addrow$(A,i,j,\beta)$, mulrow$(A,i,\beta)$ et swaprow$(A,i,j)$, qui réalisent chacune des opérations élémentaires ci-dessus sur la matrice $A$.
  2. Faire de même des fonctions addcol, mulcol, et swapcol, qui agissent sur les colonnes de $A$.

Les opérations élémentaires sous une forme matricielle

Chaque opération élémentaire sur les lignes de $A$ peut être interprétée comme le produit à gauche de $A$ par une matrice $Q$, résultat de cette opération sur la matrice identité.

Vérifions-le sur un exemple...

  >>> from sympy import *
  >>> def addrow(A,i,j,b):
  ...     A[i,:] += b*A[j,:]
  ... 
  >>> A = randMatrix(4,4)
  >>> A
  [68, 18, 58, 47]
  [84, 34, 60, 55]
  [35, 52, 10,  7]
  [34, 98, 10, 21]

Appliquons notre opération élémentaire sur une matrice identité, et multiplions cette dernière par $A$ (à gauche) :

  >>> Q = eye(4)
  >>> addrow(Q,1,2,3)
  >>> Q
  [1, 0, 0, 0]
  [0, 1, 3, 0]
  [0, 0, 1, 0]
  [0, 0, 0, 1]
  >>> Q*A
  [ 68,  18, 58, 47]
  [189, 190, 90, 76]
  [ 35,  52, 10,  7]
  [ 34,  98, 10, 21]

Vérifions que l'on trouve la même chose en appliquant directement l'opération élémentaire à $A$ :

  >>> addrow(A,1,2,3)
  >>> A
  [ 68,  18, 58, 47]
  [189, 190, 90, 76]
  [ 35,  52, 10,  7]
  [ 34,  98, 10, 21]

C'est le cas.

Décomposition LU et pivot de Gauss

Quand on fait un pivot de Gauss pour résoudre $AX=B$, on réalise des opérations élémentaires sur les lignes de $A$, pour faire en sorte que le système devienne triangulaire supérieur.

  • Le produit de ces opérations élémentaire nous donne notre matrice $L$.
  • $U$ est la partie triangulaire supérieure.

Faire une décomposition LU revient donc à faire un pivot de Gauss. Et si l'on suppose, dans une première approximation, que l'on ne rencontre pas de pivot nul, les seules opérations élémentaires que l'on exécutera seront de la forme $L_i \leftarrow L_i + \beta L_j$.

Les matrices $Q_k$ qui correspondent à ces opérations sont triangulaires inférieures à diagonale unité : elles ne diffèrent de l'unité que par l'élément $(i,j)$, avec $i>j$ (c.f. méthode du pivot).

Le produit $Q$ de ces matrices $Q_k$ reste triangulaire inférieur à diagonale unité, et on a $U=QA$ (le but de l'algorithme du pivot est d'obtenir une matrice triangulaire supérieure, que l'on note ici $U$).

Soit $L=Q^{-1}$, qui reste triangulaire inférieur à diagonale unité. Alors $A = Q^{-1} U = LU$ : on trouve notre décomposition.

Ainsi, pour obtenir la décomposition LU d'une matrice $A$, on peut lui appliquer un algorithme du pivot : le résultat (la matrice triangulaire) sera $U$, et le produit des matrices élémentaires que l'on a utilisé sera $L^{-1}$.

Reste à savoir comment obtenir $L$ à partir des opérations élémentaires. On a vu que $L^{-1} = Q = Q_1 \times ... \times Q_p$. Donc $L = I \times Q_p^{-1} \times ... \times Q_1^{-1}$, et on peut remarquer que si $Q_p$ correspond à l'opération élémentaire $L_i \leftarrow L_i + \beta L_j$, alors $Q_p^{-1}$ sera de la forme $C_j \leftarrow C_j - \beta C_i$.

Exemple

Réalisons la décomposition LU de la matrice

$A = \left( \begin{array}{ccc} 4 & 1 & 3 \\ 5 & 2 & 7 \\ 3 & 6 & 2 \end{array}\right)$

en utilisant l'algorithme du pivot de Gauss (la matrice a été choisie pour qu'il n'y ait pas de pivot nul, ce qui implique que la décomposition LU existe).

Pour ce faire, on utilisera les fonctions addrow et addcol précédemment définies, ainsi que la bibliothèque sympy.

On définit $A$, et on prend

  • pour valeur initiale de $U$ la matrice $A$,
  • pour valeur initiale de $L$ la matrice identité de taille 3...
  >>> A = Matrix([[4,1,3],[5,2,7],[3,6,2]])
>>> U = A
  >>> L = eye(3)
  >>> U
  [4, 1, 3]
  [5, 2, 7]
  [3, 6, 2]

On commence alors l'algorithme du pivot : on supprime le coefficient de la deuxième ligne, première colonne en utilisant pour pivot le coefficient en première ligne, première colonne ($L_2 \leftarrow L_2 - \frac{5}{4} L_1$) :

  >>> addrow(U,1,0,-Rational(5,4))
  >>> U
  [4,   1,    3]
  [0, 3/4, 13/4]
  [3,   6,    2]

L'opération associée est appliquée à $L$, mais sur les colonnes ($C_1 \leftarrow C_1 + \frac{5}{4} C_2$) :

  >>> addcol(L,0,1,Rational(5,4))
  >>> L
  [  1, 0, 0]
  [5/4, 1, 0]
  [  0, 0, 1]

On continue notre algorithme du pivot, pour supprimer le dernier coefficient de la première colonne, en utilisant pour pivot toujours le premier coefficient de cette colonne (soit $L_3 \leftarrow L_3 - \frac{3}{4} L_1$) :

  >>> addrow(U,2,0,-Rational(3,4))
  >>> U
  [4,    1,    3]
  [0,  3/4, 13/4]
  [0, 21/4, -1/4]

Et on reproduit l'opération associée aux colonnes de $L$ :

  >>> addcol(L,0,2,Rational(3,4))

On poursuit l'algorithme jusqu'à son terme...

  >>> addrow(U,2,1,-7)
  >>> U
  [4,   1,    3]
  [0, 3/4, 13/4]
  [0,   0,  -23]
  >>> addcol(L,1,2,7)
  >>> L
  [  1, 0, 0]
  [5/4, 1, 0]
  [3/4, 7, 1]

On a donc bien les formes espérées pour nos matrices $U$ et $L$. Reste à s'assurer que leur produit redonne bien $A$ :

  >>> L*U
  [4, 1, 3]
  [5, 2, 7]
  [3, 6, 2]

Travaux pratiques

  • Faire de même avec la matrice

$A = \left( \begin{array}{cccc} 5 & 7 & 2 & 2 \\ 4 & 3 & 1 & 6 \\ 2 & 8 & 3 & 4 \\ 4 & 9 & 1 & 5 \end{array}\right)$

  • Faire une fonction qui réalise la décomposition LU d'une matrice donnée en argument, par la méthode du pivot de Gauss.
On lèvera une exception PivotNul si un pivot nul est rencontré.

2.6  Cas général

Là où une matrice de permutation est nécessaire

On ne peut pas trouver de décomposition LU dans le cas où la matrice $A$, même inversible, n'a pas tout ses mineurs principaux non nuls.

On peut cependant montrer que, si $A$ est inversible, il existe une matrice de permutation $P$ telle que $PA$ admette une décomposition LU. Cette matrice de permutation correspond aux échanges de lignes qu'on est amené à faire dans la méthode du pivot, à chaque fois qu'un pivot est nul. On rappelle qu'on échange alors la ligne concernée avec une de celles situées en-dessous, de manière a récupérer un pivot non nul : c'est forcément possible, vu que $A$ est inversible.

L'ajout des matrices élémentaires correspondant aux $L_i \leftrightarrow L_j$ (pour $U$, et $C_j \leftrightarrow C_i$ pour $L$) aboutit à une matrice $L$ qui n'est plus triangulaire, que l'on notera $M$. Soit $P$ la matrice de permutation correspondant aux différents échanges de lignes.

On a donc $A = MU$, et $PA = PMU$. Mais alors PM redevient de la forme $L$ attendue ($PM$ contient deux fois $P$ : elles s'annulent)...

La matrice de permutation, toujours souhaitée

Il est préférable d'utiliser tout le temps le plus grand des pivots de la colonne considérée. En effet, vu que l'on divisera les coefficients de la ligne du pivot par ce dernier : plus ce pivot sera petit, plus son inverse sera grand, ce qui a pour conséquence...

  • la possibilité d'un dépassement de mémoire,
  • une petite erreur sur le pivot peut entraîner une grande erreur sur son inverse : on s'éloigne de la solution recherchée.

Pour minimiser ces effets, on intervertis donc à chaque étape la ligne courante avec la ligne de plus grand pivot : on revoit apparaître, comme ci-dessus, la matrice de permutation $P$.

Travaux pratiques

  1. Mettre en place cette permutation.
    • On renverra un message d'erreur si la matrice $A$ n'est pas inversible.
    • On construira la matrice de permutation $P$ petit à petit, à partir de la matrice identité (il suffit de lui appliquer tous les échanges de lignes rencontrés dans la méthode).
  2. Comparer l'ancienne décomposition LU (sans permutation) et la nouvelle sur la matrice

$\left( \begin{array}{ccc} 10^{-9} & 27 & 46 \\ 60 & 38 & 24 \\ 15 & 75 & 32 \end{array} \right)$

On effectuera les calculs avec numpy, pour éviter le calcul exact, le but étant de faire apparaître l'apparition d'erreurs catastrophiques, en calcul approché, quand on ne sélectionne pas systématiquement le meilleur pivot.

2.7  Algorithme de Crout

Présentation

L'algorithme de Crout :

  • applique une méthode de pivot partiel, comme dans la section précédente,
  • minimise la quantité de mémoire utilisée.

Il est donc le bienvenu dans une implémentation d'une décomposition LU d'une grande matrice.

Principes

Son principe est le suivant : comme

  • $L$ est triangulaire inférieure à diagonale unité, ses coefficients supérieurs étant nuls,
  • $U$ est triangulaire supérieure,

on se rend compte que $L$ et $U$ peuvent être stockés dans une seule matrice (les 0 de $L$ et $U$, comme les 1 sur la diagonale de $L$, peuvent ne pas être mémorisés).

De plus, il n'est pas nécessaire de former la matrice $P$ : il suffit de conserver la liste des numéros de lignes à permuter.

Les matrices utilisées dans la vraie vie ayant parfois des tailles gigantesques, ce gain en mémoire (une matrice au lieu de trois) peut être indispensable, et peut aussi entraîner un gain de temps.

Travaux pratiques

  1. Améliorer votre algorithme pour économiser de l'espace mémoire.
  2. Comparer la rapidité des différentes méthodes.

3.  La décomposition QR

3.1  Préliminaire : procédé d'orthonormalisation

Le But

Le procédé d'orthonormalisation consiste à prendre une base $(b_1, ..., b_n)$ de l'espace vectoriel considéré, et à la transformer en une base orthonormée $(e_1, ..., e_n)$ :

  • les vecteurs sont deux à deux orthogonaux,
  • ils sont de norme égale à 1.

Ce procédé se généralise immédiatement à toute famille libre $(b_1, ..., b_k)$ qui peut ainsi être transformée en une famille (libre) orthonormée $(e_1, ..., e_k)$ engendrant le même espace.

En fait, on peut montrer que $(e_1, ..., e_k)$ est l'unique famille orthonormée engendrant le même espace que $(b_1, ..., b_k)$, telle que $<e_i,b_i>>0$ (produit scalaire strictement positif, i.e. les vecteurs "pointent dans la même direction").

Finalement, la méthode de Gram-Schmidt consiste à former une famille orthonormale en "redressant", puis en normant, progressivement chacun des vecteurs fournis, pour les rendre orthogonaux aux vecteurs déjà formés.

L'Algorithme

Le principe de cette orthonormalisation se comprend aisément :

  1. On divise $b_1$ par sa norme, ce qui donne $e_1$ :

$e_1 = \frac{b_1}{||b_1||}$

Ainsi, notre premier vecteur de la nouvelle base sera bien normé.
  1. On annule la composante de $b_2$ suivant $e_1$, pour que le résultat soit orthogonal à $e_1$, puis on divise le résultat par sa norme :

$e_2 = \frac{b_2 - <b_2,e_1> e_1}{||b_2-<b_2,e_1> e_1||}$

En effet, $b_2 - <b_2,e_1> e_1$ est orthogonal à $e_1$ (prendre leur produit scalaire pour s'en convaincre) : $(e_1, e_2)$ forme bien une famille orthonormale.
  1. On annule la composante de $b_3$ suivant $e_1$ et $e_2$, pour que le vecteur résultant soit bien orthogonal à ces deux vecteurs, et on divise le résultat par sa propre norme :

$e_3 = \frac{b_3 - <b_3,e_1> e_1 - <b_3,e_2> e_2}{||b_3 - <b_3,e_1> e_1 - <b_3,e_2> e_2 e_1||}$

  1. On continue jusqu'à épuiser tous les vecteurs $b_i$...

$e_i = \frac{b_i - \sum_{j=1}^{i-1} <b_i,e_j> e_j}{||b_i - \sum_{j=1}^{i-1} <b_i,e_j> e_j||}$

Gram-Schmidt et sympy

Montrons comment réaliser le procédé d'orthonormalisation de Gram-Schmidt avec python-sympy.

Préliminaires

Commençons par définir trois vecteurs à orthonormaliser :

  >>> from sympy import *
  >>> B1 = Matrix((2,1,2)).reshape(3,1)
  >>> B2 = Matrix((12,-6,0)).reshape(3,1)
  >>> B3 = Matrix((-3,-3,18)).reshape(3,1)

Avant toutes choses, vérifions que l'on a bien une base. Nous allons former une matrice avec ces trois vecteurs, puis nous assurer que le déterminant de cette matrice est bien non nul :

  >>> M=B1.row_join(B2).row_join(B3)
  >>> M
  [2, 12, -3]
  [1, -6, -3]
  [2,  0, 18]
  >>> M.det()
  -540

Ces vecteurs forment donc bien une base !

Fonction prédéfinie

Pour orthonormaliser façon Gram-Schmidt en sympy :

  >>> GramSchmidt([B1,B2,B3], orthog = True)
  [ [2/3]
  [1/3]
  [2/3], [ 2/3]
  [-2/3]
  [-1/3], [-1/3]
  [-2/3]
  [ 2/3] ]

Si on a juste besoin d'orthogonaliser (pas normer), ne pas mettre le paramètre orthog :

  >>> GramSchmidt([B1,B2,B3])
  [ [2]
  [1]
  [2], [ 8]
  [-8]
  [-4], [ -5]
  [-10]
  [ 10] ]

Le produit scalaire entre vecteurs s'obtient avec la méthode dot, et la norme avec norm. On peut donc vérifier que l'on a bien une base orthogonale non normée :

  >>> (u,v,w) = GramSchmidt([B1,B2,B3])
  >>> u.dot(v) 
  0
  >>> u.dot(w)
  0
  >>> v.dot(w) 
  0
  >>> u.norm()
  3

A la main

On a tout ce qu'il faut pour orthonormaliser à la main. Commençons par calculer le premier vecteur $e_1$ de la nouvelle base (c'est $b_1$ divisé par sa norme) :

  >>> e1 = B1/B1.norm()
  >>> e1
  [2/3]
  [1/3]
  [2/3]
  >>> e1.norm()
  1

Continuons le procédé, en suivant l'algorithme ci-dessus. On vérifiera, à chaque étape, la norme, et le produit scalaire avec les autres vecteurs :

  >>> e2 = B2-B2.dot(e1)*e1
  >>> e2 = e2/e2.norm()
  >>> e2
  [ 2/3]
  [-2/3]
  [-1/3]
  >>> e2.norm()
  1
  >>> e1.dot(e2)
  0
  >>> e3 = B3-B3.dot(e1)*e1-B3.dot(e2)*e2
  >>> e3 = e3/e3.norm()
  >>> e3
  [-1/3]
  [-2/3]
  [ 2/3]
  >>> e3.norm()
  1
  >>> e3.dot(e1)
  0
  >>> e3.dot(e2) 
  0

On a donc réussi, en appliquant le procédé d'orthonormalisation de Gram-Schmidt, à transformer la base

$\left(\begin{array}{c} 2 \\ 1 \\ 2 \end{array}\right), \left(\begin{array}{c} 12 \\ -6 \\ 0 \end{array}\right), \left(\begin{array}{c} -3 \\ -3 \\ 18 \end{array}\right)$

en la base orthonormée :

$\left(\begin{array}{c} \frac{2}{3} \\ \frac{1}{3} \\ \frac{2}{3} \end{array}\right), \left(\begin{array}{c} \frac{2}{3} \\ -\frac{2}{3} \\ -\frac{1}{3} \end{array}\right), \left(\begin{array}{c} -\frac{1}{3} \\ -\frac{2}{3} \\ \frac{2}{3} \end{array}\right)$

Travaux pratiques

  1. Faire une fonction qui reçoit une liste de vecteurs en entrée, et qui les orthonormalise.
  2. Faire en sorte que votre fonction puisse normaliser, ou non.

3.2  Présentation de la décomposition QR

La décomposition

Soit $A$ une matrice carrée, à coefficients réels, inversible. On appelle décomposition QR de $A$ toute égalité $A = QR$

où :

  • $Q$ est une matrice orthogonale : $Q^{-1} = Q^T$,
  • $R$ est une matrice triangulaire supérieure à coefficients diagonaux strictement positifs.

Cette décomposition existe (pour $A$ inversible) et est unique.

Avec sympy

Définissons une matrice $M$ qui possède bien une décomposition QR, c'est-à-dire qui est bien inversible (déterminant non nul) :

  >>> M
  [2, 12, -3]
  [1, -6, -3]
  [2,  0, 18]
  >>> M.det()
  -540

Pour obtenir la décomposition QR de $M$, on utilise la méthode QRdecomposition :

  >>> (Q,R) = M.QRdecomposition()
  >>> Q
  [2/3,  2/3, -1/3]
  [1/3, -2/3, -2/3]
  [2/3, -1/3,  2/3]
  >>> R
  [3,  6,  9]
  [0, 12, -6]
  [0,  0, 15]

$R$ est bien une matrice triangulaire supérieure à coefficients diagonaux strictement positifs. Vérifions que $Q$ est bien orthogonale (c'est-à-dire que $Q^{-1}$ est égale à la transposée de $Q$) :

  >>> Q.inv()
  [ 2/3,  1/3,  2/3]  
  [ 2/3, -2/3, -1/3]
  [-1/3, -2/3,  2/3]
  >>> Q.T
  [ 2/3,  1/3,  2/3]
  [ 2/3, -2/3, -1/3]
  [-1/3, -2/3,  2/3]

C'est le cas.

Dernière chose : $Q \times R$ est bien égale à $M$ :

  >>> Q*R
  [2, 12, -3]
  [1, -6, -3]
  [2,  0, 18]
  >>> Q*R == M
  True

Lien avec Gram-Schmidt

L'écriture QR d'une matrice inversible $A$ est intimement liée à l'orthonormalisation de la famille des vecteurs colonnes de $A$ par Gram-Schmidt.

Soient $(a_1, ..., a_n)$ les vecteurs colonnes de $A$ ; ils forment une base, vu que $A$ est inversible. $A$ peut être interprétée comme la matrice de passage de la base canonique $(e_1, ..., e_n)$ à la base $(a_1, ..., a_n)$.

Notons $(q_1, ..., q_n)$ le résultat (la base orthonormée) du procédé de Gram-Schmidt appliqué à $(a_1, ..., a_n)$.

Soient :

  • $Q$ la matrice de passage de la base canonique à la base $(q_1, ..., q_n)$.
$Q$ est une matrice orthogonale, car c'est une matrice de passage entre deux bases orthonormées.
  • $R$ la matrice de passage de la base $(q_1, ..., q_n)$ à la base $(a_1, ..., a_n)$.
    • La matrice $R$ est triangulaire supérieure, vu que chaque $a_i$ s'exprime uniquement sur les vecteurs $q_1, ..., q_i$ (c'est une conséquence du procédé d'orthonormalisation).
    • Pour tout $i \leqslant j$, $R_{i,j}$ représente la composante de $a_j$ sur $q_i$, donc est égal à $<q_i,a_j>$, puisque la base $(q_1, ..., q_n)$ est orthonormée.
    • En particulier, $R_jj = <q_j,a_j> >0$, vu le procédé d'orthonormalisation.

L'égalité $A=Q \times R$ est une conséquence immédiate du fait que, ainsi construite, $A, Q$ et $R$ sont des matrices de passage entre bases.

Pour conclure : la décomposition QR est simplement une traduction du procédé d'orthonormalisation de Gram-Schmidt appliqué aux vecteurs de $A$.

Premier algorithme de décomposition

Les remarques ci-dessus nous permettent de construire un premier algorithme de décomposition QR :

  1. La matrice $Q$ a pour jième vecteur colonne le vecteur $q_j = \frac{a_j - \sum_{i=1}^{j-1} <a_j,q_i> q_i}{||a_j - \sum_{i=1}^{j-1} <a_j,q_i> q_i||}$
  2. Pour tout $i \leqslant j$, $R_{i,j} = <q_i,a_j>$.

Mise en oeuvre du premier algorithme avec sympy

On va montrer comment obtenir la décomposition QR à la main avec sympy, en utilisant l'algorithme basé sur Gram-Schmidt.

Considérons la matrice

$\left( \begin{array}{ccc} 2 & 12 & -3 \\ 1 & -6 & -3 \\ 2 & 0 & 18 \end{array} \right)$

et définissons-là sous python :

  >>> from sympy import *
  >>> A = Matrix((2,12,-3,1,-6,-3,2,0,18)).reshape(3,3)

Cette matrice possède une décomposition QR, vu que son déterminant est non nul :

  >>> A.det()
  -540

Occupons-nous déjà de $Q$ : sa première colonne est la première colonne de $A$, normée :

  >>> Q = zeros(3)
  >>> Q[:,0] = A[:,0]/A[:,0].norm()

Ses colonnes suivantes s'obtiennent comme dans Gram-Schmidt...

  >>> Q[:,1] = A[:,1]-A[:,1].dot(Q[:,0])*Q[:,0]
  >>> Q[:,1] = Q[:,1]/Q[:,1].norm()
  >>> Q
  [2/3,  2/3, 0]
  [1/3, -2/3, 0]
  [2/3, -1/3, 0]
  >>> Q[:,2] = A[:,2]-A[:,2].dot(Q[:,0])*Q[:,0]-A[:,2].dot(Q[:,1])*Q[:,1]
  >>> Q[:,2] = Q[:,2]/Q[:,2].norm()
  >>> Q
  [2/3,  2/3, -1/3]
  [1/3, -2/3, -2/3]
  [2/3, -1/3,  2/3]

On vérifie que $Q$ est bien orthogonale ($Q^{-1} = Q^T$) :

  >>> Q.inv()==Q.T
  True

Il reste à définir $R$, grâce au fait que $R_{i,j} = <Q_i,A_j>$ :

  >>> R = zeros(3)
  >>> for i in range(3):
  ...     for j in range(i,3):
  ...         R[i,j] = Q[:,i].dot(A[:,j])
  ... 
  >>> R
  [3,  6,  9]
  [0, 12, -6]
  [0,  0, 15]

$R$ est bien triangulaire supérieure à diagonale strictement positive. Il reste à vérifier que $Q \times R$ redonne bien $A$ :

  >>> Q*R==A
  True

Améliorations possibles

Une fois que l'on a obtenu $Q$, il n'est pas nécessaire de faire les calculs ci-dessus pour $R$. En effet, $A = QR$, donc $R = Q^{-1} A$. Mais $Q$ est orthogonale, donc $Q^{-1} = Q^T$. D'où, pour obtenir $R$, on fait $Q^T \times A$.

D'autres améliorations sont possibles pour obtenir $Q$ avec moins de calculs.

Travaux pratiques

Faire une fonction qui retourne la décomposition QR d'une matrice $A$.

3.3  Méthode de Givens

Présentation

Le calcul formel a ses limites, et l'on est souvent obligé de revenir au calcul numérique pour réaliser notre décomposition QR. Aussi, le comportement des méthodes envisagées face aux erreurs d'arrondis prend de l'importance.

La méthode de Gram-Schmidt n'est numériquement pas très stable. Aussi, on cherche d'autres méthodes, minimisant les erreurs de calculs.

La méthode de Givens, plus stable, consiste à transformer la matrice $A$ en une matrice triangulaire supérieure $R$, en lui appliquant (par produit à gauche) des matrices de rotations, qui sont orthogonales, de manière à annuler petit à petit tous les coefficients sous-diagonaux.

L'idée est la suivante :

  • on peut annuler la dernière composante du premier vecteur de la matrice $A$, en le faisant tourner autour d'un axe orthogonal au dernier vecteur de base.
  • On continue ainsi avec les autres composantes.

Les matrices de rotation

Les matrices de rotation, autour des axes canoniques de $\mathbb{R}^6$, ont la forme suivante :

$\left( \begin{array}{cccccc} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & \cos(\theta) & 0 & 0 & sin(\theta) & 0 \\ 0 & 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & -sin(\theta) & 0 & 0 & cos(\theta) & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{array}\right)$

et on peut généraliser cette forme à $\mathbb{R}^n$.

Définissons une fonction python qui retourne de telles matrices :

  >>> def G(n,i,j,c,s):
  ...     M = eye(n)
  ...     M[i,i] = M[j,j] = c
  ...     M[i,j],M[j,i] = -s, s
  ...     return M
  ... 
  >>> x = Symbol('x')
  >>> G(6,2,5,cos(x),sin(x))
  [1, 0,       0, 0, 0,      0]
  [0, 1,       0, 0, 0,      0]
  [0, 0,  cos(x), 0, 0, sin(x)]
  [0, 0,       0, 1, 0,      0]
  [0, 0,       0, 0, 1,      0]
  [0, 0, -sin(x), 0, 0, cos(x)]

Action d'une matrice de rotation

On constate que les matrices de rotation G(n,i,j,c,s) n'agissent que sur les lignes $i$ et $j$ de $A$, quand on fait le produit à gauche $G(n,i,j,c,s) \times A$ :

  >>> s,c = Symbol('s'),Symbol('c')
  >>> G(6,2,5,c,s)*randMatrix(6,6)
  [          77,           59,           57,           69,          71,            8]
  [          99,           26,           65,           10,          71,           91]
  [52*s + 50*c, 78*s + 55*c, 88*s + 59*c, 17*s + 30*c, 85*s + 6*c, 42*s + 34*c]
  [          68,           10,           26,           53,          86,           74]
  [          69,           10,           31,            6,          26,            3]
  [-50*s + 52*c, -55*s + 78*c, -59*s + 88*c, -17*c + 30*s, -6*s + 85*c,  -34*s + 42*c]

Notre but est, on le rappelle, d'obtenir une matrice triangulaire, en multipliant différentes matrices $G(n, i, j, cos(\theta), sin(\theta)$ à gauche de $A$.

On constate, sur l'exemple ci-dessus, que cela est réalisable. On peut, en effet, annuler le dernier coefficient de la première colonne, en trouvant $c$ et $s$ tels que :

  • $50*s+52*c = 0$, pour annuler ledit coefficient,
  • $c^2+s^2 = 1$, pour qu'un $\theta$ tel que $c = cos(\theta), s = sin(\theta)$ existe bien.

En posant $r = \sqrt{50^2+52^2}$, et $c = \frac{50}{r}, s = \frac{52}{r}$, on trouve la solution au système, et on annule le coefficient en bas de la première colonne.

Vers un algorithme

D'une manière générale, si on veut annuler l'élément en position $(i,j)$ dans $A$ à l'aide d'une matrice de rotation, on peut procéder ainsi :

  • Si $r = \sqrt{A_{j,j}^2 + A_{i,j}^2}$ est nul, alors $G(n,i,j,1,0)$ convient, puisque le coefficient est déjà nul.
  • Sinon, on pose $c = \frac{A_{j,j}}{r}$ et $s = \frac{A_{i,j}}{r}$. Alors $G(n,i,j,c,s)$ convient.

On conçoit donc l'idée d'un algorithme pour obtenir la décomposition QR : on annule les coefficients de la première colonne, du dernier au second, en multipliant $A$ à gauche par les bonnes matrices de rotation. Puis on annule les coefficients de la deuxième colonne, du dernier au troisième, en multipliant, dans le cas de la colonne, etc.

Le résultat de ces opérations, noté $R = G_m \times ... \times G_1 \times A$, sera par construction triangulaire supérieur.

Reste alors à retrouver $Q$. Comme $R = G_m \times ... \times G_1 \times A$, on a $A = G_m^{-1} \times G_{m-1}^{-1} ... G_1^{-1} R = Q \times R$, avec $Q = G_m^{-1} \times G_{m-1}^{-1} ... G_1^{-1}$.

Enfin, les matrices $G_1, ..., G_m$ sont des matrices de rotation, donc sont orthogonales, et leurs inverses sont leurs transposées : $Q = G_m^{T} \times G_{m-1}^{T} ... G_1^{T}$

L'algorithme

L'algorithme de Givens peut maintenant être écrit :

  $Q$ est l'identité de taille $n$
  pour j de 1 à n-1 :
      pour i de j+1 à n :
          calculer $r = \sqrt{A_{j,j}^2 + A_{i,j}^2}$.
          si $r$ n'est pas nul, alors
              calculer $G = G\left(n,i,j, \frac{A_{j,j}}{r}, \frac{A_{i,j}}{r} \right)$
              remplacer $A$ par $G\times A$,
              remplacer $Q$ par $Q \times G^T$.

L'algorithme se termine quand tous les coefficients sous-diagonaux ont tous été annulés.

Il reste cependant un dernier point à s'assurer : dans la décomposition QR, $R$ est sensée être triangulaire supérieure à coefficients diagonaux strictement positifs.

On peut prouver que, après produit à gauche par une matrice de rotation $G(n,i,j,c,s)$, avec les coefficients $c,s$ bien choisis :

  • $A_{i,j}$ s'annule (on le sait déjà, c'est le but de la manipulation),
  • $A_{i,i}$ devient strictement positif.

Ainsi, tous les coefficients diagonaux de la matrice construite deviennent strictement positifs...sauf, éventuellement le dernier, puisqu'on ne peut faire tourner notre algorithme que jusqu'à l'avant-dernière colonne (il n'y a pas de coefficient en-dessous du terme $(n,n)$).

Il convient donc, pour finir notre algorithme, de s'assurer que l'élément en position $(n,n)$ de la matrice $R$ obtenue ci-dessus est bien positif. Si tel n'était pas le cas, il faut :

  • Changer le signe de $R_{n,n}$,
  • Changer le signe de la dernière colonne de $Q$.

Réalisation avec sympy

On va réaliser une décomposition QR par méthode de Givens, avec la matrice suivante...

  >>> A = Matrix((2,12,-3,1,-6,-3,2,0,18)).reshape(3,3)
  >>> A
  [2, 12, -3]
  [1, -6, -3]
  [2,  0, 18]

On initialise notre algorithme : $R$ est égal à $A$, et $Q$ à la matrice identité :

  >>> R = Matrix((2,12,-3,1,-6,-3,2,0,18)).reshape(3,3)
  >>> Q = eye(3)

Premier passage dans la boucle. On calcule $r$, qui ici est non nul. On multiplie alors $R$ à gauche par la matrice de rotation adéquate :

  >>> r = sqrt(R[0,0]**2+R[2,0]**2)
  >>> r
  2*2**(1/2)
  >>> g = G(3,2,0,R[0,0]/r,R[2,0]/r)
  >>> R=g*R
  >>> R
  [2*2**(1/2),  6*2**(1/2), 15*2**(1/2)/2]
  [         1,          -6,            -3]
  [         0, -6*2**(1/2), 21*2**(1/2)/2]

On a bien fait apparaître un 0 là où il fallait. Reportons la matrice de rotation (ou plutôt sa transposée) dans $Q$ :

  >>> Q = Q*g.T

On recommence, pour annuler le coefficient en première colonne, deuxième ligne, par la rotation qu'il faut (qui dépend de la nullité de $r$), et on n'oublie pas de reporter ce qu'il faut dans $Q$ :

  >>> r = sqrt(R[0,0]**2+R[1,0]**2)
  >>> r
  3
  >>> g = G(3,1,0,R[0,0]/r,R[1,0]/r)  
  >>> R=g*R
  >>> R
  [3,           6,             9]
  [0, -6*2**(1/2), -9*2**(1/2)/2]
  [0, -6*2**(1/2), 21*2**(1/2)/2]
  >>> Q = Q*g.T

Le 0 est bien apparu là où il fallait, et celui qui avait précédemment été introduit n'a pas sauté ! Continuons...

  >>> r = sqrt(R[1,1]**2+R[2,1]**2)
  >>> r
  12
  >>> g = G(3,2,1,R[1,1]/r,R[2,1]/r)
  >>> R=g*R
  >>> R
  [3,  6,   9]
  [0, 12,  -6]
  [0,  0, -15]
  >>> Q = Q*g.T

On a fini notre boucle. Comme le coefficient de la dernière diagonale est négatif :

  • on change son signe,
  • on change le signe de la dernière colonne de $Q$.

Ce qui donne...

  >>> if R[2,2]<0:
  ...     R[2,2] = -R[2,2]
  ...     Q[:,2] = -Q[:,2]
  ... 

$R$ est de la forme attendue. Vérifions que $Q$ est bien orthogonale, et que le produit $Q \times R$ redonne bien $A$...

  >>> Q.inv() == Q.T
  True
  >>> Q*R == A
  True

On a réussi notre décomposition par la méthode de Givens.

Travaux pratiques

Faire une fonction qui réalise, si cela est possible, la décomposition QR d'une matrice $A$, par la méthode de Givens.

3.4  Méthode de Householder

Présentation

La méthode de Givens permet de former la décomposition QR d'une matrice inversible $A$, en lui appliquant (par produit à gauche) des matrices de rotations.

La méthode de Householder procède de la même façon, mais en utilisant des matrices de symétries vectorielles par rapport à des hyperplans.

Là où la méthode de Givens consiste à annuler un par un chacun des coefficients sous-diagonaux de la matrice $A$, la méthode de Householder procède à l'annulation simultanée de tous les coefficients sous-diagonaux d'une même colonne, à chaque étape.

Les matrices de Householder

Définition : Soit $V$ un vecteur non nul de $\mathbb{R}^n$. La matrice $H_V = I_n - \frac{2}{||V||^2} V V^T$ est appelée matrice de Householder du vecteur $V$.

$H_V$ représente, d'un point de vue géométrique, une symétrie orthogonale par rapport à l'hyperplan orthogonal à $V$.

Travaux pratiques

  1. Faire une fonction qui, à un vecteur $V$, associe sa matrice de Householder.
  2. Vérifier que la matrice de Householder d'un vecteur $V$ de votre choix est : symétrique, orthogonale, et qu'elle transforme $V$ en son opposé.

Trouver le bon hyperplan

Le but est similaire à la méthode de Givens : trouver le bon paramètre, pour annuler des coefficients sous-diagonaux.

Ici, on souhaite trouver, pour un vecteur $U$ donné, une matrice $H_V$ (c'est-à-dire un vecteur $V$) tel que $H_V(U)$ conserve les premières composantes de $U$, et annule les dernières, en trouvant la symétrie qu'il faut.

Ce vecteur $V$ existe toujours, voici comment le construire :

  • Notons $(u_1, ..., u_n)$ les composantes de $U$, et supossons que l'on souhaite annuler ses composantes à partir de la position $j$.
  • Notons $P_j = (0, ..., 0, u_j, ..., u_n)$ la projection de $U$ sur les vecteurs $(e_j, ..., e_n)$ de la base canonique, et $||P_j||$ sa norme euclidienne.
  • Alors $V = (0, ..., 0, u_j + \varepsilon ||P_j||, u_{j+1}, ..., u_n)$ convient, où $\varepsilon$ est le signe de $u_j$ (-1, 0, ou 1).

Travaux pratiques

  1. Faire une fonction qui reçoit un vecteur $U$ et un indice $i$ en entrée, et qui renvoie le vecteur $V$ définit comme ci-dessus.
  2. Définir une fonction qui reçoit une matrice $A$, et un indice $i$, et qui renvoie le produit de la matrice de Householder formée à partir de la jième colonne de $A$, par $A$.
  3. Faire la décomposition QR par la méthode de Householder :
    • On annulera d'un coup les composants de 2 à $n$ de la première colonne de $A$, en multipliant par une matrice de Householder bien choisie.
    • Puis on annulera les composants de 3 à $n$ de la deuxième colonne de la matrice résultante, toujours avec une bonne matrice de Householder.
    • etc.
    • On fera attention à bien renvoyer Q et R :
      • $R = H_{n-1} \times ... \times H_1 \times A$ est obtenue par produit sur $A$, à gauche, par les matrices de Householder.
      • $Q = I_n \times H_1 \times ... \times H_{n-1}$, vu que les $H_i$ sont orthogonales (similaire à la méthode de Givens).

Remarque finale

L'algorithme de Householder, programmé par la méthode précédente, fournit une matrice $R$ triangulaire supérieure, mais qui n'est pas forcément à diagonale strictement positive.

Pour faire en sorte qu'il en soit ainsi, on peut, en cours de méthode, lorsque l'on rencontre un $R_{k,k}<0$, procéder à des changements de signe simultanés :

  • de la ligne $k$ de $R$,
  • de la colonne $k$ de $Q$.

Travaux pratiques

  1. Adaptez votre algorithme pour qu'il prenne en considération cette ultime remarque.
  2. Retrouver les résultats de la méthode QRdecomposition de sympy, et de votre fonction Givens, sur les exemples de ce document, et sur

$\left( \begin{array}{cccc} 23 & -15 & 30 & 15 \\ -4 & -30 & -15 & 55 \\ -4 & 20 & 35 & 5 \\ -8 & 40 & -80 & -65 \end{array}\right)$

4.  Décomposition de Choleski

Terminons ce tour d'horizon des méthodes classiques de décomposition matricielle par la présentation de la décomposition de Choleski...

Soit $S$ une matrice carrée d'ordre $n$, symétrique ($S^T = S$) et définie positive (i.e. ses valeurs propres sont strictement positives).

Alors $S$ possède une unique décomposition de la forme $S = L L^T$, où $L$ est triangulaire inférieure à diagonale strictement positive : c'est la décomposition de Choleski.

Actions

Recent Changes