randomized SVD decomposes a matrix by extracting the first k singular values/vectors using k+p random projections.

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.

    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

They should be approximation to the true singular values (with suitable hypotheses, with good probability...). Intuitively, it is like trying to infer the singular values of a $n\times n$ matrix by looking at its leading $(k+p)\times (k+p)$ submatrix only and replacing the rest with zeros: it is cheaper to compute, but yields only an approximate result.

Where does the "randomized" part come into play? Well, depending on the specific matrix, the $(k+p)\times (k+p)$ submatrix may not be representative of the whole matrix: so, we conjugate everything with a random orthogonal matrix $Q$ that mixes up everything and ensures that there is "nothing special" regarding the first $k+p$ entries with respect to the rest.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.