Precept 7: Eigenvalue iterations¶

Programming tasks¶

This programming task is optional! See the PDF for proof based questions.

Task 1¶

In this precept, we will investigate power method and inverse power method for symmetric matrices.

We start by creating an example matrix. We $1000$ uniformly random points on the unit sphere by generating Gaussian points using np.random.randn and dividing by their length. Then we make a scatter plot of the points.

Your task is to compute the $n \times n$ matrix $D$ whose $(i,j)$-th entry is the distance between $X[i,:]$ and $X[j,:]$, and define $$ A(i,j) = \exp(-D(i,j)^2). $$

In [ ]:

# Task 1
import numpy as np 
import matplotlib.pyplot as plt

# n random vectors on the sphere
n = 1000
X = np.random.randn(n,3)
X = X / np.linalg.norm(X, axis =1,keepdims=True)


# Scatter plot
fig = plt.figure(figsize = (10,10))
ax = fig.add_subplot(projection='3d')
ax.scatter(X[:,0], X[:,1], X[:,2])
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
plt.show()

## Compute the matrix A :
A = np.zeros((n,n))

Task 2¶

Using the built in numpy command compute the built in eigenvalues and eigenvectors of $A$. Sort the eigenvalues is descending order and sort the columns of $V$ in a corresponding way. Check $$ \text{norm}(A - V D V^\top )/\text{norm}(A) $$ and $$ \text{norm}(I - V^\top V) $$ where $I$ is the identity matrix, to make sure that you have decomposed $A$ into $V D V^\top$, where $V$ is orthogonal and $D$ is diagonal, which is the decomposition guaranteed by the spectral theorem.

In [ ]:

# Task 2 use the following function to compute eigenvalues and eigenvectors of A
help(np.linalg.eigh)
help(np.flip) # May want to use this function to reorder eigenvalues/vectors

Task 3¶

Visualize the eigenvectors in $V$. In particular, the eigenvectors $V[:,j]$ are vectors in $\mathbb{R}^n$. Make a scatter plot of $X$ and color the scatter plot using $V_2$, the second eigenvector (equal to $V[:,1]$ in numpy), the eigenvector associated with the second largest eigenvalue These eigenvectors are very rough approximations of spherical harmonics https://en.wikipedia.org/wiki/Spherical_harmonics

In [ ]:

# Task 3: scatter plot of eigenvectors
# Scatter plot

Task 4¶

First, use the power method to find $V_1$ (the first eigenvector, associated with the greatest eigenvalue). That is, start with a random vector $v$ and keep applying the matrix and normalizing until you start to converge to $v$. You can check the error by $$ \min(\text{norm}(v - V_1),\text{norm}(v + V_1)) $$ the reason that we need to check $-v$ and $v$ is that eigenvectors are only defined up to multiplication by a constant. Even when $v$ is normalized to have length $1$ there is still a sign ambiguity.

In [ ]:

# Task 4 implement power method

def distance_to_eigenvector(vk, eigenvec):
    return np.min([np.linalg.norm(vk - eigenvec), np.linalg.norm(vk + eigenvec)])


err = np.zeros(100)

## add your code
plt.semilogy(err)

Task 5¶

Next, use the power method to find $V_2$ (the second eigenvector) using the following trick: at each iteration, before multiplying by $A$, project the vector $v$ to the space orthogonal to $V_1$. Theoretically, it would suffice to just project $v$ orthogonal to the space $V_1$ initially, but to avoid any numerical issues perform the projection each iteration. Recall that the project orthogonal to a vector $u$ which is normalized to have unit length is $$ w = v - u(u^\top v) $$ Check that you are able to get $V(:,2)$.

In [ ]:

# Task 5 implement power method for v2

def distance_to_eigenvector(vk, eigenvec):
    return np.min([np.linalg.norm(vk - eigenvec), np.linalg.norm(vk + eigenvec)])


err = np.zeros(100)

## add your code
plt.semilogy(err)

Task 6¶

Next, using the eigenvalues $d$ that you computed using the built in solver, use the inverse power method to compute $V_1$ and $V_2$

In [ ]:

# Implement task 6

# add your code

print(distance_to_eigenvector(v1, V[:,0]))

# add your code

print(distance_to_eigenvector(v2, V[:,1]))