check for zeros on the diagonal, implement wide form

Author: simonbyrne

src/GenericSVD.jl

@@ -5,11 +5,14 @@ import Base: SVD, svdvals!, svdfact!
 include("utils.jl")
 include("bidiagonalize.jl")
 
-
-
 function svdfact!(X; sorted=true, thin=true)
     m,n = size(X)
-    m >= n || error("Generic SVD requires more rows than columns.")
+    t =false
+    if m < n
+        m,n = n,m
+        X = X'
+        t = true
+    end
     B,P = bidiagonalize_tall!(X)
     U,Vt = full(P,thin=thin)
     U,S,Vt = svd!(B,U,Vt)
@@ -27,11 +30,15 @@ function svdfact!(X; sorted=true, thin=true)
         U = U[:,I]
         Vt = Vt[I,:]
     end
-    SVD(U,S,Vt)
+    t ? SVD(Vt',S,U') : SVD(U,S,Vt)
 end
 
 function svdvals!(X; sorted=true)
-    B,P = bidiagonalize_tall!(copy(X))
+    m,n = size(X)
+    if m < n
+        X = X'
+    end
+    B,P = bidiagonalize_tall!(X)
     S = svd!(B)[2]
     for i = eachindex(S)
         if signbit(S[i])
@@ -45,11 +52,16 @@ end
 """
 Tests if the B[i-1,i] element is approximately zero, using the criteria
 ```math
-    |B_{i-1,i}| < ɛ*(|B_{i-1,i-1}| + |B_{i,i}|)
+    |B_{i-1,i}| ≤ ɛ*(|B_{i-1,i-1}| + |B_{i,i}|)
 ```
 """
