Metodo del gradiente coniugato

Metodo CG §

Il metodo del gradiente coniugato (Hestenes-Stiefel 1950) costruisce la nuova direzione di discesa come combinazione lineare della direzione precedente e del residuo.

Pseudo codice

Direzione di ricerca iniziale $p_0=r_0$ for $k=0,1,\dots$ calcolo step ${\alpha }_k=\frac{r_k^T{r}_k}{p_k^{T}A{p}_k}$ update soluzione $x_{k+1}=x_k+{\alpha }_k{p}_k$ update residuo $r_{k+1}=r_k-{\alpha }_{k}A{p}_k$ calcolo step ${\beta }_{k+1}=\frac{r_{k+1}^T{r}_{k+1}}{r_k^T{r}_k}$ update direzione discesa $p_{k+1}=r_{k+1}+{\beta }_{k+1}{p}_k$ end

Osservazioni Le quantità $r_k^T{r}_k$ e $A{p}_k$ compaiono due volte, pertando conviene memorizzarle per evitare calcoli unitili (il residuo compare anche nell’iterazione successiva).

Il calcolo del passo ${\alpha }_k$ dinamico segue sempre una minimzation rule, infatti:

$${\alpha }_k=\frac{r_k^T{p}_k}{p_k^{T}A{p}_k}=\frac{r_k^T(r_k+{\beta }_k{p}_{k-1}))}{p_k^{T}A{p}_k}=\frac{r_k^T{r}_k}{p_k^{T}A{p}_k}$$

siccome per l’osservazione di ortogonalità tra i residui e le direzioni di ricerca precedenti $r_k^T{p}_{k-1}=0$.

Si può dare una caratterizzazione equivalente del passo ${\beta }_k$. Lemma Se si ha, per $k\ge 1$, $r_k^T{r}_{k+1}=0$ e, per $k>1$, $p_{k-2}^{T}A{p}_{k-1}=0$, allora vale

$${\beta }_k=-\frac{r_k^{T}A{r}_{k-1}}{p_{k-1}^{T}A{p}_{k-1}}\text{ per }k\ge 1$$

Tuttavia questa formulazione è computazionalmente meno conveniente.

Enuciamo ora un teorema fondamentale che caratterizza il metodo CG.

Teorema Siano stati compiuti $k$ passi, allora vale la seguente indentità tra sottospazi:

$$span\{p_0,\dots ,p_j\}=span\{r_0,\dots ,r_j\}\text{ per }0<j\le k$$

inoltre i residuo sono mutualmente ortogonali, e le direazioni di discesa sono $A$-coniugate (ortogonali):

$$r_k^T{r}_j=0p_k^{T}A{p}_j=0\text{ per }j<k$$

Dim La prima identità tra gli span è immediata, segue dalla definizione della direzione di discesa:

$$p_j=r_j+{\beta }_j{p}_{j-1}1\le j\le k{p}_0=r_0$$

Per dimostrare l’ortogonalità procediamo per induzione su $k$. La base $k=1$ è vera, infatti:

$$r_1^T{r}_0=r_1^T{p}_0=0$$

siccome lo step è ottimale, mentre

$$p_1^{T}A{p}_0=(r_1+{\beta }_1{p}_0{)}^{T}A{p}_0$$

e usando l’altra caratterizzazione di ${\beta }_k$ per $k=1$ risulta nullo. Assumiamo i risultato veri al passo $k-1$, con $k>1$, mostriamo l’ortgonoalità dei residui:

$$r_k^T{r}_j=(r_{k-1}-{\alpha }_{k-1}A{p}_{k-1}{)}^T{r}_j=r_{k-1}^T{r}_j-{\alpha }_{k-1}{p}_{k-1}^{T}A{r}_j=0$$

infatti per induzione, se $j<k-1$, il primo addendo è nullo, il secondo segue sempre dall’ipotesi induttiva e dall’uguaglianza degli span.

Resta il caso $j=k-1$, qui possiamo usare la definizione di ${\alpha }_{k-1}$:

$$r_{k-1}^T{r}_{k-1}-\frac{r_{k-1}^T{r}_{k-1}}{p_{k-1}^{T}A{p}_{k-1}}{p}_{k-1}^{T}A{r}_{k-1}=0$$

Infatti il denominatore si può riscrivere come

$$p_{k-1}^{T}A{p}_{k-1}=p_{k-1}^{T}A(r_{k-1}+{\beta }_{k-1}{p}_{k-2})$$

e per ipotesi induttiva l’addendo con ${\beta }_{k-1}$ è nullo.

Resta da dimostrare che le direzioni di discesa sono $A$-coniugate.

$$p_k^{T}A{p}_j=(r_k+{\beta }_k{p}_{k-1}{)}^{T}A{p}_j=r_k^{T}A{p}_j+{\beta }_k{p}_{k-1}^{T}A{p}_j$$

Se $j<k-1$ per ipotesi induttiva il secondo termine è nullo, il primo sempre per ipotesi induttiva e per l’uguaglianza degli span. Per dimostrae il caso $j=k-1$, si utilizza la definzione di ${\beta }_k$ del lemma:

$$r_k^{T}A{p}_{k-1}+{\beta }_k{p}_{k-1}^{T}A{p}_{k-1}=r_k^{T}A{p}_{k-1}-\frac{r_k^{T}A{p}_{k-1}}{p_{k-1}^{T}A{p}_{k-1}}{p}_{k-1}^{T}A{p}_{k-1}=0\square$$

Una conseguenza fondametale di questo teorema, è che (al netto di errori di calcolo) il metodo CG è un metodo diretto!

Teorema Sia $A\in {\mathbb{R}}^{n\times n}$ simmetrica e definita positiva, applicando il metodo CG al problema $Ax=b$, a partire da una qualsiasi solizione $x_0$, si ha convergenza in al più $n$ passi alla soluzione esatta.

