Randomized Krylov Subspace Methods

• Randomized Krylov subspace methods are a family of algorithms that incorporate randomness into basis construction, compression, or truncation to preserve approximation power while reducing computational and communication costs.
• They employ techniques like sketch-based orthogonalization, randomized Arnoldi, and block Krylov iterations to achieve faster convergence and lower memory pressure in solving linear systems and eigenproblems.
• These methods are applied to inverse problems, matrix functions, and low-rank approximations, demonstrating practical speedups on distributed and mixed-precision architectures.

Randomized Krylov subspace methods are Krylov algorithms in which randomness enters the construction, compression, or truncation of the approximation space. The common objective is to preserve the approximation power of deterministic Krylov procedures while reducing orthogonalization cost, communication, or memory pressure, or while enlarging the effective search space for a fixed matrix–vector budget. In the recent literature, this umbrella includes sketch-based Arnoldi and Golub–Kahan procedures, randomized block Krylov methods with Gaussian starting blocks, sketch-and-solve projected least-squares schemes, recycling and restarted methods for matrix functions, and adaptive randomized truncation estimators that remain unbiased (Balabanov et al., 2020, Musco et al., 2015, Luo et al., 4 Apr 2025).

1. Taxonomy and problem classes

A Krylov subspace method builds approximations in spaces of the form

$$\mathcal K_m(A,b)=\mathrm{span}\{b,Ab,A^2b,\dots,A^{m-1}b\}.$$

Randomization modifies this template in several distinct ways. In sketch-based methods, a random embedding replaces exact inner products by sketched inner products and yields a basis whose sketch is orthonormal. In randomized block Krylov methods, one starts from a Gaussian block $\Omega$ or $\Pi$ and works in a larger block-Krylov space. In sketch-and-solve variants, the projected least-squares problem itself is randomized. In restarted and recycled methods for matrix functions, randomization is used to control orthogonalization or to improve restarted convergence. A recent line introduces adaptive randomized truncation for unbiased Krylov estimators (Timsit et al., 2023, Tropp, 2021, Guidotti et al., 28 Mar 2025, Luo et al., 4 Apr 2025).

Paradigm	Randomization mechanism	Representative tasks
Sketch-based Arnoldi/GMRES/FOM	Oblivious subspace embedding $S$ or $\Omega$	Linear systems, eigenproblems, matrix functions (Balabanov et al., 2020, Timsit et al., 2023)
Inner-product-free sketch-and-solve	Gaussian sketch in projected least-squares	Large-scale inverse problems (Landman et al., 4 Feb 2025)
Randomized block Krylov	Gaussian start $\Omega$ or $\Pi$	Extreme eigenvalues, low-rank approximation, PCA (Tropp, 2021, Musco et al., 2015)
Recycling with randomized sketching	Sketched compression of augmented Krylov basis	Sequences of matrix-function evaluations (Burke et al., 2023)
Randomized restarting/truncation	Randomized Arnoldi inside restart or adaptive truncation probabilities	Matrix functions, unbiased estimators (Guidotti et al., 28 Mar 2025, Luo et al., 4 Apr 2025)

This taxonomy reflects a central distinction. Some methods randomize the basis construction itself, whereas others randomize the projected model, the start block, or the stopping and truncation mechanism. The literature therefore treats “randomized Krylov” not as a single algorithm but as a family of closely related design principles.

2. Sketch-based orthogonalization and randomized Arnoldi

A dominant thread replaces full-dimensional orthogonalization by orthogonalization in a smaller sketch space. The basic ingredient is an oblivious subspace embedding $S\in\mathbb R^{k\times n}$ or $\Omega\in\mathbb R^{\ell\times n}$ such that, for vectors in the relevant Krylov space,

$$(1-\varepsilon)\|x\|^2 \le \|Sx\|^2 \le (1+\varepsilon)\|Sx\|^2,$$

or equivalently the sketched inner product approximates the exact one (Balabanov et al., 2020, Timsit et al., 2023). Randomized Gram–Schmidt then orthonormalizes sketches and updates the high-dimensional vectors using the resulting small least-squares coefficients. Balabanov and Grigori’s randomized Gram–Schmidt process is designed for orthonormalization or QR factorization and is stated to be less computationally expensive than classical Gram–Schmidt while being at least as numerically stable as modified Gram–Schmidt (Balabanov et al., 2020).

