Diverging results of eigen() between Julia and Python

General Usage
1

I am converting an algorithm from python to julia that uses eigen decomposition in the “hot loop” updating the solution. Python version uses scipy eigh function, which I already learned needs to be translated as eigen(Symmetric(array, :L)) to get matching results.

I am doing the conversion in a very naive way - basically trial-and-error translation line-by-line and checking if the results match. I thought I got everything correct and running one iteration of the loop returned the same results in both languages, but the whole algorithm does not work in Julia.

I think I narrowed it down to divergence in eigen results, that is , first couple of iterations decomposition is either the same or very very close, but then it diverges so much, that the julia version does not want to converge at all (contrary to the python).

Here is an MWE using just eigen (I added PythonCall python version to see it this will change anything, but the results are the same as calling directly from python):

using LinearAlgebra
using PythonCall

sp = pyimport("scipy")
np = pyimport("numpy")

dat =  [1. 2 3; 
        4 5 6; 
        7 8 9; 
]

newdat = dat * dat'
z = eigen(Symmetric(newdat, :L))
for i in 1:10
    newz = z.vectors .+ 1
    z = eigen(Symmetric(newz, :L))
    display(z.vectors)
end

pynewdat = np.array(newdat)
q = sp.linalg.eigh(pynewdat)
for i in 1:10
    newq = q[1] + 1
    q = sp.linalg.eigh(np.array(newq))
    display(q[1])
end

It diverges for me on the third iteration.

My common sense would dictate, that the algorithm should work regardless of variation of 1e-5 range, since both languages are calling the same LAPACK routine, right? Both represent numbers as Float64. So why the difference? Can it be mitigated somehow?

2

Do Python and Julia order the eigenpairs the same way?

3

Do you mean is the order within the returned eigenvalues/vectors the same? Yes it is.
When they work in the same way, the are virtually identical.

Although on that note, I was focusing on the eigenvectors, but now I noticed that even for the first run of eigen the is actually a difference in the first eigenvalue: julia gives 3.313751296132755e-15 and python 9.71445138e-16 but this is on the verge of “epsilon” for floats in Julia (and the algorithm has a step of correcting values smaller than that).

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
5

I think you are making a comment that the MWE doesn’t make mathematical sense, right?
That definitely might be true - this is not representative of the actual code - I just replaced all other operations with the addition for brevity.

This was my attempt to see if and when some numerical differences would pile up while repeatedly calling eigen to check if what I am seeing in my actual code stems from this or did I made a mistake somewhere else.

And still, even if the results wouldn’t be the same, shouldn’t they be close to each other?
It just suprisees me that I can get virtually identical results for 2-3 iterations and then something widly different on the next one.

6

No. Different libraries — or different versions of the same library — can choose different signs for eigenvectors, leading to completely different results in an ill-posed algorithm like yours. It’s not a question of roundoff errors.

It will be about the same until it chooses a different sign, then boom.

3 Likes
7

I see, did not expect this level of indeterminancy from linear algebra, to be honest. :slight_smile:
Thanks for the explanation!

So does this sound probable to you, that because of that one version of the same algorithm would converge but the other one wouldn’t?
Is there anything to be done, to mitigate/minimize this?

8

If your algorithm depends on the overall sign choice of your eigenvectors, then your algorithm is incorrect (in any language!). Fix your algorithm.

(Think about why it depends on the sign. Can this dependency be removed? Or is there a deterministic way to pick a canonical sign in your application? I can’t give you more advice without knowing your actual problem.)

3 Likes
9

Right, so I wanted to better understand what the ICA algorithm does by translating the FastICA implementation from scikit-learn.
I started with the deflation version, which works ok for me (and it does not call eigen in the loop) and then moved to the parallel one. And here, as mentioned, when I run the algorithm on identical data, first couple of iterations are ok (which I took as a sign that I translated all the operations correctly), then it goes sideways - to the point that it does not converge any more. The eigen is used in the routine for symmetric decorrelation (which FastICA from MultivarieteStats does differently, so I had to rely only on python code).
That is why, after check other stuff for some time I thought that maybe the internal handling of calling LAPACK was is some way different etc.

I understand now that identical results are almost impossible to achieve in the long run, but now I kind of have no idea how to trace why it does not converge in the end.

10

Eigenvectors for decorrelation/whitening are perfectly fine — the eigenvector basis is still uncorrelated (still diagonalizes the correlation matrix) if you flip the signs. I don’t know what the parallel FastICA algorithm you are referring to looks like, but perhaps there is some synchronization step you are missing?

In any case, it sounds like there is some aspect of the algorithm you need to understand better in order to implement it properly.

2 Likes
11

Tip: you can have ChatGPT translate the Python code to Julia. Then go fix the errors ChatGPT made.

12

For sure. I just thought I checked all other parts a couple of times already, so this felt like the only aspect I can’t just test and had to ask for help.

To give this thread some closure - the problem was not with the implementation, but with the transition from row-major python code to a column-major. I missed a transpose in one line (and both options gave very similar results on my dummy data, so I didn’t notice).
Code is working, I know a bit more about ICA and eigen decomp., maybe I would only hope it would take me a bit less time. :slight_smile: Thanks again for clarifications!

I actually tried ChatGPT and Github Copilot at some point in the past, but the versions then struggled both with ICA and Julia in general. Probably the upcoming versions will give me a better result, but for now I’ll exercise my “meat brain”.

2 Likes