isposdef rewrite with help of cholfact

Author: fredrikekre

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]))

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,8 @@ 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{T}) where {T} =
+    isposdef!(copy_oftype(A, promote_type(typeof(chol(one(T))),Float32)))
 isposdef(x::Number) = imag(x)==0 && real(x) > 0
 
 stride1(x::Array) = 1

base/linalg/symmetric.jl

@@ -308,8 +308,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)

test/linalg/cholesky.jl

@@ -67,23 +67,33 @@ 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"
         else
             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)

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)

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
 

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]