Universal Least-Squares Solvers Overview

• Universal least-squares solvers are computational frameworks that solve diverse least-squares problems without redesigning the core algorithm for each instance.
• They employ unified approaches such as robust Gauss–Newton iterations, randomized sketching, and Kaczmarz-type methods to handle both consistent and inconsistent systems.
• These solvers enable practical applications in machine learning, statistical inference, and high-dimensional optimization through adaptable, precomputed or iterative strategies.

Universal least-squares solvers are methods, software systems, or linear operators designed to solve broad classes of least-squares problems without redesigning the core solver for each instance. In the cited literature, the term is used in more than one precise sense. In nonlinear estimation, it denotes a unified non-linear least-squares engine that transparently handles different factor-graph structures and manifold-valued variables (Grisetti et al., 2020). In randomized numerical linear algebra, it denotes algorithmic frameworks that unify single-row, block, randomized, and deterministic schemes for consistent and inconsistent least-squares problems (Shao, 2024). In operator-theoretic form, it denotes a matrix $H\in\mathbb R^{n\times m}$ such that $\hat\theta = H b$ solves $\min_{\theta\in\mathbb R^n}\|A\theta-b\|_2^2$ for every right-hand side $b\in\mathbb R^m$ (Machado et al., 4 Sep 2025).

1. Conceptual scope and mathematical formulations

Least-squares universality begins with a small number of recurring formulations. For nonlinear estimation, the canonical non-linear least-squares objective is

$F(x)=\tfrac12\sum_i r_i(x)^\top W_i r_i(x),$

or, after stacking residuals,

$F(x)=\tfrac12\,r(x)^\top W\,r(x).$

Linearization $r(x+\delta x)\approx r(x)+J(x)\delta x$ yields the Gauss–Newton normal equations $(J^\top WJ)\delta x=-J^\top Wr$, and robustification is implemented through IRLS with $\gamma_i=\rho'(u_i)/u_i$ and $W_i\leftarrow \gamma_i W_i$ (Grisetti et al., 2020).

In machine learning, the least-mean-squares template is

$\hat\theta = H b$0

with $\hat\theta = H b$1 a convex regularizer such as $\hat\theta = H b$2 or $\hat\theta = H b$3. This single form subsumes linear regression, Ridge, Lasso, Elastic-Net, SVD, and PCA in the cited treatment (Maalouf et al., 2019).

For rectangular or rank-deficient linear systems, several papers adopt the minimum-length least-squares solution $\hat\theta = H b$4, where $\hat\theta = H b$5 is the Moore–Penrose pseudoinverse. This is the unique minimum-length minimizer of $\hat\theta = H b$6, and it remains the target even when $\hat\theta = H b$7 is rank-deficient or when the system is strongly over- or under-determined (Meng et al., 2011).

A stricter meaning of universality fixes the matrix $\hat\theta = H b$8 and seeks a single solver matrix $\hat\theta = H b$9 valid for every right-hand side. In that setting, $\min_{\theta\in\mathbb R^n}\|A\theta-b\|_2^2$0 is a universal least-squares solver for $\min_{\theta\in\mathbb R^n}\|A\theta-b\|_2^2$1 if $\min_{\theta\in\mathbb R^n}\|A\theta-b\|_2^2$2 attains the minimum of $\min_{\theta\in\mathbb R^n}\|A\theta-b\|_2^2$3 for every $\min_{\theta\in\mathbb R^n}\|A\theta-b\|_2^2$4, equivalently if $\min_{\theta\in\mathbb R^n}\|A\theta-b\|_2^2$5 satisfies

$\min_{\theta\in\mathbb R^n}\|A\theta-b\|_2^2$6

These are the Moore–Penrose properties $\min_{\theta\in\mathbb R^n}\|A\theta-b\|_2^2$7 and $\min_{\theta\in\mathbb R^n}\|A\theta-b\|_2^2$8 in the cited formulation (Machado et al., 4 Sep 2025).

