Randomized Preconditioner Methods

• Randomized preconditioners are techniques that use random sampling, sketches, and projections to construct operators that enhance the numerical conditioning of linear systems.
• They employ methods like matrix probing, Nyström spectral corrections, and randomized deflation to approximate inverse operators efficiently.
• Applications include seismic inversion, image reconstruction, and data assimilation, where they reduce iterative solver iterations and overall computational cost.

A randomized preconditioner is a preconditioning construction in which random sampling, randomized sketches, randomized trial functions, or randomized subspaces are used to build, approximate, or implicitly realize an operator that improves the conditioning or spectral clustering of a target problem. In contemporary numerical linear algebra and inverse problems, the term covers several distinct mechanisms: matrix probing of inverse operators, Nyström-type low-rank spectral corrections, randomized deflation, randomized approximate factorizations, and sketch-based diagnostics for preconditioner choice. These constructions are especially prominent when the forward model is available only through operator applications, when the Hessian or normal matrix has effective low rank, or when explicit factorization is too expensive (Demanet et al., 2011, Frangella et al., 2021, Balabanov et al., 24 Sep 2025).

1. Concept and taxonomic scope

The literature does not treat randomized preconditioning as a single algorithmic template. Instead, it comprises a family of techniques in which randomness is introduced at the level of range discovery, operator probing, sampling of algebraic structure, or preconditioner assessment. In some cases the output is an explicit SPD preconditioner; in others it is a projector, a low-rank correction, a sparse factor, or an implicit Krylov-space surrogate.

Family	Random ingredient	Representative formulation
Matrix probing	Randomized trial functions in curvelet space	$H^{-1}\approx \sum_j c_j B_j$ (Demanet et al., 2011)
Sketch-based spectral correction	Gaussian sketch $\Omega$, Nyström approximation	Nyström PCG (Frangella et al., 2021)
Range deflation	Projector onto $\operatorname{range}(A_\mu \Omega)$	RandRAND (Balabanov et al., 24 Sep 2025)
Randomized factorization	Sampled fill-ins or sketched triangular factors	ParAC (Liang et al., 5 May 2025), MRCQR (Garrison et al., 16 Jun 2026)
Randomized diagnostics	Gaussian probe vectors	$S(M)=\|I-M^{-1}A\|_{\mathsf F}$ (DiPaolo et al., 2019)

This scope also includes closely related methods in which the preconditioner is not formed explicitly but is emulated by randomized subspace iteration. The block Krylov construction in “Preconditioning without a preconditioner” is the clearest example: the augmented block-CG iterate is analyzed as an implicit randomized preconditioner that can match or outperform Nyström-based deflation without explicitly constructing the preconditioner (Chen et al., 30 Jan 2025).

2. Principal construction paradigms

A first major paradigm is operator probing. In the wave-equation Hessian setting, the inverse Hessian is modeled as a pseudodifferential operator with symbol $h^{-1}(x,\xi)$, approximated by a small basis expansion $H^{-1}\approx A(c):=\sum_{j=1}^{p} c_j B_j$. The coefficients are identified by least-squares fitting from applications of the normal operator $H$ to randomized trial functions built in curvelet space, typically with i.i.d. Gaussian or Rademacher coefficients on an illumination mask. The fitting relation is $A(c)H u_i \approx u_i$, and the resulting preconditioner is cheap to apply once the symbol coefficients are known (Demanet et al., 2011).

A second paradigm is randomized low-rank spectral approximation. For SPD systems $(A+\mu I)x=b$, randomized Nyström methods form $C=A\Omega$, $\Omega$0, and $\Omega$1, then define a preconditioner from the eigendecomposition $\Omega$2. In Nyström PCG, the preconditioner is

$\Omega$3

with inverse

$\Omega$4

This form explicitly rescales the sketched dominant eigenspace while leaving its orthogonal complement untouched (Frangella et al., 2021).

A third paradigm is randomized projection or deflation. RandRAND constructs an orthogonal projector $\Omega$5 onto $\Omega$6 and then defines right, correction, or left preconditioners that replace the action of $\Omega$7 on $\Omega$8 by a controlled scalar or projector-based correction. The method is designed to avoid eigenpair computation and even low-rank approximation in the usual Nyström sense; its basic object is the randomized deflating projector itself (Balabanov et al., 24 Sep 2025).

