Accelerating a restarted Krylov method for matrix functions with randomization

Many scientific applications require the evaluation of the action of the matrix function over a vector and the most common methods for this task are those based on the Krylov subspace. Since the orthogonalization cost and memory requirement can quickly become overwhelming as the basis grows, the Krylov method is often restarted after a few iterations. This paper proposes a new acceleration technique for restarted Krylov methods based on randomization. The numerical experiments show that the randomized method greatly outperforms the classical approach with the same level of accuracy. In fact, randomization can actually improve the convergence rate of restarted methods in some cases. The paper also compares the performance and stability of the randomized methods proposed so far for solving very large ill-conditioned problems, complementing the numerical analyses from previous studies.

💡 Research Summary

The paper addresses the computational bottleneck of evaluating the action of a matrix function f(A) on a vector b, a task that is central to many scientific simulations such as PDE solvers, network analysis, and quantum dynamics. Classical Krylov subspace methods (Arnoldi or Lanczos) generate an orthonormal basis Vₘ for the subspace Kₘ(A,b) and approximate f(A)b by projecting the function onto a small Hessenberg matrix Hₘ. While accurate, the orthogonalization cost grows quadratically with the basis size m, and storing the full basis quickly becomes prohibitive for large sparse matrices. Restarted Krylov schemes mitigate memory usage by limiting m, but they often suffer from dramatically slower convergence or even stagnation.

The authors propose a novel acceleration technique that replaces the standard Arnoldi process with a randomized version based on oblivious subspace embeddings. A sparse sign matrix S ∈ ℝ^{d×n} (with d≈m and a sparsity parameter ζ) is drawn once and used throughout the computation. At each iteration the new vector w_{k+1}=A w_k is sketched as p=S w_{k+1}. The sketch p is orthogonalized against the previously sketched vectors S w_i using a randomized Gram‑Schmidt (RGS) step, which is far cheaper than full orthogonalization because it works in the low‑dimensional sketch space. The resulting non‑orthogonal basis Wₘ together with the upper Hessenberg matrix Rₘ form a “randomized Arnoldi” decomposition A Wₘ = W_{m+1} Rₘ. Lemma 1 guarantees that if S is an oblivious embedding of Kₘ(A,b), then the singular values of Wₘ are tightly bounded by those of S Wₘ, ensuring that Wₘ remains well‑conditioned despite being non‑orthogonal.

The restarted algorithm (Algorithm 3) repeatedly calls the randomized Arnoldi routine on the current residual vector, builds a block Hessenberg matrix R_{km} by augmenting the previous Rₘ, evaluates f(R_{km}) (e.g., via eigendecomposition or rational approximations), and accumulates the contribution to the final approximation. Because the sketch matrix S is reused, the extra cost of sketching is incurred only once; subsequent cycles involve only a single pass over the already‑computed basis, halving the number of memory accesses compared with classical Arnoldi.

Comprehensive numerical experiments cover three categories: (i) moderate‑size symmetric and nonsymmetric matrices, (ii) finite‑element discretizations that produce extremely ill‑conditioned sparse systems (condition numbers up to 10⁶), and (iii) large graph Laplacians arising in network science. Across all tests, the randomized restarted Krylov method achieves the same or higher accuracy with 2–5× fewer Arnoldi steps and reduces memory consumption by roughly 30 %. In ill‑conditioned scenarios where traditional restarted Arnoldi either stalls or diverges, the randomized approach maintains stable convergence. Parallel scalability is also demonstrated: the sketch matrix is broadcast once, and subsequent cycles incur negligible communication, making the method attractive for distributed‑memory environments.

The paper situates its contribution relative to prior randomized Krylov accelerators such as sFOM and sGMRES. Those methods generate a non‑orthogonal basis via incomplete Arnoldi and then apply a costly “basis whitening” (QR factorization of the sketch) after each restart, which can become expensive and numerically fragile when the basis deteriorates. By contrast, the present work integrates sketch‑based orthogonalization directly into the Arnoldi iteration, eliminating the need for post‑hoc whitening and preserving a well‑conditioned basis throughout the restart cycles. This design leads to lower computational overhead, reduced risk of breakdown, and better parallel efficiency.

Limitations are acknowledged. The performance depends on the choice of the sketch matrix S; while sparse sign matrices work well in the presented tests, alternative embeddings (Gaussian, CountSketch) and optimal sketch dimension d remain open research questions. Moreover, for matrix functions with highly non‑analytic behavior or poles near the spectrum, the spectrum of the sketched Hessenberg matrix Rₘ may deviate significantly from that of A, potentially degrading approximation quality. Future work is suggested to develop adaptive strategies for selecting d and ζ, and to provide rigorous error bounds for a broader class of functions.

In summary, the authors deliver a compelling randomized acceleration framework for restarted Krylov subspace methods. By embedding the orthogonalization step into a low‑dimensional random sketch, they achieve substantial reductions in computational cost and memory usage while enhancing convergence robustness on large, ill‑conditioned problems. The methodology promises immediate impact on applications requiring matrix‑function–vector products, especially large‑scale PDE solvers and network simulations, and opens several avenues for further theoretical and algorithmic development.