A fast solver for ill-conditioned linear systems using randomized stable solutions of its blocks

Problem Statement
The paper addresses the challenge of solving consistent linear systems of equations $Ax = b$, where $A \in \mathbb{R}^{m \times n}$ and $b \in \mathbb{R}^{m \times 1}$. The primary focus is on systems that are highly ill-conditioned, including sparse systems, dense least-squares problems, and cases that are significantly over- or under-determined. The authors highlight that standard iterative methods, such as GMRES, often struggle with these systems because minimizing the residual norm does not necessarily guarantee convergence in the relative error for ill-conditioned problems. Furthermore, traditional preconditioning techniques (e.g., ILU, Jacobi) can suffer from numerical instability, pivot breakdowns, or excessive fill-in when applied to matrices with very large condition numbers, often failing precisely when robust acceleration is most needed.

Methodology
The authors propose an enhanced algorithm termed Regularized Orthogonality and Residual based Block-Kaczmarz (ROR-BK). This method builds upon the randomized block-Kaczmarz framework but introduces three specific enhancements to improve stability and convergence:

1. Dynamic Residual-Based Aggregation: Instead of relying solely on fixed blocks, the method dynamically aggregates a subset of equations at each iteration based on the current residual vector. Specifically, it selects rows contributing the most to the squared residual to form a block, aiming to minimize the residual norm more effectively in each step.
2. Effective Mutual Orthogonality Sampling: The method generalizes the concept of orthogonality between block pairs (used in the Simple Orthogonal Block-Kaczmarz or SOBK method) to the "effective orthogonality" of a single block against all other blocks. A symmetric matrix $C$ is computed where entries represent the cosine of the effective angle between block centroids. Blocks are then sampled probabilistically based on a distribution that favors blocks with high cumulative orthogonality to the rest of the system.
3. Regularization: To ensure numerical stability, particularly when blocks are rank-deficient or ill-conditioned, the block updates utilize a regularized pseudo-inverse. Instead of the standard Moore-Penrose inverse $B^\dagger = B^T(BB^T)^{-1}$, the method uses the regularized form $B^T(BB^T + \mu I)^{-1}$, where $\mu$ is a small positive constant. This is interpreted as a stochastic proximal point algorithm.

The algorithm operates by alternating between updates using these orthogonality-sampled blocks and updates using the residual-based dynamic blocks. The authors also demonstrate that this method can serve as an inner iteration (preconditioner) for flexible GMRES (FAB-GMRES) solvers.

Key Contributions

• Algorithmic Innovation: The introduction of ROR-BK, which combines residual-based dynamic block selection, global orthogonality-aware sampling, and Tikhonov-style regularization within a block-Kaczmarz framework.
• Theoretical Analysis: The paper provides a convergence analysis (Theorems 2.3 and 2.4) showing that the method achieves stable convergence even when blocks are numerically rank-deficient. It establishes that the error reduction is bounded by the regularization parameter $\mu$ and the smallest non-zero eigenvalue of the block matrices. Theorem 2.4 further argues that selecting blocks based on maximum residuals improves the upper and lower bounds of the error.
• Preconditioning Strategy: The proposal of using ROR-BK as a pre-solver or inner iteration for GMRES to accelerate the solution of ill-conditioned systems without the need for explicit, potentially unstable, external preconditioners.

Results
The authors conducted extensive numerical experiments comparing ROR-BK against:

• Randomized Solvers: Simple Orthogonal Block-Kaczmarz (SOBK) and Tail Averaged Regularized Block Kaczmarz (TA-ReBlocK-U).
• Krylov Subspace Solvers: LSQR and various GMRES variants (AB/BA-GMRES) with preconditioners such as Jacobi, ILU, and Robust Incomplete Factorization (RIF).

Performance Findings:

• Speed and Iterations: Across a wide range of test matrices (including sparse matrices from the University of Florida collection and randomly generated dense matrices), ROR-BK consistently required fewer iterations and less CPU time than SOBK and TA-ReBlocK-U. Speed-ups ranged from 1.23x to 5.17x over SOBK and up to 49.08x over TA-ReBlocK-U in terms of CPU time.
• Ill-Conditioned Systems: In highly ill-conditioned scenarios (condition numbers up to $10^{20}$), many Krylov solvers and their preconditioned versions failed to converge within the iteration limit or produced large relative errors. ROR-BK successfully converged in these cases, often achieving lower relative errors.
• Real-World Application: The method was tested on image reconstruction problems (parallel-beam tomography and MRI) involving noisy data and severely underdetermined systems. ROR-BK demonstrated superior performance in terms of Peak Signal-to-Noise Ratio (PSNR) and Structural Similarity Index (SSIM) compared to other methods, particularly in scenarios where other solvers suffered from instability or failure.

Significance and Claims
The paper claims that ROR-BK offers a promising alternative for solving linear systems where the matrix $A$ is ill-conditioned, achieving efficient approximations without the need for explicit preconditioning. The authors argue that while preconditioners aim to reduce the condition number of the matrix, they do not necessarily reduce the condition number of the problem itself and can introduce numerical instability. By partitioning the system into blocks and utilizing effective orthogonality and regularization, ROR-BK distributes the covariant equations into potentially well-conditioned sub-problems.

The authors position the method as a robust solver for square, overdetermined, and underdetermined systems. They further suggest its utility as a pre-solver or an inner iteration for accelerating other iterative methods like GMRES, particularly in cases where standard preconditioning techniques fail due to pivot breakdowns or rank deficiency. The work emphasizes that the combination of residual-based dynamics and orthogonality sampling provides a simple yet effective mechanism for fast and stable convergence in high-condition number linear systems.