Eigenvalues and Eigenvectors
import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg as la
Definition
Let be a square matrix. A non-zero vector is an eigenvector for with eigenvalue if
Rearranging the equation, we see that is a solution of the homogeneous system of equations
where is the identity matrix of size . Non-trivial solutions exist only if the matrix is singular which means . Therefore eigenvalues of are roots of the characteristic polynomial
scipy.linalg.eig
The function scipy.linalg.eig
computes eigenvalues and eigenvectors of a square matrix .
Let's consider a simple example with a diagonal matrix:
A = np.array([[1,0],[0,-2]])
print(A)
[[ 1 0]
[ 0 -2]]
The function la.eig
returns a tuple (eigvals,eigvecs)
where eigvals
is a 1D NumPy array of complex numbers giving the eigenvalues of , and eigvecs
is a 2D NumPy array with the corresponding eigenvectors in the columns:
results = la.eig(A)
The eigenvalues of are:
print(results[0])
[ 1.+0.j -2.+0.j]
The corresponding eigenvectors are:
print(results[1])
[[1. 0.]
[0. 1.]]
We can unpack the tuple:
eigvals, eigvecs = la.eig(A)
print(eigvals)
[ 1.+0.j -2.+0.j]
print(eigvecs)
[[1. 0.]
[0. 1.]]
If we know that the eigenvalues are real numbers (ie. if is symmetric), then we can use the NumPy array method .real
to convert the array of eigenvalues to real numbers:
eigvals = eigvals.real
print(eigvals)
[ 1. -2.]
Notice that the position of an eigenvalue in the array eigvals
correspond to the column in eigvecs
with its eigenvector:
lambda1 = eigvals[1]
print(lambda1)
-2.0
v1 = eigvecs[:,1].reshape(2,1)
print(v1)
[[0.]
[1.]]
A @ v1
array([[ 0.],
[-2.]])
lambda1 * v1
array([[-0.],
[-2.]])
Examples
Symmetric Matrices
The eigenvalues of a symmetric matrix are always real and the eigenvectors are always orthogonal! Let's verify these facts with some random matrices:
n = 4
P = np.random.randint(0,10,(n,n))
print(P)
[[7 0 6 2]
[9 5 1 3]
[0 2 2 5]
[6 8 8 6]]
Create the symmetric matrix :
S = P @ P.T
print(S)
[[ 89 75 22 102]
[ 75 116 27 120]
[ 22 27 33 62]
[102 120 62 200]]
Let's unpack the eigenvalues and eigenvectors of :
evals, evecs = la.eig(S)
print(evals)
[361.75382302+0.j 42.74593101+0.j 26.33718907+0.j 7.16305691+0.j]
The eigenvalues all have zero imaginary part and so they are indeed real numbers:
evals = evals.real
print(evals)
[361.75382302 42.74593101 26.33718907 7.16305691]
The corresponding eigenvectors of are:
print(evecs)
[[-0.42552429 -0.42476765 0.76464379 -0.23199439]
[-0.50507589 -0.54267519 -0.64193252 -0.19576676]
[-0.20612674 0.54869183 -0.05515612 -0.80833585]
[-0.72203822 0.4733005 0.01415338 0.50442752]]
Let's check that the eigenvectors are orthogonal to each other:
v1 = evecs[:,0] # First column is the first eigenvector
print(v1)
[-0.42552429 -0.50507589 -0.20612674 -0.72203822]
v2 = evecs[:,1] # Second column is the second eigenvector
print(v2)
[-0.42476765 -0.54267519 0.54869183 0.4733005 ]
v1 @ v2
-1.1102230246251565e-16
The dot product of eigenvectors and is zero (the number above is very close to zero and is due to rounding errors in the computations) and so they are orthogonal!
Diagonalization
A square matrix is diagonalizable if it is similar to a diagonal matrix. In other words, is diagonalizable if there exists an invertible matrix such that is a diagonal matrix.
A beautiful result in linear algebra is that a square matrix of size is diagonalizable if and only if has independent eigevectors. Furthermore, where the columns of are the eigenvectors of and has corresponding eigenvalues along the diagonal.
Let's use this to construct a matrix with given eigenvalues , and eigenvectors .
P = np.array([[1,1],[1,-1]])
print(P)
[[ 1 1]
[ 1 -1]]
D = np.diag((3,1))
print(D)
[[3 0]
[0 1]]
M = P @ D @ la.inv(P)
print(M)
[[2. 1.]
[1. 2.]]
Let's verify that the eigenvalues of are 3 and 1:
evals, evecs = la.eig(M)
print(evals)
[3.+0.j 1.+0.j]
Verify the eigenvectors:
print(evecs)
[[ 0.70710678 -0.70710678]
[ 0.70710678 0.70710678]]
Matrix Powers
Let be a square matrix. Computing powers of by matrix multiplication
is computationally expensive. Instead, let's use diagonalization to compute more efficiently
Let's compute both ways and compare execution time.
Pinv = la.inv(P)
k = 20
%%timeit
result = M.copy()
for _ in range(1,k):
result = result @ M
42.1 µs ± 11.4 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
Let's use diagonalization to do the same computation.
%%timeit
P @ D**k @ Pinv
6.42 µs ± 1.36 µs per loop (mean ± std. dev. of 7 runs, 100000 loops each)
Diagonalization computes much faster!
Exercises
Under construction