CAMA : ma20 Convergence de Jacobi avec inertie

Cours du 11 / 05

Ajouter de l'inertie à Jacobi

La méthode de Jacobi mène au système itératif :

$${\mathbf{x}}^{k+1}=M^{-1}(N{\mathbf{x}}^k+\mathbf{b})$$

Cette méthode converge ssi la matrice

$b$ a un rayon spectral inférieur à 1 (cf ma12).

On peut agrandir le rayon de convergence en ajoutant de l'inertie:

$${\mathbf{x}}^{k+1}=w{M}^{-1}(N{\mathbf{x}}^k+\mathbf{b})+(1-w){\mathbf{x}}^k$$

• $0<w\le 1$.

• si
  $w=1$ : Jacobi classique
• si
  $w=0$ : on néglige les termes en dehors de la diagonale et
  $b$ donc ça ne marche pas

On parle d'inertie car on avance "moins vite": la nouvelle valeur de

${\mathbf{x}}^{k+1}$ est comprise entre l'ancienne valeur de

${\mathbf{x}}^{k+1}$ et

${\mathbf{x}}^k$. C'est la surrelaxation

Programmons l'inertie pour Jacobi

On commence pour un Jacobi qui diverge.

np.random.seed(799)

A = np.random.randint(10, size=(4,4))
b = A.sum(axis=1) # la solution est [1,1,1,1]

A:
 [[5 7 6 0]
 [1 7 2 5]
 [5 6 5 1]
 [0 6 3 7]] 
b:
 [18 15 17 16] 

M = np.diag(A)        # vecteur
N = np.diag(M) - A    # np.diag d'une matrice donne un vecteur, np.diag d'un vecteur donne une matrice

M:
 [[5 0 0 0]
 [0 7 0 0]
 [0 0 5 0]
 [0 0 0 7]]
N:
 [[ 0 -7 -6  0]
 [-1  0 -2 -5]
 [-5 -6  0 -1]
 [ 0 -6 -3  0]]

x0 = np.random.random(4)
x = x0
for i in range(20):
    x = (N @ x + b) / M

...
x_16 = [-4448.651 -1888.411 -4149.91  -1981.882]
x_17 = [7627.267 3238.983 7114.521 3399.456]
x_18 = [-13068.402  -5548.37  -12190.539  -5823.066]
x_19 = [22399.965  9511.401 20894.459  9982.548]

Ajoutons de l'inertie :

x = x0  # on reprend la même valeur initiale pour la comparaison
w = 0.5 # on choisit w 
for i in range(20):
    x = w * (N @ x + b) / M + (1-w) * x

...
x_17 = [1.059 0.977 0.972 1.03 ]
x_18 = [1.063 0.977 0.968 1.031]
x_19 = [1.067 0.978 0.963 1.032]

La solution est [1,1,1,1], l'inertie fonctionne.

Étudions la convergence

On trace une courbe de:

• l'erreur absolue (lorsqu'on connait la solution)
• de l'erreur relative (entre 2
  ${\mathbf{x}}^i$ successifs)
• du résidu (
  $||A{\mathbf{x}}^i-\mathbf{b}||$).

x = x0    # on reprend la même valeur initiale pour la comparaison
w = 0.5   # on choisit w 
error = [np.square(x - np.ones(4)).sum()]
for i in range(20):
    x = w * (N @ x + b) / M + (1-w) * x
    error.append(np.square(x - np.ones(4)).sum())

A l'échelle logarithmique:

Il faut toujours regarder une erreur en échelle logarithmique.

En faisant le calcul sur 200 itérations :

On s'est rapproché de la solution puis on a divergé.

Erreur relative

Regardons l'écart entre 2

${\mathbf{x}}^k$ successifs.

x = x0 # on reprend la même valeur initiale pour la comparaison
w = 0.5 # on choisit w 
error2 = []
for i in range(200):
    old_x = x
    x = w * (N @ x + b) / M + (1-w) * x
    error2.append(np.square(x - old_x).sum())

Il y a une relation entre l'écart de deux valeurs successives et l'erreur absolue.

L'écart entre 2

$\mathbf{x}$ successifs est une facon de savoir quand arrêter un algorithme itératif.

Résidu

x = x0 # on reprend la même valeur initiale pour la comparaison
w = 0.5 # on choisit w 
residu = []
for i in range(200):
    old_x = x
    x = w * (N @ x + b) / M + (1-w) * x
    residu.append(np.square(A @ x - b).sum())

Normaliser

• Si la solution est un milliard, avoir une erreur de 0.1 est très bien.
• Si la solution est 0.01, avoir une erreur de 0.1 est énorme.

On ne peut juger une erreur qu'avec une référence. Si on connait la solution exacte :

$$\frac{||{\mathbf{x}}^k-\mathbf{x}||}{||\mathbf{x}||}$$

De même, l'erreur entre 2 itérations successives doit être normalisée :

$$\frac{||{\mathbf{x}}^{k+1}-{\mathbf{x}}^k||}{||{\mathbf{x}}^k||}$$

def mk_A(seed):
    np.random.seed(seed)
    return np.random.randint(10, size=(4,4))

def plot_error(M, N, b, x0, w, n=200):
    x = x0 
    error = [np.square(x - np.ones(4)).sum()]
    error2 = []
    for i in range(n):
        old_x = x
        x = w * (N @ x + b) / M + (1-w) * x
        error.append(np.square(x - np.ones(4)).sum())
        error2.append(np.square(x - old_x).sum())

def plot_error_normalized(M, N, b, x0, w, n=200):
    x = x0 
    error = [np.square(x - np.ones(4)).sum()]
    error2 = []
    for i in range(n):
        old_x = x
        x = w * (N @ x + b) / M + (1-w) * x
        error.append((np.square(x - np.ones(4)).sum())/4) # normalisé par rapport à la solution
        error2.append((np.square(x - old_x).sum())/np.square(x).sum()) # normalisé par rapport à x

A = mk_A(799)
b = A.sum(axis=1)                    

M = np.diag(A)    
N = np.diag(M) - A 

x0 = np.random.random(4)

plot_error(M, N, b, x0, w=0.1, n=1000) 

plot_error_normalized(M, N, b, x0, w=0.1, n=1000) 

L'erreur relative normalisée se stabilise.