Meaning of universality	Representative object	Representative sources
Cross-domain optimization engine	FactorGraph + NLLS solver with interchangeable backends	(Grisetti et al., 2020)
Cross-regime iterative framework	Randomized, block, reflective, or preconditioned least-squares methods	(Shao, 2024, Zouzias et al., 2012, Meng et al., 2011, Yang et al., 2024, Sugihara et al., 14 Apr 2025)
Right-hand-side independent linear solver	Matrix $\min_{\theta\in\mathbb R^n}\|A\theta-b\|_2^2$9 such that $b\in\mathbb R^m$0 solves LS for all $b\in\mathbb R^m$1	(Machado et al., 4 Sep 2025)

This multiplicity of meanings is important. The literature does not present a single canonical universal least-squares solver; instead, it presents several distinct notions of universality tied to software architecture, algorithmic family, or generalized-inverse structure.

2. Unified nonlinear engines and probabilistic extensions

A software-centric universal solver is exemplified by the factor-graph system in “Least Squares Optimization: from Theory to Practice” (Grisetti et al., 2020). Its driver linearizes every factor, applies robust weighting if requested, assembles the global block-sparse $b\in\mathbb R^m$2 and $b\in\mathbb R^m$3, solves $b\in\mathbb R^m$4, updates $b\in\mathbb R^m$5 or $b\in\mathbb R^m$6 on manifolds, and checks convergence through cost decrease, step-size, and maximum-iteration conditions. The Hessian is assembled “per factor” in a double loop over the variables touched by that factor, accumulating only the lower-triangular blocks. Because each factor depends only on a small subset of variables, the resulting matrix is symmetric and block-sparse. To reduce fill-in, the system can apply reorderings such as AMD and COLAMD before sparse Cholesky or QR factorization (Grisetti et al., 2020).

The same system couples algorithmic generality with implementation modularity. It stores the sparse-block Hessian in a custom data structure containing only the nonzero blocks $b\in\mathbb R^m$7, supports Cholesky backends through SuiteSparse and Eigen’s SimplicialLLT, QR backends through Householder and TSQR, and iterative PCG for very large-scale dense problems such as single-pose ICP. All backends implement a common C++ interface, ISolver, and can be swapped at runtime via configuration. The core library is reported as $b\in\mathbb R^m$8 kLoC of C++17 and relies on BOSS for runtime loading and serialization, a heap-allocated FactorGraph exposing Variables and Factors, templated Variable<Dim,Type> and Factor<...> abstractions, and optional automatic differentiation via templated scalar types or analytic Jacobians (Grisetti et al., 2020).

Its extensibility model is explicitly stratified. New variables are introduced by implementing $b\in\mathbb R^m$9 and, when needed, $F(x)=\tfrac12\sum_i r_i(x)^\top W_i r_i(x),$0. New factors can be supplied as FactorBase for direct update(H,b) control, as ErrorFactor<ErrDim,Var...> with explicit error() and Jacobians(), or as ADErrorFactor implementing only error(). A correspondence-driven factor can process many data-association entries within one factor object. Robustifiers are injected as policies and can be assigned to subsets of factors, such as Huber on loop closures and Cauchy on visual odometry. The same outer-loop interface is meant to accommodate Dogleg and trust-region methods because they share the build-linear-system / solve-linear-step pattern (Grisetti et al., 2020).

The same universal engine perspective extends beyond Gaussian residual models. “Advancing Mixture Models for Least Squares Optimization” constructs an exact least-squares representation of a Gaussian mixture by extracting the dominant component $F(x)=\tfrac12\sum_i r_i(x)^\top W_i r_i(x),$1 and defining a vector residual $F(x)=\tfrac12\sum_i r_i(x)^\top W_i r_i(x),$2 of dimension $F(x)=\tfrac12\sum_i r_i(x)^\top W_i r_i(x),$3. Its first $F(x)=\tfrac12\sum_i r_i(x)^\top W_i r_i(x),$4 entries are the Mahalanobis residual of the dominant Gaussian, while the final scalar entry encodes the full log-sum-exp mixture term. The method can be inserted directly into Ceres through an AutoDiffCostFunction or into GTSAM through a custom NoiseModelFactor, with no custom loss and no changes to the core solver (Pfeifer et al., 2021).

