How to calculate Singular Value Decomposition (SVD) with Krylov subspace methods? Background : I have used Krylov subspace methods for a rather long time now, designing algorithms with the help of posing and solving mostly sparse matrix*vector 2-norm minimization problems.
One thing I have not done so far is to perform Singular Value Decomposition (SVD). It is defined for any matrix $\bf A$, to find $\bf U,V,\Sigma$  so that :
$${\bf A = U\Sigma V}^*$$ where $\bf U,V$ are orthonormal and $\bf \Sigma$ is diagonal with exclusively non-negative real entries and ${\bf M}^*$ means conjugate transpose of $\bf M$.
How could this be achieved? Not only answers describing how but also references to works, articles, even blogs will be accepted.
 A: Let me sketch the outline of a method, using the same technique I've addressed in another answer here. The key observation is that an SVD of $A$ can readily obtained from an eigenvalue computation of the related matrix
$$
B = \begin{bmatrix} 0 & A \\ A^* & 0 \end{bmatrix}.
$$
In fact the nonzero eigenvalues of $B$ are precisely $\pm$ the nonzero singular values of $A$ with eigenvectors
$$
B \begin{bmatrix} u \\ \pm v\end{bmatrix} = \pm \sigma \begin{bmatrix} u \\ \pm v\end{bmatrix}
$$
where $u$ and $v$ are right and left singular vectors of $A$! One can thus adapt any Krylov subspace method for the symmetric eigenvalue problem (e.g. Lanczos) to compute an SVD by applying this method to $B$ and recovering from it a (full or partial) SVD of $A$.
In fact, you don't even have to write the matrix $B$ down. All a Kyrlov subspace method for a symmetric eigenvalue problem requires is the ability to compute products of the form $x \mapsto Bx$. But this can be done by using subroutines to multiply by $A$ and $A^*$ by using the fact that
$$
B\begin{bmatrix} y \\ z \end{bmatrix} = \begin{bmatrix} Az \\ A^*y \end{bmatrix}.
$$
Note that multiplying by $A^*$ is absolutely essential. Any reasonable Krylov subspace method for the SVD must utilize both multiplications with $A$ and $A^*$. To see why, consider the simple case when $A = uv^*$ for vectors $u$ and $v$.Then it is very difficult to determine $v$ in $A=uv^*$ by only multiplications with $A$. This is because, once you've figured out $u$, each product $Ax = (v^*x)u$ with only tells you $v^*x$. To completely determine $v$ from inner products with other vectors $x$, you'll need many many inner products (as many as the size of $v$).
Alternately, you could use Krylov subspace methods on $A^*A = V\Sigma^2 V^*$ and $AA^* = U\Sigma^2 U^*$. This has the disadvantage of needing to compute two matvecs every iteration and needing to solve two separate problems, one for $U$ and one for $V$, but the problem does become positive semidefinite, which may be useful.
