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

Technical Summary: The Shift-and-Invert Arnoldi Method for Singular Matrix Pencils

Problem Statement
The paper addresses the challenge of solving large-scale sparse singular generalized eigenvalue problems of the form $Ax = \lambda Bx$, where the matrix pencil $A - \lambda B$ is singular (i.e., $\det(A - \lambda B) \equiv 0$). In such cases, the classical definition of eigenvalues fails because any scalar satisfies the determinant condition. The meaningful solutions are the "regular" eigenvalues, defined where the rank of $A - \lambda_0 B$ drops below the normal rank of the pencil.

While methods like the staircase method or QZ decomposition exist for small, dense problems, they are computationally prohibitive for large-scale sparse systems. Recent approaches by Hochstenbach, Mehl, and Plestenjak propose regularizing singular pencils using randomly chosen dense matrices to separate true eigenvalues from spurious ones. However, the authors argue that using dense regularization matrices destroys the sparsity structure of the original problem, leading to significant performance and memory overheads for large-scale applications. Furthermore, existing sparse direct solvers do not natively support the row/column removal or addition required to regularize these pencils without custom development.

Methodology
The authors propose a deterministic strategy to regularize singular pencils by leveraging the pivoting sequence of a sparse LU factorization. The core methodology involves the following steps:

1. Regularization via LU Pivoting: Instead of random dense matrices, the method uses a sparse LU factorization of the shifted matrix $A - \sigma B$ (where $\sigma$ is a shift). During the factorization, if a zero (or near-zero) pivot is encountered, the algorithm adds a row and column to the matrix to maintain the factorization.

   • Augmented Pencil (AGEP): The process constructs matrices $F$ and $G$ such that the augmented pencil $\begin{bmatrix} A & G \\ F^* & 0 \end{bmatrix} - \lambda \begin{bmatrix} B & 0 \\ 0 & 0 \end{bmatrix}$ becomes regular. The matrices $F$ and $G$ are derived directly from the permutation matrices and the pivot selection process of the LU factorization.
   • Projected Pencil (PGEP): Alternatively, the method constructs a projected pencil $G_\perp^* (A - \lambda B) F_\perp$ by removing the rows and columns corresponding to the zero pivots. $F_\perp$ and $G_\perp$ are derived as the identity matrix with specific rows removed.

2. Rank Detection and Correction: The LU factorization naturally attempts to reveal the normal rank of the pencil. The authors note that due to rounding errors or inappropriate tolerance settings, the detected rank ($\hat{k}$) might differ from the true normal rank ($k$).

   • If the detected rank is too low (i.e., the tolerance is too aggressive), a rank correction scheme is proposed. This involves using random orthogonal matrices to augment the detected $F$ and $G$ to the correct dimension, followed by a low-cost Schur complement update to solve the resulting system.

3. Shift-and-Invert Arnoldi: The regularized problem (either AGEP or PGEP) is solved using the two-sided shift-and-invert Arnoldi method. The authors prove an equivalence theorem showing that the Arnoldi factorization on the augmented pencil (with a specific starting vector and inner product) yields the same Ritz values and Hessenberg matrix as the projected pencil. This allows the use of the projected formulation to avoid the explicit construction of the large augmented system while maintaining the same spectral properties.

4. Rectangular Problems: The method is extended to rectangular eigenvalue problems ($A, B \in \mathbb{R}^{n \times m}$ with $n \neq m$). For full-rank rectangular problems, existing sparse solvers with row pivoting can be used directly to select a square, regular sub-problem.

Key Contributions

• Sparse Regularization: The paper introduces a method to generate sparse regularization matrices ($F, G, F_\perp, G_\perp$) directly from the pivoting sequence of sparse LU factorizations, preserving the sparsity of the original problem.
• Rank Detection: It demonstrates that standard sparse LU factorizations (using partial, rook, or complete pivoting) often successfully detect the normal rank of the pencil, facilitating the automatic construction of the regularization.
• Rank Correction: A specific algorithm is proposed to correct the regularization if the LU factorization underestimates the normal rank, ensuring the method remains robust even with imperfect rank detection.
• Theoretical Equivalence: The paper establishes a rigorous theoretical link between the augmented and projected regularized pencils within the context of the Arnoldi method, proving that they yield identical spectral approximations under specific conditions.

Results
Numerical experiments were conducted on three types of problems:

1. Finite Element Model Updating: A truss structure model with 33 degrees of freedom (resulting in a pencil of dimension ~35,937). The proposed method (using RLU pivoting) correctly identified the rank and computed eigenvalues with residuals on the order of $10^{-14}$. In contrast, the randomized approach resulted in condition numbers orders of magnitude higher ($10^6$ vs $10^3$) and larger residuals.
2. Double Eigenvalue Detection: A parametric problem involving the discretization of a 2D Laplacian. The LU-based method correctly determined the rank and achieved high accuracy. The randomized approach required constructing dense matrices, leading to significantly higher storage and computational costs.
3. Rectangular Nonlinear Eigenvalue Problem: A full-rank rectangular problem derived from a nonlinear eigenvalue problem. The sparse LU approach was significantly faster (0.05s vs 2.9s for factorization and projection) and required less memory (98MB vs 204MB) compared to the randomized dense approach.

In all cases, the sparse regularization methods preserved the sparsity of the matrices, leading to lower condition numbers for the shifted systems and faster convergence compared to randomized dense regularization.

Significance and Claims
The authors position this work as a concept rather than a complete, optimized software solution, acknowledging that robust sparse direct solvers for singular matrices are currently missing from standard software libraries. However, the paper claims that:

• Performance: Preserving sparsity through deterministic pivoting is beneficial for both the performance of the eigenvalue solver and the accuracy of the results (via better conditioning).
• Feasibility: Existing sparse direct solvers can be adapted to solve these problems, particularly for rectangular full-rank cases, without needing new, complex software infrastructure.
• Robustness: The proposed rank correction scheme ensures the method works even when the LU factorization fails to perfectly detect the normal rank.

The paper concludes that while the randomized approach is theoretically sound, the deterministic sparse LU approach offers a practical, efficient, and accurate alternative for large-scale singular eigenvalue problems, provided that the gap between singular values is sufficient for rank detection. Future work is suggested to explore solvers with global maximum volume pivoting and to exploit Kronecker structures in multiparameter problems.