import numpy as np
import numpy.linalg as lin

np.set_printoptions(precision=3, linewidth=150, suppress=True)

Conditionnement d'une matrice

Soit la matrice A suivante :

A = np.array([[10, 7, 8, 7], [7, 5, 6, 5], [8, 6, 10, 9], [7, 5, 9, 10]])
A

array([[10,  7,  8,  7],
       [ 7,  5,  6,  5],
       [ 8,  6, 10,  9],
       [ 7,  5,  9, 10]])

Une matrice symétrique qui n'a rien de méchant a priori. Son déterminant est 1.

lin.det(A)

0.9999999999999869

Construisons b de telle sorte que la solution du système matriciel A x = b soit [1,1,1,1] :

b = A.sum(axis=1)
print(b)
x = lin.solve(A, b)
x

[32 23 33 31]

array([1., 1., 1., 1.])

Perturbons légèrement b. Dans le cas d'une expérience cela s'appelle une erreur de mesure. En informatique cela peut être le résultat d'erreurs d'arrondi.

bp = [32.1, 22.9, 33.1, 30.9]
eb = lin.norm(b - bp) / lin.norm(b) # une erreur se mesure par rapport à la valeur de la donnée
eb

0.0033319453118976702

On a une erreur sur b de l'ordre de 0,3 %. On la note $ ||{\bf \delta b}|| \, / \,||{\bf b}||$.

On pourrait espérer une erreur sur le résultat du même ordre de grandeur. Regardons la solution x de notre système matriciel perturbé :

xp = lin.solve(A, bp)
xp

array([  9.2, -12.6,   4.5,  -1.1])

Cette solution n'a rien à voir avec [1, 1, 1, 1].

ex = lin.norm(x - xp) / lin.norm(x)
ex

8.19847546803699

L'erreur est de l'ordre de 8 pour 1.

L'erreur est 2460 fois plus grande que l'erreur sur b.

ex / eb

2460.567236431514

Pourquoi ?

On a

$$ \begin{align} & A ({\bf x} + {\bf \delta x}) = {\bf b} + {\bf \delta b} \quad \textrm{et donc} \\ & A \, {\bf \delta x} = {\bf \delta b} \; \textrm{ puisque } A {\bf x} = {\bf b} \quad \textrm{et finalement}\\ & {\bf \delta x} = A^{-1} \, {\bf \delta b} \end{align} $$

Comme A et son inverse sont des applications linéaires on a

$$ \begin{align} & ||{\bf b}|| \le ||A|| \, ||{\bf x}|| \quad \textrm{et} \quad ||{\bf \delta x}|| \le ||A^{-1}|| \, ||{\bf \delta b}|| \end{align} $$

donc

$$ \begin{align} \frac{||{\bf \delta x}||}{||{\bf x}||} \le ||A^{-1}|| \, \frac{||{\bf \delta b}||}{||{\bf x}||} \le ||A^{-1}|| \, ||A|| \, \frac{||{\bf \delta b}||}{||{\bf b}||} \end{align} $$

lin.norm(lin.inv(A)) * lin.norm(A)

3009.5787080586942

Le problème est là.

On appelle cela le conditionnement de A :

cond(A) = $||A^{-1}|| \, ||A||$

Une matrice mal conditionnée va générer des erreurs de calcul lors de la résolution du système matriciel.

np.linalg.cond(A) # scipy n'a pas le conditionnement mais numpy l'a. 

2984.0927016757555

Bizarre cette différence avec le calul ci-dessus qui a donné 3009.

Perturbons la matrice

Que ce passe-t-il si on perturbe A et non b ?

np.random.seed(0)

dA = 2 * np.random.random(size = A.shape) - 1
display(dA)
ea = lin.norm(dA) / lin.norm(A)
print('Erreur relative sur A :', ea)

array([[ 0.098,  0.43 ,  0.206,  0.09 ],
       [-0.153,  0.292, -0.125,  0.784],
       [ 0.927, -0.233,  0.583,  0.058],
       [ 0.136,  0.851, -0.858, -0.826]])

