Inverse of `Hermitian` complex matrix giving wrong results

[marvinlenk]

I noticed a problem with Hermitian matrices, especially the inverse of them. For any given invertible square matrix x, x / x and x \ x should be the identity matrix (or at least close to it). Invoking this on a Hermitian complex matrix y = Hermitian(x) results in a filled complex matrix, that is not even hermitian. At least y / y is the hermitian conjugate of y \ y.

Using z = hermitianpart(y) does however give correct results for z / z and z \ z. Notably, z == y and typeof(z) == typeof(y) is true. The inverse also works for Hermitian real matrices.

Minimal working example:

using LinearAlgebra

x = rand(ComplexF64,2,2)
y = Hermitian(x)
z = hermitianpart(y)

y \ y ≈ I
y / y ≈ I

z \ z ≈ I
z / z ≈ I

Output is false for

y \ y ≈ I
y / y ≈ I

but expected to be true.

Output of versioninfo():

Julia Version 1.10.0
Commit 3120989f39b (2023-12-25 18:01 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 4 × Intel(R) Core(TM) i5-6500 CPU @ 3.20GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, skylake)
  Threads: 3 on 4 virtual cores
Environment:
  JULIA_NUM_THREADS = 3

[IanButterworth]

This looks to be related to this weird bug https://github.com/JuliaLang/LinearAlgebra.jl/issues/1064

[mcabbott]

This bug seems to exist in Julia 1.0:

julia> y = [0.279982+0.988074im  0.770011+0.870555im
            0.138001+0.889728im  0.177242+0.701413im] |> Hermitian;

julia> inv(y) * y
2×2 Array{Complex{Float64},2}:
  0.648315+5.55112e-17im  0.0487625+0.0551296im  
 0.0686912-0.0776606im     0.646396-5.55112e-17im

julia> inv(collect(y)) * y
2×2 Array{Complex{Float64},2}:
          1.0-1.11022e-16im  -1.20814e-17+0.0im
 -1.56742e-17+2.77556e-17im           1.0+0.0im

julia> VERSION
v"1.0.5"

I get the same numbers on 1.12.0-DEV.320.

[dkarrasch]

We forgot to realify the main diagonal of the underlying data storage. Fix is coming.