# Computing sparse orthogonal projections
**URL:** https://discourse.julialang.org/t/computing-sparse-orthogonal-projections/62652
**Category:** Numerics
**Tags:** linearalgebra
**Created:** [June 9, 2021, 6:11pm UTC](https://discourse.julialang.org/t/computing-sparse-orthogonal-projections/62652 "2021-06-09T18:11:07Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@Thomas](https://discourse.julialang.org/u/Thomas)
#### Post date: [June 9, 2021, 6:11pm UTC](https://discourse.julialang.org/t/computing-sparse-orthogonal-projections/62652/1 "2021-06-09T18:11:07Z")
Follow up question: since this is part of my implementation of the projection matrix `A*inv(A'A)*A'` do you see a fast way to stay within SparseMatrixCSC?
Working implementation but not sparse
```julia
project(A::SparseMatrixCSC,v::SparseVectorCSC) = A*((A'*A)\(A'*Vector(v))
```
I can think of two methods:
1. The SparseMatrix library only has `lu,qr` which can take sparse input and so I think I have to implement the backslash operator for sparse inputs using these.
2. Use sparse `svd` (extra package) for `A` and reconstruct the `inv(A'A)` to stay within sparse data format.
Edit: fixed typo in `project`.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [June 10, 2021, 1:22am UTC](https://discourse.julialang.org/t/computing-sparse-orthogonal-projections/62652/2 "2021-06-10T01:22:39Z")
> [@Thomas](#):
>
> Follow up question: since this is part of my implementation of the projection matrix `A*inv(A'A)*A'` do you see a fast way to stay within SparseMatrixCSC?
Almost **never compute matrix inverses** , especially with sparse matrices. For something like this the most direct approach would be `A * ((A'A) \ (A'x))` to multiply it by a vector `x`. If you want to precompute the LU factorization you could do `LU = lu(A'A)` and then compute `A * (LU \ (A'x))`. In some cases (if you have many more rows than columns), `A'A` may be dense, in which case you might want to use `LU = lu(Matrix(A'A))` instead.
**However** , whether you are using sparse or dense matrices, you should strongly consider using the [QR factorization](https://en.wikipedia.org/wiki/QR_decomposition) for this kind of projection, in order to avoid explicitly computing A^\* A and squaring the condition number (see also the discussion [here](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/33)). In particular, if the m\times n (n \le m) full-rank matrix A=\hat{Q}\hat{R} where \hat{Q} is m\times n, \hat{Q}^\* \hat{Q} = I, and \hat{R} is square (n\times n) and invertible (i.e. the “thin” QR factorization), then A(A^\* A)^{-1} A^\* = \hat{Q}\hat{R} (\hat{R}^\* \hat{R})^{-1} \hat{R}^\* \hat{Q}^\* = \boxed{\hat{Q}\hat{Q}^\*}.
The `QR = qr(A)` function in Julia’s `LinearAlgebra` package computes the QR factorization, and `SparseArrays` provides an efficient method for sparse `A` as well. `QR.Q` represents the m\times m matrix Q of the “full” QR factorization, of which \hat{Q} is the first n columns — the matrix isn’t explicitly stored, but there are fast methods provided to multiply by Q and Q^\*.
Thus, to obtain p = \hat{Q}\hat{Q}^\* x for a vector x and a matrix A, do:
```julia
QR = qr(A)
y = QR.Q' * x
y[size(A,2)+1:end] .= 0 # project to y = Q̂'x
p = QR.Q * y
```
Unfortunately, it is somewhat more complicated in the sparse case, because in that case Julia computes a factorization P\_r A P\_c^\* = QR \Longleftrightarrow A = P\_r^{\*} Q R P\_c of a _permuted_ version of the matrix A (where P\_r and P\_c are permutation matrices for the rows and columns of A, respectively). Hence the projection becomes A(A^\* A)^{-1} A^\* = \cdots = \boxed{P\_r^{\*} \hat{Q} \hat{Q}^\* P\_r}, which is computed by
```julia
QR = qr(A)
y = QR.Q' * x[QR.prow] # compute y = Q'Pᵣx
y[size(A,2)+1:end] .= 0 # project to y = Q̂'Pᵣx
p = QR.Q * y
p[QR.prow] = p # apply Pᵣᵀ permutation
```
**Alternatively** , realize that the intermediate step (A^\*A)^{-1} A^\* x is simply the **least-squares solution**. So, your projection can also be computed by:
```julia
p = A * (A \ x)
```
since `A \ x` does least-squares (via QR) for a matrix with more rows than columns. Or, using a precomputed QR factorization:
```julia
QR = qr(A)
p = A * (QR \ x)
```
This is closely related to the solution above involving Q, but may be slightly less efficient (and less accurate?) because it involves a (cheap) additional triangular solve step using the matrix \hat{R} (though this may be offset by the reduced cost of multiplying by `A` rather than `Q`).
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [June 10, 2021, 2:20am UTC](https://discourse.julialang.org/t/computing-sparse-orthogonal-projections/62652/3 "2021-06-10T02:20:26Z")
> [@Thomas](#):
>
> Working implementation but not sparse
>
> ```julia
> project(A::SparseMatrixCSC,v::SparseVectorCSC) = A*(inv(A'A)\(A'*Vector(v))
>
> ```
Note that this code does not seem to be correct. You should _either_ use `inv` or use `\`, but not both.
> [@Thomas](#):
>
> The SparseMatrix library only has `lu,qr` which can take sparse input and so I think I have to implement the backslash operator for sparse inputs using these.
`SparseArrays` already provides sparse `\` methods, including `\` for factorization objects.
---
### Author:  [@Thomas](https://discourse.julialang.org/u/Thomas)
#### Post date: [June 10, 2021, 8:43am UTC](https://discourse.julialang.org/t/computing-sparse-orthogonal-projections/62652/4 "2021-06-10T08:43:38Z")
Sorry that was a typo; I fixed it in the question post.
Do you have this error?
```julia
using SparseArrays, LinearAlgebra
A = sprand(10,10,0.5)
b = sprand(10,0.5)
A\b
ERROR: MethodError: no method matching ldiv!(::SuiteSparse.UMFPACK.UmfpackLU{Float64, Int64}, ::SparseVector{Float64, Int64})
versioninfo()
Julia Version 1.6.1
Commit 6aaedecc44 (2021-04-23 05:59 UTC)
Platform Info:
OS: macOS (x86_64-apple-darwin18.7.0)
CPU: Intel(R) Core(TM) i7-8557U CPU @ 1.70GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-11.0.1 (ORCJIT, skylake)
```
---
### Author:  [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)
#### Post date: [June 10, 2021, 10:23am UTC](https://discourse.julialang.org/t/computing-sparse-orthogonal-projections/62652/5 "2021-06-10T10:23:52Z")
Side-note, If for some reason, your original data is of low rank, have a look at [LowRankApprox.jl](https://github.com/JuliaMatrices/LowRankApprox.jl#low-rank-factorizations) !
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [June 10, 2021, 12:13pm UTC](https://discourse.julialang.org/t/computing-sparse-orthogonal-projections/62652/6 "2021-06-10T12:13:36Z")
> [@Thomas](#):
>
> ```julia
> A = sprand(10,10,0.5)
> b = sprand(10,0.5)
> A\b
> ERROR: MethodError: no method matching ldiv!(::SuiteSparse.UMFPACK.UmfpackLU{Float64, Int64}, ::SparseVector{Float64, Int64})
>
> ```
Julia (and the underlying SuiteSparse library) only has support for solving sparse systems with dense vectors, not sparse vectors, so you need to do `A \ Vector(b)` here.
(There is not much point in supporting sparse vectors here, because the solution will almost certainly be dense in any case. Also, the storage for `A` swamps the storage for the vectors anyway.)
---
### Author:  [@Thomas](https://discourse.julialang.org/u/Thomas)
#### Post date: [June 10, 2021, 4:01pm UTC](https://discourse.julialang.org/t/computing-sparse-orthogonal-projections/62652/7 "2021-06-10T16:01:44Z")
Thanks for the clarification. I wanted to have an implementation for sparse(c) but you convinced me not to.