Same SciPy.linalg.eigh() eigenvalues are not equal

Question

I am using Python 2.7, NumPy 1.6.2 and SciPy 0.16.0 to calculate the following.

I have created a Hadamard matrix. Then I have created outer products from the 0,2,4-th vectors of the matrix and added them and made the diagonal 0. Then I have computed the eigenvalues using SciPy.linalg.eigh(). It is degenerate. First two eigenvalues are same, i.e., -3. But when I use Python to check it says they are not same. The code is given below.

from scipy import linalg as sp
import numpy
from numpy import linalg as np

def get_outer_product(vector):
    length = len(vector)

    outer_product = [[0 for x in range(length)] for x in range(length)]

    for i in range(0, length):
        for j in range(0, length):
            if i == j:
                outer_product[i][j] = 0
                continue
            outer_product[i][j] = vector[i] * vector[j]
    return outer_product

def test():
    hadamard_matrix = sp.hadamard(8)
    sum_of_outer_products = [map(sum, zip(*t)) for t in zip(get_outer_product(hadamard_matrix[0]), get_outer_product(hadamard_matrix[2]))]
    sum_of_outer_products = [map(sum, zip(*t)) for t in zip(sum_of_outer_products, get_outer_product(hadamard_matrix[4]))]

    e_vals, e_vecs = sp.eigh(sum_of_outer_products)
    print str(e_vals[0]) + " == " + str(e_vals[0]) + "?"
    print e_vals[0] == e_vals[1]

The output is:

Python 2.7 (r27:82525, Jul  4 2010, 09:01:59) [MSC v.1500 32 bit (Intel)] on win32
>>> import scipytest
>>> scipytest.test()
-3.0 == -3.0?
False

What I am doing wrong?

Can you print the actual values of the array and the types. It isn't returning a complex type or something right? Are you sure you're comparing two floats? — Brian Pendleton, Sep 08 '15 at 21:33

Salvador Dali · Answer 1 · 2016-01-09T01:27:51.690

2

Floating point problem.

When working with floats, always be careful with equal operator and remember that there might be precision problems and use safe comparison:

print str(e_vals[0]) + " == " + str(e_vals[0]) + "?"
print e_vals[0] == e_vals[1]
print numpy.isclose(e_vals[0], e_vals[1])

By the way, it returned true in both cases on my machine:

-3.0 == -3.0?
True
True

edited Jan 09 '16 at 01:27

answered Sep 08 '15 at 21:33

Salvador Dali

214,103
147
703
753

I am using NumPy 1.6.2. – Omar Shehab Sep 08 '15 at 23:04
1

@OmarShehab then you either update your library or you can implement your own version of `isclose()`, which just checks whether one floating point is approximately equal to another floating point especially seeing that documentation explains how `isclose` is implemented: `absolute(a - b) <= (atol + rtol * absolute(b))` – Salvador Dali Sep 08 '15 at 23:36
I cannot upgrade NumPy because of other dependencies. Do you have any idea why the output of print e_vals[0] == e_vals[1] is different at your end? – Omar Shehab Sep 09 '15 at 00:26
1

@OmarShehab these are just guesses, but `eig*` is using some numerical convergence algorithm and it does not always converge to absolutely the same point (may be this is -3.000000000000000018) and due to floating point inaccuracy the two numbers are not the same. It is actually good that you have noticed it now, and not in the middle of a way bigger project. I suggest never to use floating point equal operator and implement your own version of equality for floating points. – Salvador Dali Sep 09 '15 at 00:32
Does it mean that the eig() function has been modified in your version? – Omar Shehab Sep 09 '15 at 00:35
1

@OmarShehab I highly doubt that there were any changes to `eig` function from 1.6.2 to 1.7 version. It may be just converged exactly to -3. But even in this case I would not use the `e_vals[0] == e_vals[1]` code and will still use isclose or my own isclose implementation. – Salvador Dali Sep 09 '15 at 00:38

Same SciPy.linalg.eigh() eigenvalues are not equal

1 Answers1