Merge pull request #1 from simonbyrne/sb/0.5

Update for 0.5, clean up some docs

Author: simonbyrne

.travis.yml

@@ -4,8 +4,9 @@ os:
   - linux
   - osx
 julia:
-  - release
-  - nightly
+  - 0.4
+  - 0.5
+#  - nightly
 notifications:
   email: false
 # uncomment the following lines to override the default test script

REQUIRE

@@ -1 +1,2 @@
 julia 0.4
+Compat

appveyor.yml

@@ -2,8 +2,10 @@ environment:
   matrix:
   - JULIAVERSION: "julialang/bin/winnt/x86/0.4/julia-0.4-latest-win32.exe"
   - JULIAVERSION: "julialang/bin/winnt/x64/0.4/julia-0.4-latest-win64.exe"
-  - JULIAVERSION: "julianightlies/bin/winnt/x86/julia-latest-win32.exe"
-  - JULIAVERSION: "julianightlies/bin/winnt/x64/julia-latest-win64.exe"
+  - JULIAVERSION: "julialang/bin/winnt/x86/0.4/julia-0.5-latest-win32.exe"
+  - JULIAVERSION: "julialang/bin/winnt/x64/0.4/julia-0.5-latest-win64.exe"
+#  - JULIAVERSION: "julianightlies/bin/winnt/x86/julia-latest-win32.exe"
+#  - JULIAVERSION: "julianightlies/bin/winnt/x64/julia-latest-win64.exe"
 
 branches:
   only:

src/GenericSVD.jl

@@ -1,5 +1,7 @@
+__precompile__(true)
 module GenericSVD
 
+import Compat: view
 import Base: SVD
 
 include("utils.jl")
@@ -53,12 +55,15 @@ end
 
 
 """
-Tests if the B[i-1,i] element is approximately zero, using the criteria
+    is_offdiag_approx_zero!(B::Bidiagonal,i,ɛ)
+
+Tests if the element `B[i-1,i]` is approximately zero, using the criteria
 ```math
     |B_{i-1,i}| ≤ ɛ*(|B_{i-1,i-1}| + |B_{i,i}|)
 ```
+If true, sets the element to exact zero.
 """
-function offdiag_approx_zero(B::Bidiagonal,i,ɛ)
+function is_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
@@ -68,9 +73,19 @@ end
 
 
 """
-Generic SVD algorithm:
+    svd!(B::Bidiagonal [, U, Vt [, ϵ]])
+
+Compute the SVD of a bidiagonal matrix `B`, via an implicit QR algorithm with shift (known as a Golub-Kahan iterations).
+
+Optional arguments:
+
+ * `U` and `Vt`: orthogonal matrices which pre- and post-multiply `B`, for computing the SVD of a full matrix `X = U*B*Vt` from its bidiagonalized form.
 
-This finds the lowest strictly-bidiagonal submatrix, i.e. n₁, n₂ such that
+ * `ϵ`: the tolerance for testing zeros of the offdiagonal elements of `B` (see below).
+
+Algorithm:
+
+This proceeds by iteratively finding the lowest strictly-bidiagonal submatrix, i.e. n₁, n₂ such that
 ```
      [ d ?           ]
      [   d 0         ]
@@ -79,7 +94,7 @@ This finds the lowest strictly-bidiagonal submatrix, i.e. n₁, n₂ such that
   n₂ [         d 0   ]
      [           d 0 ]
 ```
-Then applies a Golub-Kahan iteration.
+then applying a Golub-Kahan QR iteration.
 """
 function svd!{T<:Real}(B::Bidiagonal{T}, U=nothing, Vt=nothing, ɛ::T = eps(T))
     n = size(B, 1)
@@ -91,28 +106,25 @@ function svd!{T<:Real}(B::Bidiagonal{T}, U=nothing, Vt=nothing, ɛ::T = eps(T))
         while true
             @label mainloop
 
-            while offdiag_approx_zero(B,n₂,ɛ)
+            while is_offdiag_approx_zero!(B,n₂,ɛ)
                 n₂ -= 1
                 if n₂ == 1
                     @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₁ = n₂
+            while true
                 n₁ -= 1
