Universal Min-Norm/Least-Squares Solvers

• Simultaneous universal min-norm/least-squares solvers are operators that, for every right-hand side, compute a least-squares solution and enforce minimum-norm behavior on the range of the matrix.
• They employ a combination of algebraic formulations, including Penrose conditions and SVD-parameterizations, to compactly characterize and optimize solver properties with sparsity promotion.
• Empirical analyses and iterative algorithms such as DRS, ADMM, and Krylov methods illustrate their efficiency and robustness across consistent, singular, and rank-deficient regimes.

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" (Machado et al., 4 Sep 2025), 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 (Machado et al., 4 Sep 2025).

1. Algebraic definition and Penrose characterization

For a matrix $A$ of rank $A\in\mathbb{R}^{m\times n}$0, a simultaneous universal solver is a matrix $A\in\mathbb{R}^{m\times n}$1 such that $A\in\mathbb{R}^{m\times n}$2 is, for every $A\in\mathbb{R}^{m\times n}$3, both a least-squares solution and, on $A\in\mathbb{R}^{m\times n}$4, the minimum-$A\in\mathbb{R}^{m\times n}$5-norm exact solution. The paper on sparse universal solvers states that this is equivalent to the three Penrose conditions

$A\in\mathbb{R}^{m\times n}$6

denoted $A\in\mathbb{R}^{m\times n}$7, $A\in\mathbb{R}^{m\times n}$8, and $A\in\mathbb{R}^{m\times n}$9 (Machado et al., 4 Sep 2025).

With the full SVD

$b$0

where $b$1 is diagonal with positive singular values, the same paper gives the parameterization

$b$2

with $b$3. 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 $b$4 spans the null-space degrees of freedom (Machado et al., 4 Sep 2025).

A key reduced characterization is the equivalence

$b$5

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 (Machado et al., 4 Sep 2025).

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 (Choi et al., 2010). 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 $b$6-norm $b$7. The paper lists six equivalent linear-constraint formulations, all with objective $b$8 (Machado et al., 4 Sep 2025).

Formulation	Constraints	Comment
$b$9	$Hb$0	Penrose-properties form
$Hb$1	$Hb$2	“PMN+P3” form
$Hb$3	$Hb$4	“PLS+PMN” form
$Hb$5	$Hb$6	Double-projector form
$Hb$7	$Hb$8	Reduced linear equation
$Hb$9	$\min_\theta \|A\theta-b\|_2^2$0	SVD-parameterization form

Each of the first five formulations can be cast as a standard-form LP in variables $\min_\theta \|A\theta-b\|_2^2$1 with $\min_\theta \|A\theta-b\|_2^2$2. The SVD form reduces the optimization to the free variable $\min_\theta \|A\theta-b\|_2^2$3, with only $\min_\theta \|A\theta-b\|_2^2$4 variables, and is therefore described as often the most compact formulation (Machado et al., 4 Sep 2025).

The same paper derives a sparsity bound for extreme points of the LP reformulation of $\min_\theta \|A\theta-b\|_2^2$5: every extreme solution has at most

$\min_\theta \|A\theta-b\|_2^2$6

nonzeros. The proof sketch given in the summary applies $\min_\theta \|A\theta-b\|_2^2$7 to obtain the single linear system

$\min_\theta \|A\theta-b\|_2^2$8

computes its matrix rank as $\min_\theta \|A\theta-b\|_2^2$9, and then invokes the standard LP fact that an extreme point over $b\in \mathrm{Range}(A)$0 linearly independent equations in $b\in \mathrm{Range}(A)$1 has at most $b\in \mathrm{Range}(A)$2 nonzeros (Machado et al., 4 Sep 2025).

This framework separates two issues that are often conflated. The Penrose constraints enforce simultaneous least-squares and minimum-norm behavior for all $b\in \mathrm{Range}(A)$3; the $b\in \mathrm{Range}(A)$4-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\in \mathrm{Range}(A)$5.

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

The computational problem in (Machado et al., 4 Sep 2025) is written in the standard “sum-of-two-convex” form

