Chapter 171 — Eigenvalue Computation

# --- Visualization: QR iteration converging to upper triangular ---
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid')

np.random.seed(42)
# Symmetric matrix for clean real eigenvalues
M = np.random.randn(4, 4)
A = M @ M.T  # positive definite symmetric

N_ITER = 25
off_diag_norms = []  # measure how "upper triangular" matrix becomes
diagonals = []       # track diagonal entries

Ak = A.copy()
for _ in range(N_ITER):
    Q, R = np.linalg.qr(Ak)
    Ak = R @ Q  # similarity transform: same eigenvalues
    # Off-diagonal Frobenius norm (strictly lower triangle)
    lower = np.tril(Ak, -1)
    off_diag_norms.append(np.linalg.norm(lower, 'fro'))
    diagonals.append(np.diag(Ak).copy())

true_eigs = sorted(np.linalg.eigvals(A).real, reverse=True)

fig, axes = plt.subplots(1, 2, figsize=(13, 5))

axes[0].semilogy(off_diag_norms, 'o-', color='steelblue', markersize=4)
axes[0].set_xlabel('QR Iteration')
axes[0].set_ylabel('Off-diagonal Frobenius norm (log scale)')
axes[0].set_title('QR algorithm convergence to upper triangular')

diagonals = np.array(diagonals)
colors = ['red', 'green', 'purple', 'orange']
for i in range(4):
    axes[1].plot(diagonals[:, i], color=colors[i], lw=1.5, label=f'diag[{i}]')
    axes[1].axhline(true_eigs[i], color=colors[i], ls='--', alpha=0.5)
axes[1].set_xlabel('QR Iteration')
axes[1].set_ylabel('Diagonal entry value')
axes[1].set_title('Diagonal entries → eigenvalues (dashed)')
axes[1].legend()

plt.tight_layout()
plt.show()

# --- Implementation: Power iteration and QR algorithm from scratch ---
import numpy as np

def power_iteration(A, n_iter=100, tol=1e-10, seed=0):
    """
    Find the dominant eigenvalue and eigenvector of A.

    Args:
        A: (n,n) square matrix
        n_iter: maximum iterations
        tol: convergence tolerance on eigenvalue change

    Returns:
        (eigenvalue, eigenvector, n_iterations)
    """
    np.random.seed(seed)
    n = A.shape[0]
    v = np.random.randn(n)
    v /= np.linalg.norm(v)

    lam_prev = 0.0
    for k in range(n_iter):
        Av = A @ v
        lam = v @ Av               # Rayleigh quotient
        v = Av / np.linalg.norm(Av)
        if abs(lam - lam_prev) < tol:
            return lam, v, k+1
        lam_prev = lam

    return lam, v, n_iter


def qr_algorithm(A, n_iter=50):
    """
    Compute all eigenvalues of A using unshifted QR iteration.

    Args:
        A: (n,n) square matrix
        n_iter: number of QR steps

    Returns:
        Approximate eigenvalues (diagonal of converged upper triangular)
    """
    Ak = A.copy().astype(float)
    for _ in range(n_iter):
        Q, R = np.linalg.qr(Ak)
        Ak = R @ Q
    return np.diag(Ak)


# Test on a known matrix
np.random.seed(5)
M = np.random.randn(5, 5)
A = M @ M.T  # symmetric positive definite

# Power iteration (dominant only)
lam_dom, v_dom, iters = power_iteration(A)
print(f"Power iteration: λ = {lam_dom:.6f} ({iters} iters)")

# QR algorithm (all eigenvalues)
our_eigs = sorted(qr_algorithm(A, n_iter=100), reverse=True)
numpy_eigs = sorted(np.linalg.eigvals(A).real, reverse=True)

print("\nQR algorithm vs NumPy:")
for i, (ours, ref) in enumerate(zip(our_eigs, numpy_eigs)):
    print(f"  λ{i+1}: ours={ours:.6f}  numpy={ref:.6f}  diff={abs(ours-ref):.2e}")

