Make isposdef return true for unsymmetric positive definite matrices

[mohamed82008]

Since this question https://discourse.julialang.org/t/isposdef-inconsistent-results-numerically-unstable/9013 got no further response, I am opening an issue hoping it won't be shot down.

[andreasnoack]

There are two definitions: a stricter which is simpler conceptually as well as easier to implement and a slightly more general definition which is less simple conceptually. Right now we use the former. What I miss is the good motivation for using the latter. If people frequently would like to test for the generalized version for non-symmetric matrices then it would make sense to consider using that version. However, in practice, the request has only come up in cases where the matrix was actually supposed to be symmetric but rounding had made it slightly non-symmetric. I don't think that is a compelling motivation. In these cases, if we changed isposdef to the more general definition, it would return true but the next thing many users would try would be cholfact which would fail because of the matrix being non-symmetric. In most cases, and specifically in your example from Discourse, the matrix is supposed to be symmetric and people would be better off using Symmetric.

[antoine-levitt]

The main confusion seems to be that isposdef returns false on non-symmetric positive definite matrices. Maybe a bit radical, but how about deprecating isposdef for general matrices (or making it be an error) and making it only work for Symmetric?

[KristofferC]

I don't think that is needed. The definition used in Julia for positive definiteness does not include the extension as described in https://en.wikipedia.org/wiki/Positive-definite_matrix#Extension_for_non-symmetric_matrices. Just add a note to the docs about it and done imo.

[mohamed82008]

Why not add an option to switch to the less restricted definition at least? If the main concern is confusing people used to or expecting the restricted definition and therefore calling cholfact on any matrix that returns true when passed to isposdef, then a default value will do the trick. Strictly speaking though, cholfact will still complain that the matrix is not symmetric which will probably trigger the user to call help on isposdef which can have the docstring explaining that the function is not checking for symmetry, so what was confusing at first will be just an educational experience. And performance-wise, I think finding the symmetric part and checking for symmetry are not too far away cost-wise, same complexity, so I think this proposal is a step towards making more people happy, at no penalty to the rest in the most modest proposal, i.e. adding an option with the current behavior as the default.

[StefanKarpinski]

If someone (perhaps you, @mohamed82008?) wants to implement the the more general definition and experiment with it, that seems like a good idea for a third-party package.

[mohamed82008]

That's possible of course, but it's unlikely that someone who wants to check for positive definiteness will want to look further than the standard library for linear algebra. If you are suggesting this because only few people may need this functionality then I understand, but I don't think it's the best way to do it.

[StefanKarpinski]

I'm suggesting that in an open source project, if someone wants to explore an idea like this, the most effective way to make that happen is to try it out themselves, demonstrate that it is useful and coherent and then use that working demonstration to convince other people that they are correct.

[mohamed82008]

If you mean to convince people that the function is theoretically correct, that doesn't need an implementation. But if you mean to make a correct implementation first to make it easier to merge later, if that's not off the table, then sure I can do that.

[Jutho]

As stated here:
https://discourse.julialang.org/t/isposdef-inconsistent-results-numerically-unstable/9013/15

A non-symmetric real matrix A cannot satisfy the definition v'*A*v > 0 if you also include complex vectors v. Symmetry (or hermiticity in case of complex A) follows directly from v'*A*v > 0 if you include complex vectors from the start.

[mohamed82008]

Note my response in the same discourse thread.

[KlausC]

https://discourse.julialang.org/t/isposdef-inconsistent-results-numerically-unstable/9013/32?u=klacru

[mohamed82008]

Can we at least stick to Householder's definition and throw an error when the matrix is not Hermitian? Quoting Householder: “A Hermitian matrix A is said to be positive definite in case x≠0→x^HAx>0; it is nonnegative semidefinite in case x^HAx≥0 for every x. For non-Hermitian matrices the concept is not defined.” I think for real matrices the fact that some authors allow unsymmetric matrices can lead to isposdef giving a silent wrong/unpredictable answer. I would rather get an error.

[KlausC]

I could imagine to make a PR with the following changes for isposdef:
isposdef(::Union{Real,Hermitian,Symmetric{<:Real}}) calculate a sensible value - no checks
isposdef(::Union{Number,Symmetric}) check by isreal and throw error if not - otherwise calculate
isposdef(::AbstractMatrix) check by ishermitian and throw error if not - otherwise calculate

But I am not sure, if @andreasnoack and @KristofferC could be convinced to approve this.

[UweAlex]

The documantation of Julia should mention that isposdef only works on a symetrical matrix

[KlausC]

It does!

help?> isposdef
search: isposdef isposdef!

  isposdef(A) -> Bool

  Test whether a matrix is positive definite (and Hermitian) by trying to perform a Cholesky
  factorization of A. See also isposdef!

For real matrices, "symmetric" is equivalent to "Hermitian" and for complex matrices, "symmetric" is not relevant for positive definiteness.