---
title: Randomized Krylov Subspace Methods
url: https://www.emergentmind.com/topics/randomized-krylov-subspace-methods
type: topic
---

# Randomized Krylov Subspace Methods

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 [2011.05090][1504.05477][2504.03914].

## 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 [2302.07466][2110.00649][2503.22631][2504.03914].

| Paradigm | Randomization mechanism | Representative tasks |
|---|---|---|
| Sketch-based Arnoldi/GMRES/FOM | Oblivious subspace embedding \(S\) or \(\Omega\) | Linear systems, eigenproblems, matrix functions [2011.05090][2302.07466] |
| Inner-product-free sketch-and-solve | Gaussian sketch in projected least-squares | Large-scale inverse problems [2502.02721] |
| Randomized block Krylov | Gaussian start \(\Omega\) or \(\Pi\) | Extreme eigenvalues, low-rank approximation, PCA [2110.00649][1504.05477] |
| Recycling with randomized sketching | Sketched compression of augmented Krylov basis | Sequences of matrix-function evaluations [2308.02290] |
| Randomized restarting/truncation | Randomized Arnoldi inside restart or adaptive truncation probabilities | Matrix functions, unbiased estimators [2503.22631][2504.03914] |

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 [2011.05090][2302.07466]. 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 [2011.05090].

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 \(u\perp^\Omega v\) when \(\langle \Omega u,\Omega v\rangle=0\), defines a sketched orthogonal projector \(P^\Omega_{V_k}\), and enforces a sketched Petrov–Galerkin condition on the residual [2302.07466]. For symmetric positive-definite \(A\), the resulting approximation satisfies an \(A\)-norm quasi-optimality bound relative to the classical orthogonal projection solution:
\[
\|x-x_k\|_A \le \frac{1+\epsilon\sqrt{\mathrm{cond}(A)}}{1-\epsilon\sqrt{\mathrm{cond}(A)}}\|x-\breve x_k\|_A,
\]
under the embedding assumption and \(\epsilon \,\mathrm{cond}(A)^{1/2}<1\) [2302.07466].

For GMRES, the sketch-based construction yields near-minimal residual behavior. In the randomized Gram–Schmidt analysis, if \(S\) embeds the Krylov space, then randomized GMRES satisfies
\[
\|Ax_k-b\| \le \kappa_2(SQ_k)\min_{v\in K_k(A,b)}\|Av-b\|
= O\!\left(\sqrt{\frac{1+\epsilon}{1-\epsilon}}\right)\min\|\cdot\|,
\]
so the method is near-optimal up to the embedding distortion [2011.05090]. The 2025 overview of randomized orthogonalization reformulates this through a “sketch-orthonormal” basis \(Q\) with \((\Omega Q)^T(\Omega Q)=I\), together with conditioning estimates of the form
\[
\mathrm{Cond}(Q)\le \sqrt{\frac{1+\epsilon}{1-\epsilon}}\mathrm{Cond}(S)=O(1),
\]
and shows how whitening via the small QR factorization of \(\Omega W\) lets one work with nonorthogonal short-recurrence bases in \(\mathbb R^n\) while retaining a well-conditioned sketched basis [2512.15455].

The algorithmic motivation is arithmetic and communication. Classical Arnoldi orthogonalization costs \(O(nk^2)\); randomized Arnoldi replaces the expensive \(O(nj)\) inner-product stage at step \(j\) by sketch application and \(O(\ell j)\) orthogonalization in \(\mathbb R^\ell\), plus a reconstruction in \(\mathbb R^n\) [2302.07466]. 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 [2011.05090][2512.15455].

## 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 \(L_k=[d_1,\dots,d_k]\) satisfying
\[
A L_k = L_{k+1} H_{k+1,k},
\]
without ever forming \(\langle \cdot,\cdot\rangle\); the process uses only coordinate evaluations and rank-one updates [2502.02721]. To obtain a minimal-residual-type projected solve, the small least-squares problem
\[
\min_y \|H_{k+1,k}y-\beta e_1\|_2
\]
is itself sketched with a Gaussian oblivious subspace embedding \(S\in\mathbb R^{\ell\times m}\), leading to
\[
y_k^{(S)}=\arg\min_y \|S(H_{k+1,k}y-\beta e_1)\|_2.
\]
With high probability, the sketched solution satisfies
\[
\|Hy_k^{(S)}-\beta e_1\|_2
\le \frac{1+\epsilon}{1-\epsilon}\min_y \|Hy-\beta e_1\|_2,
\]
so the randomized projected residual is theoretically closer to the true residual norm than in quasi-minimal residual inner-product-free methods [2502.02721]. 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 \(k\) [2502.02721].

