Chapter 169 — Eigenvectors: Intuition

# --- Visualization: Eigenvectors as direction-preserving vectors ---
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid')

# A matrix with known eigenvectors for demonstration
A = np.array([[3.0, 1.0],
              [0.0, 2.0]])

# True eigenvectors (computed analytically for this matrix)
# A is upper triangular -> eigenvalues are diagonal entries: 3, 2
# eigenvector for lambda=3: [1, 0]
# eigenvector for lambda=2: [1, -1] (normalized)
eigvals, eigvecs = np.linalg.eig(A)

# Sample random vectors on the unit circle
N_VECTORS = 16
angles = np.linspace(0, 2 * np.pi, N_VECTORS, endpoint=False)
vecs = np.stack([np.cos(angles), np.sin(angles)], axis=1)  # shape (N, 2)

# Apply transformation
transformed = (A @ vecs.T).T  # shape (N, 2)

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

for ax, (vset, title) in zip(axes, [(vecs, 'Before: Unit Circle Vectors'),
                                     (transformed, 'After: A applied')]):
    ax.set_xlim(-6, 6)
    ax.set_ylim(-6, 6)
    ax.axhline(0, color='gray', lw=0.5)
    ax.axvline(0, color='gray', lw=0.5)
    ax.set_aspect('equal')
    ax.set_title(title, fontsize=13)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    for v in vset:
        ax.annotate('', xy=v, xytext=(0, 0),
                    arrowprops=dict(arrowstyle='->', color='steelblue', lw=1.2))

# Highlight eigenvectors on both plots
for i, (ev, lam) in enumerate(zip(eigvecs.T, eigvals)):
    ev_norm = ev / np.linalg.norm(ev)
    ev_scaled = ev_norm * lam  # where it maps to
    color = ['red', 'green'][i]
    axes[0].annotate('', xy=ev_norm, xytext=(0, 0),
                     arrowprops=dict(arrowstyle='->', color=color, lw=2.5))
    axes[0].text(ev_norm[0]*1.2, ev_norm[1]*1.2, f'e{i+1}', color=color, fontsize=11)
    axes[1].annotate('', xy=ev_scaled, xytext=(0, 0),
                     arrowprops=dict(arrowstyle='->', color=color, lw=2.5))
    axes[1].text(ev_scaled[0]*1.1, ev_scaled[1]*1.1,
                 f'λ={lam:.1f}·e{i+1}', color=color, fontsize=11)

