---
title: Universal Least-Squares Solvers Overview
url: https://www.emergentmind.com/topics/universal-least-squares-solvers
type: topic
---

# Universal Least-Squares Solvers Overview

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 [2002.11051]. 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 [2407.19226]. 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\) [2509.04264].

## 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\) [2002.11051].

In machine learning, the least-mean-squares template is
\[
\min_{x\in X}\;\|Ax-b\|_2^2 + g(x),
\]
with \(g(x)\) a convex regularizer such as \(\alpha\|x\|_2^2\) or \(\alpha\|x\|_1\). This single form subsumes linear regression, Ridge, Lasso, Elastic-Net, SVD, and PCA in the cited treatment [1906.04705].

For rectangular or rank-deficient linear systems, several papers adopt the minimum-length least-squares solution \(x^*=A^\dagger b\), where \(A^\dagger\) is the Moore–Penrose pseudoinverse. This is the unique minimum-length minimizer of \(\min_x\|Ax-b\|_2\), and it remains the target even when \(A\) is rank-deficient or when the system is strongly over- or under-determined [1109.5981].

A stricter meaning of universality fixes the matrix \(A\) and seeks a single solver matrix \(H\) valid for every right-hand side. In that setting, \(H\) is a universal least-squares solver for \(A\) if \(\hat\theta=Hb\) attains the minimum of \(\min_\theta \|A\theta-b\|_2^2\) for every \(b\), equivalently if \(H\) satisfies
\[
AHA=A,\qquad (AH)^\top=AH.
\]
These are the Moore–Penrose properties \(\mathrm{(P1)}\) and \(\mathrm{(P3)}\) in the cited formulation [2509.04264].

| Meaning of universality | Representative object | Representative sources |
|---|---|---|
| Cross-domain optimization engine | FactorGraph + NLLS solver with interchangeable backends | [2002.11051] |
| Cross-regime iterative framework | Randomized, block, reflective, or preconditioned least-squares methods | [2407.19226], [1205.5770], [1109.5981], [2408.16652], [2504.09891] |
| Right-hand-side independent linear solver | Matrix \(H\) such that \(Hb\) solves LS for all \(b\) | [2509.04264] |

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” [2002.11051]. Its driver linearizes every factor, applies robust weighting if requested, assembles the global block-sparse \(H\) and \(b\), solves \((H+\Lambda)\delta x=-b\), updates \(x\leftarrow x\oplus\delta x\) or \(x\leftarrow x\boxplus\delta x\) 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 [2002.11051].

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 \(H_{ij}\), 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 \(<6\) 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 [2002.11051].

Its extensibility model is explicitly stratified. New variables are introduced by implementing \(\boxplus\) and, when needed, \(\boxminus\). 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 [2002.11051].

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 \(k=\arg\max_\ell [s_\ell\exp(e_\ell(r))]\) and defining a vector residual \(R_{\rm MSM}(r)\) of dimension \(\dim(r)+1\). Its first \(\dim(r)\) 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 [2103.02472].

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

## 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 \(A\) whose covariance is exactly the same as \(A^\top A\), so that solving the least-squares problem on the compressed instance recovers the exact solution up to machine precision [1906.04705]. The classical Carathéodory construction yields a subset of size \(\le d+1\), but its \(O(n^2d^2)\) runtime is replaced by a “booster” with complexity \(O(nd+d^4\log n)\). For large \(d\), a sparsified variant runs in \(O(nd)\) time and returns \(O(d)\) sparsified input points [1906.04705].

The universal aspect appears in the plug-in interface. After constructing a coreset \((C,y)\), one replaces
\[
\min_x\|Ax-b\|_2^2+g(x)
\quad\Longrightarrow\quad
\min_x\|Cx-y\|_2^2+g(x).
\]
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 \((d+1)^2+1\) or \(O(d)\) with exactly zero subspace-approximation error. Streaming and distributed variants are obtained by merge-and-reduce, using \(O((d^2+1)\log n)\) memory and \(O(d)\) amortized time per point in the streaming case [1906.04705].

