svds implementation

[ViralBShah]

It seems that the way octave implements svds is as follows:

(V, s) = eigs ([spzeros(m,m), b; b', spzeros(n,n)], b_k, b_sigma, b_opts);

If you have the octave repo, you can see the implementation in octave-3.6.4/scripts/sparse/svds.m.
I see some pre and post processing too.

Should we go ahead this way for our svds implementation too?

[jiahao]

I'm pretty sure this trick uses the relationship between singular values as "square roots" of the eigenvalues, and as such would be numerically unstable in that the computation works with a matrix that has a condition number which is roughly squared compared to the original matrix. (In the same vein as why normal equation methods like CGNE aren't currently favored in iterative methods)

Given that IterativeSolvers.jl now offers the Golub-Kahan-Lanczos method for computing singular values via svdvals_gkl(), I think we should close this and work on providing the full SVD factorization (i.e. return the singular vectors also) from IterativeSolvers.

[ViralBShah]

Agree.

[ViralBShah]

On second thoughts, it would be nice to have a working svds implementation, which could be replaced later.

[ViralBShah]

Related: #6610