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 for this kind of projection, in order to avoid explicitly computing A^* A and squaring the condition number (see also the discussion here). 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:

```
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

```
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:

```
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:

```
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`

).