MAT 321 - PS4

• You can (and probably should) discuss assignments with each others, and the course staff, but you must write and understand any solutions/code that you submit. You can consult any resource you want for this homework, but you should cite sources you used.
• You must upload the assignment to Gradescope before the deadline. You should submit this notebook in .ipynb format, and optinally a single PDF file with your written answers.

Q1 Diagonalization of Circulant Matrices

In these problem sets, we have used repeatedly that circulant matrices are diagonalized by the DFT to:

• derive fast mat-vecs with circulant and Toeplitz matrices
• derive fast solvers for circulant and perturbed circulant systems
• deblur large datasets of images

That is, we used the fact that:

$$C = F^{-1} \Lambda F $$

where $F$ is the DFT matrix, and $\text{diag}(\Lambda) = Fc$. The goal of this exercise is to prove this identity.

Note that this is an eigendecomposition of $C$, where $V = F^{-1}$ is the matrix of eigenvectors. Recall that the DFT matrix is

$$ F = \frac{1}{\sqrt{n}} \begin{bmatrix} 1 & 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega & \omega^2 & \omega^3 & \cdots & \omega^{n-1} \\ 1 & \omega^2 & \omega^4 & \omega^6 & \cdots & \omega^{2(n-1)} \\ 1 & \omega^3 & \omega^6 & \omega^9 & \cdots & \omega^{2(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \ddots \\ 1 & \omega^{n-1} & \omega^{2(n-1)} & \omega^{3(n-1)} & \cdots & \omega^{(n-1)(n-1)} \\ \end{bmatrix} $$

where $\omega = \exp(- \imath 2 \pi / n) $, $\imath^2 = -1$ (note that $\omega$ depends on $n$). Also recall that $F$ with this normalization is unitary, that is $F^{-1} = F^* = \bar{F}^T$, and $F$ is symmetric (but not hermitian), so $F^{-1} = \bar{F}$.

We prove this identity in 3 steps:

1. Let $D_n \in \mathbb{R}^{n \times n}$ be the circulant downward shift matrix. For $n = 5$, this matrix is:

$$ D_5 = \begin{bmatrix} & & & & 1 \\ 1 & & & & \\ & 1& & & \\ & & 1& & \\ & & & 1& \\ \end{bmatrix} $$. We will first show that the DFT diagonalizes $D_n$: $D_n = F^{-1} \Gamma F $ for some diagonal matrix $\Gamma$. Here's how to proceed:

• Show that $D_n {F}_{i+1} = \gamma_{i+1} {F}_{i+1} $ where ${F}_{i+1}$ is the $i+1$ column of ${F}$:

$$ F_{i+1} = \frac{1}{\sqrt{n}} \begin{bmatrix} 1 \\ {\omega}^{i}\\ {\omega}^{2i}\\ {\omega}^{3i}\\ \vdots \\ {\omega}^{(n-1)i}\\ \end{bmatrix} $$

and $\gamma_{i+1} = \bar{\omega}^i$.

• Argue that this implies $D_n F = F \bar{\Gamma} $, where $\bar{\Gamma}$ is a diagonal matrix with diagonal entries $\gamma_i$.
• Take the complex conjugate of previous part, and conclude that $D_n = F^{-1} \Gamma F$

2. Write any circulant matrix as an $(n-1)$ degree polynomial in $D_n$. That is, find the coefficients $\alpha_k$ so that:

$$C = \sum_{k=0}^{n-1} \alpha_k D_n^k $$

Recall that $D_n^0$ is defined as the identity. You can do this by inspection.

3. Use the previous two parts to show $C = F^{-1} \Lambda F $ for some diagonal matrix $\Lambda$. Some inspection reveals that $\text{diag}(\Lambda) = Fc$, but you do not need to prove this.

Add your answer to Q1 here

Q2 The Bartels-Stewart algorithm

In this problem, we derive an algorithm for solving a special kind of matrix equation: Sylvester equations. A Sylvester equation is of the form: $$ AX - XB = F $$ where $A,X,B,F \in \mathbb{R}^{n \times n}$. Using the Kronecker formalism from HW3, we can solve this matrix equation in $\mathcal{O}(n^6)$, but we can do better: the Bartels-Stewart algorithm is an algorithm that solves this equation in $\mathcal{O}(n^3)$. The Bartels-Stewart algorithm is based on one clever observation: if $A$ and $B$ are diagonal matrices, then this equation is easy to solve. Say $A = \text{diag}(d)$, $B = \text{diag}(p)$, then the $(i,j)$-th entry of the Sylvester equation reads: $$ d_i x_{i,j} - x_{i,j} p_j = f_{i,j} $$ i.e., the system is uncoupled and $x_{i,j}$ can be easily found as: $$ x_{i,j} = \frac{f_{i,j}}{d_i - p_j} \ . $$ From this equation, it is clear that there is a unique solution if and only if $d_i - p_j \neq 0 \leftrightarrow d_i \neq p_j$ for all $i,j$.

To solve a general Sylvester equation, Bartels-Stewart reduces a general Sylvester equation to the diagonal case using an eigenvalue decomposition.

2.1

Assume $A$ and $B$ are diagonalizable, that their eigenvalues decompositions are $$ A = U \Sigma U^{-1}, \qquad B = V\Lambda V^{-1}, $$ and that $X$ satisfies the Sylvester equation $AX - XB = F$. Show that the matrix $\hat{X} = U^{-1} X V $ satisfies the following transformed (diagonal) Sylvester equation:

$$ \Sigma \hat{X}- \hat{X} \Lambda = \hat{F} $$

for some matrix $\hat{F} $.

2.2

Use the previous part to show that if $A$ is symmetric positive definite, and $B$ is symmetric negative definite (i.e. $C := -B$ is symmetric positive definitive), then the Sylvester equation has a unique solution.

2.3

Use 1.1 to design and implement an $\mathcal{O}(n^3)$ solver for a general Sylvester equation. You may assume that the matrices are diagonalizable but not SPD.

You may use np.linalg.eig() which computes an eigenvalue decomposition of $A$ in $\mathcal{O}(n^3)$.

Reminder: the number one sin of Numerical Linear Algebra is to form the inverse of a matrix. In particular, these two commands should be avoided like the plague: np.linalg.inv(A) or np.linalg.solve(A,I) (they do the same computation). Instead, you should use: np.linalg.solve(A,F) to compute $X = A^{-1} F$ and np.linalg.solve(A.T,F.T).T to compute $X = F A^{-1} = ({A^{-1}}^T {F}^T)^T$

Add your answer to parts 1 and 2 here

import numpy as np

def Bartels_Stewart(A,B,F):
    ## ADD YOUR CODE HERE
    X = np.zeros_like(F)
    return X


# THIS CODE IS WRITTEN FOR YOU. NO NEED TO CHANGE IT
n = 10
A = np.random.randn(n,n)
B = np.random.randn(n,n)
X = np.random.randn(n,n)

F = A@X - X@(B)

Xsolve = Bartels_Stewart(A, B, F)

# Should be close to machine epsilon ~1e-14 is fine
print('error: ', np.linalg.norm(Xsolve - X )/np.linalg.norm(X))

error:  1.0

Q3 Fast PDE Solver: Poisson's Equation

In HW2, we derived a numerical scheme to solve Poisson's equation in 1D:

$$ - \frac{d^2}{dx^2} u(x) = f(x), \quad u(0) = u(1) = 0. $$

The discretized system takes the form:

$$ T u = h^2 f, $$

where

$$ T = \begin{bmatrix} 2 & -1 & & & \\ -1 & 2 & -1 & & \\ & -1 & 2 & -1 & \\ & & \ddots & \ddots & \ddots \\ & & & -1 & 2 \\ \end{bmatrix}, \quad f = \begin{bmatrix} f(h) \\ f(2h) \\ \vdots \\ f(nh) \\ \end{bmatrix}, $$

and $ h = \frac{1}{n+1} $. Since $ T \in \mathbb{R}^{n \times n} $ is tridiagonal, we can solve this system in $ O(n) $ operations using a banded LU. All variations of the 1D problem are similarly straightforward.

The 2D case is more interesting. Consider Poisson's equation on the unit square with homogeneous Dirichlet boundary conditions:

$$ - \frac{d^2}{dx^2} u(x,y) - \frac{d^2}{dy^2} u(x,y) = f(x,y), \\ u(x,0) = u(x,1) = 0, \\ u(0,y) = u(1,y) = 0. $$

Poisson's in 2D and 3D equation models diffusion processes, with applications in fluid flow, heat transfer, and chemical transport. It also serves as a foundational problem for more complex equations.

We discretize the equation on an $ n \times n $ uniform grid $ (x_i, y_j) $, where $ x_i = ih $, $ y_j = jh $, and $ h = \frac{1}{n+1} $. Let $ F_{i,j} = h^2 f(x_i, y_j) $; that is, the entry $ (i,j) $ of matrix $ F \in \mathbb{R}^{n \times n} $ is $ h^2 $ times the value of $ f $ at $ (x_i, y_j) $. We seek approximate values $ U_{i,j} \approx u(x_i, y_j) $ using finite differences.

The operator $ u \rightarrow - \frac{d^2}{dx^2} u $ with boundary conditions $ u(x,0) = u(x,1) = 0 $ can be applied independently to the columns of $ U $, resulting in $ U \rightarrow TU $. Similarly, the operator $ u \rightarrow - \frac{d^2}{dy^2} u $ with boundary conditions $ u(0,y) = u(1,y) = 0 $ can be applied independently to the rows of $ U $, yielding $ U \rightarrow UT^T = UT $.

Therefore, the discretized operator becomes $ U \rightarrow TU + UT $, leading to the matrix equation:

$$ TU + UT = F. $$

Our goal is to solve this linear system for $ U $. Using the Kronecker formalism from the previous homework, we could solve this in $ O(n^6) $ operations—not ideal. The Bartels-Stewart algorithm improves this to $ O(n^3) $ complexity, but we can do even better. A key property of the second-order finite difference matrix $ T $ is that its eigenvalue decomposition can be computed analytically:

$$ T = S^{-1} \Lambda S, $$

where $ S $ is the Discrete Sine Transform (DST) matrix with components $ S_{kj} = \sqrt{\tfrac{2}{n+1}} \sin\left(\tfrac{kj\pi}{n+1}\right) $, and the eigenvalues are $ \lambda_j = 4 \sin^2 \left( \tfrac{j \pi}{2(n+1)} \right) $. The DST is the imaginary part of the Discrete Fourier Transform (DFT). Thanks to this connection, we can compute $ Sx $ in $ \mathcal{O}(n \log n) $ time using the Fast Fourier Transform (FFT), available via scipy.fft.dst(x, type=1, norm="ortho"). The DST matrix is orthogonal and symmetric when properly normalized (norm="ortho"), so $ S = S^T = S^{-1} $.

Your task is to write a solver for Poisson's equation with $ \mathcal{O}(n^2 \log n) $ complexity using these properties of $ T $.

Programming Note: To compute $ SX $ (where $ X $ is a matrix), use scipy.fft.dst(X, type=1, norm="ortho", axis=1). To compute $ XS$, use scipy.fft.dst(X, type=1, norm="ortho", axis=0). You can use the dst(X, axis=0) function defined below to keep the syntax nicer. A function which returns the eigenvalues of $T$ is also provided for you below.

import numpy as np
import matplotlib.pyplot as plt
import scipy.fft

def fast_2d_poisson_solver(F):
    # ADD YOUR CODE HERE
    X = np.zeros_like(F)
    return X

## ALL THE CODE BELOW IS PROVIDED FOR YOU.

# A helper function you can use if you want
def dst(X,axis =0):
    return scipy.fft.dst(X, type = 1, norm = "ortho", axis =axis)

# Computes eigenvalues of the Tridiagonal matrix from known analytical formula
def get_eigenvalues_of_T(n):
    return 4 * (np.sin(np.arange(1,n+1) * np.pi / (2*( n + 1)))**2)

# Get the matrix explicitely, useful for debugging
def get_second_order_diff_matrix(n):
    T = np.diag(2 * np.ones(n)) + np.diag(-1 *  np.ones(n-1), k = 1) + np.diag(-1 * np.ones(n-1) , k = -1)
    return T

### TEST FOR CORRECTNESS

n = 100
h = 1 / (n+1)
from scipy.ndimage import gaussian_filter
F = gaussian_filter(np.random.randn(n,n) * h**2 , 3) # An interesting looking RHS function f

# O(n^3) way using your Bartles-Stewart code:
T = get_second_order_diff_matrix(n)
X = Bartels_Stewart(T, -T, F)
print('Bartels-Stewart error (if this doesnt work, fix previous question!)',np.linalg.norm(T @X + X@T - F) / np.linalg.norm(F))

# O(n^2 log(n)) way using your fast Poisson code:
X2 = fast_2d_poisson_solver(F)
print('Fast poisson solver error', np.linalg.norm(T @X2 + X2@T - F) / np.linalg.norm(F))


plt.figure()
plt.imshow(F)
plt.title("Right hand side function f(x,y)")
plt.show()


plt.figure()
plt.imshow(X2)
plt.title("Approximate solution u(x,y)")
plt.show()


### TEST FOR COMPLEXITY

### A scaling plot

import time
# Define the scaling experiment
# This one takes 1s on my laptop
sizes = [10, 20, 50, 100, 200, 300]  # Increase n gradually

# If you think your method is right, uncomment this! It takes ~10s on my laptop
# sizes = [10, 20, 50, 100, 200, 300, 500, 1000]  # Increase n gradually

bartels_times = []
fast_times = []

for n in sizes:
    h = 1 / (n + 1)
    F = np.ones((n, n)) * h**2

    # Time Bartels-Stewart solver
    T = get_second_order_diff_matrix(n)
    start = time.time()
    X_bartels = Bartels_Stewart(T, -T, F)
    bartels_times.append(time.time() - start)

    # Time the fast Poisson solver
    start = time.time()
    X_fast = fast_2d_poisson_solver(F)
    fast_times.append(time.time() - start)

# Create reference lines for O(n^2) and O(n^3) scaling
sizes_np = np.array(sizes)
o_n2_ref = (sizes_np ** 2) * (fast_times[4] / sizes_np[4] ** 2)  # Normalized to first data point
o_n3_ref = (sizes_np ** 3) * (bartels_times[4] / sizes_np[4] ** 3)  # Normalized to first data point

# Plot the scaling results with reference lines
plt.figure(figsize=(10, 6))
plt.plot(sizes, bartels_times, label="Bartels-Stewart (O(n^3))", marker="o")
plt.plot(sizes, fast_times, label="Fast Poisson Solver (O(n^2 log(n)))", marker="o")
plt.plot(sizes, o_n2_ref, '--', label="O(n^2) Reference", color='gray')
plt.plot(sizes, o_n3_ref, '--', label="O(n^3) Reference", color='black')

# Add plot labels and formatting
plt.xlabel("Matrix Size (n)")
plt.ylabel("Computation Time (seconds)")
plt.title("Scaling of 2D Poisson Solvers")
plt.legend()
plt.xscale("log")  # Set x-axis to logarithmic scale
plt.yscale("log")  # Set y-axis to logarithmic scale
plt.grid(True, which="both", linestyle="--", linewidth=0.5)
plt.show()

Bartels-Stewart error (if this doesnt work, fix previous question!) 1.0
Fast poisson solver error 1.0

Q4. Crafting the jewel

We spent two weeks deriving one of the jewels of numerical linear algebra: the shifted-QR algorithm. Time to pull all the pieces together. When applied to real symmetric matrices, the algorithm looks as follows.

$Q_0 A^{(0)} Q_0 ^T = A \qquad $ (1) (Compute a reduction to tridiagonal by Householder Reflection)

$M \leftarrow m$

for $k = 1,2, \dots$

$\mu_k = A_{M,M}^{(k-1)} \qquad $ (Compute Rayleigh quotient)

$Q^{(k)}, R^{(k)} \leftarrow A^{(k-1)} - \mu_k I \qquad $ (2) (Compute QR decomposition of the shifted matrix)

$A^{(k)} \leftarrow R^{(k)} Q^{(k)} + \mu_k I \qquad $ (3) (Recombine factors)

if $|A_{M,M-1}^{(k)} | < \text{tol}$: $\qquad $ (4) (If converged to eigenvalue)

Add $A_{M,M}^{(k)}$ to list of eigenvalues

$A^{(k)} \leftarrow A^{(k)}_{[1:M-1, 1:M-1]} \qquad $ (Deflate: ignore last row and column)

$M \leftarrow M-1$

if $ M = 1: $ (The trivial 1x1 matrix case: the entry is the eigenvalue)

Add $A_{1,1}^{(k)}$ to list of eigenvalues

terminate

end while

• Step (1) was implemented in precept 8 in $O(m^3)$

• Step (2) was implemented in HW3 $O(m^2)$ for (nonsymmetric). Use your implementation from HW2 (or better yet, modify it to make it $O(m)$ for tridiagonal (optional)!).

• Step (3) Step 3 is a matrix-matrix multiply which is $O(m^3)$ naively. However, it could be implemented in $O(m)$ since it is a matrix-matrix multiply between two very structured and sparse matrices (optional).

• Step (4) is the deflation step. We check if we have converged to an eigenvalue.

As a result of the shifting we recover the $O(\epsilon^3)$ converge of the Rayleigh quotient iteration. Thus, we expect to need only a handful of iterations per eigenvalues (3-4).

Your task is to implement the shifted-QR algorithm above and keep track of how many QR iterations you are doing. Make sure that you only need $< 4 n$ QR iterations to compute all eigenvalues.

def shifted_QR_iteration(A, tol =1e-13):
    eigenvalues = []
    n_iters = 0 
    ## ADD YOUR CODE HERE


    while ():
        ## You can replace this loop, with whatever you want, but make sure you keep track of the number of iterations 
        # (specifically, the number of QR decomposition you are computing)
        n_iters +=1 

    eigenvalues = np.zeros(A.shape[0])
    return eigenvalues,  n_iters

n = 100
# A test case: a random symmetric matrix
A = np.random.randn(n,n)
A = A + A.T 

# Compute eigenvalues with built-in:
eigvals1 = np.linalg.eigvals(A)

# Compute with your code
eigvals2, n_iters = shifted_QR_iteration(A.copy())

# Check that the number eigenvalues are correct:
print('eigenvalue error:', np.max(np.sort(eigvals1) - np.sort(eigvals2)))

# Number of n_iters divided by n:
# Should be < 4
print('number of QR iterations per eigenvalues:', n_iters/n)

eigenvalue error: 27.125159747208166
number of QR iterations per eigenvalues: 0.0

Q5. Spectral Embedding

In this problem we explore nonlinear dimensionality reduction for high-dimensional data. Our goal is to represent a complex dataset using only one parameter. As a motivating example, consider a disordered sequence of images showing a rotating duck: each image lives in a high-dimensional pixel space, yet the underlying variation can be described by a single parameter (the rotation angle).

from PIL import Image
import glob
import matplotlib.pyplot as plt
import numpy as np
import os

# Load images from the spinning_duck directory in their existing order
files = sorted(glob.glob("./spinning_duck/*.png"))

# Load the images as grayscale and stack them
frames = np.stack([np.array(Image.open(f).convert("L"), dtype=np.float64) for f in files])
N, H, W = frames.shape
print(f"Loaded {N} images from './spinning_duck' with size {H}x{W}")

# Visualize a few frames to check
fig, axes = plt.subplots(1, 5, figsize=(15, 3))
for ax, idx in zip(axes, [0, N//5, 2*N//5, 3*N//5, 4*N//5, N-1]):
    ax.imshow(frames[idx], cmap='gray')
    ax.set_title(f'Frame {idx}')
    ax.axis('off')
plt.tight_layout()
plt.show()

Loaded 47 images from './spinning_duck' with size 128x128

We will construct a one-dimensional embedding for such data using the eigenvectors of a similarity matrix.

Let $W$ be the similarity matrix, with entries

$$ w_{ij} = \exp\left(-\frac{\|y_i - y_j\|_2^2}{2\sigma^2}\right), \qquad y_i \in \mathbb{R}^{128^2}. $$ where $\{y_i\}_{i=1}^{N}$ are images and $\sigma$ is a fixed bandwidth parameter.

The only assumptions we need to apply this dimensional reduction method is that the similarity matrix $W$ is symmetric with nonnegative entries; it need not come from a Gaussian kernel specifically.

We wish to assign each image $y_i$ a scalar value $x_i \in \mathbb{R}$ such that similar images are mapped close together. To quantify this, define the embedding vector $x = (x_1, \dots, x_N)$ and energy

$$ E(x) = \sum_{i,j=1}^N w_{ij} (x_i - x_j)^2. $$

Large weights $w_{ij}$ strongly penalize separating $x_i$ and $x_j$. Our task is to find $x$ that minimizes $E(x)$.

Part A

Let $d = W\mathbf{1}$, where $\mathbf{1}$ is the all-ones vector, and define $D = \mathrm{diag}(d)$.

Find a matrix $L$ in terms of $W$ and $D$ such that

$$ E(x) = 2\, x^\top L x. $$

Prove that $L$ is symmetric and positive semidefinite.

Part B

The minimizer $x = \alpha \mathbf{1}$ (for any $\alpha$) is useless — it maps all data points to the same location. Additionally, scaling $x$ shrinks $E(x)$ arbitrarily. To avoid this, we constrain the solution:

$$ \min_{x \in \mathbb{R}^N} x^\top L x \qquad \text{s.t.} \qquad \|x\|_2 = 1, \quad \mathbf{1}^\top x = 0. $$

Introducing Lagrange multipliers gives the Lagrangian

$$ \mathcal{L}(x,\lambda,\mu) = x^\top L x - \lambda (\|x\|_2^2 - 1) - \mu (\mathbf{1}^\top x). $$

Derive the first-order optimality conditions and show that any critical point satisfies

$$ Lx = \lambda x, \qquad \mathbf{1}^\top x = 0. $$

Part C

Express the constrained minimizer in terms of the eigenvectors of $L$, assuming $L$ has a single zero eigenvalue (i.e., the nullspace of $L$ is spanned by $\mathbf{1}$).

Write your answer in here

Part A

Part B

Part C

Part D (Implementation and Visualization)

Implement your solution from Part C using the duck dataset. You may use np.linalg.eig or scipy.linalg.eigh.

def compute_spectral_embedding(W, D):
    '''
    Return the vector $x$ that minimizes $E(x)$ subject to the constraints.
    The vector $x$ should have the same length as the number of frames
    Input: 
        W: a numpy array of shape (N, N): the similarity matrix
        D: a numpy array of shape (N, N): D = diag(d) where d is the sum of the rows of W
    Output: 
        x: a numpy array of shape (N,)
    '''
    # Return the embedding vector $x$
    embedding = np.zeros(W.shape[0])
    return embedding



## This code is provided for you. No need to change it!

## Compute the similarity matrix $W$ and the diagonal matrix $D$ from the frames
def compute_W_D_from_frames(frames):
    N = frames.shape[0]
    X = frames.reshape(N, -1)
    sqn = np.sum(X * X, axis=1, keepdims=True)
    dist2 = np.maximum(sqn + sqn.T - 2 * (X @ X.T), 0.0)
    sigma = np.median(np.sqrt(dist2 + 1e-18)[~np.eye(N, dtype=bool)]) * 0.2
    W = np.exp(-dist2 / (2 * sigma**2))
    d = W.sum(axis=1)
    D = np.diag(d)
    return W, D


def compute_energy_function_naive(x, W):
    N = W.shape[0]
    energy = 0
    for i in range(N):
        for j in range(N):
            energy += W[i, j] * (x[i] - x[j])**2
    return energy


W, D = compute_W_D_from_frames(frames)


## A few tests for your code
embedding = compute_spectral_embedding(W, D)

# Is \|x\|_2 = 1?
if not np.isclose(np.linalg.norm(embedding), 1):
    print("Error: The norm of the embedding is not 1! Constraint is not satisfied!")
else:
    print("The norm of the embedding is 1")

# Is \mathbf{1}^\top x = 0?
if not np.isclose(np.sum(embedding), 0):
    print("Error: The sum of the embedding is not 0! Constraint is not satisfied!")
else:
    print("The sum of the embedding is 0")


if compute_energy_function_naive(embedding, W)  > 0.13:
    print("Error: The energy function is too large! You are not at the minimum!")
else:
    print("The energy function is small enough!")




#  Reorder the frames according to the embedding
sorted_indices = np.argsort(embedding)
frames_ordered = frames[sorted_indices]

#. Plot the sequence of images ordered by the embedding
num_show = min(10, N)  # Show up to 10 images
idxs = np.linspace(0, N-1, num_show, dtype=int)


## Plot original ordering
fig, axes = plt.subplots(1, num_show, figsize=(2*num_show, 2))
for i, ax in enumerate(axes):
    ax.imshow(frames[idxs[i]], cmap='gray')
    ax.set_title(f'Order {idxs[i]}')
    ax.axis('off')
plt.suptitle('Duck images in original order')
plt.tight_layout()
plt.show()



## Plot ordered by spectral embedding
#  It should look like a spinning duck!
fig, axes = plt.subplots(1, num_show, figsize=(2*num_show, 2))
for i, ax in enumerate(axes):
    ax.imshow(frames_ordered[idxs[i]], cmap='gray')
    ax.set_title(f'Order {idxs[i]}')
    ax.axis('off')
plt.suptitle('Duck images ordered by spectral embedding')
plt.tight_layout()
plt.show()

Error: The norm of the embedding is not 1! Constraint is not satisfied!
The sum of the embedding is 0
The energy function is small enough!

def play_duck_movie(frames, interval=0.2, ordering = ''):
    """
    Display a movie of the frames in order with the given interval between frames.
    Also shows the frame number.
    """
    import matplotlib.pyplot as plt
    import time
    from IPython.display import display, clear_output

    fig, ax = plt.subplots()
    im = ax.imshow(frames[0], cmap='gray')
    txt = ax.text(0.5, 1.05, f"Frame 0", color='black', fontsize=14, ha='center', va='bottom', transform=ax.transAxes)
    plt.axis('off')

    for i, frame in enumerate(frames):
        im.set_data(frame)
        txt.set_text(f"Frame {i} - {ordering}")
        clear_output(wait=True)
        display(fig)
        plt.pause(interval)

    plt.close(fig)

# A movie of the images in the embedding order!
play_duck_movie(frames, interval=0.2, ordering='Original ordering')
play_duck_movie(frames_ordered, interval=0.2, ordering='Spectral embedding')

Extra credit papers:

• Randomized algorithms - https://arxiv.org/pdf/2002.01387.pdf
• Divide-and-conquer algorithms for eigenvalues - https://epubs.siam.org/doi/10.1137/S0895479892241287