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.
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.