LG

Soit le système linéaire $AX=b$ avec $A$ matrice carrée $N \times N$ et $b$ et $X$ vecteurs de $\mathbb R^N$

Peut se résoudre si $det A \neq 0$.

Soit on détermine $A^{-1}$, on obtient $X = A^{-1}b$.

$$A^{-1}=\frac 1 {\det A} \,^{t}\!{\rm {com} A}$$

Dans la pratique, cette méthode est trop longue et inappliquée.

$$A = LU$$

• $L$ matrice triangulaire inférieure,
• $U$ matrice triangulaire supérieure avec des 1 sur la diagonale.

On détermine $LY=b$ puis $UX=Y$.

Lorsque L est déterminée, on ne peut calculer $y_i$ que si les $y_i$ inférieurs sont calculés :

$$y_i = \frac{b_i - \sum_{j=1}^{i-1}l_{i,j}y_j}{l_{ii}}$$

Si $l_{i,i} = 0$, deux possibilités : le numérateur de $y_i$ est différent de 0 alors, il n'y a aucune solution. Si le numérateur vaut 0, alors, il y a une infinité de résultats et on peut mettre n'importe quelle valeur à $y_i$.

Complexité : $N^2/2$ multiplications et additions, $N$ divisions.

Puis on applique le même principe avec la matrice supérieure $U$ pour déterminer le vecteur $x$.

Implémentation en python :

import numpy as np
 
def fsub(A,b):
    if len(A) != len(A[0]):
        raise ValueError('La matrice n''est pas carrée.')
    if len(A) != len(b):
        raise ValueError('La matrice A et le vecteur b n''ont pas la même dimension.')
 
    N = len(A)
    x = np.empty(N)
    sol = 1
    for i in range(0, N):
        somme = b[i] - np.sum(A[i][0:i] * x[0:i])
        if somme != 0 and A[i][i] == 0:
            return 0, np.nan
        if somme == 0 and A[i][i] == 0:
            x[i] = 0
            sol = -1
        else:
            x[i] = somme / A[i][i]
 
    return sol, x
 
def bsub(A,b):
    if len(A) != len(A[0]):
        raise ValueError('La matrice n''est pas carrée.')
    if len(A) != len(b):
        raise ValueError('La matrice A et le vecteur b n''ont pas la même dimension.')
 
    N = len(A)
    x = np.empty(N)
    sol = 1
    for i in range(N - 1, -1, -1):
        if A[i][i] == 0:
            determinant_null = True
        somme = b[i] - np.sum(A[i][i+1:N] * x[i+1:N])
# Pas de solution
        if somme != 0 and A[i][i] == 0:
            return 0, np.nan
# Infinité de solutions
        if somme == 0 and A[i][i] == 0:
            x[i] = 0
            sol = -1
        else:
            x[i] = somme / A[i][i]
 
    return sol, x
 
fsub(np.array([[1, 0, 0, 0],[1, 1, 0, 0],[1, 1, 1, 0],[1, 1, 1, 1]]), np.array([1, 2, 3, 4]))
fsub(np.array([[1, 0, 0, 0],[1, 0, 0, 0],[1, 1, 1, 0],[1, 1, 1, 1]]), np.array([1, 2, 3, 4]))
fsub(np.array([[1, 0, 0, 0],[2, 0, 0, 0],[1, 1, 1, 0],[1, 1, 1, 1]]), np.array([1, 2, 3, 4]))

Valable pour une matrice carrée définie positive.

Étape 1 : $$ \begin{align*} &j = 1 \\ &i = 1 \quad &l_{jj} &= a_{jj} \\ &\forall i \in [2,N] \quad &u_{ji} &= \frac{a_{ji}}{l_{jj}} \end{align*} $$ Étape 2 : $$ \begin{align} \forall j &\in [2,N-1] : \\ i &= 1 \quad & l_{ji} &= a_{ji} \\ \forall &i \in [2, j] \quad & l_{ji} &= a_{ji} - \sum_{k=1}^{i-1}l_{jk}u_{ki} \\ \forall &i \in [j+1, N] \quad & u_{ji} &= \frac{a_{ji} - \sum_{k=1}^{j-1}l_{jk}u_{ki}}{l_{jj}} \end{align} $$ Étape 3: $$ \begin{align} \forall i &\in [1,N] : \\ j &= N \\ l_{ji} &= a_{ji} - \sum_{k=1}^{i-1}l_{jk}u_{ki} \\ \end{align} $$ Complexité : $\frac{N^3}{3}$ multiplications, $N$ divisions.