$b\in \mathrm{Range}(A)$6

with $b\in \mathrm{Range}(A)$7 and $b\in \mathrm{Range}(A)$8, where $b\in \mathrm{Range}(A)$9 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

$Hb$0

The proximal map of $Hb$1 is elementwise soft-thresholding,

$Hb$2

while the proximal map of $Hb$3 is Euclidean projection onto the affine set $Hb$4 (Machado et al., 4 Sep 2025).

For the simultaneous $Hb$5 problem, the paper sets

$Hb$6

and gives the closed-form projection

$Hb$7

This closed form is what makes DRS practical in the full matrix variable $Hb$8 (Machado et al., 4 Sep 2025).

An ADMM variant is developed for the SVD-parameterized problem $Hb$9. Introducing $2$0 explicitly and enforcing

$2$1

the updates alternate between a least-squares step in $2$2, a soft-thresholding step in $2$3, and a dual update. In the summary, the $2$4-update is

$2$5

followed by

$2$6

A similar ADMM is stated for the related $2$7 problem (Machado et al., 4 Sep 2025).

The empirical comparison in (Machado et al., 4 Sep 2025) is explicit. DRS$2$8 is compared against direct LP via Gurobi and against the ADMM on the SVD form. On random dense $2$9 of rank $A\theta=b$0 with $A\theta=b$1 up to $A\theta=b$2, DRS$A\theta=b$3 “routinely solves in hundreds to a few thousand seconds,” Gurobi “times out beyond $A\theta=b$4,” and the ADMM is “significantly slower than DRS$A\theta=b$5 (sometimes by an order of magnitude).” A representative example with $A\theta=b$6, $A\theta=b$7, reports $A\theta=b$8 versus $A\theta=b$9’s $A$0, $A$1 versus $A$2’s $A$3, in $A$4, while Gurobi did not finish in $A$5 (Machado et al., 4 Sep 2025).

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 $A$6. Instead, an iterative method maps $A$7 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

$A$8

but, when $A$9, the final iterate $A\in\mathbb{R}^{m\times n}$00 has minimal residual without necessarily having minimal norm. The paper introduces a minimum-norm “lifting” refinement at the final step: $A\in\mathbb{R}^{m\times n}$01 and proves $A\in\mathbb{R}^{m\times n}$02. For complex-symmetric $A\in\mathbb{R}^{m\times n}$03, the refinement becomes

$A\in\mathbb{R}^{m\times n}$04

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 (Liu et al., 2023).

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 $A\in\mathbb{R}^{m\times n}$05 steps with $A\in\mathbb{R}^{m\times n}$06, the pseudoinverse solution, whether $A\in\mathbb{R}^{m\times n}$07 is nonsingular, singular compatible, or singular incompatible (Choi et al., 2010). 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 (Choi et al., 2013).

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

$A\in\mathbb{R}^{m\times n}$08

for success probability $A\in\mathbb{R}^{m\times n}$09 (Zouzias et al., 2012).

Sugihara and Hayami propose applying RRGMRES to

$A\in\mathbb{R}^{m\times n}$10

with $A\in\mathbb{R}^{m\times n}$11 symmetric positive definite or approximated by an $A\in\mathbb{R}^{m\times n}$12-step NR-SSOR inner iteration. Their theory states that if $A\in\mathbb{R}^{m\times n}$13, the method converges in at most $A\in\mathbb{R}^{m\times n}$14 steps to the unique minimum-norm solution, and if $A\in\mathbb{R}^{m\times n}$15, it converges in at most $A\in\mathbb{R}^{m\times n}$16 steps to the least-squares solution $A\in\mathbb{R}^{m\times n}$17. The reported experiments state that NR-SSOR preconditioning with $A\in\mathbb{R}^{m\times n}$18 or $A\in\mathbb{R}^{m\times n}$19 yields minimal residuals $A\in\mathbb{R}^{m\times n}$20 in $A\in\mathbb{R}^{m\times n}$21–$A\in\mathbb{R}^{m\times n}$22 as many iterations, and runs $A\in\mathbb{R}^{m\times n}$23–$A\in\mathbb{R}^{m\times n}$24 faster than the comparison methods (Sugihara et al., 14 Apr 2025).

