---
title: Universal Min-Norm/Least-Squares Solvers
url: https://www.emergentmind.com/topics/simultaneous-universal-minimum-norm-least-squares-solvers
type: topic
---

# Universal Min-Norm/Least-Squares Solvers

Searching arXiv for relevant papers on simultaneous universal minimum-norm/least-squares solvers and closely related universal solvers.
Simultaneous universal minimum-norm/least-squares solvers are operators \(H\in\mathbb{R}^{n\times m}\) associated with a fixed matrix \(A\in\mathbb{R}^{m\times n}\) such that, for every right-hand side \(b\), the vector \(Hb\) solves the ordinary least-squares problem \(\min_\theta \|A\theta-b\|_2^2\), and, when \(b\in \mathrm{Range}(A)\), \(Hb\) is also the minimum-\(2\)-norm solution of \(A\theta=b\). In the formulation studied in "On computing sparse universal solvers for key problems in statistics" [2509.04264], this is the simultaneous universal solver problem, or “P134 problem,” and it is characterized by three Penrose conditions rather than the full Moore–Penrose system. The same universal objective also appears in iterative Krylov, row-action, and factorization methods that aim to return either a least-squares minimizer or the pseudoinverse solution without case distinctions across consistent, inconsistent, singular, or rank-deficient regimes [2509.04264].

## 1. Algebraic definition and Penrose characterization

For a matrix \(A\) of rank \(r\), a simultaneous universal solver is a matrix \(H\) such that \(Hb\) is, for every \(b\), both a least-squares solution and, on \(\mathrm{Range}(A)\), the minimum-\(2\)-norm exact solution. The paper on sparse universal solvers states that this is equivalent to the three Penrose conditions

\[
AHA = A,\qquad (AH)^\top = AH,\qquad (HA)^\top = HA,
\]

denoted \(P1\), \(P3\), and \(P4\) [2509.04264].

With the full SVD
\[
A = U\Sigma V^\top = [U_1\ U_2]
\begin{bmatrix}
D & 0\\
0 & 0
\end{bmatrix}
[V_1\ V_2]^\top,
\]
where \(D\in\mathbb{R}^{r\times r}\) is diagonal with positive singular values, the same paper gives the parameterization
\[
H = V_1D^{-1}U_1^\top + V_2WU_2^\top,
\]
with \(W\in\mathbb{R}^{(n-r)\times(m-r)}\). This exhibits the entire affine family of simultaneous universal solvers. The Moore–Penrose inverse is one member of this family, but not the only one; the free block \(W\) spans the null-space degrees of freedom [2509.04264].

A key reduced characterization is the equivalence
\[
P1+P3+P4 \iff AA^\top H^\top + H^\top A^\top A = 2A^\top.
\]
This converts the simultaneous universal solver constraints into a single linear matrix equation. In the paper’s terminology, this is the “PMX” form, and it is central both for sparsity analysis and for proximal algorithms [2509.04264].

The surrounding literature clarifies why the distinction between least-squares optimality and minimum norm is necessary. For singular symmetric systems, standard MINRES may return a least-squares solution that is not the minimum-length solution, while the pseudoinverse solution requires extra structure or an explicit refinement [1003.4042]. In that sense, simultaneous universal solvers encode, at the operator level, the same separation that iterative methods confront at the algorithmic level.

## 2. Equivalent optimization formulations and sparsity structure

The sparse-solver formulation selects, among all simultaneous universal solvers, one minimizing the entrywise \(1\)-norm \(\|H\|_1\). The paper lists six equivalent linear-constraint formulations, all with objective \(\min \|H\|_1\) [2509.04264].

| Formulation | Constraints | Comment |
|---|---|---|
| \(P_{134}^1\) | \(AHA=A,\ (AH)^\top=AH,\ (HA)^\top=HA\) | Penrose-properties form |
| \(P_{\mathrm{PMN}3}^1\) | \(AA^\top H^\top=A,\ (AH)^\top=AH\) | “PMN+P3” form |
| \(P_{\mathrm{PLS\,PMN}}^1\) | \(A^\top A\,H=A^\top,\ HAA^\top=A^\top\) | “PLS+PMN” form |
| \(P_{\mathcal{RR}^\top}^1\) | \(AH=AA^\dagger,\ HA=A^\dagger A\) | Double-projector form |
| \(P_{\mathrm{PMX}}^1\) | \(AA^\top H^\top + H^\top A^\top A = 2A^\top\) | Reduced linear equation |
| \(\mathcal{P}_{134}^1\) | \(H = V_1D^{-1}U_1^\top + V_2WU_2^\top\) | SVD-parameterization form |

