Menu
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:
p.s: Here are relevant version informations:
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
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 of 1e-12 to A[5,7] and A[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:
372
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.
99
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With
