# Isposdef and eigvals do not agree

**URL:** <https://discourse.julialang.org/t/isposdef-and-eigvals-do-not-agree/118191>\
**Category:** Statistics\
**Tags:** linearalgebra, distributions, potential-bug\
**Created:** [August 14, 2024, 5:02pm UTC](https://discourse.julialang.org/t/isposdef-and-eigvals-do-not-agree/118191 "2024-08-14T17:02:19Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![fmor](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fmor/32/211337_2.png) [@fmor](https://discourse.julialang.org/u/fmor)\
**Post date:** [August 14, 2024, 5:02pm UTC](https://discourse.julialang.org/t/isposdef-and-eigvals-do-not-agree/118191/1 "2024-08-14T17:02:19Z")

</div>

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

```julia
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.

```julia
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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [August 14, 2024, 5:32pm UTC](https://discourse.julialang.org/t/isposdef-and-eigvals-do-not-agree/118191/2 "2024-08-14T17:32:15Z")

</div>

The problem is

```julia
julia> ishermitian(invA_star)
false

```

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

```julia
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

```

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [August 14, 2024, 5:34pm UTC](https://discourse.julialang.org/t/isposdef-and-eigvals-do-not-agree/118191/3 "2024-08-14T17:34:12Z")

</div>

Hmm, symmetrized inverse is recognised as positive definite:

```julia
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
julia> cond(A_star)
11.353243527890454

julia> cond(invA_star)
11.35324352789046

```

I am also eager to read some explanation here.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [August 14, 2024, 5:37pm UTC](https://discourse.julialang.org/t/isposdef-and-eigvals-do-not-agree/118191/4 "2024-08-14T17:37:51Z")

</div>

The problem is that a matrix must be hermitian (symmetric) to be positive definite (see [Definite matrix - Wikipedia](https://en.wikipedia.org/wiki/Definite_matrix)) and `inv` doesn’t preserve symmetry unless you tell it that the input is a `Hermitian` (`Symmetric`) matrix.

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [August 14, 2024, 5:59pm UTC](https://discourse.julialang.org/t/isposdef-and-eigvals-do-not-agree/118191/5 "2024-08-14T17:59:25Z")

</div>

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

```julia
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

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

```

---

<div class="post-metadata">

**Author:** ![fmor](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fmor/32/211337_2.png) [@fmor](https://discourse.julialang.org/u/fmor)\
**Post date:** [August 14, 2024, 6:55pm UTC](https://discourse.julialang.org/t/isposdef-and-eigvals-do-not-agree/118191/6 "2024-08-14T18:55:52Z")

</div>

First of all thank you everyone for the answers 🙂

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](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.

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [August 15, 2024, 4:55pm UTC](https://discourse.julialang.org/t/isposdef-and-eigvals-do-not-agree/118191/7 "2024-08-15T16:55:46Z")

</div>

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](https://en.m.wikipedia.org/wiki/Definite_matrix#Extension_for_non-Hermitian_square_matrices)

For example,

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

is unsymmetric positive definite.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [August 15, 2024, 5:01pm UTC](https://discourse.julialang.org/t/isposdef-and-eigvals-do-not-agree/118191/8 "2024-08-15T17:01:07Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![jlapeyre](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlapeyre/32/4514_2.png) [@jlapeyre](https://discourse.julialang.org/u/jlapeyre)\
**Post date:** [August 15, 2024, 10:02pm UTC](https://discourse.julialang.org/t/isposdef-and-eigvals-do-not-agree/118191/9 "2024-08-15T22:02:41Z")

</div>

[IsApprox.jl](https://github.com/jlapeyre/IsApprox.jl) supports testing for almost postitive definite.

```julia
julia> using IsApprox

julia> isposdef(invA_star, Approx())
true

julia> isposdef(invA_star, Equal())
false

```

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [August 16, 2024, 2:27pm UTC](https://discourse.julialang.org/t/isposdef-and-eigvals-do-not-agree/118191/10 "2024-08-16T14:27:33Z")

</div>

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.