On computing sparse universal solvers for key problems in statistics

Ananias Sousa Machado; Jon Lee; Marcia Fampa

REVIEW 2 major objections 5 minor 1 cited by

On computing sparse universal solvers for key problems in statistics

T0 review · 2 major / 5 minor · reviewed 2026-08-05 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

Pith's one-line read Douglas-Rachford splitting finds sparse universal solvers for least-squares problems much faster than LP-based alternatives.

desk verdict Solid math, nice new bound, but the 'performed best' claim needs the experimental details before it's believable. read the letter →

arxiv 2509.04264 v2 pith:DXIVVHOC submitted 2025-09-04 math.OC

Ananias Sousa Machado , Marcia Fampa , Jon Lee This is my paper · ORCID

classification math.OC MSC 15A0990C2590C05

keywords leastsquaresuniversalsolvergeneralizedinverseMoore-Penrosepseudoinverseminimum1-normsparsityDouglas-Rachfordsplittinglinearprogramming

verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

Least-squares problems with many different right-hand sides can be handled by one precomputed matrix, and the standard choice, the Moore-Penrose pseudoinverse, is usually dense even when the data matrix is sparse. This paper looks for sparse stand-ins that still solve least-squares, minimum-norm, or both, for every right-hand side, and induces sparsity by minimizing the sum of the absolute values of the entries. It gives several equivalent linear-programming formulations for three problem classes, proves a new upper bound on the number of nonzeros in extreme solutions of the simultaneous-solver LP, and develops ADMM and Douglas-Rachford splitting algorithms. The central computational claim is that a Douglas-Rachford variant driven by a fixed-point stopping rule solves the tested instances up to 5000-by-5000 far faster than a commercial LP solver and faster than ADMM, while producing comparably sparse solutions. If correct, practitioners could replace repeated multiplications by a dense pseudoinverse with one cheaper sparse solver computed in advance.

What carries the argument

The engine is the conversion of the four Moore-Penrose conditions into affine linear systems with closed-form projections. In SVD coordinates, a candidate H has block form V Gamma U^T; P1 fixes one block to D^{-1}, P3 zeroes another, P4 zeroes a third, leaving compact LP formulations. For the splitting algorithms, the same conditions appear as A^T A H = A^T (least squares), A^T A H = A^T plus H A A^+ = H (minimum rank), and A^T A H + H A A^T = 2A^T (simultaneous solver). Each feasible set C is affine, so the projection has a closed form, and the Douglas-Rachford iteration reduces to soft thresholding S_lambda and projection: V -> V + Pi_C(2S_lambda(V)-V)-S_lambda(V).

What would settle it

Run DRS_fp and a commercial LP solver on the same three problem classes using real or differently conditioned design matrices and several step-size/penalty settings, timing them at comparable solution quality; if the LP solver wins or the gap disappears, the headline computational claim fails. For the sparsity theorem, solve one small LP instance of the simultaneous-solver problem and count nonzeros of an extreme solution against mn-(m-r)(n-r).

Watch

Extended reading notes

Core claim

The paper claims a theoretical and a computational result. Theoretically, extreme solutions of the standard LP reformulation of the simultaneous least-squares/minimum-norm solver have at most mn-(m-r)(n-r) nonzeros; the proof computes the rank of (I_m tensor A^T A)+(A A^T tensor I_n). Computationally, Douglas-Rachford splitting applied to affine characterizations of the Moore-Penrose conditions solves instances up to 5000x5000 in minutes to tens of minutes where a commercial LP solver times out and ADMM is slower. Returned sparse solvers have about half the pseudoinverse's nonzeros; enforcing minimum rank makes them visibly denser.

Load-bearing premise

The speed comparison rests on one family of random test matrices with a fixed singular-value profile and on step-size, penalty, and tolerance values that the paper does not report; with different real data or different settings, the ranking could change, and the algorithms' convergence is inherited from standard DRS theory rather than proved for these projections.

Editorial extensions

If this is right

Sparse universal least-squares and simultaneous solvers become computable at scales where LP solvers time out: the DRS fixed-point variant handles 5000-by-5000 instances in minutes to tens of minutes rather than timing out after two hours.
On the tested matrices, minimizing the 1-norm cuts the number of nonzeros roughly in half relative to the Moore-Penrose pseudoinverse for the least-squares and simultaneous classes.
Imposing minimum rank makes the minimal-1-norm solvers noticeably denser, quantifying a real trade-off between rank and sparsity for universal solvers.
Every extreme solution of the simultaneous-solver LP obeys the new mn-(m-r)(n-r) sparsity ceiling, so LP-based 1-norm minimization cannot produce sparser universal simultaneous solvers than that bound.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

