Multilayer approximate nullspace methods for saddle point systems

Technische Samenvatting: Multilayer Approximate Nullspace Methoden voor Saddle Point Systemen

Probleemstelling
Het artikel behandelt de numerieke oplossing van grote, ijle lineaire systemen met een saddle point blokstructuur, weergegeven als:
$$W \begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} A & B \\ -C^T & 0 \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} f \\ g \end{bmatrix}$$
waarbij $A \in \mathbb{R}^{n \times n}$, $B, C \in \mathbb{R}^{n \times m}$ met $m \leq n$, en $A + A^T$ positief semidefiniet is. Deze systemen komen voor in uiteenlopende toepassingen, waaronder de tijdsintegratie van dissipatieve Hamiltoniaanse systemen, computationele vloeistofdynamica (bijv. Navier-Stokes vergelijkingen) en geoptimaliseerde beperkingen. De auteurs richten zich op de ontwikkeling van een zuiver algebraïsche, "black-box" iteratieve solver die niet afhankelijk is van probleem-specifieke fysieke informatie. De studie behandelt drie verschillende gevallen:

1. Symmetrisch: $B = C$ en $A = A^T \geq 0$.
2. Structureel symmetrisch: $B = C$ maar $A \neq A^T$ (met $A + A^T \geq 0$).
3. Algemeen: $B \neq C$ en $A \neq A^T$ (met $A + A^T \geq 0$).

Methodologie
De voorgestelde oplossing is een nieuwe klasse multilayer iteratieve schema's die fungeren als preconditioners binnen een externe Krylov-subruimte methode (specifiek, restarted flexible GMRES). De kernmethodologie vervangt klassieke exacte nullspace operaties door het gebruik van ijle benaderingen en iteratieve solvers.

1. Benadering van de Nullspace Constructie:
   In plaats van exacte nullspace bases te berekenen (die doorgaans dicht zijn), maken de auteurs gebruik van het Sparse Approximate Right Oblique Conjugation (SAROC) algoritme. Deze procedure berekent ijle benaderde bases ($\tilde{Z}$ voor de nullspace van $B^T$ en $\tilde{U}$ voor de nullspace van $C^T$) door middel van oblique conjugatie met pivoting toe te passen en drop-toleranties te gebruiken om de ijle structuur te behouden.

2. Geprojecteerde Nullspace Matrices:
   De methode construeert geprojecteerde matrices zoals $\tilde{N} = \tilde{Z}^T A \tilde{U}$ (algemeen geval) of $\tilde{N} = \tilde{Z}^T A \tilde{Z}$ (symmetrisch geval). Theoretische analyse (Lemma 3) toont aan dat, indien de benaderingsfout voldoende klein is, het symmetrische deel van deze geprojecteerde matrices de positieve definitie behoudt, wat oplosbaarheid garandeert.

3. Factorized Sparse Approximate Inverse (FSAI):
   Om de geprojecteerde systemen efficiënt te behandelen zonder expliciet dichte matrices te vormen, berekent de auteur een Factorized Sparse Approximate Inverse (FSAI) van het symmetrische deel van de geprojecteerde nullspace matrix. Deze stap genereert impliciet een preconditioner die het systeem transformeert naar een vorm die geschikt is voor shifted skew-symmetric of symmetric positive definite solvers.

4. Multilayer Iteratieve Structuur:
   Het schema werkt via een inner-outer iteratie:

   • Outer Loop: Een flexible GMRES (fGMRES) solver behandelt het hoofd saddle point systeem.
   • Inner Layers: De toepassing van de preconditioner omvat het iteratief oplossen van subproblemen:
     • Het vinden van particuliere oplossingen voor ondergedetermineerde systemen via LSQR.
     • Het oplossen van het gereduceerde nullspace systeem via Conjugate Gradient (CG) (symmetrisch geval) of fGMRES met een specifieke preconditioner voor shifted skew-symmetric systemen (structureel symmetrisch/algemeen geval).
     • Het oplossen van de binnenste shifted skew-symmetric systemen via de Minimal Residual (MRS) methode.
   • Optionele M-Orthogonalisatie: De auteurs voegen een optionele Modified Gram-Schmidt (MGS) stap toe om de nullspace basis $M$-orthogonal te maken (waarbij $M$ gerelateerd is aan $A$) om de stabiliteit te verbeteren, met name voor slecht geconditioneerde problemen.

