La méthode de Cholesky est une méthode de factorisation de matrices symétriques.
Les matrices $D+H$ et $D+V$ étant symétriques (car $D/2+H$ et $D/2+V$ le sont), on peut les factoriser par la méthode de Cholesky. Il est judicieux d'employer cette méthode pour réduire la compléxité dans la résolution des systèmes à résoudre à chaque itération. Voici une estimation du gain par rapport à résoudre ces systèmes directement par un algorithme de Gauss.
Une fois les décompositions calculées, chaque inversion coûtera de l'ordre de $n^2$ opérations (une descente et une remontée). Donc, par itération on fait de l'ordre de $2n^2$ opérations car on doit résoudre deux systèmes. Si l'on converge en $K$ itérations, cela fait une compléxité de $2Kn^2$ au cours de l'algorithme. Le calcul d'une décomposition de Cholesky coûtant de l'ordre de $n^3/6$, le coût total estimé est de $n^3/3+2Kn^2$.
Si l'on résout les systèmes linéaires directement par un algorithme de Gauss, cela coûtera de l'ordre de $2Kn^3/3$ ce qui, en fonction de la valeur de $K$, peut être beaucoup plus coûteux.
# Question 2 and 3
%pylab inline
# We build a simple symetric positive definite matrix to
# test the rest of the program
def getH(n):
H=zeros((n,n))
for i in range(n):
for j in range(n):
H[i,j]=1.+min(i,j)
return H
def Cholesky(A):
(n,m)=A.shape
if n!=m:
print('Error: not a square matrix')
return -1.
else:
B=zeros((n,n))
for j in range(n):
if j==0:
B[0,0]=sqrt(A[0,0])
B[1:n,0]=A[1:n,0]/B[0,0]
else:
B[j,j]=sqrt(A[j,j]-sum(B[j,0:j]*B[j,0:j]))
for i in range(j+1,n):
B[i,j]=(A[i,j]-sum(B[j,0:j]*B[i,0:j]))/B[j,j]
return B
def solveTriangularSystem(A,b,option):
(n,m)=A.shape
if n!=m:
print('Error: not a square matrix')
return -1.
else:
x=zeros(A.shape[0])
if option==1: # Remontée
for j in range(n-1,-1,-1):
x[j]=(b[j]-sum(A[j,j:n]*x[j:n]))/A[j,j]
else: # Descente
for j in range(n):
x[j]=(b[j]-sum(A[j,0:j]*x[0:j]))/A[j,j]
return x
def resoudreAvecCholesky(B,c):
y=solveTriangularSystem(B,c,-1)
return solveTriangularSystem(B.T,y,1)
def resolutionItertive(D,H,V,b,x0,eps,kmax):
A=D+H+V
B1=Cholesky(D+H)
B2=Cholesky(D+V)
k=0
x=x0
res=norm(dot(A,x)-b)
while(res>eps and k<kmax):
k+=1
x=resoudreAvecCholesky(B1,-dot(V,x)+b)
x=resoudreAvecCholesky(B2,-dot(H,x)+b)
res=norm(dot(A,x)-b)
return [x,res,k]
# Remarque: Le main n'était pas demandé à l'examen, mais
# il est nécessaire pour appeler les méthodes!
# Création des matrices
n=10 # Dimension de la matrice
c=2. # Création de C
D=c*eye(n)
H=getH(n) # On crée H et V qui satisfont les hypothèses
V=H
A=D+H+V
# Résolution itérative du système
b=ones(n) # Second membre
x0=zeros(n) # Vecteur d'initialisation
eps=1.E-4 # Tolérance pour le résidu
kmax=1000 # Nombre max d'itérations (non demandé à l'examen)
[x,res,k]=resolutionItertive(D,H,V,b,x0,eps,kmax)
print([res,k])