`isposdef(E::Eigen)` is incorrect

[sztal]

Hi,

The problem is as in the title. In short, isposdef method for Eigen factorization is incorrect. It does not check for hermiticity/symmetricity and only looks at the eigenvalues, which is of course not enough.

using LinearAlgebra
T = UpperTriangular(rand(3, 3) .+ 1)
isposdef(T)  # False
isposdef(eigen(T))  # True

Best!

Version info

Julia Version 1.11.3
Commit d63adeda50d (2025-01-21 19:42 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 16 × AMD Ryzen 9 5900HX with Radeon Graphics
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, znver3)
Threads: 1 default, 0 interactive, 1 GC (on 16 virtual cores)
Environment:
  JULIA_PROJECT = @.
  JULIA_PKG_PRESERVE_TIERED_INSTALLED = true
  JULIA_REVISE = manual

[dkarrasch]

I suspect this is a matter of documentation, like saying that isposdef(E::Eigen) only checks for real and positive eigenvalues, assuming that E is the eigendecomposition of a self-adjoint matrix ("RealHermOrSymComplexHerm"). I don't know how one would check for symmetry: any such check would be subject to floating point inaccuracies, so likely V'V is not the identity, and/or V'Diagonal(evals)*V is not symmetric/Hermitian. I vaguely remember, however, that there was some effort to store that piece of information (whether E is the eigendecomposition of a self-adjoint matrix).

A quick search revealed that @RalphAS, @antoine-levitt @KlausC shared interest in this.

[aravindh-krishnamoorthy]

The problem is as in the title. In short, isposdef method for Eigen factorization is incorrect. It does not check for hermiticity/symmetricity and only looks at the eigenvalues, which is of course not enough.

Could you explain why this is not enough?

[stevengj]

It depends on what definition of positive-definiteness you are using.

• The "basic" definition is that a positive-definite matrix is a Hermitian matrix with positive eigenvalues. (There are also several equivalent properties.)
• A strict generalization of this is that PD means any matrix where real(dot(x, A * x)) > 0 for all $x \ne 0$, or equivalently A + A' (which is Hermitian) has positive eigenvalues, or equivalently which implies that the eigenvalues of $A$ have positive real parts but the converse is not generally true (see below).

Up to now, isposdef has used the former, more restrictive definition — it is automatically false for any non-Hermitian matrix.

[aravindh-krishnamoorthy]

In such a case, I feel that this is more than just an issue of documentation. This is because the definition of positive definiteness changes between isposdef(A) and isposdef(eigen(A)). I'm also not sure if relying on the type RealHermOrSymComplexHerm is sufficient for the case where a symmetric A is not wrapped in Symmetric. Furthermore, looking at the definition below:

LinearAlgebra.jl/src/eigen.jl

Line 126 in 3f46f5f

isposdef(A::Union{Eigen,GeneralizedEigen}) = isreal(A.values) && all(x -> x > 0, A.values)

this issue is applicable to both real- and complex-valued matrices.

[antoine-levitt]

or equivalently the eigenvalues of A have positive real parts

You mean A+A', right?

Yes, the current behavior is an unfortunate bug, but it's not very clear what to do, except storing symmetry in Eigen, documenting the bug or just deleting the ispodef(::Eigen) method. Storing the symmetry in Eigen is probably a good idea anyway...

[antoine-levitt]

Turns out it's such a good idea that @KlausC and me did it two years ago: JuliaLang/julia#48549 JuliaLang/julia#48657. I don't remember the outcome of this but it looks like one or both were ready to merge, in which case we can just merge it and change isposdef to use this? @KlausC ?

[sztal]

or equivalently the eigenvalues of A have positive real parts

You mean A+A', right?

Yes, the current behavior is an unfortunate bug, but it's not very clear what to do, except storing symmetry in Eigen, documenting the bug or just deleting the ispodef(::Eigen) method. Storing the symmetry in Eigen is probably a good idea anyway...

Yes, I agree with that. Given the practical importance of the distinction between hermitian/non-hermitian eigendecompositions (e.g. one can be mappend to SVD while the other can't; the inverse eigenvectors are obtained for free in the hermitian case and require extra calculations in the non-hermitian case etc. etc.) it makes a lot of sense to keep this information in Eigen objects.

[aravindh-krishnamoorthy]

Thank you.

I'd actually be in favour of deprecating and eventually removing isposdef(E::Eigen).

The reasoning is that (a) given an eigen decomposition, it's rather easy for a user to check positive definiteness via ishermitian(A) && isreal(E.values) && all(x -> x > 0, E.values) and the check does not need any "internal" knowledge; (b) symmetry information doesn't really seem to be a good addition to Eigen and will not be useful for other purposes; (c) the Cholesky-decomposition-based check in isposdef(A) is not that bad in terms of computational complexity so the gains from trying to optimise it out (when an eigen decomposition is available) are not that high; (d) as @stevengj mentioned in his comment (elsewhere; and I agree), in most practical cases, positive definiteness of a matrix is usually known, e.g., A = B*B'.

[antoine-levitt]

symmetry information doesn't really seem to be a good addition to Eigen and will not be useful for other purposes

Why do you say this? It's not much bother adding it to the Eigen structure, and it enables a whole bunch of methods. Look at JuliaLang/julia#48549

[aravindh-krishnamoorthy]

Why do you say this? It's not much bother adding it to the Eigen structure, and it enables a whole bunch of methods. Look at JuliaLang/julia#48549

Thank you for the comment. I just think it's a slippery slope, as shown in the trivial example below:

julia> A = [1 1im; conj(eps()+1im) 1];
julia> B = [1 1im; conj(1im) 1];
julia> ishermitian(A)
false
julia> ishermitian(B)
true
julia> EA = eigen(A);
julia> EB = eigen(B);
julia> EA.vectors*EA.vectors' ≈ EB.vectors*EB.vectors'
true
julia> inv(EA.vectors) ≈ EA.vectors'
true
julia> inv(EB.vectors) ≈ EB.vectors'
true

On the other hand, maintaining a RealHermOrSymComplexHermEigen might work for matrices wrapped in Symmetric or Hermitian. But then again, a user is actively wrapping this so the added value is low, imho.

[antoine-levitt]

OK, point taken, my proposed change (storing the unitarity information in Eigen) would not solve this because it's not easy to know what to do when the eigenvectors are not unitary. My preference would be to error on isposdef(E::Eigen) when E.unitary_eigenvectors is false (ie, Eigen would only support those methods that make sense and are fast). That's breaking, but it could be done eventually, and it would not break the most common case (doing this on a symmetric/hermitian matrix). The alternative is the current behavior (which is just wrong) or perfoming an expensive cholesky.

[stevengj]

You mean A+A', right?

Right, sorry — if A+A' is positive definite (has positive eigenvalues), then the eigenvalues of A have positive real parts, but the converse is not necessarily true.

For example, A = [1 3; 0 2] is a 2x2 non-symmetric matrix with positive real eigenvalues 1 and 2, but A + A' is not positive definite (it has a negative eigenvalue). That is:

julia> A = [1 3; 0 2]
2×2 Matrix{Int64}:
 1  3
 0  2

julia> isposdef(A+A')
false

julia> isposdef(eigen(A))
true

which is clearly an incorrect result, even with the "generalized" definition of positive-definiteness for non-Hermitian matrices.