1

I only know of the following power iteration. But it needs to create a huge matrix A'*A when both of rows and columns are pretty large. And A is a dense matrix as well. Is there any alternative to power iteration method below? I have heard of krylov subspace method, but I am not familiar with it. In anycase I am looking for any faster method than the one mentioned below:

B = A'*A; % or B = A*A' if it is smaller
x = B(:,1); % example of starting point, x will have the largest eigenvector 
x = x/norm(x); 
for i = 1:200 
  y = B*x; 
  y = y/norm(y);
  % norm(x - y); % <- residual, you can try to use it to stop iteration
  x = y; 
end; 
n3 = sqrt(mean(B*x./x)) % translate eigenvalue of B to singular value of A
amoeba
  • 93,463
  • 28
  • 275
  • 317
user3086871
  • 583
  • 4
  • 11
  • yes I am looking for existing solutions. My matrix is roughly 8 million by 150000. so taking the smaller dimension 150k, it will be 150k by 150k taking roughly 146GB of space. – user3086871 May 23 '17 at 18:30
  • It is dense I am afraid. – user3086871 May 23 '17 at 18:35
  • 1
    The answer is to use a package, because mostly anything you try by hand will be woefully inefficient in python. For example check out: https://scicomp.stackexchange.com/questions/7566/memory-efficient-implementations-of-partial-singular-value-decompositions-svd – Alex R. May 23 '17 at 19:05
  • See https://stats.stackexchange.com/questions/2806 and https://stats.stackexchange.com/questions/41259. – amoeba May 23 '17 at 19:15
  • Possible duplicate of [How to compute SVD of a huge sparse matrix?](https://stats.stackexchange.com/questions/41259/how-to-compute-svd-of-a-huge-sparse-matrix) – amoeba May 23 '17 at 19:15
  • 1
    @amoeba are you sure that is a duplicate as the OP said it was a dense matrix? – mdewey May 24 '17 at 08:20

1 Answers1

3

There are two approaches that require only small number of multiplications of vectors by $A$ or $A'$, and so have computational complexity $O(MNk)$ for small fixed $k$. Their performance is empirically similar, too.

You absolutely do not want to implement a Lanczos-type algorithm yourself -- the rounding error is subtle and quick to anger -- but implementations are available in many packages.

You actually can implement Stochastic SVD yourself; there's a very short Matlab implementation at that link.

Thomas Lumley
  • 21,784
  • 1
  • 22
  • 73