`ldlt` failes to check whether input matrix is positive-definite

[runapp]

Julia Version 1.1.0
Commit 80516ca202 (2019-01-21 21:24 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i5-8400 CPU @ 3.40GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, znver1)

We expect both ldlt calls raise a PosDefException. However, only the second call does so.

using LinearAlgebra,SparseArrays

A=SymTridiagonal([1,0,0],[0,1])
A2=Symmetric(sparse(A))

t1=ldlt(A)
t2=ldlt(A2)

[runapp]

Further investigation shows the first ldlt is from LinearAlgebra and second from SparseSuite.CHOLMOD. The first one is pure julia and second is a wrapper for C code. Pos-def check is done at cholmod.jl:1448 while no check in LinearAlgebra.ldlt. Should we add a check for it?

[stevengj]

ldlt(A) in your example doesn't throw, but it produces Inf and NaN data, so I agree that it would be better to throw.

In principle, the LDLᵀ decomposition (unlike ordinary Cholesky) doesn't require a positive-definite matrix, but I guess we don't currently support that case? It would be nicer to just fix the code to support indefinite symmetric matrices.

[andreasnoack]

Our versions don't require the matrix to be PD, e.g.

julia> T = SymTridiagonal([1.0,-1,1], ones(2)/2)
3×3 SymTridiagonal{Float64,Array{Float64,1}}:
 1.0   0.5   ⋅
 0.5  -1.0  0.5
  ⋅    0.5  1.0

julia> eigvals(T)
3-element Array{Float64,1}:
 -1.2247448713915894
  1.0
  1.224744871391589

julia> ldlt(T)
LDLt{Float64,SymTridiagonal{Float64,Array{Float64,1}}}([1.0 0.5 0.0; 0.5 -1.25 -0.4; 0.0 -0.4 1.2])

but both versions break down for zero pivots because they don't pivot based on the values of the matrix. The sparse version from CHOLMOD does some pivoting during the symbolic factorization to reduce fill but no pivoting during the numerical factorization and the Julia version for SymTridiagonal doesn't do pivoting at all.

[stevengj]

Doesn't that also mean that ldlt is numerically unstable for indefinite matrices (arbitrarily bad roundoff errors for near-zero pivots)?

[andreasnoack]

Doesn't that also mean that ldlt is numerically unstable for indefinite matrices (arbitrarily bad roundoff errors for near-zero pivots)?

I believe it does (but that is also the case for LU with partial pivoting). It would be good to have pivoted version of these factorizations.

[stevengj]

but that is also the case for LU with partial pivoting

Not sure what you mean here — LU with partial pivoting is backwards stable. (In principle, the constant factor in the stability condition is exponentially large in the matrix size, but in practice this never seems to apply to realistic matrices. Regardless, it is technically stable.) Without pivoting, on the other hand, it is unstable (in both theory and practice).

[ViralBShah]

julia> t1=ldlt(A)
ERROR: ZeroPivotException: factorization encountered one or more zero pivots. Consider switching to a pivoted LU factorization.
Stacktrace:
 [1] ldlt!(S::SymTridiagonal{Float64, Vector{Float64}})
   @ LinearAlgebra ~/julia/usr/share/julia/stdlib/v1.8/LinearAlgebra/src/ldlt.jl:122
 [2] ldlt(M::SymTridiagonal{Int64, Vector{Int64}}; shift::Bool)
   @ LinearAlgebra ~/julia/usr/share/julia/stdlib/v1.8/LinearAlgebra/src/ldlt.jl:167
 [3] ldlt(M::SymTridiagonal{Int64, Vector{Int64}})
   @ LinearAlgebra ~/julia/usr/share/julia/stdlib/v1.8/LinearAlgebra/src/ldlt.jl:162
 [4] top-level scope
   @ REPL[33]:1

julia> t2=ldlt(A2)
ERROR: ZeroPivotException: factorization encountered one or more zero pivots. Consider switching to a pivoted LU factorization.
Stacktrace:
 [1] #ldlt!#10
   @ ~/julia/usr/share/julia/stdlib/v1.8/SuiteSparse/src/cholmod.jl:1308 [inlined]
 [2] ldlt(A::SuiteSparse.CHOLMOD.Sparse{Float64}; shift::Float64, check::Bool, perm::Nothing)
   @ SuiteSparse.CHOLMOD ~/julia/usr/share/julia/stdlib/v1.8/SuiteSparse/src/cholmod.jl:1347
 [3] ldlt
   @ ~/julia/usr/share/julia/stdlib/v1.8/SuiteSparse/src/cholmod.jl:1340 [inlined]
 [4] #ldlt#13
   @ ~/julia/usr/share/julia/stdlib/v1.8/SuiteSparse/src/cholmod.jl:1396 [inlined]
 [5] ldlt(A::Symmetric{Int64, SparseMatrixCSC{Int64, Int64}})
   @ SuiteSparse.CHOLMOD ~/julia/usr/share/julia/stdlib/v1.8/SuiteSparse/src/cholmod.jl:1396
 [6] top-level scope
   @ REPL[34]:1