Giampaolo Bottoni - 29 agosto 2008
La ricerca degli autovalori di una matrice non simmetrica è nettamente più difficile della ricerca degli autovalori di una matrice simmetrica.
L'algoritmo più solido per matrici piene ossia con quasi tutti gli elementi non nulli, è quello di Francis [ Computer J. 4, 265-271, 332-345 (1961/62) ]. Si trova citato in tutti i libri recenti ossia dagli anni 70, di analisi numerica ( tipo quello di J.Stoer e R.Bulirsch - Introduzione all'analisi numerica vol. 2 edito da Zanichelli ).
Quando la matrice è simmetrica essa può essere preliminarmente trasformata in una matrice tridiagonale simmetrica ed ovviamente le operazioni su una matrice tridiagonale sono molto più veloci che su una matrice piena.
Quando la matrice non è simmetrica sarebbe molto bello poterla trasformare in una tridiagonale non simmetrica ma.... questa trasformazione non è sicura ossia l'algoritmo, in certi casi, non può essere portato a termine ( un po' come quando si deve risolvere un sistema di equazioni lineari e la matrice risulta singolare o pericolosamente quasi singolare ).
Naturalmente si trovano in letteratura vari algoritmi per tridiagonalizzare una matrice qualsiasi e... se fossi un matematico non dilettantesco li avrei dovuti studiare bene prima di proporre il mio metodo. Non so se questo metodo è la riscoperta dell'acqua calda. Ad un certo punto avevo creduto che riuscisse sempre ad andare in porto ma poi mi sono accorto che mi sbagliavo ( o forse no, stando alla sperimentazione numerica che ho fatto). Se fosse stato un algoritmo incondizionatamente stabile sarebbe stato sicuramente nuovo perché nessun algoritmo di tridiagonalizzazione noto in precedenza poteva vantare questa dote.
Forse la sua principale dote è quella di essere concettualmente molto semplice per cui ho pensato di descriverlo pur conoscendo i suoi limiti e il rischio di descrivere, con parole mie, un algoritmo noto ai matematici professionalmente esperti di calcolo numerico.
Come premessa va ricordato che l'algoritmo di Francis ossia il metodo QR richiede che la matrice sia stata preventivamente trasformata in forma di Hessenberg ossia una triangolare superiore con anche la prima sottodiagonale non nulla. Per trasformare la matrice in forma di Hessenberg si usano matrici unitarie di Householder. Se queste matrici fossero applicate ad una matrice simmetrica il risultato sarebbe quello di ottenere una tridiagonale simmetrica e dunque la forma di Hessenberg è causata dal fatto che la matrice non è simmetrica.
Per chiarire le cose traccio qui una matrice in forma di Hessemberg con la convenzione che un punto indica un elemento nullo e una x indica un elemento solo casualmente nullo.
| x x x x x x x | | x x x x x x x | | . x x x x x x | Matrice in forma di Hessenberg : | . . x x x x x | | . . . x x x x | | . . . . x x x | | . . . . . x x |Ovviamente sarebbe molto meglio riuscire a trasformare la matrice in forma tridiagonale ossia:| x x . . . . . | | x x x . . . . | | . x x x . . . | Matrice in forma tridiagonale : | . . x x x . . | | . . . x x x . | | . . . . x x x | | . . . . . x x |Ogni matrice può essere trasformata in una altra matrice avente gli stessi autovalori solo se essa viene premoltiplicata per una matrice non singolare e postmoltiplicata per la matrice inversa di quella usata per premoltiplicare. Se A è la matrice di cui cerchiamo gli autovalori, sia T una matrice di trasformazione che ha il solo vincolo di non essere singolare ossia deve esse tale per cui esista anche l'inversa ossia T−1. Allora A ha gli stessi autovalori di T·A·T−1. Non tutte le matrici di trasformazione sono di buona qualità ossia sono tali da non produrre errori numerici tali da rendere poco affidabile la trasformazione. Se una matrice T è quasi singolare la sua inversa esiste ma è determinabile con grossa approssimazione numerica. Analogamente T rappresenta una cattiva matrice di trasformazione se i suoi autovalori sono molto differenti, in modulo, l'uno dall'altro. Le migliori matrici di trasformazione sono quelle in cui tutti gli autovalori sono identici in valore assoluto ossia tutti hanno modulo uno ( se le matrici sono reali possono avere coppie di autovalori complessi coniugati ma aventi comunque sempre modulo unitario ). Queste matrici sono dette matrici unitarie ( in due o tre dimensioni sono matrici unitarie le rotazioni; in tre dimensioni ogni rotazione ha un autovalore reale e due complessi coniugati e l'autovettore associato all'autovalore reale è il vettore diretto come l'asse di rotazione, dato che la rotazione applicata a tale vettore non lo modifica ma lo lascia inalterato ossia come se fosse stato moltiplicato per l'autovalore 1 ).
Data la loro importanza indico con U una matrice unitaria, adattissima per effettuare trasformazioni per questa sua importante caratteristica: l'inversa di una trasformazione unitaria coincide con la sua trasposta ovvero U −1 = UT. Va notato che ogni calcolo numerico provoca degli errori di arrotondamento ma questo non avviene nel caso dell'inversione delle matrici unitarie perché tutto si riduce a trasferire semplicemente dei dati ed anche queste operazioni sono relativamente poche ossia circa il quadrato dell'ordine della matrice mentre per calcolare l'inversa di una matrice qualsiasi occorre effettuare un numero di operazioni pari a circa il cubo dell'ordine della matrice.
Senza entrare nel dettaglio di come realizzare una trasformazione di Householder importa precisare che tale tipo di trasformazione corrisponde ad una matrice simmetrica ossia una matrice di Householder coincide con la sua inversa per cui... nel caso di matrici di Householder non è neppure necessario calcolare la trasposta. Altra caratteristica numerica importante è il fatto che un prodotto matrice per vettore richiede, in genere, un numero di operazioni pari al quadrato dell'ordine della matrice; nel caso delle trasformazioni di Householder invece esiste un algoritmo che consente di fare il prodotto matrice per vettore effettuando solo un numero di operazioni pari al doppio dell'ordine della matrice. Questo risparmio è fondamentale soprattutto per matrici di ordine elevato, oltre l'ordine 10.
Le matrici di Householder vengono usate per azzerare gli elementi della matrice che si cerca di semplificare. Dato un qualsiasi vettore non identicamente nullo si può trovare una trasformazione unitaria tale per cui il vettore risultato ha tutti gli elementi nulli ad eccezione di uno solo. Quando si mira a trasformare una matrice arbitraria in matrice di Hessenberg si applicano successive trasformazioni di Householder che azzerano i vettori colonna sottodiagonali in modo che resti non nullo il solo elemento della prima sottodiagonale.
Per maggiore chiarezza visualizzo cosa succede ad una matrice qualsiasi quando viene pre e postmoltiplicata per le successive trasformazioni di Householder:| x x x x x x x | | x x x x x x x | | x x x x x x x | Matrice arbitraria iniziale : | x x x x x x x | | x x x x x x x | | x x x x x x x | | x x x x x x x |Dopo aver pre e post moltiplicato per la prima trasformazione di Householder pensata per azzerare al massimo il primo vettore colonna che inizia sotto il primo elemento diagonale ottengo:| x x x x x x x | | x x x x x x x | | . x x x x x x | Matrice al passo n. 1 : | . x x x x x x | | . x x x x x x | | . x x x x x x | | . x x x x x x |Si noti che riesco a rendere nulli quasi tutti gli elementi della prima colonna sottodiagonale ma non riesco ad azzerare i simmetrici elementi della prima riga perché la matrice non è simmetrica. Se la matrice fosse simmetrica anche gli elementi simmetrici della prima riga si annullerebbero perché la pre e post moltiplicazione per la matrice di Householder conserva la simmetria. Dopo aver effettuato il secondo passo ottengo una matrice così fatta:| x x x x x x x | | x x x x x x x | | . x x x x x x | Matrice al passo n. 2 : | . . x x x x x | | . . x x x x x | | . . x x x x x | | . . x x x x x |e continuando, al passo dopo ho questa matrice:| x x x x x x x | | x x x x x x x | | . x x x x x x | Matrice al passo n. 3 : | . . x x x x x | | . . . x x x x | | . . . x x x x | | . . . x x x x |e così via fino ad ottenere una matrice in forma di Hessenberg se la matrice di partenza non era simmetrica e una tridiagonale simmetrica se la matrice di partenza era simmetrica.
Vediamo ora di formulare l'algoritmo per la tridiagonalizzazione di una matrice qualsiasi. Supponiamo di aver già effettuato vari passi con successo e di trovarci al... terzo passo da fare. La matrice da trasformare avrà la seguente forma:| x x . . . . . | | x x x . . . . | | . x x x x x x | Matrice prima del terzo passo : | . . x x x x x | | . . x x x x x | | . . x x x x x | | . . x x x x x |Come si vede, per non interferire con il lavoro fatto nei passi precedenti occorre realizzare una trasformazione ( pre e post moltiplicazione ) con una matrice che sia diagonale unitaria nelle prime tre righe e colonne e che sia una effetiva matrice di Householder dalla quarta riga e colonna in poi. Dato che stiamo supponendo di operare su una matrice non simmetrica la trasformazione adatta per manipolare una colonna non è capace di azzerare contemporaneamente la stessa riga per cui dobbiamo effettuare, in stretta successione, non una ma DUE trasformazioni di Householder. Dopo la prima trasformazione avremo una matrice così:| x x . . . . . | | x x x . . . . | | . x x x x x x | Dopo il primo semipasso : | . . x x x x x | | . . o x x x x | | . . o x x x x | | . . o x x x x |Si sono indicati con una "o" invece che con un punto gli elementi azzerati dal primo semipasso. Il secondo semipasso invece è pensato per azzerare gli elementi della corrispondente riga che, per la mancanza di simmetria, non sono stati azzerati dal primo semipasso. Dopo il secondo semipasso la matrice avrà questa struttura:| x x . . . . . | | x x x . . . . | | . x x x x o o | Dopo il secondo semipasso : | . . x x x x x | | . . o x x x x | | . . o x x x x | | . . o x x x x |Si noti questo fatto importante. Gli elementi azzerati col secondo semipasso sono uno in meno di quelli azzerati col primo semipasso. Infatti non è stato possibile usare una trasformazione di Householter di ordine uguale a quello della trasformazione usata nel primo semipasso perché, facendo così, avremmo interferito con l'altra trasformazione ossia avremmo azzerato un elemento in più sulla riga ma reso non nulli tutti gli elementi della corrispondente colonna.
Dato che ci stiamo avvicinando al punto cruciale dell'algoritmo scriviamo più estesamente la parte interessante della matrice a partire dalla terza riga e colonna distinguendo gli elementi:| a33 a34 a35 0 0 | | a43 a44 a45 a46 a47 | La parte finale della matrice in dettaglio : | 0 a54 a55 a56 a57 | | 0 a64 a65 a66 a67 | | 0 a74 a75 a76 a77 |Se a35 fosse nullo non ci sarebbe niente altro da fare perché i due semipassi avrebbero raggiunto lo scopo di rendere tridiagonale la terza riga e colonna della matrice. Se la matrice è simmetrica l'elemento a35 sarebbe stato nullo come tutti gli altri dopo a34 ossia non ci sarebbe stata la necessità di effettuare anche il secondo semipasso. Se la matrice fosse stata simmetrica sarebbero stati uguali tra loro gli elementi a34 e a43 perché tutte le trasformazioni fatte ( basate su premoltiplicazione e postmoltiplicazione con una matrice simmetrica di Householder) conservano la simmetria.
Per poter concludere il passo è indispensabile fare in modo che l'elemento a35 risulti nullo e fare anche in modo che la colonna sotto l'elemento a43 continui ad essere fatta da elementi tutti nulli. Ma questo SI PUÒ FARE ! usando questa matrice di trasformazione:| 1 0 0 0 0 | | 0 1 t 0 0 | Matrice di fine passo ( premoltiplicatrice) : | 0 0 1 0 0 | | 0 0 0 1 0 | | 0 0 0 0 1 |Per effettuare una trasformazione bisogna conoscere la matrice inversa della matrice premoltiplicatrice. Ma l'inversa della precedente matrice è banale, come è immediato verificare:| 1 0 0 0 0 | | 0 1 -t 0 0 | Matrice di fine passo ( postmoltiplicatrice) : | 0 0 1 0 0 | | 0 0 0 1 0 | | 0 0 0 0 1 |Il valore che l'elemento t deve avere è questo:
t = a35 / a34
Tutto sembrerebbe perfetto ossia, a valle di questa veloce trasformazione che interessa solo la riga 4 e la colonna 5, avremmo ottenuto il risultato di tridiagonalizzare la terza riga e colonna della nostra matrice arbitraria. Ma c'è un importante MA: il valore numerico dell'elemento t diventa ( ovvero tende all' ) infinito se l'elemento a34 vale esattamente ( ovvero tenda a ) zero.
Questa situazione non può mai verificarsi con matrici simmetriche perché se l'elemento a34 fosse nullo lo sarebbe, per simmetria, anche l'elemento a43 ma questo vorrebbe dire che la matrice è partizionabile in due blocchi distinti e tra loro indipendenti. Il fatto che t resti piccolo è una proprietà delle matrici poco asimmetriche che, a puro scopo orientativo, saranno definite come quelle matrici in cui, in tutti i passi, il valore dei vari t non supera in valore assoluto l'unità.
Dal punto di vista matematico si potrebbe aprire a questo punto un ampio settore di ricerca su come caratterizzare le matrici non simmetriche per capire se sono poco asimmetriche ovvero con tridiagonalizzazione fattibile e poco perturbabile variando leggermente il valore degli elementi della matrice di partenza.Risultati della sperimentazione numerica
Dalla sperimentazione numerica è emerso un fatto sorprendente. Una qualsiasi matrice di trasformazione T tale da non modificare il primo termine della diagonale non influisce sul fatto che il valore di t vada all'infinito. In altre parole, una trasformazione in cui la prima riga e colonna sono nulle, tranne ovviamente il termine diagonale che deve valere 1, non riesce ad evitare che t raggiunga valori altissimi. Dall'altra parte una trasformazione che coinvolge il primo termine della diagonale modifica i valori di t e dunque permette di evitare valori di t molto elevati in valore assoluto. Se dunque l'obbiettivo fosse quello di evitare che t risulti troppo grande basterebbe riapplicare l'algoritmo dopo aver applicato alla matrice originaria qualche matrice di rotazione ... a caso. La probabilità che la nuova matrice presenti gli stessi problemi di valori di t troppo alti sarebbe piuttosto bassa.
Per effettuare la sperimentazione numerica visitare questa pagina: xt2(usafiledati).xhtml.
...per modificare la matrice di test...
Matrice in costruzione:Test dell'algoritmo per la tridiagonalizzazione di una matrice
Questa è la matrice di partenza: Parte reale autovalori della matrice di test: Parte immaginaria autovalori della matrice di test: Info sul calcolo della tridiagonalizzazione: La matrice test tridiagonalizzata:La matrice tridiagonalizzata deve avere gli stessi autovalori di quella di partenza
Parte reale autovalori della matrice di test tridiagonalizzata: Parte immaginaria autovalori della matrice di test tridiagonalizzata: