Diverging results of eigen() between Julia and Python

General Usage
4

Your iteration is ill-defined.

The i-th step of your algorithm corresponds to:

\text{Symmetric}(Q^{(i-1)} + oo^T, \text{:L}) = Q^{(i)} \Lambda^{(i)} Q^{(i)}

where o is the vector of all 1’s. That is, you add 1’s to previous eigenvectors Q^{(i-1)}, symmetrize by throwing away the upper half (since Q is not symmetric! this looks suspicious to me), then find the new eigenvectors Q^{(i)}.

The problem with this iteration is that the eigenvectors are not unique: even if they are real and orthonormal and the eigenvalues are distinct, you can still multiply any eigenvector (any column of Q) by \pm 1. So, for your iteration to be well-defined, it must be invariant under sign flips of the eigenvectors, but it isn’t! It is easy to check that simply multiplying your eigenvectors on the right by Diagonal([1,-1,1]) completely changes the eigenvalues and eigenvectors after one iteration.

For example, with the simplest Q matrix of all, Q=I, you get entirely different eigenvalues/eigenvectors of Q + oo^T if you flip the sign of the middle vector:

julia> eigen(Symmetric(Q + o*o', :L))
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
3-element Vector{Float64}:
 0.9999999999999992
 1.0
 4.0
vectors:
3×3 Matrix{Float64}:
  0.751201  -0.319943  -0.57735
 -0.652679  -0.490588  -0.57735
 -0.098522   0.810531  -0.57735

julia> eigen(Symmetric(Q*Diagonal([1,-1,1]) + o*o', :L))
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
3-element Vector{Float64}:
 -0.5615528128088254
  1.0000000000000004
  3.5615528128088303
vectors:
3×3 Matrix{Float64}:
  0.260956   0.707107     -0.657192
 -0.92941    1.38778e-15  -0.369048
  0.260956  -0.707107     -0.657192

Different eigensolvers, even different versions of LAPACK/BLAS, are entitled to make different choices of signs for the eigenvectors, so your iteration can produce entirely different results. This is not a matter of roundoff errors, it’s a fundamental problem that what you are doing is ill-posed.

I’m not sure where you got this algorithm from, but probably there is a mistake in either the implementation or the derivation.

10 Likes