Numerical sensitivity of eigvecs for `≈` identical (unitful) matrices?

General Usage
, , ,
1

I came across the following discrepancy, see MWE:

using LinearAlgebra, Unitful
GHz = 1u"GHz"; wq = 6.7565*GHz; wc = 20*GHz; g = 0.1847*GHz;
Ns = 10; op = diagm(1 => sqrt.(1:Ns)); b = kron(op, I(Ns+1)); a = kron(I(Ns+1), op); 
H1 = (wq*a'*a + wc*b'*b + g *  (b' .- b) * (a' .- a )  ) ; 
H2 = (wq*a'*a + wc*b'*b - g * (b'*a + a'*b - (a'*b' + a*b))) ;
H2 ≈ H1
# true
maximum(abs, H1 .- H2)
# 2.220446049250313e-16 GHz
# Disagreement starts here: up so signs, note the second component is off by ~0.001
Es, Vs = eigen(H1 / GHz); (Vs * a' * inv(Vs))[[2,Ns+2], 1]'
# 1×2 adjoint(::Vector{Float64}) with eltype Float64:
# -0.999825  -0.0144092
 Es, Vs = eigen(H2 ./ GHz); (Vs * a' * inv(Vs))[[2,Ns+2], 1]'
# 1×2 adjoint(::Vector{Float64}) with eltype Float64:
# 0.999838  0.0134753

The nonvanishing matrix elements are not that too small, so I’m surprised such a small numerical difference between H1 and H2 skews this result. The eigenvalues seem to agree and their order is not swapped:

maximum(eigvals(H1./ GHz) .- eigvals(H2./GHz) )
# 7.105427357601002e-13

Any ideas what’s causing this disagreement? Curiously, though before calling eigvecs I divide through by GHz, setting GHz = 1 initially instead of anything unitful resolves the disagreement (they coalesce to 0.0134…).

2

Forming an inverse explicitly is almost never a good idea and I think could cause loss of precision. Try \ instead:

Es, Vs = eigen(H1 / GHz); (Vs * a' \ Vs)