Python eig for generalized eigenvalue does not return valid eigenvectors

Question

An attempt to use scipy linalg.eig to solve a generalized eigenvalue problem. Then I check the solution I get, and it does not seem to have returned my own eigenvectors. In addition, the documentation assumes that the returned vectors are normalized, and this is not the case (although this does not bother me much).

Here is an example matrix:

>>> 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 the two would be equal, but they did not ... I also tried using eigh instead without success. I would be grateful for any help, I obviously missed something.

Answer 1

You can see what happens by looking at the shape of your output. Your a*v[:,0] should give a vector, so why do you get a 3x3 array? Answer: because you are not doing matrix multiplication, you are doing array multiplication by components.

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 did you really want

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