# Solution of underdetermined systems using LAPACK

**URL:** <https://discourse.julialang.org/t/solution-of-underdetermined-systems-using-lapack/33000>\
**Category:** Performance\
**Created:** [January 5, 2020, 5:38pm UTC](https://discourse.julialang.org/t/solution-of-underdetermined-systems-using-lapack/33000 "2020-01-05T17:38:47Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![lucas711642](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lucas711642/32/12051_2.png) [@lucas711642](https://discourse.julialang.org/u/lucas711642)\
**Post date:** [January 5, 2020, 5:38pm UTC](https://discourse.julialang.org/t/solution-of-underdetermined-systems-using-lapack/33000/1 "2020-01-05T17:38:47Z")

</div>

I need to get any solution of a bunch of underdetermined systems `A*x = b` fastly and cheaply. My first thought is to compute a minimum norm solution using the LQ decomposition of `A` (which can be used as workspace to minimize allocations) with something like this:

```julia
using LinearAlgebra
m = 5; n = 10
A = randn(m, n)
b = randn(m)
x = similar(A, n)
ldiv!(x, lq!(A), b)

```

but this results in a `DimensionMismatch` error, which I think is due to the square matrix `Q` of the LQ decomposition.

Since I want to make this operation as cheap as possible, I considered calling `LAPACK.gels!` directly, but the documentation says that `b` is updated with the solution. How does this work in the minimum norm case, when the size of `b` is less than the size of the solution `x`?

Right now my best clue is to iterate over the LQ decomposition of `A` (allocating more memory than necessary) and perform the matrix products manually:

```julia
using LinearAlgebra
m = 5; n = 10
A = randn(m, n)
b = randn(m)
L, Q = lq!(A)
ldiv!(LowerTriangular(L), b)
Q'b

```

Note: using a preallocated `x` and trying `mul!(x, Q', b)` also produces an `DimensionMismatch` error.

Is there any better way than this?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 5, 2020, 6:38pm UTC](https://discourse.julialang.org/t/solution-of-underdetermined-systems-using-lapack/33000/2 "2020-01-05T18:38:55Z")

</div>

> [@lucas711642](#):
>
> solution of a bunch of underdetermined systems

Since they have an infinite number of solutions (the subspace), do you need the one with the minimum norm, or would an arbitrary one do?

---

<div class="post-metadata">

**Author:** ![lucas711642](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lucas711642/32/12051_2.png) [@lucas711642](https://discourse.julialang.org/u/lucas711642)\
**Post date:** [January 5, 2020, 6:58pm UTC](https://discourse.julialang.org/t/solution-of-underdetermined-systems-using-lapack/33000/3 "2020-01-05T18:58:10Z")

</div>

It can be any solution. I chose the minimum norm for convenience with the LAPACK functions.

I assume `A` has full rank, but I don’t have which columns of `A` form a nonsingular submatrix of `A`.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 5, 2020, 7:10pm UTC](https://discourse.julialang.org/t/solution-of-underdetermined-systems-using-lapack/33000/4 "2020-01-05T19:10:14Z")

</div>

> [@lucas711642](#):
>
> I need to get any solution of a bunch of underdetermined systems `A*x = b` fastly and cheaply.

Just do `x = A \ b` — it gives you the minimum-norm solution by default when `A` has more columns than rows. If you need to solve a sequence of right-hand sides for the same `A`, use the pivoted QR factorization:

```julia
QR = qr(A, Val(true))
x = QR \ b

```

If you [search discourse for “minimum norm”](https://discourse.julialang.org/search?q=minimum%20norm) you will find more discussions of these topics.

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [January 6, 2020, 6:11pm UTC](https://discourse.julialang.org/t/solution-of-underdetermined-systems-using-lapack/33000/5 "2020-01-06T18:11:57Z")

</div>

> [@stevengj](#):
>
> Just do `x = A \ b` — it gives you the minimum-norm solution by default when `A` has more columns than rows.

It is nice to see that.  
I was under the impression Julia’s `\` is identical to MATLAB.  
Yet MATLAB’s solution isn’t guaranteed to return the Least Norm solution for Under Determined Systems.  
I look on the documentation, both use `QR` yet Julia states it uses Pivoted QR. Is that the trick for that?  
Does it cost more (Performance wise) than regular `QR`?

---

<div class="post-metadata">

**Author:** ![lucas711642](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lucas711642/32/12051_2.png) [@lucas711642](https://discourse.julialang.org/u/lucas711642)\
**Post date:** [January 10, 2020, 10:45pm UTC](https://discourse.julialang.org/t/solution-of-underdetermined-systems-using-lapack/33000/6 "2020-01-10T22:45:20Z")

</div>

I was looking for a solution that minimizes memory allocation because I’m solving several systems in a loop. Since my systems have full rank, the cheaper LQ factorization can achieve the same result as pivoted QR.

I was wondering if there is any function available which can use the LQ factorization itself, without constructing the L and Q matrices explicitly.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 10, 2020, 11:36pm UTC](https://discourse.julialang.org/t/solution-of-underdetermined-systems-using-lapack/33000/7 "2020-01-10T23:36:10Z")

</div>

> [@lucas711642](#):
>
> but this results in a `DimensionMismatch` error, which I think is due to the square matrix `Q` of the LQ decomposition.

Seems like it is just a bug / missing feature. `lq(A) \ b` should be perfectly well defined.

> <https://github.com/JuliaLang/julia/issues/34348>
>
> As described in the \[LAPACK documentation\](https://www.netlib.org/lapack/lug/nod…e41.html), the LQ factorization can be used to compute a minimum-norm solution for (full-rank) underdetermined systems (more columns than rows). However \`lq(A) \\ b\` currently calls \`checksquare\`, and the \`ldiv!\` routine also doesn't handle the case where the output solution has more rows than the input right-hand-side (\[discourse\](https://discourse.julialang.org/t/solution-of-underdetermined-systems-using-lapack/33000)).
> 
> Should be straightforward to update the code to eliminate this restriction, since it doesn't require any new solver.

---

<div class="post-metadata">

**Author:** ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)\
**Post date:** [January 13, 2020, 2:28pm UTC](https://discourse.julialang.org/t/solution-of-underdetermined-systems-using-lapack/33000/8 "2020-01-13T14:28:32Z")

</div>

The dense QR and LQ solvers in LinearAlgebra all seem to materialize a new triangular matrix, against Lucas’ wish to avoid allocation. IIUC the corresponding LAPACK routines use views into the packed factorizations, so why don’t the Julia versions?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 13, 2020, 3:52pm UTC](https://discourse.julialang.org/t/solution-of-underdetermined-systems-using-lapack/33000/9 "2020-01-13T15:52:49Z")

</div>

> [@Ralph\_Smith](#):
>
> IIUC the corresponding LAPACK routines use views into the packed factorizations, so why don’t the Julia versions?

The [new LQ solver](https://github.com/JuliaLang/julia/pull/34350) works in-place when you call `ldiv!`.

The pivoted-QR solver needs to [call `tzrzf!`](https://github.com/JuliaLang/julia/blob/ade471359849cb2767eed69921892ab4280c3b81/stdlib/LinearAlgebra/src/qr.jl#L791), which requires a new array. (Maybe the subsequent line constructing an `UpperTriangular` matrix could use a view, though? Better yet, maybe this factorization could be cached in the `QRPivoted` for subsequent calls?)

(The pivoted-QR `ldiv!` implementation is fairly old — it derives from the `/` function in [#3315](https://github.com/JuliaLang/julia/pull/3315) based on LAPACK’s [xgelsy](https://www.netlib.org/lapack/lug/node27.html), which was ported to an “in-place” function in [#5526](https://github.com/JuliaLang/julia/pull/5526) without changing the fact that it allocated new matrices. It would probably be useful to revisit it, cc @andreasnoack.)

---

<div class="post-metadata">

**Author:** ![andreasnoack](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasnoack/32/27_2.png) [@andreasnoack](https://discourse.julialang.org/u/andreasnoack)\
**Post date:** [January 14, 2020, 9:19am UTC](https://discourse.julialang.org/t/solution-of-underdetermined-systems-using-lapack/33000/10 "2020-01-14T09:19:12Z")

</div>

It’s been a while but I guess I couldn’t figure out a way to avoid allocations. I think the LAPACK version uses work space for this step. However, it might very well be possible to reduce the number of allocations here.