A related but distinct form of sketch-based universality is preconditioning by random normal projection. LSRN treats
\[
\min_{x\in\mathbb R^n}\|Ax-b\|_2
\]
when \(m\gg n\), \(m\ll n\), or \(A\) 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 \(s=\lceil\gamma n\rceil\), draws \(G\in\mathbb R^{s\times m}\) with i.i.d. \(N(0,1)\) entries, forms \(\widetilde A=GA\), computes its SVD, and builds a right preconditioner \(N=\widetilde V\Sigma^{-1}\). The resulting spectrum of \(AN\) depends on a Gaussian matrix rather than on the original conditioning of \(A\), and if \(s=2r\) then \(\kappa(AN)<6\) with overwhelming probability [1109.5981].

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 \(A\) to be dense, sparse, or a linear operator, since only matrix-matrix and matrix-vector multiplications are required [1109.5981].

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 [1906.04705] [1109.5981].

## 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 \(Ax=b\) by projecting onto hyperplanes \(A_i^\top x=b_i\), and randomized Kaczmarz chooses row \(i_k\) with probability \(\|A_{i_k}\|^2/\|A\|_F^2\), yielding the expected linear convergence
\[
\mathbb E\big[\|x^k-x^*\|^2\big]\le
\Bigl(1-\frac{\sigma_{\min}^2(A)}{\|A\|_F^2}\Bigr)^k\|x^0-x^*\|^2.
\]
Block Kaczmarz replaces one hyperplane by the intersection associated with a block \(Z\), updating \(x^{k+1}=x^k-A_Z^+(A_Zx^k-b_Z)\) [2407.19226].

Reflective block Kaczmarz replaces projection by reflection. For a single row, the reflection operator is
\[
R_i=I-2\,\frac{A_iA_i^\top}{\|A_i\|^2},
\]
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
\[
R_Z=I_n-2A_Z^\top(A_ZA_Z^\top)^{-1}A_Z.
\]
The averaged iterate \(\bar x^N=\tfrac1N\sum_{k=0}^{N-1}x^k\) is the returned estimate [2407.19226].

For least-squares problems, the reflective framework analyzes the inconsistent decomposition \(b=Ax^*+c\) with \(A^\top c=0\), where \(x^*=A^+b\). The cited bounds show that the averaged reflective iterate converges with an \(O(1/N)\) term plus an error floor depending on the inconsistency. In the block-consistent setting, the contraction improves through a block-conditioning parameter \(\gamma\). Numerical experiments on \(300\times100\) Gaussian matrices report that reflective methods often reach a prescribed \(10^{-3}\) error in 30–50% fewer steps than standard Kaczmarz, with an optimal runtime trade-off near block size \(q\approx10\) for one tested configuration [2407.19226].

Randomized Extended Kaczmarz addresses the minimum-Euclidean-norm least-squares problem directly. It interleaves row and column sampling: one update denoises \(b\) by projecting an auxiliary vector \(z^{(k)}\) onto \(\mathrm{null}(a_{j_k}^\top)\), and the other performs a Kaczmarz step on the corrected system \(Ax=b-z^{(k)}\). In exact arithmetic,
\[
E\|x^{(T)}-x^*\|_2^2
\le
\Bigl(1-\frac1{\kappa^2(A)}\Bigr)^{\lfloor T/2\rfloor}
\Bigl(1+2\|A\|_F^2/\sigma_{\min}^2\Bigr)\|b\|_2^2,
\]
and the total arithmetic cost to reach error \(\varepsilon\) is
\[
O\bigl(\mathrm{nnz}(A)\,\kappa^2(A)\,\ln(1/\varepsilon)\bigr).
\]
The method is described as sparse, streaming-friendly, and naturally complementary to direct factorizations, Krylov methods, and sketch-and-solve approaches [1205.5770].

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 [2407.19226] [1205.5770].

## 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 \(\min\|Ax-b\|_2\) and, when the system is singular, computes the unique minimum-length solution \(A^+b\). Its defining feature is a QLP factorization of the Lanczos tridiagonal subproblem: standard MINRES reduces the projected system to an upper-triangular \(R_k\), while MINRES-QLP applies additional right reflectors to obtain a lower-triangular \(L_k\), whose diagonal entries reveal numerical rank more faithfully. The method uses one \(A\)-product and one preconditioner solve \(Mq=z\) per iteration, preserves short recurrences, and is intended precisely for symmetric or Hermitian indefinite, singular, and least-squares settings [1301.2707].