Because the tests use one random matrix family with a fixed singular-value profile, the natural next check is DRS_fp against the same LP baseline on ill-conditioned or application-specific design matrices; the ranking may not transfer.
The rank computation behind the new bound ties sparsity to LP degeneracy: if the bound is loose, the improvement must come from showing every basic feasible solution has forced zero variables, which also suggests a route to tightening the minimum-rank bound.
The DRS fixed-point map is exactly the kind of map Anderson acceleration targets, so the already fast variant may be made faster; the paper names this as future work.
For the open NP-hardness question about sparsity maximization under P1+P3+P4, the new bound does not settle complexity, but it constrains the structure any hard instances would need to have.

Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Desk Editor's Note

read the letter

Plain take: the core mathematics is in good shape, the experimental claim is not yet fully supported. The new Douglas-Rachford algorithms with closed-form projections are a real step beyond the authors' prior ADMM work, and Theorem 19 — the mn − (m−r)(n−r) sparsity bound for extreme LP solutions of the simultaneous solver problem — is a clean, genuinely new result. The SVD-block proofs of Theorems 9, 17, and 19 check out, and the final projection formulas in Propositions 7, 11, and 20 are standard; I spot-checked them and they behave.

The soft spots are in presentation and reproducibility. The printed derivations of those three projections contain intermediate algebra that doesn't parse — dimension mismatches and dropped terms — so the proofs as written should not appear in a final version. I suspect the results are right, but a referee will need to have the authors rewrite those steps carefully.

The bigger concern is the computational section. The headline speedups over Gurobi and the prior ADMM come from experiments where the DRS step size λ, the ADMM penalty ρ, and the tolerances ϵ_abs and ϵ_rel are never reported. The final tables show a single instance per matrix size, with no seeds, no error bars, and no performance profiles. The instance-generation description in §4.4 is ambiguous: the singular-value vector is written as 'M × (α1,...,αr)' but α is then given as a scalar. No code or data is shipped. So the 'performed best' statement is plausible — I'd expect a first-order method with closed-form projections to beat Gurobi on large LPs of this kind — but the paper as written doesn't allow an independent check. The appendix does average over five instances for the preliminary Gurobi comparison, which helps a little, but the main tables don't.

The self-citation is heavy but not circular; the new algorithms and bound stand on their own, and the comparisons are against independent baselines. The paper is a solid, incremental contribution to a program the authors have been developing for years. People working on sparse generalized inverses, universal solvers, or DRS/ADMM applications will get value from it. I'd send it to peer review, with the requested revisions being: rewrite the projection proofs, report all hyperparameters and seeds, and either publish code/data or add reproducibility details.

Referee Report

2 major / 5 minor

Summary. The paper studies the computation of sparse universal solvers—generalized inverses H that return optimal solutions for every right-hand side b—for three problem classes: least-squares (characterized by P1+P3), minimum-rank least-squares (P1+P2+P3), and simultaneous least-squares/minimum-2-norm (P1+P3+P4). For each class the authors present several LP-reformulable formulations, develop proximal algorithms (DRS for the first and third classes and for the minimum-rank class, plus an ADMM for P1+P3+P4), prove a new sparsity bound for extreme points of the simultaneous-solver LP (Theorem 19), and report computational comparisons on random instances. The headline claim is that the new Douglas-Rachford variants, especially DRSfp, outperform both Gurobi on the LP formulations and the existing ADMM.

Significance. The manuscript has a solid theoretical core: Theorem 19's bound mn-(m-r)(n-r) is a new and nontrivial extension of the earlier bounds, and the SVD-based characterizations in Theorems 9 and 17 appear correct. The paper also usefully catalogs many LP formulations and shows how solver performance varies dramatically with formulation. The final projection formulas in Propositions 7, 11, and 20 pass direct feasibility checks, and the algorithms are natural applications of standard DRS/ADMM machinery. The main weakness is empirical: the speed comparison that supports the title claim is not yet reproducible because key hyperparameters, seeds, and instance-generation details are missing, and the main tables report single runs. With a reproducibility appendix and corrected proofs, the paper would be a valuable contribution.

major comments (2)

