---
title: VISI601_CMI (4c) Algèbre linéaire numérique (III)
type: slide
slideOptions:
transition: slide
progress: true
slideNumber: true
---
# Algèbre linéaire numérique (III)
> [name=Jacques-Olivier Lachaud][time=January 2022][color=#907bf7]
> (Les images peuvent être soumises à des droits d'auteur. Elles sont utilisées ici exclusivement dans un but pédagogique)
###### tags: `visi601`
Retour à [VISI601_CMI (Main) Algorithmique numérique](https://codimd.math.cnrs.fr/s/IWTaBkA9m)
---
## Décomposition en matrices triangulaires
On résoud facilement les systèmes triangulaires $\Rightarrow$ on factorise $A$ avec ces matrices
:::info
0. Elimination de Gauss ou Gauss-Jordan
(principe)
2. Décomposition $A=LU$
(didactique, pas utilisé en pratique, instable numériquement)
3. Décomposition $A=PLU$
(bonne façon de résoudre $A\mathbf{x}=\mathbf{b}$ si $A$ matrice pleine, non symétrique)
5. Decomposition $A=LL^T$ ou $A=LDL^T$ : factorisation de Cholesky
(bonne façon de résoudre $A\mathbf{x}=\mathbf{b}$ si $A$ matrice pleine et Sym.Def.Pos.)
:::
---
## Elimination de Gauss ou Gauss-Jordan
(Chine, vers 180, Newton vers 1670, puis Gauss)
* Transformer une matrice en une matrice en forme échelonnée réduite par:
* Échange de deux lignes
* Multiplication d'une ligne par un scalaire non nul
* Ajout du multiple d'une ligne à une autre ligne
* Calcul du rang, du déterminant, solution $Ax=b$, de l'inverse $A^{-1}$
$$
[{\mathrm {A}}|{\mathrm {I}}]=\left({\begin{array}{rrr|rrr}2&-1&0&1&0&0\\-1&2&-1&0&1&0\\0&-1&2&0&0&1\end{array}}\right)\quad
[{\mathrm {I}}|{\mathrm {B}}]=\left({\begin{array}{rrr|rrr}1&0&0&{\frac {3}{4}}&{\frac {1}{2}}&{\frac {1}{4}}\\[3pt]0&1&0&{\frac {1}{2}}&1&{\frac {1}{2}}\\[3pt]0&0&1&{\frac {1}{4}}&{\frac {1}{2}}&{\frac {3}{4}}\end{array}}\right)
$$
[Voir Wikipedia](https://fr.wikipedia.org/wiki/Élimination_de_Gauss-Jordan)
---
## Décomposition A=LU
$$
A=LU=\begin{bmatrix}
1& 0 & \cdots & 0\\
l_{21}&1& \cdots & 0\\
\vdots&\vdots & \ddots & \vdots\\
l_{n1}&l_{n2}& \cdots & 1\\
\end{bmatrix}
\begin{bmatrix}
u_{11}& u_{12} & \cdots & u_{1n}\\
0 & u_{22} & \cdots & u_{2n}\\
\vdots&\vdots & \ddots & \vdots\\
0 & 0 & \cdots & u_{nn}\\
\end{bmatrix}
$$
par élimination de Gauss successives
Résoudre $A\mathbf{x}=\mathbf{b}$ en 3 étapes:
1. Trouvez $A=LU$ ==environ $\frac{2}{3}n^3$ flops==
2. Résoudre $L\mathbf{y}=\mathbf{b}$ (facile, sys. tri., $\mathbf{y}$ inconnu) ==environ $n^2$ flops==
3. Résoudre $U\mathbf{x}=\mathbf{y}$ (facile, sys. tri., $\mathbf{x}$ inconnu) ==environ $n^2$ flops==
---
$$
\underbrace{\begin{bmatrix}
1 & 0 & \cdots & 0\\
-\frac{a_{21}}{a_{11}} & 1 & \cdots & 0\\
\vdots&\vdots & \ddots & \vdots\\
-\frac{a_{n1}}{a_{11}} & 0 & \cdots & 1\\
\end{bmatrix}}_{M_1}
\underbrace{\begin{bmatrix}
a_{11}&a_{12}& \cdots & a_{1n}\\
a_{21}&a_{22}& \cdots & a_{2n}\\
\vdots&\vdots & \ddots & \vdots\\
a_{m1}&a_{m2}& \cdots & a_{mn}\\
\end{bmatrix}}_{A^{(1)}=A}=
\underbrace{\begin{bmatrix}
a_{11}&a_{12}& \cdots & a_{1n}\\
0 &a^{(2)}_{22}& \cdots & a^{(2)}_{2n}\\
\vdots&\vdots & \ddots & \vdots\\
0 &a^{(2)}_{m2}& \cdots & a^{(2)}_{mn}\\
\end{bmatrix}}_{A^{(2)}}
$$
On forme successivement des matrices $M_i$ qui modifient une partie $n-i \times n-i$ de $A_i$
$M_n \cdots M_1 A = A^{(n)} = U$ est triangulaire supérieure
Et $(M_n \cdots M_1)^{-1}=M_1^{-1} \cdots M_n^{-1}=\begin{bmatrix}
1 & 0 & \cdots & 0\\
\frac{a_{21}}{a_{11}} & 1 & \cdots & 0\\
\vdots&\vdots & \ddots & \vdots\\
\frac{a_{n1}}{a_{11}} & \frac{a^{(2)}_{n2}}{a^{(2)}_{22}} & \cdots & 1\\
\end{bmatrix}$ est triangulaire inférieure
---
## Décomposition A=PLU ou PA=LU
* $A$ peut être inversible et ne pas avoir de décomposition $LU$
* instabilité numérique de $LU$ si les pivots $a^{(k)}_{kk}$ proches de 0
* Matrices de permutation $P_k$ où on met le + grand pivot sur la diagonale
${\displaystyle P_{1}A=\left({\begin{array}{c|ccc}a&&w^{\textsf {T}}&\\\hline &&&\\v&&A'&\\&&&\end{array}}\right)}=
\left({\begin{array}{c|ccc}1&&0&\\\hline &&&\\ \frac{1}{a}v&&I_{n-1}&\\&&&\end{array}}\right)\left({\begin{array}{c|c}a&w^{\textsf {T}}\\\hline &\\0&A'-\frac{1}{a}vw^{\textsf {T}}\\&\end{array}}\right)$
Récursivement, si $A^{(2)}=A'-\frac{1}{a}vw^{\textsf {T}}$, on résoud $P'A^{(2)}=L'U'$ et
${\displaystyle \left({\begin{array}{c|ccc}1&&0&\\\hline &&&\\0&&P'&\\&&&\end{array}}\right)P_{1}A=\left({\begin{array}{c|ccc}1&&0&\\\hline &&&\\ \frac{1}{a}P'v&&L'&\\&&&\end{array}}\right)\left({\begin{array}{c|ccc}a&&w^{\textsf {T}}&\\\hline &&&\\0&&U'&\\&&&\end{array}}\right),}$
---
## Factorisation LLT de Cholesky
:::success
Si $A$ est sym. déf. pos., alors il existe une matrice tri. inf. $L$,
de termes diagonaux $>0$, telle que $A=LL^T$.
:::
$$
l_{11} =\sqrt{a_{11}} \quad\text{et pour $j=2,\ldots, n$ :} \quad l_{j1} = a_{1j}/h_{11}
$$
et pour tout $i=2, \ldots, n$ :
$$
l_{ii} = \left( a_{ii} - \sum_{k=1}^{i-1} l^2_{ik} \right)^{\frac{1}{2}}, \quad \text{et pour $j=i+1, \ldots, n$ : }\quad l_{ji} = \left( a_{ij} - \sum_{k=1}^{i-1} l_{ik} l_{jk} \right) / l_{ii}.
$$
---
## Factorisation LDLT de Cholesky
:::success
Variante Cholesky $A=LL^T=(L'D^{1/2})(L'D^{1/2})^T=L'DL'^T$, $D$ diagonale,
avec $L'$ avec des 1 sur la diagonale. (Evite les calculs de $\sqrt{\cdot}$ )
:::
$$
{\displaystyle D_{j}=A_{jj}-\sum _{k=1}^{j-1}L_{jk}^{2}D_{k},}\quad
{\displaystyle L_{ij}={\frac {1}{D_{j}}}\left(A_{ij}-\sum _{k=1}^{j-1}L_{ik}L_{jk}D_{k}\right)\quad {\text{for }}i>j.}
$$
---
## Stabilité numérique
* Inverse **exact** connu $A^{-1}$ puis multiplication:
$\widehat{\mathbf{x}}=A^{-1}\mathbf{b}$, erreur $\| \mathbf{x}-\widehat{\mathbf{x}}\| \le K(A) \| \mathbf{x} \| O(\frak{u})$
* Encore pire si on calcule l'inverse $\widehat{A^{-1}}$
* décomposition PLU, trouve $\widehat{\mathbf{x}}$ tel que $\widehat{\mathbf{b}}=A\widehat{\mathbf{x}}$ et :
erreur amont $\| \mathbf{b} - A \mathbf{\widehat{x}} \| \le (\| A \| +\| U \|)\| \mathbf{x} \| O(\frak{u})$
==toujours plus favorable==
* factorisation de Cholesky:
erreur amont $\| \mathbf{b} - A \mathbf{\widehat{x}} \| \le 8n(n+1)\| A \|_2 \| \mathbf{x} \| \frak{u}$
---
## Considérations algorithmiques
[Wikipedia Cholesky](https://en.wikipedia.org/wiki/Cholesky_decomposition)
* algorithme Cholesky–Banachiewicz calcule ligne par ligne (en lisant $A$ par ligne)
* algorithme Cholesky–Crout calcule colonne par colonne (en lisant $A$ par colonne)
* Cholesky $\approx \frac{1}{3}n^3$ flops vs PLU $\approx \frac{2}{3}n^3$ flops