CAMA : ma42 Gradient conjugué

Cours du 25/05

Programmer le gradient conjugué

A partir de ce cours sur le gradient conjugué programmez en Python + Numpy le gradient conjugué en exploitant les astuces mathématiques indiquées pour optimiser
votre code.

• Effectuez des tests pour valider votre code.
• Comparez la vitesse de convergence à celle du gradient avec μ optimal. Tracez des courbes de convergence (cf la feuille qui en parle)
• Comparez les temps de calcul.

Note : Veuillez écrire des fonctions les plus propres possibles, en particulier qui n'utilisent pas des variables globales comme c'est le cas dans ma correction du gradient (ma33).

Solution

import numpy as np
import scipy.linalg as lin
import matplotlib.pylab as plt

%matplotlib inline
%config InlineBackend.figure_format = 'retina'

def make_system(N):
    A = 1.0 * np.random.randint(-10, 10, size=(N,N))
    A[np.diag_indices(N)] = 0.1 + np.abs(A).sum(axis=0)    # diag dominante
    A = A + A.T                                            # symétrique
    A = A / np.abs(A).sum(axis=0).mean()
    b = 1.0 * np.random.randint(-10,10,size=(N))
    return A, b

A, b = make_system(2)
print(A, "\n\n", b)

[[ 0.65174129 -0.32338308]
 [-0.32338308  0.70149254]] 

 [ 6. -1.]

def gradient_conjugué(A, x0, b, error=1E-9, convergence=False):
    x = x0.copy()  # je ne veux pas modifier les paramètres qu'on me donne
    e2 = error**2
    r = A @ x - b  # le gradient mais aussi le résidu
    r2 = r @ r
    p = -r
    if convergence:
        conv = [np.sqrt(r2)]
    while r2 > e2:
        alpha = (r @ r) / np.dot(A @ p, p)
        x += alpha * p
        r += alpha * (A @ p)
        beta = (r @ r) / r2
        p = -r + beta * p
        r2 = r @ r
        if convergence:
            conv.append(np.sqrt(r2))
    return np.array(conv) if convergence else x

gradient_conjugué(A, np.array([0.,0.]), b, convergence=True)

array([6.08276253e+00, 2.51964707e+00, 2.22044605e-16])

def compute_error(N, method=gradient_conjugué):
    A, b = make_system(N)
    x = method(A, np.zeros(N), b)
    err = A @ x - b
    return np.sqrt(err @ err)

compute_error(10)

4.165926057296536e-15

Comparons avec le gradient simple

def gradient(A, x0, b, e = 1E-9, convergence=False, max_iterations=1000):
    x = x0.copy()
    e2 = e**2
    k = 0
    gradJ = A @ x - b
    g2 = gradJ @ gradJ
    divergence_limite = 1E6 * g2
    if convergence:
        conv = [np.sqrt(g2)]
    while g2 > e2:
        µ = np.dot(gradJ, gradJ) / np.dot(A @ gradJ, gradJ)
        x -= µ * gradJ
        gradJ = A @ x - b
        g2 = gradJ @ gradJ
        if convergence:
            conv.append(np.sqrt(g2))

        # la suite n'est que des tests pour se protéger
        if g2 > divergence_limite:  # au cas où on diverge
            print("DIVERGE")
            break
        k += 1
        if k > max_iterations:  # c'est trop long, je crains la boucle infinie
            print('Trop long, boucle infinie ?')
            break
    return np.array(conv) if convergence else x

# vérifions que ca marche

compute_error(10, method=gradient)

6.767792116739432e-10

Perfs

# comparons les performances

seed = 123
np.random.seed(seed)

%timeit compute_error(1000, method=gradient)

34 ms ± 4.23 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

seed = 123
np.random.seed(seed)

%timeit compute_error(1000, method=gradient_conjugué)

32.1 ms ± 847 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Le gain n'est pas clair…

Nombre d'iteration dans les 2 cas

N = 1000
A,b = make_system(N)
x0 = np.zeros(N)

Pour le gradient simple

err = gradient(A, x0, b, convergence=True)

plt.plot(np.arange(err.shape[0]), err)
plt.title('Convergence du gradient')
plt.semilogy();

Pour le gradient conjuge

err = gradient_conjugué(A, x0, b, convergence=True)

plt.plot(np.arange(err.shape[0]), err)
plt.title('Convergence du gradient conjugué')
plt.semilogy();

Argh, le gradient conjugué n'est pas la révolution prédite !

Un cas réel

Logiquement vous devriez être décu aussi on va tester avec un problème réel qui correspond à cet exemple de l'équation de Poisson. Le système matriciel de ce problème est téléchargeable ici. Une fois le fichier sauvé, pour récupérer A et b faites :

npz = np.load('/tmp/Ab.npz')
A = npz['A']
b = npz['b']

• Faites une étude rapide de A, indiquez quel pourcentage des valeurs de A est différent de 0. Afficher l'image de la matrice avec plt.imshow(A) (faire une grande image pour voir quelque chose).
• Refaites la comparaison entre les deux méthodes avec ce système matriciel.
• Regardez la documentation de lin.solve (en particulier les options) et comparer lin.solve à vos deux algorithmes.

Solution

print(A.min(), A.max())

-1.5693731138089555 4.357203686821435

diff0 = (A != 0).sum() / (A.shape[0] * A.shape[1])
print(f"Pourcentage de valeurs != 0 : {100 * diff0:.3} %")

Pourcentage de valeurs != 0 : 0.339 %

plt.figure(figsize=(15,15))
plt.imshow(A)

Comparaison gradient simple et conjugué

%time err = gradient_conjugué(A, np.zeros(len(A)), b, convergence=True)

CPU times: user 1.75 s, sys: 155 ms, total: 1.91 s
Wall time: 521 ms

plt.plot(np.arange(err.shape[0]), err)
plt.title('Convergence du gradient conjugué')
plt.semilogy();

# le gradient simple

%time err = gradient(A, np.zeros(len(A)), b, convergence=True, max_iterations=10000)

CPU times: user 1min 11s, sys: 4.3 s, total: 1min 15s
Wall time: 19.4 s

plt.plot(np.arange(err.shape[0]), err)
plt.title('Convergence du gradient')
plt.semilogy();

On voit la supériorité du gradient conjugué tant en nombre d'itérations (175 contre 7800) qu'en temps de calcul (0,5 s contre 20 s).

Comparaison avec lin.solve

?lin.solve

%time x = lin.solve(A, b, assume_a='pos')

CPU times: user 170 ms, sys: 85.1 ms, total: 255 ms
Wall time: 94.4 ms

r = A @ x - b
r @ r

2.0359257929180678e-25

On note aussi lin.solve est plus rapide et sa solution est nettement meilleure… lin.solve utilise une méthode directe ici. Cela est dû au fait que Scipy utilise la bibliothèque Lapack (qui est imbatable).

Le gradient conjugué de Scipy (avec Lapack)

Le gradient conjugué à tout son sens pour les matrices creuses aussi il est dans la partie "sparse" de Scipy. On a vu que notre matrice à plus de 99 % de valeur nulles ce qui en fait bien une matrice creuse. Aussi je la charge dans le format COO qui ne stocke que les valeurs non nulles et:

import scipy.sparse as sparse
from scipy.sparse.linalg import cg

Ac = sparse.load_npz('/tmp/Acoo.npz')
%time x = cg(Ac, b)

CPU times: user 24.8 ms, sys: 15.1 ms, total: 39.9 ms
Wall time: 10.8 ms

On gagne un facteur 10 !