This construction addresses a common misconception that Gaussian mixtures must be approximated before they can be used in least-squares solvers. In the cited formulation, the model is exact at the log-likelihood level and retains an almost-linear local structure through the dominant mode. The reported experiments include 100 000 Monte Carlo runs in which MSM achieved 100% success across all 1D/2D settings considered, while Max-Mixture and Sum-Mixture exhibited failure modes on asymmetric or higher-dimensional cases. In point-set registration, MSM is reported to deliver the best trade-off among exact GMM modeling, robust convergence, small RMSE, credible covariances, and modest overhead over MM (Pfeifer et al., 2021).

3. Compression, sketching, and solver wrapping

A second route to universality compresses the data rather than generalizing the optimizer itself. “Fast and Accurate Least-Mean-Squares Solvers” constructs a weighted subset of rows of $F(x)=\tfrac12\sum_i r_i(x)^\top W_i r_i(x),$5 whose covariance is exactly the same as $F(x)=\tfrac12\sum_i r_i(x)^\top W_i r_i(x),$6, so that solving the least-squares problem on the compressed instance recovers the exact solution up to machine precision (Maalouf et al., 2019). The classical Carathéodory construction yields a subset of size $F(x)=\tfrac12\sum_i r_i(x)^\top W_i r_i(x),$7, but its $F(x)=\tfrac12\sum_i r_i(x)^\top W_i r_i(x),$8 runtime is replaced by a “booster” with complexity $F(x)=\tfrac12\sum_i r_i(x)^\top W_i r_i(x),$9. For large $F(x)=\tfrac12\,r(x)^\top W\,r(x).$0, a sparsified variant runs in $F(x)=\tfrac12\,r(x)^\top W\,r(x).$1 time and returns $F(x)=\tfrac12\,r(x)^\top W\,r(x).$2 sparsified input points (Maalouf et al., 2019).

The universal aspect appears in the plug-in interface. After constructing a coreset $F(x)=\tfrac12\,r(x)^\top W\,r(x).$3, one replaces

$F(x)=\tfrac12\,r(x)^\top W\,r(x).$4

The paper states that the solution is identical up to numerical precision and presents wrappers for LinearRegression, RidgeCV, LassoCV, and ElasticNetCV. It also states that SVD/PCA admits a coreset of size $F(x)=\tfrac12\,r(x)^\top W\,r(x).$5 or $F(x)=\tfrac12\,r(x)^\top W\,r(x).$6 with exactly zero subspace-approximation error. Streaming and distributed variants are obtained by merge-and-reduce, using $F(x)=\tfrac12\,r(x)^\top W\,r(x).$7 memory and $F(x)=\tfrac12\,r(x)^\top W\,r(x).$8 amortized time per point in the streaming case (Maalouf et al., 2019).

A related but distinct form of sketch-based universality is preconditioning by random normal projection. LSRN treats

$F(x)=\tfrac12\,r(x)^\top W\,r(x).$9

