Add bidiagonal cholesky version (#1409)

Co-authored-by: Mateus Araújo <[email protected]>

Author: dkarrasch

src/bidiag.jl

@@ -1467,6 +1467,35 @@ function inv(B::Bidiagonal{T}) where T
     return B.uplo == 'U' ? UpperTriangular(dest) : LowerTriangular(dest)
 end
 
+# cholesky-version for (sym)tridiagonal matrices
+for (T, uplo) in ((:UpperTriangular, :(:U)), (:LowerTriangular, :(:L)))
+    @eval function _chol!(A::Bidiagonal, ::Type{$T})
+        dv = real(A.dv)
+        ev = A.ev
+        n = length(dv)
+        @inbounds for i in 1:n-1
+            iszero(dv[i]) && throw(ZeroPivotException(i))
+            ev[i] /= dv[i]
+            dv[i+1] -= abs2(ev[i])*dv[i]
+            Akk = dv[i]
+            Akk, info = _chol!(Akk, UpperTriangular)
+            if info != 0
+                return $T(A), convert(BlasInt, i)
+            end
+            dv[i] = Akk
+            ev[i] *= dv[i]
+        end
+        Akk = dv[n]
+        Akk, info = _chol!(Akk, $T)
+        if info != 0
+            return $T(A), convert(BlasInt, n)
+        end
+        dv[n] = Akk
+        B = Bidiagonal(dv, ev, $uplo)
+        return $T(B), convert(BlasInt, 0)
+    end
+end
+
 # Eigensystems
 eigvals(M::Bidiagonal) = copy(M.dv)
 function eigvecs(M::Bidiagonal{T}) where T

src/special.jl

@@ -605,8 +605,8 @@ function _hvcat(rows::Tuple{Vararg{Int}}, A::Union{AbstractArray,AbstractQ,Unifo
     end
 end
 
-# factorizations
-function cholesky(S::RealHermSymComplexHerm{<:Real,<:SymTridiagonal}, ::NoPivot = NoPivot(); check::Bool = true)
+# tridiagonal cholesky factorization
+function cholesky(S::RealSymHermitian{<:BiTriSym}, ::NoPivot = NoPivot(); check::Bool = true)
     T = choltype(S)
     B = Bidiagonal{T}(diag(S, 0), diag(S, S.uplo == 'U' ? 1 : -1), sym_uplo(S.uplo))
     cholesky!(Hermitian(B, sym_uplo(S.uplo)), NoPivot(); check = check)

src/tridiag.jl

@@ -1039,7 +1039,7 @@ function dot(x::AbstractVector, A::Tridiagonal, y::AbstractVector)
     return r
 end
 
-function cholesky(S::SymTridiagonal, ::NoPivot = NoPivot(); check::Bool = true)
+function cholesky(S::Union{SymTridiagonal,Tridiagonal}, ::NoPivot = NoPivot(); check::Bool = true)
     if !ishermitian(S)
         check && checkpositivedefinite(-1)
         return Cholesky(S, 'U', convert(BlasInt, -1))

test/cholesky.jl

@@ -492,9 +492,41 @@ end
 end
 
 @testset "Cholesky for AbstractMatrix" begin
-    S = SymTridiagonal(fill(2.0, 4), ones(3))
-    C = cholesky(S)
-    @test C.L * C.U ≈ S
+    for T in (identity, big), M in (SymTridiagonal(fill(T(2.0), 4), ones(3)),
+        Symmetric(SymTridiagonal(fill(T(2.0), 4), ones(3)), :U),
+        Symmetric(SymTridiagonal(fill(T(2.0), 4), ones(3)), :L),
+        Tridiagonal(ones(3), fill(T(2.0), 4), ones(3)),
+        Hermitian(Tridiagonal(ones(3), fill(T(2.0), 4), ones(3)), :U),
+        Hermitian(Bidiagonal(fill(T(2.0), 4), ones(3), :U), :U),
+        Hermitian(Bidiagonal(fill(T(2.0), 4), ones(3), :U), :L),
+        Hermitian(Bidiagonal(fill(T(2.0), 4), ones(3), :L), :L),
+        )
+        C = cholesky(M)
+        @test C.L * C.U ≈ M
+        @test parent(C.U) isa Bidiagonal
+        M[1,1] *= -1
+        @test_throws PosDefException cholesky(M)
+        C = cholesky(M, check=false)
+        @test C.info > 0
+        M[1,1] *= -1
+        M[end,end] *= -1
+        @test_throws PosDefException cholesky(M)
+        C = cholesky(M, check=false)
+        @test C.info > 0
+    end
+    # test LowerTriangular version
+    M = Hermitian(Bidiagonal(fill(2.0, 4), im * ones(3), :L), :L)
+    C = cholesky!(copy(M))
+    @test C.L * C.U ≈ M
+    # non-(RealOrComplex) eltype
+    A = Tridiagonal(randn(Quaternion{Float64}, 4, 4) |> t -> t't)
+    C = cholesky(A)
+    @test C.L * C.U ≈ A
+    @test parent(C.U) isa Bidiagonal
+    A = Hermitian(Tridiagonal(randn(Quaternion{Float64}, 4, 4) |> t -> t't), :L)
+    C = cholesky(A)
+    @test C.L * C.U ≈ A
+    @test parent(C.U) isa Bidiagonal
 end
 
 @testset "constructor with non-BlasInt arguments" begin