Factorisation de Cholesky

La factorisation de Cholesky, nommée d'après André-Louis Cholesky, consiste, pour une matrice symétrique définie positive A, à déterminer une matrice triangulaire inférieure L telle que : A = LL^T.

La matrice L est en quelque sorte une « racine carrée » de A. Cette décomposition permet notamment de calculer la matrice inverse A⁻¹, de calculer le déterminant de A (égal au carré du produit des éléments diagonaux de L) ou encore de simuler une loi multinormale. Elle est aussi utilisée en chimie quantique pour accélérer les calculs (voir Décomposition de Cholesky (chimie quantique)).

Exemple

[modifier | modifier le code ]

La matrice symétrique A :

${\begin{pmatrix}1&1&1&1\1円&5&5&5\1円&5&14&14\1円&5&14&15\\\end{pmatrix}}$ {\displaystyle {\begin{pmatrix}1&1&1&1\1円&5&5&5\1円&5&14&14\1円&5&14&15\\\end{pmatrix}}}

est égale au produit de la matrice triangulaire L :

${\begin{pmatrix}1&0&0&0\1円&2&0&0\1円&2&3&0\1円&2&3&1\\\end{pmatrix}}$ {\displaystyle {\begin{pmatrix}1&0&0&0\1円&2&0&0\1円&2&3&0\1円&2&3&1\\\end{pmatrix}}}

avec à sa droite sa transposée L^T :

${\begin{pmatrix}1&1&1&1\0円&2&2&2\0円&0&3&3\0円&0&0&1\\\end{pmatrix}}$ {\displaystyle {\begin{pmatrix}1&1&1&1\0円&2&2&2\0円&0&3&3\0円&0&0&1\\\end{pmatrix}}}

Théorème

[modifier | modifier le code ]

Théorème — Si A est une matrice symétrique définie positive, il existe une matrice réelle triangulaire inférieure L telle que^[1]:

A = LL^T.

De plus cette décomposition est unique si l'on impose à L d'avoir des coefficients diagonaux strictement positifs.

Démonstration

Notons $B=(e_{1},...e_{n})$ {\displaystyle B=(e_{1},...e_{n})} la base canonique de $\mathbb {R} ^{n}$ {\displaystyle \mathbb {R} ^{n}}.

A est symétrique définie positive, donc elle représente un certain produit scalaire dans B.

Soit $B'=(e_{1}',...e_{n}')$ {\displaystyle B'=(e_{1}',...e_{n}')} une base orthonormée pour ce produit scalaire. On note P la matrice de passage de B' vers B.

La formule de changement de base pour une application bilinéaire donne $A=P^{T}\mathrm {I} _{n}P=P^{T}P$ {\displaystyle A=P^{T}\mathrm {I} _{n}P=P^{T}P}

On va chercher à bien choisir B' de sorte que P soit triangulaire supérieure à coefficients diagonaux positifs.

Pour cela, on cherche B' orthonormée telle que :