[§1 (Experimental setup), §2.2, §3.3, §4.4, Appendix] The computational claim that DRSfp 'performed best' is the paper's headline, but the experiments are not reproducible as reported. The DRS step size λ and ADMM penalty ρ are never specified; the tolerances ϵ_abs, ϵ_rel in the stopping criteria of §2.1.3, §3.1 and §4.2 are never given numerical values. The final tables (Tables 1, 3, 5) show a single row per instance size, with no seeds, no error bars, and no performance profiles; the Appendix averages five instances, but the main tables do not. Instance generation is also ambiguous: §4.4 defines rc as M(α1,...,αr) with M=2 and then sets α=(1/M)(2/(r+1)), so the intended decreasing sequence is not specified. These omissions must be repaired (and preferably code/data released) before the speed comparison can be assessed.
[§2.1.1, §3.2, §4.3 (Propositions 7, 11, 20)] The three closed-form projection formulas are central to the DRS implementations. Each formula can be checked by direct substitution and is consistent with standard projections onto affine subspaces, but the printed proofs contain invalid algebra. In Proposition 7, after deriving A^T = A^T A V - (1/2) A^T A(A^T AΛ), the proof multiplies by A†^T and replaces A†^T A^T A A^T AΛ by A A^T AΛ; this is not a valid simplification and the equation sizes do not match. Analogous dimensionally inconsistent manipulations appear in Propositions 11 and 20. Since these proofs are the only derivations of the projections used by the algorithms, they must be rewritten carefully.

minor comments (5)

[Theorem 17 / §4.1] Calling P1+P3+P4 'reflexive ah-ha-symmetric' is not correct: P2 does not follow from P1+P3+P4 (e.g., A=diag(1,0), H=I2). The equivalence with PMX is fine; please adjust the terminology and any related statements.
[Table 6] The caption defines β134 := mn + (m-r)(n-r), but Theorem 19 proves the bound mn - (m-r)(n-r). The numerical ratios in the table are consistent with the minus version, so this appears to be a typo; please correct.
[§3.1] The output of the ADMM for P^1_123 is written as V1D^{-1}U1^T + V2Z^{k+1}; the second term should be V2Z^{k+1}U1^T. Please fix.
[Lemma 18] Lemma 18 is quoted from [XFLP21] without proof; a one-line proof (rank of [B,-B] equals rank B) would make Theorem 19 self-contained.
[Figures 1 and 2] The plots lack axis labels and legend details in the text; state explicitly what is plotted and for which instance/run. Also, the convergence plots are for a single instance; label them as such.

Circularity Check

0 steps flagged · score 0.0 of 10

No significant circularity: the new DRS/ADMM algorithms and Theorem 19 are self-contained; self-citations are standard prior results or external facts.

full rationale

I examined the derivation chain for each claimed contribution. The optimization formulations (P1_13, P1_123, P1_134 and their variants) are equivalent reformulations of the same feasible sets, established by theorems (Theorems 2–5, 9, 16–17) whose proofs are given or cited to external standard references. The DRS algorithms are derived from the reformulations via proximal operators and projections; the projection formulas and residuals are proved in the paper (Propositions 7, 11, 20, 21) rather than assumed. The new sparsity bound (Theorem 19) uses Lemma 18, cited to [XFLP21], but Lemma 18 is a standard elementary LP fact about extreme points having at most rank(B) nonzeros; its use does not assume the target result, and the rank computation in Theorem 19 is an independent algebraic calculation. The numerical comparisons are against independent baselines (Gurobi, A†, the earlier ADMM), and the DRSfp output is not fed back as an input into the optimization problems. The stopping protocol that matches the 1-norm of DRSfp is a legitimate performance comparison, not a fitted-input prediction. The paper does self-cite frequently and reuses prior formulations and bounds, but these citations are either external results (BG74, Penrose), standard facts, or clearly identified previous work that does not smuggle in the conclusions. Missing hyperparameters (λ, ρ, tolerances), unreported seeds, and single-run experiments are reproducibility concerns, not circularity. Therefore no circular step can be exhibited, and the score is 0.

Assumptions & free parameters 4 free parameters · 4 assumptions · 0 invented entities

The mathematical claims rest on standard linear algebra and convex optimization results. The only hand-set quantities are algorithm hyperparameters and the 0-norm threshold, none of which are reported with values except the threshold. No new physical or mathematical entities are introduced.

free parameters (4)

DRS step size lambda
Appears in prox_{lambda ||.||_1} and residual definitions in Sections 2.1, 3.2, and 4.3; the value is not reported.
ADMM penalty rho
Augmented Lagrangian penalty in Sections 3.1 and 4.2; the value is not reported.
Stopping tolerances epsilon_abs and epsilon_rel
Used in DRS_res, DRS_fp, and ADMM stopping criteria; values are not given in the paper.
0-norm threshold = 1e-5
All sparsity counts treat entries below 1e-5 as zero; this hand-chosen threshold directly defines the 0-norm tables.

assumptions (4)