Erreur relative sur A : 0.06868857112100454

Ap = A + dA
Ap

array([[10.098,  7.43 ,  8.206,  7.09 ],
       [ 6.847,  5.292,  5.875,  5.784],
       [ 8.927,  5.767, 10.583,  9.058],
       [ 7.136,  5.851,  8.142,  9.174]])

xp = lin.solve(Ap, b)
xp

array([-12.365,  15.574,  10.146,  -5.94 ])

ex = lin.norm(xp - x) / lin.norm(x)
ex

11.432687335993894

ex / ea

166.44235204505293

On note que l'erreur est nettement moins grande. La raison est qu'on n'a pas trouvé l'erreur sur A qui perturbera le plus possible le résultat. En fait ce n'est pas que le conditionnement de A qui compte, l'erreur est aussi importante. Deux erreurs de même norme pertuberont différemment le résultat.

Notons aussi que dans ce cas les maths sont un peu différente mais on retrouve le conditionnement de A :

$$ \begin{align} & (A + \Delta A) \, ({\bf x} + {\bf \delta x}) = {\bf b} \quad \textrm{et donc} \\ & A \, {\bf \delta x} + \Delta A \, ({\bf x} + {\bf \delta x}) = 0 \; \textrm{ puisque } A {\bf x} = {\bf b} \quad \textrm{et finalement}\\ & {\bf \delta x} = -A^{-1} \,\Delta A \, ({\bf x} + {\bf \delta x}) \quad \textrm{et} \\ & ||{\bf \delta x}|| \le ||A^{-1}|| \, ||\Delta A|| \, ||{\bf x} + {\bf \delta x}|| \end{align} $$

Ainsi

$$ \begin{align} \frac{||{\bf \delta x}||}{||{\bf x} + {\bf \delta x}||} \le ||A^{-1}|| \, ||\Delta A|| = ||A^{-1}|| \, ||A|| \, \frac{||\Delta A||}{||A||} \end{align} $$

ou

$$ \begin{align} \frac{||{\bf \delta x}||}{||{\bf x} + {\bf \delta x}||} \le cond(A) \, \frac{||\Delta A||}{||A||} \end{align} $$

L'erreur ne se mesure pas par rapport à x mais par rapport à ${\bf x} + {\bf \delta x}$.

Propriétés

• $cond(A) \ge 1$ car $Id = A\, A^{-1}$ et donc $1 \le ||A||\, ||A^{-1}||$ (pour les normes subordonnées car $||Id||_F = \sqrt{n}$)
• $cond(A) = cond(A^{-1})$ par définition du conditionnement
• $\displaystyle cond_2(A) = \frac{\max_i |\lambda_i|}{\min_i |\lambda_i|}$ si la matrice est normale où le 2 indique qu'on utilise la norme 2 et $\lambda_i$ sont les valeurs propres de A

vp = lin.eigvals(A)
vp.max() / vp.min()  # et voici la différence avec le calcul ||A|| ||inv(A)|| ci-dessus
                     # ca sent l'erreur numérique

(2984.092701676269+0j)

• si A est unitaire ou orthogonale alors $cond_2(A) = 1$ (une rotation ou symétrie ne va pas agrandir l'erreur)
• le conditionnement n'est pas modifié par transformation unitaire

Préconditionnement

Si le conditionnement n'est pas modifié par une transformation unitaire, il l'est par d'autres transformations. Ainsi

$$ \forall A, \exists B \; \textrm{appelée matrice de préconditionnement t.q.} \quad cond(B\, A) \le cond(A) $$

aussi on lieu de résoudre $A {\bf x} = {\bf b}$ on résoud $B\, A {\bf x} = B\, {\bf b}$

Toute la difficulté consiste à trouver une matrice de préconditionnement B qui soit simple à calculer.

Exercice 11.1

Quelle est la meilleure matrice de préconditionnement de A ? Pourquoi ?