ishermitian test's A[i,j] != ctranspose(A[j,i]) which causes problems with numerical accuracy

[dhoegh]

I have a function where i determine the hessian with finite difference using the hessian function in the Calculus package. Then I test if the hessian is positive definite with isposdef this always return false because, ishermitian is called at:

julia/base/linalg/dense.jl

Lines 27 to 28 in 5e2a617

	isposdef!{T<:BlasFloat}(A::StridedMatrix{T}, UL::Symbol) = LAPACK.potrf!(char_uplo(UL), A)[2] == 0
	isposdef!(A::StridedMatrix) = ishermitian(A) && isposdef!(A, :U)

before the LAPACK.potrf! is called. Due to numerical inaccuracy will ishermitian always return false because A[i,j] != ctranspose(A[j,i]) in

julia/base/linalg/generic.jl

Line 282 in 5e2a617

if A[i,j] != ctranspose(A[j,i])

, will return true for some i,j due to numerical inaccuracy. If I change the line to !isapprox(A[i,j],ctranspose(Aj,i])) the ishermitian returns true and hence also call the LAPACK.potrf! to check if it is positive definite.

So from my point of view, should ishermitian test approximately equality.

[andreasnoack]

I think the better solution here is to define isposdef! for Hermitian and Symmetric types because you probably know that the matrix is theoretically Hermitian, right?

[dhoegh]

Yes I do, but then I would need to construct the special matrix myself? If it isn’t changed a note in the documentation about the problem would probably be a good idea.

[andreasnoack]

You would only have to write isposdef(Hermitian(A)). Wouldn't that be okay? It would probably also be much faster because you can avoid the check for Hermitianity.

[dhoegh]

When you know that is the trick, it is a good solution. On the other side a newcomer to the language could get really stumped by the fact isposdef(A) do not return true as it should and do not realize that he should do isposdef(Hermitian(A)), and I must admit that my linear algebra lessons where too far away to remember what Hermitian meant before I wikied it. I only found the problem because I explored the source and found that LAPACK.potrf! returned what I expected.

[jiahao]

The only reasonable thing I can think of doing here is to relax the exact inequality test to an approximate one with user-specifiable threshold.

[andreasnoack]

I don't think the two solutions exclude each other: having a floating point issym and ishermitian for StridedMatrix with approximate testing, but also isposdef! methods for Hermitian and Symmetric.

[dhoegh]

Sounds like a good idea to test for approximate equality for StridedMatrix and for performance isposdef(Hermitian(A)) if the user supply the info. Maybe isposdef for StridedMatrix should throw an argument error if the matrix is not approximately Hermitian to alert the user.

[stevengj]

Note that Matlab's ishermitian function seems to check for exact equality. @alanedelman, do you have an opinion here?

