numpy - Eigenvalues in Python: A Bug?

Question

Here are two assumptions about eigenvectors and eigenvalues of square matrices. I believe that both are true:

1. If a matrix is symmetric and contains only real values, then it is a Hermitian matrix, and then all eigenvalues should be real numbers and all components of all eigenvectors should also be real numbers. No complex numbers should appear in the results when you calculate eigenvectors and eigenvalues from Hermitian matrices.

2. The eigenvector of a given eigenvalue, calculated from a given matrix should always point into a direction that is determined only by the matrix and the eigenvalue. The algorithm used to calculate it has no influence on the result, as long as the algorithm is implemented correctly.

But both assumptions do not hold when you use standard libraries in Python to calculate eigenvectors and eigenvalues. Do those methods contain bugs?

There are four different methods to calculate eigenvalues and eigenvectors from Hermitian matrices:

1. numpy.linalg.eig
2. scipy.linalg.eig
3. numpy.linalg.eigh
4. scipy.linalg.eigh

#1 and #2 can be used for any square matrix (including Hermitian matrices).
#3 and #4 are made for Hermitian matrices only. As far as I did understand their purpose is just that they run faster, but the results should be the same (as long as the input is really Hermitian).

But the four methods deliver three different results for the very same input. Here is the program that I used to test all four methods:

#!/usr/bin/env python3

import numpy as np
import scipy.linalg as la

A = [
    [19, -1, -1, -1, -1, -1, -1, -1],
    [-1, 19, -1, -1, -1, -1, -1, -1],
    [-1, -1, 19, -1, -1, -1, -1, -1],
    [-1, -1, -1, 19, -1, -1, -1, -1],
    [-1, -1, -1, -1, 19, -1, -1, -1],
    [-1, -1, -1, -1, -1, 19, -1, -1],
    [-1, -1, -1, -1, -1, -1, 19, -1],
    [-1, -1, -1, -1, -1, -1, -1, 19]
]

A = np.array(A, dtype=np.float64)

delta = 1e-12
A[5,7] += delta
A[7,5] += delta

if np.array_equal(A, A.T):
    print('input is symmetric')
else:
    print('input is NOT symmetric')

methods = {
    'np.linalg.eig'  : np.linalg.eig,
    'la.eig'         : la.eig,
    'np.linalg.eigh' : np.linalg.eigh,
    'la.eigh'        : la.eigh
}

for name, method in methods.items():

    print('============================================================')
    print(name)
    print()

    eigenValues, eigenVectors = method(A)
    
    for i in range(len(eigenValues)):
        print('{0:6.3f}{1:+6.3f}i '.format(eigenValues[i].real, eigenValues[i].imag), end=' |  ')
        line = eigenVectors[i]
        for item in line:
            print('{0:6.3f}{1:+6.3f}i '.format(item.real, item.imag), end='')
        print()

    print('---------------------')

    for i in range(len(eigenValues)):
        if eigenValues[i].imag == 0:
            print('real    ', end=' |  ')
        else:
            print('COMPLEX ', end=' |  ')
        line = eigenVectors[i]
        for item in line:
            if item.imag == 0:
                print('real    ', end='')
            else:
                print('COMPLEX ', end='')
        print()

    print()

And here is the output it produces:

input is symmetric
============================================================
np.linalg.eig

12.000+0.000i  |  -0.354+0.000i  0.913+0.000i  0.204+0.000i -0.013+0.016i -0.013-0.016i  0.160+0.000i -0.000+0.000i  0.130+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.183+0.000i  0.208+0.000i  0.379-0.171i  0.379+0.171i -0.607+0.000i  0.000+0.000i -0.138+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.182+0.000i  0.203+0.000i -0.468-0.048i -0.468+0.048i  0.153+0.000i  0.001+0.000i -0.271+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.182+0.000i  0.203+0.000i  0.657+0.000i  0.657-0.000i  0.672+0.000i -0.001+0.000i  0.617+0.000i 
20.000-0.000i  |  -0.354+0.000i -0.182+0.000i  0.203+0.000i -0.276+0.101i -0.276-0.101i -0.361+0.000i  0.001+0.000i -0.644+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.001+0.000i -0.612+0.000i -0.001+0.000i -0.001-0.000i  0.001+0.000i  0.706+0.000i -0.000+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.182+0.000i  0.203+0.000i -0.276+0.101i -0.276-0.101i -0.018+0.000i -0.000+0.000i  0.306+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.001+0.000i -0.612+0.000i -0.001+0.000i -0.001-0.000i  0.001+0.000i -0.708+0.000i  0.000+0.000i 
---------------------
real     |  real    real    real    COMPLEX COMPLEX real    real    real    
real     |  real    real    real    COMPLEX COMPLEX real    real    real    
real     |  real    real    real    COMPLEX COMPLEX real    real    real    
COMPLEX  |  real    real    real    real    real    real    real    real    
COMPLEX  |  real    real    real    COMPLEX COMPLEX real    real    real    
real     |  real    real    real    COMPLEX COMPLEX real    real    real    
real     |  real    real    real    COMPLEX COMPLEX real    real    real    
real     |  real    real    real    COMPLEX COMPLEX real    real    real    

