# 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:** 1
**Showing post:** 2

<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: [June 10, 2021, 1:22am UTC](https://discourse.julialang.org/t/computing-sparse-orthogonal-projections/62652/2 "2021-06-10T01:22:39Z")

</div>

> [@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`).

---

_[View the full topic](https://discourse.julialang.org/t/computing-sparse-orthogonal-projections/62652)._