When inserted into Arnoldi, the result is a randomized Arnoldi relation in which the basis is orthonormal in the sketched inner product rather than in the Euclidean one. In the randomized orthogonal projection framework, one writes $\Omega$0 when $\Omega$1, defines a sketched orthogonal projector $\Omega$2, and enforces a sketched Petrov–Galerkin condition on the residual (Timsit et al., 2023). For symmetric positive-definite $\Omega$3, the resulting approximation satisfies an $\Omega$4-norm quasi-optimality bound relative to the classical orthogonal projection solution: $\Omega$5 under the embedding assumption and $\Omega$6 (Timsit et al., 2023).

For GMRES, the sketch-based construction yields near-minimal residual behavior. In the randomized Gram–Schmidt analysis, if $\Omega$7 embeds the Krylov space, then randomized GMRES satisfies

$\Omega$8

so the method is near-optimal up to the embedding distortion (Balabanov et al., 2020). The 2025 overview of randomized orthogonalization reformulates this through a “sketch-orthonormal” basis $\Omega$9 with $\Pi$0, together with conditioning estimates of the form

$\Pi$1

and shows how whitening via the small QR factorization of $\Pi$2 lets one work with nonorthogonal short-recurrence bases in $\Pi$3 while retaining a well-conditioned sketched basis (Damas et al., 17 Dec 2025).

The algorithmic motivation is arithmetic and communication. Classical Arnoldi orthogonalization costs $\Pi$4; randomized Arnoldi replaces the expensive $\Pi$5 inner-product stage at step $\Pi$6 by sketch application and $\Pi$7 orthogonalization in $\Pi$8, plus a reconstruction in $\Pi$9 (Timsit et al., 2023). The literature repeatedly emphasizes that this is especially attractive on distributed-memory and mixed-precision architectures, because global dot products are replaced by smaller dense kernels and local updates (Balabanov et al., 2020, Damas et al., 17 Dec 2025).

3. Randomized Krylov methods for linear systems and inverse problems

For inverse problems, the randomized Krylov literature splits into two complementary directions. One keeps the projection framework but removes inner products from basis generation. The other extends sketch-based projections from square to rectangular operators.

The inner-product-free line is developed through a Hessenberg process with partial pivoting that generates a nonorthogonal Krylov basis $S$0 satisfying

$S$1

without ever forming $S$2; the process uses only coordinate evaluations and rank-one updates (Landman et al., 4 Feb 2025). To obtain a minimal-residual-type projected solve, the small least-squares problem

$S$3

is itself sketched with a Gaussian oblivious subspace embedding $S$4, leading to

$S$5

With high probability, the sketched solution satisfies

$S$6

so the randomized projected residual is theoretically closer to the true residual norm than in quasi-minimal residual inner-product-free methods (Landman et al., 4 Feb 2025). Numerically, the resulting sCMRH and sLSLU methods achieve residual norms and reconstruction errors nearly identical to GMRES and LSQR, while removing all inner products and reducing global communication by an order proportional to the iteration count $S$7 (Landman et al., 4 Feb 2025).

The rectangular-matrix line is based on randomized Golub–Kahan bidiagonalization. Here one uses two sketches, $S$8 and $S$9, and builds sketched-orthonormal left and right bases $\Omega$0 and $\Omega$1 satisfying randomized analogues of

$\Omega$2

This gives randomized LSQR-, CGLS-, and LSMR-type solvers together with hybrid projection methods that combine randomized Arnoldi or randomized Golub–Kahan factorizations with Tikhonov regularization, and select regularization parameters automatically during the iteration via discrepancy principle, GCV, or WGCV (Chung et al., 27 Aug 2025). Theoretical guarantees state that, with high probability, randomized residuals remain within $\Omega$3 of their classical counterparts, while experiments in image deblurring and seismic tomography report comparable solution quality with reduced orthogonalization cost (Chung et al., 27 Aug 2025).

A common structural feature of both directions is that the randomized component acts only on the small projected geometry. The large-scale information still enters through matrix–vector products with $\Omega$4 and, when needed, $\Omega$5. This suggests that the randomized inverse-problem solvers remain recognizably Krylov in approximation mechanism, even when orthogonality is substantially relaxed.

4. Matrix functions, recycling, restarting, and residual-based control

The evaluation of $\Omega$6 has become a major testing ground for randomized Krylov ideas. For a general matrix $\Omega$7, classical Arnoldi computes an orthonormal basis $\Omega$8 and forms

