Repository: JuliaLinearAlgebra/GenericLinearAlgebra.jl (julia)

Issue:
svd not supported for lower bidiagonal matrices
As reported in https://github.com/JuliaLinearAlgebra/TSVD.jl/issues/33 and fixed by #156. However, I think this should be done with a pass of Givens rotations to flip the Bidiagonal. Currently, the error is

```julia
julia> GenericLinearAlgebra._svd!(Bidiagonal(randn(4), randn(3), :L))
ERROR: SVD for lower triangular bidiagonal matrices isn't implemented yet.
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] __svd!(B::Bidiagonal{Float64, Vector{Float64}}, U::Matrix{Float64}, Vᴴ::Matrix{Float64}; tol::Float64)
   @ GenericLinearAlgebra ~/.julia/dev/GenericLinearAlgebra/src/svd.jl:252
 [3] __svd!
   @ ~/.julia/dev/GenericLinearAlgebra/src/svd.jl:156 [inlined]
 [4] #_svd!#40
   @ ~/.julia/dev/GenericLinearAlgebra/src/svd.jl:292 [inlined]
 [5] _svd!(B::Bidiagonal{Float64, Vector{Float64}})
   @ GenericLinearAlgebra ~/.julia/dev/GenericLinearAlgebra/src/svd.jl:286
 [6] top-level scope
   @ REPL[17]:1
```

Pull request:
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.

Supersedes #156. Closes #157.

Patch:
diff --git a/Project.toml b/Project.toml
index baffda4..a5ff553 100644
--- a/Project.toml
+++ b/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"
diff --git a/src/svd.jl b/src/svd.jl
index 8af3949..67259a3 100644
--- a/src/svd.jl
+++ b/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 patch:
diff --git a/test/svd.jl b/test/svd.jl
index bf211b4..9bcecfe 100644
--- a/test/svd.jl
+++ 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