5

I am finding that scipy.linalg.eig sometimes gives inconsistent results. But not every time.

>>> import numpy as np
>>> import scipy.linalg as lin
>>> modmat=np.random.random((150,150))
>>> modmat=modmat+modmat.T  # the data i am interested in is described by real symmetric matrices
>>> d,v=lin.eig(modmat)
>>> dx=d.copy()
>>> vx=v.copy()
>>> d,v=lin.eig(modmat)
>>> np.all(d==dx)
False
>>> np.all(v==vx)
False
>>> e,w=lin.eigh(modmat)
>>> ex=e.copy()
>>> wx=w.copy()
>>> e,w=lin.eigh(modmat)
>>> np.all(e==ex)
True
>>> e,w=lin.eigh(modmat)
>>> np.all(e==ex)
False

While I am not the greatest linear algebra wizard, I do understand that the eigendecomposition is inherently subject to weird rounding errors, but I don't understand why repeating the computation would result in a different value. But my results and reproducibility are varying.

What exactly is the nature of the problem -- well, sometimes the results are acceptably different, and sometimes they aren't. Here are some examples:

>>> d[1]
(9.8986888573772465+0j)
>>> dx[1]
(9.8986888573772092+0j)

The difference above of ~3e-13 does not seem like an enormously big deal. Instead, the real problem (at least for my present project) is that some of the eigenvalues cannot seem to agree on the proper sign.

>>> np.all(np.sign(d)==np.sign(dx))
False
>>> np.nonzero(np.sign(d)!=np.sign(dx))
(array([ 38,  39,  40,  41,  42,  45,  46,  47,  79,  80,  81,  82,  83,
    84, 109, 112]),)
>>> d[38]
(-6.4011617320002525+0j)
>>> dx[38]
(6.1888785138080209+0j)

Similar code in MATLAB does not seem to have this problem.

aestrivex
  • 5,170
  • 2
  • 27
  • 44
  • I tried to reproduce this using NumPy 1.6.1 / SciPy 0.10.1, but couldn't. – NPE Feb 04 '13 at 20:50
  • I am using numpy 1.6.1 and scipy 0.10.0. Also, when I didn't use copy() I wasn't able to produce this error (but it does exist in my larger application where something similar to copy() is going on). Which doesn't mean much, because it is mindbogglingly inconsistent. – aestrivex Feb 04 '13 at 20:53

2 Answers2

5

The eigenvalue decompositions satisfy A V = V Lambda, which is all what is guaranteed --- for instance the order of the eigenvalues is not.

Answer to the second part of your question:

Modern compilers/linear algebra libraries produce/contain code that does different things depending on whether the data is aligned in memory on (e.g.) 16-byte boundaries. This affects rounding error in computations, as floating point operations are done in a different order. Small changes in rounding error can then affect things such as ordering of the eigenvalues if the algorithm (here, LAPACK/xGEEV) does not guarantee numerical stability in this respect.

(If your code is sensitive to things like this, it is incorrect! Running e.g. it on a different platform or different library version would lead to a similar problem.)

The results usually are quasi-deterministic --- for instance you get one of 2 possible results, depending if the array happens to be aligned in memory or not. If you are curious about alignment, check A.__array_interface__['data'][0] % 16.

See http://www.nccs.nasa.gov/images/FloatingPoint_consistency.pdf for more

pv.
  • 33,875
  • 8
  • 55
  • 49
  • It would be extremely useful to include a solution the problem here in addition to the explanation. I don't know how/where to use `A.__array_interface__['data'][0] % 16`, and basically, given that this is my problem, how do I change my code so that it is not sensitive in this way? – Aaron Bramson Jan 12 '16 at 09:47
2

I think your problem is you are expecting the eigenvalues to be returned in a particular order, and they don't always come out the same. Sort them, and you'll be on your way. If I run your code to generate d and dx with eig I get the following:

>>> np.max(d - dx)
(19.275224236664116+0j)

But...

>>> d_i = np.argsort(d)
>>> dx_i = np.argsort(dx)
>>> np.max(d[d_i] - dx[dx_i])
(1.1368683772161603e-13+0j)
Jaime
  • 65,696
  • 17
  • 124
  • 159
  • Ah. That's not quite the right answer, but it's close enough that it pointed me in the right direction. What was going on in my program is that somewhere else in my code I was indexing the eigenvectors at locations where they weren't. It did help me fix my bug (so thank you) but it doesn't explain the question I posed, which is why this calculation provides different answers -- or a different ordering -- when run multiple times. That is, okay, the eigenvalues are not guaranteed to be in any particular order, but if you give it the same matrix as input why would they change? – aestrivex Feb 04 '13 at 22:50