Metodo del gradiente coniugato

In analisi numerica, il metodo del gradiente coniugato (spesso abbreviato in CG, dall'inglese conjugate gradient) è un algoritmo per la risoluzione numerica di un sistema lineare la cui matrice sia simmetrica e definita positiva.

Il metodo è stato inizialmente proposto nel 1952 dai matematici Magnus Hestenes e Eduard Stiefel^[1] e costituisce una variante del metodo del gradiente in cui la direzione di discesa a ogni passo è scelta in modo tale da garantire la convergenza del metodo in un numero di iterazioni pari al più alla dimensione del sistema da risolvere.

Il metodo del gradiente biconiugato ne fornisce una generalizzazione al caso di matrici non simmetriche.

Si voglia calcolare la soluzione ${\displaystyle \mathbf {x} \in \mathbb {R} ^{n}}$ del sistema lineare

${\displaystyle A\mathbf {x} =\mathbf {b} }$

dove ${\displaystyle A\in \mathbb {R} ^{n\times n}}$ è una matrice simmetrica e definita positiva a coefficienti reali e ${\displaystyle \mathbf {b} \in \mathbb {R} ^{n}}$ è il termine noto.

La matrice ${\displaystyle A}$, grazie alle sue proprietà, induce un prodotto scalare ${\displaystyle \langle \cdot ,\cdot \rangle _{A}:\mathbb {R} ^{n}\times \mathbb {R} ^{n}\to \mathbb {R} }$ definito da

${\displaystyle \langle \mathbf {u} ,\mathbf {v} \rangle _{\mathbf {A} }:=\mathbf {u} ^{T}\mathbf {A} \mathbf {v} ,\qquad \forall \mathbf {u} ,\mathbf {v} \in \mathbb {R} ^{n}.}$

Una coppia di vettori ${\displaystyle \mathbf {u} ,\mathbf {v} }$ che soddisfa ${\displaystyle \langle \mathbf {u} ,\mathbf {v} \rangle _{A}=0}$, cioè ortogonale rispetto a questo prodotto scalare, si dice ${\displaystyle A}$-coniugata.

Inoltre la soluzione ${\displaystyle \mathbf {x} }$ del sistema lineare precedente corrisponde al punto di minimo della forma quadratica

${\displaystyle Q(\mathbf {x} )={\frac {1}{2}}\mathbf {x} ^{T}A\mathbf {x} -\mathbf {x} ^{T}\mathbf {b} .}$

Infatti:

${\displaystyle \nabla Q(\mathbf {x} )=A\mathbf {x} -\mathbf {b} }$

da cui

${\displaystyle \nabla Q(\mathbf {x} )=0\Leftrightarrow A\mathbf {x} =\mathbf {b} .}$

Questo suggerisce di procedere iterativamente, partendo da una data soluzione iniziale ${\displaystyle \mathbf {x} _{0}}$ e muovendosi lungo direzioni ${\displaystyle \left\{\mathbf {p} _{k}\right\}_{k=0}^{n}}$ che minimizzano la forma quadratica ${\displaystyle Q(\mathbf {x} ).}$ A differenza del metodo del gradiente, in cui la direzione di discesa ${\displaystyle \mathbf {p} _{k}}$ al ${\displaystyle k}$-esimo passo è scelta pari a ${\displaystyle \mathbf {p} _{k}=-\nabla Q(\mathbf {x} _{k})}$, nel caso del gradiente coniugato essa viene scelta in modo che risulti ${\displaystyle A}$-ortogonale alle direzioni precedenti, cioè ${\displaystyle \langle \mathbf {p} _{j},\mathbf {p} _{k}\rangle _{A}=0,\,\forall j=0,\dots ,k-1.}$ Il significato geometrico di tale scelta è mostrato nella figura a lato, da cui emerge in particolare il vantaggio di scegliere direzioni ${\displaystyle A}$-ortogonali e non semplicemente ortogonali alle linee di livello della funzione ${\displaystyle Q(\mathbf {x} )}$.

Alla ${\displaystyle k}$-esima iterazione la soluzione viene dunque aggiornata nel modo seguente:

${\displaystyle \mathbf {x} _{k+1}=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k},}$

