Search code examples
numpyeigenvalueeigenvector

Sanity Checking Output of Python numpy eig() function

I have a question on numpy.linalg.eig().

Here's my covariance matrix after normalizing /standardizing the data.

lr_cov = np.cov(lr_norm, rowvar = False, ddof = 0)
lr_cov

array([[ 0.95454545,  0.88156287,  0.8601369 ],
       [ 0.88156287,  0.95454545,  0.87367031],
       [ 0.8601369 ,  0.87367031,  0.95454545]])

I use the eig() function as below -- no problems here.

eig_val, eig_vec = np.linalg.eig(lr_cov)

eig_vec

array([[-0.57694452, -0.6184592 ,  0.53351967],
       [-0.57990975, -0.14982268, -0.80078577],
       [-0.57518668,  0.77140222,  0.27221115]])

eig_val

array([ 2.69815538,  0.09525935,  0.07022164])

But when I proceed to sanity check that (Covariance Matrix)*(Eigen vector) = (Eigen Value)*(Eigen Vector), the LHS and RHS in this case don't match up.

lr_cov*eig_vec

array([[-0.55071977, -0.54521067,  0.45889996],
       [-0.5112269 , -0.14301256, -0.69962276],
       [-0.49473928,  0.67395122,  0.25983791]])

eig_val*eig_vec

array([[-1.55668595, -0.05891402,  0.03746463],
       [-1.5646866 , -0.01427201, -0.05623249],
       [-1.55194302,  0.07348327,  0.01911511]])

What I am doing incorrectly?

Solution

  • Two points:

    • * is element-wise multipication. Use the dot() method for matrix multiplication.
    • eig_val is a 1d array. Convert it to a 2d square diagonal array with np.diag(eig_val).

    Example:

    In [70]: cov
Out[70]: 
array([[ 0.95454545,  0.88156287,  0.8601369 ],
       [ 0.88156287,  0.95454545,  0.87367031],
       [ 0.8601369 ,  0.87367031,  0.95454545]])

In [71]: eig_val, eig_vec = np.linalg.eig(cov)

In [72]: cov.dot(eig_vec)
Out[72]: 
array([[-1.55668595, -0.05891401,  0.03746463],
       [-1.56468659, -0.01427202, -0.05623249],
       [-1.55194302,  0.07348327,  0.01911511]])

In [73]: eig_vec.dot(np.diag(eig_val))
Out[73]: 
array([[-1.55668595, -0.05891401,  0.03746463],
       [-1.56468659, -0.01427202, -0.05623249],
       [-1.55194302,  0.07348327,  0.01911511]])

    In the last line, np.diag(eig_val) is on the right in order to multiply each column of eig_vec by the corresponding eigenvalue.

    If you take advantage of numpy's broadcasting, you don't have to use np.diag(eig_val), and you can use element-wise multiplication (in either order, since element-wise multiplication is commutative):

    In [75]: eig_vec * eig_val  # element-wise multiplication with broadcasting

Out[75]: 
array([[-1.55668595, -0.05891401,  0.03746463],
       [-1.56468659, -0.01427202, -0.05623249],
       [-1.55194302,  0.07348327,  0.01911511]])

In [76]: eig_val * eig_vec
Out[76]: 
array([[-1.55668595, -0.05891401,  0.03746463],
       [-1.56468659, -0.01427202, -0.05623249],
       [-1.55194302,  0.07348327,  0.01911511]])