Each of the first five formulations can be cast as a standard-form LP in variables \(H^+,H^-\ge 0\) with \(H=H^+-H^-\). The SVD form reduces the optimization to the free variable \(W\), with only \((n-r)(m-r)\) variables, and is therefore described as often the most compact formulation [2509.04264].

The same paper derives a sparsity bound for extreme points of the LP reformulation of \(P_{134}^1\): every extreme solution has at most
\[
mn-(m-r)(n-r)
\]
nonzeros. The proof sketch given in the summary applies \(\mathrm{vec}(\cdot)\) to obtain the single linear system
\[
(I_m\otimes A^\top A + AA^\top\otimes I_n)\,\mathrm{vec}(H)=2\,\mathrm{vec}(A^\top),
\]
computes its matrix rank as \(mn-(m-r)(n-r)\), and then invokes the standard LP fact that an extreme point over \(p\) linearly independent equations in \(\mathbb{R}^N\) has at most \(p\) nonzeros [2509.04264].

This framework separates two issues that are often conflated. The Penrose constraints enforce simultaneous least-squares and minimum-norm behavior for all \(b\); the \(1\)-norm objective is a secondary design criterion used to induce sparsity. A plausible implication is that universal solvers can be tailored to implementation goals—such as sparsity or communication reduction—without changing the target mapping \(b\mapsto Hb\).

## 3. Proximal-point and splitting algorithms for sparse simultaneous solvers

The computational problem in [2509.04264] is written in the standard “sum-of-two-convex” form
\[
\min_H f(H)+g(H),
\]
with \(f(H)=\|H\|_1\) and \(g(H)=I_C(H)\), where \(C\) is the affine constraint set of one of the equivalent formulations. Two algorithmic routes are emphasized: Douglas–Rachford Splitting (DRS) and ADMM-like methods.

For DRS, the iterates are
\[
\begin{aligned}
H^{k+\tfrac12} &:= \mathrm{prox}_{\lambda f}(V^k),\\
V^{k+\tfrac12} &:= 2H^{k+\tfrac12}-V^k,\\
H^{k+1} &:= \mathrm{prox}_{\lambda g}(V^{k+\tfrac12}),\\
V^{k+1} &:= V^k + H^{k+1}-H^{k+\tfrac12}.
\end{aligned}
\]
The proximal map of \(f\) is elementwise soft-thresholding,
\[
\mathrm{prox}_{\lambda f}(V)=S_\lambda(V),\qquad
S_\lambda(v)=\mathrm{sign}(v)\max\{|v|-\lambda,0\},
\]
while the proximal map of \(g\) is Euclidean projection onto the affine set \(C\) [2509.04264].

For the simultaneous \(P_{134}^1\) problem, the paper sets
\[
C = \{H : A^\top A\,H = A^\top,\ HAA^\top = A^\top\},
\]
and gives the closed-form projection
\[
\Pi_C(V)
=
V
-A^\dagger A\,V
+A^\dagger
-VAA^\dagger
+A^\dagger A\,VAA^\dagger.
\]
This closed form is what makes DRS practical in the full matrix variable \(H\) [2509.04264].

An ADMM variant is developed for the SVD-parameterized problem \(\mathcal{P}_{134}^1\). Introducing \(H\) explicitly and enforcing
\[
H=V_1D^{-1}U_1^\top+V_2WU_2^\top,
\]
the updates alternate between a least-squares step in \(W\), a soft-thresholding step in \(H\), and a dual update. In the summary, the \(W\)-update is
\[
W^{k+1}=V_2^\top V^kU_2,
\]
followed by
\[
H^{k+1}=S_{1/\rho}\!\left(V_1D^{-1}U_1^\top+V_2W^{k+1}U_2^\top+\Lambda^k\right).
\]
A similar ADMM is stated for the related \(\mathcal{P}_{123}^1\) problem [2509.04264].

The empirical comparison in [2509.04264] is explicit. DRS\(_{fp}\) is compared against direct LP via Gurobi and against the ADMM on the SVD form. On random dense \(A\in\mathbb{R}^{m\times n}\) of rank \(r\approx 0.25m\) with \(m=n\) up to \(5000\), DRS\(_{fp}\) “routinely solves in hundreds to a few thousand seconds,” Gurobi “times out beyond \(m\approx 100\),” and the ADMM is “significantly slower than DRS\(_{fp}\) (sometimes by an order of magnitude).” A representative example with \(m=n=2000\), \(r=500\), reports \(\|H\|_0\approx 1.0\times 10^7\) versus \(A^\dagger\)’s \(\sim 2.0\times 10^7\), \(\|H\|_1\approx 1.3\times 10^3\) versus \(A^\dagger\)’s \(\sim 1.3\times 10^4\), in \(\approx 300\,\mathrm{s}\), while Gurobi did not finish in \(7200\,\mathrm{s}\) [2509.04264].