dove ${\displaystyle \alpha _{k}\in \mathbb {R} ^{+}}$ corrisponde alla lunghezza del passo di discesa. È possibile dimostrare (si veda ad esempio Soluzione di sistemi lineari) che la scelta ottimale per ${\displaystyle \alpha _{k}}$, che porta cioè al minimo di ${\displaystyle Q(\mathbf {x} _{k+1})}$, è

${\displaystyle \alpha _{k}={\frac {\mathbf {p} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}A\mathbf {p} _{k}}},}$

dove

${\displaystyle \mathbf {r} _{k}=\mathbf {b} -A\mathbf {x} _{k}}$

è il residuo del sistema.

Un metodo per calcolare direzioni di discesa ${\displaystyle A}$-ortogonali alle precedenti è il seguente^[2]:

${\displaystyle \mathbf {p} _{k+1}=\mathbf {r} _{k+1}-\beta _{k}\mathbf {p} _{k},}$

con ${\displaystyle \mathbf {p} _{0}=\mathbf {r} _{0}}$; la scelta ottimale per ${\displaystyle \beta _{k}}$ è

${\displaystyle \beta _{k}={\frac {\mathbf {p} _{k}^{T}A\mathbf {r} _{k+1}}{\mathbf {p} _{k}^{T}A\mathbf {p} _{k}}}.}$

Lo schema generale per la soluzione mediante metodo del gradiente coniugato è il seguente:

${\displaystyle {\begin{aligned}&\mathbf {x} _{0}=vettore\ iniziale\ arbitrario\\&\mathbf {r} _{0}=\mathbf {b} -A\mathbf {x} _{0}\\&\mathbf {p} _{0}=\mathbf {r} _{0}\\&{\text{for }}k=0,\dots ,n\\&\qquad \alpha _{k}={\frac {\mathbf {p} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}A\mathbf {p} _{k}}}\\&\qquad \mathbf {x} _{k+1}=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}\\&\qquad \mathbf {r} _{k+1}=\mathbf {b} -A\mathbf {x} _{k+1}\\&\qquad \beta _{k}={\frac {\mathbf {p} _{k}^{T}A\mathbf {r} _{k+1}}{\mathbf {p} _{k}^{T}A\mathbf {p} _{k}}}\\&\qquad \mathbf {p} _{k+1}=\mathbf {r} _{k+1}-\beta _{k}\mathbf {p} _{k}\\&\qquad k=k+1\\&{\hbox{end}}.\end{aligned}}}$

L'eventuale implementazione dell'algoritmo in aritmetica floating point, in cui la convergenza in al più ${\displaystyle n}$ passi non è garantita, il ciclo for può essere sostituito da un ciclo while che verrà eseguito finché la norma del residuo ${\displaystyle \|\mathbf {r} _{k}\|}$ non sia più piccola di una tolleranza impostata dall'utente.

In molti casi è possibile accelerare ulteriormente la velocità di convergenza dell'algoritmo migliorando le proprietà di condizionamento della matrice ${\displaystyle A}$. Si introduca a tal fine una matrice di precondizionamento ${\displaystyle P}$ simmetrica e definita positiva. L'algoritmo corrispondente al metodo del gradiente coniugato precondizionato (spesso abbreviato in PCG, dall'inglese preconditioned conjugate gradient) si ottiene applicando la versione senza precondizionamento per trovare la soluzione ${\displaystyle {\hat {\mathbf {x} }}}$ del seguente sistema:

${\displaystyle R^{-1}AR^{-1}{\hat {\mathbf {x} }}=R^{-1}\mathbf {b} }$,

dove ${\displaystyle R}$ è la radice quadrata di ${\displaystyle P}$ e ${\displaystyle {\hat {\mathbf {x} }}=R\mathbf {x} }$.

Lo schema risolutivo in questo caso diventa^[2]:

${\displaystyle {\begin{aligned}&{\text{for }}k=0,\dots ,n\\&\qquad \alpha _{k}={\frac {\mathbf {p} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}A\mathbf {p} _{k}}}\\&\qquad \mathbf {x} _{k+1}=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}\\&\qquad \mathbf {r} _{k+1}=\mathbf {b} -A\mathbf {x} _{k+1}\\&\qquad {\text{trovare la soluzione }}\mathbf {z} _{k+1}{\text{ del sistema }}P\mathbf {z} _{k+1}=\mathbf {r} _{k+1}\\&\qquad \beta _{k}={\frac {\mathbf {p} _{k}^{T}A\mathbf {z} _{k+1}}{\mathbf {p} _{k}^{T}A\mathbf {p} _{k}}}\\&\qquad \mathbf {p} _{k+1}=\mathbf {z} _{k+1}-\beta _{k}\mathbf {p} _{k}\\&\qquad k=k+1\\&{\hbox{end}}.\end{aligned}}}$

