Chapter 170 — Eigenvalues: Intuition

# --- Visualization: Unit circle → ellipse, eigenvalues as axis lengths ---
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid')

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

eigvals, eigvecs = np.linalg.eig(A)

# Unit circle
theta = np.linspace(0, 2*np.pi, 300)
circle = np.stack([np.cos(theta), np.sin(theta)])  # (2, 300)
ellipse = A @ circle  # (2, 300)

fig, ax = plt.subplots(figsize=(8, 8))
ax.plot(circle[0], circle[1], 'steelblue', lw=1.5, label='Unit circle (input)')
ax.plot(ellipse[0], ellipse[1], 'tomato', lw=1.5, label='Ellipse (output = A·circle)')

# Draw eigenvectors scaled by their eigenvalues
colors = ['green', 'purple']
for i in range(2):
    ev = eigvecs[:, i].real
    lam = eigvals[i].real
    ev_unit = ev / np.linalg.norm(ev)
    ax.annotate('', xy=ev_unit, xytext=(0, 0),
                arrowprops=dict(arrowstyle='->', color=colors[i], lw=2, alpha=0.7))
    ax.annotate('', xy=ev_unit * lam, xytext=(0, 0),
                arrowprops=dict(arrowstyle='->', color=colors[i], lw=3))
    ax.text(ev_unit[0]*lam*1.1, ev_unit[1]*lam*1.1,
            f'λ{i+1}={lam:.2f}', color=colors[i], fontsize=12)

ax.set_xlim(-5, 5)
ax.set_ylim(-5, 5)
ax.set_aspect('equal')
ax.axhline(0, color='gray', lw=0.5)
ax.axvline(0, color='gray', lw=0.5)
ax.set_title('Eigenvalues = scaling factors along eigenvector directions', fontsize=13)
ax.set_xlabel('x'); ax.set_ylabel('y')
ax.legend()
plt.tight_layout()
plt.show()
print(f"det(A) = {np.linalg.det(A):.4f}, λ1*λ2 = {np.prod(eigvals.real):.4f}")
print(f"trace(A) = {np.trace(A):.4f}, λ1+λ2 = {np.sum(eigvals.real):.4f}")

# Characteristic polynomial from scratch (2x2 case)
import numpy as np

def eigenvalues_2x2(A):
    """
    Compute eigenvalues of a 2x2 matrix via the characteristic polynomial.

    Args:
        A: (2,2) numpy array

    Returns:
        Tuple of two eigenvalues (may be complex)
    """
    tr = A[0,0] + A[1,1]          # trace
    det = A[0,0]*A[1,1] - A[0,1]*A[1,0]  # determinant
    discriminant = tr**2 - 4*det
    sqrt_disc = np.sqrt(complex(discriminant))
    lam1 = (tr + sqrt_disc) / 2
    lam2 = (tr - sqrt_disc) / 2
    return lam1, lam2

# Test
A = np.array([[3.0, 1.0], [0.5, 2.0]])
lam1, lam2 = eigenvalues_2x2(A)
numpy_vals = np.linalg.eigvals(A)

print(f"Our eigenvalues:   {lam1:.4f}, {lam2:.4f}")
print(f"NumPy eigenvalues: {sorted(numpy_vals.real, reverse=True)}")

# Verify trace and det relations
print(f"\ntrace(A) = {np.trace(A):.4f}, λ1+λ2 = {(lam1+lam2).real:.4f}")
print(f"det(A)   = {np.linalg.det(A):.4f}, λ1*λ2 = {(lam1*lam2).real:.4f}")

# --- Implementation: Visualize the characteristic polynomial ---
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid')

def char_poly_2x2(A, lam):
    """
    Evaluate det(A - lambda*I) for a 2x2 matrix at given lambda values.

    Args:
        A: (2,2) numpy array
        lam: scalar or array of lambda values

    Returns:
        det(A - lambda*I) for each lambda
    """
    lam = np.atleast_1d(lam)
    results = []
    for l in lam:
        M = A - l * np.eye(2)
        results.append(M[0,0]*M[1,1] - M[0,1]*M[1,0])
    return np.array(results)

A = np.array([[3.0, 1.0], [0.5, 2.0]])
lam_range = np.linspace(-1, 6, 500)
p = char_poly_2x2(A, lam_range)

true_eigs = np.linalg.eigvals(A).real

