CODEWITHSUNDEEP
pythonnumpyfloating-pointlinear-algebra
jarondl
jarondl

Reputation: 1663

Precision of numpy's eigenvaluesh

First I find the eigenvalues of a (4000x4000) matrix by using numpy.linalg.eigvalsh. Then, I change the boundary conditions, expecting only a minor change in the eigenvalues.

Subtracting the eigenvalues is vulnerable to floating point errors, so I've used some relative tolerance.

Now say I have a eigenvalue A = 1.0001e-10, and another B = 1.0050e-10. According to my humble knowledge of floating-point arithmetic, A - B != 0. The problem is, that these numbers come from linear algebra calculations involving many orders of magnitude. Other eigenvalues might for example be of order 1.

The question is, what is the precision of eigenvalues calculated using numpy.linalg.eigvalsh? Is this precision relative to the value (A * eps), or is it relative to the largest eigenvalue? or perhaps relative to elements of the original matrix?

For example, this matrix:

1      1e-20
1e-20  3

gives the same eigenvalues as this:

1     1e-5
1e-5  3

Upvotes: 3

Views: 5015

Answers (2)

tiago
tiago

Reputation: 23502

First, the solver is not exact. Second, your example matrices are poorly conditioned: the diagonal elements are orders of magnitude larger than the off-diagonal ones. This is always going to cause numerical issues.

From simple algebra, the determinant of the second matrix is (1 * 3) - (1e5 * 1e5) = 3 - 1e-10. You can already see that the precision problem is actually twice as large the precision of your smallest element. (The same applies for the eigenvalues.) Even though linalg is using double precision, because the solver is approximate you get the same answer. If you change the small values to 1e-3 you start seeing a difference, because now the precision is on the order of the numerical approximation.

This specific problem has been asked before. You can see in this answer how to use sympy to solve the eigenvalues with arbitrary precision.

Upvotes: 0

aka.nice
aka.nice

Reputation: 9512

I'm not sure if Lapack is used underneath eigvalsh, but this might be of interest:

Lapack error bounds for the Symmetric/Unsymmetric eigenproblem:

http://www.netlib.org/lapack/lug/node89.html

http://www.netlib.org/lapack/lug/node91.html

Upvotes: 3

Related Questions