I was doing some matrix calculations and wanted to calculate the eigenvalues and eigenvectors of this particular matrix:
I found its eigenvalues and eigenvectors analytically and wanted to confirm my answer using numpy.linalg.eigh
, since this matrix is symmetric. Here is the problem: I find the expected eigenvalues, but the corresponding eigenvectors appear to be not eigenvectors at all
Here is the little piece of code I used:
import numpy as n
def createA():
#create the matrix A
m=3
T = n.diag(n.ones(m-1.),-1.) + n.diag(n.ones(m)*-4.) +\
n.diag(n.ones(m-1.),1.)
I = n.identity(m)
A = n.zeros([m*m,m*m])
for i in range(m):
a, b, c = i*m, (i+1)*m, (i+2)*m
A[a:b, a:b] = T
if i < m - 1:
A[b:c, a:b] = A[a:b, b:c] = I
return A
A = createA()
ev,vecs = n.linalg.eigh(A)
print vecs[0]
print n.dot(A,vecs[0])/ev[0]
So for the first eigenvalue/eigenvector pair, this yields:
[ 2.50000000e-01 5.00000000e-01 -5.42230975e-17 -4.66157689e-01
3.03192985e-01 2.56458619e-01 -7.84539156e-17 -5.00000000e-01
2.50000000e-01]
[ 0.14149052 0.21187998 -0.1107808 -0.35408209 0.20831606 0.06921674
0.14149052 -0.37390646 0.18211242]
In my understanding of the Eigenvalue problem, it appears that this vector doesn't suffice the equation A.vec = ev.vec, and that therefore this vector is no eigenvalue at all.
I am pretty sure the matrix A itself is correctly implemented and that there is a correct eigenvector. For example, my analytically derived eigenvector:
rvec = [0.25,-0.35355339,0.25,-0.35355339,0.5,-0.35355339,0.25,
-0.35355339,0.25]
b = n.dot(A,rvec)/ev[0]
print n.allclose(real,b)
yields True
.
Can anyone, by any means, explain this strange behaviour? Am I misunderstanding the Eigenvalue problem? Might numpy be erroneous?
(As this is my first post here: my apologies for any unconventionalities in my question. Thanks you in advance for your patience.)