Multilayer approximate nullspace methods for saddle point systems

Technical Summary: Multilayer Approximate Nullspace Methods for Saddle Point Systems

Problem Statement
The paper addresses the numerical solution of large, sparse linear systems with a saddle point block structure, represented as:
$$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}$$
where $A \in \mathbb{R}^{n \times n}$, $B, C \in \mathbb{R}^{n \times m}$ with $m \leq n$, and $A + A^T$ is positive semidefinite. These systems arise in diverse applications, including the time integration of dissipative Hamiltonian systems, computational fluid dynamics (e.g., Navier-Stokes equations), and constrained optimization. The authors focus on developing a purely algebraic, "black-box" iterative solver that does not rely on problem-specific physical information. The study covers three distinct cases:

1. Symmetric: $B = C$ and $A = A^T \geq 0$.
2. Structurally Symmetric: $B = C$ but $A \neq A^T$ (with $A + A^T \geq 0$).
3. General: $B \neq C$ and $A \neq A^T$ (with $A + A^T \geq 0$).

Methodology
The proposed solution is a new class of multilayer iterative schemes that function as preconditioners within an outer Krylov subspace method (specifically, restarted flexible GMRES). The core methodology modifies the classical exact nullspace method by replacing exact operations with sparse approximations and iterative solvers.

1. Approximate Nullspace Construction:
   Instead of computing exact nullspace bases (which are typically dense), the authors employ the Sparse Approximate Right Oblique Conjugation (SAROC) algorithm. This procedure computes sparse approximate bases ($\tilde{Z}$ for the nullspace of $B^T$ and $\tilde{U}$ for the nullspace of $C^T$) by performing oblique conjugation with pivoting and applying drop tolerances to maintain sparsity.

2. Projected Nullspace Matrices:
   The method constructs projected matrices such as $\tilde{N} = \tilde{Z}^T A \tilde{U}$ (general case) or $\tilde{N} = \tilde{Z}^T A \tilde{Z}$ (symmetric case). Theoretical analysis (Lemma 3) demonstrates that if the approximation error is sufficiently small, the symmetric part of these projected matrices retains positive definiteness, ensuring solvability.

3. Factorized Sparse Approximate Inverse (FSAI):
   To handle the projected systems efficiently without explicitly forming dense matrices, the authors compute a Factorized Sparse Approximate Inverse (FSAI) of the symmetric part of the projected nullspace matrix. This step implicitly generates a preconditioner that transforms the system into a form suitable for shifted skew-symmetric or symmetric positive definite solvers.

4. Multilayer Iterative Structure:
   The scheme operates as an inner-outer iteration:

   • Outer Loop: A flexible GMRES (fGMRES) solver handles the main saddle point system.
   • Inner Layers: The preconditioner application involves solving sub-problems iteratively:
     • Finding particular solutions for underdetermined systems via LSQR.
     • Solving the reduced nullspace system via Conjugate Gradient (CG) (symmetric case) or fGMRES with a specific preconditioner for shifted skew-symmetric systems (structurally symmetric/general cases).
     • Solving the innermost shifted skew-symmetric systems via the Minimal Residual (MRS) method.
   • Optional M-Orthogonalization: The authors include an optional Modified Gram-Schmidt (MGS) step to $M$-orthogonalize the nullspace basis (where $M$ is related to $A$) to improve stability, particularly for ill-conditioned problems.

Key Contributions

• Generalization of Nullspace Methods: The paper extends the nullspace approach to the general case where $B \neq C$ and $A$ is non-symmetric, a scenario often neglected due to computational costs.
• Sparse Approximation Strategy: By introducing SAROC and FSAI, the method avoids the density issues associated with exact nullspace bases, making the approach scalable for large sparse problems.
• Theoretical Analysis: The authors provide a rigorous theoretical framework (Theorem 1, Corollary 1, Lemma 3) establishing conditions under which the approximate projected matrices remain positive definite, even when $A$ is indefinite or the nullspace bases are approximate.
• Robust Solver Class: The development of a unified solver class that adapts to symmetric, structurally symmetric, and general saddle point structures using a consistent algebraic framework.

Experimental Results
The authors evaluated the new scheme against a baseline of restarted GMRES preconditioned with ILUTP (Incomplete LU with pivoting) and, for symmetric cases, an incomplete LDL$^T$ factorization. Tests were conducted on matrices from the SuiteSparse collection and generated from IFISS (fluid dynamics) and random generators.

• Robustness: The new scheme demonstrated superior robustness. In the symmetric and general cases, the baseline ILUTP frequently failed due to zero pivots or divergence, whereas the new scheme converged for all test problems (with minor exceptions in the general case under very loose tolerances).
• Memory Efficiency: In terms of the number of nonzeros (a proxy for memory usage), the new scheme generally required fewer nonzeros than the ILUTP baseline, particularly in the symmetric and general cases.
• Iteration Counts: The new scheme consistently converged within the maximum iteration limit (1,000), while the baseline often hit this limit without satisfying the residual tolerance.
• Impact of Orthogonalization: For the structurally symmetric and general cases, M-orthogonalization did not significantly alter convergence rates. However, in the symmetric case, it proved beneficial for highly ill-conditioned problems (e.g., the Goddard rocket problem), reducing outer iterations significantly.

Significance and Claims
The paper claims that the proposed multilayer iterative schemes offer a robust and effective alternative to classical preconditioning methods for saddle point systems, particularly when problem-specific information is unavailable. The authors emphasize that the method successfully addresses the computational and memory disadvantages of exact nullspace methods by utilizing sparse approximations.

The significance lies in the ability to treat a wide range of saddle point problems (symmetric to general) with a single, purely algebraic framework. The authors modestly note that while the current implementation is a "black-box" approach, performance could be further enhanced by incorporating additional preconditioning for the internal least squares solvers. They also identify parallel implementation for large-scale problems as a direction for future work. The method is presented as particularly relevant for implicit discretization of time-varying or nonlinear dissipative Hamiltonian systems where sequences of slightly varying systems must be solved.