Accelerating a restarted Krylov method for matrix functions with randomization

1. Problem Statement

Many scientific applications (e.g., solving Partial Differential Equations, analyzing complex networks, lattice quantum chromodynamics) require computing the action of a matrix function $f(A)$ on a vector $b$, denoted as $f(A)b$.

• The Challenge: For large-scale matrices, explicitly forming $f(A)$ is computationally infeasible. The standard approach uses Krylov subspace methods (specifically the Arnoldi iteration) to approximate $f(A)b$.
• The Bottleneck: The Arnoldi process requires storing and orthogonalizing a growing basis of vectors. As the basis size $m$ increases:
  1. Memory usage becomes prohibitive ($O(nm)$).
  2. Orthogonalization cost becomes overwhelming, growing quadratically ($O(m^2n)$) for non-Hermitian matrices.
• Current Mitigation: Restarted Krylov methods are used to limit the basis size to a fixed $m$, discarding old vectors and restarting with the current approximation. However, restarting often leads to slower convergence, stagnation, or divergence because the method loses information about the global spectrum of $A$.
• The Gap: While randomized algorithms (random sketching) have been proposed to accelerate standard Krylov methods, their application to restarted methods, particularly for ill-conditioned large-scale problems, remains under-explored.

2. Methodology

The authors propose a Randomized Restarted Krylov Method that replaces the standard orthogonalization step with a Randomized Gram-Schmidt (RGS) process using oblivious subspace embeddings.

Core Components:

1. Random Sketching: Instead of full orthogonalization, the method uses a random matrix $S$ (specifically a sparse sign matrix) to project the Krylov basis into a lower-dimensional space.

   • The sketch $S$ satisfies the Johnson-Lindenstrauss-like property: $(1-\epsilon)\|x\|^2 \leq \|Sx\|^2 \leq (1+\epsilon)\|x\|^2$.
   • Orthogonalization is performed on the sketched vectors ($SW_m$) rather than the full vectors ($W_m$). This ensures the resulting basis $W_m$ is well-conditioned without the $O(m^2n)$ cost of full orthogonalization.

2. Randomized Arnoldi Iteration (Algorithm 2):

   • Constructs a non-orthogonal but well-conditioned basis $W_m$.
   • Computes the projection matrix $R_m$ (analogous to the Hessenberg matrix $H_m$ in standard Arnoldi) using the sketched inner products.
   • Complexity: Reduces the cost of orthogonalization from $O(m^2n)$ to $O(m^2 + \zeta nm)$, where $\zeta$ is the sparsity of the sketch.

3. Restarting Procedure (Algorithm 3):

   • The method employs a sequence of fixed-size Krylov subspaces.
   • At each restart cycle $k$, the process starts with the last vector of the previous cycle ($w_{m+1}^{(k-1)}$).
   • The global approximation is updated by combining the block-triangular structure of the accumulated projection matrices.
   • Crucially, the method avoids storing the full basis from previous cycles, maintaining low memory usage.

3. Key Contributions

1. Novel Algorithm: Introduction of a Randomized Restarted Krylov framework that integrates RGS-based Arnoldi iterations into the restarting cycle.
2. Acceleration of Restarted Methods: Demonstration that randomization not only reduces computational cost but can improve the convergence rate of restarted methods compared to classical restarted Arnoldi.
3. Analysis of Ill-Conditioned Problems: A comprehensive study of randomized methods on large-scale, ill-conditioned matrices (arising from Finite Element Method discretizations and Graph Laplacians), a regime where previous randomized studies were limited to well-conditioned or small problems.
4. Ritz Value Distribution Insight: The paper provides numerical evidence that randomization perturbs the Ritz values (eigenvalues of the projected matrix) in a way that prevents the "clustering" phenomenon often seen in classical restarted methods. This leads to a better representation of the matrix spectrum and faster convergence.
5. Scalability Analysis: Investigation into the parallel scalability of the proposed method, showing that randomization does not negatively impact strong scaling.

4. Numerical Results

The authors tested their method on three distinct problem classes:

• Convection-Diffusion Equation (3D FEM): Large, sparse, ill-conditioned matrices.
• Circular Membrane Vibration: Oscillatory problem involving the matrix cosine.
• Graph Laplacian Diffusion: Massive complex networks (Kronecker, Orkut, uk-2002).

Key Findings:

• Convergence: The Randomized Restarted (restart-rand) method consistently outperformed the classical Restarted (restart) method. In the Convection-Diffusion test, it converged significantly faster, requiring fewer restart cycles to reach the same accuracy.
• Stability:
  • Classical incomplete orthogonalization often diverged or stagnated.
  • The sFOM (sketched FOM) method showed instability in oscillatory problems.
  • The proposed restart-rand maintained stability and accuracy comparable to the full Arnoldi method but with lower cost.
• Performance (Speed):
  • Non-restarted: The randomized Arnoldi (rand) was up to 2.92x faster than classical Arnoldi due to reduced orthogonalization costs and better cache utilization (sketches fit in L3 cache).
  • Restarted: The speedup was more modest (1.29x to 1.56x) because the overhead of sketching became significant relative to the short restart cycles. However, because restart-rand converged in fewer cycles, the total time to solution was still lower.
• Parameter Sensitivity:
  • Sketching Dimension ($d$): Too small ($d < 160$) led to erratic convergence or divergence. Optimal performance was found around $d \approx 320$ for restart-rand.
  • Sparsity ($\zeta$): Lower sparsity ($\zeta=1$) generally offered better performance than $\zeta=4$ due to lower sketching costs, provided the basis remained well-conditioned.

5. Significance

• Practical Impact: The method offers a viable solution for computing matrix functions on extremely large, ill-conditioned systems (e.g., $N \approx 10^7$) where memory constraints force the use of restarting, but classical restarting is too slow.
• Theoretical Insight: It challenges the assumption that randomization is only a "preconditioning" tool for small problems, proving its efficacy in the high-dimensional, ill-conditioned regime.
• Future Directions: The paper highlights that while empirical results are strong, a formal theoretical analysis of how randomization affects Ritz value distribution and convergence rates is still needed. It also suggests combining this approach with spectral deflation or thick restarting for further gains.

In summary, the paper successfully demonstrates that random sketching can be effectively integrated into restarted Krylov methods to simultaneously reduce memory usage, lower computational cost, and accelerate convergence for difficult, large-scale scientific computing problems.