Numpy eigenvectors aren't eigenvectors?
Asked Answered
D

1

7

I was doing some matrix calculations and wanted to calculate the eigenvalues and eigenvectors of this particular matrix:

enter image description here

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.)

Dari answered 22/4, 2015 at 13:28 Comment(0)
D
9

The eigen vectors are stored as column vectors as described here. So you have to use vecs[:,0] instead vecs[0]

For example this here works for me (I use eig because A is not symmetric)

import numpy as np
import numpy.linalg as LA
import numpy.random

A = numpy.random.randint(10,size=(4,4))
# array([[4, 7, 7, 7],
#        [4, 1, 9, 1],
#        [7, 3, 7, 7],
#        [6, 4, 6, 5]])

eval,evec = LA.eig(A)

evec[:,0]
# array([ 0.55545073+0.j,  0.37209887+0.j,  0.56357432+0.j,  0.48518131+0.j])

np.dot(A,evec[:,0]) / eval[0]
# array([ 0.55545073+0.j,  0.37209887+0.j,  0.56357432+0.j,  0.48518131+0.j])
Dateless answered 22/4, 2015 at 13:37 Comment(2)
That's it. -I feel like quite the fool- Thank you a lot, that solves everything!Dari
@Dari Don't feel like a fool, the numpy/scipy eigenvector output as columns instead of rows is highly irregular and not well documented...having run into this myself and explored the internet on this topic it's a very common mistake.Arnst

© 2022 - 2024 — McMap. All rights reserved.