La méthode du gradient conjugué

La méthode du gradient consiste essentiellement en deux phases: choisir une direction de descente (celle du résidu) et un point où \(\Phi\) atteint un minimum local le long de cette direction. La seconde phase est indépendante de la première. En effet, étant donné une direction \(\mathbf{p}^{(k)}\), on peut choisir \(\alpha_k\) comme étant la valeur du paramètre \(\alpha\) telle que \(\Phi(\mathbf{x}^{(k)} + \alpha \mathbf{p}^{(k)})\) est minimum. En écrivant que la dérivée par rapport à \(\alpha\) s’annule au point où la fonction admet un minimum local, on obtient

\[\label{alphap} \alpha_k = \frac{ ({\mathbf{p}^{(k)}},\mathbf{r}^{(k)} )}{ ({\mathbf{p}^{(k)}}, A \mathbf{p}^{(k)} )}\]
La question est de savoir comment choisir \(\mathbf{p}^{(k)}\).

Une alternative à l’approche qui nous a conduit à identifier \(\mathbf{p}^{(k)}\) avec \(\mathbf{r}^{(k)}\) est proposée dans la suite.

Définition: Vecteur optimal

Un vecteur \(\mathbf{x}^{(k)}\) est dit optimal par rapport à une direction \(\mathbf{p} \neq \mathbf{0}\) si

\[\label{star2} \Phi ( {\bf x}^{(k)} ) \leq \Phi ( {\bf x}^{(k)} + \lambda {\bf p} ), \qquad \forall \lambda \in \mathbb{R}\]

Si \(\mathbf{x}^{(k)}\) est optimal par rapport à n’importe quelle direction d’un espace vectoriel \(V\), on dit que \(\mathbf{x}^{(k)}\) est optimal par rapport à \(V\).

Il découle de la définition de l’optimalité que \(\mathbf{p}\) doit être orthogonal au résidu \(\mathbf{r}^{(k)}\). En effet, on déduit de [eq18] que \(\Phi\) admet un minimum local le long de \(\mathbf{p}\) pour \(\lambda=0\). La dérivée de \(\Phi\) par rapport à \(\lambda\) doit donc s’annuler pour \(\lambda=0\). Or

\[\frac{d}{d \lambda} \Phi( \mathbf{x}^{(k)} + \lambda \mathbf{p} ) = (\mathbf{p}, A \mathbf{x}^{(k)} - \mathbf{b}) + \lambda (\mathbf{p},A\mathbf{p})\]

On a donc

\[\left.\frac{d}{d \lambda} \Phi( \mathbf{x}^{(k)} + \lambda \mathbf{p}) \right|_{\lambda = 0} = 0 \quad \Longleftrightarrow \quad (\mathbf{p},\mathbf{r}^{(k)}) = 0\]

c’est-à-dire, \(\mathbf{p} \perp \mathbf{r}^{(k)}\). Remarquer que le vecteur \(\mathbf{x}^{(k+1)}\) de la méthode du gradient (non préconditionné) est optimal par rapport à \(\mathbf{r}^{(k)}\) puisque, d’après le choix de \(\alpha_k\), on a \(\mathbf{r}^{(k+1)} \perp \mathbf{r}^{(k)}\), mais cette propriété n’est plus vraie pour l’itérée suivante \(\mathbf{x}^{(k+2)}\). On est donc amené à se demander s’il existe des directions de descente qui conservent l’optimalité des vecteurs à chaque itération. Soit donc

\[\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + \mathbf{q}\]

et supposons \(\mathbf{x}^{(k)}\) optimal par rapport à une direction \(\mathbf{p}\) (ainsi \(\mathbf{r}^{(k)} \perp \mathbf{p}\)). Imposons que \(\mathbf{x}^{(k+1)}\) demeure optimal par rapport à \(\mathbf{p}\), c’est-à-dire \(\mathbf{r}^{(k+1)} \perp \mathbf{p}\). On obtient alors

\[0 = (\mathbf{p},\mathbf{r}^{(k+1)}) = (\mathbf{p},\mathbf{r}^{(k)} - A\mathbf{q} ) = - (\mathbf{p},A\mathbf{q})\]

Ainsi, pour préserver l’optimalité entre les itérées successives, les directions de descentes doivent être mutuellement A-orthogonales (on dit aussi A-conjuguées), i.e.

\[(\mathbf{p},A\mathbf{q}) = 0\]

Une méthode qui utilise des directions de descente A-conjuguées est dite conjuguée. Expliquons maintenant comment générer automatiquement une suite de directions conjuguées. Soit \(\mathbf{p}^{(0)} = \mathbf{r}^{(0)}\) et cherchons une direction de la forme

