Richiami dello schema del GCM

Il metodo del Gradiente Coniugato Modificato (GCM) si basa sulla tecnica del Gradiente Coniugato (GC) sviluppata nei primi anni '50 da Hestenes e Stiefel per la soluzione di sistemi lineari con matrice simmetrica e definita positiva. Secondo il metodo del GC, la soluzione del generico sistema:

$$A \bx = \bb$$ (1)

si ottiene costruendo una successione di approssimazioni successive $\bx_k$ a partire da un vettore iniziale arbitrario $\bx_0$ secondo la formula ricorrente:

$$\bx_{k+1} = \bx_k + \alpha_k \bp_k$$ (2)

Lo scalare $\alpha_k$ viene scelto in modo da minimizzare una opportuna misura $\Phi$ dell'errore commesso assumendo $\bx_{k+1}$ come soluzione del sistema (1). Adottando come misura dell'errore la seguente forma quadratica:

$$\Phi (\bx_{k+1}) = \be_{k+1}^T A \be_{k+1} = \left( \bx - \bx_{k+1} \right)^T A \left( \bx - \bx_{k+1} \right)$$ (3)

si ottiene:

$$\alpha_k = \frac{\bp_k^T \br_k}{\bp_k^T A \bp_k}$$ (4)

dove $\br_k = \bb - A\bx_k$ è il residuo associato alla $k$-esima iterazione. Il nuovo residuo relativo alla soluzione $\bx_{k+1}$ si ricava facilmente dalla (2):

$$\br_{k+1} = \bb - A \left( \bx_k + \alpha_k \bp_k \right) = \br_k - \alpha_k A \bp_k$$ (5)

La direzione di ricerca $\bp_{k+1}$ utilizzata per calcolare la nuova approssimazione della soluzione viene determinata $A$-ortogonalizzando $\bp_{k+1}$ rispetto a $\bp_k$ secondo la formula ricorrente:

$$\bp_{k+1} = \br_{k+1} + \beta_k \bp_k$$ (6)

Imponendo la condizione di $A$-ortogonalità fra i vettori $\bp_k$ e $\bp_{k+1}$:

$$\bp_{k+1}^T A \bp_k = 0$$ (7)

si ottiene la relazione che definisce lo scalare $\beta_k$:

$$\beta_k = - \frac{\br_{k+1}^T A \bp_k}{\bp_k^T A \bp_k}$$ (8)

Lo schema del GC va inizializzato scegliendo un vettore soluzione arbitrario $\bx_0$, il corrispondente residuo $\br_0=\bb-A\bx_0$ e ponendo $\bp_0=\br_0$. Si applicano quindi in successione le formule ricorrenti (2), (5) e (6) mediante le (4) e (8) finché una opportuna norma del residuo, ad esempio la norma euclidea, risulta inferiore ad una prefissata tolleranza. Si osservi che è generalmente molto più conveniente fissare la tolleranza in funzione della norma del termine noto. Il residuo, quindi, su cui viene effettivamente eseguito il controllo sarà quello relativo:

$$r_r = \frac{\left\Vert \br_{k+1} \right\Vert}{\left\Vert \bb \right\Vert} < {\mbox {toll}}$$ (9)

rendendo adimensionale la tolleranza prescelta. L'ultima approssimazione $\bx^*$ del vettore $\bx$ così ottenuta viene assunta come soluzione del sistema lineare (1). Va osservato che, specialmente in problemi mal-condizionati, gli errori di arrotondamento che si accumulano nel calcolo del residuo mediante la formula ricorrente (5) possono far sì che quest'ultimo sia in realtà anche molto diverso dal residuo vero e che, di conseguenza, $\bx^*$ sia lontano dalla soluzione esatta del sistema (1). Per evitare tale inconveniente, è opportuno calcolare al termine della procedura iterativa del GC il residuo relativo vero $r_r^*$:

$$r_r^* = \frac{\left\Vert \bb-A\bx^* \right\Vert}{\left\Vert \bb \right\Vert} < {\mbox{toll}}$$ (10)

e verificare che sia effettivamente inferiore alla tolleranza prefissata.

Si può dimostrare che il vettore $\br_{k+1}$ risulta ortogonale ai precedenti $\br_0, \ldots, \br_k$ e che, pertanto, il residuo $n$-esimo, avendo indicato con $n$ la dimensione del sistema (1), deve necessariamente essere nullo. A causa degli errori di arrotondamento del calcolatore, tuttavia, ciò non si verifica ed il GC, pur essendo teoricamente un metodo diretto, va catalogato fra le tecniche iterative. In conseguenza a questo aspetto, che incide negativamente sulla performance del solutore, il GC, se non opportunamente accelerato, non offre alcun vantaggio rispetto ai metodi diretti classici.

