Computing sparse orthogonal projections

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

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.

1 Like

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

13 Likes

Note that this code does not seem to be correct. You should either use inv or use \, but not both.

SparseArrays already provides sparse \ methods, including \ for factorization objects.

1 Like

Sorry that was a typo; I fixed it in the question post.

Do you have this error?

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)

Side-note, If for some reason, your original data is of low rank, have a look at LowRankApprox.jl !

1 Like

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

3 Likes

Thanks for the clarification. I wanted to have an implementation for sparse(c) but you convinced me not to.