Looking for incomplete SVD -- need only some of the singular values

[Erik Schnetter]

I am looking for a routine that calculate the SVD (singular value decomposition) of a square, complex, dense, non-symmetric matrix. I am aware of Julia's SVD routine (and the respective LAPACK routines that one could call directly). However, I don't need all of the singular values -- I need only the largest one (or some of the largest ones), as well as the associated entries of U. Is there such a routine? I didn't find one in Julia.

I've looked for other packages that could be wrapped, but couldn't find any that offers this feature. The only thing I found is a description of the "svds" Matlab routine, which is apparently based on "eigs".

-erik

--
Erik Schnetter <schn...@gmail.com>
http://www.perimeterinstitute.ca/personal/eschnetter/

My email is as private as my paper mail. I therefore support encrypting
and signing email messages. Get my PGP key from https://sks-keyservers.net.

[Christian Peel]

What's the reason that you don't want to use Julia's SVD directly?   Is your matrix especially large?   Some links that may be useful:

Approximating the SVD using Julia:
   http://beowulf.lcs.mit.edu/18.337/projects/Turner-Presentation_SVD-Julia.pdf
   https://github.com/alexjturner/SVDapprox
Approximating the SVD using the power method
   http://www.sci.ccny.cuny.edu/~szlam/npisvdnipsshort.pdf
   https://en.wikipedia.org/wiki/Power_iteration

--

chris...@ieee.org

[Erik Schnetter]

I am looking for an algorithm that produces only some of the singular values, not all of them. I assume that would be cheaper. I didn't see such a function in Julia's documentation.

-erik

[Viral Shah]

There's svds, which uses Arnoldi iterations.

-viral

[Jiahao Chen]

> I've looked for other packages that could be wrapped, but couldn't find any that offers this feature. The only thing I found is a description of the "svds" Matlab routine, which is apparently based on "eigs". 

If you want a canned routine, you're best off with svds(A) in base Julia (wraps ARPACK) or tsvd(A) in https://github.com/andreasnoack/PROPACK.jl (work in progress to wrap PROPACK). The other package I'm aware of is SLEPc, but it requires the entire PETSc stack.

Note that the singular values of A are essentially the square root of the magnitude of the eigenvalues of AA', or A'A, or [0 A; A' 0], so eigenvalue routines can be used to solve the singular value problem.

If you only need the largest singular value and its left singular vector, a simple quick and dirty DIY scheme is to do power iteration on AA':

A = randn(m, n)

u = randn(n) #some guess

for i=1:10 #some number of iterations

   u = A*(A'u)

   scale!(u, 1/norm(u))

end

σ = norm(A'u) #σ, u are the largest singular value and its left singular vector

There is a simple generalization to the largest k singular values and their left singular vectors.

[Andreas Noack]

I've tried to make the package that Jiahao mentioned usable. I think it works, but it probably still has some rough edges. You can find it here

https://github.com/andreasnoack/TSVD.jl

and there is a help entry for tsvd that explains the arguments.

For a 2000x2000 dense non-symmetric complex matrix my tsvd runs in

julia> @time tsvd(A);

elapsed time: 0.354467696 seconds (16 MB allocated, 0.64% gc time in 1 pauses with 0 full sweep)

The svds in base runs a bit slower

julia> @time svds(A);
elapsed time: 2.034591908 seconds (134 MB allocated, 0.25% gc time in 6 pauses with 0 full sweep)

[paul analyst]

eigs(A[, B], ; nev=6, which=”LM”, tol=0.0, maxiter=1000, sigma=nothing, ritzvec=true, v0=zeros((0, )))
-> (d[, v ], nconv, niter, nmult, resid)
eigs computes eigenvalues d of A using Lanczos or Arnoldi iterations for real symmetric or general nonsymmetric matrices

nev= number of eigens...
Paul