How to do eigendecomposition on giant dense PSD matrices? I need to perform eigendecomposition on giant matrices (at least with dimensions of 600K by 600K). I need both eigenvalues and eigenvectors, however, only top k of them, e.g. k=100. In addition, the matrices are known to be positive definite (if that helps in any way). However, unfortunately the matrices are dense and thus sparsity-based approaches do not apply.
I think the only hope is to perform some form of random sampling of the giant matrix (e.g. rows or columns, or small sub-matrices), then do some computation on it, and repeat this in a loop in a way that each iteration can progressively produce a better estimate of the top k eigenvalues and eigenvectors.
Is this possible? If so, could you please advise me about the method?
Thank you!
Golabi
 A: If you have the ability to compute matrix-vector products, then a Krylov subspace methods such as the Lanczos method are available to you. These methods are essentially the "power method on steroids". The highly respected ARPACK software has the ability to compute the dominant eigenvalues and eigenvectors of a symmetric matrix $A$ when only provided with a black box function to compute $x \mapsto Ax$.
Another family of methods that have received a lot of attention in recent years are randomized algorithms. The key insight is that if $A$ is multiplied by a random vector $\omega$, the largest components of $A\omega$ will be in the direction of the dominant eigenvectors. Based on this idea, if $\Omega$ is an $n\times \ell$ random matrix with IID Gaussian elements, the columns of $Y = A\Omega$ are very likely to contain within their span the first few eigenvectors of $A$. From here, a good low-rank approximation of $A$ is the Nystrom approximation $A \approx A\Omega (\Omega^* A \Omega)^{-1} (A\Omega)^*$. (Note that this approximation is exactly if $\ell = n$ so that $\Omega$ is square and invertible with probability $1$.) A numerically sound way of implementing this randomized idea is presented in Algorithm 16 of this recent review, which also contains analysis for this method. This kind of method is best for the eigenvalues of $A$ rapidly decline rather than slowly declining. This idea can be combined with power iteration (see Section 14.5 of the same review), but implementing this method is very delicate numerically. One attractive feature about these algorithms is that since matrix-multiplication is a linear operation, the product $Y = A\Omega$ can be computed as $Y = A_1\Omega + \cdots + A_n\Omega$, where $A_j$ is the matrix $A$ with all but the $j$th row zeroed out. Thus, this randomized algorithm can be used as a "streaming method" where the matrix $A$ is never stored all at once but is generated one row at a time. Once a row is processed, it is discarded.