For general sparse least-squares problems, Flexible Modified LSMR replaces the two inner linear solves of right-preconditioned LSMR by a single solve
\[
M\,\tilde v_{k+1}=\tilde p_{k+1},
\]
and then imports the flexible-GMRES idea by allowing this solve to be performed by a different effective preconditioner \(M_k^{-1}\) at each outer step. The stopping test remains the norm-wise backward-error criterion
\[
\|A^T(b-Ax_k)\|_2 \le \epsilon \|A\|_2(\|b\|_2+\|A\|_2\|x_k\|_2).
\]
Across eight large sparse test matrices, FMLSMR converged on all eight problems to the \(10^{-12}\) normal-residual tolerance, whereas standard LSMR failed on three and Flexible LSMR failed on five. On the \(1850\times712\) “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” \((65\,536\times65\,536)\), only FMLSMR is reported to finish, in approximately 1306 s [2408.16652].

A further universal construction treats arbitrary singular and rectangular least-squares problems through RRGMRES applied to
\[
A\,C\,A^\top z=b,\qquad C\succ0,
\]
with recovery \(x=C A^\top z\). Because \(A C A^\top\) is symmetric and therefore range-symmetric, RRGMRES on this transformed system is stated to be breakdown-free for arbitrary \(A\), \(b\), and \(x_0\), with convergence in at most \(\mathrm{rank}(A)\) steps. The NR-SSOR inner-iteration right preconditioner generates \(C^{(\ell)}\succ0\), and the nonzero eigenvalues of \(A C^{(\ell)}A^\top\) are clustered in \([1-\rho(H)^\ell,\,1+\rho(H)^\ell]\cup\{0\}\). In the reported experiments, NR-SSOR-preconditioned RRGMRES is 3–6× faster than the unpreconditioned variant and attains \(\|A^T r_k\|\) values around \(10^{-14\text{–}-8}\) in the tested settings [2504.09891].

Not all universal solvers are Krylov methods. The MinCos line develops a geometrical inverse approximation for least-squares by minimizing \(1-\cos(X(A^\top A),I)\), with \(M=A^\top A\succ0\) in the full-column-rank case. The simplified gradient-type scheme MinCos produces a convergent sequence \(X^k\to (A^\top A)^{-1}\), 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 [1902.08388].

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 \((A^\top A)^{-1}\).

## 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” [2509.04264]. For fixed \(A\in\mathbb R^{m\times n}\), the set of all universal LS solvers is
\[
\{H\in\mathbb R^{n\times m}: AHA=A,\ (AH)^\top=AH\}.
\]
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-\(\ell_1\) universal LS solver minimizes \(\|H\|_1\) under \(\mathrm{(P1)}\) and \(\mathrm{(P3)}\). The minimum-rank universal LS solver adds reflexivity \(\mathrm{(P2)}: HAH=H\), and the simultaneous universal LS-and-minimum-norm solver adds \(\mathrm{(P4)}:(HA)^\top=HA\). The corresponding sparsity bounds for extreme solutions of the LP reformulations are stated as at most \(mr\) nonzeros for \((P_{13}^1)\), at most \(mr+(m-r)(n-r)\) nonzeros for \((P_{123}^1)\), and at most \(mn-(m-r)(n-r)\) nonzeros for \((P_{134}^1)\), where \(r=\mathrm{rank}(A)\) [2509.04264].

The computational framework is first-order rather than factorization-based. The \(\ell_1\)-objective is split from the affine Moore–Penrose constraint set, the proximal map of \(\|H\|_1\) is soft-thresholding, and the projection onto the constraint set \(\mathcal C\) is available in closed form through SVD or pseudoinverse manipulations. The Douglas–Rachford iteration alternates between \(S_\lambda(V^k)\), reflection, projection \(\Pi_{\mathcal C}\), and relaxation. In the reported experiments on random dense matrices with sizes up to \(m=5000\), DRS\(_{\rm fp}\) is the fastest method and scales easily to \(m=5000\), whereas Gurobi cannot solve beyond \(m\approx200\) within 2 h. DRS\(_{\rm fp}\) and ADMM achieve virtually identical \(\ell_1\)-values, but DRS is 5–10× faster. Relative to the Moore–Penrose inverse \(A^\dagger\), the \(\ell_1\)-solutions reduce \(\|H\|_1\) by up to 90% and \(\|H\|_0\) by up to 60% [2509.04264].

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

Source: https://www.emergentmind.com/topics/universal-least-squares-solvers