plt.suptitle('Eigenvectors preserve direction under transformation', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()
print(f"Eigenvalues: {eigvals}")
print(f"Eigenvectors (columns):\n{eigvecs}")

# Worked numeric example
import numpy as np

A = np.array([[3.0, 1.0],
              [0.0, 2.0]])

# Verify eigenvector by hand: v = [1, 0], lambda = 3
v1 = np.array([1.0, 0.0])
lam1 = 3.0

lhs = A @ v1       # A v
rhs = lam1 * v1    # λ v

print("A @ v1 =", lhs)
print("lam1 * v1 =", rhs)
print("Equal?", np.allclose(lhs, rhs))

# Verify second eigenvector: v = [1, -1], lambda = 2
v2 = np.array([1.0, -1.0])
lam2 = 2.0
print("\nA @ v2 =", A @ v2)
print("lam2 * v2 =", lam2 * v2)
print("Equal?", np.allclose(A @ v2, lam2 * v2))

# --- Implementation: Visualizing direction change for all unit vectors ---
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid')

def angle_change(A, n_samples=360):
    """
    For each unit vector at angle theta, compute how much its direction
    changes after applying matrix A.

    Args:
        A: (2,2) matrix
        n_samples: number of angles to sample

    Returns:
        angles: input angles in degrees
        delta: direction change in degrees after applying A
    """
    angles = np.linspace(0, 360, n_samples, endpoint=False)
    deltas = []
    for theta in angles:
        rad = np.radians(theta)
        v = np.array([np.cos(rad), np.sin(rad)])
        Av = A @ v
        # Angle of transformed vector
        out_angle = np.degrees(np.arctan2(Av[1], Av[0])) % 360
        delta = abs((out_angle - theta + 180) % 360 - 180)  # signed wrap
        deltas.append(delta)
    return angles, np.array(deltas)

A = np.array([[3.0, 1.0], [0.0, 2.0]])
angles, deltas = angle_change(A)

# Find angles where direction change ≈ 0 (eigenvector directions)
eigen_angles = angles[deltas < 2.0]

plt.figure(figsize=(10, 4))
plt.plot(angles, deltas, 'steelblue', lw=1.5, label='Direction change (degrees)')
for ea in eigen_angles:
    plt.axvline(ea, color='red', ls='--', alpha=0.6, lw=1)
plt.xlabel('Input vector angle (degrees)')
plt.ylabel('Direction change after A (degrees)')
plt.title('Direction change under A — zeros = eigenvector directions')
plt.legend()
plt.tight_layout()
plt.show()
print(f"Approximate eigenvector directions: {eigen_angles}°")

# --- Experiment 1: Rotation matrix has no real eigenvectors ---
# Hypothesis: A pure rotation has no direction that is preserved
# Try changing: THETA to 0 or 180 degrees and observe what happens

import numpy as np
THETA = 45  # degrees -- try: 0, 90, 180, 45

rad = np.radians(THETA)
R = np.array([[np.cos(rad), -np.sin(rad)],
              [np.sin(rad),  np.cos(rad)]])

eigvals = np.linalg.eigvals(R)
print(f"Rotation by {THETA}°")
print(f"Eigenvalues: {eigvals}")
print(f"Real? {np.all(np.isreal(eigvals))}")
print("Note: complex eigenvalues mean no real eigenvectors")

# --- Experiment 2: Symmetric matrix always has real eigenvectors ---
# Hypothesis: Symmetric matrices always have real, orthogonal eigenvectors
# Try changing: the off-diagonal entries

import numpy as np

OFF_DIAG = 2.0  # try: 0, 1, 3, -2

S = np.array([[4.0, OFF_DIAG],
              [OFF_DIAG, 1.0]])  # symmetric by construction

eigvals, eigvecs = np.linalg.eig(S)
print(f"Matrix S with off-diagonal = {OFF_DIAG}:")
print(S)
print(f"Eigenvalues: {eigvals}")
print(f"Eigenvectors (columns):\n{eigvecs}")
print(f"Dot product of eigenvectors (should be 0): {eigvecs[:,0] @ eigvecs[:,1]:.6f}")

# --- Experiment 3: Power iteration finds dominant eigenvector ---
# Hypothesis: Repeatedly applying A and normalizing converges to the
#             eigenvector with the largest eigenvalue
# Try changing: N_ITER and the initial vector

import numpy as np

N_ITER = 20  # try: 5, 10, 50

A = np.array([[3.0, 1.0], [0.0, 2.0]])
v = np.array([0.5, 0.5])  # arbitrary starting vector
v = v / np.linalg.norm(v)

history = [v.copy()]
for _ in range(N_ITER):
    v = A @ v
    v = v / np.linalg.norm(v)  # normalize to stay on unit sphere
    history.append(v.copy())

true_evec = np.linalg.eig(A)[1][:, 0]  # dominant eigenvector
true_evec /= np.linalg.norm(true_evec)

print(f"After {N_ITER} iterations: {v}")
print(f"True dominant eigenvector: {true_evec}")
print(f"Angle between them: {np.degrees(np.arccos(abs(v @ true_evec))):.4f}°")

# --- Mini Project: Visualize power iteration convergence ---
# Problem: Animate how a random vector converges to the dominant eigenvector
#          of a 2x2 matrix under repeated normalized multiplication.
# Task: Complete the animation and measure convergence rate.

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

# TODO: Try modifying this matrix
A = np.array([[4.0, 2.0],
              [1.0, 3.0]])

eigvals, eigvecs = np.linalg.eig(A)
dominant_idx = np.argmax(np.abs(eigvals))
true_evec = eigvecs[:, dominant_idx]
true_evec /= np.linalg.norm(true_evec)

# Run power iteration
np.random.seed(42)
v = np.random.randn(2)
v /= np.linalg.norm(v)

N_STEPS = 15
history = [v.copy()]
errors = []

for _ in range(N_STEPS):
    v = A @ v
    v /= np.linalg.norm(v)
    history.append(v.copy())
    # Angle error (account for sign ambiguity)
    err = np.degrees(np.arccos(min(1.0, abs(v @ true_evec))))
    errors.append(err)

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

# Plot trajectory of iterates
ax = axes[0]
ax.set_xlim(-1.5, 1.5)
ax.set_ylim(-1.5, 1.5)
ax.set_aspect('equal')
theta = np.linspace(0, 2*np.pi, 200)
ax.plot(np.cos(theta), np.sin(theta), 'gray', lw=0.5, alpha=0.4)
for i, h in enumerate(history):
    alpha = 0.2 + 0.8 * i / len(history)
    ax.annotate('', xy=h, xytext=(0,0),
                arrowprops=dict(arrowstyle='->', color='steelblue', lw=1.5, alpha=alpha))
ax.annotate('', xy=true_evec, xytext=(0,0),
            arrowprops=dict(arrowstyle='->', color='red', lw=2.5))
ax.text(true_evec[0]+0.05, true_evec[1]+0.05, 'True eigenvector', color='red')
ax.set_title('Power iteration trajectory')
ax.set_xlabel('x'); ax.set_ylabel('y')

# Plot error over iterations
axes[1].semilogy(range(1, N_STEPS+1), np.array(errors)+1e-10, 'o-', color='tomato')
axes[1].set_xlabel('Iteration')
axes[1].set_ylabel('Angular error (degrees, log scale)')
axes[1].set_title('Convergence rate of power iteration')

plt.tight_layout()
plt.show()
print(f"Eigenvalues: {eigvals}")
print(f"Ratio |λ2/λ1| = {min(abs(eigvals))/max(abs(eigvals)):.3f}")
print("Convergence rate is determined by this ratio")