$\Omega$9

The randomized alternative emphasizes that exact orthonormality is not essential if one can approximate the compression $\Omega$0 efficiently (Cortinovis et al., 2022). Two constructions are highlighted: sketch-and-orthonormalize and truncated orthogonalization. The former maintains an orthonormal sketched basis $\Omega$1, and the latter orthogonalizes only against the last $\Omega$2 vectors, optionally using a sketch to monitor $\Omega$3 and trigger whitening (Cortinovis et al., 2022). With an accurate least-squares correction, the randomized approximation is as accurate as Arnoldi in exact arithmetic for analytic $\Omega$4 on a convex domain containing the numerical range; a cheaper variant omits the least-squares correction and can be much faster, but may converge erratically if the last correction column matters (Cortinovis et al., 2022).

Recycling changes the objective from a single matrix-function evaluation to a sequence $\Omega$5. Burke and Güttel derive a closed-form augmented FOM approximation

$\Omega$6

with $\Omega$7 combining a recycling space $\Omega$8 and a fresh Krylov basis $\Omega$9 (Burke et al., 2023). The sketched version replaces full orthogonalization by a subspace embedding and yields

$\Pi$0

together with a whitening step based on a thin QR factorization of $\Pi$1 (Burke et al., 2023). In this formulation, recycling and randomized sketching are fused: spectral information is carried across problems through $\Pi$2, while the orthogonalization bottleneck is shifted into the sketch space. The paper states that this is the first method to combine subspace recycling with randomized sketching for matrix functions and reports reduced orthogonalization cost and storage relative to classical FOM-based recycling (Burke et al., 2023).

Restarted methods add another layer. Since orthogonalization cost and memory grow with the basis, matrix-function solvers are often restarted after a fixed number of Arnoldi steps. Guidotti, Martinsson, Acebrón, and Monteiro replace standard Arnoldi inside each restart cycle by randomized Arnoldi, producing a method denoted restart-rand (Guidotti et al., 28 Mar 2025). Their numerical experiments on very large finite element problems show that randomization can improve the convergence rate of restarted methods in some cases, with reported speedups of 15–50% over classical restart and a better spread of Ritz values over the spectral interval (Guidotti et al., 28 Mar 2025). The paper explicitly interprets the improved behavior through a de-clustering of Ritz values in restarted cycles, although it also states that a formal spectral analysis is still lacking (Guidotti et al., 28 Mar 2025).

For exponential integration and related matrix-function ODEs, residual-based stopping has recently been extended to the sketched setting. A general framework for Krylov ODE residuals yields a sketched residual norm $\Pi$3 computable from the last projected block and leads to an a posteriori error estimate of variation-of-constants type; in the standard Euclidean setting,

$\Pi$4

and in the sketched setting the bound is modified by the distortion factor $\Pi$5 (Krieger et al., 20 Oct 2025). The corresponding residual-based sketched Arnoldi method is reported to achieve speedups of $\Pi$6–$\Pi$7 over classical FOM on large-scale ODE models while preserving a reliable stopping criterion (Krieger et al., 20 Oct 2025).

A separate 2025 direction introduces adaptive randomized truncation estimators for Krylov subspace methods. The estimator remains unbiased while optimizing the trade-off between solution variance and computational cost, computes truncation probabilities on the fly through a constrained optimization problem, has a closed-form solution under a diminishing returns property, and proves that the optimal adaptive truncation distribution is impossible to obtain in the general case (Luo et al., 4 Apr 2025). This moves randomized Krylov methods beyond sketching alone and into unbiased stochastic estimation.

5. Randomized block Krylov methods for eigenvalues, singular values, and low-rank approximation

Randomized block Krylov methods use a random starting block rather than a random sketch of the basis. For a symmetric matrix $\Pi$8, one draws $\Pi$9 with i.i.d. Gaussian entries and forms

$S\in\mathbb R^{k\times n}$0

The projected Rayleigh–Ritz estimate of the top eigenvalue is then computed on $S\in\mathbb R^{k\times n}$1, where $S\in\mathbb R^{k\times n}$2 (Tropp, 2021). The analysis of randomized block Krylov methods for extreme eigenvalues derives gap-free and spectral-gap bounds in terms of a stable-rank quantity $S\in\mathbb R^{k\times n}$3, and shows that for matrices with polynomial spectral decay an accurate spectral norm estimate can be obtained using only a constant number of steps, with the dependence on block size $S\in\mathbb R^{k\times n}$4 entering delicately through the probability tails and expectation bounds (Tropp, 2021). Numerical evidence in the same paper shows that even $S\in\mathbb R^{k\times n}$5 or $S\in\mathbb R^{k\times n}$6 can significantly reduce misconvergence probability relative to $S\in\mathbb R^{k\times n}$7 (Tropp, 2021).

