Is scipy.linalg.eig giving the correct left eigenvectors?

Question

I have a question regarding the way how scipy.linalg.eig computes left and right eigenvectors. Maybe I misunderstood everything, but things seem not to be right to me...

From the beginning. To get eigenvalues and both eigenvectors I used the following:

ev, left_v, right_v = scipy.linalg.eig(A, left=True)

According to the manual, after setting left=True while calling the function I should expect to get left eigenvectors as columns of left_v where the ith column refers to the ith eigenvalue. However, the results were not what i anticipated so I did a simple check.

I computed right and left eigenvectors invoking the function twice (look here for details):

right_ev, right_v_2 = scipy.linalg.eig(A)
left_ev, left_v_2 = scipy.linalg.eig(A.T)

where columns ofleft_v_2 are eigenvectors associated with corresponding values in left_ev. Worth stressing that both right_ev_2 and left_ev_2 give the same eigenvalues, however they are in different order, which needs to be accounted for.

Comparing left_ev and left_ev_2 (after reordering with respect to eigenvalues) one can quickly spot that the former is the conjugate of the latter and therefore left_ev obtained from scipy.linalg.eig with left=True is not a valid left eigenvector.

Another check on the validity of the eigenvectors can be done based on the fact that for an arbitrary real square matrix left and right eigenvectors are biorthogonal, i.e.:

left_v.T.dot(right_v) should give a diagonal matrix, but it doesn't, until i change it to: left_v.T.conj().dot(right_v),

while:

left_v_2.T.dot(right_v_2) gives an anticipated diagonal matrix.

Did anyone encounter similar problem before? Am I right with what I say? Is the sciPy manual a bit imprecise while describing eig? Can you give any advice?

Many thanks!

Answer 1

About vl, the eig docstring says:

a.H vl[:,i] = w[i].conj() b.H vl[:,i]

Or, taking the conjugate transpose (i.e. Hermitian transpose) of both sides (which is what .H means), and assuming b is the identity,

vl[:,i].H a = w[i] vl[:,i].H

So the rows of the conjugate transpose of vl are the actual left eigenvectors of a.

Numpy arrays don't actually have the .H attribute, so you must use .conj().T.

Here's a script to verify the calculation:

import numpy as np
from scipy.linalg import eig

# This only affects the printed output.
np.set_printoptions(precision=4)

a = np.array([[6, 2],
              [-1, 4]])

w, vl, vr = eig(a, left=True)

print "eigenvalues:", w
print

# check the left eigenvectors one-by-one:
for k in range(a.shape[0]):
    val = w[k]
    # Use a slice to maintain shape; vec is a 2x1 array.
    # That allows a meaningful transpose using .T.
    vec = vl[:, k:k+1]
    # rowvec is 1x2; it is the conjugate transpose of vec.
    # This should be the left eigenvector.
    rowvec = vec.conj().T
    # Verify that rowvec is a left eigenvector
    lhs = rowvec.dot(a)
    rhs = val * rowvec
    print "Compare", lhs, "to", rhs
    print rowvec, "is",
    if not np.allclose(lhs, rhs):
        print "*NOT*",
    print "a left eigenvector for eigenvalue", val

print
print "Matrix version:"
print "This"
print vl.conj().T.dot(a)
print "should equal this"
print np.diag(w).dot(vl.conj().T)

Output:

eigenvalues: [ 5.+1.j  5.-1.j]

Compare [[ 1.6330+2.4495j  4.0825+0.8165j]] to [[ 1.6330+2.4495j  4.0825+0.8165j]]
[[ 0.4082+0.4082j  0.8165-0.j    ]] is a left eigenvector for eigenvalue (5+1j)
Compare [[ 1.6330-2.4495j  4.0825-0.8165j]] to [[ 1.6330-2.4495j  4.0825-0.8165j]]
[[ 0.4082-0.4082j  0.8165+0.j    ]] is a left eigenvector for eigenvalue (5-1j)

Matrix version:
This
[[ 1.6330+2.4495j  4.0825+0.8165j]
 [ 1.6330-2.4495j  4.0825-0.8165j]]
should equal this
[[ 1.6330+2.4495j  4.0825+0.8165j]
 [ 1.6330-2.4495j  4.0825-0.8165j]]

Now, the eig docstring also says in the description of the return values:

vl : double or complex ndarray
    The normalized left eigenvector corresponding to the eigenvalue
    ``w[i]`` is the column v[:,i]. Only returned if ``left=True``.
    Of shape ``(M, M)``.

and that is potentially misleading, since the conventional definition of a left eigenvector (e.g. http://mathworld.wolfram.com/LeftEigenvector.html or http://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors#Left_and_right_eigenvectors) is a row vector, so it is the conjugate transpose of the column of vl that is actually the left eigenvector.