\[\label{zerocg} \mathbf{p}^{(k+1)} = \mathbf{r}^{(k+1)} - \beta_k \mathbf{p}^{(k)}, \quad k = 0, 1, \ldots\]

où \(\beta_k \in \mathbb{R}\) est à déterminer de manière à ce que

\[\label{unocg} (A \mathbf{p}^{(j)},\mathbf{p}^{(k+1)}) = 0, \quad j = 0, 1, \ldots, k\]

En écrivant [eq20] pour \(j=k\), on obtient à partir de [eq19]

\[\label{betak} \beta_k = \frac{( A\mathbf{p}^{(k)}, \mathbf{r}^{(k+1)})} {( A\mathbf{p}^{(k)},\mathbf{p}^{(k)})}, \quad k = 0, 1, \ldots\]

On peut montrer (par recurrence) que l’égalité [eq20] est aussi vérifiée pour \(j = 0, 1, \ldots, k-1\).

La méthode du gradient conjugué (notée GC) est donc la méthode obtenue en choisissant comme directions de descente les vecteurs \({\bf p}^{(k)}\) donnés par [eq19] et comme paramètres d’accélération les \(\alpha_k\) définis en [eq17]. Par conséquent, en posant \({\bf r}^{(0)} = {\bf b} - {\rm A} {\bf x}^{(0)}\) et \({\bf p}^{(0)} = {\bf r}^{(0)}\), la \(k\)-ième itération de la méthode du gradient conjugué s’écrit

\[\begin{array}{l} \displaystyle{\alpha_k = \frac{ ({\bf p}^{(k)},{\bf r}^{(k)}) }{ ({\bf p}^{(k)}, {\rm A} {\bf p}^{(k)})} } \\[5pt] {\bf x}^{(k+1)} = {\bf x}^{(k)} + \alpha_k {\bf p}^{(k)} \\[5pt] {\bf r}^{(k+1)} = {\bf r}^{(k)} - \alpha_k {\rm A} {\bf p}^{(k)} \\[5pt] \displaystyle{\beta_{k} = \frac{( {\rm A}{\bf p}^{(k)},{\bf r}^{(k+1)} ) } { ({\bf p}^{(k)}, {\rm A}{\bf p}^{(k)}) } } \\[5pt] {\bf p}^{(k+1)} = {\bf r}^{(k+1)} - \beta_{k} {\bf p}^{(k)} \end{array}\]

Concernant la convergence de la méthode du gradient conjugué, on a le résultat suivant:

Théorème: convergence de la méthode du gradient conjugué

Soit A une matrice symétrique définie positive d’ordre \(n\). Toute méthode qui utilise des directions conjuguées pour résoudre \(A\mathbf{x}=\mathbf{b}\) conduit à la solution exacte en au plus \(n\) itérations (en arithmétique exacte!).

Preuve

Les directions \({\bf p}^{(0)}, {\bf p}^{(1)}, \ldots, {\bf p}^{(n-1)}\) forment une base A-orthogonale de \(\mathbb{R}^n\). De plus, puisque \({\bf x}^{(n)}\) est optimal par rapport à toutes les directions \({\bf p}^{(j)}\), \(j=0, \ldots, n-1\), le vecteur \({\bf r}^{(n)}\) est orthogonal à l’espace \(S_{n-1} = \langle {\bf p}^{(0)}, {\bf p}^{(1)},\ldots, {\bf p}^{(n-1)} \rangle\). Par conséquent, \({\bf r}^{(n)} \perp S_{n-1} = \mathbb{R}^n\) et donc \({\bf r}^{(n)} = {\bf 0}\) ce qui implique \({\bf x}^{(n)} = {\bf x}\). \(\blacksquare\)

On démontre aussi le théorème suivant:

Théorème: estimation de l’erreur de la méthode du gradient conjugué

Soit \(A\) une matrice symétrique définie positive. La méthode du gradient conjugué pour la résolution de \(A\mathbf{x}=\mathbf{b}\) converge après au plus \(n\) étapes. De plus, l’erreur \({\bf e}^{(k)}\) à la \(k\)-ème itération (avec \(k < n\)) est orthogonale à \({\bf p}^{(j)}\), pour \(j=0, \ldots, k-1\) et

\[\label{errcg} \| {\bf e}^{(k)} \|_{\rm A} \leq \frac{2 c^k}{1+c^{2k}} \| {\bf e}^{(0)} \|_{\rm A}, \quad \mbox{avec } c = \frac{\sqrt{K_2 ({\rm A}) } - 1 } {\sqrt{K_2 ( {\rm A}) } + 1}.\]

La \(k\)-ième itération de la méthode du gradient conjugué n’est bien définie que si la direction de descente \({\bf p}^{(k)}\) est non nulle. En outre, si \({\bf p}^{(k)} = {\bf 0}\), alors l’itérée \({\bf x}^{(k)}\) doit coı̈ncider avec la solution x du système.

