# Solving sparse underdetermined linear systems

**URL:** <https://fortran-lang.discourse.group/t/solving-sparse-underdetermined-linear-systems/6472>\
**Category:** Help\
**Created:** [September 5, 2023, 9:08pm UTC](https://fortran-lang.discourse.group/t/solving-sparse-underdetermined-linear-systems/6472 "2023-09-05T21:08:13Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\
**Post date:** [September 5, 2023, 9:08pm UTC](https://fortran-lang.discourse.group/t/solving-sparse-underdetermined-linear-systems/6472/1 "2023-09-05T21:08:13Z")

</div>

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](https://netlib.org/lapack/explore-html/d7/d3b/group__double_g_esolve_ga225c8efde208eaf246882df48e590eac.html). 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](https://github.com/jacobwilliams/LSQR) (sometimes works, sometimes doesn’t converge) and [lusol](https://github.com/nwh/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?

---

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [September 5, 2023, 9:17pm UTC](https://fortran-lang.discourse.group/t/solving-sparse-underdetermined-linear-systems/6472/2 "2023-09-05T21:17:23Z")

</div>

Maybe one of these could be of use?

- [qr\_mumps / qr\_mumps · GitLab](https://gitlab.com/qr_mumps/qr_mumps)
- [https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2023-2/sparse-qr-routines.html](https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2023-2/sparse-qr-routines.html)

---

<div class="post-metadata">

**Author:** ![everythingfunctional](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/everythingfunctional/32/176_2.png) [@everythingfunctional](https://fortran-lang.discourse.group/u/everythingfunctional)\
**Post date:** [September 5, 2023, 9:18pm UTC](https://fortran-lang.discourse.group/t/solving-sparse-underdetermined-linear-systems/6472/3 "2023-09-05T21:18:51Z")

</div>

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

> **[Parallel Cholesky Factorization for Banded Matrices using OpenMP Tasks](https://arxiv.org/abs/2305.04635)**
>
> Cholesky factorization is a widely used method for solving linear systems
> involving symmetric, positive-definite matrices, and can be an attractive
> choice in applications where a high degree of numerical stability is needed.
> One such application is...

Might be some inspiration in there.

---

<div class="post-metadata">

**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\
**Post date:** [September 5, 2023, 11:06pm UTC](https://fortran-lang.discourse.group/t/solving-sparse-underdetermined-linear-systems/6472/4 "2023-09-05T23:06:33Z")

</div>

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"

---

<div class="post-metadata">

**Author:** ![FedericoPerini](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/federicoperini/32/1750_2.png) [@FedericoPerini](https://fortran-lang.discourse.group/u/FedericoPerini)\
**Post date:** [September 6, 2023, 6:41am UTC](https://fortran-lang.discourse.group/t/solving-sparse-underdetermined-linear-systems/6472/5 "2023-09-06T06:41:47Z")

</div>

I’ve never solved such problems with sparse matrices, however it does look like SciPy provides both `LSQR` and [`LSMR`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.lsmr.html).

The [`LSMR` paper](https://arxiv.org/pdf/1006.0758.pdf) 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](https://web.stanford.edu/group/SOL/software/lsmr/), there are many implementations including Fortran, may well be worth a try.

---

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [September 13, 2023, 5:26pm UTC](https://fortran-lang.discourse.group/t/solving-sparse-underdetermined-linear-systems/6472/6 "2023-09-13T17:26:50Z")

</div>

There is the R package [sparsenet: Fit Sparse Linear Regression Models via Nonconvex Optimization](https://cran.r-project.org/web/packages/sparsenet/index.html) by Mazumder, Friedman and Hastie, with underlying Fortran (originally [mortran](https://www.slac.stanford.edu/exp/e871/documentation/offline/mortran.html)) code [sparsenet.f](https://github.com/cran/sparsenet/blob/master/src/sparsenet.f), with the algorithm explained in a [paper](https://hastie.su.domains/public/Papers/Sparsenet/Mazumder-SparseNetCoordinateDescent-2011.pdf).

---

<div class="post-metadata">

**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\
**Post date:** [February 19, 2024, 10:13pm UTC](https://fortran-lang.discourse.group/t/solving-sparse-underdetermined-linear-systems/6472/7 "2024-02-19T22:13:24Z")

</div>

Follow up: it does appear the `qr_mumps` is exactly what I want. But I can’t get it to work, at least on my Mac laptop. Opened a ticket here: [Apple Silicon Segmentation fault (#2) · Issues · qr\_mumps / qr\_mumps · GitLab](https://gitlab.com/qr_mumps/qr_mumps/-/issues/2)

I’m curious to know if anyone here has used this library?

---

<div class="post-metadata">

**Author:** ![rwmsu](https://avatars.discourse-cdn.com/v4/letter/r/48db29/32.png) [@rwmsu](https://fortran-lang.discourse.group/u/rwmsu)\
**Post date:** [February 21, 2024, 2:52pm UTC](https://fortran-lang.discourse.group/t/solving-sparse-underdetermined-linear-systems/6472/8 "2024-02-21T14:52:49Z")

</div>

@jacobwilliams, just a wild guess based on looking over the qr\_mumps and StarPU documentation but you might try installing hwloc as suggested in the StarPU docs if you don’t already have it installed. Your error messages suggest its something to do with how StarPU schedules work on processors.

See

[https://www.open-mpi.org/projects/hwloc/](https://www.open-mpi.org/projects/hwloc/)

Also

StarPU runtime is here

> **[GitHub - starpu-runtime/starpu: This is a mirror of...](https://github.com/starpu-runtime/starpu)**
>
> This is a mirror of https://gitlab.inria.fr/starpu/starpu where our development happens, but contributions are welcome here too! - starpu-runtime/starpu

You might also look at the SVD package PROPACK. You would have to write your own code to solve the LS problem using SVD generated by PROPACK but thats just a few matrix multiplies. I recall reading somewhere that some folks suggested using SVD for underdetermined systems (I think thats what LAPACK does). The problem with SVD though is its high computational cost. It looks like PROPACK uses openMP so it might not be so bad.

PROPACK link:

[http://sun.stanford.edu/~rmunk/PROPACK/](http://sun.stanford.edu/~rmunk/PROPACK/)

---

<div class="post-metadata">

**Author:** ![MarDie](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mardie/32/54_2.png) [@MarDie](https://fortran-lang.discourse.group/u/MarDie)\
**Post date:** [February 21, 2024, 6:03pm UTC](https://fortran-lang.discourse.group/t/solving-sparse-underdetermined-linear-systems/6472/9 "2024-02-21T18:03:26Z")

</div>

I maintain the Arch Linux package and have it installed. I can test your code on x86\_64 if it is easy to install

---

<div class="post-metadata">

**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\
**Post date:** [February 25, 2024, 9:34pm UTC](https://fortran-lang.discourse.group/t/solving-sparse-underdetermined-linear-systems/6472/10 "2024-02-25T21:34:50Z")

</div>

Oh cool!

Yeah, if you are bored one day and want to try this example, that would be great!:

> **[GitHub - jacobwilliams/minnorm\_test at qrm](https://github.com/jacobwilliams/minnorm_test/tree/qrm)**
>
> qrm

(note it’s the `qrm` branch). All this is doing is reading a sparse matrix (about 24000x28000) from a data file and computing a minimum norm solution. first with Lapack (`dgels`) and then with `qr_mumps`. The `qr_mumps` one is totally untested, but I based it on one of the examples, so hopefully it will work. the Lapack one takes a few minutes on my laptop.

You can compile with FPM: `fpm run --profile release --flag "-lblas -llapack"` but would need to add whatever extra flags you need to get `qr_mumps` to link.

Would really be interested in seeing the comparison.

---

<div class="post-metadata">

**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\
**Post date:** [February 25, 2024, 9:36pm UTC](https://fortran-lang.discourse.group/t/solving-sparse-underdetermined-linear-systems/6472/11 "2024-02-25T21:36:55Z")

</div>

Yeah I have StarPU compiled and hwloc installed. I’m not sure if any of the developers of any of these libraries have ever tried it on an Apple Silicon Mac, so it’s hard for me to say what the problem could be.

Thanks, I’ll look at PROPACK. But glancing at one of the pdfs, its seems to say that it requires `m>=n`, so that doesn’t work for my case.

---

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [February 25, 2024, 10:04pm UTC](https://fortran-lang.discourse.group/t/solving-sparse-underdetermined-linear-systems/6472/12 "2024-02-25T22:04:41Z")

</div>

Maybe you can find something of use in the book - [Sparse Solutions of Underdetermined Linear Systems and Their Applications, 2021, SIAM](https://epubs.siam.org/doi/book/10.1137/1.9781611976519). It’s supposed to contain 64 algorithms (in MATLAB). I only had a glance at the front matter (which is openly available).

---

<div class="post-metadata">

**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\
**Post date:** [February 25, 2024, 11:15pm UTC](https://fortran-lang.discourse.group/t/solving-sparse-underdetermined-linear-systems/6472/13 "2024-02-25T23:15:16Z")

</div>

Interesting, I’ve never seen this book before.

They seen to be solving something slightly different though based on this from the intro:

> The aim of the study is to find a sparse solution in the sense that the number of nonzero entries of the solution is the smallest.

Which isn’t really the same thing as a minimum norm solution of a sparse system. But I’ll look at it. Thanks!

---

<div class="post-metadata">

**Author:** ![rwmsu](https://avatars.discourse-cdn.com/v4/letter/r/48db29/32.png) [@rwmsu](https://fortran-lang.discourse.group/u/rwmsu)\
**Post date:** [February 25, 2024, 11:32pm UTC](https://fortran-lang.discourse.group/t/solving-sparse-underdetermined-linear-systems/6472/14 "2024-02-25T23:32:51Z")

</div>

Jacob, if you dig deep into the hwloc documention they mention it supporting Apple M1. There is an xml file that has what I guess is a “profile” for an M1 along with several others. Again, the error message you posted makes me think you problem is with StarPU and not qr\_mumps. I didn’t look at the qr\_mumps build files but you might have to build StarPU and hwloc outside of qr\_mumps first. Also, I was looking at the source code and didn’t see any comment about restriction to overdetermined or square matrices. Still might be worth a try just to see if it works. Also, PROPACK looks like a good candidate for modernizing and adding FPM support. Might give that a go and see where it leads.

---

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [February 25, 2024, 11:34pm UTC](https://fortran-lang.discourse.group/t/solving-sparse-underdetermined-linear-systems/6472/15 "2024-02-25T23:34:33Z")

</div>

Oh right, the preface repeats this,

> linear system of equations consists of a known matrix A and a known vector b such that Ax=b for an unknown vector x, where A is of size m×n. When m=n and m \> n, methods of solving for x are well known and are discussed in standard numerical analysis textbooks. However, when m \< n, Ax=b is called an underdetermined linear system. Such a linear system has become a subject of research in the last 15 years as part of an extremely active study in the community of compressive sensing. Mathematically, an underdetermined linear system has many solutions. _The aim of the study is to find a sparse solution in the sense that the number of nonzero entries of the solution is the smallest._ [emphasis added] Although there exists a sparse solution, the difficulty is that it may take years to find such an exact solution using our current computer power.Therefore, we must look for an approximation of the sparse solution which can be found within reasonable time and tolerance.

An under-determined system has an infinite number of solutions. Judging by the description of [compressed sensing](https://en.wikipedia.org/wiki/Compressed_sensing) on Wikipedia, a solution with a small number of non-zeros is different type of constraint from the minimization of the L\_2-norm (the least squares solution), and is more suitable for some types of problems.

---

<div class="post-metadata">

**Author:** ![JohnCampbell](https://avatars.discourse-cdn.com/v4/letter/j/5daacb/32.png) [@JohnCampbell](https://fortran-lang.discourse.group/u/JohnCampbell)\
**Post date:** [February 26, 2024, 1:34am UTC](https://fortran-lang.discourse.group/t/solving-sparse-underdetermined-linear-systems/6472/16 "2024-02-26T01:34:21Z")

</div>

Could linear programming strategies be adapted for optimised solutions to this problem ?

---

<div class="post-metadata">

**Author:** ![MarDie](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mardie/32/54_2.png) [@MarDie](https://fortran-lang.discourse.group/u/MarDie)\
**Post date:** [February 29, 2024, 12:26pm UTC](https://fortran-lang.discourse.group/t/solving-sparse-underdetermined-linear-systems/6472/17 "2024-02-29T12:26:36Z")

</div>

took a little bit longer because I mixed up MUMPS and QR\_MUMPS. But installing QR\_MUMPS is very straight forward:

```sh
cmake -B build -DARITH=d -DCMAKE_INSTALL_PREFIX=${PWD} -DBUILD_SHARED_LIBS=ON -DQRM_ORDERING_METIS=ON
cmake --build build --parallel --target install  

```

I have installed METIS as a package. After exporting the `LD_LIBRARY_PATH`,

```auto
fpm run --profile release --flag "-lblas -llapack -L/tmp/qr_mumps/lib -ldqrm -lqrm_common -I/tmp/qr_mumps/include"

```

I get

```sh
Project is up to date
 reading data/matrix_test.json
 check_solution_in_file
 b error: 6.3514977426607910E-013
 solve_with_qrm
 error: 4.6937255843139300E-012
 time: 0.25215055500000000 sec
 solve_with_lapack
 allocate dense matrix
 call dgels
 info = 0
 error: 4.1530253854489664E-012
 time: 634.22027651799999 sec

```

---

<div class="post-metadata">

**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\
**Post date:** [February 29, 2024, 3:22pm UTC](https://fortran-lang.discourse.group/t/solving-sparse-underdetermined-linear-systems/6472/18 "2024-02-29T15:22:18Z")

</div>

0.2 sec! Wow, thanks! I guess I’m going to try Linux. How many cores does the system you ran this on have?

Also did you not have to build/install StarPU?

---

<div class="post-metadata">

**Author:** ![MarDie](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mardie/32/54_2.png) [@MarDie](https://fortran-lang.discourse.group/u/MarDie)\
**Post date:** [February 29, 2024, 3:34pm UTC](https://fortran-lang.discourse.group/t/solving-sparse-underdetermined-linear-systems/6472/19 "2024-02-29T15:34:18Z")

</div>

I did not install StarPU and use a Intel® Core™ i7-8665U (4 cores with hyperthreading). I could not detect any traces of parallel execution, `/usr/bin/time` reports 99% use for the whole fpm comman.d

---

<div class="post-metadata">

**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\
**Post date:** [February 29, 2024, 3:43pm UTC](https://fortran-lang.discourse.group/t/solving-sparse-underdetermined-linear-systems/6472/20 "2024-02-29T15:43:14Z")

</div>

Oh interesting! So you think this was just the serial execution? So the speed up wrt lapack is just from using the sparsity, rather than doing anything in parallel?

[Next page](https://fortran-lang.discourse.group/t/solving-sparse-underdetermined-linear-systems/6472.md?page=2)