For low-rank approximation and PCA, Musco and Musco’s randomized block Krylov iteration forms

$S\in\mathbb R^{k\times n}$8

from a Gaussian start $S\in\mathbb R^{k\times n}$9, orthonormalizes $\Omega\in\mathbb R^{\ell\times n}$0, and extracts the best rank-$\Omega\in\mathbb R^{\ell\times n}$1 component in that space (Musco et al., 2015). The central result is a gap-independent improvement over simultaneous power iteration: a $\Omega\in\mathbb R^{\ell\times n}$2 spectral norm low-rank approximation is obtained in $\Omega\in\mathbb R^{\ell\times n}$3 iterations instead of $\Omega\in\mathbb R^{\ell\times n}$4, and the method also yields nearly optimal PCA guarantees for each principal component (Musco et al., 2015). The key proof idea is polynomial filtering: block Krylov spaces can approximate Chebyshev-type polynomials that separate large and small singular values more rapidly than monomials.

The same block-Krylov philosophy has been transferred to matrix functions. For symmetric $\Omega\in\mathbb R^{\ell\times n}$5, a Krylov-aware randomized low-rank approximation of $\Omega\in\mathbb R^{\ell\times n}$6 runs one block-Lanczos process of depth $\Omega\in\mathbb R^{\ell\times n}$7 on a Gaussian start $\Omega\in\mathbb R^{\ell\times n}$8, forms the small projected matrix $\Omega\in\mathbb R^{\ell\times n}$9, and returns

$$(1-\varepsilon)\|x\|^2 \le \|Sx\|^2 \le (1+\varepsilon)\|Sx\|^2,$$0

or its rank-$$(1-\varepsilon)\|x\|^2 \le \|Sx\|^2 \le (1+\varepsilon)\|Sx\|^2,$$1 truncation (Persson et al., 3 Feb 2025). The method “opens” the inner Krylov solver instead of treating $$(1-\varepsilon)\|x\|^2 \le \|Sx\|^2 \le (1+\varepsilon)\|Sx\|^2,$$2 as a black-box matrix–vector oracle. It uses exactly $$(1-\varepsilon)\|x\|^2 \le \|Sx\|^2 \le (1+\varepsilon)\|Sx\|^2,$$3 matrix–vector products with $$(1-\varepsilon)\|x\|^2 \le \|Sx\|^2 \le (1+\varepsilon)\|Sx\|^2,$$4, the same asymptotic cost as the naive outer randomized SVD with inner Krylov solves, but the effective projection space has dimension $$(1-\varepsilon)\|x\|^2 \le \|Sx\|^2 \le (1+\varepsilon)\|Sx\|^2,$$5 rather than $$(1-\varepsilon)\|x\|^2 \le \|Sx\|^2 \le (1+\varepsilon)\|Sx\|^2,$$6 (Persson et al., 3 Feb 2025). The resulting Frobenius-norm bounds are correspondingly sharper.

Tensor extensions show that randomized Krylov ideas are not confined to matrices. A randomized block Krylov method for truncated tensor SVD constructs a tensor block-Krylov space by repeated application of $$(1-\varepsilon)\|x\|^2 \le \|Sx\|^2 \le (1+\varepsilon)\|Sx\|^2,$$7, computes a t-orthonormal basis via t-QR, and projects to a small tensor for T-SVD (Kooshkghazi et al., 7 Apr 2025). Separately, randomized tensor global GMRES based on a sketched Einstein inner product replaces Frobenius inner products by mode-wise random projections and derives residual bounds and projected Tikhonov regularization for image and video restoration (Badahmane, 28 Feb 2026). These developments indicate that randomized Krylov methodology has become a template rather than a matrix-specific trick.

6. Hybridizations with randomized iterative methods and block conjugate gradients

Another recent development links randomized iterative methods directly to Krylov subspaces. In the iterative-sketching framework of “Connecting randomized iterative methods with Krylov subspaces,” the update direction is the sketched normal vector