============================================================
la.eig

12.000+0.000i  |  -0.354+0.000i  0.913+0.000i  0.204+0.000i -0.013+0.016i -0.013-0.016i  0.160+0.000i -0.000+0.000i  0.130+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.183+0.000i  0.208+0.000i  0.379-0.171i  0.379+0.171i -0.607+0.000i  0.000+0.000i -0.138+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.182+0.000i  0.203+0.000i -0.468-0.048i -0.468+0.048i  0.153+0.000i  0.001+0.000i -0.271+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.182+0.000i  0.203+0.000i  0.657+0.000i  0.657-0.000i  0.672+0.000i -0.001+0.000i  0.617+0.000i 
20.000-0.000i  |  -0.354+0.000i -0.182+0.000i  0.203+0.000i -0.276+0.101i -0.276-0.101i -0.361+0.000i  0.001+0.000i -0.644+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.001+0.000i -0.612+0.000i -0.001+0.000i -0.001-0.000i  0.001+0.000i  0.706+0.000i -0.000+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.182+0.000i  0.203+0.000i -0.276+0.101i -0.276-0.101i -0.018+0.000i -0.000+0.000i  0.306+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.001+0.000i -0.612+0.000i -0.001+0.000i -0.001-0.000i  0.001+0.000i -0.708+0.000i  0.000+0.000i 
---------------------
real     |  real    real    real    COMPLEX COMPLEX real    real    real    
real     |  real    real    real    COMPLEX COMPLEX real    real    real    
real     |  real    real    real    COMPLEX COMPLEX real    real    real    
COMPLEX  |  real    real    real    real    real    real    real    real    
COMPLEX  |  real    real    real    COMPLEX COMPLEX real    real    real    
real     |  real    real    real    COMPLEX COMPLEX real    real    real    
real     |  real    real    real    COMPLEX COMPLEX real    real    real    
real     |  real    real    real    COMPLEX COMPLEX real    real    real    

============================================================
np.linalg.eigh

12.000+0.000i  |  -0.354+0.000i  0.000+0.000i  0.000+0.000i -0.086+0.000i  0.905+0.000i -0.025+0.000i  0.073+0.000i  0.205+0.000i 
20.000+0.000i  |  -0.354+0.000i  0.000+0.000i -0.374+0.000i  0.149+0.000i -0.236+0.000i -0.388+0.000i  0.682+0.000i  0.206+0.000i 
20.000+0.000i  |  -0.354+0.000i  0.001+0.000i  0.551+0.000i  0.136+0.000i -0.180+0.000i  0.616+0.000i  0.317+0.000i  0.201+0.000i 
20.000+0.000i  |  -0.354+0.000i  0.001+0.000i -0.149+0.000i  0.719+0.000i -0.074+0.000i -0.042+0.000i -0.534+0.000i  0.207+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.005+0.000i  0.505+0.000i -0.386+0.000i -0.214+0.000i -0.556+0.000i -0.274+0.000i  0.203+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.707+0.000i -0.004+0.000i  0.002+0.000i  0.001+0.000i  0.002+0.000i -0.000+0.000i -0.612+0.000i 
20.000+0.000i  |  -0.354+0.000i  0.003+0.000i -0.529+0.000i -0.535+0.000i -0.203+0.000i  0.398+0.000i -0.262+0.000i  0.203+0.000i 
20.000+0.000i  |  -0.354+0.000i  0.707+0.000i  0.001+0.000i  0.001+0.000i  0.000+0.000i -0.005+0.000i -0.001+0.000i -0.612+0.000i 
---------------------
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    

