Python eig for generalized eigenvalue does not return correct eigenvectors

Question

Trying to use scipy's linalg.eig to solve a generalized eigenvalue problem. I then check the solution I get and it does not seem like proper eigenvectors were returned. Also, the documentation suggests vectors returned are normalized and this is not the case (though that doesn't bother me that much).

Here are sample matrices:

>>> a
array([[  2.05630374e-01,   8.89584493e-10,  -1.46171715e-06],
       [  8.89584493e-10,   2.38374743e-02,   9.43440334e-06],
       [ -1.46171715e-06,   9.43440334e-06,   1.39685787e-02]])
>>> b
array([[ 0.22501692, -0.07509864, -0.05774453],
       [-0.07509864,  0.02569336,  0.01976284],
       [-0.05774453,  0.01976284,  0.01524993]])

Running eig I get:

>>> w,v = linalg.eig(a,b)
>>> w
array([  3.08431414e-01+0.j,   5.31170281e+01+0.j,   6.06298605e+02+0.j])
>>> v
array([[-0.26014092, -0.46277857, -0.0224057 ],
       [ 0.76112351, -0.59384527, -0.83594841],
       [ 1.        , -1.        ,  1.        ]])

And then testing the result:

>>> a*v[:,0]
array([[ -5.34928750e-02,   6.77083674e-10,  -1.46171715e-06],
       [ -2.31417329e-10,   1.81432622e-02,   9.43440334e-06],
       [  3.80252446e-07,   7.18074620e-06,   1.39685787e-02]])
>>> w[0]*b*v[:,0]
array([[-0.01805437+0.j, -0.01762974+0.j, -0.01781023+0.j],
       [ 0.00602559-0.j,  0.00603163+0.j,  0.00609548+0.j],
       [ 0.00463317-0.j,  0.00463941+0.j,  0.00470356+0.j]])

I thought those two will be equal but they are not... I also tried using eigh instead with no success. Would appreciate any help, I'm obviously missing something.

Answer 1

You can see what's going on by looking at the shape of your output. Your a*v[:,0] should give a vector, so why are you getting a 3x3 array? Answer: because you're not doing matrix multiplication, you're doing component-wise array multiplication.

IOW, you did

>>> a * v[:,0]
array([[ -5.34928759e-02,   6.77083679e-10,  -1.46171715e-06],
       [ -2.31417334e-10,   1.81432623e-02,   9.43440334e-06],
       [  3.80252453e-07,   7.18074626e-06,   1.39685787e-02]])
>>> w[0] * b * v[:,0]
array([[-0.01805437+0.j, -0.01762974+0.j, -0.01781023+0.j],
       [ 0.00602559-0.j,  0.00603163+0.j,  0.00609548+0.j],
       [ 0.00463317-0.j,  0.00463941+0.j,  0.00470356+0.j]])

when you really wanted

>>> a.dot(v[:,0])
array([-0.05349434,  0.0181527 ,  0.01397614])
>>> w[0] * b.dot(v[:,0])
array([-0.05349434+0.j,  0.01815270+0.j,  0.01397614+0.j])

or

>>> matrix(a)*matrix(v[:,0]).T
matrix([[-0.05349434],
        [ 0.0181527 ],
        [ 0.01397614]])
>>> w[0]*matrix(b)*matrix(v[:,0]).T
matrix([[-0.05349434+0.j],
        [ 0.01815270+0.j],
        [ 0.01397614+0.j]])