Precision in numpy: issues while comparing numbers

Question

A bit of background first. I am finding the eigenvalues and eigenvectors of a real symmetric matrix, in which rows sum to 0. More specifically, once I find an eigenvector, I use $argsort$ to find the permutation that sorts one of the eigenvalues and apply the permutation to the original matrix.

Now, I implemented the code in python, using the numpy package. The code itself is recursive, and if it finds a set of values in the eigenvector that are equal, it extracts the symmetric submatrix corresponding to the indices for which we have equal values, and applies the algorithm all over again on this matrix.

While this is all very well, and mostly grunt work, I was caught by surprise when a bunch of indices that should have corresponded to equal entries in the eigenvector were not recognized as having equal values. The problem was that the values were computed to machine precision by some algorithm (possibly Lanczos, but I am not entirely familiar with numpy). This is a sample output, in which I explicitly check the difference between two entries in the eigenvector:

    >>> T=spectral.seriation(A,index)

    columns [ 0  1  2  3  4  5  6  7  8  9 10 11]

    [  3.30289130e-01  -2.75240941e-01  -2.75240941e-01   3.30289130e-01
    -2.75240941e-01   3.30289130e-01  -2.75240941e-01   3.30289130e-01
    3.30289130e-01  -2.75240941e-01  -1.69794463e-16  -2.75240941e-01]

    [ 4  6  9  1  2 11 10  0  5  7  8  3]

    difference   -5.55111512313e-17

The routine seriation() is a recursive function. The array of floats is the eigenvector under consideration, and the array below that gives the sorted order of the columns. Note that columns [4,6,9,1,2,11] have the same value. However, eigenvector and eigenvalue calculations are always approximations, and indeed, when I output the difference between the entry in column 9 and column 2, it is non-zero. Where the algorithm should group [4,6,9,1,2,11], it only groups [4,6,9], and puts the rest in another group, throwing a wrench into the works.

So the question is this: is there a method to perform arbitrary precision calculations in numpy? Failing this, what would be a `good' workaround to this problem?

Also, I should probably mention that it can be mathematically proven that these entries must be equal. That is a property of the matrix, but hopefully not germane to the question.

Answer 1

Doubles are not exactly real numbers [not even rational]. There are infinite number of rationals in every range [well, every range with at least two elements, to be exact], but only a finite number of bits to represent them.
Thus, you should expect some rounding errors for "exact" calculations.

For more inforamtion, you might want to read what every computer scientist should know about floating-point arithmetic

Answer 2

Check numpy.allclose and numpy.isclose functions for testing equality within a tolerance.