CODEWITHSUNDEEP
pythonnumpy
rocksNwaves
rocksNwaves

Reputation: 6163

numpy always gets complex eigenvalues, and wrong eigenvectors

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.

Upvotes: 3

Views: 3617

Answers (2)

Richard Nemeth
Richard Nemeth

Reputation: 1874

What you notice here is the relative approximation error.

Numpy is not computing the eigenvalues and eigenvectors exactly, instead it is utilising powerful and optimised numerical algorithms that result not in the exact answers, but in fairly close ones instead. For these ones, it is utilising an older LAPACK library dggev: http://www.netlib.org/lapack/double/dggev.f

If you are looking for an exact answer, then Numpy is not a solution for you. Instead have a look at Sympy. A similar problem is also being discussed here: how are numpy eigenvalues and eigenvectors computed

Upvotes: 1

salt-die
salt-die

Reputation: 854

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.

Upvotes: 2

Related Questions