# Result changes depending on whether matrix is sparse or not

**URL:** <https://discourse.julialang.org/t/result-changes-depending-on-whether-matrix-is-sparse-or-not/30669>\
**Category:** General Usage\
**Tags:** linearalgebra\
**Created:** [November 3, 2019, 11:01pm UTC](https://discourse.julialang.org/t/result-changes-depending-on-whether-matrix-is-sparse-or-not/30669 "2019-11-03T23:01:29Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![gmouts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gmouts/32/9896_2.png) [@gmouts](https://discourse.julialang.org/u/gmouts)\
**Post date:** [November 3, 2019, 11:01pm UTC](https://discourse.julialang.org/t/result-changes-depending-on-whether-matrix-is-sparse-or-not/30669/1 "2019-11-03T23:01:29Z")

</div>

I get this weird bug (I assume) but I cannot reproduce it with another matrix. So this is the setup:

```julia
A = sparse([1,2,1,3,3,4,5], [3,3,4,4,5,6,6,], [1,1,1,1,1,1,1], 6, 6 )

A = convert(SparseMatrixCSC{Float64,Int64},A)

d = vec(sum(A,dims=1))

M = spdiagm(0 => d)-transpose(A)

dropzeros!(M)

```

After this this happens:

```julia
julia> qr(M)\d
6-element Array{Float64,1}:
 -2.0
 -0.0
  0.0
  0.0
  1.0
  1.5

julia> qr(Matrix(M),Val(true))\d
6-element Array{Float64,1}:
 -1.4178403755868538
 -0.9389671361502343
 -0.17840375586854448
  0.20187793427230036
  0.8215962441314553
  1.5117370892018778

```

The second one is the correct result. Is this a bug? Am I doing something wrong here?

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [November 4, 2019, 2:40am UTC](https://discourse.julialang.org/t/result-changes-depending-on-whether-matrix-is-sparse-or-not/30669/3 "2019-11-04T02:40:37Z")

</div>

The matrix is singular. Using Ax=b terminology, the right-hand-side vector b is in the range of A, so solutions to Ax=b exist, but they are not unique. Different algorithms will produce different answers x, each of which satisfies Ax=b

```julia
julia> x2 = qr(Matrix(M),Val(true))\d
6-element Array{Float64,1}:
 -1.4178403755868538 
 -0.9389671361502343 
 -0.17840375586854448
  0.20187793427230036
  0.8215962441314553 
  1.5117370892018778 

julia> x1 = qr(M)\d
6-element Array{Float64,1}:
 -2.0
 -0.0
  0.0
  0.0
  1.0
  1.5

julia> M*x1
6-element Array{Float64,1}:
 0.0
 0.0
 2.0
 2.0
 1.0
 2.0

julia> M*x2
6-element Array{Float64,1}:
 0.0               
 0.0               
 1.9999999999999991
 1.9999999999999991
 0.9999999999999998
 2.0               

```

The first two rows of M are zero

```julia
6×6 Array{Float64,2}:
  0.0 0.0 0.0 0.0 0.0 0.0
  0.0 0.0 0.0 0.0 0.0 0.0
 -1.0 -1.0 2.0 0.0 0.0 0.0
 -1.0 0.0 -1.0 2.0 0.0 0.0
  0.0 0.0 -1.0 0.0 1.0 0.0
  0.0 0.0 0.0 -1.0 -1.0 2.0

```

so M has at best rank 4. The SVD shows that this is the case: two singular values are zero (to floating-point precision).

```julia
julia> (U,sigma,V) = svd(Matrix(M));

julia> sigma
6-element Array{Float64,1}:
 2.9775294768839364   
 2.5083437595312112   
 1.9623854744780629   
 0.9957776096426464   
 1.506335424828049e-16
 8.938809282980296e-18

```

---

<div class="post-metadata">

**Author:** ![gmouts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gmouts/32/9896_2.png) [@gmouts](https://discourse.julialang.org/u/gmouts)\
**Post date:** [November 4, 2019, 2:49am UTC](https://discourse.julialang.org/t/result-changes-depending-on-whether-matrix-is-sparse-or-not/30669/4 "2019-11-04T02:49:10Z")

</div>

Yes, the matrix is singular but the result of `qr(M)\d` should be the pseudo-inverse of `M` applied on `d`. The pseudo-inverse is unique and the result should be unique.

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [November 4, 2019, 3:04am UTC](https://discourse.julialang.org/t/result-changes-depending-on-whether-matrix-is-sparse-or-not/30669/5 "2019-11-04T03:04:44Z")

</div>

I wouldn’t make any stronger assumption than that QR backslash provides a minimizer of ||Ax-b||, and a non-unique minimizer for singular A. If you want the pseudo inverse, use `pinv`.

```julia
julia> pinv(Matrix(M))*d
6-element Array{Float64,1}:
 -1.4178403755868545 
 -0.9389671361502349 
 -0.17840375586854493
  0.2018779342723    
  0.8215962441314556 
  1.511737089201878  

```

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [November 4, 2019, 3:40am UTC](https://discourse.julialang.org/t/result-changes-depending-on-whether-matrix-is-sparse-or-not/30669/6 "2019-11-04T03:40:13Z")

</div>

Hmmm, my suggestion to use `pinv` is not a practical one. It’s never a good idea to compute an inverse to solve an Ax=b problem.

Is there any special reason in your application why one solution for x is better than another?

Regarding pseudo-inverse and QR, this thread is relevant.

> [@Moore-Penrose Generalized Inverse of Sparse Matrix](https://discourse.julialang.org/t/moore-penrose-generalized-inverse-of-sparse-matrix/17414/3):
>
> You almost never compute the inverse (pseudo or otherwise) of a sparse matrix because the inverse is generally dense. On the other hand, you do compute the application of the inverse to a vector, and you often precompute factorizations that let you apply the inverse more quickly. In the case of the ordinary inverse, you can apply a sparse inverse with A \ b and compute the factorization with lu(A) etcetera. In the case of the pseudoinverse, applying it computes the least-squares solution, and …

---

<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:** [November 4, 2019, 4:22am UTC](https://discourse.julialang.org/t/result-changes-depending-on-whether-matrix-is-sparse-or-not/30669/7 "2019-11-04T04:22:12Z")

</div>

Why do you think it is not a good idea?  
Computationally wise it might be better use direct methods to find the solution.  
But accuracy wise and certainly if one needs the Least Norm solution of the Least Squares there is nothing wrong about it.

Could it be that the Myth about “Don’t compute the inverse” is overblown?  
See [How Accurate is inv(A)\*b](https://arxiv.org/abs/1201.6035)?

By the way, MATLAB has the function [`lsqminnorm()`](https://www.mathworks.com/help/matlab/ref/lsqminnorm.html) which uses the complete orthogonal decomposition (COD) to compute the minimum norm least squares solution.  
If it is available in Julia one could use it as well as it should be faster than `pinv()` which uses the SVD.

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [November 4, 2019, 4:29am UTC](https://discourse.julialang.org/t/result-changes-depending-on-whether-matrix-is-sparse-or-not/30669/8 "2019-11-04T04:29:00Z")

</div>

SuiteSparseQR, which is used when A is sparse, returns a basic least-squares solutions, which means one with many zeros (hence the zeros that you observe). LAPACK, used in the dense case, returns the minimum-norm solution.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [November 4, 2019, 5:00am UTC](https://discourse.julialang.org/t/result-changes-depending-on-whether-matrix-is-sparse-or-not/30669/9 "2019-11-04T05:00:06Z")

</div>

@RoyiAvital it tends to be less about accuracy than speed. ([Don’t invert that matrix](https://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/)). For ill conditioned matrices, you very well may want to be using some sort of regularized method anyway ([Tikhonov regularization - Wikipedia](https://en.wikipedia.org/wiki/Tikhonov_regularization))

---

<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:** [November 4, 2019, 7:45am UTC](https://discourse.julialang.org/t/result-changes-depending-on-whether-matrix-is-sparse-or-not/30669/10 "2019-11-04T07:45:08Z")

</div>

> [@gmouts](#):
>
> Yes, the matrix is singular but the result of `qr(M)\d` should be the pseudo-inverse of `M` applied on `d` .

I am not sure about this in general (only for full rank matrices, and your matrix isn’t).

---

<div class="post-metadata">

**Author:** ![gmouts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gmouts/32/9896_2.png) [@gmouts](https://discourse.julialang.org/u/gmouts)\
**Post date:** [November 4, 2019, 9:49am UTC](https://discourse.julialang.org/t/result-changes-depending-on-whether-matrix-is-sparse-or-not/30669/11 "2019-11-04T09:49:37Z")

</div>

I do think that the result of `qr(M)\d` should be the minimum norm least squares solutions. If nothing else, just for consistency.

But in any case I have the following problem:

I have a singular sparse matrix `M` whose dimension can be very large and a vector `d` that may or may not be in the range of `M` and I want to compute `pinv(M)*d`. What is the most efficient way to do that?

---

<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:** [November 4, 2019, 9:54am UTC](https://discourse.julialang.org/t/result-changes-depending-on-whether-matrix-is-sparse-or-not/30669/12 "2019-11-04T09:54:07Z")

</div>

> [@Oscar\_Smith](#):
>
> @RoyiAvital it tends to be less about accuracy than speed. ([Don’t invert that matrix](https://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/)).

Great, so we agree that accuracy wise there is no advantage.  
Regarding speed, well, I don’t fully agree with the analysis you linked to.

He argues one can solve the linear system taking advantage of the special properties of the matrix. Yet this can be done for calculation of the inverse as well.  
Let’s take Tridiagonal Matrix for instance. You can solve directly using known methods. But you can also do the same tricks for calculating the inverse.

---

<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:** [November 4, 2019, 1:24pm UTC](https://discourse.julialang.org/t/result-changes-depending-on-whether-matrix-is-sparse-or-not/30669/13 "2019-11-04T13:24:54Z")

</div>

> [@RoyiAvital](#):
>
> Let’s take Tridiagonal Matrix for instance. You can solve directly using known methods. But you can also do the same tricks for calculating the inverse.

Solving a general N×N tridiagonal system Ax=b takes O(N) time and memory. Computing the inverse of A takes O(N²) time and memory because there are N² nonzero elements of the inverse to compute. This is better than O(N³), but still far worse than not computing the inverse at all.

For a _general_ dense matrix, it’s true that directly computing the inverse is not _terrible_, maybe a factor of 2 more time than just computing the LU factorization and within a small constant factor of the accuracy. But there’s usually no compelling reason to do it, and it’s arguably a bad habit because it will lead you into trouble for sparse and structured matrices.

---

<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:** [November 4, 2019, 2:17pm UTC](https://discourse.julialang.org/t/result-changes-depending-on-whether-matrix-is-sparse-or-not/30669/14 "2019-11-04T14:17:36Z")

</div>

> [@stevengj](#):
>
> > [@RoyiAvital](#):
> >
> > Let’s take Tridiagonal Matrix for instance. You can solve directly using known methods. But you can also do the same tricks for calculating the inverse.
> 
> Solving a general N×N tridiagonal system Ax=b takes O(N) time and memory. Computing the inverse of A takes O(N²) time and memory because there are N² nonzero elements of the inverse to compute. This is better than O(N³), but still far worse than not computing the inverse at all.
> 
> For a _general_ dense matrix, it’s true that directly computing the inverse is not _terrible_, maybe a factor of 2 more time than just computing the LU factorization and within a small constant factor of the accuracy. But there’s usually no compelling reason to do it, and it’s arguably a bad habit because it will lead you into trouble for sparse and structured matrices.

I totally agree with you.  
But in the case of the OP he had no choice but using the Pseudo Inverse.  
Then he was warned not to do it as usually most of us are warned in their first numerical course not to invert matrices for solving Linear Equations. Usually we are hinted it has something to do with accuracy and not only speed.

I was just pointing this is something that might be exaggerated.

---

<div class="post-metadata">

**Author:** ![gmouts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gmouts/32/9896_2.png) [@gmouts](https://discourse.julialang.org/u/gmouts)\
**Post date:** [November 4, 2019, 2:19pm UTC](https://discourse.julialang.org/t/result-changes-depending-on-whether-matrix-is-sparse-or-not/30669/15 "2019-11-04T14:19:30Z")

</div>

Oh, I just noticed your question. Yes, I need the solution with minimal norm. The matrix I use is the transpose of the Laplacian of a graph and it can be very large but also very sparse, so I’m convinced there has to be a better way of doing that.

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [November 4, 2019, 2:55pm UTC](https://discourse.julialang.org/t/result-changes-depending-on-whether-matrix-is-sparse-or-not/30669/16 "2019-11-04T14:55:07Z")

</div>

In general, you compute the minimum-norm least-squares solution using a complete orthogonal decomposition (COD). That’s what LAPACK does, but not SuiteSparseQR. It consists in computing a second QR factorization (of R’). There’s a Matlab package named Factorize that, among other things, calls SuiteSparseQR to compute a COD: [https://www.mathworks.com/matlabcentral/fileexchange/24119-don-t-let-that-inv-go-past-your-eyes-to-solve-that-system-factorize](https://www.mathworks.com/matlabcentral/fileexchange/24119-don-t-let-that-inv-go-past-your-eyes-to-solve-that-system-factorize). As far as I know, that hasn’t been ported to Julia.

---

<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:** [November 4, 2019, 4:22pm UTC](https://discourse.julialang.org/t/result-changes-depending-on-whether-matrix-is-sparse-or-not/30669/17 "2019-11-04T16:22:23Z")

</div>

> [@dpo](#):
>
> That’s what LAPACK does, but not SuiteSparseQR. It consists in computing a second QR factorization (of R’). There’s a Matlab package named Factorize that, among other things, calls SuiteSparseQR to compute a COD:

As I understand it, the latest algorithms to compute minimum-norm solutions using SuiteSparseQR are found in the [`spqr_rank` package](https://github.com/DrTimothyAldenDavis/SuiteSparse/tree/master/MATLAB_Tools/spqr_rank) described in [this paper](https://dl.acm.org/citation.cfm?id=2513116&picked=formats). In particular, section 2.7 of the paper says that the COD-based algorithm typically requires more time and memory (about a factor of 2 from fig. 6) than the newer algorithm. Both `spqr_rank` and the older `Factorize` package [are BSD-licensed](https://github.com/DrTimothyAldenDavis/SuiteSparse/blob/389d5df669fa7c48d14197c97272abb6155e3ee8/LICENSE.txt#L638-L673) by Tim Davis and are part of SuiteSparse. Would be nice to port `spqr_rank` to Julia. (See also [#20409](https://github.com/JuliaLang/julia/issues/20409).)

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [November 4, 2019, 4:51pm UTC](https://discourse.julialang.org/t/result-changes-depending-on-whether-matrix-is-sparse-or-not/30669/18 "2019-11-04T16:51:51Z")

</div>

Thanks. Wouldn’t it be awesome to have all that in pure Julia?