plt.figure(figsize=(9, 5))
plt.plot(lam_range, p, 'steelblue', lw=2, label='det(A - λI)')
plt.axhline(0, color='black', lw=0.8)
for ev in true_eigs:
    plt.axvline(ev, color='red', ls='--', lw=1.5)
    plt.scatter([ev], [0], color='red', zorder=5, s=60)
    plt.text(ev+0.05, 0.5, f'λ={ev:.2f}', color='red', fontsize=10)
plt.xlabel('λ')
plt.ylabel('det(A - λI)')
plt.title('Characteristic polynomial — roots are eigenvalues')
plt.legend()
plt.tight_layout()
plt.show()

# --- Experiment 1: Effect of eigenvalue sign on dynamics ---
# Hypothesis: Negative eigenvalue causes oscillating behavior under repeated application
# Try changing: LAMBDA to positive, zero, or complex

import numpy as np

# Construct matrix with desired eigenvalues
LAM1 = -0.8  # try: 1.2, 0.5, 0.0, -0.8
LAM2 = 0.6

# Build A = V * diag(lam) * V^-1 with random V
np.random.seed(7)
V = np.random.randn(2, 2)
A = V @ np.diag([LAM1, LAM2]) @ np.linalg.inv(V)

# Iterate a vector
x = np.array([1.0, 0.5])
N_ITER = 20
norms = [np.linalg.norm(x)]
for _ in range(N_ITER):
    x = A @ x
    norms.append(np.linalg.norm(x))

import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid')
plt.figure(figsize=(8, 4))
plt.plot(norms, 'o-', color='steelblue')
plt.xlabel('Iteration')
plt.ylabel('||A^n x||')
plt.title(f'Vector norm evolution (λ1={LAM1}, λ2={LAM2})')
plt.tight_layout()
plt.show()

# --- Experiment 2: Eigenvalues of a rank-1 matrix ---
# Hypothesis: A rank-1 matrix has one nonzero eigenvalue and (n-1) zeros
# Try changing: the vectors u and v

import numpy as np

u = np.array([1.0, 2.0, 3.0])  # try: [1,0,0], [1,1,1]
v = np.array([4.0, 1.0, 2.0])

# Rank-1 matrix: outer product
A = np.outer(u, v)  # shape (3, 3), rank 1

eigvals = np.linalg.eigvals(A)
print(f"Rank of A: {np.linalg.matrix_rank(A)}")
print(f"Eigenvalues: {eigvals.real}")
print(f"Nonzero eigenvalue = u·v = {u @ v:.4f}")
print(f"trace(A) = {np.trace(A):.4f} = sum of eigenvalues = {eigvals.real.sum():.4f}")

# --- Mini Project: Map eigenvalue type to geometric behavior ---
# Problem: Build a visualization that shows how changing a 2x2 matrix's
#          entries changes the eigenvalue type (real/complex/repeated)
#          and what each type looks like geometrically.

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

def classify_eigenvalues(A):
    """
    Classify a 2x2 matrix's eigenvalues.

    Returns: 'two_real', 'repeated', or 'complex'
    """
    tr = np.trace(A)
    det = np.linalg.det(A)
    disc = tr**2 - 4*det
    if abs(disc) < 1e-10:
        return 'repeated'
    elif disc > 0:
        return 'two_real'
    else:
        return 'complex'

# Sweep off-diagonal entries
K_VALUES = np.linspace(-3, 3, 200)
results = {'two_real': [], 'repeated': [], 'complex': []}

for k in K_VALUES:
    A = np.array([[2.0, k], [-k, 2.0]])
    etype = classify_eigenvalues(A)
    results[etype].append(k)

# Visualize characteristic polynomial for three cases
cases = [
    ('two_real', 0.5),
    ('repeated', 0.0),   # k=0 gives repeated λ=2
    ('complex', 1.5)
]

fig, axes = plt.subplots(1, 3, figsize=(15, 4))
lam_range = np.linspace(-2, 6, 400)

for ax, (label, k) in zip(axes, cases):
    A = np.array([[2.0, k], [-k, 2.0]])
    p = [(A[0,0]-l)*(A[1,1]-l) - A[0,1]*A[1,0] for l in lam_range]
    ax.plot(lam_range, p, 'steelblue', lw=2)
    ax.axhline(0, color='black', lw=0.8)
    eigs = np.linalg.eigvals(A)
    for ev in eigs.real:
        ax.axvline(ev, color='red', ls='--', alpha=0.7)
    ax.set_title(f'{label} (k={k})')
    ax.set_xlabel('λ')
    ax.set_ylabel('det(A-λI)')
    ax.set_ylim(-5, 10)

plt.suptitle('Eigenvalue type from characteristic polynomial shape', fontsize=13)
plt.tight_layout()
plt.show()