============================================================
la.eigh

12.000+0.000i  |  -0.354+0.000i  0.000+0.000i  0.000+0.000i -0.225+0.000i  0.882+0.000i  0.000+0.000i  0.065+0.000i -0.205+0.000i 
20.000+0.000i  |  -0.354+0.000i  0.000+0.000i -0.395+0.000i  0.332+0.000i -0.156+0.000i  0.227+0.000i  0.701+0.000i -0.205+0.000i 
20.000+0.000i  |  -0.354+0.000i  0.001+0.000i  0.612+0.000i  0.011+0.000i -0.204+0.000i -0.597+0.000i  0.250+0.000i -0.200+0.000i 
20.000+0.000i  |  -0.354+0.000i  0.001+0.000i -0.086+0.000i  0.689+0.000i  0.030+0.000i -0.054+0.000i -0.589+0.000i -0.205+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.005+0.000i  0.413+0.000i -0.264+0.000i -0.245+0.000i  0.711+0.000i -0.165+0.000i -0.205+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.707+0.000i -0.004+0.000i -0.000+0.000i  0.001+0.000i -0.002+0.000i -0.001+0.000i  0.612+0.000i 
20.000+0.000i  |  -0.354+0.000i  0.003+0.000i -0.540+0.000i -0.542+0.000i -0.309+0.000i -0.290+0.000i -0.261+0.000i -0.205+0.000i 
20.000+0.000i  |  -0.354+0.000i  0.707+0.000i  0.001+0.000i -0.000+0.000i  0.001+0.000i  0.005+0.000i -0.001+0.000i  0.612+0.000i 
---------------------
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real 

As you can see, numpy.linalg.eig and scipy.linalg.eig produce complex numbers in their output, but they shouldn't. This could be accepted as some kind of rounding error, if the magnitude of the imaginary part would by tiny compared to the magnitude of the real part. But this is not the case. One of the numbers that are produced is -0.013+0.016i. Here the imaginary part has an even higher magnitude than the real part.

Even worse: The four methods produce three different results.

All four methods calculate only once an eigenvalue of 12 and 7 times an eigenvalue of 20. And all eigenvectors always have the length 1. This means, all four methods should produce the very same eigenvector for eigenvalue 12. But only numpy.linalg.eig and scipy.linalg.eig produce the same output.

Here are the components of the eigenvector for eigenvalue 12. Have a closer look to the lines

Answer 1

As far as I know, assumption 1 is correct, but assumption 2 is not.

A Real Symmetric matrix produces eigenvalues and eigenvectors that are real only.

However, for a given eigenvalue, the associated eigenvector isn't necessarily unique.

Furthermore, round-off error shouldn't be so significant for a matrix that actually isn't that big, or contain numbers that aren't very small.

For comparison, I ran your test matrix through a JavaScript version of RG.F (Real General, from the EISPACK Library): Eigenvalues and Eigenvectors Calculator

Here is the output:

Eigenvalues:

   20
   12
   20
   20
   20
   20
   20
   20

Eigenvectors:

 0.9354143466934854     0.35355339059327395     -0.021596710639534     -0.021596710639534     -0.021596710639534     -0.021596710639534     -0.021596710639533997     -0.021596710639533997
-0.1336306209562122     0.3535533905932738     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797
-0.1336306209562122     0.3535533905932738     0.9286585574999623     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797
-0.1336306209562122     0.3535533905932738     -0.15117697447673797     0.9286585574999623     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797
-0.1336306209562122     0.3535533905932738     -0.15117697447673797     -0.15117697447673797     0.9286585574999623     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797
-0.1336306209562122     0.3535533905932738     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     0.9286585574999623     -0.15117697447673797     -0.15117697447673797
-0.1336306209562122     0.3535533905932738     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     0.9286585574999622     -0.15117697447673797
-0.1336306209562122     0.3535533905932738     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     0.9286585574999622

No imaginary components.

To confirm, or deny, the validity of results, you could always write a small program that plugs the results back into the original equation. Simple matrix and vector multiplication. Then you'd know for sure whether or not the outputs are correct. Or, if they are wrong, how far away from correct answers they are.