From 3d9121c644443e69d49d2344eeb35f5c57d7fcc0 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Mon, 12 Jun 2017 17:00:32 -0400 Subject: [PATCH] Delay throwing when CHOLMOD factorizations fail Introduce issuccess function to test if a Factorization succeeded. Fixes #22335 --- base/linalg/bunchkaufman.jl | 17 +++++++--- base/linalg/cholesky.jl | 12 +++++-- base/linalg/factorization.jl | 20 +++++++++++ base/linalg/linalg.jl | 1 + base/linalg/lu.jl | 35 +++++++++++--------- base/sparse/cholmod.jl | 64 +++++++++++++++++++++--------------- base/sparse/linalg.jl | 22 +++++++------ doc/src/stdlib/linalg.md | 1 + test/linalg/bunchkaufman.jl | 3 ++ test/linalg/cholesky.jl | 6 ++-- test/linalg/lu.jl | 8 +++-- test/sparse/cholmod.jl | 19 ++++++++--- 12 files changed, 140 insertions(+), 68 deletions(-) diff --git a/base/linalg/bunchkaufman.jl b/base/linalg/bunchkaufman.jl index 2161af6469bdc..3b674df501d44 100644 --- a/base/linalg/bunchkaufman.jl +++ b/base/linalg/bunchkaufman.jl @@ -88,8 +88,15 @@ size(B::BunchKaufman, d::Integer) = size(B.LD, d) issymmetric(B::BunchKaufman) = B.symmetric ishermitian(B::BunchKaufman) = !B.symmetric +issuccess(B::BunchKaufman) = B.info == 0 + +function Base.show(io::IO, F::BunchKaufman) + println(io, summary(F)) + print(io, "successful: $(issuccess(F))") +end + function inv(B::BunchKaufman{<:BlasReal}) - if B.info > 0 + if !issuccess(B) throw(SingularException(B.info)) end @@ -101,7 +108,7 @@ function inv(B::BunchKaufman{<:BlasReal}) end function inv(B::BunchKaufman{<:BlasComplex}) - if B.info > 0 + if !issuccess(B) throw(SingularException(B.info)) end @@ -121,7 +128,7 @@ function inv(B::BunchKaufman{<:BlasComplex}) end function A_ldiv_B!(B::BunchKaufman{T}, R::StridedVecOrMat{T}) where T<:BlasReal - if B.info > 0 + if !issuccess(B) throw(SingularException(B.info)) end @@ -132,7 +139,7 @@ function A_ldiv_B!(B::BunchKaufman{T}, R::StridedVecOrMat{T}) where T<:BlasReal end end function A_ldiv_B!(B::BunchKaufman{T}, R::StridedVecOrMat{T}) where T<:BlasComplex - if B.info > 0 + if !issuccess(B) throw(SingularException(B.info)) end @@ -161,7 +168,7 @@ function logabsdet(F::BunchKaufman) p = F.ipiv n = size(F.LD, 1) - if F.info > 0 + if !issuccess(F) return eltype(F)(-Inf), zero(eltype(F)) end s = one(real(eltype(F))) diff --git a/base/linalg/cholesky.jl b/base/linalg/cholesky.jl index 434a7fdbe5a5c..36de4a5007289 100644 --- a/base/linalg/cholesky.jl +++ b/base/linalg/cholesky.jl @@ -410,8 +410,13 @@ function getindex(C::CholeskyPivoted{T}, d::Symbol) where T<:BlasFloat throw(KeyError(d)) end -show(io::IO, C::Cholesky{<:Any,<:AbstractMatrix}) = - (println(io, "$(typeof(C)) with factor:");show(io,C[:UL])) +issuccess(C::Cholesky) = C.info == 0 + +function show(io::IO, C::Cholesky{<:Any,<:AbstractMatrix}) + println(io, "$(typeof(C)) with factor:") + show(io, C[:UL]) + print(io, "\nsuccessful: $(issuccess(C))") +end A_ldiv_B!(C::Cholesky{T,<:AbstractMatrix}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = @assertposdef LAPACK.potrs!(C.uplo, C.factors, B) C.info @@ -464,11 +469,12 @@ end isposdef(C::Cholesky) = C.info == 0 function det(C::Cholesky) + isposdef(C) || throw(PosDefException(C.info)) dd = one(real(eltype(C))) @inbounds for i in 1:size(C.factors,1) dd *= real(C.factors[i,i])^2 end - @assertposdef dd C.info + return dd end function logdet(C::Cholesky) diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl index ef5b69ad1575f..11efddb8191f9 100644 --- a/base/linalg/factorization.jl +++ b/base/linalg/factorization.jl @@ -16,6 +16,26 @@ macro assertnonsingular(A, info) :($(esc(info)) == 0 ? $(esc(A)) : throw(SingularException($(esc(info))))) end +""" + + issuccess(F::Factorization) + +Test that a factorization of a matrix succeeded. + +```jldoctest +julia> cholfact([1 0; 0 1]) +Base.LinAlg.Cholesky{Float64,Array{Float64,2}} with factor: +[1.0 0.0; 0.0 1.0] +successful: true + +julia> cholfact([1 0; 0 -1]) +Base.LinAlg.Cholesky{Float64,Array{Float64,2}} with factor: +[1.0 0.0; 0.0 -1.0] +successful: false +``` +""" +issuccess(F::Factorization) + function logdet(F::Factorization) d, s = logabsdet(F) return d + log(s) diff --git a/base/linalg/linalg.jl b/base/linalg/linalg.jl index ffc97a3dd78c4..1cc45399ccae5 100644 --- a/base/linalg/linalg.jl +++ b/base/linalg/linalg.jl @@ -91,6 +91,7 @@ export ishermitian, isposdef, isposdef!, + issuccess, issymmetric, istril, istriu, diff --git a/base/linalg/lu.jl b/base/linalg/lu.jl index 8d2558b9cb9a2..d4b7e7af9238b 100644 --- a/base/linalg/lu.jl +++ b/base/linalg/lu.jl @@ -147,7 +147,7 @@ function lufact(A::AbstractMatrix{T}) where T AA = similar(A, S, size(A)) copy!(AA, A) F = lufact!(AA, Val{false}) - if F.info == 0 + if issuccess(F) return F else AA = similar(A, S, size(A)) @@ -227,11 +227,14 @@ function getindex(F::LU{T,<:StridedMatrix}, d::Symbol) where T end end -function show(io::IO, C::LU) - println(io, "$(typeof(C)) with factors L and U:") - show(io, C[:L]) +issuccess(F::LU) = F.info == 0 + +function show(io::IO, F::LU) + println(io, "$(typeof(F)) with factors L and U:") + show(io, F[:L]) println(io) - show(io, C[:U]) + show(io, F[:U]) + print(io, "\nsuccessful: $(issuccess(F))") end A_ldiv_B!(A::LU{T,<:StridedMatrix}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = @@ -271,31 +274,31 @@ Ac_ldiv_Bc(A::LU{T,<:StridedMatrix}, B::StridedVecOrMat{T}) where {T<:BlasComple @assertnonsingular LAPACK.getrs!('C', A.factors, A.ipiv, ctranspose(B)) A.info Ac_ldiv_Bc(A::LU, B::StridedVecOrMat) = Ac_ldiv_B(A, ctranspose(B)) -function det(A::LU{T}) where T - n = checksquare(A) - A.info > 0 && return zero(T) +function det(F::LU{T}) where T + n = checksquare(F) + issuccess(F) || return zero(T) P = one(T) c = 0 @inbounds for i = 1:n - P *= A.factors[i,i] - if A.ipiv[i] != i - c += 1 + P *= F.factors[i,i] + if F.ipiv[i] != i + c += 1 end end s = (isodd(c) ? -one(T) : one(T)) return P * s end -function logabsdet(A::LU{T}) where T # return log(abs(det)) and sign(det) - n = checksquare(A) - A.info > 0 && return log(zero(real(T))), log(one(T)) +function logabsdet(F::LU{T}) where T # return log(abs(det)) and sign(det) + n = checksquare(F) + issuccess(F) || return log(zero(real(T))), log(one(T)) c = 0 P = one(T) abs_det = zero(real(T)) @inbounds for i = 1:n - dg_ii = A.factors[i,i] + dg_ii = F.factors[i,i] P *= sign(dg_ii) - if A.ipiv[i] != i + if F.ipiv[i] != i c += 1 end abs_det += log(abs(dg_ii)) diff --git a/base/sparse/cholmod.jl b/base/sparse/cholmod.jl index 08eec5cdac939..cfd67a1595575 100644 --- a/base/sparse/cholmod.jl +++ b/base/sparse/cholmod.jl @@ -7,7 +7,7 @@ import Base: (*), convert, copy, eltype, get, getindex, show, size, import Base.LinAlg: (\), A_mul_Bc, A_mul_Bt, Ac_ldiv_B, Ac_mul_B, At_ldiv_B, At_mul_B, cholfact, cholfact!, det, diag, ishermitian, isposdef, - issymmetric, ldltfact, ldltfact!, logdet + issuccess, issymmetric, ldltfact, ldltfact!, logdet importall ..SparseArrays @@ -1195,6 +1195,7 @@ function showfactor(io::IO, F::Factor) @printf(io, "method: %10s\n", s.is_super!=0 ? "supernodal" : "simplicial") @printf(io, "maxnnz: %10d\n", Int(s.nzmax)) @printf(io, "nnz: %13d\n", nnz(F)) + @printf(io, "success: %9s\n", "$(s.minor == size(F, 1))") end # getindex not defined for these, so don't use the normal array printer @@ -1356,8 +1357,6 @@ function cholfact!(F::Factor{Tv}, A::Sparse{Tv}; shift::Real=0.0) where Tv # Compute the numerical factorization factorize_p!(A, shift, F, cm) - s = unsafe_load(get(F.p)) - s.minor < size(A, 1) && throw(Base.LinAlg.PosDefException(s.minor)) return F end @@ -1397,8 +1396,6 @@ function cholfact(A::Sparse; shift::Real=0.0, # Compute the numerical factorization cholfact!(F, A; shift = shift) - s = unsafe_load(get(F.p)) - s.minor < size(A, 1) && throw(Base.LinAlg.PosDefException(s.minor)) return F end @@ -1443,13 +1440,17 @@ cholfact(A::Union{SparseMatrixCSC{T}, SparseMatrixCSC{Complex{T}}, function ldltfact!(F::Factor{Tv}, A::Sparse{Tv}; shift::Real=0.0) where Tv - cm = common() + cm = defaults(common()) + set_print_level(cm, 0) + + # Makes it an LDLt + unsafe_store!(common_final_ll, 0) + # Really make sure it's an LDLt by avoiding supernodal factorization + unsafe_store!(common_supernodal, 0) # Compute the numerical factorization factorize_p!(A, shift, F, cm) - s = unsafe_load(get(F.p)) - s.minor < size(A, 1) && throw(Base.LinAlg.ArgumentError("matrix has one or more zero pivots")) return F end @@ -1494,10 +1495,6 @@ function ldltfact(A::Sparse; shift::Real=0.0, # Compute the numerical factorization ldltfact!(F, A; shift = shift) - s = unsafe_load(get(F.p)) - if s.minor < size(A, 1) - throw(Base.LinAlg.ArgumentError("matrix has one or more zero pivots")) - end return F end @@ -1699,11 +1696,16 @@ for f in (:\, :Ac_ldiv_B) @eval function ($f)(A::Union{Symmetric{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}}, Hermitian{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}}, Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}}, B::StridedVecOrMat) - try - return ($f)(cholfact(A), B) - catch e - isa(e, LinAlg.PosDefException) || rethrow(e) - return ($f)(ldltfact(A) , B) + F = cholfact(A) + if issuccess(F) + return ($f)(F, B) + else + ldltfact!(F, A) + if issuccess(F) + return ($f)(F, B) + else + return ($f)(lufact(convert(SparseMatrixCSC{eltype(A), SuiteSparse_long}, A)), B) + end end end end @@ -1750,17 +1752,27 @@ end det(L::Factor) = exp(logdet(L)) -function isposdef(A::SparseMatrixCSC{<:VTypes,SuiteSparse_long}) - if !ishermitian(A) - return false - end - try - f = cholfact(A) - catch e - isa(e, LinAlg.PosDefException) || rethrow(e) +function issuccess(F::Factor) + s = unsafe_load(F.p) + return s.minor == size(F, 1) +end + +function isposdef(F::Factor) + if issuccess(F) + s = unsafe_load(F.p) + if s.is_ll == 1 + return true + else + # try conversion to LLt + change_factor!(eltype(F), true, s.is_super, true, s.is_monotonic, F) + b = issuccess(F) + # convert back + change_factor!(eltype(F), false, s.is_super, true, s.is_monotonic, F) + return b + end + else return false end - true end function ishermitian(A::Sparse{Float64}) diff --git a/base/sparse/linalg.jl b/base/sparse/linalg.jl index 23590f216bf60..030724d4dd791 100644 --- a/base/sparse/linalg.jl +++ b/base/sparse/linalg.jl @@ -904,19 +904,21 @@ function factorize(A::SparseMatrixCSC) end function factorize(A::Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}) where Ti - try - return cholfact(A) - catch e - isa(e, PosDefException) || rethrow(e) - return ldltfact(A) + F = cholfact(A) + if LinAlg.issuccess(F) + return F + else + ldltfact!(F, A) + return F end end function factorize(A::Hermitian{Complex{Float64}, SparseMatrixCSC{Complex{Float64},Ti}}) where Ti - try - return cholfact(A) - catch e - isa(e, PosDefException) || rethrow(e) - return ldltfact(A) + F = cholfact(A) + if LinAlg.issuccess(F) + return F + else + ldltfact!(F, A) + return F end end diff --git a/doc/src/stdlib/linalg.md b/doc/src/stdlib/linalg.md index 66dde6330d3e8..dcf9d47449f56 100644 --- a/doc/src/stdlib/linalg.md +++ b/doc/src/stdlib/linalg.md @@ -95,6 +95,7 @@ Base.LinAlg.logm Base.LinAlg.sqrtm Base.LinAlg.lyap Base.LinAlg.sylvester +Base.LinAlg.issuccess Base.LinAlg.issymmetric Base.LinAlg.isposdef Base.LinAlg.isposdef! diff --git a/test/linalg/bunchkaufman.jl b/test/linalg/bunchkaufman.jl index 233447c99f66a..42289b4a7732f 100644 --- a/test/linalg/bunchkaufman.jl +++ b/test/linalg/bunchkaufman.jl @@ -50,6 +50,7 @@ bimg = randn(n,2)/2 @testset "(Automatic) Bunch-Kaufman factor of indefinite matrix" begin bc1 = factorize(asym) + @test LinAlg.issuccess(bc1) @test logabsdet(bc1)[1] ≈ log(abs(det(bc1))) if eltya <: Real @test logabsdet(bc1)[2] == sign(det(bc1)) @@ -72,6 +73,7 @@ bimg = randn(n,2)/2 @testset "Bunch-Kaufman factors of a pos-def matrix" begin @testset for rook in (false, true) bc2 = bkfact(apd, :U, issymmetric(apd), rook) + @test LinAlg.issuccess(bc2) @test logdet(bc2) ≈ log(det(bc2)) @test logabsdet(bc2)[1] ≈ log(abs(det(bc2))) @test logabsdet(bc2)[2] == sign(det(bc2)) @@ -103,6 +105,7 @@ end @testset for rook in (false, true) F = bkfact(As, :U, issymmetric(As), rook) + @test !LinAlg.issuccess(F) @test det(F) == 0 @test_throws LinAlg.SingularException inv(F) @test_throws LinAlg.SingularException F \ ones(size(As, 1)) diff --git a/test/linalg/cholesky.jl b/test/linalg/cholesky.jl index 3152ca219475d..dce15ce0d776c 100644 --- a/test/linalg/cholesky.jl +++ b/test/linalg/cholesky.jl @@ -51,6 +51,7 @@ using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted, PosDefException @test E[i,j] <= (n+1)ε/(1-(n+1)ε)*real(sqrt(apd[i,i]*apd[j,j])) end @test apd*inv(capd) ≈ eye(n) + @test LinAlg.issuccess(capd) @test abs((det(capd) - det(apd))/det(capd)) <= ε*κ*n # Ad hoc, but statistically verified, revisit @test @inferred(logdet(capd)) ≈ log(det(capd)) # logdet is less likely to overflow @@ -77,7 +78,7 @@ using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted, PosDefException @test isposdef(capds) end ulstring = sprint(show,capds[:UL]) - @test sprint(show,capds) == "$(typeof(capds)) with factor:\n$ulstring" + @test sprint(show,capds) == "$(typeof(capds)) with factor:\n$ulstring\nsuccessful: true" else capdh = cholfact(apdh) @test inv(capdh)*apdh ≈ eye(n) @@ -95,7 +96,7 @@ using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted, PosDefException @test logdet(capdh) ≈ log(det(capdh)) @test isposdef(capdh) ulstring = sprint(show,capdh[:UL]) - @test sprint(show,capdh) == "$(typeof(capdh)) with factor:\n$ulstring" + @test sprint(show,capdh) == "$(typeof(capdh)) with factor:\n$ulstring\nsuccessful: true" end # test chol of 2x2 Strang matrix @@ -281,6 +282,7 @@ end A = T[1 2; 2 1]; B = T[1, 1] C = cholfact(A) @test !isposdef(C) + @test !LinAlg.issuccess(C) @test_throws PosDefException C\B @test_throws PosDefException det(C) @test_throws PosDefException logdet(C) diff --git a/test/linalg/lu.jl b/test/linalg/lu.jl index fd6a59621a9c0..e05b755909192 100644 --- a/test/linalg/lu.jl +++ b/test/linalg/lu.jl @@ -67,7 +67,7 @@ debug && println("(Automatic) Square LU decomposition. eltya: $eltya, eltyb: $el lstring = sprint(show,l) ustring = sprint(show,u) - @test sprint(show,lua) == "$(typeof(lua)) with factors L and U:\n$lstring\n$ustring" + @test sprint(show,lua) == "$(typeof(lua)) with factors L and U:\n$lstring\n$ustring\nsuccessful: true" let Bs = b, Cs = c for atype in ("Array", "SubArray") if atype == "Array" @@ -97,6 +97,7 @@ debug && println("(Automatic) Square LU decomposition. eltya: $eltya, eltyb: $el debug && println("Tridiagonal LU") κd = cond(Array(d),1) lud = lufact(d) + @test LinAlg.issuccess(lud) @test lufact(lud) == lud @test_throws KeyError lud[:Z] @test lud[:L]*lud[:U] ≈ lud[:P]*Array(d) @@ -139,7 +140,7 @@ debug && println("Tridiagonal LU") du[1] = zero(eltya) dl[1] = zero(eltya) zT = Tridiagonal(dl,dd,du) - @test lufact(zT).info == 1 + @test !LinAlg.issuccess(lufact(zT)) end debug && println("Thin LU") @@ -211,3 +212,6 @@ end # Issue 21453. @test_throws ArgumentError LinAlg._cond1Inf(lufact(randn(5,5)), 2, 2.0) + +# Singular LU +@test !LinAlg.issuccess(lufact(zeros(3,3))) diff --git a/test/sparse/cholmod.jl b/test/sparse/cholmod.jl index 014fd6317b38c..da632383f629d 100644 --- a/test/sparse/cholmod.jl +++ b/test/sparse/cholmod.jl @@ -348,14 +348,14 @@ for elty in (Float64, Complex{Float64}) # Factor @test_throws ArgumentError cholfact(A1) - @test_throws Base.LinAlg.PosDefException cholfact(A1 + A1' - 2eigmax(Array(A1 + A1'))I) - @test_throws Base.LinAlg.PosDefException cholfact(A1 + A1', shift=-2eigmax(Array(A1 + A1'))) - @test_throws ArgumentError ldltfact(A1 + A1' - 2real(A1[1,1])I) - @test_throws ArgumentError ldltfact(A1 + A1', shift=-2real(A1[1,1])) @test_throws ArgumentError cholfact(A1) @test_throws ArgumentError cholfact(A1, shift=1.0) @test_throws ArgumentError ldltfact(A1) @test_throws ArgumentError ldltfact(A1, shift=1.0) + @test !isposdef(cholfact(A1 + A1' - 2eigmax(Array(A1 + A1'))*I)) + @test !isposdef(cholfact(A1 + A1', shift=-2eigmax(Array(A1 + A1')))) + @test !LinAlg.issuccess(ldltfact(A1 + A1' - 2real(A1[1,1])*I)) + @test !LinAlg.issuccess(ldltfact(A1 + A1', shift=-2real(A1[1,1]))) F = cholfact(A1pd) tmp = IOBuffer() show(tmp, F) @@ -738,3 +738,14 @@ for F in (cholfact(AtA), cholfact(AtA, perm=1:5), ldltfact(AtA), ldltfact(AtA, p @test C == Ctest #Make sure C didn't change end end + +@testset "Issue #22335" begin + A = speye(3) + @test LinAlg.issuccess(cholfact(A)) + A[3, 3] = -1 + F = cholfact(A) + @test !LinAlg.issuccess(F) + @test LinAlg.issuccess(ldltfact!(F, A)) + @test A[:, 3:-1:1]\ones(3) == [-1, 1, 1] +end +