Algoritmo degli autovalori di Jacobi

Questa voce è orfana, ovvero priva di collegamenti in entrata da altre voci.

Inseriscine almeno uno pertinente e utile e rimuovi l'avviso. Segui i suggerimenti del progetto di riferimento.

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,^[1] ma divenne ampiamente utilizzato solo negli anni '50 con l'avvento dei computer.^[2]

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

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

è simmetrica e simile a ${\displaystyle S}$ .

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

${\displaystyle {\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 ${\displaystyle s=\sin(\theta )}$ e ${\displaystyle c=\cos(\theta )}$.

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

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

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

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

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

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

Per ottimizzare questo effetto, ${\displaystyle 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.

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

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

Di conseguenza ${\displaystyle \Gamma (S^{J})^{2}=\Gamma (S)^{2}-2p^{2}}$ . Questo implica ${\displaystyle \Gamma (S^{J})^{2}\leq (1-1/N)\Gamma (S)^{2}}$ oppure ${\displaystyle \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 ${\displaystyle (1-1/N)^{1/2}}$ ad una matrice diagonale.

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

${\displaystyle \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 ≈ ${\displaystyle e^{1/2}}$ .

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

${\displaystyle 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 ${\displaystyle S^{s}}$ denota il risultato, allora

${\displaystyle \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 ${\displaystyle \Gamma (S)<{\frac {d}{2+{\sqrt {{\frac {n}{2}}-1}}}}}$

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 ${\displaystyle m_{1},\,\dots \,,\,m_{n-1}}$, con la proprietà che ${\displaystyle 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 ${\displaystyle (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 ${\displaystyle m_{i}}$ non è né k né l, è sufficiente confrontare il vecchio massimo in ${\displaystyle m_{i}}$ con i nuovi elementi, aggiornando ${\displaystyle m_{i}}$ se necessario. Se ${\displaystyle 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 ${\displaystyle 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 ${\displaystyle N_{S}<N}$ .

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 ${\displaystyle e_{i}}$ è un autovalore e ogni colonna ${\displaystyle E_{i}}$ un autovettore ortonormale per ${\displaystyle e_{i}}$, io = 1, ..., n .

1. L'array logico changed mantiene lo stato di ciascun autovalore. Se il valore numerico di ${\displaystyle e_{k}}$ o ${\displaystyle 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 ${\displaystyle 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.

Sia ${\displaystyle 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) :

${\displaystyle e_{1}=2585.25381092892231}$

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

${\displaystyle e_{2}=37.1014913651276582}$

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

${\displaystyle e_{3}=1.4780548447781369}$

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

${\displaystyle e_{4}=0.1666428611718905}$

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

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). ${\displaystyle A}$ sono le radici quadrate degli autovalori (non negativi) di ${\displaystyle A^{T}A}$ . Nel caso di una matrice simmetrica ${\displaystyle S}$ si ha che ${\displaystyle S^{T}S=S^{2}}$, quindi i valori singolari di ${\displaystyle S}$ sono i valori assoluti degli autovalori di ${\displaystyle 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 ${\displaystyle \|Ax\|_{2}}$ quando x attraversa tutti i vettori con ${\displaystyle \|x\|_{2}=1}$ . È il valore singolare più grande di ${\displaystyle 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 ${\displaystyle A}$ è definito come ${\displaystyle {\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 ${\displaystyle A}$ ha rango ${\displaystyle r}$ se ha ${\displaystyle r}$ colonne linearmente indipendenti, mentre le restanti colonne sono linearmente dipendenti da queste. Equivalentemente, ${\displaystyle r}$ è la dimensione dell'intervallo di ${\displaystyle 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 ${\displaystyle A}$ è la matrice unica ${\displaystyle X=A^{+}}$ per cui ${\displaystyle AX}$ E ${\displaystyle XA}$ sono simmetriche e per cui la condizione${\displaystyle AXA=A,XAX=X}$ è verificata.

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

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

Soluzione dei minimi quadrati

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

Matrice esponenziale

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

Equazioni differenziali lineari

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

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