## 4. Iterative universal solvers: Krylov, row-action, and projection methods

In numerical linear algebra, the same universal objective is often realized without forming a global operator \(H\). Instead, an iterative method maps \(b\) to either a least-squares minimizer or the pseudoinverse solution, depending on consistency and rank structure.

For Hermitian or complex-symmetric problems, "Obtaining Pseudo-inverse Solutions With MINRES" shows that standard MINRES solves
\[
x_t=\arg\min_{x\in K_t(A,b)}\|b-Ax\|
\]
but, when \(b\notin \mathrm{range}(A)\), the final iterate \(x_g\) has minimal residual without necessarily having minimal norm. The paper introduces a minimum-norm “lifting” refinement at the final step:
\[
\hat x^* = x_g - \frac{\langle r_g,x_g\rangle}{\|r_g\|^2}r_g,
\qquad r_g=b-Ax_g,
\]
and proves \(\hat x^*=A^+b\). For complex-symmetric \(A\), the refinement becomes
\[
\hat x^* = x_g - \frac{\langle \overline r_g,x_g\rangle}{\|r_g\|^2}\,\overline r_g.
\]
The added work is one dot product, one saxpy, and one norm, so the extra cost is negligible and no additional Krylov basis storage is needed [2309.17096].

MINRES-QLP takes a different route. It augments the MINRES QR treatment of the Lanczos tridiagonal matrix by a QLP decomposition, and in exact arithmetic terminates in at most \(\ell\) steps with \(x_\ell=x^\dagger\), the pseudoinverse solution, whether \(A\) is nonsingular, singular compatible, or singular incompatible [1003.4042]. The algorithmic paper emphasizes that if the system is singular, MINRES-QLP computes the unique minimum-length solution, which generally eludes MINRES, and does so while preserving short recurrences and allowing a positive-definite preconditioner [1301.2707].

For arbitrary matrices, Zouzias–Freris’ Randomized Extended Kaczmarz interleaves a column-projection loop for the inconsistent component with randomized Kaczmarz steps on the evolving right-hand side. The result is an algorithm that “exponentially converges in expectation to the minimum Euclidean norm least squares solution,” with overall expected flop count
\[
O\!\Bigl(\mathrm{nnz}(A)\,\kappa^2(A)\,\ln\frac{\kappa^2(A)}{\eta}\Bigr)
\]
for success probability \(1-\eta\) [1205.5770].

Sugihara and Hayami propose applying RRGMRES to
\[
ACA^\top z=b,\qquad x=CA^\top z,
\]
with \(C\) symmetric positive definite or approximated by an \(\ell\)-step NR-SSOR inner iteration. Their theory states that if \(b\in\mathrm{Range}(A)\), the method converges in at most \(r\) steps to the unique minimum-norm solution, and if \(b\notin\mathrm{Range}(A)\), it converges in at most \(r\) steps to the least-squares solution \(x^*\). The reported experiments state that NR-SSOR preconditioning with \(\ell=1\) or \(\ell=4\) yields minimal residuals \(\approx 10^{-14}\) in \(1/3\)–\(1/2\) as many iterations, and runs \(5\times\)–\(10\times\) faster than the comparison methods [2504.09891].

Kalantari’s family of iteration functions
\[
F_t(r)=r-\sum_{i=1}^t \alpha_{t,i}(r)H^ir
\]
gives another unified construction. When applied to \(H=AA^\top\) or directly to the normal equations, it yields \(\varepsilon\)-approximate minimum-norm solutions for consistent systems or \(\varepsilon\)-approximate least-squares solutions for inconsistent ones; if \(s\) is the degree of the minimal polynomial of the residual with respect to \(H\), then \(F_s(r)\) gives the minimum-norm solution of \(Ax=b\) or an exact solution of \(A^\top Ax=A^\top b\) in \(O(sN+s^3)\) operations [2304.04940].

These methods differ in recurrence structure, storage, and admissible matrix classes. This suggests that “universality” is not a single algorithmic design but a shared target condition: exact or asymptotic recovery of the least-squares minimizer together with minimum-norm selection when the null space is nontrivial.

## 5. Direct factorizations and large-scale complete decompositions

