L'idée de la décomposition d'une matrice $\boldsymbol{A}$ comme le produit de deux matrices $ \boldsymbol{L}$ et $\boldsymbol{U}$ remonte à \citet{Turing48}\sindex[aut]{Turing A.}. Elle a été anticipée sur un cas particulier par la méthode de Cholesky (voir \citet{CommandantBenoit24}\sindex[aut]{Commandant Benoit}). Elle consiste à décomposer une matrice en deux sous matrices triangulaires, la première ($\boldsymbol{L}$ comme lower) étant triangulaire inférieure et la seconde ($\boldsymbol{U}$ comme upper) étant triangulaire supérieure. Par exemple, pour une matrice $4\times3$, on aura
\[ \boldsymbol{A}=\begin{bmatrix} a_{11} &a_{12}&a_{13}& a_{14}\\ a_{21} &a_{22}&a_{23}& a_{24}\\ a_{31} &a_{32}&a_{33}& a_{34}\\ a_{41} &a_{42}&a_{43}& a_{44} \end{bmatrix}=\begin{bmatrix} 1 &0&0& 0\\ \alpha_{21} &1&0& 0\\ \alpha_{31} &\alpha_{32}&1& 0\\ \alpha_{41} &\alpha_{42}&\alpha_{43}& 1 \end{bmatrix} \begin{bmatrix} \beta_{11} &\beta_{12}&\beta_{13}& \beta_{14}\\ 0 &\beta_{22}&a_{23}& \beta_{24}\\ 0 &0&\beta_{33}& \beta_{34}\\ 0 &0&0& \beta_{44} \end{bmatrix}=\boldsymbol{L}\boldsymbol{U} \]Il n'est pas toujours possible de décomposer directement la matrice $\boldsymbol{A}$ sous la forme $\boldsymbol{L}\boldsymbol{U}$. C'est, par exemple, le cas de la matrice :
\[ \boldsymbol{A}=\begin{bmatrix} 0 &a_{12}&a_{13}& a_{14}\\ a_{21} &a_{22}&a_{23}& a_{24}\\ a_{31} &a_{32}&a_{33}& a_{34}\\ a_{41} &a_{42}&a_{43}& a_{44} \end{bmatrix} \]simplement du fait de la présence du zéro. Cependant, il est possible, par permutation, de transformer $\boldsymbol{A}$ de telle sorte que l'opération puisse être validée. Pour ce faire, on introduit une matrice de permutation $\boldsymbol{P}$ telle que :
\[ \boldsymbol{P}^{-1} \boldsymbol{A} = \boldsymbol{L} \boldsymbol{U}\hspace{.25cm}\Longleftrightarrow\hspace{.25cm} \boldsymbol{A} = \boldsymbol{P}^{-1}\boldsymbol{L}\boldsymbol{U} \]où $\boldsymbol{P}$ est une matrice qui échange les lignes d'une autre matrice. Par exemple :
\[ \boldsymbol{P} = \begin{bmatrix} 0 & 1 & 0 \\ 1 & 0 & 0\\ 0 & 0 & 1 \end{bmatrix} \]est une matrice de permutation\footnote{Précisément, une matrice de permutation est composée de $0$ et de $1$ ayant un unique $1$ par ligne et par colonne.}.
Si on applique cette matrice à une matrice $3\times 3$ quelconque, par exemple la matrice
\[ \boldsymbol{A} = \begin{bmatrix} a & a & a \\ b & b & b\\ c & c & c \end{bmatrix} \]on obtient
\[ \boldsymbol{P}\boldsymbol{A} = \begin{bmatrix} 0 & 1 & 0 \\ 1 & 0 & 0\\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} a & a & a \\ b & b & b\\ c & c & c \end{bmatrix}= \begin{bmatrix} b & b & b \\ a & a & a\\ c & c & c \end{bmatrix} \]L'intérêt de la décomposition\footnote{On parle aussi de factorisation.} $\boldsymbol{L} \boldsymbol{U}$ est encore une fois de diminuer la complexité liée à l'inversion d'une matrice nécessaire pour résoudre un système d'équations. Supposons que $\boldsymbol{P} = \boldsymbol{I}$, alors, si le système :
\[ \boldsymbol{A}\boldsymbol{x}= \boldsymbol{b} \]est décomposé, on a :
\[ \boldsymbol{A}\boldsymbol{x}= \boldsymbol{L}\boldsymbol{U}\boldsymbol{x}= \boldsymbol{b} \]on peut, dans un premier temps, le résoudre en introduisant le vecteur $\boldsymbol{y} =\boldsymbol{U}\boldsymbol{x}$, ce qui donne :
\[ \boldsymbol{L}\boldsymbol{y} = \boldsymbol{x} \]nouveau système qui est solvable de manière récursive vers le bas (on parle de descente) étant donné la structure de $\boldsymbol{L}$ ce qui est moins coûteux que d'inverser $\boldsymbol{L}$. Enfin, de $\boldsymbol{y}$, il est élémentaire de trouver $\boldsymbol{x}$ par des opérations du même type. On notera que cette seconde étape se fait par remontée, étant donnée la structure de $\boldsymbol{U}$. Objectivement, la décomposition $\boldsymbol{L} \boldsymbol{U}$ n'est autre qu'une forme de l'élimination gaussienne. Il demeure néanmoins à exposer un algorithme qui permet de l'obtenir directement. Voici une méthode qui évite de réaliser un pivotage.
On écrit les matrices comme des matrices par block, i.e. :
\[ \boldsymbol{A} = \left[\begin{array}{c|c}a_{11} & \boldsymbol{A}^\top_{12}\\\hline \boldsymbol{A}_{21} & \boldsymbol{A}_{22}\end{array} \right], \hspace{.25cm}\boldsymbol{L} = \left[\begin{array}{c|c}1 & \boldsymbol{0}_{n}^\top\\\hline \boldsymbol{L}_{21} & \boldsymbol{L}_{22}\end{array} \right] \hspace{.25cm}\text{et}\hspace{.25cm}\boldsymbol{U} = \left[\begin{array}{c|c}u_{11} & \boldsymbol{U}^\top_{12}\\\hline \boldsymbol{0}_{m} & \boldsymbol{U}_{22}\end{array} \right] \]où $a_{11}$ sont des scalaires et les $\boldsymbol{0}_k$ ---~$k= m \text{ ou }n$ des vecteurs de $0$. On détermine $\boldsymbol{L}$ et $\boldsymbol{U}$ de l'égalité $\boldsymbol{A} = \boldsymbol{L}\boldsymbol{U}$, c'est-à -dire de :
\[ \begin{bmatrix} a_{11} & \boldsymbol{A}^\top_{12}\\ \boldsymbol{A}_{21} & \boldsymbol{A}_{22} \end{bmatrix}= \begin{bmatrix} 1 & \boldsymbol{0}^\top_n\\ \boldsymbol{L}_{21} & \boldsymbol{L}_{22} \end{bmatrix}\begin{bmatrix} u_{11} & \boldsymbol{U}^\top_{12}\\ \boldsymbol{0}_{m} & \boldsymbol{U}_{22} \end{bmatrix}= \begin{bmatrix} u_{11} & \boldsymbol{U}^\top_{12}\\ u_{11} \boldsymbol{L}_{21}&\boldsymbol{L}_{21}\boldsymbol{U}^\top_{12}+\boldsymbol{L}_{22}\boldsymbol{U}_{22} \end{bmatrix} \]On commence par déterminer tout d'abord les lignes de $\boldsymbol{U}$ et les colonnes de $\boldsymbol{L}$ par identification : \[ u_{11} = a_{11}, \hspace{.25cm}\boldsymbol{U}_{12} = \boldsymbol{A}_{12}\hspace{.25cm}\text{et}\hspace{.25cm}\boldsymbol{L}_{21} = \frac{\boldsymbol{A}_{21}}{a_{11}} \]
Il ne reste plus qu'à factoriser la matrice $(n-1) \times (n-1)$ :
\[ \boldsymbol{A}_{22}- \boldsymbol{L}_{21}\boldsymbol{U}^\top_{12}=\boldsymbol{A}_{22}- \frac{\boldsymbol{A}^\top_{21}}{a_{11}}\boldsymbol{A}_{12}=\boldsymbol{L}_{22}\boldsymbol{U}_{22} \]Un exemple permet d'apprécier la démarche. Considérons la matrice :
\[ \boldsymbol{A} = \left[\begin{array}{c|ccc} 3 & 2 & 5 & 1\\\hline 2 & 4 & 6 & 2\\ 1 & 1 & 2 & 0\\ 0 & 0 & 1 & 1 \end{array}\right]= \left[\begin{array}{c|ccc} 1 & 0 & 0 &0\\\hline l_{21} & 1 & 0 & 0\\ l_{31} & l_{32} & 1 & 0\\ l_{41} & l_{42} & l_{43} & 1\\ \end{array}\right]\left[\begin{array}{c|ccc} u_{11} & u_{12} & u_{13} &u_{14}\\\hline 0 & u_{22} & u_{23} & u_{24}\\ 0 & 0 & u_{33} & u_{34}\\ 0 & 0 & 0 & u_{44}\\ \end{array}\right] \]On applique les différentes formules pour obtenir :
\[ u_{11} = a_{11}=3, \hspace{.25cm} \boldsymbol{U}^\top_{12} = \boldsymbol{A}^\top_{12}= \begin{bmatrix} 2 & 5 & 1 \end{bmatrix}\hspace{.25cm}\text{et}\hspace{.25cm} \boldsymbol{L}_{21} = \frac{\boldsymbol{A}_{21}}{a_{11}}=\begin{bmatrix} 2/3\\1/3\\0 \end{bmatrix} \]ce qui signifie que :
\[ \boldsymbol{A} = \left[\begin{array}{c|ccc} 3 & 2 & 5 & 1\\\hline 2 & 4 & 6 & 2\\ 1 & 1 & 2 & 0\\ 0 & 0 & 1 & 1 \end{array}\right]= \left[\begin{array}{c|ccc} 1 & 0 & 0 &0\\\hline 2/3 & 1 & 0 & 0\\ 1/3 & l_{32} & 1 & 0\\ 0 & l_{42} & l_{43} & 1\\ \end{array}\right]\left[\begin{array}{c|ccc} 3 & 2 & 5 &1\\\hline 0 & u_{22} & u_{23} & u_{24}\\ 0 & 0 & u_{33} & u_{34}\\ 0 & 0 & 0 & u_{44}\\ \end{array}\right] \]Pour finir
\begin{eqnarray*} \boldsymbol{A}_{22}- \frac{\boldsymbol{A}^\top_{21}}{a_{11}}\boldsymbol{A}_{12}&=& \begin{bmatrix} 4 & 6 & 2\\ 1 & 2 & 0\\ 0 & 1 & 1 \end{bmatrix}- \begin{bmatrix} 2/3\\ 1/3\\ 0 \end{bmatrix} \begin{bmatrix} 2 & 5 & 1 \end{bmatrix} = \begin{bmatrix} 4 & 6 & 2\\ 1 & 2 & 0\\ 0 & 1 & 1 \end{bmatrix}- \begin{bmatrix} 4/3 & 10/3 & 2/3\\ 2/3 & 5/3 & 0\\ 0 & 0 & 0 \end{bmatrix}\\ &=&\begin{bmatrix} 8/3 & 10/3 & 8/3\\ 2/3 & 0 & 2/3\\ 2/3 & 5/3 & 0 \end{bmatrix}=\begin{bmatrix} 1 & 0 & 0\\ l_{32} & 1 & 0\\ l_{42} & l_{43} & 1 \end{bmatrix}\begin{bmatrix} u_{22} & u_{23} & u_{24}\\ 0 & u_{33} & u_{34}\\ 0 & 0 & u_{44} \end{bmatrix}\\&=& \begin{bmatrix} u_{22} & u_{23} & u_{24}\\ l_{32}u_{22} & l_{32}u_{23}+u_{33} & l_{32}u_{24}+u_{34}\\ l_{42}u_{22} & l_{42}u_{23}+ l_{43}u_{33}& l_{42}u_{34}+l_{43}u_{34}+u_{44} \end{bmatrix}\\ &=& \begin{bmatrix} 8/3 & 10/3 & 8/3\\ l_{32}(8/3) & l_{32}(10/3)+u_{33} & l_{32}(8/3)+u_{34}\\ l_{42}u_{22} & l_{42}u_{23}+ l_{43}u_{33}& l_{42}u_{34}+l_{43}u_{34}+u_{44} \end{bmatrix} \end{eqnarray*}La seconde ligne donne naissance à trois égalités. La première $2/3 = l_{32}(8/3)$ donne $l_{32} =1/4$. Après substitution dans les deux autres équations, on a :
\[ \begin{array}{ccc} 0 = (1/4)(10/3) +u_{33} &&u_{33} = -(5/6)\\ & \Longrightarrow&\\ 2/3=(1/4)(8/3)+u_{34}&&u_{34}= -(2/3) \end{array} \]ce qui montre que la méthode est récursive et enfin, comme $0 = l_{42}(8/3)$, ce qui implique que $l_{42} = 0$ et :
\[ \begin{array}{ccc} 5/3 = -l_{43}(5/6) &&l_{43} = -2\\ & \Longrightarrow&\\ 0=2(2/3)+u_{44}&&u_{44}= -(4/3) \end{array} \]Ainsi :
\[ \boldsymbol{L}= \begin{bmatrix} 1 & 0 & 0 &0\\ 2/3 & 1 & 0 & 0\\ 1/3 & 1/4 & 1 & 0\\ 0 & 0 & -2 & 1\\ \end{bmatrix}\hspace{.25cm}\text{et}\hspace{.25cm}\boldsymbol{U}=\begin{bmatrix} 3 & 2 & 5 &1\\ 0 & 8/3 & 10/3 & 8/3\\ 0 & 0 & -(5/6) & -(2/3)\\ 0 & 0 & 0 & -(4/3)\\ \end{bmatrix} \]Un dernier point : la décomposition $\boldsymbol{L}\boldsymbol{U}$ des matrices n'est pas nécessairement unique.
$\texttt{SageMath}$ a une commande qui permet d'obtenir la décomposition $\boldsymbol{L}\boldsymbol{U}$ d'une matrice. Comme cette commande possède des options, il est nécessaire d'être précis.
| $\mathtt{option1} :$ | stratégie de pivotement |
| $\mathtt{auto}$ (default) - regarde si les entrées de la matrice sont ordonnées auquel cas on adopte une stratégie de pivot partiel. Dans le cas contraire, on revient à une stratégie non zéro. C'est le meilleur choix pour les algorithmes généraux qui peuvent appeler des matrices dont les entrées sont de types variés. | |
| $\mathtt{partial}$ - recherche dans chaque colonne l'élément qui possède la plus grande valeur absolue et on permute la ligne contenant cet élément. | |
| $\mathtt{nonzero}$ - trouve le premier élément non nul dans une colonne et on utilise la ligne correspondant à cet élément. | |
| $\mathtt{option2} :$ | format |
| $\mathtt{plu}$ (default) - un triplet de matrices $\boldsymbol{P}$, $\boldsymbol{L}$ et $\boldsymbol{U}$ tels que $\boldsymbol{A} = \boldsymbol{P}*\boldsymbol{L}*\boldsymbol{U}$. | |
| $\mathtt{compact}$ - une paire ; la permutation des lignes donnée comme un tuple et les matrices $\boldsymbol{L}$ and $\boldsymbol{U}$ combined en une seule matrice. |
On notera que les matrices peuvent être différentes selon la méthode retenue comme le montrent les exemples d'application suivants (empruntés à la documentation). On commence ave l'option $\mathtt{partial}$ (ne pas oublier les guillemets pour passer une chaîne de caractères).
Puis, on continue avec l'option $\mathtt{nonzero}$
Puis, on continue avec l'option $\mathtt{compact}$
On note que le format $\mathtt{compact}$ n'est pas très utile pour la suite du cours. Néanmoins, si on veut comprendre le lien entre les formats, on peut consulter la doc à la page Base class for matrices, part 2
Donner une décomposition $\boldsymbol{L}\boldsymbol{U}$ des matrices suivantes :