Méthode du gradient conjugué

La méthode du gradient conjugué est une méthode numérique pour obtenir la solution d’un système d’équations linéaires $A.x = b$ où $A$ est une matrice symétrique définie positive. Il s’agit généralement d’une méthode itérative, qui converge vers un vecteur solution $x$. La convergence de cette méthode utilise les propriétés du produit scalaire associé à la matrice $A$. La page Wikipedia propose l’implémentation suivante en Matlab :

function x = conjgrad(A, b, x)
    r = b - A * x;
    p = r;
    rsold = r' * r;

    for i = 1:length(b)
        Ap = A * p;
        alpha = rsold / (p' * Ap);
        x = x + alpha * p;
        r = r - alpha * Ap;
        rsnew = r' * r;
        if sqrt(rsnew) < 1e-10
            break
        end
        p = r + (rsnew / rsold) * p;
        rsold = rsnew;
    end
end

Le paramètre $x$ sert d’initialisation à l’algorithme, mais peut être choisi comme égal au vecteur nul.

  • Questions :
    1. Par quels aspects cette implémentation ne respecte t’elle pas des standards de codage très sains ?
    2. Par quel aspect cette méthode est-elle différente de la méthode décrite mathématiquement dans la section précédente ? Quel est le gain de complexité occasionné ? Est-il systématique ?
    3. Implémenter la méthode du gradient conjugué et la tester.
    4. Il est possible d’utiliser cette méthode avec un préconditionneur, comme évoqué dans le paragraphe précédent. La même page propose un algorithme avec préconditionneur (mais sans pseudo-code).
    5. Implémenter la méthode du gradient conjugué avec préconditionneur et la tester.

On n’applique pas le préconditionneur $M$ sur la matrice $A$ mais sur le vecteur $A.p$ !

$M$ étant factorisée, on ne calcule pas $M^{-1}.b$ mais on résoud le système linéaire $M.x=b$, c.f. remarques précédentes.