\forall k\in \{1,...,n\},{\begin{cases}{\rm {Vect}}(e_{1},...,e_{k})={\rm {Vect}}(e_{1}',...,e_{k}')\\e_{k}.e_{k}'>0\end{cases}}

{\displaystyle \forall k\in \{1,...,n\},{\begin{cases}{\rm {Vect}}(e_{1},...,e_{k})={\rm {Vect}}(e_{1}',...,e_{k}')\\e_{k}.e_{k}'>0\end{cases}}}

L'algorithme de Gram-Schmidt garantit l'existence et l'unicité d'une telle base, d'où l'existence et l'unicité de P. Enfin on pose $L=P^{T}$ {\displaystyle L=P^{T}} dont on déduit la forme voulue.

Algorithme

[modifier | modifier le code ]

On cherche la matrice :

L={\begin{bmatrix}l_{11}\\l_{21}&l_{22}\\\vdots &\vdots &\ddots \\l_{n1}&l_{n2}&\cdots &l_{nn}\end{bmatrix}}

{\displaystyle L={\begin{bmatrix}l_{11}\\l_{21}&l_{22}\\\vdots &\vdots &\ddots \\l_{n1}&l_{n2}&\cdots &l_{nn}\end{bmatrix}}}

De l'égalité A = LL^T on déduit :

a_{ij}=\left(LL^{T}\right)_{ij}={\sum _{k=1}^{n}l_{ik}l_{jk}}={\sum _{k=1}^{\min \left\{i,j\right\}}l_{ik}l_{jk}},\;1\leq i,j\leq n

{\displaystyle a_{ij}=\left(LL^{T}\right)_{ij}={\sum _{k=1}^{n}l_{ik}l_{jk}}={\sum _{k=1}^{\min \left\{i,j\right\}}l_{ik}l_{jk}},\;1\leq i,j\leq n}

puisque l_pq=0 si 1 ≤ p < q ≤ n.

La matrice A étant symétrique, il suffit que les relations ci-dessus soient vérifiées pour i ≤ j, c'est-à-dire que les éléments l_ij de la matrice L doivent satisfaire :

a_{ij}={\sum _{k=1}^{i}l_{ik}l_{jk}},\;1\leq i\leq j\leq n

{\displaystyle a_{ij}={\sum _{k=1}^{i}l_{ik}l_{jk}},\;1\leq i\leq j\leq n}

Pour i = 1, on détermine la première colonne de L :

a_{11}=l_{11}l_{11}

{\displaystyle a_{11}=l_{11}l_{11}} d'où

l_{11}={\sqrt {a_{11}}}

{\displaystyle l_{11}={\sqrt {a_{11}}}}

a_{j1}=l_{11}l_{j1}

{\displaystyle a_{j1}=l_{11}l_{j1}} d'où

l_{j1}={\frac {a_{j1}}{l_{11}}},\;\;\;2\leq j\leq n

{\displaystyle l_{j1}={\frac {a_{j1}}{l_{11}}},\;\;\;2\leq j\leq n}

On détermine la i-ème colonne de L 2 ≤ i ≤ n, après avoir calculé les (i – 1) premières colonnes :

a_{ii}=l_{i1}l_{i1}+\ldots +l_{ii}l_{ii}

{\displaystyle a_{ii}=l_{i1}l_{i1}+\ldots +l_{ii}l_{ii}} d'où

l_{ii}={\sqrt {a_{ii}-{\sum _{k=1}^{i-1}l_{ik}^{2}}}}

{\displaystyle l_{ii}={\sqrt {a_{ii}-{\sum _{k=1}^{i-1}l_{ik}^{2}}}}}

a_{ij}=l_{i1}l_{j1}+\ldots +l_{ii}l_{ji}

{\displaystyle a_{ij}=l_{i1}l_{j1}+\ldots +l_{ii}l_{ji}} d'où

l_{ji}={\frac {a_{ij}-{\sum _{k=1}^{i-1}l_{ik}l_{jk}}}{l_{ii}}},\;\;\;i+1\leq j\leq n

{\displaystyle l_{ji}={\frac {a_{ij}-{\sum _{k=1}^{i-1}l_{ik}l_{jk}}}{l_{ii}}},\;\;\;i+1\leq j\leq n}

Il résulte du théorème précédent qu'il est possible de choisir tous les éléments l_ii > 0 en assurant que toutes les quantités

a_{11},\ldots ,a_{ii}-{\sum _{k=1}^{i-1}l_{ik}^{2}},\ldots

{\displaystyle a_{11},\ldots ,a_{ii}-{\sum _{k=1}^{i-1}l_{ik}^{2}},\ldots }

sont positives.

Décomposition de Cholesky alternative ou décomposition de Crout

[modifier | modifier le code ]

La décomposition de Cholesky alternative permet d'éviter l'utilisation des racines carrées au sein des sommes, source potentielle de problème en calcul numérique. Elle n'impose plus non plus l'obligation que la matrice A soit définie positive et se calcule de la façon suivante^[2]^,^[3]:

A=LDL^{\mathrm {T} }

{\displaystyle A=LDL^{\mathrm {T} }}

où D est une matrice diagonale, et L une matrice triangulaire inférieure avec des 1 sur sa diagonale.

D_{jj}=A_{jj}-\sum _{k=1}^{j-1}L_{jk}^{2}D_{kk}

{\displaystyle D_{jj}=A_{jj}-\sum _{k=1}^{j-1}L_{jk}^{2}D_{kk}}

L_{ij}={\frac {1}{D_{jj}}}\left(A_{ij}-\sum _{k=1}^{j-1}L_{ik}L_{jk}D_{kk}\right),\qquad {\text{pour }}i>j.

{\displaystyle L_{ij}={\frac {1}{D_{jj}}}\left(A_{ij}-\sum _{k=1}^{j-1}L_{ik}L_{jk}D_{kk}\right),\qquad {\text{pour }}i>j.}

Les remarque suivantes n'ont d'intérêt qu'en mathématique pure, mais pas pour la résolution de système d'équations par la méthode LDL^T :

Les factorisations LDL^T et LL^T (il faut noter que si les noms des méthodes sont ¿ proches/voisins ? , la matrice L est différente dans les deux cas) sont liées :

A=LDL^{\mathrm {T} }=LD^{\frac {1}{2}}D^{\frac {1}{2}}L^{\mathrm {T} }=LD^{\frac {1}{2}}(D^{\frac {1}{2}})^{\mathrm {T} }L^{\mathrm {T} }=LD^{\frac {1}{2}}(LD^{\frac {1}{2}})^{\mathrm {T} }

{\displaystyle A=LDL^{\mathrm {T} }=LD^{\frac {1}{2}}D^{\frac {1}{2}}L^{\mathrm {T} }=LD^{\frac {1}{2}}(D^{\frac {1}{2}})^{\mathrm {T} }L^{\mathrm {T} }=LD^{\frac {1}{2}}(LD^{\frac {1}{2}})^{\mathrm {T} }}

La dernière expression étant le produit d'une matrice triangulaire et de sa transposée, de la même manière que dans la factorisation LL^T.

On remarquera que la racine carrée d'une matrice diagonale (ici, D^1/2) se calcule trivialement en prenant les racines carrées de chacun de ses éléments.

Résolution d'un système d'équation par la méthode de Cholesky alternative

[modifier | modifier le code ]

Soit ${\textstyle b=AX}$ {\textstyle b=AX} un système d'équation linéaire à matrice symétrique. La méthode de résolution du système se décompose en quatre étapes :

Calcul des éléments des matrice L et D à l'aide des formules de la section précédente
Résolution du système intermédiaire ${\textstyle b=LY'}$ {\textstyle b=LY'} (voir méthode LU pour la résolution d'un système d'équation sous forme d'une matrice triangulaire inférieure)
Division des valeurs de Y par les coefficients de D : ${\textstyle Y'=DY\Rightarrow Y=D^{-1}Y'}$ {\textstyle Y'=DY\Rightarrow Y=D^{-1}Y'} (D étant diagonale ${\textstyle D^{-1}}$ {\textstyle D^{-1}} contient l'inverse des éléments de la diagonale de D)
Résolution du système intermédiaire ${\textstyle Y=L^{T}X}$ {\textstyle Y=L^{T}X} (voir méthode LU pour la résolution d'un système d'équation sous forme d'une matrice triangulaire supérieure)

Histoire

[modifier | modifier le code ]

La décomposition porte le nom d'André-Louis Cholesky un officier et ingénieur français. Elle figure dans le manuscrit intitulé « Sur la résolution numérique des systèmes d'équations linéaires », manuscrit porté en 2005 aux Archives de l'École Polytechnique. Daté du 2 décembre 1910, son contenu n'était auparavant connu que par une publication du commandant Benoît, qui décrivit la méthode de Cholesky en 1924, soit plusieurs années après sa mort^[4]. Il est probable que Cholesky ait découvert cette méthode en 1902^[4].

La méthode, définie pour un problème de topographie, resta longtemps inconnue des mathématiciens^[4]. Elle fut remise en avant par John Todd (en) en 1946 dans son cours d'analyse numérique au King's College de Londres ^[4].

Cette méthode est aujourd'hui centrale en analyse numérique.

Note

[modifier | modifier le code ]

↑ Xavier Gourdon, Les Maths en tête: Algèbre, page 247
↑ (en) D. Watkins, Fundamentals of Matrix Computations, John Wiley & Sons, 2002, 84 p. (ISBN 9780471213949, DOI 10.1002/0471249718 ).
↑ « Décomposition de Cholesky et de Crout », sur www.bibmath.net (consulté le 21 juillet 2024)
↑ ^{a b c et d} Claude Brezinski et Dominique Tournès, « André–Louis Cholesky (1875-1918), mathématicien, topographe, enseignant et officier », sur Images des maths , 15 juillet 2015.

Voir aussi

[modifier | modifier le code ]

Articles connexes

[modifier | modifier le code ]

Bibliographie

[modifier | modifier le code ]

La méthode de Cholesky est essentielle en analyse numérique. Il existe donc une multitude de références, parmi lesquelles :

Philippe Ciarlet, Introduction à l'analyse numérique matricielle et à l'optimisation, 1985 (rééd. 2001), éd. Masson, coll. Math. Appl. pour la Maîtrise (ISBN 2-225-68893-1)
Patrick Lascaux et Raymond Theodor, Analyse numérique matricielle appliquée à l'art de l'ingénieur - Volume 1 Méthodes directes, Dunod, 2000

Liens externes

[modifier | modifier le code ]

Sur la résolution numérique des systèmes linéaires , manuscrit de Cholesky en ligne et commenté sur le site BibNum.

v · m

Matrices

Forme

Transformée

Relation

Propriété

Famille

Associée

Résultats

Décompositions	LU QR Cholesky Schur Polaire Valeurs singulières
Diagonalisation Trigonalisation Réduction de Jordan Facteurs invariants

v · m Analyse numérique
Recherche de zéro	Méthode de Josephy-Newton Méthode de la sécante Méthode de Newton Point fixe
Transformations de matrice	Matrice de Hessenberg Décomposition LU Factorisation de Cholesky
Résolutions de systèmes	Méthode de Gauss-Seidel Méthode de surrelaxation successive (SOR) Méthode de Jacobi Décomposition QR Décomposition LU
Intégration numérique	Méthode du point médian Méthode des trapèzes Méthode de Simpson Méthodes de quadrature de Gauss Formule de Newton-Cotes Méthode de Romberg Méthode de Monte-Carlo
Équations différentielles	Méthode d'Euler (et semi-implicite) Méthode de Heun Méthodes de Runge-Kutta Méthode de Runge–Kutta symplectique d’ordre 4 (RKS4) Intégration de Verlet Leapfrog Méthode linéaire à pas multiples (Adams-Bashforth, backward differentiation formula (en))
Interpolation numérique	Courbe de Bézier Surface de Bézier Spline B-spline NURBS Interpolation polynomiale Interpolation lagrangienne Interpolation newtonienne Interpolation d'Hermite Suite de polynômes orthogonaux Interpolation bilinéaire Interpolation bicubique