Julia: isposdef() is incorrect

[wilburtownsend]

isposdef() is claiming that some (though not all) positive-semi definite matrices are positive definite when they are not:

julia> A = [1.69  2.21; 2.21  2.89]
2×2 Array{Float64,2}:
 1.69  2.21
 2.21  2.89
julia> det(A)
0.0
julia> isposdef(A)
true
julia> isposdef(ones(2,2))
false

I'm not sure what's causing this behaviour -- it looks like isposdef() attempts to make a Cholesky decomposition of the matrix, and all positive semi-definite matrices have a Cholesky decomposition (though IIRC there is no general algorithm for finding them). It seems to me that it would be simpler to just implement Sylvestor's criterion with something like

sylvestor(A) = issymmetric(A) & all([det(A[1:j, 1:j]) > 0 for j = 1:size(A)[1]])

Version info:

julia> versioninfo()
Julia Version 0.5.0
Commit 3c9d753 (2016-09-19 18:14 UTC)
Platform Info:
  System: NT (x86_64-w64-mingw32)
  CPU: Intel(R) Xeon(R) CPU E3-1240 v3 @ 3.40GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.7.1 (ORCJIT, haswell)

[andreasnoack]

The rounding during the computation of the Cholesky factorization can cause this. Generally, it is not possible to reliably test singularity of floating point matrices. The matrix you provided is actually not singular in its binary floating point representation

julia> det(big(A))
2.13162820728030047872728228063319445882714347172137703267935663215914838016073e-16

Determinants are generally quite unreliable in numerical computations because the terms grow so fast so Sylverster's criterion is not really practical. In numerical computations, the test via Cholesky is in most cases what you want and it is fairly cheap so I don't think we'll change it. However, it might be a good idea to extend the documentation to explain the method and maybe also use the example here to show the limitations. I'll add the documentation label to the issue.

[akaysh]

@andreasnoack Can this be the summary of the situation here. I will add it to the docs along with the example after approval.
isposdef(A) → Bool
Test whether a matrix is positive definite using Cholesky decomposition as Sylvester's criterion becomes impractical because of unreliable determinant computations of large numbers. In case of singular (non-invertible) positive semi-definite floating point matrices the function might not work as desired because of rounding involved in Cholesky factorisation and the non singular binary floating point representation of the matrix.

[tkelman]

I don't think it's necessary to mention Sylvester's criterion.