È possibile dimostrare che l'errore commesso alla ${\displaystyle k}$-esima iterazione del metodo del gradiente coniugato soddisfa la seguente stima^[2]:

${\displaystyle \|\mathbf {e} _{k}\|_{A}\leq {\frac {2c^{k}}{1+c^{2k}}}\|\mathbf {e} _{0}\|_{A},}$

dove

${\displaystyle c={\frac {{\sqrt {\kappa (A)}}-1}{{\sqrt {\kappa (A)}}+1}},}$

${\displaystyle \kappa (A)}$ il numero di condizionamento in norma ${\displaystyle 2}$ di ${\displaystyle A}$ e ${\displaystyle \|\mathbf {x} \|_{A}:=\langle \mathbf {x} ,\mathbf {x} \rangle _{A}}$ è la norma indotta da ${\displaystyle A}$.

Nel caso precondizionato vale la stessa stima con

${\displaystyle c={\frac {{\sqrt {\kappa (P^{-1}A)}}-1}{{\sqrt {\kappa (P^{-1}A)}}+1}}.}$

Si riporta un esempio di possibile implementazione del metodo del gradiente coniugato non precondizionato compatibile con i linguaggi di programmazione Octave e MATLAB.

function [xk, iter] = gradiente_coniugato(A, b, x0, toll, nmax)
    xk = x0;        
    rk = b - A * xk;
    pk = rk;
    iter = 0;
    while (norm(rk) >= toll*norm(b))
        alphak = (pk' * rk) / (pk' * A * pk);
        xk = xk + alphak * pk;
        rk = b - A * xk;
        betak = (pk' * A * rk) / (pk' * A * pk);
        pk = rk - betak * pk;
        iter = iter+1;
      if (iter == nmax && norm(rk) > toll*norm(b)) 
        disp(['warning: Convergenza non raggiunta in ' num2str(iter) ' iterazioni!']);
        break
      end
    end
end

La funzione che implementa il metodo del gradiente coniugato precondizionato è già salvata in MATLAB nel comando pcg().
Esempio:

x=pcg(A,b) 
%determina la soluzione x del sistema lineare Ax=b di una matrice simmetrica e definita positiva mediante il metodo del gradiente coniugato a partire dal vettore iniziale x0 nullo.

x=pcg(A,b,tol,nmax)
%determina la soluzione x imponendo come criterio d'arresto la tolleranza e il numero di iterazioni.

• A. Quarteroni, R. Sacco e F. Saleri, Matematica numerica, 4ª ed., Milano, Springer Verlag, 13 marzo 2014, ISBN 978-88-470-5643-5.

• Metodo del gradiente
• Algoritmo di Gauss-Newton

• Il metodo del Gradiente Coniugato (PDF), su ing.unibs.it. URL consultato il 17 giugno 2016.
• Metodo del gradiente coniugato per la soluzione di sistemi lineari (PDF), su bugs.unica.it. URL consultato il 17 giugno 2016.
• (EN) An Introduction to the Conjugate Gradient Method Without the Agonizing Pain (PDF), su cs.cmu.edu. URL consultato il 17 giugno 2016.

Portale Matematica: accedi alle voci di Teknopedia che trattano di matematica