# Julia Linear System Solver (\`mldivide()\`)

**URL:** <https://discourse.julialang.org/t/julia-linear-system-solver-mldivide/33447>\
**Category:** Numerics\
**Tags:** linearalgebra\
**Created:** [January 16, 2020, 6:09pm UTC](https://discourse.julialang.org/t/julia-linear-system-solver-mldivide/33447 "2020-01-16T18:09:30Z")\
**Posts on this page:** 7\
**Page:** 1

<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 16, 2020, 6:09pm UTC](https://discourse.julialang.org/t/julia-linear-system-solver-mldivide/33447/1 "2020-01-16T18:09:30Z")

</div>

In the question [Solution of Underdetermined Systems Using LAPACK](https://discourse.julialang.org/t/33000) @stevengj [mentioned that the `\` operator in Julia returns the Minimum Norm Solution](https://discourse.julialang.org/t/solution-of-underdetermined-systems-using-lapack/33000/4) for [Under determined System](https://en.wikipedia.org/wiki/Underdetermined_system).

[The Matrix Division operator, `\` according to Julia’s documentation](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#Base.:%5C-Tuple%7BAbstractArray%7BT,2%7D%20where%20T,Union%7BAbstractArray%7BT,1%7D,%20AbstractArray%7BT,2%7D%7D%20where%20T%7D) uses, for this case, the Pivoted QR algorithm to return the Minimum Norm Solution.

In [MATLAB’s documentation for the algorithms](https://www.mathworks.com/help/matlab/ref/mldivide.html#bt4jslc-6) it says it uses also the QR solver, yet it doesn’t say what kind. Yet MATLAB doesn’t guarantee to return the Minimum Norm solution and indeed in most cases doesn’t.

My questions are:

1. Does Julia return the Minimum Norm solution since it uses the Pivoted QR or is there an additional step on top of that?
2. Performance wise, does it have any effect on the speed of the operation? Could it be faster if there was no requirements for the minimum norm solution?

---

<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 16, 2020, 6:43pm UTC](https://discourse.julialang.org/t/julia-linear-system-solver-mldivide/33447/2 "2020-01-16T18:43:28Z")

</div>

Since Matlab is closed-source, it’s difficult to know for certain what algorithm it is using. For non-square `x = A \ b` my understanding is that it calls [`lscov`](https://www.mathworks.com/help/matlab/ref/lscov.html), which uses QR but does _not_ find the minimum-norm solution. Instead, it finds a “basic” solution in which `x` is as _sparse_ as possible (by some algorithm that they don’t document, but probably is similar to the one [described here](http://www.netlib.org/lapack/lug/node42.html)).

Julia also uses QR, but in a different way (dating back to [#3315](https://github.com/JuliaLang/julia/pull/3315) in 2013): it implements the algorithm from LAPACK’s [`xgelsy`](http://www.netlib.org/lapack/explore-html/d7/d3b/group__double_g_esolve_ga385713b8bcdf85663ff9a45926fac423.html) to obtain the minimum-norm solution from the pivoted QR factorization.

So, basically, the difference lies mainly in the post-processing, i.e. how you _use_ the factorization. I expect that the costs are pretty similar, since in both cases most of the time is spent: (i) performing the QR factorization, (ii) doing a triangular solve, and (iii) multiplying by Q’. However, in the minimum-norm case there is some additional work in determining the rank and calling [`xtzrzf`](http://www.netlib.org/lapack/lapack-3.1.1/html/dtzrzf.f.html), so it is probably slower by a small factor.

If you have a “wide” matrix that is not rank-deficient and you want the minimum-norm solution, it will soon be faster (\> 2× for large `A`) to do `lq(A) \ b` ([#34350](https://github.com/JuliaLang/julia/pull/34350)). The xgelsy algorithm with pivoted QR is more general in that it handles rank-deficient `A`, however.

---

<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 16, 2020, 6:55pm UTC](https://discourse.julialang.org/t/julia-linear-system-solver-mldivide/33447/3 "2020-01-16T18:55:08Z")

</div>

(I have little doubt that Julia’s [pivoted-QR solver](https://github.com/JuliaLang/julia/blob/60b60f101e545bba849d639cbba34cd1c0b485cf/stdlib/LinearAlgebra/src/qr.jl#L763-L797) could be further optimized, though, since it hasn’t really changed since the first draft in 2013. e.g. it doesn’t have any special-case code for the common situation of a full-rank `A`, which could probably use a faster method.)

---

<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 17, 2020, 6:37am UTC](https://discourse.julialang.org/t/julia-linear-system-solver-mldivide/33447/4 "2020-01-17T06:37:15Z")

</div>

@stevengj, I appreciate your through answer.

If I want a solution, any solution, with the highest speed, which function should I use?

---

<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 17, 2020, 12:06pm UTC](https://discourse.julialang.org/t/julia-linear-system-solver-mldivide/33447/5 "2020-01-17T12:06:09Z")

</div>

> [@RoyiAvital](#):
>
> If I want a solution, any solution, with the highest speed, which function should I use?

For what kind of problem?

---

<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 17, 2020, 12:18pm UTC](https://discourse.julialang.org/t/julia-linear-system-solver-mldivide/33447/6 "2020-01-17T12:18:00Z")

</div>

I guess for:

1. `A` is Square, no other property.
2. `A` is not square, over determined. No other property.
3. `A` is not square, under determined. No other property.
4. `A` is not square, over determined, Full rank. No other property.

For the case `A` is under determined, Full rank. No other property you wrote the optimal will be `lq(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 17, 2020, 12:37pm UTC](https://discourse.julialang.org/t/julia-linear-system-solver-mldivide/33447/7 "2020-01-17T12:37:19Z")

</div>

If by “no other property” you mean that you don’t know whether it is rank-deficient, that is hard to answer in general. A typical default would be using pivoted QR, but you might need something fancier like an SVD-based solution to regularize your problem.

As I mentioned above, I’m guessing the minimum-norm solution from pivoted QR is slightly slower than getting a basic solution from pivoted QR, but we don’t currently have a function to give us a basic solution as far as I know, so you’d have to implement it yourself.

> `A` is not square, over determined, Full rank. No other property.

If it is far from rank-deficient, i.e. it is very well conditioned so that you are okay with losing twice as many digits to roundoff errors, you could use the normal equations as I mentioned in another approach. Otherwise probably `qr(A) \ b` (i.e. non-pivoted QR).

For very small matrices, of course, matters change because you don’t want the overhead of LAPACK. (I’m not sure where the crossover point occurs, but certainly \< 10x10 you want simple triple loops or even something unrolled like StaticArrays.)