Dim Segue direttamente dal teorema precedente, infatti $r_n=0$ siccome $span\{r_0,\dots ,r_{n-1}\}={\mathbb{R}}^n$ e $r_n^T{r}_j=0$ per $j<n$. Quindi $x_n=x$, ovvero la tesi. $\square$

Velocità di convergenza §

Si può dimostrare che il fattore di abbattimento dell’errore del metodo CG è

$$\frac{2{\left(\frac{\sqrt{k_2(A)}-1}{\sqrt{k_2(A)}+1}\right)}^{k+1}}{1+{\left(\frac{\sqrt{k_2(A)}-1}{\sqrt{k_2(A)}+1}\right)}^{2(k+1)}}\approx 2{\left(\frac{\sqrt{k_2(A)}-1}{\sqrt{k_2(A)}+1}\right)}^{k}k>>1$$

quindi migliore di quello di Steepest descent.

Metodo precondizionato PCG §

Per migliorare il condizionamento del sistema $Ax=b$ possiamo usare le tecniche di precondizionamento che preservino la proprietà di $A$ di simmetria e definita positiva, ad esempio con la fattorizzazione di Cholesky incompleta:

$$L_{inc}^{-1}A{L}_{inc}^{-T}x=L_{inc}^{-1}b$$

Tuttavia quando si implementa, non si calcola esplicitamente il nuovo sistema lineare precondizionato, ma si lavoro su di esso implicitamente. Per illustrare questa strategia, ragioniamo con il prodotto scalare indotto da un precondizionarore $P$ simmetrico e definitio positivo; valgono le seguenti osservazioni: Osservazione Se $A$ è simmetrica e definita positiva, allora $P^{-1}A$ risulta simmetrica e definita positiva rispetto al $P$-prodotto scalare. Dim Segue dalla definizione:

$$(P^{-1}Ax,y{)}_P=x^T{A}^T{P}^{-T}Py=x^{T}Ay=x^{T}P{P}^{-1}Ay=(x,P^{-1}Ay{)}_P$$ $$(P^{-1}Ax,x{)}_P=x^{T}Ax>0\forall x\ne 0$$

Inoltre si sta minimizzando il medesimo funzionale:

$${\Phi }_P(x)=\frac{1}{2}{x}^{T}P(P^{-1}A)x-x^{T}P{P}^{-1}b=\Phi (x)$$

Quindi si lavoro con il $P$-prodotto scalare sul sistema $P^{-1}A$.

Come si modifica lo pseudocodice?

Calcolo il residuo del sistema precondizionato $P^{-1}Ax=P^{-1}b$ $z=P^{-1}x$. Se usiamo $P=L_{inc}{L}_{inc}^T$ possiamo risolverlo usando forward e backward.

Calcolo dello step ottimale ${\alpha }_k=\frac{(z_k,z_k{)}_P}{(p_k,P^{-1}A{p}_k{)}_P}$ svolgendo i conti viene fuori:

$${\alpha }_k=\frac{r_k^T{z}_k}{p_k^{T}A{p}_k}$$

calcolo step per dir discesa: Anche qui risulta

$${\beta }_{k+1}=\frac{r_{k+1}^T{z}_{k+1}}{r_k^T{z}_k}$$

Direzione di ricerca iniziale $p_0=z_0$

for $k=0,1,\dots$ calcolo step ${\alpha }_k=\frac{r_k^T{r}_k}{p_k^{T}A{p}_k}$ update soluzione $x_{k+1}=x_k+{\alpha }_k{p}_k$ update residuo $r_{k+1}=r_k-{\alpha }_{k}A{p}_k$ calcolo step ${\beta }_{k+1}=\frac{r_{k+1}^T{r}_{k+1}}{r_k^T{r}_k}$ update direzione discesa $p_{k+1}=r_{k+1}+{\beta }_{k+1}{p}_k$ end

Criterio di arresto §

Generalmente si usa il rapporto tra la norma euclidea del residuo con la norma euclidea del termine noto

$$\frac{\Vert r_k{\Vert }_2}{\Vert b{\Vert }_2}<\text{ tol }$$

non il residuo precondizionato! Non vale la pena calcolare la norma dell’energia del residuo ad ogni passo.

Interpretazione come metodo di Krylov §

Il metodo CG può essere visto come un metodo di Krylov in cui la soluzione al passo $k$ è scelta per minimizzare la norma dell’energia dell’errore.

Teorema Sia $A\in {\mathbb{R}}^{n\times n}$, simmetrica definita positiva (SPD), applicando il metodo CG con $x_0=0$, la soluzione approssimata $x_k$ risulta l’unico vettore che rende minima la norma dell’energia dell’errore, ovvero:

$$\Vert e_k{\Vert }_A=\underset{v\in K_k(A,r_0)}{\min }\Vert x-v{\Vert }_A$$

Dim Supponiamo che $\hat{x}=x_k-\Delta x$ sia la soluzione al problema di minimizzazione, Facciamo vedere che $\Delta x=0$. Calcoliamo la $A$-norma dell’errore:

$$\Vert \hat{e}{\Vert }_A^2=(e_k+\Delta x{)}^{T}A(e_k+\Delta x)=e_{k}TA{e}_k+\Delta x^{T}A\Delta x+2e_k^{T}A\Delta x$$

ma $A{e}_k=r_k$ ed è ortogonale a $\Delta x\in K_k(A,r_0)$. Ma $A$ è definita positiva, dunque:

$$\Vert \hat{e}{\Vert }_A^2\ge \Vert e_k{\Vert }_A^2$$

con uguaglianza se e solo se $\Delta x=0$, ovvero $\hat{x}=x_k$. $\square$