Solving sparse underdetermined linear systems

Fortranners,

I’m wondering if anyone knows of a good Fortran library for solving a sparse underdetermined system of equations in the minimum-norm sense?

That is, solve Ax = b, where A is m x n, with m<n (so, more columns than rows). Return the solution that minimizes the norm of x.

I can solve these equations with the LAPACK routine dgels. That works very well, but is assuming a dense system. My system is large (say n around 30,000) and sparse (only a small precent of the elements are non-zero, mostly along the diagonal), and this routine is quite slow for those cases.

I’ve been experimenting with lsqr (sometimes works, sometimes doesn’t converge) and lusol (can’t seem to get to work, claims systems are singular that dgels solves just fine). I’m not a linear algebra guru though.

What else is out there?

1 Like

Maybe one of these could be of use?

2 Likes

Not a direct answer, but I just saw a talk about something somewhat related:

Might be some inspiration in there.

2 Likes

Oooo… qr_mumps looks interesting! I’ll try it out. Unfortunately looks like it has a complicated build process involving several other libs that I also have to figure out how to build and install. (like the pre-FPM bad old days).

the MKL one says " NOTE: The underdetermined systems of equations are not supported"

I’ve never solved such problems with sparse matrices, however it does look like SciPy provides both LSQR and LSMR.

The LSMR paper shows it’s based on the GMRES method, so I would bet it’s more robust than conjugate gradient for general real-world cases. It’s page is here, there are many implementations including Fortran, may well be worth a try.

3 Likes

There is the R package sparsenet: Fit Sparse Linear Regression Models via Nonconvex Optimization by Mazumder, Friedman and Hastie, with underlying Fortran (originally mortran) code sparsenet.f, with the algorithm explained in a paper.