A direct-factorization analogue appears in the 2025 symmetric-indefinite algorithm of Coria, Urkullu, Uriarte, and Fernández-de-Bustos. The method factors
\[
A = MLDL^\top M^\top,
\]
with \(M\) a product of symmetric permutation–rotation matrices, \(L\) unit lower-triangular, and \(D\) block-diagonal, using Jacobi rotations and Rook’s pivoting. From the rank-revealing structure, it extracts a fundamental null basis
\[
Z = M\begin{bmatrix}
-L_{11}^{-\top}L_{21}^\top\\
I_{n-r}
\end{bmatrix},
\]
writes any least-squares solution as \(x_{\mathrm{LS}}=x_{\mathrm{part}}+Zy\), and then selects the minimum-norm solution through
\[
x_{\min}=x_{\mathrm{part}}-Z(Z^\top Z)^{-1}Z^\top x_{\mathrm{part}}.
\]
The summary states \(O(n^3)\) leading cost, “same \(1/3\,n^3\) as Bunch–Kaufman,” and reports that for determinate compatible systems the error is “approximately 50 % smaller” than Bunch–Kaufman, while in minimal least squares with minimum norm problems the computational cost is “at least 20 % smaller” than Complete Orthogonal Decomposition [2501.17696].

For matrices too large to fit in main memory, the randUTV-based out-of-core methods in [2408.05238] use complete orthogonal decompositions that “guarantee that both conditions of a least squares solution are met, regardless of the rank properties of the matrix.” With
\[
A=UTV^\top,
\]
where \(T\) is trapezoidal and rank revealing, the solver first computes a least-squares minimizer and then, in “safe” mode, performs an extra small COD/SVD on \([T_{11}\ T_{12}]\) so that
\[
x_{\mathrm{minnorm}} = V\begin{bmatrix}\tilde x\\0\end{bmatrix},
\qquad
\tilde x=S_{11}^{-1}(U_1^\top b).
\]
The performance study reports that \(\|Ax_{\mathrm{minnorm}}-b\|_2\) and \(\|x_{\mathrm{minnorm}}\|_2\) agree to within round-off with Intel MKL’s dGELSS and dGELSY solvers, with “typical residuals” \(10^{-12}\)–\(10^{-6}\), and gives explicit out-of-core timings such as \(\approx 6640\,\mathrm{s}\) on CPU and \(\approx 4900\,\mathrm{s}\) on GPU for a \(92\,160\times 92\,160\) rank-\(90000\) least-squares problem [2408.05238].

Both papers show that universal minimum-norm/least-squares behavior can be embedded in direct solvers rather than only in iterative ones. One does so through an explicit null-basis correction after \(LDL^\top\); the other does so through a complete orthogonal decomposition coupled to a final minimum-norm enforcement step.

## 6. Scope, misconceptions, and open problems

A recurring misconception is that a least-squares minimizer is automatically the minimum-norm minimizer. The MINRES refinement paper states the general solution of
\[
\min_x \|b-Ax\|^2
\]
as
\[
x=A^+b+(I-A^+A)y,
\]
with \(y\) arbitrary. Thus the least-squares set is affine along \(\mathrm{Null}(A)\), and minimum norm requires an additional selection principle [2309.17096]. This same distinction underlies the simultaneous Penrose constraints \(P1,P3,P4\), where least-squares and minimum-norm requirements are imposed together but are not redundant [2509.04264].

A second misconception is that universal behavior requires a dense SVD or a complete pseudoinverse construction. The literature summarized here shows otherwise. Universal or near-universal behavior is realized by a one-step lifting on top of MINRES, by a QLP correction to Lanczos tridiagonalization, by randomized row-and-column projections, by right-preconditioned RRGMRES on \(ACA^\top\), by null-basis corrections after \(LDL^\top\), and by out-of-core complete orthogonal decomposition [2309.17096]. A plausible implication is that the universal property is best understood as a target specification, not as a commitment to any one numerical primitive.

The current open problems are stated most explicitly in the MINRES lifting work: “dynamic stopping criteria that balance residual vs. norm,” “the design of optimal singular sub-preconditioners \(S\),” and “extensions to other Krylov methods (e.g. GMRES)” [2309.17096]. The sparse-solver work suggests an additional direction: computing simultaneous universal solvers that are not merely correct but structurally economical, with DRS currently the strongest performer among the compared convex-optimization approaches [2509.04264].

Across these strands, the subject has a common organizing principle. A simultaneous universal minimum-norm/least-squares solver must act correctly on every right-hand side, including incompatible and rank-deficient cases, and must do so while distinguishing residual minimization from norm minimization. Whether that goal is expressed through Penrose equations, Krylov recurrences, row-action iterations, or rank-revealing factorizations, the central problem is the same: constructing a map that reproduces least-squares behavior globally and the pseudoinverse solution whenever exact solvability lies in \(\mathrm{Range}(A)\).

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