Isposdef and eigvals do not agree

Specific Domains Statistics
, ,
1

Hello, I was dealing with some matrices (for a statistical model), when I came across this case which raised an error during the fitting. I was trying to compute

beta_iter[t] = rand(MvNormal(b_star, inv(A_star)))

but a matrix A_star raised the error “LinearAlgebra.PosDefException(-1)”.

However, I don’t understand why, since all its eigenvalues are strictly positive, which would mean that the matrix, and also its inverse inside the multivariate normal, are positive definite. But instead the isposdef() command seems to disagree about the inverse.

A_star = [0.5695141296842557 -0.5958323106341032; -0.5958323106341032 5.607107634272117]
invA_star = inv(A_star)
# invA_star = [1.9755086272843259 0.20992496432168561; 0.20992496432168564 0.20065248429953192]

println(eigvals(A_star)) # [0.5000000000009021, 5.67662176395547] all >0
println(eigvals(invA_star)) # [0.1761611115874664, 1.9999999999963913] still all >0

println(isposdef(A_star)) # true
println(isposdef(invA_star)) # false

The erorr occured at the 26th iteration, meaning that all the before cases were correct, so for now I was thinking that it is a numerical problem, especially also because the definition of A_star is the sum of positive definite matrices so it should preserve the positive definitiveness:

A_star = I(p)/s2_beta
for j in 1:n
	X_jt = Xlk_covariates[j,:,t]
	A_star += (X_jt * X_jt') / sig2h_iter[Si_iter[j],t]
end

The packages involved should be just LinearAlgebra and Distributions.

1 Like
2

The problem is

julia> ishermitian(invA_star)
false

If you instead make A_star Symmetric, it will be faster and work as expected.

julia> A_star = Symmetric([0.5695141296842557 -0.5958323106341032; -0.5958323106341032 5.607107634272117])
2×2 Symmetric{Float64, Matrix{Float64}}:
  0.569514  -0.595832
 -0.595832   5.60711

julia> invA_star = inv(A_star)
2×2 Symmetric{Float64, Matrix{Float64}}:
 1.97551   0.209925
 0.209925  0.200652

julia> isposdef(invA_star)
true
3 Likes
3

Hmm, symmetrized inverse is recognised as positive definite:

julia> invA_star
2×2 Matrix{Float64}:
 1.97551   0.209925
 0.209925  0.200652

julia> invA_star_sym = (invA_star+invA_star')/2
2×2 Matrix{Float64}:
 1.97551   0.209925
 0.209925  0.200652

julia> norm(invA_star-invA_star_sym)
7.850462293418876e-17

julia> isposdef(invA_star_sym)
true

Could the conditioning of the matrix be the culprit? It is not that bad, is it?

julia> cond(A_star)
11.353243527890454

julia> cond(invA_star)
11.35324352789046

I am also eager to read some explanation here.

4

The problem is that a matrix must be hermitian (symmetric) to be positive definite (see Definite matrix - Wikipedia) and inv doesn’t preserve symmetry unless you tell it that the input is a Hermitian (Symmetric) matrix.

1 Like
5

I see. So the matrix must be perfectly symmetric for the isposdef function, not just numerically so-so symmetric:

julia> isapprox(invA_star,invA_star',rtol=10e-17)
true

julia> isapprox(invA_star,invA_star',rtol=10e-18)
false

Indeed, the source code for the isposdef reads

isposdef(A::AbstractMatrix) =
    ishermitian(A) && isposdef(cholesky(Hermitian(A); check = false))
6

First of all thank you everyone for the answers :slight_smile:

Then in the end I also ended up to this other question (that I didn’t found when composing mine): https://discourse.julialang.org/t/all-positive-eigen-values-but-false-in-isposdef/107393 so yes the problem seems to be solvable by forcing the matrix to be Symmetric or Hermitian.

1 Like
7

That isn’t accurate. The Wikipedia page says “A real symmetric matrix is positive definite if…”. Unsymmetric square matrices can be positive definite as well. There’s a section about that on the same page: Definite matrix - Wikipedia

For example,

\begin{bmatrix} 1 & 0 & -2 \\ 0 & 2 & 0 \\ 2 & 0 & 1 \end{bmatrix}

is unsymmetric positive definite.

8

Ah, the Julia isposdef isn’t aware (or actively doesn’t support, I’m not fully sure) this extension of the definition. isposdef is documented as

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

9

IsApprox.jl supports testing for almost postitive definite.

julia> using IsApprox

julia> isposdef(invA_star, Approx())
true

julia> isposdef(invA_star, Equal())
false
2 Likes
10

I think it would make sense for isposdef to require the input to be Hermitian. Cholesky won’t be well defined for unsymmetric matrices, even one that is positive definite according to the extended definition.