computing the truncated SVD, one singular value/vector at a time

Question

Is there a truncated SVD algorithm that computes the singular values one at a time?

My problem: I would like to compute the first $k$ singular values (and singular vectors) of a large dense matrix $M$, but I don't know what an appropriate value of $k$ would be. $M$ is large, so for efficiency reasons, I would rather not evaluate the full SVD only to truncate off the smallest SV's afterwards.

Ideally, there would be a way to compute the singular values $\sigma_1, \sigma_2,\ldots$ serially, from largest ($\sigma_1$) to smallest ($\sigma_n$). That way, I could simply halt the computation after computing the $k$th singular value if $\sigma_k/\sigma_1$ falls below some threshold.

Does such an algorithm exist (preferably with a Python implementation)? In my googling around, I've only found truncated SVD functions that take k as a parameter, thus forcing you to guess it a priori.

Answer 1

There are a couple of options available if you want an approximate rank-k factorizations.

1. Strongly rank-revealing QR factorizations
2. Interpolative decomposition (ID) and other randomized techniques.

Generally speaking, they provide a factorization of the form \begin{equation}\| A - MN^T\| \leq \text{factor}\times \sigma_{k+1}(A) := \epsilon \end{equation}

An approximate factorization of the above form can be converted into a standard decomposition like QR or SVD using standard techniques. A good review is available in the paper by Halko, Martinsson and Tropp "Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions"

In terms of software an interface to ID algorithms is available in scipy (scipy.linalg.interpolative) http://docs.scipy.org/doc/scipy-dev/reference/linalg.interpolative.html that allows you the user to specify $\epsilon$.

Answer 2

(Edited, because I misread the question at first; you already know that there are routines available to calculate the first $k$ singular values.)

If you exclude the approach of calculating the whole SVD, partial SVD algorithms reduce to using iterative methods to solve a related Hermitian eigenvalue problem. So, one strategy you could take would be to hand-code this sort of thing yourself, and keep solving for the largest remaining unsolved singular value until you want to stop, using something like a shift-and-invert strategy. There may be elegant ways of doing this sort of thing in sophisticated packages like SLEPc.

Another strategy would be the following:

• Calculate the largest singular value $s_{1}$.
• Set the absolute tolerance of the sparse SVD routine to $\tau \cdot s_{1} \cdot f$, where $\tau$ is your threshold, and $0 < f \leq 1$ is some safety factor to determine how many possibly extraneous singular values you want to calculate.
• Call the sparse SVD routine.

If the sparse SVD routine calculates a thin SVD (and I can't see why it wouldn't), then this strategy gives you all of the singular values you want (plus possibly some extra ones), because values below the absolute tolerance will be treated as zero. In that case, you can use scipy.sparse.linalg.svds, noting that $k$ is an optional parameter, and that you don't have to specify it a priori.