numpy always gets complex eigenvalues, and wrong eigenvectors

Question

I'm working on doing simple linear algebra manipulations with numpy. Everything has been really great until now, when I take simple 2x2 matrices whose eigenvalues and vectors I know, and test numpy on them.

For example the example matrix below, there is a single eigenvalue e=1, and a single associated eigenvector [-3, 1]:

A = np.array([[-2, -9],
              [ 1,  4]])


vals, vects = np.linalg.eig(A)
vals2, vects2 = np.linalg.eigh(A)
vals3 = np.linalg.eigvals(A)

print("Eigenvalues (linalg.eig): \n", vals)
print("Eigenvalues (linalg.eigvals): \n", vals3)
print("Eigenvectors: \n", vects)

results in the following :

Eigenvalues (linalg.eig):
 [1.+1.89953279e-08j 1.-1.89953279e-08j]

Eigenvalues (linalg.eigvals):
  [1.+1.89953279e-08j 1.-1.89953279e-08j]

Eigenvectors:
 [[ 0.9486833 +0.00000000e+00j  0.9486833 -0.00000000e+00j]
 [-0.31622777-2.00228337e-09j -0.31622777+2.00228337e-09j]] 

I realize the eigenvectors are in column format. If you neglect the small imaginary parts, both vectors are ALMOST scalar multiples of the single correct eigenvector.

My matrix is not symmetric or conjugate symmetric, and therefore linalg.eigh should not work, but I tried it anyway. It gives me real values, but they are totally bogus.

I have also searched other similar questions, and have found no satisfying answers. The closest I have come is one answer that suggest writing a function to strip the eigenvalues of their small imaginary parts. This works, but still leaves incorrect eigenvectors.

What's going on here? How can I fix it simply? It seems a bit much to write my own functions to correct such a well established library in order to make it get simple calculations correct.

Answer 1

Numpy is returning a normalized eigenvector for each eigenvalue; as the eigenvalue here has multiplicity 2, it returns the same eigenvector twice. You can use np.real_if_close to reconvert to real:

In [156]: A = np.array([[-2, -9],
     ...:               [ 1,  4]])
     ...: 
     ...: w, v = np.linalg.eig(A)
     ...: v = np.real_if_close(v, tol=1)
     ...: v
Out[156]: 
array([[ 0.9486833 ,  0.9486833 ],
       [-0.31622777, -0.31622777]])

Sanity check:

In [157]: v * 10**.5
Out[157]: 
array([[ 3.,  3.],
       [-1., -1.]])

If you'd like exact solutions instead of these prone to machine error, may I suggest sympy:

In [159]: from sympy import Matrix
     ...: m = Matrix(A)
     ...: m.eigenvects()
Out[159]: 
[(1,
  2,
  [Matrix([
   [-3],
   [ 1]])])]

Where we get the eigenvalue 1, its multiplicity 2, and the eigenvector associated with it.