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)
3

In this case, H1 and H2 are symmetric and the matrices of eigenvectors have condition number very close to one, as expected. The eigenvalues are reasonably well separated, so the eigenspaces for each eigenvalue are also well conditioned. Using inv is a good thing to avoid, but it should be harmless here. The difference in sign is not surprising, since eigenvectors can certainly be scaled by \pm 1. I think the real problem is that scaling the eigenvectors by \pm 1 does not correspond to just a change of sign in Vs * a' * inv(Vs). If D is a matrix such that Vs*D has columns multiplied by \pm 1, then the matrix of interest becomes Vs * D * a' * inv(Vs * D), which does not just change signs. Scaling the rows of Vs by \pm 1 would just change signs in that matrix, but it doesn’t work for columns (eigenvectors).

To do a quick check, with a quick change to align the signs, I get the expected results for:

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)
Es1, V1s = eigen(H1 / GHz)
@show (Vs1 * a' * inv(Vs1))[[2,Ns+2], 1]'
Es2, Vs2 = eigen(H2 ./ GHz)
for j in 1:length(Es2)
    Vs2[:,j] = Vs2[:,j]*sign(Vs2[:,j]'*Vs1[:,j])
end
@show  (Vs2 * a' * inv(Vs2))[[2,Ns+2], 1]'

which gives

((Vs1 * a' * inv(Vs1))[[2, Ns + 2], 1])' = adjoint([0.9998247036167949, 0.01440923431692306])
((Vs2 * a' * inv(Vs2))[[2, Ns + 2], 1])' = adjoint([0.9998247036167945, 0.014409234316921217])

Of course, this is all a bit trickier if this is just a simple MWE and the real matrix has ill conditioned eigenspaces.

4 Likes
4

This looks like the right answer! These are indeed my matrices of interest, so I’m not using poorer-conditioned cases. I think I had implicitly assumed that my object of interest Vs * a' * inv(Vs) is “gauge invariant” and would not depend on the phases of the eigenvectors. The “correct” sign convention would be that diag(Vs) are all approximately +1 and not minuses, which I think in this case actually corresponds to fixing Vs1 to match Vs2 not vice versa like you have :slight_smile: