Here are two assumptions about eigenvectors and eigenvalues of square matrices. I believe that both are true:
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.
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 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 marked with an arrow (<==). Here you find three different values, but the values should be exactly equal. And if you have a second look, you will see, that only in the 1st line all three values are equal. In all other lines you will find 2 or 3 different values.
numpy.linalg.eig | |
scipy.linalg.eig | numpy.linalg.eigh | scipy.linalg.eigh
------------------+---------------------+-------------------
-0.354+0.000i | -0.354+0.000i | -0.354+0.000i
0.913+0.000i | 0.000+0.000i | 0.000+0.000i
0.204+0.000i | 0.000+0.000i | 0.000+0.000i
-0.013+0.016i | -0.086+0.000i | -0.225+0.000i <===
-0.013-0.016i | 0.905+0.000i | 0.882+0.000i <===
0.160+0.000i | -0.025+0.000i | 0.000+0.000i <===
-0.000+0.000i | 0.073+0.000i | 0.065+0.000i <===
0.130+0.000i | 0.205+0.000i | -0.205+0.000i
Here are my questions:
- How is this possible?
- Are these bugs?
- Is one of the results correct?
- If there is a method that delivers correct results: Which is it?
p.s: Here are relevant version informations:
- I did run this code on an iMac (macOS Catalina Version 10.15.7)
- The python version is 3.8.5
- The version of numpy is 1.19.5
- The version of scipy is 1.6.0
This is the output of numpy.show_config()
(as requested in a comment):
blas_mkl_info:
NOT AVAILABLE
blis_info:
NOT AVAILABLE
openblas_info:
libraries = ['openblas', 'openblas']
library_dirs = ['/usr/local/lib']
language = c
define_macros = [('HAVE_CBLAS', None)]
blas_opt_info:
libraries = ['openblas', 'openblas']
library_dirs = ['/usr/local/lib']
language = c
define_macros = [('HAVE_CBLAS', None)]
lapack_mkl_info:
NOT AVAILABLE
openblas_lapack_info:
libraries = ['openblas', 'openblas']
library_dirs = ['/usr/local/lib']
language = c
define_macros = [('HAVE_CBLAS', None)]
lapack_opt_info:
libraries = ['openblas', 'openblas']
library_dirs = ['/usr/local/lib']
language = c
define_macros = [('HAVE_CBLAS', None)]
None
ADDENDUM (added 1 day after the question was asked)
Reaction to comments:
complex eigenvectors of real symmetric matrices
@Rethipher: Thank you! I did read and understand the question you linked to (Can a real symmetric matrix have complex eigenvectors?), and I also did read the answers, but I didn’t understand them. Did they say “yes” or “no”? (rhetoric question, no need to answer, see next line)
@Mark Dickinson & @bnaecker: Thank you for making clear, that my assumption was wrong.real symmetric matrices vs. Hermitian matrices
@bnaecker: The set of real numbers is a subset of the set of complex numbers. Those complex numbers which are equal to their own complex conjugate are called real. So, the set of real symmetric matrices is a subset of Hermitian matrices. This is important, because numpy.linalg.eigh and scipy.linalg.eigh are designed to handle Hermitian matrices. And because every real symmetric matrix is a Hermitian matrix, those modules also can be used for my purposes.mixing up rows and columns
@Mark Dickinson & @bnaecker: Thank you, I think you are right. Also the documentations says so, I should have read it more carefully. But even if you compare columns instead of rows you will still find that the 4 methods produce 3 different results. But if the result contains a 7-dimensional subspace that can be described with 7 real basis vectors only, I still find it strange, that an algorithm produces a complex basis.“a bug would be surprising”
@bnaecker: This is true, but surprising bugs do exist. (Like Heartbleed and some others.) So, this is not really an argument.“I get reals” - “your sample matrix doesn't contain floats”
@Stef & @JohanC: Sorry, you didn’t read my program carefully enough. I added a value of1e-12toA[5,7]andA[7,5]to simulate tiny rounding errors that appear inevitably in my real app before it comes to the calculation of eigenvalues and eigenvectors. (What I’ve posted here is just a tiny test program, just big enough to demonstrate the issue.)
And you are right, Stef: Without adding this tiny noise, I also get real results. But only a tiny change of one millionth of one millionth makes such a big difference, and I can't understand why.
Reaction to DavidB2013’s answer :
I tried the tool you suggested, and I got different results. I think you also forgot to add that little noise of 1e-12 to A[5,7] and A[7,5]. However, all results are still real. I did get these eigenvalues:
12.000000000000249
20
20.00000000000075
19.999999999999
20
20
20
20
and these eigenvectors:
0.3535533905932847 0.9128505045937204 0.20252576206455747 0.002673672081814904 -0.09302397289286794 -0.09302397289286794 -0.09302397289286794 -0.09302397289286794
0.3535533905932848 -0.18259457246238117 0.20444330131542393 -0.00009386949436945406 -0.20415317121194954 -0.20415317121194954 -0.20415317121194954 -0.20415317121194954
0.3535533905932848 -0.18259457246238117 0.20444330131542393 -0.00009386949436945406 -0.20415317121194954 -0.20415317121194954 -0.20415317121194954 0.9080920678356449
0.3535533905932848 -0.18259457246238117 0.20444330131542393 -0.00009386949436945406 -0.20415317121194954 0.9080920678356449 -0.20415317121194954 -0.20415317121194954
0.3535533905932848 -0.18259457246238117 0.20444330131542393 -0.00009386949436945406 0.9080920678356449 -0.20415317121194954 -0.20415317121194954 -0.20415317121194954
0.35355339059324065 -0.00011103276380548543 -0.6116010247648269 0.7060012169461334 0.0005790869815273477 0.0005790869815273477 0.0005790869815273477 0.0005790869815273477
0.3535533905932848 -0.18259457246238117 0.20444330131542393 -0.00009386949436945406 -0.20415317121194954 -0.20415317121194954 0.9080920678356449 -0.20415317121194954
0.35355339059324054 0.0002333904819895115 -0.6131412438770024 -0.7082055415560993 0.0009655029234935232 0.0009655029234935232 0.0009655029234935232 0.0009655029234935232
Only the vector for eigenvalue 12 has the same values as your calculation: (There is a difference of approx. 1.1e-14 in 6 dimensions and 3.3e-14 in the two other dimensions, but I count this as rounding error.) All other vectors are significantly different (the smallest differences are of the size of 0.02). It puzzles me, that a tiny rounding error of 1e-12 in just 2 elements of the input matrix can produce so big differences.
I calculated the eigenvalues with another method (with the help of https://www.wolframalpha.com), and when I didn’t add the tiny delta values, which should simulate rounding errors, I only get two different eigenvalues which are 12 and 20.
The characteristic polynomial of the given matrix is:
(20 - λ)^7 * (12 - λ)
So, it has one root at λ=12 and 7 roots at λ=20 and these 8 roots are the 8 eigenvalues. All of them real numbers.
When I add the tiny delta values, I get this characteristic polynomial:
(20 - λ)^5 * (19999999999999/1000000000000 - λ) * (1000000000000 λ^2 - 32000000000001 λ + 240000000000014)/1000000000000
It has these roots:
λ=12.00000000000024999999999998 (rounded)
λ=19.999999999999 (exact value)
λ=20 (exact value)
λ=20 (exact value)
λ=20 (exact value)
λ=20 (exact value)
λ=20 (exact value)
λ=20.00000000000075000000000002 (rounded)
And again all 8 eigenvalues are real numbers.
Then I calculated the eigenvectors. Without adding 1e-12 I get this results:
Vector for eigenvalue 12:
v = (1,1,1,1,1,1,1,1)
The length of this vector is sqrt(8), and if you multiply the vector with 1/sqrt(8), you get exactly the result from the other calculations (0.35355339 in each dimension).
But the seven eigenvectors for eigenvalue 20 are very different. They are:
(-1,1,0,0,0,0,0,0)
(-1,0,1,0,0,0,0,0)
(-1,0,0,1,0,0,0,0)
(-1,0,0,0,1,0,0,0)
(-1,0,0,0,0,1,0,0)
(-1,0,0,0,0,0,1,0)
(-1,0,0,0,0,0,0,1)
Even if you bring them to the length 1, they are different from all other results and it is very easy to see that they are correct. The other results are also correct, but I would prefer these simple results.
I also calculated the eigenvalues for the version with the tiny noise. All 8 vectors are so close to the noise-less results, that even Wolfram Alpha rounded them to exactly the same values as before. And this is exactly the behavior that I would expect from an algorithm that calculates eigenvalues and eigenvectors:
- Tiny variations in the input should - when ever it is possible - return tiny variations in the results.
eighresults. The eigenvector for12is[-0.35355339 -0.35355339 -0.35355339 -0.35355339 -0.35355339 -0.35355339 -0.35355339 -0.35355339](the first column of the second output ofeigh). The eigenspace for20of the originalA(without the delta addition) is of course 7-dimensional, so there are many ways to choose a basis of eigenvectors. – Ruffnereigdoesn't sort its output based on eigenvalue, whileeighdoes. And the routines are designed for different things, so slightly different results should not be surprising. LAPACK is one of the most widely-used and heavily-tested libraries on the planet; a bug would be surprising indeed. – DesignateAhas an eigenvalue of exactly20with eigenspace of dimension five, along with two distinct eigenvalues that are close to but not exactly equal to20, and an eigenvalue close to (but again not exactly equal to)12.) – Ruffnerreals in all four cases (version 1.19.4, 1.5.4, cblas) – Choicblasand Apple's own BLAS/LINALG libraries (which I presume NumPy on macOS links against)? To the OP: can you show the output ofnumpy.show_config()? – Ruffnersympy.Matrix(A).eigenvects()gives[-1, 1, 0, 0, 0, 0, 0, 0]; [-1, 0, 1, 0, 0, 0, 0, 0]; [-1, 0, 0, 1, 0, 0, 0, 0]; [-1, 0, 0, 0, 1, 0, 0, 0]; [-1, 0, 0, 0, 0, 1, 0, 0]; [-1, 0, 0, 0, 0, 0, 1, 0]; [-1, 0, 0, 0, 0, 0, 0, 1]as the eigenvectors corresponding to 20. – Latton1e-12toA[5,7]andA[7,5]– Choi