# Eigenvalues and Eigenvectors

import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg as la
%matplotlib inline


## Definition

Let $A$ be a square matrix. A non-zero vector $\mathbf{v}$ is an eigenvector for $A$ with eigenvalue $\lambda$ if

$$A\mathbf{v} = \lambda \mathbf{v}$$

Rearranging the equation, we see that $\mathbf{v}$ is a solution of the homogeneous system of equations

$$\left( A - \lambda I \right) \mathbf{v} = \mathbf{0}$$

where $I$ is the identity matrix of size $n$. Non-trivial solutions exist only if the matrix $A - \lambda I$ is singular which means $\mathrm{det}(A - \lambda I) = 0$. Therefore eigenvalues of $A$ are roots of the characteristic polynomial

$$p(\lambda) = \mathrm{det}(A - \lambda I)$$

## scipy.linalg.eig

The function scipy.linalg.eig computes eigenvalues and eigenvectors of a square matrix $A$.

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 $A$, and eigvecs is a 2D NumPy array with the corresponding eigenvectors in the columns:

results = la.eig(A)


The eigenvalues of $A$ 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 $A$ 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)

[[9 4 7 9]
[4 3 9 5]
[7 3 6 8]
[0 5 1 1]]


Create the symmetric matrix $S = P + P^T$:

S = P @ P.T
print(S)

[[227 156 189  36]
[156 131 131  29]
[189 131 158  29]
[ 36  29  29  27]]


Let's unpack the eigenvalues and eigenvectors of $S$:

evals, evecs = la.eig(S)
print(evals)

[5.03989629e+02+0.j 3.13953244e-01+0.j 1.61350731e+01+0.j
2.25613442e+01+0.j]


The eigenvalues all have zero imaginary part and so they are indeed real numbers:

evals = evals.real
print(evals)

[5.03989629e+02 3.13953244e-01 1.61350731e+01 2.25613442e+01]


The corresponding eigenvectors of $A$ are:

print(evecs)

[[ 0.66646882  0.62743094  0.31665684 -0.24875323]
[ 0.48300604  0.03934429 -0.77317651  0.4090908 ]
[ 0.55645551 -0.77634479  0.18643326 -0.2299754 ]
[ 0.11349778 -0.04551065  0.5168841   0.84727673]]


Let's check that the eigenvectors are orthogonal to each other:

v1 = evecs[:,0] # First column is the first eigenvector
print(v1)

[0.66646882 0.48300604 0.55645551 0.11349778]

v2 = evecs[:,1] # Second column is the second eigenvector
print(v2)

[ 0.62743094  0.03934429 -0.77634479 -0.04551065]

v1 @ v2

-2.7755575615628914e-16


The dot product of eigenvectors $\mathbf{v}_1$ and $\mathbf{v}_2$ 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 $M$ is diagonalizable if it is similar to a diagonal matrix. In other words, $M$ is diagonalizable if there exists an invertible matrix $P$ such that $D = P^{-1}MP$ is a diagonal matrix.

A beautiful result in linear algebra is that a square matrix $M$ of size $n$ is diagonalizable if and only if $M$ has $n$ independent eigevectors. Furthermore, $M = PDP^{-1}$ where the columns of $P$ are the eigenvectors of $M$ and $D$ has corresponding eigenvalues along the diagonal.

Let's use this to construct a matrix with given eigenvalues $\lambda_1 = 3, \lambda_2 = 1$, and eigenvectors $v_1 = [1,1]^T, v_2 = [1,-1]^T$.

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 $M$ 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 $M$ be a square matrix. Computing powers of $M$ by matrix multiplication

$$M^k = \underbrace{M M \cdots M}_k$$

is computationally expensive. Instead, let's use diagonalization to compute $M^k$ more efficiently

$$M^k = \left( P D P^{-1} \right)^k = \underbrace{P D P^{-1} P D P^{-1} \cdots P D P^{-1}}_k = P D^k P^{-1}$$

Let's compute $M^20$ both ways and compare execution time.

Pinv = la.inv(P)

from numpy.linalg import matrix_power as mpow

k = 20

%%timeit
mpow(M,k)

12.9 µs ± 1.65 µs per loop (mean ± std. dev. of 7 runs, 100000 loops each)


Let's use diagonalization to do the same computation.

%%timeit
P @ D**k @ Pinv

6.94 µs ± 495 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)


Using diagonalization computes $M^{20}$ about twice as a fast!

## Exercises

Under construction