Soit \(A\) la matrice \[ A = \begin{pmatrix} 2 & -1 & 0\\ 4 & -1 & 2\\ -6& 2 & 0 \end{pmatrix} \]
La décomposition \(LU\) de la matrice \(A\) est de la forme \[ \begin{pmatrix} 2 & -1 & 0\\ 4 & -1 & 2\\ -6& 2 & 0 \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0\\ L_{21} & 1 & 0\\ L_{31} & L_{32} & 1 \end{pmatrix} \begin{pmatrix} U_{11} & U_{12} & U_{13}\\ 0 & U_{22} & U_{23}\\ 0 & 0 & U_{33} \end{pmatrix} \] Pour \(j=1\) à \(n: U_{1j}=a_{1j}\), i.e., \(U_{11}=2, U_{12}=-1\) et \(U_{13}=0\).
1ère colonne de \(L\) : \(L_{21}=\frac{a_{21}}{U_{11}}=\frac{4}{2}=2\) et \(L_{31}=\frac{a_{31}}{U_{11}}=\frac{-6}{2}=-3\).
2ème ligne de \(U\) : \(U_{22}=a_{22}-L_{21}U_{12}=-1-2(-1)=1\) et \(U_{23}=a_{23}-L_{21}U_{13}=2-2(0)=2\).
2ème colonne de \(L\) : \(L_{32}=\frac{1}{U_{22}}(a_{32}-L_{31}U_{12})=(2-(-3)(-1))/1=-1\).
3ème ligne de \(U\) : \(U_{33}=a_{33}-(L_{31}U_{13}+L_{32}U_{23})=0-((-3)(0)+(-1)2)=2\).
La décomposition \(LU\) de \(A\) est donc \[ L= \begin{pmatrix} 1 & 0 & 0\\ 2 & 1 & 0\\ -3 & -1 & 1 \end{pmatrix},\quad U= \begin{pmatrix} 2 & -1 & 0\\ 0 & 1 & 2\\ 0 & 0 & 2 \end{pmatrix}. \] Comme \(A=LU\) et que le déterminant d’une matrice triangulaire (supérieure ou inférieure) est égale au produit des éléments diagonaux, on a \[ \det (A)=\det(LU)=\det(L)\det(U)=(1\cdot 1 \cdot 1) (2\cdot 1 \cdot 2)=4 \]
Comme \(A=LU\), la résolution du système revient alors à résoudre deux systèmes triangulaires, i.e, \[\begin{align*} A\vec{x}=\vec{b} \Leftrightarrow LU\vec{x}=\vec{b} \Leftrightarrow \begin{cases} L\vec{y}=\vec{b}\\ U\vec{x}=\vec{y} \end{cases} \end{align*}\] Résolution de \(L\vec{y}=\vec{b}\) (descente) : \[ \begin{pmatrix} 1 & 0 & 0\\ 2 & 1 & 0\\ -3 & -1 & 1 \end{pmatrix} \begin{pmatrix} y_1\\ y_2\\ y_3 \end{pmatrix} = \begin{pmatrix} 2\\ 14\\ 12 \end{pmatrix} \] \(y1=2, y_2=14-2y_1=10\) et \(y_3=12+3y_1+y_2=28\).
Résolution de \(U\vec{x}=\vec{y}\) (remontée) : \[ \begin{pmatrix} 2 & -1 & 0\\ 0 & 1 & 2\\ 0 & 0 & 2 \end{pmatrix} \begin{pmatrix} x_1\\ x_2\\ x_3 \end{pmatrix} = \begin{pmatrix} 2\\ 10\\ 28 \end{pmatrix} \] \(x_1=-8, x_2=-18\) et \(x_3=14\).
\[\begin{align*} A^2\vec{x}=\vec{b} \Leftrightarrow AA\vec{x}=\vec{b} \Leftrightarrow \begin{cases} A\vec{y}=\vec{b}\\ A\vec{x}=\vec{y} \end{cases} \Leftrightarrow \begin{cases} A\vec{y}=\vec{b}\\ L\vec{z}=\vec{y}\\ U\vec{x}=\vec{z} \end{cases}. \end{align*}\] D’après la question précédente, on \(\vec{y}=(-8,-18,14)^T\), la résolution des deux systèmes restants conduit à \(\vec{z}=(-8,-2,-12)^T\) et \(\vec{x}=(1,10,-6)^T\).
Soit \(A\) la matrice \[ A = \begin{pmatrix} 2 & 6 & -1\\ 1 & 3 & 1\\ 1& 1 & 1 \end{pmatrix} \] 1. \(\det(A)=6\neq 0\), \(A\) est inversible.
Une matrice \(A\) admet une factorisation \(LU\) ssi : toutes les sous matrices principales \(A_k\) de la matrice \(A\) sont inversibles.
Les sous-matrices principales de \(A\) sont \(A_3=A\), \(A_2=\begin{pmatrix}2 & 6 \\1 & 3 \end{pmatrix}\) et \(A_1=2\). \(A_2\) n’est pas inversible car \(\det(A_2)=0\). Par conséquent, même si \(A\) est inversible, elle n’admet pas de décomposition \(LU\).
trisup <- function(A, b){
# trisup : resolution system triangulaire superieur
#
# Entrees : - A une matrice carre de taille nxn
# - b un vecteur de R^n
# Sortie : x solution du sytem Ax=b
#
if (ncol(A)!=nrow(A)){stop("A pas carre")}
n <- ncol(A)
x <- rep(NA, n)
x[n] <- b[n]/A[n,n]
for (i in (n-1):1){
j <- (i+1):n
x[i] <- (b[i]-sum(A[i,j]*x[j]))/A[i,i]
}
return(x)
}triinf <- function(A, b){
n <- ncol(A)
x <- rep(NA, n)
x[1] <- b[1]/A[1,1]
for (i in 2:n){
j <- 1:(i-1)
x[i] <- (b[i]-sum(A[i,j]*x[j]))/A[i,i]
}
x
}Ti <- matrix(c(1,0,0,
2,3,0,
1,4,-1), ncol = 3, byrow = T)
bi <- c(1, 8, 10)
xi <- triinf(Ti,bi)
xi## [1] 1 2 -1solve(Ti, bi)## [1] 1 2 -1Ts <- matrix(c(1,2,3,
0,4,8,
0,0,5), ncol = 3, byrow = T)
bs <- c(6, 16, 15)
xs <- trisup(Ts,bs)
xs## [1] 1 -2 3solve(Ts, bs)## [1] 1 -2 3Pour déterminier le nombre d’opérations nécessaires dans la résolution d’un système triangulaire supérieur ou inférieur, on rappelle l’algorithme :
Triangulaire supérieur (Remonter) : \[ x_n = \frac{b_n}{a_{nn}}, \] pour \(i=n-1\) à \(1\), \[ x_i=(b_i-\sum_{j=i+1}^n a_{ij}x_j)/a_{ii}. \]
L’étape \(n\) requiert 1 division. Les étapes suivantes nécessitent \((n - i)\) multiplications et \((n - i - 1)\) additions pour chaque terme de la somme, puis une soustraction et une division.
On note \(N_{MD}\) le nombre de multiplications/divisions intervenant dans l’algorithme et \(N_{SD}\) le nombre de soustractions/additions. On a donc
\[ N_{MD}=1+\sum_{i=1}^{n-1}((n-i)+1)=1+\left(\sum_{i=1}^{n-1}(n-i)\right)+n-1=n+\sum_{i=1}^{n-1}i=\frac{n^2+n}{2}, \] et \[ N_{SD}=\sum_{i=1}^{n-1}((n-i-1)+1)=\sum_{i=1}^{n-1}i=\frac{n^2-n}{2} \] Par conséquent, le nombre total d’opérations nécessaires à la résolution d’un système triangulaire est \(n^2\). Par le même raisonnement ont obtient également \(n^2\) opérations arithmétiques dans le cas d’un système triangulaire inférieur.
if(abs(prod(diag(A))) < .Machine$double.eps){
stop("akk doit etre non 0")
}sysGauss <- function(A, b){
# Mise sous forme echelonne d'un systeme Ax=b
# Entree :
# - A une matrice carre de taille n
# - b un vecteur de Rn (le second membre)
# Sortie :
# - x la solution du systeme Ax=b
if (ncol(A)!=nrow(A)){stop("A pas carre")}
n <- ncol(A)
for (k in 1:(n-1)){
for (i in (k+1):n){
mik <- A[i,k]/A[k,k]
for (j in k:n){
A[i,j] <- A[i,j] - mik*A[k,j]
}
b[i] <- b[i] - mik*b[k]
}
}
return(list("A"=A, "b"=b))
}uplowdiag <- function(x, k){
n <- length(x)
m <- n + abs(k)
y <- matrix(0, nrow=m, ncol=m)
y[col(y) == row(y) + k] <- x
y
}
n <- 3
sA <- rep(-1, n-1)
A <- 2*diag(n) + uplowdiag(sA, 1) + uplowdiag(sA, -1)
# autre option
#A <- 2*diag(n) +
# rbind(cbind(rep(0,n-1), diag(sA)), rep(0,n)) +
# rbind(rep(0,n), cbind(diag(sA), rep(0,n-1)))
b <- c(1,rep(0,n-2),1)
sys <- sysGauss(A,b)
trisup(sys$A,sys$b)## [1] 1 1 1triGauss <- sysGauss(A, b)
triGauss## $A
## [,1] [,2] [,3]
## [1,] 2 -1.0 0.000000
## [2,] 0 1.5 -1.000000
## [3,] 0 0.0 1.333333
##
## $b
## [1] 1.000000 0.500000 1.333333x <- trisup(triGauss$A, triGauss$b)
x## [1] 1 1 1solve(A,b)## [1] 1 1 1LUdec <- function(A, b){
if (ncol(A)!=nrow(A)){stop("A pas carre")}
n <- ncol(A)
L <- diag(n)
for (k in 1:(n-1)){
for (i in (k+1):n){
L[i, k] <- A[i,k]/A[k,k]
A[i, ] <- A[i, ] - L[i, k]*A[k, ]
}
}
return(list("U"=A, "L"=L))
}A5 <- matrix(c(2,4,-6,-1,-1,2,0,2,0),3,3)
A5## [,1] [,2] [,3]
## [1,] 2 -1 0
## [2,] 4 -1 2
## [3,] -6 2 0LUA5 <- LUdec(A5)
LUA5$L## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 2 1 0
## [3,] -3 -1 1LUA5$U## [,1] [,2] [,3]
## [1,] 2 -1 0
## [2,] 0 1 2
## [3,] 0 0 2LUA5$L%*%LUA5$U## [,1] [,2] [,3]
## [1,] 2 -1 0
## [2,] 4 -1 2
## [3,] -6 2 0La factorisation de Cholesky constitue une extension de la méthode \(LU\) pour les matrices symétriques définies positives. On cherche une décomposition de la forme \(A=BB^T\) où \(B\) est une matrice triangulaire inférieure dont les composantes diagonales sont positives (permettant ainsi de garantir l’unicité de la décomposition).
L’algorithme de décomposition repose sur une identification successive des éléments colonne par colonne, où, pour chaque colonne \(j\), on détermine l’élément de la ligne \(i\) pour \(i=j\) à \(n\) à l’aide des éléments précédement calculés. Un algorithme procédant sur ligne par ligne est également possible. Le choix entre les deux dépend essentiellement de la manière dont la matrice est stockée en machine (i.e, par ligne ou par colonne).
(l’algorithme est dans le cours).
facto_cholesky <- function(A){
if (all(A==t(A))!=TRUE){stop("A pas sym")}
if (all(eigen(A,only.values = T)$values >0)!=TRUE){stop("A pas > 0")}
n <- nrow(A)
L <- matrix(0,n,n)
L[1,1] <- sqrt(A[1,1])
for (j in 2:n){
L[j,1] <- A[j,1]/L[1,1]
}
for (i in 2:(n-1)){
k <- 1:(i-1)
L[i,i] <- sqrt(A[i,i]-sum(L[i,k]^2))
for (j in (i+1):n){
L[j,i] <- ( A[j,i] - sum(A[j,k]*L[i,k]) )/L[i,i]
}
}
k <- 1:(n-1)
L[n,n] <- sqrt(A[n,n]-sum(L[n,k]^2))
return(L)
}
L <- facto_cholesky(A)
L ## [,1] [,2] [,3]
## [1,] 1.4142136 0.0000000 0.000000
## [2,] -0.7071068 1.2247449 0.000000
## [3,] 0.0000000 -0.8164966 1.154701L%*%t(L)## [,1] [,2] [,3]
## [1,] 2 -1 0
## [2,] -1 2 -1
## [3,] 0 -1 2Hypothèses Cholesky :
\(A\) symétrique ;
\(A\) définie positive.
But : factoriser \(A\) sous la forme \[ A=LDL^T \] où \[ L = \begin{pmatrix} 1 & 0& \ldots & 0\\ L_{21} & \ddots & \ddots & \vdots\\ \vdots & \ddots & \ddots &0\\ L_{n1} & \ldots & L_{n-11}& 1 \end{pmatrix}, \quad \text{et}\quad D = \begin{pmatrix} D_{11}& 0 & \ldots & 0\\ 0 & \ddots & \ddots& \vdots\\ \vdots & \ddots & \ddots& \vdots\\ 0 & \ldots &0 & D_{nn} \end{pmatrix} \] Gain en terme de stockage (i.e., uniquement \(L\) et pas \(L\) et \(U\) par rapport à \(LU\)).
Inconvénient Cholesky : présence de la racine carrée (i.e, extraire une racine prend du temps).
Éviter le calcul de la racine en passant par \(LDL^T\).
Cholesky : \(A=\tilde{L}\tilde{L}^T\) avec \[ L = \begin{pmatrix} \tilde{L}_{11} & 0& \ldots & 0\\ \tilde{L}_{21} & \ddots & \ddots & \vdots\\ \vdots & \ddots & \ddots &0\\ \tilde{L}_{n1} & \ldots & \tilde{L}_{n-11}& \tilde{L}_{nn}. \end{pmatrix} \] On a \[\begin{align*} \det(A) &= \det(\tilde{L}\tilde{L}^T) = \det(\tilde{L})\det(\tilde{L}^T)\\ & = \det(\tilde{L})^2= \left(\prod_{i=1}^n \tilde{L}_{ii}\right)^2=\prod_{i=1}^n \tilde{L}_{ii}^2 \\ &\neq 0 \quad (A\in S_n^{++}(\mathbb{R})\Longrightarrow A~\text{inversible}~\Leftrightarrow \det(A)\neq 0). \end{align*}\] Soit \(\tilde{L}_k\) une sous-matrice d’ordre \(k\) de \(\tilde{L}\). Pour \(i=1,\ldots,n\) \(\tilde{L}_{ii}\neq 0 \Longrightarrow det(\tilde{L}_k)=\prod_{i=1}^n \tilde{L}_{kk}\neq 0\). Donc toutes les sous-matrices de \(\tilde{L}\) sont inversibles (\(\Longrightarrow\) LU ok).
On applique la décomposition \(LU\) à \(\tilde{L}\). On pose \(\tilde{L}=LU\).
Par suite, \(A=\tilde{L}\tilde{L}^T=LU(LU)^T=LUU^TL^T=LU^2L^T\).
En posant \(D=U^2\), on en déduit que \(A=LDL^T\).
\(L_{ij}\) pour \(i=2\) à \(n\) et \(j=1\) à \(n-1\) ;
\(D_{ii}\) pour \(i=1\) à \(n\).
Démarche : \(A=LDL^T \Rightarrow A_{ij}=(LDL^T)_{ij}=\sum_{k=1}^n (LD)_{ij} \cdot L_{kj}^T\).
Donc \[\begin{align*} A_{ij} &= \sum_{k=1}^n\left[\sum_{m=1}^n L_{im}D_{mk}\right]L_{kj}^T\\ & = \sum_{k=1}^n\left[\sum_{m=1}^n L_{im}D_{mk}\right]L_{kj}^T\\ & = \sum_{k=1}^n L_{ik}D_{kk}L_{jk} \quad (D_{mk}=0, D~\text{diag})\\ & = \sum_{k=1}^i L_{ik}D_{kk}L_{jk} \quad (L~\text{tri inf}) \end{align*}\] On distingue deux cas : 1er cas \(i=j \rightarrow D_{ii}\), 2e cas \(i> j \rightarrow L\).
1er cas \(i=j\) \[\begin{align*} A_{ii} &= \sum_{k=1}^{i-1} (L_{ik})^2D_{kk} + (L_ii)D_{ii}\\ & = \sum_{k=1}^{i-1} (L_{ik})^2D_{kk} + D_{ii}. \end{align*}\] D’où \[ D_{ii} = A_{ii} - \sum_{k=1}^{i-1} (L_{ik})^2D_{kk}. \]
2e cas \(i> j\) \[\begin{align*} A_{ij} &= \sum_{k=1}^{i-1} L_{ik}D_{kk}L_{jk} + L_{ii}D_{ii}L_{ji}\\ &= \sum_{k=1}^{i-1} L_{ik}D_{kk}L_{jk} + D_{ii}L_{ji}. \end{align*}\] D’où \[ L_{ji} = \left[A_{ij} - \sum_{k=1}^{i-1} L_{ik}D_{kk}L_{jk}\right]/D_{ii}. \]
On obtient l’algorithme suivant :
Pour \(i=1\) à \(n\) faire
\(D_{ii}=A_{ii} - \sum_{k=1}^{i-1} (L_{ik})^2D_{kk}\)
\(L_{ii}=1\)
Pour \(j=i+1\) à \(n\)
\(L_{ji} = [A_{ij} - \sum_{k=1}^{i-1} L_{ik}D_{kk}L_{jk}]/D_{ii}\)
Comme \(A=LDL^T\), la résolution de \(A\vec{x}=\vec{b}\) revient à résoudre deux systèmes triangulaires et un système diagonale. \[\begin{align*} A\vec{x}=\vec{b} \Leftrightarrow LDL^T\vec{x}=\vec{b} \Leftrightarrow \begin{cases} L\vec{y}=\vec{b}\\ D\vec{z}=\vec{y}\\ L^T\vec{x}=\vec{z} \end{cases} \end{align*}\] …
HR
\[\begin{align*} f_{n+1}f_{n+3} - f^2_{n+3} &= f_{n+1}(f_{n+2}+f_{n+1})-f_{n+2}(f_{n+1}+f_n)\\ &=f_{n+1}^2-f_{n+2}f_n+f_{n+1}f_{n+2}-f_{n+2}f_{n+1}\\ &=-(-1)^{n+1}\\ &=(-1)^{n+2} \end{align*}\] Donc \(\forall n\in\mathbb{N}\), \(f_nf_{n+2}-f^2_{n+1}=(-1)^{n+1}\).
\[ x_1 = \frac{ \begin{vmatrix} f_{n+2} & f_{n+1}\\ f_{n+3} & f_{n+2} \end{vmatrix}}{ \begin{vmatrix} f_n & f_{n+1}\\ f_{n+1} & f_{n+2} \end{vmatrix} } =\frac{(f_{n+1})^2-f_{n+1}f_{n+3}}{f_nf_{n+2}-f^2_{n+1}}=\frac{-(-1)^{n+2}}{(-1)^{n+1}}=1 \] et \[ x_2 = \frac{ \begin{vmatrix} f_{n} & f_{n+1}\\ f_{n+1} & f_{n+2} \end{vmatrix}}{ \begin{vmatrix} f_n & f_{n+1}\\ f_{n+1} & f_{n+2} \end{vmatrix} } =\frac{f_nf_{n+3}-f_{n+1}f_{n+2}}{f_nf_{n+2}-f^2_{n+1}}=\frac{f_n(f_{n+2}+f_{n+1})-f_{n+1}(f_n+f_{n+1})}{f_nf_{n+2}-f^2_{n+1}}=\frac{f_nf_{n+2}-f_{n+1}^2}{f_nf_{n+2}-f_{n+1}^2}=1. \] Donc \(\forall n\in \mathbb{N}\), \(x_1=x_2=1\).
Comme la suite \((f_n)_n\) est strictement croissante, on a \[\begin{align*} \|A_n\|_1&=\max \{|f_n|+|f_{n+1}| ; |f_{n+1}| + |f_{n+2}|\}\\ &=\max \{f_n+f_{n+1} ; f_{n+1} + f_{n+2}\}\\ &=\max \{f_{n+2} ; f_{n+3}\}\\ &= f_{n+3}. \end{align*}\]
De plus, on a \[ A^{-1}_n = \frac{1}{\det A_n} \begin{pmatrix} f_{n+2} & -f_{n+1}\\ -f_{n+1} & f_{n} \end{pmatrix} =(-1)^{n+1} \begin{pmatrix} f_{n+2} & -f_{n+1}\\ -f_{n+1} & f_{n} \end{pmatrix} \]
\[\begin{align*} \|A_n^{-1}\|_1&=\max \{f_{n+2}+f_{n+1} ; f_{n+1}+f_n\}=\max\{f_{n+3}; f_{n+2}\}=f_{n+3}\\ &= f_{n+3}. \end{align*}\]
Donc \(\mathrm{cond_1}(A_n)=(f_{n+3})^2\).
Comme \(f_n \rightarrow \infty\) quand \(n\rightarrow \infty\), on en déduit que \(\mathrm{cond_1}(A_n)\rightarrow \infty\) quand \(n\rightarrow \infty\).
\(A_n\) est un exemple de matrice mal conditionnée pour \(n\) grand.
library(expm) # pour calculer la puissance d'une matrice
An <- function(n){
f0 <- 0
f1 <- 1
f2 <- 1
A0 <- matrix (c(f0 , f1, f1, f2) ,2, 2)
Anp1 <- A0%^%n
#bnp1 <- c(Anp1[2,2], Anp1[2,2]+Anp1[2,1])
#solve(Anp1,bnp1)
#invAnp1 <- solve(Anp1)
#condAn <- norm(Anp1,"1")*norm(invAnp1,"1")
condAn <- (Anp1[2,2]+Anp1[2,1])^2
cat("cond_1(A",n-1,")=",condAn,sep = "")
}
An(1)## cond_1(A0)=4An(10)## cond_1(A9)=20736#An(50)