Belangrijkste Bijdragen

• Generalisering van Nullspace Methoden: Het artikel breidt de nullspace benadering uit naar het algemene geval waar $B \neq C$ en $A$ niet-symmetrisch is, een scenario dat vaak wordt genegeerd vanwege de computationele kosten.
• Strategie voor IJle Benadering: Door SAROC en FSAI te introduceren, vermijdt de methode de dichtheidsproblemen die geassocieerd worden met exacte nullspace bases, waardoor de aanpak schaalbaar is voor grote ijle problemen.
• Theoretische Analyse: De auteurs bieden een rigoureus theoretisch kader (Theorem 1, Corollary 1, Lemma 3) dat de voorwaarden vaststelt waaronder de benaderde geprojecteerde matrices positief definit blijven, zelfs wanneer $A$ indefinit of de nullspace bases benaderd zijn.
• Robuuste Solver Class: De ontwikkeling van een verenigde solver class die zich aanpast aan symmetrische, structureel symmetrische en algemene saddle point structuren met behulp van een consistente algebraïsche framework.

Experimentele Resultaten
De auteurs evalueerden de nieuwe methode tegenover een baseline van restarted GMRES gepreconditioneerd met ILUTP (Incomplete LU met pivoting) en, voor symmetrische gevallen, een incomplete LDL$^T$ factorisatie. Tests werden uitgevoerd op matrices uit de SuiteSparse collectie en gegenereerd vanuit IFISS (vloeistofdynamica) en willekeurige generatoren.

• Robuustheid: De nieuwe methode toonde superieure robuustheid. In de symmetrische en algemene gevallen faalde de baseline ILUTP regelmatig door nul-pivots of divergentie, terwijl de nieuwe methode convergeerde voor alle testproblemen (met uitzondering van kleine incidenten in het algemene geval onder zeer ruime toleranties).
• Geheugenefficiëntie: Wat betreft het aantal nonzeros (een proxy voor geheugengebruik), vereiste de nieuwe methode over het algemeen minder nonzeros dan de ILUTP baseline, met name in de symmetrische en algemene gevallen.
• Iteratie-aantallen: De nieuwe methode convergeerde consistent binnen de maximale iteratielimiet (1.000), terwijl de baseline vaak deze limiet bereikte zonder de residu-tolerantie te voldoen.
• Impact van Orthogonalisatie: Voor de structureel symmetrische en algemene gevallen veranderde M-orthogonalisatie de convergentiesnelheden niet significant. Echter, in het symmetrische geval bleek het nuttig voor zeer slecht geconditioneerde problemen (bijv. het Goddard rocket probleem), waarbij het het aantal externe iteraties aanzienlijk verminderde.

Betekenis en Claims
Het artikel claimt dat de voorgestelde multilayer iteratieve schema's een robuust en effectief alternatief bieden voor klassieke preconditioning methoden voor saddle point systemen, vooral wanneer probleem-specifieke informatie niet beschikbaar is. De auteurs benadrukken dat de methode erin slaagt de computationele en geheugenmatige nadelen van exacte nullspace methoden aan te pakken door gebruik te maken van ijle benaderingen.

De betekenis ligt in het vermogen om een breed scala aan saddle point problemen (van symmetrisch tot algemeen) te behandelen met één enkel, zuiver algebraïsch framework. De auteurs merken bescheiden op dat hoewel de huidige implementatie een "black-box" benadering is, de prestaties verder verbeterd kunnen worden door extra preconditioning voor de interne least squares solvers toe te voegen. Zij identificeren ook parallelle implementatie voor grootschalige problemen als een richting voor toekomstig werk. De methode is bijzonder relevant voor de impliciete discretisering van tijdvariërende of niet-lineaire dissipatieve Hamiltoniaanse systemen waar reeksen van licht variërende systemen opgelost moeten worden.