Lien de la note Hackmd
Cours du 04/05
Programmation vectorielle
Le but des exercices est
- d’avoir un programme qui donne la bonne réponse
- qui soit le plus rapide possible (et pour cela on utilise massivement Numpy)
En règle général si vous avez des for imbriqués c’est mauvais signe.
Méthode du pivot de Gauss partiel
L’ennoncé est dans le cours
Solution
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
def solve_gauss_partial(A, b): # on prend le max dans la colonne i parmi les lignes en dessous (plus facile)
for i in range(len(A)-1):
pivot = np.argmax(np.abs(A[i:, i])) # il n'y a que 3 lignes à ajouter pour échanger les lignes
A[[i, pivot]] = A[[pivot, i]]
b[[i, pivot]] = b[[pivot, i]]
E = np.diag(np.array([1.,] * len(A), dtype=A.dtype))
coefs = - A[i+1:,i] / A[i,i]
E[i+1:,i] = coefs
A[i:, i:] = E[i:,i:] @ A[i:,i:]
b[i+1:] += coefs * b[i] # multiplication terme à terme
# A est maintenant triangulaire surpérieur
res = np.zeros(len(b), dtype=b.dtype)
res[-1] = b[-1] / A[-1,-1]
for i in range(len(A)-1)[::-1]:
res[i] = (b[i] - A[i,i+1:] @ res[i+1:]) / A[i,i]
return res
1
2
3
4
5
6
7
e = 1E-6
A = np.array([[e, 1], [1, 2]], dtype='float32')
b = np.array([1., 3.], dtype='float32')
print(f"A\n {A} \nb\n {b}\n")
x = solve_gauss_partial(A, b)
print('solution : ',x)
print('vérification\n', A@x)
1
2
3
4
5
6
7
8
9
A
[[0.000001 1. ]
[1. 2. ]]
b
[1. 3.]
solution : [1.0000019 0.99999905]
vérification
[3. 0.999997]
Factorisation de Choleski
Il s’agit de décomposer A en $A = B\, B^T$ avec B une matrice triangulaire inférieure. Cela n’est possible que si la matrice A est symétrique et définie positive (c’est d’ailleurs un facon de vérifier qu’une matrice est définie positive).
Écrire l’algorithme de Choleski qui prend A et retourne B (pour deviner l’algorithme, essayez de trouver les coefficients de B à partir des coefficients de A sur une matrice A 4x4).
Solution
\(A = B\, B^T = \begin{bmatrix} b_{11} & 0 & \dots & 0\\ b_{21} & b_{22} & \dots & 0\\ & \vdots&\\ b_{n1} & b_{n2} & \dots& b_{n,n} \end{bmatrix} \begin{bmatrix} b_{11} & b_{21} & \dots & b_{n1}\\ 0 & b_{22} & \dots & b_{n2}\\ & \vdots&\\ b_{n1} & b_{n2} & \dots& b_{n,n} \end{bmatrix}= \begin{bmatrix} b_{11}^2 & b_{11}b_{21} & \dots & b_{11}b_{n1}\\ x & \sum_{i=1}^2b_{2i}^2 & \dots & \sum_{i=1}^2b_{2i}b_{ni}\\ & & \vdots&\\ x & x & \dots& \sum_{i=1}^2b_{n,i}^2 \end{bmatrix}\) avec $x$ la même valeur que de l’autre coté de la diagonale
On voit que $b_{11} = \sqrt{a_{11}}$ et maintenant qu’on a $b_{11}$ on peut trouver toute la première ligne de $B^T$ : $b_{j1}=a_{1j}/b_{11}$.
Une fois qu’on connait la première ligne de $B^T$ , on s’attaque à la deuxième en commencant par trouver $b_{22}$ puis ensuite tous les autres éléments de la ligne comme on a fait pour la première ligne.
On a donc dans le cas général :
- $b_{ii} = \sqrt{a_{ii} - \sum_{k=1}^{i-1}b_{ik}^2}$
- $b_{ji} = a_{ij} - \sum_{k=1}^{i-1}b_{ik}b_{jk}/b_{ii} = a_{ij} - \sum_{k=1}^{i-1}b_{ik}b_{kj}^T/b_{ii} \space\forall j\gt i$
1
2
3
4
5
6
def Choleski(A):
B = np.zeros(A.shape)
for i in range(len(A)):
B[i,i] = np.sqrt(A[i,i] - np.sum(np.square(B[i, :i]))) # garanti ok car A est def positive
B[i+1:, i] = (A[i, i+1:] - B[i, :i] @ B.T[:i, i+1:]) / B[i,i] # les ∑ sous forme de prod. scalaire
return B
Matrice symetrique
Rappel : pas de boucles for imbriquées mais des opérations vectorielles et matricielles (sur des sous-matrices).
Créer une matrice symétrique définie positive est vérifier que votre programme marche bien.
Solution
1
2
3
4
5
6
7
A = np.random.randint(10, size=(4,4))
A = A + A.T # symmétrique
A = A + np.diag(A.sum(axis=0)) # diagonale dominante
print('A:\n', A)
B = Choleski(A)
print('B\n', B)
print('vérification\n', B @ B.T)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
A:
[[55 8 18 5]
[ 8 33 7 10]
[18 7 54 9]
[ 5 10 9 28]]
B
[[7.4161984871 0. 0. 0. ]
[1.0787197799 5.6423721639 0. 0. ]
[2.4271195049 0.7765914857 6.8924593995 0. ]
[0.6741998625 1.6434093681 0.8831939788 4.9055711788]]
vérification
[[55. 8. 18. 5.]
[ 8. 33. 7. 10.]
[18. 7. 54. 9.]
[ 5. 10. 9. 28.]]
Améliorer Jacobi
Lorsqu’on écrit une itération de la méthode de Jacobi avec l’ensemble des coefficients, on constate que si on calcule la nouvelle valeur de x élément par élement alors lorsqu’on veut mettre à jour x[1], on connait déjà x[0]. Idem lorsqu’on met à jour x[2] on connait déjà x[0] et x[1], etc.
L’idée de la version amméliorée de Jacobi est d’utiliser la nouvelle valeur de x[0] et non pas l’ancienne comme c’est le cas dans l’algorithme du cours. Ainsi en utilisant des valeurs mise à jour on peut espérer converger plus vite.
Dans ce chaque itération demande un calcul ligne par ligne et donc une boucle for dans la boucle for sur les itérations.
Test d’arrêt
On ajoutera un argument error à la fonction qui indique la précision désirée du résultat. Par défaut sa valeur est de 1E-6 et pour offrir une bonne garantie on arrête l’algorithme lorsque $||x_{t+1} - x_t|| < \texttt{error}\, / \, 10$.
Solution
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
def Jacobi(A, b, error=1E-6, verbose=False):
L = np.tril(A)
U = -np.triu(A, k=1)
if verbose:
print(f"L:\n {L}\nU\n {U}\n")
previous_x = np.zeros(len(b))
x = np.random.random(len(b))
err = (error / 10) ** 2
while np.sum(np.square(x - previous_x)) > err:
previous_x = x.copy()
if verbose:
print(f"x = {x}")
# on résoud L x = U x + b avec L matrice triangulaire inférieure
y = U @ x + b
x[0] = y[0] / L[0,0]
for i in range(1,len(L)):
x[i] = (y[i] - L[i,:i] @ x[:i]) / L[i,i]
return x
1
2
3
4
5
6
A = np.random.randint(10, size=(4,4))
A = A + np.diag(A.sum(axis=0))
b = A.sum(axis=1) # ainsi la solution est [1,1,1,1]
print('A:\n', A, "\nb:\n", b, "\n")
Jacobi(A,b, verbose=True)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
A:
[[24 2 1 7]
[ 5 19 4 6]
[ 9 2 20 9]
[ 2 9 9 32]]
b:
[34 34 40 52]
L:
[[24 0 0 0]
[ 5 19 0 0]
[ 9 2 20 0]
[ 2 9 9 32]]
U
[[ 0 -2 -1 -7]
[ 0 0 -4 -6]
[ 0 0 0 -9]
[ 0 0 0 0]]
x = [0.8870874823 0.8448958895 0.2146829205 0.8281640711]
x = [1.0957657001 1.194392389 1.0147923641 0.9351814319]
x = [1.0020897014 1.0168049182 1.0265474982 0.9876765266]
x = [1.0010877908 0.9980164155 1.0052544156 0.9990120918]
x = [1.0002345046 0.9991440665 1.000424625 1.000106649 ]
x = [1.0000225291 0.9998709979 0.9999547701 1.0000475947]
x = [0.999998753 0.9999948204 0.9999796615 1.0000072549]
x = [0.9999991631 1.000002211 0.9999968908 1.0000003049]
x = [0.9999998564 1.0000005961 0.9999998678 0.9999998785]
x = [0.9999999913 1.0000000685 1.0000000518 0.9999999667]
1
array([1.0000000018, 0.9999999991, 1.0000000142, 0.9999999961])