What would be the right criterion? norm(A - A') ≪ norm(A) in some cheap norm? But this might be deceptive if you care mainly about the smallest-magnitude eigenvalues.

[dhoegh]

If we can agree on something I would be happy to prepare a pull. As @stevengj pointed out Matlab do test for exact equallity for both ishermitian and issymmetric which might be a reason to keep doing this. However in order to prevent the isposdef to return a false negative to the user, which do not realize this is due to the matrix fails ishermitian due to numerical noise we could make isposdef raise an error if the matrix fails ishermitian. This error could be turnoff by having an extra argument raise=true which defaults to true.
The second part of the fix will be to implementing isposdef! for Hermitian and Symmetric types.

Do this sound like a viable path of solving this issue?

[andreasnoack]

I think we should use a tolerance for matrices with floating point elements even though MATLAB has chosen a different solution. I don't see why you'd be interested in exact symmetry when working with floating point values. I asked @alanedelman and he is fine with such a change. So basically there should be two versions of issym and ishermitian: one for floating point elements and one for everything else where the former uses approximate and the latter uses exact equality.

It would be great if you could prepare a pull request with something like this.

[dhoegh]

Should istril and istriu also be tested if they are approximate zero for floatingpoints? Should the approximate test for ishermitian be performed with isapprox with default rtol,atol and then they can be changed if the user supply them.

[stevengj]

@andreasnoack, even if you accept that these tests should be approximate, it still seems to be an open question what those tests should be. (The test in #10369 is clearly wrong.) As I mentioned above, possibly it should be something like somenorm(A - A') <= thresh*somenorm(A), but it's not obvious to me what norm should be used and what its impact would be.

[andreasnoack]

I'm pretty sure we won't do this at this point. One thing we could consider is to remove the Hermitian wrapper in

julia/base/linalg/dense.jl

Line 92 in 5132d3e

isposdef!(A::AbstractMatrix) = ishermitian(A) && isposdef(cholfact!(Hermitian(A)))

If we do that then isposdef will throw for non-Hermitian matrices instead of returning false. This will force users to explicitly wrap their matrices in Hermitian if they want to suppress the small floating point errors. The only problem with that is that #17367 is still open.

[fredrikekre]

It is weird if isposdef throws though, it should just answer the question with true or false. I imagine the ishermitian check is there for this particular reason.

[andreasnoack]

I can see that but it will mainly throw in the cases people have complained about it returning false, i.e. when the matrix is almost Hermitian. For most people, the right solution would be to have the matrix wrapped in Hermitian when passing it. In that case, it won't throw.

[fredrikekre]

Wouldn't we end up here then:

A = [2. 1.; 0. 1.]
isposdef(A)            # errors because A is not Hermitian, error message suggests using Hermitian(A)
isposdef(Hermitian(A)) # Success! But clearly A is not positive definite.

[fredrikekre]

Similarly, this error message is not very good either:

julia> A = [2. 1.; 0. 1.];

julia> cholfact(A)
ERROR: ArgumentError: matrix is not symmetric/Hermitian. This error can be avoided by calling cholfact(Hermitian(A)) which will ignore either the upper or lower triangle of the matrix.
Stacktrace:
 [1] cholfact(::Array{Float64,2}, ::Val{false}) at ./linalg/cholesky.jl:318 (repeats 2 times)

julia> cholfact(Hermitian(A))
Base.LinAlg.Cholesky{Float64,Array{Float64,2}} with factor:
[1.41421 0.707107; 0.0 0.707107]
successful: true

Hermitian(A) is clearly not the matrix I tried to factorize, instead I should obtain a factorization with successful: false. Now that we delay throwing for cholfact, we should probably remove that error, no? Currently we treat lufact and cholfact differently here.

[andreasnoack]

I think the first error message is pretty clear about what will happen if you wrap your matrix in Hermitian. The Hermitian wrapper enforces that the matrix is Hermitian so the Cholesky should only fail when the matrix, as it looks like through the Hermitian wrapper, is not PD so I think the current behavior is reasonable.

The point regarding the isposdef is that people mainly tests with isposdof to figure out the eigenvalues are positive, e.g. that to ensure that it is a valid covariance matrix. They already know/assume that the matrix is symmetric/Hermitian.

Currently we treat lufact and cholfact differently here.

Could you explain what you mean here. I'm not sure I follow.

[stevengj]

There is also a more general meaning of "positive-definite" that applies to arbitrary matrices: whether the Hermitian part (A+A')/2 is positive-definite, or equivalently whether eigvals(A) have positive real parts.

[KristofferC]

Could you explain what you mean here. I'm not sure I follow.

cholfact throws on non posdef input, lufact returns a sucess = false factorization on singular inputs and defers throwing to when the factorization object is actually used.

[fredrikekre]

Could you explain what you mean here. I'm not sure I follow.

cholfact throws on non posdef input, lufact returns a sucess = false factorization on singular inputs and defers throwing to when the factorization object is actually used.

Exactly, we now treat bad input differently; for cholfact we throw if the input is non-hermitian, but not if the input is bad in the sense of being non-positive definite. This changed with #21976, pre #21976 we threw for both cases of bad input to cholfact now we throw for one case but not the other. I think we should never throw, like lufact. Getting of topic here maybe...

[andreasnoack]

Okay, I get it. It makes sense now that we set the flag. This was also kind of a deprecation of the old behavior. Since it happened in 0.6 we could change to just setting the flag now.

[andreasnoack]

I don't there is anything left to do here.

[GiggleLiu]

I think the behavior of ishermitian is not consistent with other julia methods such as eigen.
As shown in #28885, This can be quite non-intuitive. Please at least consider using keyword parameteratol.