when $r(x+\delta x)\approx r(x)+J(x)\delta x$0, $r(x+\delta x)\approx r(x)+J(x)\delta x$1, or $r(x+\delta x)\approx r(x)+J(x)\delta x$2 is rank-deficient, and it permits Tikhonov regularization by rewriting the problem as an unregularized least-squares system. In the over-determined case it chooses $r(x+\delta x)\approx r(x)+J(x)\delta x$3, draws $r(x+\delta x)\approx r(x)+J(x)\delta x$4 with i.i.d. $r(x+\delta x)\approx r(x)+J(x)\delta x$5 entries, forms $r(x+\delta x)\approx r(x)+J(x)\delta x$6, computes its SVD, and builds a right preconditioner $r(x+\delta x)\approx r(x)+J(x)\delta x$7. The resulting spectrum of $r(x+\delta x)\approx r(x)+J(x)\delta x$8 depends on a Gaussian matrix rather than on the original conditioning of $r(x+\delta x)\approx r(x)+J(x)\delta x$9, and if $(J^\top WJ)\delta x=-J^\top Wr$0 then $(J^\top WJ)\delta x=-J^\top Wr$1 with overwhelming probability (Meng et al., 2011).

This makes the iteration count predictable for LSQR or for the Chebyshev semi-iterative method. The latter is emphasized as particularly efficient on clusters with high communication cost because it avoids the per-iteration global reductions of LSQR. The same architecture permits $(J^\top WJ)\delta x=-J^\top Wr$2 to be dense, sparse, or a linear operator, since only matrix-matrix and matrix-vector multiplications are required (Meng et al., 2011).

A common misconception is that exactness and compression are incompatible. The coreset-based LMS framework explicitly states exact covariance preservation and exact LMS outputs, while LSRN shows that randomized sketches can be used not as approximate solvers but as preconditioners for high-precision iterative refinement (Maalouf et al., 2019, Meng et al., 2011).

4. Projection, reflection, and Kaczmarz-type universality

Projection methods form another universal family because they naturally admit consistent, inconsistent, block, and randomized variants. The classical Kaczmarz method solves $(J^\top WJ)\delta x=-J^\top Wr$3 by projecting onto hyperplanes $(J^\top WJ)\delta x=-J^\top Wr$4, and randomized Kaczmarz chooses row $(J^\top WJ)\delta x=-J^\top Wr$5 with probability $(J^\top WJ)\delta x=-J^\top Wr$6, yielding the expected linear convergence

$(J^\top WJ)\delta x=-J^\top Wr$7

Block Kaczmarz replaces one hyperplane by the intersection associated with a block $(J^\top WJ)\delta x=-J^\top Wr$8, updating $(J^\top WJ)\delta x=-J^\top Wr$9 (Shao, 2024).

Reflective block Kaczmarz replaces projection by reflection. For a single row, the reflection operator is

$\gamma_i=\rho'(u_i)/u_i$0

and the iterate is updated through a Householder reflection. In the consistent case, the iterates lie on a high-dimensional sphere centered at the solution. In the deterministic cyclic version, the even and odd subsequences lie on two spheres whose centers coincide exactly at the true solution. Block reflections generalize this geometry to intersections of multiple hyperplanes via

$\gamma_i=\rho'(u_i)/u_i$1

The averaged iterate $\gamma_i=\rho'(u_i)/u_i$2 is the returned estimate (Shao, 2024).

For least-squares problems, the reflective framework analyzes the inconsistent decomposition $\gamma_i=\rho'(u_i)/u_i$3 with $\gamma_i=\rho'(u_i)/u_i$4, where $\gamma_i=\rho'(u_i)/u_i$5. The cited bounds show that the averaged reflective iterate converges with an $\gamma_i=\rho'(u_i)/u_i$6 term plus an error floor depending on the inconsistency. In the block-consistent setting, the contraction improves through a block-conditioning parameter $\gamma_i=\rho'(u_i)/u_i$7. Numerical experiments on $\gamma_i=\rho'(u_i)/u_i$8 Gaussian matrices report that reflective methods often reach a prescribed $\gamma_i=\rho'(u_i)/u_i$9 error in 30–50% fewer steps than standard Kaczmarz, with an optimal runtime trade-off near block size $W_i\leftarrow \gamma_i W_i$0 for one tested configuration (Shao, 2024).