standard math Penrose characterization of the Moore-Penrose pseudoinverse and the SVD block characterization of the four Moore-Penrose properties
Used throughout to define feasible sets and to derive equivalent formulations; quoted from Pen55 and BG74/PFLX24b.
standard math Universal-solver characterizations: P1+P3 iff least-squares solver, P1+P4 iff minimum-norm solver, AH=AA†, PLS, PMN, PLSr, PMX
These equivalences from BG74 and the paper's Theorems 5, 9, 16, 17 are the basis for the LP and splitting formulations.
standard math Convergence of Douglas-Rachford splitting and ADMM for convex closed proper functions, including l1 and indicator functions of affine sets
The algorithms inherit convergence from BC17, BPC+11, PB14, and FZB20; the paper does not prove convergence of the specific fixed-point maps.
domain assumption Random sprand instances with prescribed rank and singular-value profile are representative of real least-squares problems
All benchmarking conclusions rest on this distribution; no real-world data or cross-validation is used.

how reviews work

0 comments

Cite this review

Pith. "Pith review of On computing sparse universal solvers for key problems in statistics." pith.science (2026). https://pith.science/paper/DXIVVHOC

@misc{pith2026250904264,
  author       = {Pith},
  title        = {Pith review of: On computing sparse universal solvers for key problems in statistics},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/DXIVVHOC}},
  note         = {Machine review of arXiv:2509.04264}
}

read the original abstract

We give sparsity results and present algorithms for calculating minimum (vector) 1-norm universal solvers connected to least-squares problems. In particular, besides universal least-squares solvers, we consider minimum-rank universal least-squares solvers, and simultaneous universal minimum-norm/least-squares solvers. For all of these, we present and compare several new alternative linear-programming formulations and very effective proximal-point algorithms. Overall, we found that our new Douglas-Rachford splitting algorithms for these problems performed best.

Discussion (0). Continue with ORCID to comment.

Forward citations

Cited by 1 Pith paper

Reviewed papers in the Pith corpus that reference this work. Sorted by Pith novelty score. Full citation record

Sparse symmetric generalized inverses for sparse symmetric matrices
math.OC 2026-05 unverdicted novelty 6.0 of 10

A Douglas-Rachford splitting algorithm with closed-form projection computes substantially sparser symmetric generalized inverses than the Moore-Penrose pseudoinverse for sparse symmetric matrices.

Reference graph

Works this paper leans on

16 extracted references · 7 canonical work pages · cited by 1 Pith paper

[10]

https://doi.org/10.1007/s10898-021-01087-y . FZB20. Anqi Fu, Junzi Zhang, and Stephen Boyd. Anderson accelerated Douglas– Rachford splitting. SIAM Journal on Scientific Computing, 42(6):A3560–A3583, January

work page doi:10.1007/s10898-021-01087-y
[11]

https://doi.org/10.1137/19M1290097. PB14. Neal Parikh and Stephen Boyd. Proximal algorithms. Foundations and Trends in Optimization , 1(3):127–239,

doi:10.1137/19m1290097
[13]

https://arxiv.org/ abs/2401.17540. PFLX24b. Gabriel Ponte, Marcia Fampa, Jon Lee, and Luze Xu. On computing sparse generalized inverses. Operations Research Letters, 52:107058,

arXiv
[14]

https:// doi.org/10.1016/j.orl.2023.107058. 34 A.S. Machado, M. Fampa & J. Lee XFL20. Luze Xu, Marcia Fampa, and Jon Lee. 1-norm minimization and minimum- rank structured sparsity for symmetric and ah-symmetric generalized inverses: rank one and two,

work page arXiv 2023
[15]

To appear in: Fields Institute Communications vol- ume on Data Science and Optimization, S. Dang, A. Deza, S. Gupta, P.D. McNicholas, S. Pokutta, M. Sugiyama, eds. https://link.springer.com/book/ 9783032038432, https://arxiv.org/abs/2010.11406. XFLP21. Luze Xu, Marcia Fampa, Jon Lee, and Gabriel Ponte. Approximate 1-norm minimization and minimum-rank stru...

arXiv 2010
[16]

https://doi.org/10.1137/19M1281514

doi:10.1137/19m1281514
[1974]

https://link.springer.com/book/10.1007/b97366. BPC+11. Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein. Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations and Trends in Machine Learning , 3(1):1– 122,

doi:10.1007/b97366
[2011]

https://doi.org/10.1561/2200000016. DG17a. Ivan Dokmani´ c and R´ emi Gribonval. Beyond Moore-Penrose Part I: generalized inverses that minimize matrix norms,

doi:10.1561/2200000016

Show all 16 references

[2013]

