# How to project a vector onto a very large, non-orthogonal subspace

I have a difficult problem.

I have a very large, non-orthogonal matrix $A$ and need to project the vector $y$ onto the subspace spanning the columns of $A$. If this were a small matrix, I would use Gram-Schmidt or just compute $A(A^TA)^{-1}A^Ty$. Unfortunately, $A$ is just too big to do this in a feasible amount of time. Is there a good trick to do this projection?

EDIT2: Some background info: I work in signal processing and work with Fusion Frames quite a bit. For those unfamiliar, it is the problem of representing a signal with a union of subspaces from an overcomplete dictionary, typically using a block-sparse solver (in my case, a program called SPGL1). The problem is, I need to project my signal vector onto this dictionary before performing the inversion. The dictionary is coherent so under normal circumstances I need to use the formula above or a method like Gram-Schmidt to get an orthonormal basis.

Right now, my working solution is to approximate this projection by using the truncated SVD and using the $U$ matrix as an approximate orthonormal basis (a la PCA). But, it would be nice if anyone knew any iterative projection algorithms. I know this isn't a scientific computing site, but I figured there may be someone here who knows some linear algebra tricks.

• The problem is that you should be able to decide whether the columns span the entire space at the very least (i.e., whether the projection is the vector itself). I see no way to do it in fewer than $n^3$ operations but maybe I'm just too stupid. – fedja Oct 24 '13 at 1:31
• if you are willing to put up with approximation, why not just solve a least squares problem using an iterative method? – Suvrit Oct 24 '13 at 2:47
• @suv....rit: do you know any such algorithms? I have looked in journal articles for projection algorithms but haven't gotten very far. – rlbond Oct 24 '13 at 14:07
• OP has indicated that he/she is a graduate student, that this is a serious research problem arising in thesis work, and avers that the question was closed without due consideration. Would the closers mind explaining their decision to close? I've opened a meta discussion: meta.mathoverflow.net/questions/1085/… – Todd Trimble Oct 24 '13 at 17:00
• rlbond, now that the question has been reopened, and seeing that several people in the meta discussion wished that more background had been given, would you be willing to provide some? Also, it might help draw the right people to your question if other tags were added (numerical analysis?) -- please consider. – Todd Trimble Oct 24 '13 at 19:38

Depending on how many rows, maybe you could use some statistics results.

When I was doing Principal Component Analysis I ended up having a matrix $X_{75\times 375000}$ and I could still find a good vector space that described well these vectors (and projections).

That's the whole point of Principal Component Analysis. Since the amount of information you have is just too much, you can get rid of some redundancy and reduce the dimension (really, my case was a 75 dimension space with that huge matrix).

Maybe these articles can give you the basic idea. http://www.cs.otago.ac.nz/cosc453/student_tutorials/principal_components.pdf http://www.stat.cmu.edu/~cshalizi/490/pca/pca-handout.pdf

It's a well known algorithm, mainly for computer vision (Eigenfaces). Maybe it won't be useful at all, since you want the span of these vectors and this is an approximation method. Anyway, I always consider this algorithm when some gigantic vectors need to be analyzed.

• I ended up trying something similar to this and it works pretty well actually. In MATLAB I did [U,S] = svds(A, 150); y_proj = U * (U' * y); and it worked pretty well. I always forget about simple PCA... – rlbond Oct 25 '13 at 1:11

It makes a considerable difference if you need to project a vector just once or repeatedly inside some loop of another algorithm. If you would be in the latter case than it would be indeed a good idea to follow suv...rits suggestion and use an iterative algorithm to solve the normal equations (conjugate gradients, probably with preconditioning, suggest themselves). In this case you can often just use a small and fixed number of CG iterations and restart the next CG iteration with what you've got from the previous run. If the point you want to project does not change very much in the outer loop, CG should be converging eventually.

However, your EDIT2 suggests that your only project you data once. Without knowing your precise approximation problem, it could be possible to make a reformulation of the problem: If $y$ is your signal but you need to put $z = A(A^T A)^{-1}A^Ty$ in your algorithm, you could just put in $z$ and add the constraints $z = Aw$ and $A^T Aw = A^T y$, i.e. you introduce two new variables and two new linear constraints. Then use some algorithm for your new inversion which can handle linear constraints. Without knowing more details I can't be more concrete here…