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).