From 6a3364125346400d8d0e5a43c1d189a501ae7bab Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Wed, 7 Jun 2017 23:38:13 +0200 Subject: [PATCH] isposdef rewrite with help of cholfact --- NEWS.md | 4 ++++ base/deprecated.jl | 4 ++++ base/linalg/cholesky.jl | 10 ++++++---- base/linalg/dense.jl | 22 ++++++++-------------- base/linalg/symmetric.jl | 2 -- test/linalg/cholesky.jl | 11 +++++++++++ test/linalg/dense.jl | 7 ++++++- test/linalg/diagonal.jl | 1 + test/linalg/eigen.jl | 3 ++- 9 files changed, 42 insertions(+), 22 deletions(-) diff --git a/NEWS.md b/NEWS.md index 0f62292db74c1..bda76f483ccd5 100644 --- a/NEWS.md +++ b/NEWS.md @@ -77,6 +77,10 @@ Deprecated or removed * The `cholfact`/`cholfact!` methods that accepted an `uplo` symbol have been deprecated in favor of using `Hermitian` (or `Symmetric`) views ([#22187], [#22188]). + * `isposdef(A::AbstractMatrix, UL::Symbol)` and `isposdef!(A::AbstractMatrix, UL::Symbol)` + have been deprecated in favor of `isposdef(Hermitian(A, UL))` and `isposdef!(Hermitian(A, UL))` + respectively ([#22245]). + * The function `current_module` is deprecated and replaced with `@__MODULE__` ([#22064]). This caused the deprecation of some reflection methods (such as `macroexpand` and `isconst`), which now require a module argument. diff --git a/base/deprecated.jl b/base/deprecated.jl index 555998ca6f546..4defa517f7c34 100644 --- a/base/deprecated.jl +++ b/base/deprecated.jl @@ -1361,6 +1361,10 @@ end @deprecate cholfact!(A::StridedMatrix, uplo::Symbol, ::Type{Val{true}}; tol = 0.0) cholfact!(Hermitian(A, uplo), Val{true}, tol = tol) @deprecate cholfact(A::StridedMatrix, uplo::Symbol, ::Type{Val{true}}; tol = 0.0) cholfact(Hermitian(A, uplo), Val{true}, tol = tol) +# PR #22245 +@deprecate isposdef(A::AbstractMatrix, UL::Symbol) isposdef(Hermitian(A, UL)) +@deprecate isposdef!(A::StridedMatrix, UL::Symbol) isposdef!(Hermitian(A, UL)) + # also remove all support machinery in src for current_module when removing this deprecation # and make Base.include an error _current_module() = ccall(:jl_get_current_module, Ref{Module}, ()) diff --git a/base/linalg/cholesky.jl b/base/linalg/cholesky.jl index ca9f7c2be51c3..434a7fdbe5a5c 100644 --- a/base/linalg/cholesky.jl +++ b/base/linalg/cholesky.jl @@ -288,7 +288,7 @@ and return a `Cholesky` factorization. The matrix `A` can either be a [`Symmetri `StridedMatrix` or a *perfectly* symmetric or Hermitian `StridedMatrix`. The triangular Cholesky factor can be obtained from the factorization `F` with: `F[:L]` and `F[:U]`. The following functions are available for `Cholesky` objects: [`size`](@ref), [`\\`](@ref), -[`inv`](@ref), and [`det`](@ref). +[`inv`](@ref), [`det`](@ref), [`logdet`](@ref) and [`isposdef`](@ref). # Example @@ -461,17 +461,19 @@ function A_ldiv_B!(C::CholeskyPivoted, B::StridedMatrix) end end +isposdef(C::Cholesky) = C.info == 0 + function det(C::Cholesky) - C.info == 0 || 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 - dd + @assertposdef dd C.info end function logdet(C::Cholesky) - C.info == 0 || throw(PosDefException(C.info)) + # need to check first, or log will throw DomainError + isposdef(C) || throw(PosDefException(C.info)) dd = zero(real(eltype(C))) @inbounds for i in 1:size(C.factors,1) dd += log(real(C.factors[i,i])) diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl index c1dc43965e13b..c58d9afe6fe10 100644 --- a/base/linalg/dense.jl +++ b/base/linalg/dense.jl @@ -29,13 +29,12 @@ function scale!(X::Array{T}, s::Real) where T<:BlasComplex X end -# Test whether a matrix is positive-definite -isposdef!(A::StridedMatrix{<:BlasFloat}, UL::Symbol) = LAPACK.potrf!(char_uplo(UL), A)[2] == 0 - """ isposdef!(A) -> Bool -Test whether a matrix is positive definite, overwriting `A` in the process. +Test whether a matrix is positive definite by trying to perform a +Cholesky factorization of `A`, overwriting `A` in the process. +See also [`isposdef`](@ref). # Example @@ -51,16 +50,14 @@ julia> A 2.0 6.78233 ``` """ -isposdef!(A::StridedMatrix) = ishermitian(A) && isposdef!(A, :U) +isposdef!(A::AbstractMatrix) = ishermitian(A) && isposdef(cholfact!(Hermitian(A))) -function isposdef(A::AbstractMatrix{T}, UL::Symbol) where T - S = typeof(sqrt(one(T))) - isposdef!(S == T ? copy(A) : convert(AbstractMatrix{S}, A), UL) -end """ isposdef(A) -> Bool -Test whether a matrix is positive definite. +Test whether a matrix is positive definite by trying to perform a +Cholesky factorization of `A`. +See also [`isposdef!`](@ref) # Example @@ -74,10 +71,7 @@ julia> isposdef(A) true ``` """ -function isposdef(A::AbstractMatrix{T}) where T - S = typeof(sqrt(one(T))) - isposdef!(S == T ? copy(A) : convert(AbstractMatrix{S}, A)) -end +isposdef(A::AbstractMatrix) = ishermitian(A) && isposdef(cholfact(Hermitian(A))) isposdef(x::Number) = imag(x)==0 && real(x) > 0 stride1(x::Array) = 1 diff --git a/base/linalg/symmetric.jl b/base/linalg/symmetric.jl index a3c2aa39c435e..4538617411db3 100644 --- a/base/linalg/symmetric.jl +++ b/base/linalg/symmetric.jl @@ -314,8 +314,6 @@ det(A::Symmetric) = det(bkfact(A)) inv(A::Hermitian{T,S}) where {T<:BlasFloat,S<:StridedMatrix} = Hermitian{T,S}(inv(bkfact(A)), A.uplo) inv(A::Symmetric{T,S}) where {T<:BlasFloat,S<:StridedMatrix} = Symmetric{T,S}(inv(bkfact(A)), A.uplo) -isposdef!(A::HermOrSym{<:BlasFloat,<:StridedMatrix}) = ishermitian(A) && LAPACK.potrf!(A.uplo, A.data)[2] == 0 - eigfact!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}) = Eigen(LAPACK.syevr!('V', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)...) function eigfact(A::RealHermSymComplexHerm) diff --git a/test/linalg/cholesky.jl b/test/linalg/cholesky.jl index 531f19ee15230..3152ca219475d 100644 --- a/test/linalg/cholesky.jl +++ b/test/linalg/cholesky.jl @@ -67,10 +67,14 @@ using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted, PosDefException capds = cholfact(apds) @test inv(capds)*apds ≈ eye(n) @test abs((det(capds) - det(apd))/det(capds)) <= ε*κ*n + @test logdet(capds) ≈ log(det(capds)) + @test isposdef(capds) if eltya <: BlasReal capds = cholfact!(copy(apds)) @test inv(capds)*apds ≈ eye(n) @test abs((det(capds) - det(apd))/det(capds)) <= ε*κ*n + @test logdet(capds) ≈ log(det(capds)) + @test isposdef(capds) end ulstring = sprint(show,capds[:UL]) @test sprint(show,capds) == "$(typeof(capds)) with factor:\n$ulstring" @@ -78,12 +82,18 @@ using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted, PosDefException capdh = cholfact(apdh) @test inv(capdh)*apdh ≈ eye(n) @test abs((det(capdh) - det(apd))/det(capdh)) <= ε*κ*n + @test logdet(capdh) ≈ log(det(capdh)) + @test isposdef(capdh) capdh = cholfact!(copy(apdh)) @test inv(capdh)*apdh ≈ eye(n) @test abs((det(capdh) - det(apd))/det(capdh)) <= ε*κ*n + @test logdet(capdh) ≈ log(det(capdh)) + @test isposdef(capdh) capdh = cholfact!(copy(apd)) @test inv(capdh)*apdh ≈ eye(n) @test abs((det(capdh) - det(apd))/det(capdh)) <= ε*κ*n + @test logdet(capdh) ≈ log(det(capdh)) + @test isposdef(capdh) ulstring = sprint(show,capdh[:UL]) @test sprint(show,capdh) == "$(typeof(capdh)) with factor:\n$ulstring" end @@ -270,6 +280,7 @@ end for T in (Float32, Float64, Complex64, Complex128) A = T[1 2; 2 1]; B = T[1, 1] C = cholfact(A) + @test !isposdef(C) @test_throws PosDefException C\B @test_throws PosDefException det(C) @test_throws PosDefException logdet(C) diff --git a/test/linalg/dense.jl b/test/linalg/dense.jl index 3f1b316ba3905..74d26f0601c93 100644 --- a/test/linalg/dense.jl +++ b/test/linalg/dense.jl @@ -49,7 +49,12 @@ bimg = randn(n,2)/2 apd = ainit'*ainit # symmetric positive-definite @testset "Positive definiteness" begin - @test isposdef(apd,:U) + @test !isposdef(ainit) + @test isposdef(apd) + if eltya != Int # cannot perform cholfact! for Matrix{Int} + @test !isposdef!(copy(ainit)) + @test isposdef!(copy(apd)) + end end @testset "For b containing $eltyb" for eltyb in (Float32, Float64, Complex64, Complex128, Int) binit = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex.(breal, bimg) : breal) diff --git a/test/linalg/diagonal.jl b/test/linalg/diagonal.jl index 6c2f8b2977bcb..d551b24548db0 100644 --- a/test/linalg/diagonal.jl +++ b/test/linalg/diagonal.jl @@ -235,6 +235,7 @@ end end @testset "isposdef" begin + @test isposdef(Diagonal(1.0 + rand(n))) @test !isposdef(Diagonal(-1.0 * rand(n))) end diff --git a/test/linalg/eigen.jl b/test/linalg/eigen.jl index c44563e5676f2..5e6ec48db16f2 100644 --- a/test/linalg/eigen.jl +++ b/test/linalg/eigen.jl @@ -35,14 +35,15 @@ aimg = randn(n,n)/2 @test eab[1] == eigvals(fill(α,1,1),fill(β,1,1)) @test eab[2] == eigvecs(fill(α,1,1),fill(β,1,1)) - d,v = eig(a) @testset "non-symmetric eigen decomposition" begin + d, v = eig(a) for i in 1:size(a,2) @test a*v[:,i] ≈ d[i]*v[:,i] end f = eigfact(a) @test det(a) ≈ det(f) @test inv(a) ≈ inv(f) + @test isposdef(a) == isposdef(f) @test eigvals(f) === f[:values] @test eigvecs(f) === f[:vectors]