import numpy as np
from numpy import linalg as la
import pandas as pd
import matplotlib.pyplot as plt
import requests
# Read in Visualization code from Github (requires bokeh module)
exec(requests.get('https://raw.githubusercontent.com/edunford/ppol564/master/lectures/visualization_library/visualize.py').content)
vla = LinearAlgebra # assign class to an simplier naming convention.
Say we have some transformation matrix $\textbf{A}$.
A = np.array([[1,2],
[2,1]])
A
Let's transform some vector $\vec{x}$ by $\textbf{A}$.
x = np.array([-.3,.1])
print('Original vector x')
print(x)
# What happens to as we transform x by matrix A?
print("\nAx")
print(A.dot(x))
print("\nAAx")
print(A.dot(A.dot(x)))
print("\nAAAx")
print(A.dot(A.dot(A.dot(x))))
Let's do the same thing again with vector $\vec{v}$
v = np.array([.7,.7])
print('Original vector v')
print(v)
# What happens to as we transform v by matrix A?
print("\nAv")
print(A.dot(v))
print("\nAAv")
print(A.dot(A.dot(v)))
print("\nAAAv")
print(A.dot(A.dot(A.dot(v))))
Let's visualize this...
plot = vla()
plot.graph(extent=5)
plot.vector(x,add_color="black")
plot.vector(v,add_color="blue")
plot.show()
# Plot vector x and all its transformations
plot.vector(A.dot(x),add_color="black",alpha=.2)
plot.vector(A.dot(A.dot(x)),add_color="black",alpha=.2)
plot.vector(A.dot(A.dot(A.dot(x))),add_color="black",alpha=.2)
# Plot vector v and all its transformations
plot.vector(A.dot(v),add_color="blue",alpha=.2)
plot.vector(A.dot(A.dot(v)),add_color="blue",alpha=.2)
plot.vector(A.dot(A.dot(A.dot(v))),add_color="blue",alpha=.2)
plot.show()
What is going on here?
$\vec{v}$ isn't being transformed like $\vec{x}$ is each time; rather, it's being scaled!
print(v)
print(v*3)
print(v*(3**2))
print(v*(3**3))
Why?
$\vec{v}$ is an Eigenvector!
An Eigen vector is a vector that only expands or shrinks when undergoing a linear transformation. An eigenvector (or the characteristic vector of a linear transformation) is a non-zero vector that changes by only a scalar factor when that linear transformation is applied to it.
The factor by which that Eigenvector shrinks or expands is called its Eigenvalue, which we denote with $\lambda$.
The equation for locating an Eigenvalue and Eigenvector is as follows.
\begin{equation} \textbf{A}\vec{v} = \lambda \vec{v} \end{equation}
We want to locate some $\vec{v}$ where the transformation of that vector simply yields a scaled version of that original vector.
Think about this: our transformation matrix $\textbf{A}$ collapses to a single scalar value when applied on $\vec{v}$
First, let's do a little re-arranging
We're going from matrix multiplication to scalar multiplication, which is odd. To get around this, we can turn a scalar into a matrix by taking the scalar product of the identity matrix.
$$ \lambda \textbf{I}$$
$$ \lambda \begin{bmatrix}1 & 0 \\ 0 & 1 \end{bmatrix}$$
$$ \begin{bmatrix}\lambda & 0 \\ 0 & \lambda \end{bmatrix}$$
Also, $\vec{v}$ is on both sides of the equation, so let's re-arrange things.
$$ \textbf{A}\vec{v} = \lambda \textbf{I} \vec{v} $$
$$ \textbf{A}\vec{v} - \lambda \textbf{I} \vec{v} = 0$$
$$ (\textbf{A} - \lambda \textbf{I})\vec{v} = 0$$
Assuming that $\textbf{A}$ is a $2 \times 2$ matrix
$$ ( \begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix} - \begin{bmatrix}\lambda & 0 \\ 0 & \lambda \end{bmatrix})\vec{v} = 0$$
$$ \begin{bmatrix} a_{11} - \lambda & a_{12} \\ a_{21} & a_{22} - \lambda \end{bmatrix}\vec{v} = 0$$
Recall that what we are aiming for is locating an instance when $ \begin{bmatrix} a_{11} - \lambda & a_{12} \\ a_{21} & a_{22} - \lambda \end{bmatrix}$ collapses down into a singular value.
Put differently, we are looking for values of $\lambda$ that result in the determinant of our A matrix equaling zero.
$$ det(\textbf{A} - \lambda \textbf{I}) = 0$$
$$det(\begin{bmatrix} a_{11} - \lambda & a_{12} \\ a_{21} & a_{22} - \lambda \end{bmatrix}) = 0$$
So, we merely need to plug in values for $\lambda$ where this condition holds.
Let's try doing this computationally by just brute force plugging in values and having the computer do the work for us.
def find_lambdas(A,start,end):
'''
Calculates the determinant across a specified range of lamba entries.
'''
store = []
for lamb in np.arange(start,end+1,.5):
det = la.det(A - (lamb*np.eye(2)))
store.append([lamb,det])
return np.array(store)
# Scan between -5 and 5
vals = find_lambdas(A,-5,5)
# Let's plot this
plt.figure()
plt.ylabel('Determinant')
plt.xlabel('Lambda')
plt.axhline(0,color="black")
plt.plot(vals[:,0],vals[:,1],
linewidth=3,color="orange",alpha=.5)
# Plot instances when det == 0
plt.scatter(vals[np.where(vals[:,1]==0),0],
vals[np.where(vals[:,1]==0),1],
color="black")
plt.show()
# What are these values where the det == 0??
vals[np.where(vals[:,1]==0),0]
Naturally, we could also derive these values algebraically.
Recall the composition of matrix A.
\begin{equation} A = \begin{bmatrix} 1 & 2 \\ 2 & 1 \end{bmatrix} \end{equation}
\begin{equation} det(\begin{bmatrix} 1 - \lambda & 2 \\ 2 & 1 - \lambda \end{bmatrix}) = 0 \end{equation}
\begin{equation} (1 - \lambda)(1 - \lambda) - (2)(2)= 0 \end{equation}
\begin{equation} \lambda^2 - 2\lambda - 3 = 0 \end{equation}
The above equation is known as the "characteristic polynomial".
\begin{equation} (\lambda + 1)(\lambda - 3) = 0 \end{equation}
\begin{equation} \lambda = - 1 \\ \lambda = 3 \end{equation}
Let's now find the eigenvector that corresponds with each eigenvalue. We can do this easily by
plugging in our respective eigenvalues, and
Solving for the values of $\vec{v}$ in equation 1
\begin{equation} \begin{bmatrix} 1 - \lambda & 2 \\ 2 & 1 - \lambda \end{bmatrix} \begin{bmatrix}v_1\\v_2\end{bmatrix} = \begin{bmatrix}0\\0\end{bmatrix} \end{equation}
Solving for the eigenvector when $\lambda = 3$, using reduced row echelon form to solve the system of equations. The aim is to use row-wise manipulations to convert the matrix to an identity matrix (see lecture 14). Since the our matrix isn't full rank (by definition of what the eigen values are doing), we'll only be able to recover one of the unit vectors. In other words, our aim is to convert first vector column into the $hat{i}$ unit vector.
\begin{equation} \left| \begin{array}{cc|c} 1 - (3) & 2 & 0 \\ 2 & 1 -(3) & 0 \\ \end{array} \right| \end{equation}
\begin{equation} \left| \begin{array}{cc|c} -2 & 2 & 0 \\ 2 & -2 & 0 \\ \end{array} \right| \end{equation}
Add the first row to the second.
\begin{equation} \left| \begin{array}{cc|c} -2 & 2 & 0 \\ 0 & 0 & 0 \\ \end{array} \right| \end{equation}
Then divide the first row by -2.
\begin{equation} \left| \begin{array}{cc|c} 1 & -1 & 0 \\ 0 & 0 & 0 \\ \end{array} \right| \end{equation}
Let's re-represent this as the original matrix format
\begin{equation} \begin{bmatrix} 1 & -1 \\ 0 & 0 \end{bmatrix}\begin{bmatrix} v_1 \\ v_2 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix} \end{equation}
Multiple $\vec{v}$ by the reduced form equation.
\begin{equation} v_1(1) - v_2(1) = 0 \end{equation}
\begin{equation} v_1 = v_2 \\ v_1 = 1 \\ v_2 = 1 \end{equation}
Thus, the eigenvector $\vec{e}_{\lambda = 3}$ is all vectors that span the vector $\begin{bmatrix} 1 \\1 \end{bmatrix}$
Now, Solving for the eigenvector when $\lambda = -1$
\begin{equation} \left| \begin{array}{cc|c} 1 - (-1) & 2 & 0 \\ 2 & 1 -(-1) & 0 \\ \end{array} \right| \end{equation}
\begin{equation} \left| \begin{array}{cc|c} 2 & 2 & 0 \\ 2 & 2 & 0 \\ \end{array} \right| \end{equation}
Let's subtract the first row from the second and divide the first row by 2.
\begin{equation} \left| \begin{array}{cc|c} 1 & 1 & 0 \\ 0 & 0 & 0 \\ \end{array} \right| \end{equation}
\begin{equation} \left| \begin{array}{cc|c} 1 & 1 & 0 \\ 0 & 0 & 0 \\ \end{array} \right| \end{equation}
\begin{equation} \begin{bmatrix} 1 & 1 \\ 0 & 0 \end{bmatrix}\begin{bmatrix} v_1 \\ v_2 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix} \end{equation}
\begin{equation} v_1 + v_2 = 0 \\ v_1 = -v_2 \\ v_1 = -1 \\ v_2 = 1 \end{equation}
Thus, the eigenvector $\vec{e}_{\lambda = -1}$ is all vectors that span the vector $\begin{bmatrix} -1 \\1 \end{bmatrix}$
Let's visualize these two vectors to see how we did.
e1 = np.array([1,1])
e2 = np.array([1,-1])
# Generate Graphic
plot.clear().graph(20)
plot.vector(A.dot(e1),add_color="blue",alpha=.2)
plot.vector(A.dot(A.dot(e1)),add_color="blue",alpha=.2)
plot.vector(A.dot(A.dot(A.dot(e1))),add_color="blue",alpha=.2)
plot.vector(A.dot(e2),add_color="red",alpha=.2)
plot.vector(A.dot(A.dot(e2)),add_color="red",alpha=.2)
plot.vector(A.dot(A.dot(A.dot(e2))),add_color="red",alpha=.2)
plot.show()
What sorts of conclusions can we draw from these two eigenvectors and their corresponding eigenvalues?
With each transformation, eigenvector that corresponds with the eigen value 3 grows, whereas the eigenvector that corresponds with the eigen value -1 is stable (it just flips signs but never really changes.
This information is really useful when studying a system of linear equations to know if some process is shrinking, growing, or at equilibrium (see the reading for examples!)
numpy
¶evals, evecs = la.eig(A)
evals
evecs
Why do these eigenvectors differ slightly from the ones we located manually?
e1 = np.array([1,1])
e2 = np.array([-1,1])
print("\nOur Eigenvectors normalized\n")
print('E1 =',e1/la.norm(e1))
print('E2 =',e2/la.norm(e2))
Building off the above where we derived each eigenvalue and eigenvector individually, we can express some matrix $\textbf{A}$ as three matrices. This process of breaking down a matrix into the product of other matrices is known as factorization or a decomposition.
Matrix decompositions come in many flavors in linear algebra and you might run into different ones to help resolve different types of problems.
If $\textbf{A}$ is an $n \times n$ matrix (with $n$ linearly independent eigenvectors), then we can express the decomposition as follows.
\begin{equation} \textbf{A}\textbf{V} = \textbf{V}\Lambda \end{equation}
\begin{equation} \textbf{A} = \textbf{V}\Lambda\textbf{V}^{-1} \end{equation}
where
L = np.diag(evals)
V = evecs
V_inv = la.inv(V)
print('Original Matrix A\n')
print(A)
print('\nDecomposed Parts\n')
print('Eigenvectors V\n')
print(V.round(2))
print('\nEigenvalues L\n')
print(L.round(2))
print('\nEigenvectors Inverse V^-1\n')
print(V_inv.round(2))
print('\nRecomposed Matrix VLV^-1\n')
print(V.dot(L).dot(V_inv))
Eigen decompositions have some very interesting properties, one of the most interesting features is that if we want to apply a transformation multiple times, we merely need to take a polynomial of our eigenvalues! In way, we already saw this in the initial example above.
\begin{equation} \textbf{A} = \textbf{V}\Lambda\textbf{V}^{-1} \end{equation}
\begin{equation} \textbf{A}\textbf{A} = \textbf{A}^2 = \textbf{V}\Lambda(\textbf{V}^{-1} \textbf{V})\Lambda\textbf{V}^{-1} = \textbf{V}\Lambda^2\textbf{V}^{-1} \end{equation}
\begin{equation} \textbf{A}^n = \textbf{V}\Lambda^n\textbf{V}^{-1} \end{equation}
This can greatly ease some computations, such as Markov chains!
A.dot(A.dot(A.dot(A)))
V.dot(L**4).dot(V_inv)
So we needed to transform $\vec{x}$ by $\textbf{A}$ 25 times we could just do this.
A_25 = V.dot(L**25).dot(V_inv)
A_25.dot(x)
la.det(A)
evals.prod()
A.trace()
sum(evals)
evals
la.eigvals(la.inv(A))
This means each eigenvector describes a unique dimension!
evecs[:,0].dot(evecs[:,1])