# Differences in A \\ b for sparse and nonsparse rank-deficient A
**URL:** https://discourse.julialang.org/t/differences-in-a-b-for-sparse-and-nonsparse-rank-deficient-a/66917
**Category:** General Usage
**Tags:** linearalgebra, sparse
**Created:** [August 24, 2021, 3:18pm UTC](https://discourse.julialang.org/t/differences-in-a-b-for-sparse-and-nonsparse-rank-deficient-a/66917 "2021-08-24T15:18:33Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@Thomas](https://discourse.julialang.org/u/Thomas)
#### Post date: [August 24, 2021, 3:18pm UTC](https://discourse.julialang.org/t/differences-in-a-b-for-sparse-and-nonsparse-rank-deficient-a/66917/1 "2021-08-24T15:18:33Z")
I have a question regarding the operator `\` for SpareseMatrixCSC. Thanks for your help.
```julia
using LinearAlgebra, SparseArrays
Is = [9, 8, 9, 12, 2, 5, 6, 11, 8, 12]
Js = [1, 3, 4, 4, 5, 5, 7, 7, 10, 10]
Vs = [0.046668312772398135, 0.0732172552606527, 0.07228019735547542, 0.27500506619596665, 0.8821260519923395, 0.188438924730004, 0.22883463110234215, 0.17705262933291666, 0.297278962545636, 0.07770904672896628]
A = sparse(Is, Js, Vs)
b = [0.0
0.26743848063303416
0.0
0.0
0.057129952809003536
0.20824640373847988
0.0
0.1580262484804872
0.09667497046894227
0.0
0.16112322314768995
0.27723584012581337]
qr(A) \ b ≈ A \ b # True
A \ b ≈ pinv(Matrix(A))*b # False
qr(Matrix(A), Val(true)) \ b ≈ pinv(Matrix(A))*b # True
```
Is this behavior expected?
- The first comparison is true because [https://github.com/JuliaLang/julia/blob/1b93d53fc4bb59350ada898038ed4de2994cce33/stdlib/SparseArrays/src/linalg.jl#L1551](https://github.com/JuliaLang/julia/blob/1b93d53fc4bb59350ada898038ed4de2994cce33/stdlib/SparseArrays/src/linalg.jl#L1551)
- The second comparison is questionable given the discussion (admittedly for dense matrices) [Efficient way of doing linear regression - #33 by stevengj](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/33)
- The last comparison makes sense because one needs a pivoted QR for rank-deficient case.
Perhaps there is a kwarg in `qr(::SparseMatrixCSC)` that I can use to get the same result as the dense fully pivoted case and without intermediate conversion of `A` to dense.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [August 24, 2021, 4:30pm UTC](https://discourse.julialang.org/t/differences-in-a-b-for-sparse-and-nonsparse-rank-deficient-a/66917/3 "2021-08-24T16:30:20Z")
> [@Thomas](#):
>
> Is this behavior expected?
You have a 12 \times 10 matrix A of rank 5, so the least-square problem \min\_x \Vert Ax - b\Vert\_2 (which is what `A \ b` and your other expressions solve) has infinitely many solutions.
For dense matrices, Julia’s `A \ b` (along with `qr(A) \ b` and `pinv(A) * b`) computes the _minimum-norm_ solution to the least-squares problem.
However, for sparse matrices, Julia’s (SuiteSparse’s) `A \ b` (and `qr(A) \ b`) computes the _sparsest_ solution, also called the “basic” solution (the solution with the most zero entries).
In your case, you’ll notice that `Matrix(A) \ b` has 4 zero entries and norm `≈ 1.57`, whereas `A \ b` has 5 zero entries and norm `≈ 2.62`.
This is documented if you look at the documentation for `\`, for example.
---
### Author:  [@Thomas](https://discourse.julialang.org/u/Thomas)
#### Post date: [August 24, 2021, 7:17pm UTC](https://discourse.julialang.org/t/differences-in-a-b-for-sparse-and-nonsparse-rank-deficient-a/66917/4 "2021-08-24T19:17:26Z")
Many thanks! I should have read the documentation _very_ carefully. Could you have a look at the logic of the following? I did some partial testing and it seems to return the correct solution.
```julia
# Compute the pseudoinverse solution for least squares problem
function pinvsol(A::SparseMatrixCSC, b::AbstractVector)
m, n = size(A)
F = qr(A) # SuiteSparse.SPQR.QRSparse
rnk = rank(F)
if rnk == n # full column rank
return qr(A) \ b
end
# rank-deficient
R11 = F.R[Base.OneTo(rnk),Base.OneTo(rnk)]
R12 = F.R[Base.OneTo(rnk),rnk+1:end]
S = sparse(R11 \ R12)
xb = R11 \ (F.Q[:,Base.OneTo(rnk)]'*b[F.prow])
p = n - rnk
x2 = [S; -sparse(I,p,p)] \ [xb; zeros(p)]
x1 = xb - S*x2
x = [x1; x2]
x[F.pcol] = x
return x # least squares least norm solution for sparse A
end
```
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [August 24, 2021, 8:00pm UTC](https://discourse.julialang.org/t/differences-in-a-b-for-sparse-and-nonsparse-rank-deficient-a/66917/5 "2021-08-24T20:00:47Z")
> [@Thomas](#):
>
> `F.Q[:,Base.OneTo(rnk)]`
Unfortunately, this creates a dense matrix. (`F.Q` is not stored explicitly, but as a composition of sparse or low-rank operations so that you can compute `F.Q` times a vector efficiently. Trying to pull out a subset of the columns of `Q` discards this structure.)
Instead, I would compute `(F.Q' * b[F.prow])[1:rnk]`, which should be equivalent and efficient without constructing a dense matrix.
> [@Thomas](#):
>
> `S = sparse(R11 \ R12)`
Note that this is also not a sparse matrix in general, since `inv(R11)` is in general not sparse. In general you want to compute the action of matrices like this on vectors rather than computing the matrix itself. But it may be okay for your purposes since maybe you are interested only in the case where the rank is very low, so a dense `R11` matrix inverse is acceptable?
I’ve just skimmed your code looking for common sparse-matrix gotchas; I haven’t really looked at the structure of your algorithm. Note that there are some papers on minimum-norm solutions with sparse matrices that may be helpful; for example, [Solution of Sparse Underdetermined Systems of Linear Equations (1984)](https://epubs.siam.org/doi/abs/10.1137/0905068?journalCode=sijcd4)
See also [Finding the least norm solution of least squares problem with a sparse matrix](https://discourse.julialang.org/t/finding-the-least-norm-solution-of-least-squares-problem-with-a-sparse-matrix/30702) for iterative methods.
---
### Author:  [@Thomas](https://discourse.julialang.org/u/Thomas)
#### Post date: [August 25, 2021, 12:28pm UTC](https://discourse.julialang.org/t/differences-in-a-b-for-sparse-and-nonsparse-rank-deficient-a/66917/6 "2021-08-25T12:28:05Z")
Many thanks for your help. The logic of this implementation follows from the book by Bjorck (section 2.7.4).
> [@stevengj](#):
>
> But it may be okay for your purposes since maybe you are interested only in the case where the rank is very low, so a dense `R11` matrix inverse is acceptable?
This is actually an overlook on my part. Generally I would avoid this if possible.
I had a look at the paper you sent and it appears the complete orthogonal decomposition is a good candidate to implement. Could you outline what are the “best” (simplicity vs speed) approaches for this problem that people have come up after the book of Bjorck?