Merge pull request #9425 from jaak-s/svds

initial svds support based on eigs
JuliaLang · Jan 11, 2015 · 2acaf34 · 2acaf34
2 parents 02c0501 + b9d7a30
commit 2acaf34
Show file tree

Hide file tree

Showing 5 changed files with 92 additions and 1 deletion.
diff --git a/base/exports.jl b/base/exports.jl
@@ -695,6 +695,7 @@ export
     svd,
     svdfact!,
     svdfact,
+    svds,
     svdvals!,
     svdvals,
     sylvester,

diff --git a/base/linalg.jl b/base/linalg.jl
@@ -102,6 +102,7 @@ export
     svd,
     svdfact!,
     svdfact,
+    svds,
     svdvals!,
     svdvals,
     sylvester,

diff --git a/base/linalg/arnoldi.jl b/base/linalg/arnoldi.jl
@@ -106,3 +106,36 @@ function eigs(A, B;
                                  resid, ncv, v, ldv, sigma, iparam, ipntr, workd, workl, lworkl, rwork)
 
 end
+
+
+## svds
+
+type SVDOperator{T,S} <: AbstractArray{T, 2}
+    X::S
+    m::Int
+    n::Int
+    SVDOperator(X::S) = new(X, size(X,1), size(X,2))
+end
+
+## v = [ left_singular_vector; right_singular_vector ]
+*{T,S}(s::SVDOperator{T,S}, v::Vector{T}) = [s.X * v[s.m+1:end]; s.X' * v[1:s.m]]
+size(s::SVDOperator)  = s.m + s.n, s.m + s.n
+issym(s::SVDOperator) = true
+
+function svds{S}(X::S; nsv::Int = 6, ritzvec::Bool = true, tol::Float64 = 0.0, maxiter::Int = 1000)
+  if nsv > minimum(size(X)); error("nsv($nsv) should be at most $(minimum(size(X)))"); end
+
+  otype = eltype(X)
+  ex    = eigs(SVDOperator{otype,S}(X), I; ritzvec = ritzvec, nev = 2*nsv, tol = tol, maxiter = maxiter)
+  ind   = [1:2:nsv*2]
+  sval  = abs(ex[1][ind])
+
+  if ! ritzvec
+    return sval, ex[2], ex[3], ex[4], ex[5]
+  end
+
+  ## calculating singular vectors
+  left_sv  = sqrt(2) * ex[2][ 1:size(X,1),     ind ] .* sign(ex[1][ind]')
+  right_sv = sqrt(2) * ex[2][ size(X,1)+1:end, ind ]
+  return left_sv, sval, right_sv, ex[3], ex[4], ex[5], ex[6]
+end
diff --git a/doc/stdlib/linalg.rst b/doc/stdlib/linalg.rst
@@ -568,7 +568,7 @@ Linear algebra functions in Julia are largely implemented by calling functions f
 
    Conjugate transpose array ``src`` and store the result in the preallocated array ``dest``, which should have a size corresponding to ``(size(src,2),size(src,1))``. No in-place transposition is supported and unexpected results will happen if `src` and `dest` have overlapping memory regions.
 
-.. function:: 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)
+.. function:: eigs(A, [B,]; nev=6, which="LM", tol=0.0, maxiter=300, 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 respectively. If ``B`` is provided, the generalized eigen-problem is solved.  The following keyword arguments are supported:
     * ``nev``: Number of eigenvalues
@@ -604,6 +604,23 @@ Linear algebra functions in Julia are largely implemented by calling functions f
       real or complex inverse with level shift ``sigma`` :math:`(A - \sigma I )^{-1}`
       =============== ================================== ==================================
 
+.. function:: svds(A; ritzvec=true, args...) -> (left_sv, s, right_sv, nconv, niter, nmult, resid)
+
+   ``svds`` computes singular values ``s`` of ``A`` using Lanczos or Arnoldi iterations. Uses ``eigs`` underneath so following keyword arguments are supported:
+    * ``nev``: Number of singular values.
+    * ``ncv``: Number of Krylov vectors used in the computation; see ``eigs`` manual.
+    * ``ritzvec``: Whether to return the left and right singular vectors ``left_sv`` and ``right_sv``, default is ``true``. If ``false`` the singular vectors are omitted from the output.
+    * ``which``: type of singular values (and vectors) to compute, default is largest values. See ``eigs`` manual.
+    * ``tol``: tolerance, see ``eigs``.
+    * ``maxiter``: Maximum number of iterations, see ``eigs``.
+    * ``sigma``: See ``eigs``.
+    * ``v0``: starting vector of right singular vector from which to start the iterations.
+
+   **Example**::
+
+      X = sprand(10, 5, 0.2)
+      svds(X, nev = 2)
+
 .. function:: peakflops(n; parallel=false)
 
    ``peakflops`` computes the peak flop rate of the computer by using double precision :func:`Base.LinAlg.BLAS.gemm!`. By default, if no arguments are specified, it multiplies a matrix of size ``n x n``, where ``n = 2000``. If the underlying BLAS is using multiple threads, higher flop rates are realized. The number of BLAS threads can be set with ``blas_set_num_threads(n)``.

diff --git a/test/linalg/arnoldi.jl b/test/linalg/arnoldi.jl
@@ -0,0 +1,39 @@
+using Base.Test
+
+let # svds test
+    A = sparse([1, 1, 2, 3, 4], [2, 1, 1, 3, 1], [2.0, -1.0, 6.1, 7.0, 1.5])
+    S1 = svds(A, nsv = 2)
+    S2 = svd(full(A))
+
+    ## singular values match:
+    @test_approx_eq S1[2] S2[2][1:2]
+
+    ## 1st left singular vector
+    s1_left = sign(S1[1][3,1]) * S1[1][:,1]
+    s2_left = sign(S2[1][3,1]) * S2[1][:,1]
+    @test_approx_eq s1_left s2_left
+
+    ## 1st right singular vector
+    s1_right = sign(S1[3][3,1]) * S1[3][:,1]
+    s2_right = sign(S2[3][3,1]) * S2[3][:,1]
+    @test_approx_eq s1_right s2_right
+end
+
+let # complex svds test
+    A = sparse([1, 1, 2, 3, 4], [2, 1, 1, 3, 1], exp(im*[2.0:2:10]))
+    S1 = svds(A, nsv = 2)
+    S2 = svd(full(A))
+
+    ## singular values match:
+    @test_approx_eq S1[2] S2[2][1:2]
+
+    ## left singular vectors
+    s1_left = abs(S1[1][:,1:2])
+    s2_left = abs(S2[1][:,1:2])
+    @test_approx_eq s1_left s2_left
+
+    ## right singular vectors
+    s1_right = abs(S1[3][:,1:2])
+    s2_right = abs(S2[3][:,1:2])
+    @test_approx_eq s1_right s2_right
+end