En l’absence d’erreur d’arrondi, la méthode du gradient conjugué peut donc être vue comme une méthode directe puisqu’elle converge en un nombre fini d’étapes Néanmoins, pour les matrices de grandes tailles, elle est généralement utilisée comme méthode itérative puisqu’on l’interrompt dès que l’erreur devient inférieure à une tolérance fixée. De ce point de vue, la dépendance de l’erreur par rapport au conditionnement de la matrice est plus favorable que pour la méthode du gradient. Signalons également que l’estimation [eq22] est souvent beaucoup trop pessimiste et ne prend pas en compte le fait que dans cette méthode, et contrairement à la méthode du gradient, la convergence est influencée par la totalité du spectre de A, et pas seulement par les valeurs propres extrémales.
Effet des erreurs d’arrondi

La méthode du gradient conjugué ne converge en un nombre fini d’étapes qu’en arithmétique exacte. L’accumulation des erreurs d’arrondi détruit l’A-orthogonalité des directions de descente et peut même provoquer des divisions par zéro lors du calcul des coefficients \(\alpha_k\) et \(\beta_k\). Ce dernier phénomène (appelé breakdown dans la littérature anglo-saxonne) peut être évité grâce à des procédés de stabilisation; on parle alors de méthodes de gradient stabilisées. Il se peut, malgré ces stratégies, que GC ne converge pas (en arithmétique finie) après \(n\) itérations. Dans ce cas, la seule possibilité raisonnable est de redémarrer les itérations avec le dernier résidu calculé. On obtient alors la méthode du gradient conjugué cyclique ou méthode du gradient conjugué avec redémarrage, qui ne possède pas les propriétés de convergence de GC.

1. La méthode du gradient conjugué préconditionné (Sec. 4.3.5 du livre)

Si \(P\) est une matrice de préconditionnement symétrique définie positive, la méthode du gradient conjugué préconditionné consiste à appliquer la méthode du gradient conjugué au système préconditionné

\[\begin{aligned} {\rm P}^{-1/2} {\rm A} {\rm P}^{-1/2} {\bf y} = {\rm P}^{-1/2} {\bf b}, \qquad \mbox{ avec } {\bf y} = {\rm P}^{1/2} {\bf x} \end{aligned}\]

En pratique, la méthode est implémentée sans calculer explicitement \(P^{1/2}\) ou \(P^{-1/2}\). Après un peu d’algèbre, on obtient le schéma suivant:

on se donne \({\bf x}^{(0)}\) et on pose \({\bf r}^{(0)} = {\bf b} - {\rm A}{\bf x}^{(0)}\), \({\bf z}^{(0)} = {\rm P}^{-1} {\bf r}^{(0)}\) et \({\bf p}^{(0)} = {\bf z}^{(0)}\), la \(k\)-ième itération s’écrit alors

\[\begin{array}{l} \displaystyle{\alpha_k = \frac{ ({\bf p}^{(k)}, {\bf r}^{(k)}) }{ ({\bf p}^{(k)},{\rm A} {\bf p}^{(k)})}}\\[5pt] {\bf x}^{(k+1)} = {\bf x}^{(k)} + \alpha_k {\bf p}^{(k)} \\[5pt] {\bf r}^{(k+1)} = {\bf r}^{(k)} - \alpha_k {\rm A} {\bf p}^{(k)} \\[5pt] {\rm P} {\bf z}^{(k+1)} = {\bf r}^{(k+1)} \\[5pt] \displaystyle{\beta_{k} = \frac{ ( {\rm A}{\bf p}^{(k)}, {\bf z}^{(k+1)}) } { ({\bf p}^{(k)}, {\rm A}{\bf p}^{(k)}) } } \\[5pt] {\bf p}^{(k+1)} = {\bf z}^{(k+1)} - \beta_{k} {\bf p}^{(k)} \end{array}\]
Le coût du calcul est plus élevé que pour GC puisqu’on doit résoudre à chaque itération le système linéaire \({\rm P} {\bf z}^{(k+1)} = {\bf r}^{(k+1)}\). L’estimation de l’erreur est la même que pour la méthode non préconditionnée, à condition de remplacer A par \({\rm P}^{-1} {\rm A}\).

2. Implémentation

Voici une implémentation de la méthode du gradient conjugué préconditionné pour la résolution d’un système linéaire \(A\mathbf{x}=\mathbf{b}\) où \(A\) est une matrice symétrique définie positive et \(P\) est une matrice de préconditionnement symétrique définie positive.

import numpy as np