[1] §2.2 L'algorithme de Gauss.

À propos des méthodes de décomposition de type GAUSS Archive 14/​10/​2018

Si l'un des $l_{jj}$ vaut 0, il faut appliquer un pivot.

Résolution de systèmes linéaires par des méthodes directes : Gauss, LU Archive 14/10/2018

Exemple :

Université Laval avec une matrice $4*4$ Archive 14/10/2018

Valable pour une matrice carrée symétrique définie positive.

$$A = C D {}^t C$$

• $C$ matrice triangulaire inférieure ayant des 1 sur la diagonale.
• $D$ matrice diagonale

Soit $\forall i \in [1,N]$ $$ d_i = a_{ii} - \sum_{j=1}^{i-1} l_{ij}^2 d_j \\ \forall j \in [i+1,N]\\ l_{ji} = (a_{ji} - \sum_{k=1}^{i-1} l_{jk} l_{ik} d_k)/d_i $$

Exemple : Résolution de systèmes linéaires par des méthodes directes : LDL’ et Choleski en page 14 Archive 14/10/2018

Valable pour une matrice carrée symétrique définie positive.

$$A = H {}^t H$$

• $H$ matrice triangulaire inférieure.

Méthode :

The Cholesky algorithm Archive 17/10/2018

Soit $$A^i= \begin{pmatrix} I_{i-1} & 0 & 0 \\ 0 & a_{i,i} & b_i^T \\ 0 & b_i & B^i \end{pmatrix}$$

$$L_i:= \begin{pmatrix} I_{i-1} & 0 & 0 \\ 0 & \sqrt{a_{i,i}} & 0 \\ 0 & \frac{b_i}{\sqrt{a_{i,i}}} & I_{n-i} \end{pmatrix}$$

$$ A^{i+1}= \begin{pmatrix} I_{i-1} & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & B^i - \frac{b_i b_i^T}{a_{i,i}} \end{pmatrix} $$

Exemple :

$$ A = A^1= \begin{pmatrix} 4 & 12 & -16 \\ 12 & 37 & -43 \\ -16 & -43 & 98 \end{pmatrix} $$

$$ L_1= \begin{pmatrix} \sqrt4 & 0 & 0 \\ \frac{12}{\sqrt4} & 1 & 0 \\ \frac{-16}{\sqrt4} & 0 & 1 \end{pmatrix}= \begin{pmatrix} 2 & 0 & 0 \\ 6 & 1 & 0 \\ -8 & 0 & 1 \end{pmatrix} $$

$$ A^2= \begin{pmatrix} 1 & 0 & 0 \\ 0 & B^1 - \frac{b_1 b_1^T}{a_{1,1}} & \cdots \\ 0 & \vdots & \ddots \end{pmatrix} \quad B^1= \begin{pmatrix} 37 & -43 \\ -43 & 98 \end{pmatrix} \quad b_1=\begin{pmatrix} 12 \\ -16 \end{pmatrix} \quad A^2= \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 5 \\ 0 & 5 & 34 \end{pmatrix} $$

$$ L_2= \begin{pmatrix} 1 & 0 & 0 \\ 0 & \sqrt1 & 0 \\ 0 & \frac{5}{\sqrt1} & 1 \end{pmatrix}= \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 5 & 1 \end{pmatrix} $$

$$ A^3= \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & B^2 - b_2 b_2^T \end{pmatrix} \quad B^2= \begin{matrix} 34 \end{matrix} \quad b_2=\begin{matrix} 5 \end{matrix} \quad A^3= \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 9 \end{pmatrix} $$

$$ L_3= \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & \sqrt9 \end{pmatrix}= \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 3 \end{pmatrix}$$

$$ L = L_1 L_2 L_3 = \begin{pmatrix} 2 & 0 & 0 \\ 6 & 1 & 0 \\ -8 & 5 & 3 \end{pmatrix} $$