$$(1-\varepsilon)\|x\|^2 \le \|Sx\|^2 \le (1+\varepsilon)\|Sx\|^2,$$8

and the new iterate is chosen in the affine subspace

$$(1-\varepsilon)\|x\|^2 \le \|Sx\|^2 \le (1+\varepsilon)\|Sx\|^2,$$9

with the coefficients determined by an optimal projection mechanism (Sun et al., 27 May 2025). As the memory parameter $\Omega$00 and $\Omega$01, the method recovers the usual normal-equation Krylov subspace. At $\Omega$02, it recovers stochastic gradient or randomized Kaczmarz with an adaptive Polyak step size (Sun et al., 27 May 2025). The method is proved to converge linearly in expectation, thereby giving a literal bridge between low-memory randomized iterations and full-memory Krylov schemes (Sun et al., 27 May 2025).

For symmetric positive-definite systems, a related perspective appears in randomized block conjugate gradient. Here one solves a single system $\Omega$03 but augments $\Omega$04 with $\Omega$05 Gaussian vectors and runs block Lanczos on the resulting block $\Omega$06 (Chen et al., 30 Jan 2025). The key approximation is

$\Omega$07

and the same basis serves all regularization parameters $\Omega$08 because $\Omega$09 (Chen et al., 30 Jan 2025). The paper proves that this randomized block CG method can outperform a broad class of Nyström-preconditioned CG methods without explicitly constructing a preconditioner, and uses the same block-Lanczos run to compute the entire ridge-regression regularization path and to generate multiple Gaussian samples (Chen et al., 30 Jan 2025).

These two lines sharpen the meaning of “randomized Krylov.” Randomization need not only approximate orthogonality; it can also define a projection-based update rule whose limiting cases interpolate between stochastic first-order methods and Krylov acceleration, or it can embed preconditioning implicitly inside a block Krylov process.

7. Convergence irregularities, similarity restoration, and open issues

The strongest recurring qualification in the literature concerns irregular convergence. In randomized orthogonal projection methods, numerical experiments show convergence rates similar to standard orthogonal projection methods but with occasional spikes, especially early in the iteration sequence or near spectral cluster transitions; the sharper quasi-optimality bound in terms of computable factors $\Omega$10 and $\Omega$11 is introduced precisely to capture these spikes (Timsit et al., 2023). For restarted randomized methods for matrix functions, the papers likewise report occasional error spikes and note that too small a sketch dimension $\Omega$12 may lead to erratic behavior or stagnation (Guidotti et al., 28 Mar 2025).

A structural explanation is that the Hessenberg matrix produced by randomized Arnoldi is generally not similar to the one produced by standard Arnoldi. This lack of similarity can delay convergence or create spike-like irregularities in eigenvalue and matrix-function computations (Grigori et al., 15 Jan 2026). The remedy proposed in 2026 is a similarity-restoring randomized Arnoldi process: after the usual randomized Arnoldi cycle, one solves one additional least-squares problem,

$\Omega$13

forms

$\Omega$14

and thereby enforces Euclidean orthogonality of the last Arnoldi vector to the current subspace (Grigori et al., 15 Jan 2026). In exact arithmetic, the resulting projected matrix is similar to the standard Arnoldi Hessenberg, and for eigenvalue problems and matrix functions the approximations become identical to those of standard Arnoldi while retaining the speed of randomized orthogonalization (Grigori et al., 15 Jan 2026).

Several limitations remain explicit. For the randomized matrix-function schemes based on fast nonorthonormal bases, theoretical forward stability in finite precision is identified as an open problem (Cortinovis et al., 2022). For restarted randomized Krylov methods, a formal account of how random sketching perturbs Ritz values and improves approximation is also listed as lacking (Guidotti et al., 28 Mar 2025). The adaptive randomized truncation work adds a different kind of obstruction: it proves that obtaining the optimal adaptive truncation distribution is impossible in the general case, and only under a diminishing returns condition does a closed-form optimum exist (Luo et al., 4 Apr 2025).

These issues define the current research frontier. The field has already established quasi-optimality, near-minimal residuality, gap-independent spectral approximation, and practical speedups across linear systems, inverse problems, matrix functions, ODEs, and tensor problems. At the same time, similarity preservation, finite-precision behavior, restart stability, and the statistical design of unbiased truncation mechanisms remain active technical questions rather than settled doctrine (Grigori et al., 15 Jan 2026, Luo et al., 4 Apr 2025).