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

**URL:** <https://discourse.julialang.org/t/numerical-sensitivity-of-eigvecs-for-identical-unitful-matrices/134846>\
**Category:** General Usage\
**Tags:** linearalgebra, unitful, eigenvectors, accuracy\
**Created:** [January 2, 2026, 7:27pm UTC](https://discourse.julialang.org/t/numerical-sensitivity-of-eigvecs-for-identical-unitful-matrices/134846 "2026-01-02T19:27:42Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![bremez](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bremez/32/38777_2.png) [@bremez](https://discourse.julialang.org/u/bremez)\
**Post date:** [January 2, 2026, 7:27pm UTC](https://discourse.julialang.org/t/numerical-sensitivity-of-eigvecs-for-identical-unitful-matrices/134846/1 "2026-01-02T19:27:42Z")

</div>

I came across the following discrepancy, see MWE:

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

```julia-auto
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…).

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [January 2, 2026, 8:00pm UTC](https://discourse.julialang.org/t/numerical-sensitivity-of-eigvecs-for-identical-unitful-matrices/134846/2 "2026-01-02T20:00:12Z")

</div>

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

```julia-auto
Es, Vs = eigen(H1 / GHz); (Vs * a' \ Vs)

```

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [January 3, 2026, 1:00am UTC](https://discourse.julialang.org/t/numerical-sensitivity-of-eigvecs-for-identical-unitful-matrices/134846/3 "2026-01-03T01:00:04Z")

</div>

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:

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

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

---

<div class="post-metadata">

**Author:** ![bremez](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bremez/32/38777_2.png) [@bremez](https://discourse.julialang.org/u/bremez)\
**Post date:** [January 3, 2026, 2:11am UTC](https://discourse.julialang.org/t/numerical-sensitivity-of-eigvecs-for-identical-unitful-matrices/134846/4 "2026-01-03T02:11:09Z")

</div>

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 🙂