Randomized Extended Kaczmarz addresses the minimum-Euclidean-norm least-squares problem directly. It interleaves row and column sampling: one update denoises $W_i\leftarrow \gamma_i W_i$1 by projecting an auxiliary vector $W_i\leftarrow \gamma_i W_i$2 onto $W_i\leftarrow \gamma_i W_i$3, and the other performs a Kaczmarz step on the corrected system $W_i\leftarrow \gamma_i W_i$4. In exact arithmetic,

$W_i\leftarrow \gamma_i W_i$5

and the total arithmetic cost to reach error $W_i\leftarrow \gamma_i W_i$6 is

$W_i\leftarrow \gamma_i W_i$7

The method is described as sparse, streaming-friendly, and naturally complementary to direct factorizations, Krylov methods, and sketch-and-solve approaches (Zouzias et al., 2012).

These results directly counter the assumption that Kaczmarz-type methods are restricted to consistent systems or to feasibility formulations. In the cited literature, both reflective Kaczmarz and REK are formulated as least-squares solvers, including explicit treatment of inconsistency, minimum-norm solutions, and block extensions (Shao, 2024, Zouzias et al., 2012).

5. Krylov, flexible preconditioning, and inverse approximation

Krylov methods supply universality through matrix structure and preconditioning rather than through data compression. For symmetric or Hermitian problems, MINRES-QLP solves $W_i\leftarrow \gamma_i W_i$8 and, when the system is singular, computes the unique minimum-length solution $W_i\leftarrow \gamma_i W_i$9. Its defining feature is a QLP factorization of the Lanczos tridiagonal subproblem: standard MINRES reduces the projected system to an upper-triangular $\hat\theta = H b$00, while MINRES-QLP applies additional right reflectors to obtain a lower-triangular $\hat\theta = H b$01, whose diagonal entries reveal numerical rank more faithfully. The method uses one $\hat\theta = H b$02-product and one preconditioner solve $\hat\theta = H b$03 per iteration, preserves short recurrences, and is intended precisely for symmetric or Hermitian indefinite, singular, and least-squares settings (Choi et al., 2013).

For general sparse least-squares problems, Flexible Modified LSMR replaces the two inner linear solves of right-preconditioned LSMR by a single solve

$\hat\theta = H b$04

and then imports the flexible-GMRES idea by allowing this solve to be performed by a different effective preconditioner $\hat\theta = H b$05 at each outer step. The stopping test remains the norm-wise backward-error criterion

$\hat\theta = H b$06

Across eight large sparse test matrices, FMLSMR converged on all eight problems to the $\hat\theta = H b$07 normal-residual tolerance, whereas standard LSMR failed on three and Flexible LSMR failed on five. On the $\hat\theta = H b$08 “well1850” problem, the reported outer iteration counts are 463 for LSMR, 167 for FLSMR, and 117 for FMLSMR with 8 inner CG steps. On “delaunay_n16” $\hat\theta = H b$09, only FMLSMR is reported to finish, in approximately 1306 s (Yang et al., 2024).

A further universal construction treats arbitrary singular and rectangular least-squares problems through RRGMRES applied to

$\hat\theta = H b$10

with recovery $\hat\theta = H b$11. Because $\hat\theta = H b$12 is symmetric and therefore range-symmetric, RRGMRES on this transformed system is stated to be breakdown-free for arbitrary $\hat\theta = H b$13, $\hat\theta = H b$14, and $\hat\theta = H b$15, with convergence in at most $\hat\theta = H b$16 steps. The NR-SSOR inner-iteration right preconditioner generates $\hat\theta = H b$17, and the nonzero eigenvalues of $\hat\theta = H b$18 are clustered in $\hat\theta = H b$19. In the reported experiments, NR-SSOR-preconditioned RRGMRES is 3–6× faster than the unpreconditioned variant and attains $\hat\theta = H b$20 values around $\hat\theta = H b$21 in the tested settings (Sugihara et al., 14 Apr 2025).

