CAMA : ma11 Conditionnement d'une matrice

Cours du 27 / 04

Soit la matrice symétrique

$A$ suivante :

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

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

lin.det(A) # calcul son determinant

0.9999999999999869  # on peut arrondir a 1

Construisons

$b$ tel que

$A\mathbf{x}=b$ et

$\mathbf{x}\mathbf{=}\mathbf{[}\mathbf{1}\mathbf{,}\mathbf{1}\mathbf{,}\mathbf{1}\mathbf{,}\mathbf{1}\mathbf{]}$

b = A.sum(axis=1)

[32 23 33 31]

x = lin.solve(A, b)

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

Perturbons

$b$, comme s'il y avait une erreur de mesure ou 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

0.0033319453118976702

On a une erreur de l'ordre de

$0,3\mathrm{\%}$.

On note l'erreur :

$$\frac{||\delta \mathbf{b}||}{||\mathbf{b}||}$$

Regardons la solution

$\mathbf{x}$ de notre système matriciel perturbé:

xp = lin.solve(A, bp)

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

La solution n'a rien n'a voir avec

$[1,1,1,1]$

ex = lin.norm(x - xp) / lin.norm(x) #mesure de l'erreur

8.19847546803699

L'erreur est de l'ordre de 8.

ex / eb

2460.567236431514

C'est

$2460$ fois plus que l'erreur sur

$b$.

Pourquoi ?

$$\begin{array}{rl}A(\mathbf{x}+\delta \mathbf{x})& =\mathbf{b}+\delta \mathbf{b}\text{et donc}\\ A\delta \mathbf{x}& =\delta \mathbf{b}\text{ puisque }A\mathbf{x}=\mathbf{b}\text{et finalement}\\ \delta \mathbf{x}& =A^{-1}\delta \mathbf{b}\end{array}$$
Comme

$A$ et son inverse sont des applications linéaires :

$$\begin{array}{rl}& ||\mathbf{b}||\le ||A||||\mathbf{x}||\text{et}||\delta \mathbf{x}||\le ||A^{-1}||||\delta \mathbf{b}||\end{array}$$
donc :

$$\begin{array}{r}\frac{||\delta \mathbf{x}||}{||\mathbf{x}||}\le ||A^{-1}||\frac{||\delta \mathbf{b}||}{||\mathbf{x}||}\le ||A^{-1}||||A||\frac{||\delta \mathbf{b}||}{||\mathbf{b}||}\end{array}$$

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

3009.5787080586942

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 # different de 3009

Perturbons la matrice

np.random.seed(0)

dA = 2 * np.random.random(size = A.shape) - 1

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]])

ea = lin.norm(dA) / lin.norm(A) # erreur relative sur A

0.06868857112100454

Ap = A + dA

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)

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

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

11.432687335993894

ex / ea # valeur de l'erreur

166.44235204505293

L'erreur est moins grande.

Une erreur peut fortement perturber

$A$, le conditionnement et l'erreur sont tous les deux importants.

Pour retrouver le conditionnement de

$A$ dans ce cas :

$$\begin{array}{rl}& (A+\mathrm{\Delta }A)(\mathbf{x}+\delta \mathbf{x})=\mathbf{b}\text{et donc}\\ & A\delta \mathbf{x}+\mathrm{\Delta }A(\mathbf{x}+\delta \mathbf{x})=0\text{ puisque }A\mathbf{x}=\mathbf{b}\text{et finalement}\\ & \delta \mathbf{x}=-A^{-1}\mathrm{\Delta }A(\mathbf{x}+\delta \mathbf{x})\text{et}\\ & ||\delta \mathbf{x}||\le ||A^{-1}||||\mathrm{\Delta }A||||\mathbf{x}+\delta \mathbf{x}||\end{array}$$
Donc

$$\begin{array}{r}\frac{||\delta \mathbf{x}||}{||\mathbf{x}+\delta \mathbf{x}||}\le ||A^{-1}||||\mathrm{\Delta }A||=||A^{-1}||||A||\frac{||\mathrm{\Delta }A||}{||A||}\end{array}$$
et

$$\begin{array}{r}\frac{||\delta \mathbf{x}||}{||\mathbf{x}+\delta \mathbf{x}||}\le cond(A)\frac{||\mathrm{\Delta }A||}{||A||}\end{array}$$

Propriétés

• $cond(A)\ge 1$ car
  $Id=A{A}^{-1}$ et donc
  $1\le ||A||||A^{-1}||$
• $cond(A)=cond(A^{-1})$
• ${\displaystyle con{d}_2(A)=\frac{\underset{i}{max}|{\lambda }_i|}{\underset{i}{min}|{\lambda }_i|}}$ si la matrice est réelle
  • 2 indique qu'on utilise la norme 2
  • ${\lambda }_i$ sont les valeurs propres de A

vp = lin.eigvals(A)
vp.max() / vp.min()

• si A est unitaire ou orthogonale :
  $con{d}_2(A)=1$
• le conditionnement n'est pas modifié par transformation unitaire

Préconditionnement

Le conditionnement peut etre tranformé :

$$\mathrm{\forall }A,\mathrm{\exists }B\text{appelée matrice de préconditionnement t.q.}cond(BA)<cond(A)$$

On résoud

$BA\mathbf{x}=B\mathbf{b}$.