# --- Experiment 1: Convergence rate of power iteration ---
# Hypothesis: Convergence slows when eigenvalue gap is small
# Try changing: GAP (difference between the two largest eigenvalues)

import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid')

GAP = 0.5  # try: 0.1, 0.5, 2.0, 0.01
LAM1 = 5.0
LAM2 = LAM1 - GAP

np.random.seed(3)
V = np.random.randn(2, 2)
A = V @ np.diag([LAM1, LAM2]) @ np.linalg.inv(V)

# Track error in eigenvalue estimate
v = np.random.randn(2); v /= np.linalg.norm(v)
errors = []
for _ in range(80):
    Av = A @ v
    lam_est = v @ Av
    errors.append(abs(lam_est - LAM1))
    v = Av / np.linalg.norm(Av)

plt.figure(figsize=(8, 4))
plt.semilogy(errors, 'steelblue')
plt.xlabel('Iteration')
plt.ylabel('|λ_est - λ_true| (log)')
plt.title(f'Power iteration convergence: gap = {GAP}, ratio = {LAM2/LAM1:.3f}')
plt.tight_layout()
plt.show()

# --- Experiment 2: Condition number and eigenvalue sensitivity ---
# Hypothesis: Adding small noise to a nearly-singular matrix causes large eigenvalue changes
# Try changing: EPSILON

import numpy as np

EPSILON = 1e-6  # try: 1e-1, 1e-3, 1e-6, 1e-10

# Nearly singular matrix
A = np.array([[1.0, 1.0],
              [1.0, 1.0 + EPSILON]])

eigs = np.linalg.eigvals(A)
cond = np.linalg.cond(A)
print(f"Epsilon = {EPSILON:.1e}")
print(f"Eigenvalues: {sorted(eigs, reverse=True)}")
print(f"Condition number: {cond:.2e}")
print(f"Smallest eigenvalue: {min(abs(eigs)):.2e}")

# --- Mini Project: Eigenvalues of an adjacency matrix ---
# Problem: The eigenvalues of a graph's adjacency matrix encode its structure.
#          The largest eigenvalue relates to connectivity; the second-smallest
#          Laplacian eigenvalue reveals cluster structure.
# Task: Build a small graph, compute its Laplacian, find eigenvalues,
#       and use them to detect clusters.

import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid')

# Adjacency matrix: two clusters (nodes 0-2 and 3-5) with a weak bridge
n = 6
A_adj = np.zeros((n, n))

# Dense connections within clusters
for i in range(3):
    for j in range(i+1, 3):
        A_adj[i, j] = A_adj[j, i] = 1
for i in range(3, 6):
    for j in range(i+1, 6):
        A_adj[i, j] = A_adj[j, i] = 1

# Weak bridge between clusters
A_adj[2, 3] = A_adj[3, 2] = 0.1  # TODO: try 0, 0.5, 1.0

# Laplacian: L = D - A where D = degree matrix
D = np.diag(A_adj.sum(axis=1))
L = D - A_adj

eigvals, eigvecs = np.linalg.eigh(L)  # eigh for symmetric

print("Laplacian eigenvalues:")
for i, ev in enumerate(eigvals):
    print(f"  λ{i} = {ev:.4f}")

# The Fiedler vector: eigenvector of second-smallest eigenvalue
fiedler_vec = eigvecs[:, 1]

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].bar(range(n), fiedler_vec, color=['steelblue' if f > 0 else 'tomato' for f in fiedler_vec])
axes[0].axhline(0, color='black', lw=0.8)
axes[0].set_xlabel('Node')
axes[0].set_ylabel('Fiedler vector value')
axes[0].set_title('Fiedler vector: positive=cluster A, negative=cluster B')

axes[1].bar(range(n), eigvals, color='steelblue')
axes[1].set_xlabel('Eigenvalue index')
axes[1].set_ylabel('λ')
axes[1].set_title('Laplacian spectrum')

plt.tight_layout()
plt.show()
print(f"\nFiedler value (λ1={eigvals[1]:.4f}): connectivity measure")
print("Signs of Fiedler vector indicate cluster membership:")
print(['+' if f > 0 else '-' for f in fiedler_vec])