Not all universal solvers are Krylov methods. The MinCos line develops a geometrical inverse approximation for least-squares by minimizing $\hat\theta = H b$22, with $\hat\theta = H b$23 in the full-column-rank case. The simplified gradient-type scheme MinCos produces a convergent sequence $\hat\theta = H b$24, and the paper studies three accelerations: STEA2 matrix extrapolation, ABBmin deterministic step modification, and random relaxation. The reported guidance is explicit: use STEA2 for moderate sizes, Random MinCos for very large or ill-conditioned systems, and LSQR or CG only when specialized Krylov structure is required (Chehab et al., 2019).

Taken together, these methods show that universality in Krylov and inverse-approximation settings is not tied to one recurrence. It may mean short-recurrence robustness under singularity, flexible inner solves, communication avoidance, spectral clustering through preconditioning, or direct approximation of $\hat\theta = H b$25.

6. Sparse universal solver matrices, structural optimization, and interpretation

The operator-theoretic definition of a universal least-squares solver is sharpened in “On computing sparse universal solvers for key problems in statistics” (Machado et al., 4 Sep 2025). For fixed $\hat\theta = H b$26, the set of all universal LS solvers is

$\hat\theta = H b$27

The paper then asks for solvers that are sparse, low-rank, or simultaneously least-squares and minimum-norm.

Three optimization problems are distinguished. The minimum-$\hat\theta = H b$28 universal LS solver minimizes $\hat\theta = H b$29 under $\hat\theta = H b$30 and $\hat\theta = H b$31. The minimum-rank universal LS solver adds reflexivity $\hat\theta = H b$32, and the simultaneous universal LS-and-minimum-norm solver adds $\hat\theta = H b$33. The corresponding sparsity bounds for extreme solutions of the LP reformulations are stated as at most $\hat\theta = H b$34 nonzeros for $\hat\theta = H b$35, at most $\hat\theta = H b$36 nonzeros for $\hat\theta = H b$37, and at most $\hat\theta = H b$38 nonzeros for $\hat\theta = H b$39, where $\hat\theta = H b$40 (Machado et al., 4 Sep 2025).

The computational framework is first-order rather than factorization-based. The $\hat\theta = H b$41-objective is split from the affine Moore–Penrose constraint set, the proximal map of $\hat\theta = H b$42 is soft-thresholding, and the projection onto the constraint set $\hat\theta = H b$43 is available in closed form through SVD or pseudoinverse manipulations. The Douglas–Rachford iteration alternates between $\hat\theta = H b$44, reflection, projection $\hat\theta = H b$45, and relaxation. In the reported experiments on random dense matrices with sizes up to $\hat\theta = H b$46, DRS$\hat\theta = H b$47 is the fastest method and scales easily to $\hat\theta = H b$48, whereas Gurobi cannot solve beyond $\hat\theta = H b$49 within 2 h. DRS$\hat\theta = H b$50 and ADMM achieve virtually identical $\hat\theta = H b$51-values, but DRS is 5–10× faster. Relative to the Moore–Penrose inverse $\hat\theta = H b$52, the $\hat\theta = H b$53-solutions reduce $\hat\theta = H b$54 by up to 90% and $\hat\theta = H b$55 by up to 60% (Machado et al., 4 Sep 2025).

This operator view clarifies a final ambiguity in the subject. A universal least-squares solver need not be an iterative procedure that acts separately on each right-hand side; it may be a precomputed matrix that encodes the least-squares map for all right-hand sides at once. Conversely, the software and algorithmic literatures use “universal” to mean broad applicability across structures, residual models, or matrix regimes rather than a single closed-form operator. The cited record therefore supports a plural interpretation: universal least-squares solving is a family of design goals—cross-domain modeling, cross-regime iteration, and right-hand-side independent solution maps—rather than a single algorithmic doctrine (Machado et al., 4 Sep 2025).