Algoritmo degli autovalori di Jacobi

In algebra lineare, l'algoritmo degli autovalori di Jacobi è un metodo iterativo per il calcolo degli autovalori e degli autovettori di una matrice simmetrica reale (processo noto come diagonalizzazione). Prende il nome da Carl Gustav Jacob Jacobi, che per primo propose il metodo nel 1846, ma divenne ampiamente utilizzato solo negli anni '50 con l'avvento dei computer.

Descrizione

Sia $S$ una matrice simmetrica e $G=G(i,j,\theta )$ una matrice di rotazione di Givens. Si dimostra che la matrice

$S'=GSG^{\top }\,$

è simmetrica e simile a $S$ .

Inoltre, $S^{\prime }$ ha i seguenti valori:

${\begin{aligned}S'_{ii}&=c^{2}\,S_{ii}-2\,sc\,S_{ij}+s^{2}\,S_{jj}\\S'_{jj}&=s^{2}\,S_{ii}+2sc\,S_{ij}+c^{2}\,S_{jj}\\S'_{ij}&=S'_{ji}=(c^{2}-s^{2})\,S_{ij}+sc\,(S_{ii}-S_{jj})\\S'_{ik}&=S'_{ki}=c\,S_{ik}-s\,S_{jk}&k\neq i,j\\S'_{jk}&=S'_{kj}=s\,S_{ik}+c\,S_{jk}&k\neq i,j\\S'_{kl}&=S_{kl}&k,l\neq i,j\end{aligned}}$

dove $s=\sin(\theta )$ e $c=\cos(\theta )$.

Essendo $G$ ortogonale, $S$ E $S^{\prime }$ hanno la stessa norma di Frobenius $||\cdot ||_{F}$ (la somma delle radici quadrate dei quadrati di tutti gli elementi). In ogni caso, è possibile scegliere $\theta$ tale che $S_{ij}^{\prime }=0$, in tal caso $S^{\prime }$ ha una somma maggiore dei quadrati sulla diagonale:

$S'_{ij}=\cos(2\theta )S_{ij}+{\tfrac {1}{2}}\sin(2\theta )(S_{ii}-S_{jj})$

Imponendo $S'_{ij}$ uguale a 0, e riorganizzando i termini, si ottiene:

$\tan(2\theta )={\frac {2S_{ij}}{S_{jj}-S_{ii}}}$

Se $S_{jj}=S_{ii}$, allora

$\theta ={\frac {\pi }{4}}$

Per ottimizzare questo effetto, $S_{ij}$ dovrebbe essere l'elemento fuori diagonale con il valore assoluto più grande, chiamato pivot .

Il metodo degli autovalori di Jacobi esegue ripetutamente rotazioni finché la matrice non diventa quasi diagonale. Iterando questo processo, gli elementi nella diagonale sono approssimazioni degli autovalori (reali) di S.

Convergenza

Se $p=S_{kl}$ è un elemento pivot, allora per definizione $|S_{ij}|\leq |p|$ per $1\leq i,j\leq n,i\neq j$ .

Si definisce $\Gamma (S)^{2}$ come la somma dei quadrati di tutti gli elementi fuori diagonale di $S$ . Essendo che $S$ ha esattamente $2N:=n(n-1)$ elementi fuori diagonale, si ha che $p^{2}\leq \Gamma (S)^{2}\leq 2Np^{2}$ O $2p^{2}\geq \Gamma (S)^{2}/N$ .

Di conseguenza $\Gamma (S^{J})^{2}=\Gamma (S)^{2}-2p^{2}$ . Questo implica $\Gamma (S^{J})^{2}\leq (1-1/N)\Gamma (S)^{2}$ oppure $\Gamma (S^{J})\leq (1-1/N)^{1/2}\Gamma (S)$; in altre parole, la sequenza delle rotazioni di Jacobi converge almeno linearmente di un fattore $(1-1/N)^{1/2}$ ad una matrice diagonale.

Un numero di $N$ rotazioni di Jacobi sono chiamate spazzate; si definisce con $S^{\sigma }$ il risultato delle $N$ spazzate. La stima precedente rende

$\Gamma (S^{\sigma })\leq \left(1-{\frac {1}{N}}\right)^{N/2}\Gamma (S)$ ;

Ovvero, la sequenza di spazzate converge almeno linearmente con un fattore ≈ $e^{1/2}$ .

Tuttavia il seguente risultato di Schönhage fornisce una convergenza localmente quadratica. A tal fine sia S la presenza di m autovalori distinti $\lambda _{1},...,\lambda _{m}$ con molteplicità $\nu _{1},...,\nu _{m}$, e sia d > 0 la distanza più piccola tra due autovalori diversi. Chiamiamo un numero di

$N_{S}:={\frac {n(n-1)}{2}}-\sum _{\mu =1}^{m}{\frac {1}{2}}\nu _{\mu }(\nu _{\mu }-1)\leq N$

rotazioni di Jacobi una spazzata di Schönhage. Se $S^{s}$ denota il risultato, allora

$\Gamma (S^{s})\leq {\sqrt {{\frac {n}{2}}-1}}\left({\frac {\gamma ^{2}}{d-2\gamma }}\right),\quad \gamma :=\Gamma (S)$ .

Quindi la convergenza diventa quadratica non appena $\Gamma (S)<{\frac {d}{2+{\sqrt {{\frac {n}{2}}-1}}}}$

Costo

Ogni rotazione di Jacobi può essere eseguita in O(n) passi quando l'elemento pivot p è noto. Tuttavia la ricerca di p richiede l'ispezione di tutti gli N ≈ 1/2 n ² elementi fuori diagonale. Possiamo ridurre anche questo procedimento alla complessità O(n) se introduciamo un array di indici aggiuntivo $m_{1},\,\dots \,,\,m_{n-1}$, con la proprietà che $m_{i}$ è l'indice dell'elemento più grande nella riga i, (i = 1, ..., n − 1) dell'attuale S. Allora, gli indici del pivot ( k, l ) devono essere una delle coppie $(i,m_{i})$.

Anche l'aggiornamento dell'array dell'indice può essere eseguito con complessità nel caso medio O(n): in primo luogo, l'elemento maggiore nelle righe aggiornate k e l può essere trovato nei passaggi O(n). Nelle altre righe i cambiano solo gli elementi nelle colonne k e l. Iterando su queste righe, se $m_{i}$ non è né k né l, è sufficiente confrontare il vecchio massimo in $m_{i}$ con i nuovi elementi, aggiornando $m_{i}$ se necessario. Se $m_{i}$ dovrebbe essere uguale a k o ad l, e l'elemento corrispondente è diminuito durante l'aggiornamento, il massimo sulla riga i deve essere trovato da zero nella complessità O(n). Tuttavia, ciò avverrà in media solo una volta per rotazione. Pertanto, ogni rotazione ha complessità nel caso medio O( n ), e ogni spazzata ha complessità media O (n ³), che equivale a una moltiplicazione di matrici. Inoltre il $m_{i}$ deve essere inizializzato prima dell'avvio del processo, operazione che può essere eseguita in n ² passaggi.

Tipicamente, il metodo di Jacobi converge alla precisione numerica dopo un piccolo numero di passaggi. Si noti che più autovalori riducono il numero di iterazioni, in quanto $N_{S}<N$ .

Algoritmo

Il seguente algoritmo è una descrizione del metodo Jacobi in notazione di tipo matematico. Questo algoritmo calcola un vettore e che contiene gli autovalori, e una matrice E che contiene i corrispondenti autovettori; di conseguenza, ogni $e_{i}$ è un autovalore e ogni colonna $E_{i}$ un autovettore ortonormale per $e_{i}$, io = 1, ..., n .

Appunti

1. L'array logico changed mantiene lo stato di ciascun autovalore. Se il valore numerico di $e_{k}$ o $e_{l}$ viene modificato durante un'iterazione, il componente corrispondente di changed viene impostato su true, altrimenti su false. La variabile intera state conta il numero di componenti di changed che hanno valore true . L'iterazione si interrompe non appena state è uguale a 0. Ciò significa che nessuna delle approssimazioni $e_{1},\,...\,,e_{n}$ ha recentemente cambiato il suo valore, e quindi non è molto probabile che ciò accada se l'iterazione continua. Qui si presuppone che le operazioni in virgola mobile vengano arrotondate in modo ottimale al numero in virgola mobile più vicino.

2. Il triangolo superiore della matrice S viene distrutto, mentre il triangolo inferiore e la diagonale rimangono invariati. Pertanto, è possibile ripristinare S se necessario, seguendo il seguente algoritmo:

3. Gli autovalori non sono necessariamente in ordine decrescente. Ciò può essere ottenuto mediante un semplice algoritmo di ordinamento.

4. L'algoritmo è scritto utilizzando la notazione matriciale (array basati su 1 invece che su 0).

5. Quando si implementa l'algoritmo, la parte specificata utilizzando la notazione matriciale deve essere eseguita simultaneamente.

6. Questa implementazione non tiene conto correttamente del caso in cui una dimensione è un sottospazio indipendente. Ad esempio, nel caso in cui questo algoritmo analizzi una matrice diagonale, l'implementazione mostrata in precedenza non terminerà mai, poiché nessuno degli autovalori cambierà. Pertanto, nelle implementazioni reali, è necessario aggiungere ulteriore logica per tenere conto di questo caso.

Esempio

Sia $S={\begin{pmatrix}4&-30&60&-35\\-30&300&-675&420\\60&-675&1620&-1050\\-35&420&-1050&700\end{pmatrix}}$

Quindi jacobi produce i seguenti autovalori e autovettori dopo 3 scansioni (19 iterazioni) :

$e_{1}=2585.25381092892231$

$E_{1}={\begin{pmatrix}0.0291933231647860588\\-0.328712055763188997\\0.791411145833126331\\-0.514552749997152907\end{pmatrix}}$

$e_{2}=37.1014913651276582$

$E_{2}={\begin{pmatrix}-0.179186290535454826\\0.741917790628453435\\-0.100228136947192199\\-0.638282528193614892\end{pmatrix}}$

$e_{3}=1.4780548447781369$

$E_{3}={\begin{pmatrix}-0.582075699497237650\\0.370502185067093058\\0.509578634501799626\\0.514048272222164294\end{pmatrix}}$

$e_{4}=0.1666428611718905$

$E_{4}={\begin{pmatrix}0.792608291163763585\\0.451923120901599794\\0.322416398581824992\\0.252161169688241933\end{pmatrix}}$

Applicazioni per matrici reali simmetriche

Quando si conoscono gli autovalori (e gli autovettori) di una matrice simmetrica, si calcolano facilmente i seguenti valori.

Valori singolari

I valori singolari di una matrice (quadrata). $A$ sono le radici quadrate degli autovalori (non negativi) di $A^{T}A$ . Nel caso di una matrice simmetrica $S$ si ha che $S^{T}S=S^{2}$, quindi i valori singolari di $S$ sono i valori assoluti degli autovalori di $S$.

Norma di ordine 2 e raggio spettrale

La norma di ordine 2 di una matrice A è la norma basata sulla norma vettoriale euclidea; cioè, il valore più grande $\|Ax\|_{2}$ quando x attraversa tutti i vettori con $\|x\|_{2}=1$ . È il valore singolare più grande di $A$ . Nel caso di una matrice simmetrica, è anche il valore assoluto più grande dei suoi autovettori, e quindi uguale al suo raggio spettrale.

Numero di condizione

Il numero di condizione di una matrice non singolare $A$ è definito come ${\mbox{cond}}(A)=\|A\|_{2}\|A^{-1}\|_{2}$ . Nel caso di una matrice simmetrica, è il valore assoluto del quoziente tra l'autovalore maggiore e quello minore. Matrici con numeri di condizione elevati possono causare risultati numericamente instabili: piccole perturbazioni possono causare grandi errori. Le matrici di Hilbert sono le matrici mal condizionate più famose. Ad esempio, la matrice di Hilbert del quarto ordine ha una condizione pari a 15514, mentre per l'ordine 8 è 2,7 × 10 ⁸ .

Rango

Una matrice $A$ ha rango $r$ se ha $r$ colonne linearmente indipendenti, mentre le restanti colonne sono linearmente dipendenti da queste. Equivalentemente, $r$ è la dimensione dell'intervallo di $A$ . Inoltre è il numero di valori singolari diversi da zero.

Nel caso di una matrice simmetrica, r è il numero di autovalori diversi da zero. Sfortunatamente, a causa di errori di arrotondamento, le approssimazioni numeriche di autovalori nulli potrebbero non essere zero (può anche succedere che un'approssimazione numerica sia nulla, mentre il valore reale non sia effettivamente 0). Pertanto è possibile calcolare il rango numerico solo decidendo quali autovalori sono sufficientemente vicini a zero.

Pseudo-inversa della matrice

La pseudo-inversa di una matrice $A$ è la matrice unica $X=A^{+}$ per cui $AX$ E $XA$ sono simmetriche e per cui la condizione$AXA=A,XAX=X$ è verificata.

Se $A$ non è singolare, quindi $A^{+}=A^{-1}$ .

Quando viene chiamata la procedura jacobi(S, e, E), allora la relazione $S=E^{T}{\mbox{Diag}}(e)E$ vale dove Diag(e) denota la matrice diagonale con il vettore e sulla diagonale. Si definisce con $e^{+}$ il vettore dove $e_{i}$ è sostituito da $1/e_{i}$ se $e_{i}\leq 0$, e da 0 se $e_{i}$ è (numericamente vicino a) zero. Poiché la matrice E è ortogonale, ne consegue che la pseudo-inversa di S è data da $S^{+}=E^{T}{\mbox{Diag}}(e^{+})E$ .

Soluzione dei minimi quadrati

Se matrice $A$ non ha rango completo, potrebbe non esserci una soluzione del sistema lineare $Ax=b$ . Tuttavia, si può cercare un vettore x per il quale $\|Ax-b\|_{2}$ è minimo. La soluzione è $x=A^{+}b$ . Nel caso di una matrice simmetrica S come prima, si ha $x=S^{+}b=E^{T}{\mbox{Diag}}(e^{+})Eb$ .

Matrice esponenziale

Da $S=E^{T}{\mbox{Diag}}(e)E$ si trova $\exp S=E^{T}{\mbox{Diag}}(\exp e)E$ dove exp $e$ è il vettore dove $e_{i}$ è sostituito da $\exp e_{i}$ . Allo stesso modo, $f(S)$ può essere calcolato in modo ovvio per qualsiasi funzione (analitica) $f$ .

Equazioni differenziali lineari

L'equazione differenziale $x'=Ax,x(0)=a$ ha la soluzione $x(t)=\exp(tA)$ . Per una matrice simmetrica $S$, ne consegue che $x(t)=E^{T}{\mbox{Diag}}(\exp te)Ea$ . Se $a=\sum _{i=1}^{n}a_{i}E_{i}$ è l'espansione di $a$ dagli autovettori di $S$, Poi $x(t)=\sum _{i=1}^{n}a_{i}\exp(te_{i})E_{i}$ .

Permettere $W^{s}$ essere lo spazio vettoriale attraversato dagli autovettori di $S$ che corrispondono ad un autovalore negativo e $W^{u}$ analogamente per gli autovalori positivi. Se $a\in W^{s}$ Poi ${\mbox{lim}}_{t\rightarrow \infty }x(t)=0$ ; cioè, il punto di equilibrio 0 è attrattivo $x(t)$ . Se $a\in W^{u}$ Poi ${\mbox{lim}}_{t\rightarrow \infty }x(t)=\infty$ ; cioè, 0 è ripugnante $x(t)$ . $W^{s}$ E $W^{u}$ sono chiamate varietà stabili e instabili per $S$ . Se $a$ ha componenti in entrambe le varietà, quindi un componente viene attratto e l'altro viene respinto. Quindi $x(t)$ approcci $W^{u}$ COME $t\to \infty$ .