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.
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.
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.
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
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.
On peut installer sympy ainsi sous une debian GNU/Linux, ou sous une Ubuntu :
$ aptitude install python-sympy
dans un shell, en étant administrateur.
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.
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
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.
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
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
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}
Soit $A$ une matrice carrée d'ordre $n$ et inversible. On appelle décomposition LU de $A$ toute égalité $A=LU$, où
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
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).
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.
On peut prouver que :
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)$
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.
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 ... >>>
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).
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]
$ 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)$
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 :
On trouve les mêmes opérations sur les colonnes.
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.
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.
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$.
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
>>> A = Matrix([[4,1,3],[5,2,7],[3,6,2]])
>>> 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]
$A = \left( \begin{array}{cccc} 5 & 7 & 2 & 2 \\ 4 & 3 & 1 & 6 \\ 2 & 8 & 3 & 4 \\ 4 & 9 & 1 & 5 \end{array}\right)$
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)...
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...
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$.
$\left( \begin{array}{ccc} 10^{-9} & 27 & 46 \\ 60 & 38 & 24 \\ 15 & 75 & 32 \end{array} \right)$
L'algorithme de Crout :
Il est donc le bienvenu dans une implémentation d'une décomposition LU d'une grande matrice.
Son principe est le suivant : comme
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.
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)$ :
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.
Le principe de cette orthonormalisation se comprend aisément :
$e_1 = \frac{b_1}{||b_1||}$
$e_2 = \frac{b_2 - <b_2,e_1> e_1}{||b_2-<b_2,e_1> e_1||}$
$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||}$
$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||}$
Montrons comment réaliser le procédé d'orthonormalisation de Gram-Schmidt avec python-sympy.
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 !
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
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)$
Soit $A$ une matrice carrée, à coefficients réels, inversible. On appelle décomposition QR de $A$ toute égalité $A = QR$
où :
Cette décomposition existe (pour $A$ inversible) et est unique.
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
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 :
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$.
Les remarques ci-dessus nous permettent de construire un premier algorithme de décomposition QR :
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
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.
Faire une fonction qui retourne la décomposition QR d'une matrice $A$.
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 :
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)]
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 :
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.
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 :
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 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 :
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 :
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 :
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.
Faire une fonction qui réalise, si cela est possible, la décomposition QR d'une matrice $A$, par la méthode de Givens.
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.
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$.
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 :
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 :
$\left( \begin{array}{cccc} 23 & -15 & 30 & 15 \\ -4 & -30 & -15 & 55 \\ -4 & 20 & 35 & 5 \\ -8 & 40 & -80 & -65 \end{array}\right)$
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.