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

import numpy as np
import numpy.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)

def gradient_conjugué(A, x0, b, error=1E-9, convergence=False):

# vérfions que ca marche

# combien faut-il d'itérations

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)

Comparons avec le gradient simple

def gradient(A, x0, b, e = 1E-9, convergence=False, max_iterations=1000):

# vérifions que ca marche

compute_error(10, method=gradient)

Perfs

# comparons les performances

seed = 123
np.random.seed(seed)

%timeit compute_error(1000, method=gradient)

seed = 123
np.random.seed(seed)

%timeit compute_error(1000, method=gradient_conjugué)

Nombre d'iteration dans les 2 cas

On teste avec N=1000 puis N=10000.

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

Pour le gradient simple

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

#plot

Pour le gradient conjugué

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

#plot

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',allow_pickle=True)
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.

Comparaison gradient simple et conjugué

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

#plot

# le gradient simple

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

# plot

On voit la supériorité du gradient conjugué tant en nombre d'itérations (175 contre 7800) qu'en temps de calcul (environ 40 fois plus rapide).

Comparaison avec lin.solve de Scipy

Le solveur de Scipy a des option que celui de Numpy n'a pas.

# regarder la doc pour avoir les options optimales

# calculer le résidu

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 (téléchargez Acoo.npz) 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)

On gagne un ordre de grandeur.