Randomized SVD and singular values

Question

Randomized SVD decomposes a matrix by extracting the first k singular values/vectors using k+p random projections. This works surprisingly well for large matrices.

My question concerns the singular values that are output from the algorithm. Why aren't the values equal to the first k-singular values if you do the full SVD?

Below I have a simple implementation in R. Any suggestions on improving the performance would be appreciated.

 rsvd = function(A, k=10, p=5){
       n = nrow(A)
       y = A %*% matrix(rnorm(n * (k+p)), nrow=n)
       q = qr.Q(qr(y))
       b = t(q) %*% A
       svd = svd(b)
       list(u=q %*% svd$u, d=svd$d, v=svd$v)
    }

    > set.seed(10)

    > A <- matrix(rnorm(500*500),500,500)

    > svd(A)$d[1:15]
     [1] 44.94307 44.48235 43.78984 43.44626 43.27146 43.15066 42.79720 42.54440 42.27439 42.21873 41.79763 41.51349 41.48338 41.35024 41.18068

    > rsvd.o(A,10,5)$d
     [1] 34.83741 33.83411 33.09522 32.65761 32.34326 31.80868 31.38253 30.96395 30.79063 30.34387 30.04538 29.56061 29.24128 29.12612 27.61804

    B <- matrix(rnorm(500*50),500,500)  # rank 50

> rsvd(B,10,5)$d
 [1] 86.48035 83.02114 81.03988 80.04358 77.24979 76.10945 74.47357 74.08382
 [9] 72.85898 72.06897 69.59526 67.70750 66.53867 62.96446 61.50838

> svd(B)$d[1:15]
 [1] 92.44779 91.47689 88.71948 88.08170 87.24533 85.13312 84.14741 83.71757
 [9] 82.80832 81.43005 80.73903 79.92959 78.87421 78.33509 77.38431

As Joris pointed out, I have this posted on stackoverflow as well. you can find the revelant conversation here

https://stackoverflow.com/questions/4224031/randomized-svd-singular-values

Also see the relevant paper by Martinsson et al: A randomized algorithm for the decomposition of matrices

Answer 1

We have implemented this (along with a power iteration refinement) in the scikit-learn python package.

Our implementation is able to find the exact same singular values and vectors if k + p > rank(M) as demonstrated in the tests.

If you cut (k + p) before reaching near zero singular values (i.e. in the k + p < rank(M) case) then the singular vectors are indeed different from the ones you get with the un-truncated version but they might still be very useful in practice for features extraction in machine learning: for instance 'truncated' eigenfaces at 150 work as good for face recognition task with SVM as the top 150 first singular vectors of the full decomposition even though the rank of my faces dataset seems to be much higher.

This randomized / truncated SVD method looks really interesting in practice: it can really cut down the computation time as shown in this benchmark:

Answer 2

I do not think the singular values should match those of the full matrix. You are computing an approximation of the input matrix by projection onto $k+p$ random vectors. For a rank $k+p$ matrix to approximate a rank $n \gg k+p$ matrix, the trace should probably be the same, but then if the first $k$ singular values are to overlap, you have to push a lot of 'variance' to the last $p$ singular values of the approximation (probably so many that they are no longer the least significant singular values).

Another way of looking at this is one is approximating $A = U\Sigma V'$ by another decomposition, $T \Gamma W'$. We should not expect a fast randomized algorithm to magically work such that $T$ is the first $k$ columns of $U$, $\Gamma$ is a submatrix of $\Sigma$, etc.