Kalantari’s family of iteration functions

$A\in\mathbb{R}^{m\times n}$25

gives another unified construction. When applied to $A\in\mathbb{R}^{m\times n}$26 or directly to the normal equations, it yields $A\in\mathbb{R}^{m\times n}$27-approximate minimum-norm solutions for consistent systems or $A\in\mathbb{R}^{m\times n}$28-approximate least-squares solutions for inconsistent ones; if $A\in\mathbb{R}^{m\times n}$29 is the degree of the minimal polynomial of the residual with respect to $A\in\mathbb{R}^{m\times n}$30, then $A\in\mathbb{R}^{m\times n}$31 gives the minimum-norm solution of $A\in\mathbb{R}^{m\times n}$32 or an exact solution of $A\in\mathbb{R}^{m\times n}$33 in $A\in\mathbb{R}^{m\times n}$34 operations (Kalantari, 2023).

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\in\mathbb{R}^{m\times n}$35

with $A\in\mathbb{R}^{m\times n}$36 a product of symmetric permutation–rotation matrices, $A\in\mathbb{R}^{m\times n}$37 unit lower-triangular, and $A\in\mathbb{R}^{m\times n}$38 block-diagonal, using Jacobi rotations and Rook’s pivoting. From the rank-revealing structure, it extracts a fundamental null basis

$A\in\mathbb{R}^{m\times n}$39

writes any least-squares solution as $A\in\mathbb{R}^{m\times n}$40, and then selects the minimum-norm solution through

$A\in\mathbb{R}^{m\times n}$41

The summary states $A\in\mathbb{R}^{m\times n}$42 leading cost, “same $A\in\mathbb{R}^{m\times n}$43 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 (Coria et al., 29 Jan 2025).

For matrices too large to fit in main memory, the randUTV-based out-of-core methods in (Chillarón et al., 2024) 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\in\mathbb{R}^{m\times n}$44

where $A\in\mathbb{R}^{m\times n}$45 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 $A\in\mathbb{R}^{m\times n}$46 so that

$A\in\mathbb{R}^{m\times n}$47

The performance study reports that $A\in\mathbb{R}^{m\times n}$48 and $A\in\mathbb{R}^{m\times n}$49 agree to within round-off with Intel MKL’s dGELSS and dGELSY solvers, with “typical residuals” $A\in\mathbb{R}^{m\times n}$50–$A\in\mathbb{R}^{m\times n}$51, and gives explicit out-of-core timings such as $A\in\mathbb{R}^{m\times n}$52 on CPU and $A\in\mathbb{R}^{m\times n}$53 on GPU for a $A\in\mathbb{R}^{m\times n}$54 rank-$A\in\mathbb{R}^{m\times n}$55 least-squares problem (Chillarón et al., 2024).

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 $A\in\mathbb{R}^{m\times n}$56; 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

$A\in\mathbb{R}^{m\times n}$57

as

$A\in\mathbb{R}^{m\times n}$58

with $A\in\mathbb{R}^{m\times n}$59 arbitrary. Thus the least-squares set is affine along $A\in\mathbb{R}^{m\times n}$60, and minimum norm requires an additional selection principle (Liu et al., 2023). This same distinction underlies the simultaneous Penrose constraints $A\in\mathbb{R}^{m\times n}$61, where least-squares and minimum-norm requirements are imposed together but are not redundant (Machado et al., 4 Sep 2025).

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 $A\in\mathbb{R}^{m\times n}$62, by null-basis corrections after $A\in\mathbb{R}^{m\times n}$63, and by out-of-core complete orthogonal decomposition (Liu et al., 2023). 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 $A\in\mathbb{R}^{m\times n}$64,” and “extensions to other Krylov methods (e.g. GMRES)” (Liu et al., 2023). 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 (Machado et al., 4 Sep 2025).

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 $A\in\mathbb{R}^{m\times n}$65.