Convert lower to upper bidgiagonal in _svd! and _svdvals! (#158)

This is done by running a single pass of givens rotations from the
left while also updating the U matrix. For svdvals!, the uplo field
is simply flipped as the singular values of the transpose are the
same.

Author: andreasnoack

Project.toml

@@ -1,6 +1,6 @@
 name = "GenericLinearAlgebra"
 uuid = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
-version = "0.3.17"
+version = "0.3.18"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/svd.jl

@@ -71,7 +71,7 @@ function svdIter!(
         d[n2] = -G.s * ei1 + G.c * di1
 
     else
-        throw(ArgumentError("lower bidiagonal version not implemented yet"))
+        throw(ArgumentError("please convert to upper bidiagonal"))
     end
 
     return B
@@ -254,6 +254,7 @@ function __svd!(
 end
 
 function _svdvals!(B::Bidiagonal{T}; tol = eps(T)) where {T<:Real}
+    B = Bidiagonal(B.dv, B.ev, :U)
     __svd!(B, tol = tol)
     return sort(abs.(diag(B)), rev = true)
 end
@@ -283,12 +284,31 @@ function _sort_and_adjust!(U, s, Vᴴ)
     return nothing
 end
 
+# Convert a lower Bidiagonal to an upper by running a single pass
+# of givens rotations from the left. A noop if the Bidiagnoal is
+# already upper.
+function _lower_to_upper!(B::Bidiagonal, U::Matrix)
+    if istriu(B)
+        return B
+    end
+    for i in 1:length(B.ev)
+        G, r = givens(B.dv[i], B.ev[i], i, i + 1)
+        B.dv[i] = r
+        B.ev[i] = G.s * B.dv[i + 1]
+        B.dv[i + 1] = G.c * B.dv[i + 1]
+        rmul!(U, G')
+    end
+    return Bidiagonal(B.dv, B.ev, :U)
+end
+
 function _svd!(B::Bidiagonal{T}; tol = eps(T)) where {T<:Real}
     n = size(B, 1)
 
     U = Matrix{T}(I, n, n)
     Vᴴ = Matrix{T}(I, n, n)
 
+    B = _lower_to_upper!(B, U)
+
     __svd!(B, U, Vᴴ, tol = tol)
 
     s = diag(B)

test/svd.jl

@@ -266,4 +266,12 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions, DoubleFloats
         @test svdvals(BigFloat[0 0; 1 -1]) ≈ [sqrt(2), 0]
         @test svdvals(BigFloat[1 0 0; 0 0 0; 0 1 -1]) ≈ [sqrt(2), 1, 0]
     end
+
+    @testset "Lower Bidiagonal matrices. Issue 157" begin
+        B = Bidiagonal(randn(4), randn(3), :L)
+        @test GenericLinearAlgebra._svdvals!(copy(B)) ≈ GenericLinearAlgebra._svdvals!(copy(B'))
+        U, s, V = GenericLinearAlgebra._svd!(copy(B))
+        @test B ≈ U * Diagonal(s) * V'
+        @test_throws ArgumentError("please convert to upper bidiagonal") GenericLinearAlgebra.svdIter!(B, 1, size(B, 1), 1.0)
+    end
 end