Basic implementation of GK SVD algorithm

Author: simonbyrne

LICENSE.md

@@ -1,6 +1,6 @@
 The GenericSVD.jl package is licensed under the MIT "Expat" License:
 
-> Copyright (c) 2016: Simon Byrne.
+> Copyright (c) 2016: Simon Byrne, Andreas Noack.
 >
 > Permission is hereby granted, free of charge, to any person obtaining
 > a copy of this software and associated documentation files (the

README.md

@@ -1,3 +1,7 @@
-# GenericSVD
+# GenericSVD.jl
 
 [![Build Status](https://travis-ci.org/simonbyrne/GenericSVD.jl.svg?branch=master)](https://travis-ci.org/simonbyrne/GenericSVD.jl)
+
+Implements Singular Value Decomposition for generic number types, such as `BigFloat` and `Complex{BigFloat}`. It internally overloads several Base functions such that existing methods (`svd`, `svdfact` and `svdvals`) should work directly.
+
+It uses a Golub-Kahan 2-stage algorithm of bidiagonalization with Householder reflections, followed by an implicit QR with shift.

src/GenericSVD.jl

@@ -1,5 +1,245 @@
 module GenericSVD
 
-# package code goes here
+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.")
+    B,P = bidiagonalize_tall!(X)
+    U,Vt = full(P,thin=thin)
+    U,S,Vt = svd!(B,U,Vt)
+    for i = 1:n
+        if signbit(S[i])
+            S[i] = -S[i]
+            for j = 1:n
+                Vt[i,j] = -Vt[i,j]
+            end
+        end
+    end
+    if sorted
+        I = sortperm(S,rev=true)
+        S = S[I]
+        U = U[:,I]
+        Vt = Vt[I,:]
+    end
+    SVD(U,S,Vt)
+end
+
+function svdvals!(X; sorted=true)
+    B,P = bidiagonalize_tall!(copy(X))
+    S = svd!(B)[2]
+    for i = eachindex(S)
+        if signbit(S[i])
+            S[i] = -S[i]
+        end
+    end
+    sorted ? sort!(S,rev=true) : S
+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}|)
+```
+"""
+offdiag_approx_zero(B::Bidiagonal,i,ɛ) =
+    abs(B.ev[i-1]) < ɛ*(abs(B.dv[i-1]) + abs(B.dv[i]))
+
+
+"""
+This finds the lowest strictly-bidiagonal submatrix, i.e. n₁, n₂ such that
+
+```
+     [ d ?           ]
+     [   d 0         ]
+  n₁ [     d e       ]
+     [       d e     ]
+  n₂ [         d 0   ]
+     [           d 0 ]
+```
+"""
+function svd!{T<:Real}(B::Bidiagonal{T}, U=nothing, Vt=nothing, ɛ::T = eps(T))
+    n = size(B, 1)
+    n₂ = n
+
+    if istriu(B)
+        while true
+            while offdiag_approx_zero(B,n₂,ɛ)
+                n₂ -= 1
+                if n₂ == 1
+                    return U,B.dv,Vt
+                end
+            end
+            n₁ = n₂ - 1
+            while n₁ > 1 && !offdiag_approx_zero(B,n₁,ɛ)
+                n₁ -= 1
+            end
+
+
+
+            # TODO: check for diagonal zeros
+            # if n₂ == n₁ + 1 # 2x2 block
+                
+            #     B.dv[n₁] = s₂
+            #     B.dv[n₂] = s₁
+            #     B.eb[n₁] = 0
+                
+            #     # TODO: apply op
+            # end
+
+            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
+            h = hypot(d₂,e)
+            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
+end
+
+
+
+"""
+Applies a Golub-Kahan SVD step.
+
+A Givens rotation is applied to the top 2x2 matrix, and the resulting "bulge" is "chased" down the diagonal to the bottom of the matrix.
+"""
+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
+        #  [b ,d₂]   [0  ,d₂′]
+
+
+        d₁ =  d₁′*G.c + e₁′*G.s
+        e₁ = -d₁′*G.s + e₁′*G.c
+        b  =  d₂′*G.s
+        d₂ =  d₂′*G.c
+
+        for i = n₁:n₂-2
+
+            #  [. ,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
+            
+            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
+            #  [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
+            
+            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₁] 
+        #  [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.ev[n₂-1] =  G.c*e₁ + G.s*d₂
+        B.dv[n₂]   = -G.s*e₁ + G.c*d₂
+    else
+        throw(ArgumentError("lower bidiagonal version not implemented yet"))
+    end
+
+    return B,U,Vt
+end
+
+
+"""
+The singular values of the matrix
+```
+B = [ f g ;
+      0 h ]
+```
+(i.e. the sqrt-eigenvalues of `B'*B`).
+
+This is a direct translation of LAPACK [DLAS2](http://www.netlib.org/lapack/explore-html/d8/dfd/dlas2_8f.html).
+"""
+function svdvals2x2(f, h, g)
+    fa = abs(f)
+    ga = abs(g)
+    ha = abs(h)
+
+    fhmin = min(fa,ha)
+    fhmax = max(fa,ha)
+    
+    if fhmin == 0
+        ssmin = zero(f)
+        if fhmax == 0
+            ssmax = zero(f)
+        else
+            ssmax = max(fhmax,ga)*sqrt(1+(min(fhmax,ga)/max(fhmax,ga))^2)
+        end
+    else
+        if ga < fhmax
+            as = 1 + fhmin/fhmax
+            at = (fhmax-fhmin)/fhmax
+            au = (ga/fhmax)^2
+            c = 2/(sqrt(as^2 + au) + sqrt(at^2+au))
+            ssmin = fhmin*c
+            ssmax = fhmax/c
+        else
+            au = fhmax / ga
+            if au == 0
+                ssmin = (fhmin*fhmax)/ga
+                ssmax = ga
+            else
+                as = 1+fhmin/fhmax
+                at = (fhmax-fhmin)/fhmax
+                c = 1/(sqrt(1 + (as*au)^2) + sqrt(1 + (at*au)^2))
+                ssmin = 2((fhmin*c)*au)
+                ssmax = ga/(2c)
+            end
+        end
+    end
+    ssmin,ssmax
+end
+
+
+
 
 end # module

src/bidiagonalize.jl

@@ -0,0 +1,68 @@
+"""
+Packed storage of bidiagonalizing QR reflectors.
+"""
+immutable PackedUVt{T}
+    A::Matrix{T}
+end
+
+
+"""
+Bidiagonalize a tall matrix `A` into `B`. Both arguments are overwritten.
+"""
+function bidiagonalize_tall!{T}(A::Matrix{T},B::Bidiagonal)
+    m, n = size(A)
+    # tall case: assumes m >= n
+    # upper bidiagonal
+
+    for i = 1:n
+        x = slice(A, i:m, i)
+        τi = LinAlg.reflector!(x)
+        B.dv[i] = real(A[i,i])
+        LinAlg.reflectorApply!(x, τi, slice(A, i:m, i+1:n))
+        A[i,i] = τi # store reflector in diagonal coefficient
+
+        # for Real, this only needs to be n-2
+        # needed for Complex to ensure superdiagonal is Real
+        if i <= n-1
+            x = slice(A, i, i+1:n)
+            conj!(x)
+            τi = LinAlg.reflector!(x)
+            B.ev[i] = real(A[i,i+1])
+            LinAlg.reflectorApply!(slice(A, i+1:m, i+1:n), x, τi)
+            A[i,i+1] = τi
+        end
+    end
+    B.isupper = true
+
+    B, PackedUVt(A)
+end
+
+function bidiagonalize_tall!{T}(A::Matrix{T})
+    m,n = size(A)
+    R = real(T)
+    B = Bidiagonal(Array(R,n),Array(R,n-1),true)
+    bidiagonalize_tall!(A,B)
+end
+
+function Base.full{T}(P::PackedUVt{T};thin=true)
+    A = P.A
+    m,n = size(A)
+
+    # U = Q_1 ... Q_n I_{m,n}
+    w = thin ? n : m
+    U = eye(T,m,w)
+    for i = n:-1:1
+        τi = A[i,i]
+        x = slice(A, i:m, i)
+        LinAlg.reflectorApply!(x, τi', slice(U, i:m, i:w))
+    end
+
+    # Vt = P_{n_2} ... P_1
+    Vt = eye(T,n,n)
+    for i = n-1:-1:1
+        τi = A[i,i+1]
+        x = slice(A, i, i+1:n)
+        LinAlg.reflectorApply!(slice(Vt, i:n, i+1:n), x, τi')
+    end
+    U,Vt
+end

src/utils.jl

@@ -0,0 +1,33 @@
+import Base.LinAlg: reflectorApply!
+@inline function reflectorApply!(A::StridedMatrix, x::AbstractVector, τ::Number) # apply conjugate transpose reflector from right.
+    m, n = size(A)
+    if length(x) != n
+        throw(DimensionMismatch("reflector must have same length as second dimension of matrix"))
+    end
+    @inbounds begin
+        for i = 1:m
+            Aiv = A[i, 1]
+            for j = 2:n
+                Aiv += A[i, j]*x[j]
+            end
+            Aiv = Aiv*τ
+            A[i, 1] -= Aiv
+            for j = 2:n
+                A[i, j] -= Aiv*x[j]'
+            end
+        end
+    end
+    return A
+end
+
+
+import Base: A_mul_B!, A_mul_Bc!, A_ldiv_B!
+
+A_mul_B!(G::LinAlg.Givens, ::Void) = nothing
+A_mul_Bc!(::Void, G::LinAlg.Givens) = nothing
+
+function A_ldiv_B!{Ta,Tb}(A::SVD{Ta}, B::StridedVecOrMat{Tb})
+    T = promote_type(Ta,Tb)
+    k = length(find(A.S .> eps(real(T))*maximum(A.S)))
+    A.Vt[1:k,:]' * (A.S[1:k] .\ (A.U[:,1:k]' * B))
+end