+
                 # check for diagonal zeros
                 if abs(B.dv[n₁]) ≤ ɛ*maxB
                     svd_zerodiag_row!(U,B,n₁,n₂)
                     @goto mainloop
                 end
+                if n₁ == 1 || is_offdiag_approx_zero!(B,n₁,ɛ)
+                    break
+                end
             end
 
             if abs(B.dv[n₂]) ≤ ɛ*maxB
@@ -140,7 +152,9 @@ end
 
 
 """
-Sets B[n₁,n₁] to zero, then zeros out row n₁ by applying sequential row (left) Givens rotations up to n₂.
+    svd_zerodiag_row!(U,B,n₁,n₂)
+
+Sets `B[n₁,n₁]` to zero, then zeros out row `n₁` by applying sequential row (left) Givens rotations up to `n₂`, and the corresponding inverse rotations to `U` (preserveing `U*B`.
 """
 function svd_zerodiag_row!(U,B,n₁,n₂)
     e = B.ev[n₁]
@@ -167,7 +181,9 @@ end
 
 
 """
-Sets B[n₂,n₂] to zero, then zeros out column n₂ by applying sequential column (right) Givens rotations up to n₁.
+    svd_zerodiag_col!(B::Bidiagonal,Vt,n₁,n₂)
+
+Sets `B[n₂,n₂]` to zero, then zeros out column `n₂` by applying sequential column (right) Givens rotations up to `n₁`, and the corresponding inverse rotations to `Vt` (preserving `B*Vt`).
 """
 function svd_zerodiag_col!(B,Vt,n₁,n₂)
     e = B.ev[n₂-1]
@@ -187,7 +203,7 @@ function svd_zerodiag_col!(B,Vt,n₁,n₂)
 
         if n₁ < i
             eᵢ = B.ev[i-1]
-            e       = -G.s*eᵢ
+            e         = -G.s*eᵢ
             B.ev[i-1] = G.c*eᵢ
         end
     end
@@ -196,7 +212,9 @@ end
 
 
 """
-Applies a Golub-Kahan SVD step.
+    svd_gk!{T<:Real}(B::Bidiagonal{T},U,Vt,n₁,n₂,shift)
+
+Applies a Golub-Kahan SVD step (an implicit QR with shift) to the submatrix `B[n₁:n₂,n₁:n₂]`, applying the inverse transformations to `U` and `Vt` (preserving `U*B*Vt`).
 
 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.
 """
@@ -275,6 +293,8 @@ end
 
 
 """
+    svdvals2x2(f, h, g)
+
 The singular values of the matrix
 ```
 B = [ f g ;
@@ -289,8 +309,7 @@ function svdvals2x2(f, h, g)
     ga = abs(g)
     ha = abs(h)
 
-    fhmin = min(fa,ha)
-    fhmax = max(fa,ha)
+    fhmin, fhmax = minmax(fa,ha)
 
     if fhmin == 0
         ssmin = zero(f)

src/bidiagonalize.jl

@@ -15,20 +15,20 @@ function bidiagonalize_tall!{T}(A::Matrix{T},B::Bidiagonal)
     # upper bidiagonal
 
     for i = 1:n
-        x = slice(A, i:m, i)
+        x = view(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))
+        LinAlg.reflectorApply!(x, τi, view(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)
+            x = view(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)
+            LinAlg.reflectorApply!(view(A, i+1:m, i+1:n), x, τi)
             A[i,i+1] = τi
         end
     end
@@ -53,16 +53,16 @@ function Base.full{T}(P::PackedUVt{T};thin=true)
     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))
+        x = view(A, i:m, i)
+        LinAlg.reflectorApply!(x, τi', view(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')
+        x = view(A, i, i+1:n)
+        LinAlg.reflectorApply!(view(Vt, i:n, i+1:n), x, τi')
     end
     U,Vt
 end

src/utils.jl

@@ -26,9 +26,11 @@ 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})
-    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))
+if VERSION < v"0.5.0-"
+    function A_ldiv_B!{Ta,Tb}(A::SVD{Ta}, B::StridedVecOrMat{Tb})
+        k = searchsortedlast(A.S, eps(real(Ta))*A.S[1], rev=true)
+        view(A.Vt,1:k,:)' * (view(A.S,1:k) .\ (view(A.U,:,1:k)' * B))
+    end
 end