A fourth paradigm randomizes algebraic structure. For Laplacian systems, ParAC performs approximate Cholesky elimination in which dense clique updates are replaced by sampled sparse updates so that clique contributions are preserved in expectation while fill is dramatically reduced. In tall-skinny QR, randomized XR-based and SRFT-based constructions compute a sketched triangular factor $\Omega$9 or $\operatorname{range}(A_\mu \Omega)$0 and use right preconditioning $\operatorname{range}(A_\mu \Omega)$1 or $\operatorname{range}(A_\mu \Omega)$2 to stabilize subsequent CholeskyQR steps (Liang et al., 5 May 2025, Fan et al., 2021, Garrison et al., 16 Jun 2026).

3. Probabilistic analysis and conditioning theory

The analytical core of randomized preconditioning is typically a high-probability or in-expectation statement about spectrum, stability, or subspace capture. One line of work studies preconditioner quality through the stability proxy

$\operatorname{range}(A_\mu \Omega)$3

If $\operatorname{range}(A_\mu \Omega)$4 and $\operatorname{range}(A_\mu \Omega)$5, then $\operatorname{range}(A_\mu \Omega)$6, which leads to a dimension-free Gaussian sketch estimator. The same paper gives a tight tail bound and a union-bound guarantee for selecting the minimum-stability preconditioner among $\operatorname{range}(A_\mu \Omega)$7 candidates with work comparable to about $\operatorname{range}(A_\mu \Omega)$8 conjugate-gradient iterations (DiPaolo et al., 2019).

A more structural line of analysis studies randomized preconditioning without a deterministic sandwich $\operatorname{range}(A_\mu \Omega)$9. “Preconditioning in Expectation” shows that sampled Laplacian preconditioners can decrease energy in expectation even when the standard deterministic requirement fails. For randomized Richardson iteration, the paper proves a constant-factor contraction in expectation and combines this with recursive preconditioning to obtain algorithms for SDD systems and electrical flow with expected time close to $S(M)=\|I-M^{-1}A\|_{\mathsf F}$0 (Cohen et al., 2014).

For operator probing, the key random object is the least-squares design matrix. In matrix probing of the wave-equation Hessian, invertibility is controlled by a probabilistic condition of the form

$S(M)=\|I-M^{-1}A\|_{\mathsf F}$1

with high probability, where $S(M)=\|I-M^{-1}A\|_{\mathsf F}$2 is the rank of the informative subspace and $S(M)=\|I-M^{-1}A\|_{\mathsf F}$3 is the number of fitted symbol parameters. This formalizes the empirical observation that richer phase-space coverage permits larger symbol expansions (Demanet et al., 2011).

For Nyström preconditioning, the dominant quantity is the effective dimension

$S(M)=\|I-M^{-1}A\|_{\mathsf F}$4

With Gaussian $S(M)=\|I-M^{-1}A\|_{\mathsf F}$5 and $S(M)=\|I-M^{-1}A\|_{\mathsf F}$6, Nyström PCG proves

$S(M)=\|I-M^{-1}A\|_{\mathsf F}$7

and therefore constant-iteration CG convergence once the sketch rank is comparable to $S(M)=\|I-M^{-1}A\|_{\mathsf F}$8 (Frangella et al., 2021).

RandRAND pushes this perspective further. Its condition number bounds depend only weakly on the problem size and reduce to a small constant when the dimension of the deflated subspace is comparable to the effective spectral dimension. The emphasis is on projection error $S(M)=\|I-M^{-1}A\|_{\mathsf F}$9, rather than on explicit low-rank approximation error, which differentiates it from standard Nyström analyses (Balabanov et al., 24 Sep 2025).

4. Computational realization and implementation regimes