def pcg(A, b, P=None, x0=None, tol=1e-6, maxiter=100):
    """
    Calcule la solution du système linéaire Ax = b en utilisant la méthode du gradient conjugué préconditionné.
    
    Paramètres :
    ----------
    A : ndarray, shape (n, n)
        Matrice des coefficients.
    
    b : ndarray, shape (n,)
        Vecteur de droite.
    
    P : ndarray, shape (n, n), optionnel
        Matrice de préconditionnement. Par défaut, il s'agit de la matrice d'identité.
    
    x0 : ndarray, shape (n,), optionnel
        Estimation initiale de la solution. Par défaut, il s'agit d'un vecteur nul.
    
    maxiter : int, optionnel
        Nombre maximal d'itérations pour le solveur. La valeur par défaut est 100.
    
    tol : float, optionnel
        Tolérance pour le critère de convergence de la norme résiduelle. La valeur par défaut est 1e-6.
    
    Retourne :
    -------
    x : ndarray, shape (n,)
        Solution approximative de Ax = b.
    
    niter : int
        Nombre d'itérations effectuées.

    inc : liste
        L'incrément à chaque itération.

    res : list
        La norme résiduelle à chaque itération.

    Notes :
    -----
    La méthode du gradient conjugué préconditionné utilise la matrice P pour accélérer la convergence. La méthode itère
    tant que la norme résiduelle est supérieure à la tolérance spécifiée ou jusqu'à ce que le nombre maximal d'itérations soit atteint.
    """
    if P is None:
        P = np.eye(len(A))
    if x0 is None:
        x0 = np.zeros(len(A))
    niter = 0
    x = x0
    r = b - A @ x
    r0 = r
    z = np.linalg.solve(P, r)
    p = z
    res = [1]
    inc = []
    #while res[-1] > tol and niter < maxiter:
    while niter < maxiter:
        niter += 1
        pold = p
        zold = z
        rold = r
        xold = x
        alpha = np.dot(p, rold) / np.dot(p, A @ p)
        x = xold + alpha * p
        r = rold - alpha * A @ p
        z = np.linalg.solve(P, r)
        beta = np.dot(z, A @ pold ) / np.dot(pold, A @ pold)
        p =  z - beta * pold
        res.append(np.linalg.norm(r)/np.linalg.norm(r0))
        inc.append(np.linalg.norm(x - xold))
        if inc[-1] < tol:
            break
    return x, niter, inc, res

3. Exemple

Reprenons la matrice symétrique définie positive qui nous a servie à illustrer l’algorithme du gradient.

import numpy as np
# pick a 4x4 matrix symétrique définie positive
A=np.array([[2.,-1.,0.,0.],[-1.,2.,-1.,0.],[0.,-1.,2.,-1.],[0.,0.,-1.,2.]])
print(f"A = {A}")
# on calcule les valeurs propres de A
print(f"Valeurs propres de A: {np.linalg.eigvals(A)}")
print(f"Conditionnement de A: {np.linalg.cond(A)}")
Results 
A = [[ 2. -1.  0.  0.]
 [-1.  2. -1.  0.]
 [ 0. -1.  2. -1.]
 [ 0.  0. -1.  2.]]
Valeurs propres de A: [3.61803399 2.61803399 0.38196601 1.38196601]
Conditionnement de A: 9.472135954999578

On peut maintenant résoudre le système \(A\mathbf{x}=\mathbf{b}\) avec la méthode du gradient conjugué préconditionné.

from tan.syslin.pcg import pcg
[x_pcg,niter_pcg,res_pcg,inc_pcg]=pcg( A, np.array([1.,0.,1.,0]), P=np.diag(np.diag(A)) )
print(f"Solution: x = {x_pcg}")
print(f"Nombre d'itérations: {niter_pcg}")
Results 
Solution: x = [1.2 1.4 1.6 0.8]
Nombre d'itérations: 5

Nous pouvons tracer la convergence du résidu ainsi que celle de l’incrément avec Plotly

import plotly.graph_objects as go
fig = go.Figure()
fig.add_trace(go.Scatter(x=np.arange(0,niter_pcg+1), y=inc_pcg, mode='lines+markers', name='Incrément'))
fig.add_trace(go.Scatter(x=np.arange(0,niter_pcg+1), y=res_pcg, mode='lines+markers', name='Résidu'))
fig.update_layout(title='Convergence de la méthode du gradient conjugué préconditionné', xaxis_title='Itération', yaxis_title='Incrément et Résidu', yaxis_type='log')
fig.show()
Results
nous observons que le résidu et l’incrément atteignent le 0 machine en 4 itérations. Cela est dû au fait que la matrice \(A\) est de taille 4 et que la méthode du gradient conjugué converge en au plus 4 itérations pour une matrice de taille 4, voir cette remarque.