Determinazione dell'autovalore minimo

La ricerca dell'autovalore minimo di $A$ e dell'autovettore ad esso collegato si può ricondurre alla minimizzazione del quoziente di Rayleigh in $\Re^n$. Se $A$ è una matrice simmetrica e definita positiva, il campo $q$ è una forma quadratica definita positiva del tutto simile a quella utilizzata nella ricerca della soluzione di un sistema lineare. Si può, pertanto, pensare di applicare lo schema del gradiente coniugato per minimizzare $q$.

Anche in questo caso l'esperienza insegna che lo schema del gradiente coniugato viene notevolmente accelerato se si utilizza una opportuna matrice di precondizionamento $K^{-1}$ per $A$. I criteri di scelta di $K^{-1}$ sono del tutto analoghi a quelli adottati nella soluzione di un sistema lineare. La pratica dimostra che l'uso della decomposta incompleta di Cholesky come precondizionatore consente di aumentare notevolmente la velocità di convergenza dello schema.

Il vettore $\bz$ che minimizza $q$ viene cercato costruendo una successione $\bz_k$ mediante la formula ricorrente:

$$\bz_{k+1} = \bz_k + \alpha_k \bp_k$$ (13)

Affinchè la successione dei vettori così generata converga effettivamente a $\bv_n$, il coefficiente scalare $\alpha_k$ viene calcolato in modo tale da minimizzare $q \left( \bz_{k+1} \right)$:

