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.