-offdiag_approx_zero(B::Bidiagonal,i,ɛ) =
-    abs(B.ev[i-1]) < ɛ*(abs(B.dv[i-1]) + abs(B.dv[i]))
+function offdiag_approx_zero(B::Bidiagonal,i,ɛ)
+    iszero = abs(B.ev[i-1]) ≤ ɛ*(abs(B.dv[i-1]) + abs(B.dv[i]))
+    if iszero
+        B.ev[i-1] = 0
+    end
+    iszero
+end
 
 
 """
@@ -70,34 +82,112 @@ function svd!{T<:Real}(B::Bidiagonal{T}, U=nothing, Vt=nothing, ɛ::T = eps(T))
     n = size(B, 1)
     n₂ = n
 
+    maxB = max(maxabs(B.dv),maxabs(B.ev))
+
     if istriu(B)
         while true
+            @label mainloop
+
             while offdiag_approx_zero(B,n₂,ɛ)
                 n₂ -= 1
                 if n₂ == 1
-                    return U,B.dv,Vt
+                    @goto done
                 end
             end
+
+
+
             n₁ = n₂ - 1
+            # check for diagonal zeros
+            if abs(B.dv[n₁]) ≤ ɛ*maxB
+                svd_zerodiag_row!(U,B,n₁,n₂)
+                @goto mainloop
+            end
             while n₁ > 1 && !offdiag_approx_zero(B,n₁,ɛ)
                 n₁ -= 1
+                # check for diagonal zeros
+                if abs(B.dv[n₁]) ≤ ɛ*maxB
+                    svd_zerodiag_row!(U,B,n₁,n₂)
+                    @goto mainloop
+                end
+            end
+
+            if abs(B.dv[n₂]) ≤ ɛ*maxB
+                svd_zerodiag_col!(B,Vt,n₁,n₂)
+                @goto mainloop
             end
 
-            # TODO: check for diagonal zeros
 
             d₁ = B.dv[n₂-1]
             d₂ = B.dv[n₂]
             e  = B.ev[n₂-1]
-            
+
             s₁, s₂ = svdvals2x2(d₁, d₂, e)
-            # use singular value closest to
+            # use singular value closest to sqrt of final element of B'*B
             h = hypot(d₂,e)
-            shift = abs(s₁-h) < abs(s₂-h) ? s₁ : s₂            
+            shift = abs(s₁-h) < abs(s₂-h) ? s₁ : s₂
             svd_gk!(B, U, Vt, n₁, n₂, shift)
         end
     else
         throw(ArgumentError("lower bidiagonal version not implemented yet"))
     end
+    @label done
+    U, B.dv, Vt
+end
+
+
+"""
+Sets B[n₁,n₁] to zero, then zeros out row n₁ by applying sequential row (left) Givens rotations up to n₂.
+"""
+function svd_zerodiag_row!(U,B,n₁,n₂)
+    e = B.ev[n₁]
+    B.dv[n₁] = 0 # set to zero
+    B.ev[n₁] = 0
+
+    for i = n₁+1:n₂
+        # n₁ [0 ,e ] = G * [e ,0 ]
+        #    [ ... ]       [ ... ]
+        # i  [dᵢ,eᵢ]       [dᵢ,eᵢ]
+        dᵢ = B.dv[i]
+
+        G,r = givens(dᵢ,e,i,n₁)
+        A_mul_Bc!(U,G)
+        B.dv[i] = r # -G.s*e + G.c*dᵢ
+
+        if i < n₂
+            eᵢ = B.ev[i]
+            e       = G.s*eᵢ
+            B.ev[i] = G.c*eᵢ
+        end
+    end
+end
+
+
+"""
+Sets B[n₂,n₂] to zero, then zeros out column n₂ by applying sequential column (right) Givens rotations up to n₁.
+"""
+function svd_zerodiag_col!(B,Vt,n₁,n₂)
+    e = B.ev[n₂-1]
+    B.dv[n₂] = 0 # set to zero
+    B.ev[n₂-1] = 0
+
+    for i = n₂-1:-1:n₁
+        #   i      n₂     i      n₂
+        #  [eᵢ,...,e ] = [eᵢ,...,0 ] * G'
+        #  [dᵢ,...,0 ]   [dᵢ,...,e ]
+        dᵢ = B.dv[i]
+
+        G,r = givens(dᵢ,e,i,n₂)
+        A_mul_B!(G,Vt)
+
+        B.dv[i] = r # G.c*dᵢ + G.s*e
+
+        if n₁ < i
+            eᵢ = B.ev[i-1]
+            e       = -G.s*eᵢ
+            B.ev[i-1] = G.c*eᵢ
+        end
+    end
 end
 
 
@@ -110,15 +200,15 @@ A Givens rotation is applied to the top 2x2 matrix, and the resulting "bulge" is
 function svd_gk!{T<:Real}(B::Bidiagonal{T},U,Vt,n₁,n₂,shift)
 
     if istriu(B)
-        
+
         d₁′ = B.dv[n₁]
         e₁′ = B.ev[n₁]
         d₂′ = B.dv[n₁+1]
 
         G, r = givens(d₁′ - abs2(shift)/d₁′, e₁′, n₁, n₁+1)
         A_mul_B!(G, Vt)
 
-        #  [d₁,e₁] = [d₁′,e₁′] * G
+        #  [d₁,e₁] = [d₁′,e₁′] * G'
         #  [b ,d₂]   [0  ,d₂′]
 
 
@@ -129,48 +219,48 @@ function svd_gk!{T<:Real}(B::Bidiagonal{T},U,Vt,n₁,n₂,shift)
 
         for i = n₁:n₂-2
 
-            #  [. ,e₁′,b′ ] = G * [d₁,e₁,0 ] 
+            #  [. ,e₁′,b′ ] = G * [d₁,e₁,0 ]
             #  [0 ,d₂′,e₂′]       [b ,d₂,e₂]
 
             e₂ = B.ev[i+1]
 
             G, r = givens(d₁, b, i, i+1)
             A_mul_Bc!(U, G)
 
-            B.dv[i] =  G.c*d₁ + G.s*b
-            
+            B.dv[i] =  r # G.c*d₁ + G.s*b
+
             e₁′ =  G.c*e₁ + G.s*d₂
             d₂′ = -G.s*e₁ + G.c*d₂
-            
+
             b′  =  G.s*e₂
             e₂′ =  G.c*e₂
 
-            #  [. ,0 ] = [e₁′,b′ ] * G
+            #  [. ,0 ] = [e₁′,b′ ] * G'
             #  [d₁,e₁]   [d₂′,e₂′]
             #  [b ,d₂]   [0  ,d₃′]
 
             d₃′ = B.dv[i+2]
 
             G, r = givens(e₁′, b′, i+1, i+2)
             A_mul_B!(G, Vt)
-            
-            B.ev[i] = e₁′*G.c + b′*G.s
-            
+
+            B.ev[i] = r # e₁′*G.c + b′*G.s
+
             d₁ =  d₂′*G.c + e₂′*G.s
             e₁ = -d₂′*G.s + e₂′*G.c
 
             b  = d₃′*G.s
             d₂ = d₃′*G.c
         end
 
-        #  [. ,.] = G * [d₁,e₁] 
+        #  [. ,.] = G * [d₁,e₁]
         #  [0 ,.]       [b ,d₂]
 
         G, r = givens(d₁,b,n₂-1,n₂)
         A_mul_Bc!(U, G)
-        
-        B.dv[n₂-1] =  G.c*d₁ + G.s*b
-        
+
+        B.dv[n₂-1] =  r # G.c*d₁ + G.s*b
+
         B.ev[n₂-1] =  G.c*e₁ + G.s*d₂
         B.dv[n₂]   = -G.s*e₁ + G.c*d₂
     else
@@ -198,7 +288,7 @@ function svdvals2x2(f, h, g)
 
     fhmin = min(fa,ha)
     fhmax = max(fa,ha)
-    
+
     if fhmin == 0
         ssmin = zero(f)
         if fhmax == 0

src/utils.jl

@@ -30,3 +30,17 @@ function A_ldiv_B!{Ta,Tb}(A::SVD{Ta}, B::StridedVecOrMat{Tb})
     k = searchsortedlast(A.S, eps(real(Ta))*A.S[1], rev=true)
     sub(A.Vt,1:k,:)' * (sub(A.S,1:k) .\ (sub(A.U,:,1:k)' * B))
 end
+
+
+# we have to define our own givens function due to ordering restriction in Base (#14936)
+function givens{T}(f::T, g::T, i1::Integer, i2::Integer)
+    if i1 == i2
+        throw(ArgumentError("Indices must be distinct."))
+    end
+    c, s, r = Base.LinAlg.givensAlgorithm(f, g)
+    if i1 > i2
+        s = -s
+        i1,i2 = i2,i1
+    end
+    Base.LinAlg.Givens(i1, i2, convert(T, c), convert(T, s)), r
+end

test/runtests.jl

@@ -12,21 +12,79 @@ x,y = GenericSVD.svdvals2x2(a,b,c)
 @test sort(sqrt(eigvals(U'*U))) ≈ [x,y]
 @test sort(svdvals(U)) ≈ [x,y]
 
+U = eye(3)
+B = Bidiagonal([0.0,1.0,2.0],[3.0,4.0],true)
+B1 = copy(B)
+
+GenericSVD.svd_zerodiag_row!(U,B,1,3)
+@test B[1,1] == 0
+@test B[1,2] == 0
+@test U*full(B) ≈ B1
+
+Vt = eye(3)
+B = Bidiagonal([1.0,2.0,0.0],[3.0,4.0],true)
+B1 = copy(B)
+
+GenericSVD.svd_zerodiag_col!(B,Vt,1,3)
+@test B[3,3] == 0
+@test B[2,3] == 0
+@test full(B)*Vt ≈ B1
+
+
+
 n,m = 100,20
 
 X = randn(n,m)
 bX = big(X)
 bS = svdfact(bX)
 @test isapprox(full(bS), bX, rtol=1e3*eps(BigFloat))
 @test isapprox(svdvals(bX), svdvals(X), rtol=1e3*eps())
-
+@test bX == X # check we didn't modify the input
 
 bY = big(randn(n))
 @test isapprox(qrfact(bX,Val{false}) \ bY, bS \ bY, rtol=1e3*eps(BigFloat))
+@test bX == X # check we didn't modify the input
+
+bXt = bX'
+bSt = svdfact(bXt)
+@test isapprox(full(bSt), bXt, rtol=1e3*eps(BigFloat))
+@test isapprox(svdvals(bXt), svdvals(X), rtol=1e3*eps())
+@test bXt == X' # check we didn't modify the input
+
+X = Float64[1 2 0; 0 1 2; 0 0 0]
+bX = big(X)
+bS = svdfact(bX)
+@test isapprox(full(bS), bX, rtol=1e3*eps(BigFloat))
+@test isapprox(svdvals(bX), svdvals(X), rtol=1e3*eps())
+@test bX == X # check we didn't modify the input
+
+X = Float64[0 2 0; 0 1 2; 0 0 1]
+bX = big(X)
+bS = svdfact(bX)
+@test isapprox(full(bS), bX, rtol=1e3*eps(BigFloat))
+@test isapprox(svdvals(bX), svdvals(X), rtol=1e3*eps())
+@test bX == X # check we didn't modify the input
+
+
+bD = big(randn(m))
+bX = diagm(bD)
+bS = svdfact(bX)
+
+@test isapprox(full(bS), bX, rtol=1e3*eps(BigFloat))
+@test bS.S == sort(abs(bD),rev=true)
+
+
+
 
 X = randn(n,m)+im*randn(n,m)
 bX = big(X)
 bS = svdfact(bX)
 @test isapprox(full(bS), bX, rtol=1e3*eps(BigFloat))
 @test isapprox(svdvals(bX), svdvals(X), rtol=1e3*eps())
+@test bX == X # check we didn't modify the input
 
+bXt = bX'
+bSt = svdfact(bXt)
+@test isapprox(full(bSt), bXt, rtol=1e3*eps(BigFloat))
+@test isapprox(svdvals(bXt), svdvals(X), rtol=1e3*eps())
+@test bXt == X' # check we didn't modify the input