$$ \begin{align*} \forall i &\in [1,N] \\ j &= i \quad & L_{j,j} &= \sqrt{ A_{j,j} - \sum_{k=1}^{j-1} L_{j,k}^2 } \\ \forall j &\in [1,i-1] \quad & L_{i,j} &= \frac{1}{L_{j,j}} \left( A_{i,j} - \sum_{k=1}^{j-1} L_{i,k} L_{j,k} \right) \end{align*} $$ Exemple : $A = \begin{pmatrix} 4 & 12 & -16
12 & 37 & -43
-16 & -43 & 98
\end{pmatrix}$

$L = \begin{pmatrix} \sqrt 4 = 2 & 0 & 0 \\ \frac{12}{2} = 6 & \sqrt{37-6^2} = 1 & 0 \\ \frac{-16}{2} = -8 & \frac{-43-(-8)*6}{1} = 5 & \sqrt{98-8^2-5^2} = 3 \end{pmatrix}$

$QR = A$ avec R matrice diagonale supérieure et Q matrice orthogonale.

Soit $A_0 = A$

$\alpha_i = A_{i-1}[1]$ avec $A[1]$ correspondant à la première colonne de $A$ et $i \in [1,N-1]$.

$v_i = \alpha_i \pm \|\alpha_i\| e_1 $ avec $e_1$ vecteur nul sauf la première ligne vaut 1. Le signe à prendre est le même que la première ligne de $A[1]$ pour éviter une valeur nulle.

$$Q_i = I - 2 \frac{v_i v_i^t }{\|v_i\|_2^2}$$

$$A_i = \prod_{k=i-1}^{1}Q_i A$$

$$A_i= \begin{pmatrix} U_{i-1} & B \\ 0 & A_{i+1} \\ \end{pmatrix}$$ avec $U_x$ : matrice triangulaire supérieure de dimension $x$, et $b$ vecteur.

$$R = A_{N-1}$$

$$Q = \prod_{k=1}^{N-1}Q_i$$

Décomposition QR Archive du 06/02/2019

Factorisation QR et systèmes surdéterminés Archive du 06/02/2019

Exemple :

$$A=A_0= \begin{pmatrix} 12 & -51 & 4 \\ 6 & 167 & -68 \\ -4 & 24 & -41 \\ \end{pmatrix}$$

$$\alpha_1 = \begin{pmatrix} 12 \\ 6 \\ -4 \\ \end{pmatrix}$$

$$\|\alpha_1\| = 14$$

$$v_1 = \alpha_1 \pm \|\alpha_1\| e_1 = \begin{pmatrix} 26 \\ 6 \\ -4 \\ \end{pmatrix}$$

$$Q_1 = I - 2 \frac{v_1 v_1^T}{\|v_1\|_2^2} = \begin{pmatrix} -6/7 & -3/7 & 2/7 \\ -3/7 & 82/91 & 6/91 \\ 2/7 & 6/91 & 87/91 \\ \end{pmatrix}$$

$$A_1 = Q_1 A_0 = \begin{pmatrix} -14 & -21 & 14 \\ 0 & 2261/13 & -854/13 \\ 0 & 252/13 & -553/13 \\ \end{pmatrix}$$

$$\alpha_2 = \begin{pmatrix} 2261/13 \\ 252/13 \\ \end{pmatrix}$$

$$\|\alpha_2\| = 175$$

$$v_2 = \alpha_2 \pm \|\alpha_2\| e_1 = \begin{pmatrix} 4536/13 \\ 252/13 \\ \end{pmatrix}$$

$$Q_2 = I - 2 \frac{v_2 v_2^T}{\|v_2\|_2^2} = \begin{pmatrix} -323/325 & -36/325 \\ -36/325 & 323/325 \\ \end{pmatrix}$$

$$R = A_2 = Q_2 A_1 = \begin{pmatrix} -14 & -21 & 14 \\ 0 & -175 & 70 \\ 0 & 0 & -35 \\ \end{pmatrix}$$

$$Q = \begin{pmatrix} -6/7 & -3/7 & 2/7 \\ -3/7 & 82/91 & 6/91 \\ 2/7 & 6/91 & 87/91 \\ \end{pmatrix} \begin{pmatrix} 1 & 0 & 0 \\ 0 & -323/325 & -36/325 \\ 0 & -36/325 & 323/325 \\ \end{pmatrix} = \begin{pmatrix} -6/7 & 69/175 & 58/175 \\ -3/7 & -158/175 & -6/175 \\ 2/7 & -6/35 & 33/35 \\ \end{pmatrix}$$