https://doi.org/10.1109/ICASSP.2013.6638923. FFL16. Victor K. Fuentes, Marcia Fampa, and Jon Lee. Sparse pseudoinverses via LP and SDP relaxations of Moore-Penrose. In CLAIO 2016 (18th Latin-Iberian- American Conference on Operations Research) , pages 343–350,

2013
[2014]

https://doi.org/10.1561/2400000003. Pen55. Roger Penrose. A generalized inverse for matrices. Mathematical Proceedings of the Cambridge Philosophical Society , 51:406–413,

doi
[2016]

https: //marciafampa.com/pdf/CLAIO_2016_Proceedings_FuentesFampaLee.pdf. FFL20. Victor Fuentes, Marcia Fampa, and Jon Lee. Diving for sparse partially-reflexive generalized inverses. In Le Thi, H.A., et al., editor, WCGO 2019 , pages 89–98,

2019
[2017]

https://doi.org/10.1007/978-3-319-48311-5 . BG74. Adi Ben-Israel and Thomas N.E. Greville. Generalized Inverses: Theory and Ap- plications. Springer,

doi
[2018]

https://doi.org/10.1016/j.orl.2018. 09.005. FLP21. Marcia Fampa, Jon Lee, and Gabriel Ponte. Trading off 1-norm and sparsity against rank for linear models using mathematical optimization: 1-norm mini- mizing partially reflexive ah-symmetric generalized inverses. Open Journal ...

2018 doi
[2020]

https://doi.org/10.1007/978-3-030-21803-4_9 . FL18. Marcia Fampa and Jon Lee. On sparse reflexive generalized inverses. Operations Research Letters, 46(6):605–610,

doi
[2021]

https://doi.org/10.5802/ojmo.6. FLPX21. Marcia Fampa, Jon Lee, Gabriel Ponte, and Luze Xu. Experimental analysis of local searches for sparse reflexive generalized inverses. Journal of Global Opti- mization, 81:1057–1093,

doi
[2024]

https://arxiv.org/abs/2410.19506. DKV13. Ivan Dokmani´ c, Mihailo Kolundˇ zija, and Martin Vetterli. Beyond Moore- Penrose: sparse pseudoinverse. In ICASSP 2013 (38th International Confer- ence on Acoustics, Speech, and Signal Processing), pp. 6526–6530. IEEE,

2013 arXiv

Pith tools

Reviewed August 5, 2026 · model on record in the stance chip above.

[1] [10]

https://doi.org/10.1007/s10898-021-01087-y . FZB20. Anqi Fu, Junzi Zhang, and Stephen Boyd. Anderson accelerated Douglas– Rachford splitting. SIAM Journal on Scientific Computing, 42(6):A3560–A3583, January

work page doi:10.1007/s10898-021-01087-y

[2] [11]

https://doi.org/10.1137/19M1290097. PB14. Neal Parikh and Stephen Boyd. Proximal algorithms. Foundations and Trends in Optimization , 1(3):127–239,

doi:10.1137/19m1290097

[3] [13]

https://arxiv.org/ abs/2401.17540. PFLX24b. Gabriel Ponte, Marcia Fampa, Jon Lee, and Luze Xu. On computing sparse generalized inverses. Operations Research Letters, 52:107058,

arXiv

[4] [14]

https:// doi.org/10.1016/j.orl.2023.107058. 34 A.S. Machado, M. Fampa & J. Lee XFL20. Luze Xu, Marcia Fampa, and Jon Lee. 1-norm minimization and minimum- rank structured sparsity for symmetric and ah-symmetric generalized inverses: rank one and two,

work page arXiv 2023

[5] [15]

To appear in: Fields Institute Communications vol- ume on Data Science and Optimization, S. Dang, A. Deza, S. Gupta, P.D. McNicholas, S. Pokutta, M. Sugiyama, eds. https://link.springer.com/book/ 9783032038432, https://arxiv.org/abs/2010.11406. XFLP21. Luze Xu, Marcia Fampa, Jon Lee, and Gabriel Ponte. Approximate 1-norm minimization and minimum-rank stru...

arXiv 2010

[6] [16]

https://doi.org/10.1137/19M1281514

doi:10.1137/19m1281514

[7] [1974]

https://link.springer.com/book/10.1007/b97366. BPC+11. Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein. Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations and Trends in Machine Learning , 3(1):1– 122,

doi:10.1007/b97366

[8] [2011]

https://doi.org/10.1561/2200000016. DG17a. Ivan Dokmani´ c and R´ emi Gribonval. Beyond Moore-Penrose Part I: generalized inverses that minimize matrix norms,

doi:10.1561/2200000016