The rectangular-matrix line is based on randomized Golub–Kahan bidiagonalization. Here one uses two sketches, \(S_m\in\mathbb R^{s_m\times m}\) and \(S_n\in\mathbb R^{s_n\times n}\), and builds sketched-orthonormal left and right bases \(U_{k+1}\) and \(V_{k+1}\) satisfying randomized analogues of
\[
A V_k = U_{k+1} B_k,\qquad A^T U_k = V_{k+1} B_k^T.
\]
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 [2508.20269]. Theoretical guarantees state that, with high probability, randomized residuals remain within \((1\pm\epsilon)\) of their classical counterparts, while experiments in image deblurring and seismic tomography report comparable solution quality with reduced orthogonalization cost [2508.20269].

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 \(A\) and, when needed, \(A^T\). 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 \(f(A)b\) has become a major testing ground for randomized Krylov ideas. For a general matrix \(A\), classical Arnoldi computes an orthonormal basis \(U_m\) and forms
\[
f(A)b \approx \|b\|_2 U_m f(K_m)e_1.
\]
The randomized alternative emphasizes that exact orthonormality is not essential if one can approximate the compression \(V_m^\dagger A V_m\) efficiently [2212.12758]. Two constructions are highlighted: sketch-and-orthonormalize and truncated orthogonalization. The former maintains an orthonormal sketched basis \(S_i=\Theta V_i\), and the latter orthogonalizes only against the last \(k\) vectors, optionally using a sketch to monitor \(\mathrm{cond}(\Theta V_i)\) and trigger whitening [2212.12758]. With an accurate least-squares correction, the randomized approximation is as accurate as Arnoldi in exact arithmetic for analytic \(f\) 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 [2212.12758].

Recycling changes the objective from a single matrix-function evaluation to a sequence \(f(A^{(i)})b^{(i)}\). Burke and Güttel derive a closed-form augmented FOM approximation
\[
F_m=\underline V_m\, f\!\left[\underline V_m^\dagger A \underline V_m\right]\underline V_m^\dagger,
\]
with \(\underline V_m=[U,V_m]\) combining a recycling space \(U\) and a fresh Krylov basis \(V_m\) [2308.02290]. The sketched version replaces full orthogonalization by a subspace embedding and yields
\[
\hat F_m=\underline V_m\, f\!\left[(S_m)^\dagger(SA_m)\right](S_m)^\dagger S,
\]
together with a whitening step based on a thin QR factorization of \(S_m=S\underline V_m\) [2308.02290]. In this formulation, recycling and randomized sketching are fused: spectral information is carried across problems through \(U\), 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 [2308.02290].

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 [2503.22631]. 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 [2503.22631]. 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 [2503.22631].

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 \(\|r_m(t)\|_*\) computable from the last projected block and leads to an a posteriori error estimate of variation-of-constants type; in the standard Euclidean setting,
\[
\|e_m(t)\|\le C_1\, t\, \phi_1(-\omega_1 t)\,\max_{0\le \theta\le t}\|r_m(\theta)\|,
\]
and in the sketched setting the bound is modified by the distortion factor \(1/\sqrt{1-\epsilon}\) [2510.17538]. The corresponding residual-based sketched Arnoldi method is reported to achieve speedups of \(2\times\)–\(5\times\) over classical FOM on large-scale ODE models while preserving a reliable stopping criterion [2510.17538].

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 [2504.03914]. 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 \(A\), one draws \(\Omega\in\mathbb R^{n\times \ell}\) with i.i.d. Gaussian entries and forms
\[
K_q(A;\Omega)=\mathrm{span}\{\Omega,A\Omega,A^2\Omega,\dots,A^q\Omega\}.
\]
The projected Rayleigh–Ritz estimate of the top eigenvalue is then computed on \(\mathrm{range}(S_q)\), where \(S_q=[\Omega\mid A\Omega\mid\cdots\mid A^q\Omega]\) [2110.00649]. The analysis of randomized block Krylov methods for extreme eigenvalues derives gap-free and spectral-gap bounds in terms of a stable-rank quantity \(\mathrm{srk}(\nu)\), 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 \(\ell\) entering delicately through the probability tails and expectation bounds [2110.00649]. Numerical evidence in the same paper shows that even \(\ell=2\) or \(3\) can significantly reduce misconvergence probability relative to \(\ell=1\) [2110.00649].

