# Diverging results of eigen() between Julia and Python

**URL:** <https://discourse.julialang.org/t/diverging-results-of-eigen-between-julia-and-python/97258>\
**Category:** General Usage\
**Tags:** linearalgebra\
**Created:** [April 8, 2023, 8:53pm UTC](https://discourse.julialang.org/t/diverging-results-of-eigen-between-julia-and-python/97258 "2023-04-08T20:53:15Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![mkoculak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkoculak/32/28310_2.png) [@mkoculak](https://discourse.julialang.org/u/mkoculak)\
**Post date:** [April 8, 2023, 8:53pm UTC](https://discourse.julialang.org/t/diverging-results-of-eigen-between-julia-and-python/97258/1 "2023-04-08T20:53:15Z")

</div>

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

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

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [April 8, 2023, 9:04pm UTC](https://discourse.julialang.org/t/diverging-results-of-eigen-between-julia-and-python/97258/2 "2023-04-08T21:04:18Z")

</div>

Do Python and Julia order the eigenpairs the same way?

---

<div class="post-metadata">

**Author:** ![mkoculak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkoculak/32/28310_2.png) [@mkoculak](https://discourse.julialang.org/u/mkoculak)\
**Post date:** [April 8, 2023, 9:14pm UTC](https://discourse.julialang.org/t/diverging-results-of-eigen-between-julia-and-python/97258/3 "2023-04-08T21:14:36Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 8, 2023, 9:20pm UTC](https://discourse.julialang.org/t/diverging-results-of-eigen-between-julia-and-python/97258/4 "2023-04-08T21:20:30Z")

</div>

> [@mkoculak](#):
>
> My common sense would dictate, that the algorithm should work regardless of variation of 1e-5 range

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
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.

---

<div class="post-metadata">

**Author:** ![mkoculak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkoculak/32/28310_2.png) [@mkoculak](https://discourse.julialang.org/u/mkoculak)\
**Post date:** [April 8, 2023, 9:35pm UTC](https://discourse.julialang.org/t/diverging-results-of-eigen-between-julia-and-python/97258/5 "2023-04-08T21:35:07Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 8, 2023, 9:37pm UTC](https://discourse.julialang.org/t/diverging-results-of-eigen-between-julia-and-python/97258/6 "2023-04-08T21:37:11Z")

</div>

> [@mkoculak](#):
>
> And still, even if the results wouldn’t be the same, shouldn’t they be close to each other?

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.

> [@mkoculak](#):
>
> It just surprises me that I can get virtually identical results for 2-3 iterations and then something widly different on the next one.

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

---

<div class="post-metadata">

**Author:** ![mkoculak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkoculak/32/28310_2.png) [@mkoculak](https://discourse.julialang.org/u/mkoculak)\
**Post date:** [April 8, 2023, 9:43pm UTC](https://discourse.julialang.org/t/diverging-results-of-eigen-between-julia-and-python/97258/7 "2023-04-08T21:43:18Z")

</div>

I see, did not expect this level of indeterminancy from linear algebra, to be honest. 🙂  
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?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 8, 2023, 9:48pm UTC](https://discourse.julialang.org/t/diverging-results-of-eigen-between-julia-and-python/97258/8 "2023-04-08T21:48:06Z")

</div>

> [@mkoculak](#):
>
> Is there anything to be done, to mitigate/minimize this?

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

---

<div class="post-metadata">

**Author:** ![mkoculak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkoculak/32/28310_2.png) [@mkoculak](https://discourse.julialang.org/u/mkoculak)\
**Post date:** [April 8, 2023, 10:02pm UTC](https://discourse.julialang.org/t/diverging-results-of-eigen-between-julia-and-python/97258/9 "2023-04-08T22:02:17Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 8, 2023, 10:20pm UTC](https://discourse.julialang.org/t/diverging-results-of-eigen-between-julia-and-python/97258/10 "2023-04-08T22:20:03Z")

</div>

> [@mkoculak](#):
>
> 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).

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.

---

<div class="post-metadata">

**Author:** ![blackeneth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/blackeneth/32/10353_2.png) [@blackeneth](https://discourse.julialang.org/u/blackeneth)\
**Post date:** [April 9, 2023, 6:59pm UTC](https://discourse.julialang.org/t/diverging-results-of-eigen-between-julia-and-python/97258/11 "2023-04-09T18:59:01Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![mkoculak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkoculak/32/28310_2.png) [@mkoculak](https://discourse.julialang.org/u/mkoculak)\
**Post date:** [April 11, 2023, 12:58am UTC](https://discourse.julialang.org/t/diverging-results-of-eigen-between-julia-and-python/97258/12 "2023-04-11T00:58:50Z")

</div>

> [@stevengj](#):
>
> In any case, it sounds like there is some aspect of the algorithm you need to understand better in order to implement it properly.

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. 🙂 Thanks again for clarifications!

> [@blackeneth](#):
>
> Tip: you can have ChatGPT translate the Python code to Julia. Then go fix the errors ChatGPT made.

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”.