Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

isposdef() is incorrect #20004

Closed
wilburtownsend opened this issue Jan 12, 2017 · 3 comments · Fixed by #22245
Closed

isposdef() is incorrect #20004

wilburtownsend opened this issue Jan 12, 2017 · 3 comments · Fixed by #22245
Labels
docs This change adds or pertains to documentation linear algebra Linear algebra

Comments

@wilburtownsend
Copy link

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)
@tkelman tkelman added the linear algebra Linear algebra label Jan 12, 2017
@andreasnoack
Copy link
Member

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.

@andreasnoack andreasnoack added the docs This change adds or pertains to documentation label Jan 13, 2017
@akaysh
Copy link
Contributor

akaysh commented Jan 13, 2017

@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
Copy link
Contributor

tkelman commented Jan 14, 2017

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
docs This change adds or pertains to documentation linear algebra Linear algebra
Projects
None yet
Development

Successfully merging a pull request may close this issue.

4 participants