For low-rank approximation and PCA, Musco and Musco’s randomized block Krylov iteration forms
\[
K=[A\Pi,(AA^T)A\Pi,(AA^T)^2A\Pi,\dots,(AA^T)^qA\Pi]
\]
from a Gaussian start \(\Pi\), orthonormalizes \(K\), and extracts the best rank-\(k\) component in that space [1504.05477]. The central result is a gap-independent improvement over simultaneous power iteration: a \((1+\epsilon)\) spectral norm low-rank approximation is obtained in \(\tilde O(1/\sqrt{\epsilon})\) iterations instead of \(\tilde O(1/\epsilon)\), and the method also yields nearly optimal PCA guarantees for each principal component [1504.05477]. 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 \(A\), a Krylov-aware randomized low-rank approximation of \(f(A)\) runs one block-Lanczos process of depth \(s+r\) on a Gaussian start \(\Omega\), forms the small projected matrix \(T_{s+r}\), and returns
\[
Q_s\,[\,f(T_{s+r})\,]_{1:s\ell,\,1:s\ell}\,Q_s^T
\]
or its rank-\(k\) truncation [2502.01888]. The method “opens” the inner Krylov solver instead of treating \(f(A)\) as a black-box matrix–vector oracle. It uses exactly \((s+r)\ell\) matrix–vector products with \(A\), the same asymptotic cost as the naive outer randomized SVD with inner Krylov solves, but the effective projection space has dimension \(s\ell\) rather than \(\ell\) [2502.01888]. 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 \(\mathcal X*\mathcal X^\top\), computes a t-orthonormal basis via t-QR, and projects to a small tensor for T-SVD [2504.04989]. 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 [2603.00839]. 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
\[
d_k=-A^TS_kS_k^T(Ax^k-b),
\]
and the new iterate is chosen in the affine subspace
\[
\Pi_k=\mathrm{aff}\{x^{j_k},\dots,x^k,x^k+d_k\},
\]
with the coefficients determined by an optimal projection mechanism [2505.20602]. As the memory parameter \(\ell\to\infty\) and \(S_k\equiv I\), the method recovers the usual normal-equation Krylov subspace. At \(\ell=1\), it recovers stochastic gradient or randomized Kaczmarz with an adaptive Polyak step size [2505.20602]. The method is proved to converge linearly in expectation, thereby giving a literal bridge between low-memory randomized iterations and full-memory Krylov schemes [2505.20602].

For symmetric positive-definite systems, a related perspective appears in randomized block conjugate gradient. Here one solves a single system \((A+\mu I)x=b\) but augments \(b\) with \(m-1\) Gaussian vectors and runs block Lanczos on the resulting block \(B=[b\mid \Omega]\) [2501.18717]. The key approximation is
\[
x_t(\mu)=Q_t (T_t+\mu I)^{-1}E_1 B_{-1},
\]
and the same basis serves all regularization parameters \(\mu\) because \(\mathcal K_t(A_\mu,B)=\mathcal K_t(A,B)\) [2501.18717]. 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 [2501.18717].

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 \(\alpha_k\) and \(\beta_k\) is introduced precisely to capture these spikes [2302.07466]. For restarted randomized methods for matrix functions, the papers likewise report occasional error spikes and note that too small a sketch dimension \(d\) may lead to erratic behavior or stagnation [2503.22631].

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 [2601.10248]. 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,
\[
\widehat h_m=(U_m^HU_m)^{-1}(U_m^Hu_{m+1}),
\]
forms
\[
\widehat u_{m+1}=u_{m+1}-U_m\widehat h_m,\qquad
\widehat H_m=H_m+\widehat h_m c_m^H,
\]
and thereby enforces Euclidean orthogonality of the last Arnoldi vector to the current subspace [2601.10248]. 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 [2601.10248].

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 [2212.12758]. For restarted randomized Krylov methods, a formal account of how random sketching perturbs Ritz values and improves approximation is also listed as lacking [2503.22631]. 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 [2504.03914].

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 [2601.10248][2504.03914].

Source: https://www.emergentmind.com/topics/randomized-krylov-subspace-methods