Randomized preconditioners are often motivated by access models in which the matrix is unavailable but operator products are cheap. In variational image reconstruction, randomized Nyström preconditioning is built entirely from operator calls to $h^{-1}(x,\xi)$0, $h^{-1}(x,\xi)$1, $h^{-1}(x,\xi)$2, and $h^{-1}(x,\xi)$3, using a stabilized Nyström pipeline $h^{-1}(x,\xi)$4, $h^{-1}(x,\xi)$5, Cholesky of $h^{-1}(x,\xi)$6, and an SVD of a skinny matrix $h^{-1}(x,\xi)$7. The paper emphasizes on-the-fly construction on GPUs and reports that moderate sketch ranks, such as $h^{-1}(x,\xi)$8–$h^{-1}(x,\xi)$9, were sufficient in its experiments (Hong et al., 2024).

A closely related matrix-free regime appears in strong-constraint 4D-Var. There, randomized preconditioners are assembled from TLM and adjoint products. RandSVD and Nyström require separate TLM and adjoint phases, while SingleView computes $H^{-1}\approx A(c):=\sum_{j=1}^{p} c_j B_j$0 and $H^{-1}\approx A(c):=\sum_{j=1}^{p} c_j B_j$1 concurrently; the paper reports peak parallel efficiency of about $H^{-1}\approx A(c):=\sum_{j=1}^{p} c_j B_j$2 for RandSVD/Nyström and about $H^{-1}\approx A(c):=\sum_{j=1}^{p} c_j B_j$3 for SingleView (Subrahmanya et al., 2024).

Randomization is also compatible with mixed precision. MRCQR constructs a right preconditioner $H^{-1}\approx A(c):=\sum_{j=1}^{p} c_j B_j$4 from a subsampled randomized trigonometric transform and then applies double-precision CholeskyQR to $H^{-1}\approx A(c):=\sum_{j=1}^{p} c_j B_j$5. Its perturbation analysis isolates the preconditioner from the final accuracy requirement and yields the practical thresholds $H^{-1}\approx A(c):=\sum_{j=1}^{p} c_j B_j$6 for an FP32 sketch and $H^{-1}\approx A(c):=\sum_{j=1}^{p} c_j B_j$7 for an FP16 sketch, while still achieving $H^{-1}\approx A(c):=\sum_{j=1}^{p} c_j B_j$8 in double precision for condition numbers up to $H^{-1}\approx A(c):=\sum_{j=1}^{p} c_j B_j$9 (Garrison et al., 16 Jun 2026). The mixed-precision normal-equations analysis of Garrison and Ipsen reaches a similar conclusion: the conditioning depends only mildly on the quality of the randomized preconditioner, and automatic precision selection can be driven by a fast condition estimator (Garrison et al., 17 Mar 2026).

When the target matrix is a grounded Laplacian, randomized factorization becomes a parallel sparse-linear-algebra problem. ParAC dynamically discovers dependencies created by sampled fill-ins, uses on-the-fly scheduling on CPUs and GPUs, and constructs $H$0 without the heavy symbolic preprocessing typical of classical incomplete factorizations (Liang et al., 5 May 2025).

5. Domain-specific realizations

In wave-based imaging, randomized preconditioning is strongly tied to microlocal structure. Matrix probing of the wave-equation Hessian treats the normal operator as pseudodifferential, uses curvelets to localize in position, scale, and direction, and fits a compressed inverse symbol from randomized illuminated probes. The resulting preconditioner corrects illumination and resolution in a dip-dependent manner and was proposed specifically for linearized seismic inversion (Demanet et al., 2011).

In data assimilation, the dominant use-case is the Gauss–Newton or inner-loop Hessian. For weak-constraint 4D-Var, randomized low-rank eigenspace approximations from the current Hessian are used to build limited-memory spectral preconditioners that are more robust than recycling Ritz information from previous inner loops (Daužickaitė et al., 2021). In strong-constraint 4D-Var, three sketching methods—RandSVD, Nyström, and SingleView—construct scalable preconditioners for the Gauss–Newton system and substantially reduce both PCG iterations and total TLM/adjoint counts (Subrahmanya et al., 2024). In ensemble variational data assimilation, the preconditioner is constructed from a low-rank approximation of the control Hessian using a sketching matrix built from differences of the right-hand sides across ensemble members, then reused across the entire ensemble (Daužickaitė et al., 22 May 2026).