$\forall i \in [1,N]$

$a_i$ les vecteurs colonnes de A.

$u_i = a_i - \sum_{j=1}^{i-1} (a_i \cdot e_j) e_j$

$e_i = \frac{u_i}{\|u_i\|}$

$Q = \left[e_1 | … | e_N \right]$

$ \begin{align} & \forall i > j & R_{i,j} = 0 \\ & \forall i \le j & R_{i,j} = a_j e_i \\ \end{align} $

QR Decomposition with Gram-Schmidt Archive du 16/02/2019

Exemple :

$$A= \begin{pmatrix} 12 & -51 & 4 \\ 6 & 167 & -68 \\ -4 & 24 & -41 \\ \end{pmatrix}$$

$$u_1 = \begin{pmatrix} 12 \\ 6 \\ -4 \\ \end{pmatrix}$$

$$\|u_1\| = 14$$

$$e_1 = \frac{u_1}{\|u_1\|} = \begin{pmatrix} 6/7 \\ 3/7 \\ -2/7 \\ \end{pmatrix}$$

$$a_2 \cdot e_1 = (-51)*6/7 + 167 * 3/7 + 24 * (-2) / 7 = 21$$

$$u_2 = a_2 - (a_2 \cdot e_1) e_1 = \begin{pmatrix} -51 - 21 * 6 / 7 \\ 167 - 21 * 3 / 7 \\ 24 - 21 * (-2) / 7 \\ \end{pmatrix} = \begin{pmatrix} -69 \\ 158 \\ 30 \\ \end{pmatrix} $$

$$\|u_2\| = 175$$

$$e_2 = \frac{u_2}{\|u_2\|}$$

$$a_3 \cdot e_1 = 4 * 6 / 7 + (-68) * 3 / 7 + (-41) * (-2) / 7 = -14$$

$$a_3 \cdot e_2 = 4 * (-69) / 175 + (-68) * 158 / 175 + (-41) * 30 / 175 = -70$$

$$u_3 = a_3 - (a_3 \cdot e_1) e_1 - (a_3 \cdot e_2) e_2 = \begin{pmatrix} 4 - (-14) * 6 / 7 - (-70) * (-69) / 175 \\ -68 - (-14) * 3 / 7 - (-70) * 158 / 175 \\ -41 - (-14) * (-2) / 7 - (-70) * 30 / 175 \\ \end{pmatrix} = \begin{pmatrix} -58 / 5 \\ 6/5 \\ -33 \\ \end{pmatrix} $$

$$\|u_3\| = 35$$

$$e_3 = \frac{u_3}{\|u_3\|}$$

$$Q = \begin{pmatrix} 6/7 & -69/175 & -58/175 \\ 3/7 & 158/175 & 6/175 \\ -2/7 & 6/ 35 & -33/ 35 \\ \end{pmatrix} $$

$$R = \begin{pmatrix} 14 & 21 & -14 \\ 0 & 175 & -70 \\ 0 & 0 & 35 \\ \end{pmatrix} $$

$\forall j \in [1,N-1]$ avec j désignant les colonnes.

$\forall i \in [j+1,N]$ par ordre décroissant avec désignant les lignes.

On cherche donc à annuler chaque terme sur le triangle inférieur.

Au départ $R_{1,N}$ vaut A.

$r_{i,j} = \sqrt{\left ( R^{ij}_{i,j} \right ) ^ 2 + \left ( R^{ij}_{i-1,j} \right ) ^ 2}$

On calcule la matrice de rotation $G^{i,j}$ : vaut la matrice identité sauf : $$G_{i,i} = \frac{R_{i-1,j}}{r_{i,j}}$$ $$G_{i-1,i-1} = \frac{R_{i-1,j}}{r_{i,j}}$$ $$G_{i,i-1} = -\frac{R_{i,j}}{r_{i,j}}$$ $$G_{i-1,i} = \frac{R_{i,j}}{r_{i,j}}$$

Le $R^{(n+1)}$ suivant vaut $G*R^{(n)}$.

An Example of QR Decomposition Archive du 17/02/2019

The QR Reduction Archive du 17/02/2019

