The shift-and-invert Arnoldi method for singular matrix pencils

Sintesi Tecnica: Il Metodo di Arnoldi Shift-and-Invert per Fasci di Matrici Singolari

Enunciato del Problema
Il lavoro affronta la sfida di risolvere problemi agli autovalori generalizzati singolari, sparsi e su larga scala della forma $Ax = \lambda Bx$, dove il fascio di matrici $A - \lambda B$ è singolare (cioè, $\det(A - \lambda B) \equiv 0$). In tali casi, la definizione classica di autovalori fallisce perché qualsiasi scalare soddisfa la condizione del determinante. Le soluzioni significative sono gli autovalori "regolari", definiti nei punti in cui il rango di $A - \lambda_0 B$ scende al di sotto del rango normale del fascio.

Sebbene esistano metodi come il metodo a gradini (staircase method) o la decomposizione QZ per problemi piccoli e densi, questi sono computazionalmente proibitivi per sistemi sparsi su larga scala. Approcci recenti proposti da Hochstenbach, Mehl e Plestenjak suggeriscono di regolarizzare fasci singolari utilizzando matrici dense scelte casualmente per separare gli autovalori reali da quelli spurii. Tuttavia, gli autori sostengono che l'uso di matrici di regolarizzazione dense distrugge la struttura sparsa del problema originale, portando a significativi sovraccarichi di prestazioni e memoria per applicazioni su larga scala. Inoltre, i solutori diretti sparsi esistenti non supportano nativamente la rimozione o l'aggiunta di righe/colonne necessaria per regolarizzare questi fasci senza uno sviluppo personalizzato.

Metodologia
Gli autori propongono una strategia deterministica per regolarizzare fasci singolari sfruttando la sequenza di pivot di una fattorizzazione LU sparsa. La metodologia di base comprende i seguenti passaggi:

1. Regolarizzazione tramite Pivot LU: Invece di matrici dense casuali, il metodo utilizza una fattorizzazione LU sparsa della matrice spostata $A - \sigma B$ (dove $\sigma$ è uno spostamento). Durante la fattorizzazione, se viene incontrato un pivot nullo (o vicino allo zero), l'algoritmo aggiunge una riga e una colonna alla matrice per mantenere la fattorizzazione.

   • Fascio Augmentato (AGEP): Il processo costruisce le matrici $F$ e $G$ in modo che il fascio augmentato $\begin{bmatrix} A & G \\ F^* & 0 \end{bmatrix} - \lambda \begin{bmatrix} B & 0 \\ 0 & 0 \end{bmatrix}$ diventi regolare. Le matrici $F$ e $G$ sono derivate direttamente dalle matrici di permutazione e dal processo di selezione dei pivot della fattorizzazione LU.
   • Fascio Proiettato (PGEP): In alternativa, il metodo costruisce un fascio proiettato $G_\perp^* (A - \lambda B) F_\perp$ rimuovendo le righe e le colonne corrispondenti ai pivot nulli. $F_\perp$ e $G_\perp$ sono derivati come la matrice identità con specifiche righe rimosse.

2. Rilevamento e Correzione del Rango: La fattorizzazione LU tenta naturalmente di rivelare il rango normale del fascio. Gli autori notano che, a causa di errori di arrotondamento o impostazioni di tolleranza inappropriate, il rango rilevato ($\hat{k}$) potrebbe differire dal vero rango normale ($k$).

   • Se il rango rilevato è troppo basso (cioè, la tolleranza è troppo aggressiva), viene proposto uno schema di correzione del rango. Questo prevede l'uso di matrici ortogonali casuali per aumentare le $F$ e $G$ rilevate alla dimensione corretta, seguito da un aggiornamento a basso costo del complemento di Schur per risolvere il sistema risultante.

3. Arnoldi Shift-and-Invert: Il problema regolarizzato (sia AGEP che PGEP) viene risolto utilizzando il metodo di Arnoldi shift-and-invert a due lati. Gli autori dimostrano un teorema di equivalenza che mostra come la fattorizzazione di Arnoldi sul fascio augmentato (con un vettore di partenza specifico e un prodotto scalare interno) produca gli stessi valori di Ritz e la stessa matrice di Hessenberg del fascio proiettato. Ciò consente l'uso della formulazione proiettata per evitare la costruzione esplicita del grande sistema augmentato, mantenendo le stesse proprietà spettrali.