Inverse problems in imaging supply another large class of examples. In $H$1–$H$2 iteratively reweighted reconstruction, randomized Nyström preconditioners target the weighted Hessian $H$3 inside PCG. In $H$4–$H$5 reconstruction, the same low-rank metric is integrated into weighted accelerated proximal-gradient schemes, with efficient weighted proximal maps derived from the diagonal-plus-low-rank structure (Hong et al., 2024).

Least-squares and kernel methods have produced several distinct randomized preconditioning models. Randomized preconditioned normal equations use DCT-based smoothing plus row sampling to construct an upper-triangular preconditioner $H$6, and both the symmetrically preconditioned and half-preconditioned normal equations are shown to be almost as accurate as Matlab’s QR-based mldivide when the preconditioner is effective (Ipsen, 24 Jul 2025). A separate line uses randomized preconditioner selection: among candidate block-diagonal or low-rank-plus-pinch preconditioners, Gaussian stability sketches can select a preconditioner that never uses more iterations than the unpreconditioned solver in the kernel-regression tests reported in the paper (DiPaolo et al., 2019). The block-Krylov approach of “preconditioning without a preconditioner” extends this spectrum-deflation viewpoint to ridge-regression paths and Gaussian sampling, with a single randomized Krylov basis reused across all shifts $H$7 (Chen et al., 30 Jan 2025).

The concept also extends beyond linear solves. In CCA and FDA on generalized Stiefel manifolds, randomized sketches are used to build Riemannian metrics $H$8, thereby reducing the local Hessian condition number of orthogonality-constrained optimization problems (Shustin et al., 2021). In tall-matrix nullspace computation, RLOBPCG uses a sketch-and-precondition construction $H$9 inside LOBPCG to compute smallest singular triplets of $A(c)H u_i \approx u_i$0 (Epperly et al., 18 Feb 2026).

6. Limitations, misconceptions, and active directions

Randomized preconditioning is not a universal shortcut. Its strongest guarantees generally require either effective low rank, rapid spectral decay, or a well-resolved visible subspace. Nyström PCG is explicitly limited when $A(c)H u_i \approx u_i$1, when $A(c)H u_i \approx u_i$2 is very small, or when the matrix lacks low effective rank (Frangella et al., 2021). In matrix probing, severe heterogeneity, anisotropy, or very limited aperture can shrink the visible set $A(c)H u_i \approx u_i$3, reduce the effective rank $A(c)H u_i \approx u_i$4, and make coefficient fitting ill-conditioned (Demanet et al., 2011).

A second misconception is that every randomized criterion directly predicts Krylov performance. The stability proxy $A(c)H u_i \approx u_i$5 correlates strongly with iterations in the experiments of Frangella, Tropp, and Udell, but the same paper explicitly states that Frobenius-norm closeness to the identity may not perfectly align with CG or GMRES ordering in pathological cases. It also evaluates candidate preconditioners rather than constructing them; randomized selection and randomized preconditioning are adjacent but not identical notions (DiPaolo et al., 2019).

A third distinction concerns where the randomness actually enters. In overlapping Schwarz preconditioners for randomized neural networks, the hidden layers are randomly initialized, but the additive and restricted additive Schwarz preconditioners are deterministic once the least-squares matrix $A(c)H u_i \approx u_i$6 is assembled. The paper states this explicitly: stochastic variability arises solely from the randomized bases that form $A(c)H u_i \approx u_i$7, while the Schwarz operators themselves are deterministic (Shang et al., 2024).

Current directions are correspondingly technical rather than merely heuristic. They include adaptive rank selection and error-based doubling for Nyström PCG (Frangella et al., 2021), adaptive bases and cross-validation strategies for inverse-Hessian probing (Demanet et al., 2011), multi-round filtering for candidate selection (DiPaolo et al., 2019), Q-less projector implementations and factorization-free range-deflation variants (Balabanov et al., 24 Sep 2025), and mixed-precision precision selection driven by fast condition estimation (Garrison et al., 17 Mar 2026, Garrison et al., 16 Jun 2026). This suggests that the modern meaning of “randomized preconditioner” is less a single algorithm than a design principle: exploit randomized access to informative spectral or phase-space structure, then convert that structure into a cheap operator that improves the linear algebra of the original problem.