$$\frac{\partial}{\partial \alpha_k} \left[ \frac{\left( \bz_k + \a... ...] = \frac{\partial}{\partial \alpha_k} \left( \frac{\Phi_1}{\Phi_2} \right) = 0$$ (14)

Il calcolo esplicito della (14) risulta abbastanza laborioso. Conviene determinare separatamente le derivate rispetto ad $\alpha_k$ di $\Phi_1$ e $\Phi_2$:

$$\frac{\partial \Phi_1}{\partial \alpha_k} = \frac{\partial \Phi_1}{\partial \bz_{k+1}} \cdot \frac{\partial \bz_{k+1}}{\partial \alpha_k} = 2 \bp_k^T A \left( \bz_k + \alpha_k \bp_k \right)$$ (15)

$$\frac{\partial \Phi_2}{\partial \alpha_k} = \frac{\partial \Phi_2}{\partial \bz_{k+1}} \cdot \frac{\partial \bz_{k+1}}{\partial \alpha_k} = 2 \bp_k^T \left( \bz_k + \alpha_k \bp_k \right)$$ (16)

le quali possono essere scritte, per agevolare i calcoli, nel seguente modo:

$$\frac{\partial \Phi_1}{\partial \alpha_k} = 2 \left( a + \alpha_k b \right)$$ (17)

$$\frac{\partial \Phi_2}{\partial \alpha_k} = 2 \left( c + \alpha_k d \right)$$ (18)

mediante l'introduzione degli scalari $a$, $b$, $c$ e $d$ definiti come:

$$a = \bp_k^T A \bz_k$$ (19)

$$b = \bp_k^T A \bp_k$$ (20)

$$c = \bp_k^T \bz_k$$ (21)

$$d = \bp_k^T \bp_k$$ (22)

Si scriva ora in modo esplicito $\Phi_1$ e $\Phi_2$:

$$\Phi_1 = \bz_k^T A \bz_k + 2 \alpha_k \bp_k^T A \bz_k + \alpha_k^2 \bp_k^T A \bp_k$$ (23)

$$\Phi_2 = \bz_k^T \bz_k + 2 \alpha_k \bp_k^T \bz_k + \alpha_k^2 \bp_k^T \bp_k$$ (24)

che, grazie agli scalari:

$$e = \bz_k^T A \bz_k$$ (25)

$$m = \bz_k^T \bz_k$$ (26)

diventano:

$$\Phi_1 = e + 2 \alpha_k a + \alpha_k^2 b$$ (27)

$$\Phi_2 = m + 2 \alpha_k c + \alpha_k^2 d$$ (28)

Sviluppando ora formalmente l'equazione (14), si ottiene:

$$\frac{\Phi_2 \displaystyle{\frac{\partial \Phi_1}{\partial \alpha... ... \Phi_1 \displaystyle{\frac{\partial \Phi_2}{\partial \alpha_k}}}{\Phi_2^2} = 0$$ (29)

Introducendo nella (29) le (17), (18), (27) e (28), calcolate con l'uso delle (19)-(22), (25) e (26), si ricava con semplici calcoli:

$$\alpha_k^2 \left( bc - ad \right) + \alpha_k \left( bm - de \right) + \left( am - ce \right) = 0$$ (30)

da cui si ha infine lo scalare $\alpha_k$:

$${\alpha_k}_{1,2} = \frac{\left( de - bm \right) \pm \sqrt{\Delta}}{2 \left( bc - ad \right)}$$ (31)

avendo posto:

$$\Delta = \left( de - bm \right)^2 - 4 \left( bc - ad \right) \left( am - ce \right)$$ (32)

Si può facilmente verificare che per minimizzare $\Phi_1/\Phi_2$ è necessario scegliere nella (31) la soluzione con il segno ``+'' davanti alla radice quadrata.

Dopo aver calcolato la nuova approssimazione $\bz_{k+1}$ dell'autovettore $\bv_n$ mediante la (13), si può verificare se la convergenza dello schema è stata raggiunta. In base alla definizione (1) di autovalori ed autovettori, il vettore residuo $\br_{k+1}$ del GCM applicato nella minimizzazione di $q \left( \bz \right)$ si definisce come:

$$\br_{k+1} = A \bz_{k+1} - q \left( \bz_{k+1} \right) \bz_{k+1}$$ (33)

Come già osservato nella verifica di convergenza del GCM per la soluzione di sistemi lineari, conviene adimensionalizzare il problema effettuando il controllo sul residuo relativo anziché su quello assoluto definito in (33). Pertanto, lo schema va terminato quando risulta verificata la condizione:

$$r_r = \frac{\left\vert \br_{k+1} \right\vert}{\left\vert A \bz_{k+1} \right\vert} <$$ (34)

dove toll è un'opportuna tolleranza fissata arbitrariamente (ad esempio, $10^{-6}$).

Qualora la condizione (34) non sia verificata, è necessario aggiornare la direzione di ricerca $\bp_k$ e proseguire con lo schema iterativo. Come nel GCM per la soluzione di un sistema lineare il vettore $\bp_k$ era collegato al gradiente della forma quadratica da minimizzare (in quel caso coincidente con il residuo $\br_{k}$), così ora la nuova direzione di ricerca $\bp_{k+1}$ sarà funzione del gradiente $\bg_{k+1}$ di $q$ calcolato in $\bz_{k+1}$. La relazione ricorrente che definisce $\bp_{k+1}$ risulta pertanto:

$$\bp_{k+1} = K^{-1} \bg_{k+1} + \beta_{k} \bp_k$$ (35)

dove $K^{-1}$ è la matrice di precondizionamento (ad esempio, la decomposta incompleta di Cholesky), $\bg_{k+1}$ in virtù della (11) e della (33) viene calcolato come:

$$\bg_{k+1} = \frac{2 \br_{k+1}}{\bz_{k+1}^T \bz_{k+1}}$$ (36)

e $\beta_k$ è uno scalare determinato imponendo che il vettore $\bp_{k+1}$ sia $A$-ortogonale a $\bp_k$:

$$\beta_k = - \frac{\bg_{k+1}^T K^{-1} A \bp_k}{\bp_k^T A \bp_k}$$ (37)

Lo schema del GCM per determinare l'autovalore minimo e l'autovettore ad esso associato in una matrice simmetrica e definita positiva va inizializzato scegliendo un vettore arbitrario $\bz_0$, calcolando $\br_0$, $\bg_0$ e quindi $\bp_0 = K^{-1} \bg_0$. La nuova approssimazione dell'autovettore $\bv_n$ si calcola mediante la (13) con $\alpha_k$ definito dalla (31) in cui si prende il segno ``+'' davanti alla radice quadrata. Se è verificata la condizione (34) lo schema ha raggiunto la convergenza e pertanto $\bv_n = \bz_{k+1}$ e $\lambda_n = q \left( \bz_{k+1} \right)$. Altrimenti si calcola la nuova direzione di ricerca (35) mediante la (36) e la (37), e si torna alla (13) fino a che la condizione (34) non è verificata. Si può notare che il costo di una singola iterazione della procedura appena indicata è sostanzialmente dovuto a 3 prodotti matrice-vettore: $A \bp_k$, $A \bz_{k+1}$ e $K^{-1} \bg_{k+1}$ (il prodotto $A \bz_k$ è infatti disponibile dall'iterazione precedente), per l'implementazione dei quali si rimanda allo schema del GCM per la soluzione di un sistema lineare.

Si osservi che, se il vettore iniziale $\bz_0$ possiede componente nulla lungo $\bv_n$, vale a dire è ortogonale a $\bv_n$, la successione dei $\bz_k$ non può convergere a $\bv_n$ e quindi il $\lambda_n$ così determinato non è corretto. In pratica, gli errori di arrotondamento commessi nel corso della procedura fanno acquisire gradualmente ai vettori $\bz_k$ una componente diversa da 0 lungo $\bv_n$, permettendo quindi ugualmente la convergenza, anche se in tempi più lunghi, all'autovettore cercato. La probabilità che tale inconveniente si ripeta cambiando il vettore iniziale $\bz_0$ è praticamente nulla.