La convergenza del GC risulta particolarmente accelerata qualora gli autovalori della matrice $A$ siano raggruppati attorno ad unico valore, il che corrisponde ad interpretare geometricamente l'iper-ellissoide associato alla matrice $A$ come quasi-sferico. Si può, pertanto, cercare di risolvere un nuovo sistema:

$$B \by = \bc$$ (11)

equivalente al precedente, in cui, tuttavia, la matrice $B$ presenta i propri autovalori raggruppati attorno all'unità. Il nuovo sistema (11), detto sistema precondizionato, è collegato al sistema originale (1) mediante le seguenti relazioni:

$$B = X^{-1} A X^{-1} \;\;\;\;\;\;\; \by = X \bx \;\;\;\;\;\;\; \bc = X^{-1} \bb$$ (12)

La matrice $K^{-1} = X^{-1} X^{-1}$ è detta matrice di precondizionamento. Se si sceglie $K^{-1} = A^{-1}$ si osserva facilmente che la matrice $B$ coincide con l'identità e la soluzione di (11) è immediata. È ovvio che tale opportunità è del tutto teorica, in quanto se si conoscesse già l'inversa di $A$ non servirebbe ricorrere ad alcun metodo per risolvere il sistema (1). Tuttavia, questa osservazione ci permette di concludere che la matrice di precondizionamento prescelta sarà tanto più efficace quanto più il prodotto $A K^{-1}$ si avvicinerà in qualche modo all'identità.

Riscrivendo le equazioni del GC per il sistema (11) e successivamente ripristinando le variabili originali, si ottiene lo schema del Gradiente Coniugato Modificato (GCM):

$$\bx_{k+1} = \bx_k + \alpha_k \bp_k$$

$$\br_{k+1} = \br_k - \alpha_k A \bp_k$$ (13)

$$\bp_{k+1} = K^{-1} \br_{k+1} + \beta_k \bp_k$$

dove:

$$\alpha_k = \frac{\bp_k^T \br_k}{\bp_k^T A \bp_k} \;\;\;\;\; \beta_k = - \frac{\br_{k+1}^T K^{-1} A \bp_k}{\bp_k^T A \bp_k}$$ (14)

Lo schema del GCM va inizializzato partendo da un vettore arbitrario $\bx_0$, con il relativo residuo $\br_0=\bb-A\bx_0$, e dalla direzione $\bp_0 = K^{-1} \br_0$. La soluzione iniziale può essere migliorata con alcune (poche) iterazioni effettuate con lo schema delle Correzioni Residue (CR), che coincide con lo schema di Jacobi applicato al sistema (11) ove si assuma la diagonale di $B$ unitaria:

$$\by_0 = \bc$$

$$\by_1 = \left( I - B \right) \by_0 + \bc$$

$$\vdots$$ (15)

$$\by_{k+1} = \left( I - B \right) \by_k + \bc$$

Espresso nelle variabili originali, si verifica facilmente che lo schema (15) risulta:

$$\bx_0 = K^{-1} \bb$$

$$\bx_1 = \bx_0 + K^{-1} \br_0$$

$$\vdots$$ (16)

$$\bx_{k+1} = \bx_k + K^{-1} \br_k$$

Lo schema del GCM si colloca, allo stato attuale, fra i metodi più moderni ed efficienti per la soluzione di sistemi lineari sparsi, simmetrici e definiti positivi, quali tipicamente si possono derivare dall'applicazione del metodo degli elementi finiti in svariati problemi di ingegneria civile, ambientale, meccanica, ecc. Uno dei punti vincenti del GCM sta nella semplicità di implementazione e nel ridotto costo computazionale rischiesto ad ogni iterazione, essendo le operazioni più dispendiose limitate a due prodotti matrice-vettore ($A \bp_k$ e $K^{-1} \br_{k+1}$). Tuttavia, per sfruttare appieno da un punto di vista computazionale le potenzialità di tale tecnica è opportuno memorizzare la matrice del sistema non in forma tabellare, bensì in forma compatta. Sarà pertanto necessario ri-programmare l'operazione di prodotto matrice-vettore secondo questa nuova notazione, come vedremo in seguito. A tali costi, va infine aggiunto quello per il calcolo iniziale della matrice di precondizionamento $K^{-1}$ e della sua applicazione nello schema (13)-(14).