Merge pull request #32 from lostella/lower-rank-fix

Low rank fix

src/Arpack.jl

@@ -8,7 +8,7 @@ Arnoldi and Lanczos iteration for computing eigenvalues
 module Arpack
 
 # To make Pkg aware that this dependency
-# will be injected by BinaryProvider. 
+# will be injected by BinaryProvider.
 using Libdl
 
 const depsfile = joinpath(@__DIR__, "..", "deps", "deps.jl")
@@ -326,20 +326,26 @@ function _svds(X; nsv::Int = 6, ritzvec::Bool = true, tol::Float64 = 0.0, maxite
     # ind   = [1:2:ncv;]
     # sval  = abs.(ex[1][ind])
 
-    svals = sqrt.(real.(ex[1]))
+    realex1 = real.(ex[1])
+    threshold = max(eps(real(otype))*realex1[1], eps(real(otype)))
+    firstzero = findfirst(v -> v <= threshold, realex1)
+    r = firstzero === nothing ? nsv : firstzero-1 # rank of the decomposition
+    realex1[r+1:end] .= zero(real(otype))
+    svals = sqrt.(realex1)
+
     if ritzvec
         # calculating singular vectors
         # left_sv  = sqrt(2) * ex[2][ 1:size(X,1),     ind ] .* sign.(ex[1][ind]')
         if size(X, 1) >= size(X, 2)
             V = ex[2]
             # We cannot assume that X*V is a Matrix even though V is. This is not
             # the case for e.g. LinearMaps.jl so we convert to Matrix explicitly
-            U = Array(qr(rmul!(convert(Matrix, X*V), Diagonal(inv.(svals)))).Q)
+            U = Array(qr(rmul!(convert(Matrix, X*V), Diagonal([inv.(svals[1:r]); ones(nsv-r)]))).Q)
         else
             U = ex[2]
             # We cannot assume that X'U is a Matrix even though U is. This is not
             # the case for e.g. LinearMaps.jl so we convert to Matrix explicitly
-            V = Array(qr(rmul!(convert(Matrix, X'U), Diagonal(inv.(svals)))).Q)
+            V = Array(qr(rmul!(convert(Matrix, X'U), Diagonal([inv.(svals[1:r]); ones(nsv-r)]))).Q)
         end
 
         # right_sv = sqrt(2) * ex[2][ size(X,1)+1:end, ind ]

test/runtests.jl

@@ -293,3 +293,33 @@ LinearAlgebra.adjoint(A::MyOp) = MyOp(adjoint(A.mat))
     A = MyOp(randn(10, 9))
     @test svds(A, v0 = ones(9))[1].S == svds(A.mat, v0 = ones(9))[1].S
 end
+
+@testset "low rank" begin
+    Random.seed!(123)
+    @testset "$T coefficients" for T in [Float64, Complex{Float64}]
+        @testset "rank $r" for r in [2, 5, 10]
+            m, n = 3*r, 4*r
+            nsv = 2*r
+
+            FU = qr(randn(T, m, r))
+            U = Matrix(FU.Q)
+            S = 0.1 .+ sort(rand(r), rev=true)
+            FV = qr(randn(T, n, r))
+            V = Matrix(FV.Q)
+
+            A = U*Diagonal(S)*V'
+            F = svds(A, nsv=nsv)[1]
+
+            @test F.S[1:r] ≈ S
+            if T == Complex{Float64}
+                # This test fails since ARPACK does not have an Hermitian solver
+                # for the complex case. This problem occurs for U in the "fat"
+                # case. In the "tall" case the same may happen for V instead.
+                @test_broken F.U'*F.U ≈ Matrix{T}(I, nsv, nsv)
+            else
+                @test F.U'*F.U ≈ Matrix{T}(I, nsv, nsv)
+            end
+            @test F.V'*F.V ≈ Matrix{T}(I, nsv, nsv)
+        end
+    end
+end