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 not defined for symmetric matrices #422

Closed
mcopik opened this issue Apr 26, 2017 · 5 comments · Fixed by JuliaLang/julia#21561
Closed

isposdef not defined for symmetric matrices #422

mcopik opened this issue Apr 26, 2017 · 5 comments · Fixed by JuliaLang/julia#21561

Comments

@mcopik
Copy link

mcopik commented Apr 26, 2017

The standard linalg function isposdef does not work for symmetric matrices (recent git master):

x = rand(5, 5)
x = x*x'
isposdef(x)
isposdef(Symmetric(x))

The last line fails with an error:

ERROR: MethodError: no method matching isposdef!(::Symmetric{Float64,Array{Float64,2}})

which I guess is caused by a missing definition. One workaround is to access underlying matrix directly but it is inconvenient:

isposdef(Symmetric(x).data)

I think this behavior is not correct and surprising for each user, I have not found any mention of such limitation in the documentation. There is an additional implementation for Diagonal matrices and this works fine:

isposdef(Diagonal(rand(1, 5)))
@fredrikekre
Copy link
Member

The same applies for Hermitian. Something like this seems to be sufficient (for dense matrices at least):

isposdef!(A::HermOrSym{<:BlasFloat, <:StridedMatrix}) = ishermitian(A) && LAPACK.potrf!(A.uplo, A.data)[2] == 0

mirroring the definitions here (and here). I will submit a PR.

@mcopik
Copy link
Author

mcopik commented Apr 26, 2017

@fredrikekre Why is ishermitian(A) necessary here? You're just passing a Hermitian/Symmetric matrix, right?

@fredrikekre
Copy link
Member

For cases like this:

julia> A = complex.(rand(2, 2), rand(2, 2))
2×2 Array{Complex{Float64},2}:
 0.127036+0.875987im  0.702463+0.796754im
 0.631628+0.488688im  0.214832+0.825713im

julia> ishermitian(Symmetric(A))
false

@mcopik
Copy link
Author

mcopik commented Apr 26, 2017

@fredrikekre I see, thanks. I forgot that one can make a complex symmetric matrix.

@andreasnoack
Copy link
Member

I think the post 0.6 solution should be to store info in the Cholesky type and delay throwing. That is what we do for LU and BunchKaufman and then we could just have isposdef(F::Cholesky) = F.info == 0 and something like isposdef(A::AbstractMatrix) = ishermitian(A) && isposdef(cholfact(A)).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

4 participants