Julia thinks inverse of a Hermitian matrix is not Hermitian

General Usage
1

Here is my code:

using LinearAlgebra

a = randn(50,50);
a = (a+a')/2;
println(ishermitian(a))
println(ishermitian(inv(a)))

The output is

True
False

However, the inverse of a Hermitian matrix is supposed to be Hermitian.

Turns out that the non-Hermitian part of inv(a) is small, so this is a machine precision error. The number maximum(abs.(inv(a)-inv(a)')) is \sim10^{-14}.

I am using inverse of a particular Hermitian matrix to generate a quantum Hamiltonian, and it needs to be Hermitian. As a check, I am verifying the Hamiltonian is hermitian before running the code.

However, due to the numerical error I mentioned above, there is no way to distinguish between non-hermiticity generated by machine precision errors, and non-hermiticity generated by typo/bugs in my code.

Is there a solution to this problem?

2

You can use a Hermitian wrapper. The function inv will return another instance of Hermitian.

2 Likes
3

Can you please provide a small example of how to use the wrapper?

4 
julia> using LinearAlgebra

julia> a = Hermitian(randn(3,3))
3×3 Hermitian{Float64, Matrix{Float64}}:
  0.155291  -1.07756   -0.185214
 -1.07756   -0.241846  -0.642284
 -0.185214  -0.642284  -0.152622

julia> inv(a)
3×3 Hermitian{Float64, Matrix{Float64}}:
  2.90744    0.352187  -5.01043
  0.352187   0.448985  -2.31687
 -5.01043   -2.31687    9.27839
5 Likes
5

note also that giving Julia this extra info will make the inverse faster since you don’t have to look at the whole matrix.

Also as per usual, you almost never should be using inv. factorization and linear solving is typically better.

4 Likes