Exemple :

$$R^{(0)} = A= \begin{pmatrix} 12 & -51 & 4 \\ 6 & 167 & -68 \\ -4 & 24 & -41 \\ \end{pmatrix}$$

$i=3, j=1$

$r_{3,1} = \sqrt{(-4)^2+6^2} = 2 \sqrt{13}$

$$G^{(1)} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \frac{G^{(0)}_{i-1,j}}{r_{3,1}} & \frac{G^{(0)_{i,j}}}{r_{3,1}} \\ 0 & -\frac{G^{(0)}_{i,j}}{r_{3,1}} & \frac{G^{(0)}_{i-1,j}}{r_{3,1}} \\ \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \frac{6}{2 \sqrt{13}} & \frac{-4}{2 \sqrt{13}} \\ 0 & -\frac{-4}{2 \sqrt{13}} & \frac{6}{2 \sqrt{13}} \\ \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \frac{3}{\sqrt{13}} & \frac{-2}{\sqrt{13}} \\ 0 & \frac{2}{\sqrt{13}} & \frac{3}{\sqrt{13}} \\ \end{pmatrix} $$

$$R^{(1)} = G^{(1)} R^{(0)} = \begin{pmatrix} 12 & -51 & 4 \\ 2 \sqrt{13} & 453/\sqrt{13} & -122/\sqrt{13} \\ 0 & 406/\sqrt{13} & -259/\sqrt{13} \\ \end{pmatrix} $$

$i=2, j=1$

$r_{2,1} = \sqrt{(12)^2+(2 \sqrt{13})^2} = 14$

$$G^{(2)} = \begin{pmatrix} \frac{G^{(0)}_{i-1,j}}{r_{2,1}} & \frac{G^{(0)_{i,j}}}{r_{2,1}} & 0 \\ -\frac{G^{(0)}_{i,j}}{r_{2,1}} & \frac{G^{(0)}_{i-1,j}}{r_{2,1}} & 0 \\ 0 & 0 & 1 \\ \end{pmatrix} = \begin{pmatrix} \frac{12}{14} & \frac{2 \sqrt{13}}{14} & 0 \\ -\frac{2 \sqrt{13}}{14} & \frac{12}{14} & 0 \\ 0 & 0 & 1 \\ \end{pmatrix} = \begin{pmatrix} \frac{6}{7} & \frac{\sqrt{13}}{7} & 0 \\ -\frac{\sqrt{13}}{7} & \frac{6}{7} & 0 \\ 0 & 0 & 1 \\ \end{pmatrix} $$

$$R^{(2)} = G^{(2)} R^{(1)} = \begin{pmatrix} 14 & 21 & -14 \\ 0 & 483/\sqrt{13} & -122/\sqrt{13} \\ 0 & 406/\sqrt{13} & -259/\sqrt{13} \\ \end{pmatrix} $$

$i=3, j=2$

$r_{3,2} = \sqrt{(483/\sqrt{13})^2+(406/\sqrt{13})^2} = 175$

$$G^{(3)} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \frac{G^{(0)}_{i-1,j}}{r_{3,2}} & \frac{G^{(0)_{i,j}}}{r_{3,2}} \\ 0 & -\frac{G^{(0)}_{i,j}}{r_{3,2}} & \frac{G^{(0)}_{i-1,j}}{r_{3,2}} \\ \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \frac{483/\sqrt{13}}{175} & \frac{406/\sqrt{13}}{175} \\ 0 & -\frac{406/\sqrt{13}}{175} & \frac{483/\sqrt{13}}{175} \\ \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \frac{69}{25\sqrt{13}} & \frac{58}{25\sqrt{13}} \\ 0 & -\frac{58}{25\sqrt{13}} & \frac{69}{25\sqrt{13}} \\ \end{pmatrix} $$

$$R^{(2)} = G^{(2)} R^{(1)} = \begin{pmatrix} 14 & 21 & -14 \\ 0 & 175 & -70 \\ 0 & 0 & -35 \\ \end{pmatrix} $$

[1] : DESTUYNDER, Philippe. Méthodes numériques pour l'ingénieur. Lavoisier, 2010, 248 pages.

[2] : Timothy, Vismor. Matrix Algorithms. January 30, 2015, 69 pages. Archive 11/11/2018