# Obtaining a projection onto the range of a matrix through its QR decomposition

**URL:** <https://discourse.julialang.org/t/obtaining-a-projection-onto-the-range-of-a-matrix-through-its-qr-decomposition/117334>\
**Category:** Numerics\
**Created:** [July 22, 2024, 10:56am UTC](https://discourse.julialang.org/t/obtaining-a-projection-onto-the-range-of-a-matrix-through-its-qr-decomposition/117334 "2024-07-22T10:56:21Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![cgois](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgois/32/37096_2.png) [@cgois](https://discourse.julialang.org/u/cgois)\
**Post date:** [July 22, 2024, 10:56am UTC](https://discourse.julialang.org/t/obtaining-a-projection-onto-the-range-of-a-matrix-through-its-qr-decomposition/117334/1 "2024-07-22T10:56:21Z")

</div>

Given a matrix A, I want to compute a matrix with rank(A) column vectors that span the range of A. With an SVD, this can be done with:

```julia
import LinearAlgebra as LA
function ortr_svd(A::AbstractMatrix{T}) where {T<:Number}
    dec = LA.svd(A)
    tol = maximum(dec.S) * maximum(size(A)) * eps(T)
    rank = sum(dec.S .> tol)
    dec.U[:, 1:rank] # The first `rank` columns of U span A.
end

```

Because this does not work with sparse matrices, I want to find an alternative using a QR decomposition. My attempt is:

```julia
function ortr_qr(A::AbstractMatrix{T}) where {T<:Number}
    dec = LA.qr(A, LA.ColumnNorm())
    tol = maximum(abs.(dec.R)) * maximum(size(A)) * eps(T)
    rank = sum(abs.(LA.Diagonal(dec.R)) .> tol)
    dec.Q[:, 1:rank]
end

```

This seems to work when A is dense, but otherwise I get the error

```julia
ERROR: Pivoting Strategies are not supported by `SparseMatrixCSC`s

```

Is there any workaround? For instance, can I do it without any pivoting, or can I make pivoting work on sparse arrays?

Thanks for any help!

---

<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:** [July 22, 2024, 12:17pm UTC](https://discourse.julialang.org/t/obtaining-a-projection-onto-the-range-of-a-matrix-through-its-qr-decomposition/117334/2 "2024-07-22T12:17:47Z")

</div>

> [@cgois](#):
>
> Because this does not work with sparse matrices, I want to find an alternative using a QR decomposition.

See also: [How to find the linearly independent columns (rows) of a matrix - #14 by stevengj](https://discourse.julialang.org/t/how-to-find-the-linearly-independent-columns-rows-of-a-matrix/109772/14) (although this wasn’t written with sparse QR in mind). This returns a list of which of the columns of A are independent (up to some tolerance), and you can then use those columns as a basis — this might be a good approach to adapt for sparse A since its columns are sparse (albeit not orthonormal).

What sizes and ranks of sparse matrices are you interested in?

---

<div class="post-metadata">

**Author:** ![abulak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abulak/32/28314_2.png) [@abulak](https://discourse.julialang.org/u/abulak)\
**Post date:** [July 22, 2024, 12:25pm UTC](https://discourse.julialang.org/t/obtaining-a-projection-onto-the-range-of-a-matrix-through-its-qr-decomposition/117334/3 "2024-07-22T12:25:13Z")

</div>

this is what I did with svd and qr:

> <https://github.com/kalmarek/SymbolicWedderburn.jl/blob/master/src/image_basis.jl#L82-L103>

Incidentally I do know the mathematical rank beforehand and **for my examples** this numerical `rank(F)` has always agreed.

Matrices I deal with are of sizes up to 60k (of course sparse ;).

---

<div class="post-metadata">

**Author:** ![cgois](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgois/32/37096_2.png) [@cgois](https://discourse.julialang.org/u/cgois)\
**Post date:** [July 23, 2024, 10:09am UTC](https://discourse.julialang.org/t/obtaining-a-projection-onto-the-range-of-a-matrix-through-its-qr-decomposition/117334/4 "2024-07-23T10:09:34Z")

</div>

@stevengj In fact I had read your post before, but got stuck in adapting it to sparse matrices because the `qr` does not allow for pivoting if the matrix is sparse. Together with the answer by @abulak, this can be solved by getting the permutations from the `rpivinv` field of a `SparseArrays.SPQR.QRSparse`.

I expect usage with matrices of quite varying sizes (order thousands to tens of thousands) and the rank should often be around 1/5 of the dimension.

@abulak I had previously attempted something quite similar, but I was using `F.prow` instead of `F.rpivinv` for the rows permutation vector, and this somehow was not working. Indeed, these two fields are not equal, but I cannot seem to find the difference between them. Would you happen to know this?

Thanks a lot to both!

---

<div class="post-metadata">

**Author:** ![abulak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abulak/32/28314_2.png) [@abulak](https://discourse.julialang.org/u/abulak)\
**Post date:** [July 29, 2024, 8:26am UTC](https://discourse.julialang.org/t/obtaining-a-projection-onto-the-range-of-a-matrix-through-its-qr-decomposition/117334/5 "2024-07-29T08:26:42Z")

</div>

These are inverses of each other – the difference is on which side you’re permuting rows. `F.prow` is permuting rows of `A`, its inverse will be permuting rows of `Q`. For solving linear equations with `\` the `prow` is more useful. For our purpose we need the correct rows of `Q`.

Btw. every time when you ask for `F.prow` you’re inverting `F.rpivinv` (which you can also access like this, no need to `getfield(F, :rpivinv)`.