4. Problemi Rettangolari: Il metodo è esteso ai problemi agli autovalori rettangolari ($A, B \in \mathbb{R}^{n \times m}$ con $n \neq m$). Per problemi rettangolari a rango pieno, i solutori sparsi esistenti con pivotazione sulle righe possono essere utilizzati direttamente per selezionare un sottoproblema quadrato e regolare.

Contributi Chiave

• Regolarizzazione Sparsa: Il lavoro introduce un metodo per generare matrici di regolarizzazione sparse ($F, G, F_\perp, G_\perp$) direttamente dalla sequenza di pivot delle fattorizzazioni LU sparse, preservando la sparsità del problema originale.
• Rilevamento del Rango: Dimostra che le fattorizzazioni LU sparse standard (utilizzando pivotazione parziale, a scacchiera o completa) spesso rilevano con successo il rango normale del fascio, facilitando la costruzione automatica della regolarizzazione.
• Correzione del Rango: Viene proposto un algoritmo specifico per correggere la regolarizzazione se la fattorizzazione LU sottostima il rango normale, garantendo che il metodo rimanga robusto anche con un rilevamento del rango imperfetto.
• Equivalenza Teorica: Il lavoro stabilisce un legame teorico rigoroso tra i fasci regolarizzati augmentati e proiettati nel contesto del metodo di Arnoldi, dimostrando che producono approssimazioni spettrali identiche in condizioni specifiche.

Risultati
Sono stati condotti esperimenti numerici su tre tipi di problemi:

1. Aggiornamento del Modello ad Elementi Finiti: Un modello di struttura reticolare con 33 gradi di libertà (risultante in un fascio di dimensione ~35.937). Il metodo proposto (utilizzando pivotazione RLU) ha identificato correttamente il rango e calcolato gli autovalori con residui dell'ordine di $10^{-14}$. Al contrario, l'approccio randomizzato ha portato a numeri di condizione ordini di grandezza superiori ($10^6$ contro $10^3$) e residui più grandi.
2. Rilevamento di Autovalori Doppi: Un problema parametrico che coinvolge la discretizzazione di un Laplaciano 2D. Il metodo basato su LU ha determinato correttamente il rango e raggiunto un'alta accuratezza. L'approccio randomizzato ha richiesto la costruzione di matrici dense, portando a costi di archiviazione e computazionali significativamente più elevati.
3. Problema Non Lineare agli Autovalori Rettangolare: Un problema rettangolare a rango pieno derivato da un problema agli autovalori non lineare. L'approccio LU sparsa è stato significativamente più veloce (0,05s contro 2,9s per fattorizzazione e proiezione) e ha richiesto meno memoria (98MB contro 204MB) rispetto all'approccio randomizzato denso.

In tutti i casi, i metodi di regolarizzazione sparsa hanno preservato la sparsità delle matrici, portando a numeri di condizione inferiori per i sistemi spostati e a una convergenza più rapida rispetto alla regolarizzazione densa randomizzata.

Significato e Affermazioni
Gli autori posizionano questo lavoro come un concetto piuttosto che come una soluzione software completa e ottimizzata, riconoscendo che i solutori diretti sparsi robusti per matrici singolari mancano attualmente nelle librerie software standard. Tuttavia, il lavoro afferma che:

• Prestazioni: Preservare la sparsità attraverso la pivotazione deterministica è vantaggioso sia per le prestazioni del risolutore agli autovalori che per l'accuratezza dei risultati (tramite un migliore condizionamento).
• Fattibilità: I solutori diretti sparsi esistenti possono essere adattati per risolvere questi problemi, in particolare per i casi rettangolari a rango pieno, senza la necessità di nuove infrastrutture software complesse.
• Robustezza: Lo schema di correzione del rango proposto garantisce che il metodo funzioni anche quando la fattorizzazione LU non rileva perfettamente il rango normale.

Il lavoro conclude che, sebbene l'approccio randomizzato sia teoricamente valido, l'approccio deterministico basato su LU sparsa offre un'alternativa pratica, efficiente e accurata per problemi agli autovalori singolari su larga scala, a condizione che il divario tra i valori singolari sia sufficiente per il rilevamento del rango. Si suggerisce un lavoro futuro per esplorare solutori con pivotazione globale a volume massimo e per sfruttare le strutture di Kronecker nei problemi a più parametri.