Eigenvalues and Eigenvectors

Think of a shadow. You hold an odd-shaped object under a light. As you rotate the object, its shadow changes shape — sometimes long and thin, sometimes short and wide. There is one specific angle where the shadow is longest. That angle is like an eigenvector — the direction that reveals the most about the object's structure. The length of that longest shadow is like the eigenvalue.

PCA (Principal Component Analysis) uses exactly this idea. It finds the directions in your data where the shadow is longest — where the data varies the most — and uses those as the new axes. Module 33 covers PCA in full. This module gives you the mathematical foundation to understand how it works.

A vector e is an eigenvector of matrix A if multiplying A by e gives back e scaled by a number λ:

import numpy as np

# Define a simple 2×2 matrix
A = np.array([
    [3, 1],
    [1, 3],
])

# ── Manually verify the eigenvector definition ─────────────────────────
e1 = np.array([1, 1])   # eigenvector 1
e2 = np.array([1, -1])  # eigenvector 2

print("Verifying A · e = λ · e for eigenvector 1:")
Ae1  = A @ e1
lam1 = 4
print(f"  A · e1        = {Ae1}")
print(f"  λ₁ × e1 (×4) = {lam1 * e1}")
print(f"  Are they equal? {np.allclose(Ae1, lam1 * e1)}")  # True

print("\nVerifying A · e = λ · e for eigenvector 2:")
Ae2  = A @ e2
lam2 = 2
print(f"  A · e2        = {Ae2}")
print(f"  λ₂ × e2 (×2) = {lam2 * e2}")
print(f"  Are they equal? {np.allclose(Ae2, lam2 * e2)}")  # True

# ── NumPy computes eigenvalues and eigenvectors in one call ────────────
eigenvalues, eigenvectors = np.linalg.eig(A)

print(f"\nNumPy computed eigenvalues:  {eigenvalues}")
print(f"NumPy computed eigenvectors:\n{eigenvectors}")
# Each COLUMN of the eigenvectors matrix is one eigenvector
# Column 0 corresponds to eigenvalues[0], Column 1 to eigenvalues[1]

print(f"\nEigenvector 0 (column 0): {eigenvectors[:, 0].round(4)}")
print(f"Eigenvector 1 (column 1): {eigenvectors[:, 1].round(4)}")
# Note: NumPy returns unit vectors (length=1), so [1,1] becomes [0.707, 0.707]
# They point in the same direction — same eigenvector, just normalised

# ── A random vector is NOT an eigenvector ─────────────────────────────
random_vec = np.array([3, 1])   # not an eigenvector
print(f"\nRandom vector [3, 1]:")
print(f"  A · v = {A @ random_vec}  — direction changed, not an eigenvector")

import numpy as np

np.random.seed(42)

# ── Simulate a DoorDash customer dataset: 1000 customers, 5 features ────
# Features: order_frequency, avg_order_value, delivery_distance,
#           night_orders_pct, weekend_orders_pct
# Features are correlated — frequent orderers tend to spend more

n = 1000
order_freq      = np.random.normal(15, 5, n)
avg_value       = 150 + 8 * order_freq + np.random.normal(0, 30, n)  # correlated
delivery_dist   = np.random.normal(4, 1.5, n)
night_pct       = np.random.normal(0.3, 0.1, n).clip(0, 1)
weekend_pct     = np.random.normal(0.4, 0.1, n).clip(0, 1)

X = np.column_stack([order_freq, avg_value, delivery_dist,
                     night_pct, weekend_pct])

# ── Step 1: Standardise the data (zero mean, unit variance) ───────────
# This is CRITICAL before computing eigenvalues
# Without it, features with larger values (like avg_value in dollars)
# will dominate simply because their numbers are bigger
X_mean = X.mean(axis=0)
X_std  = X.std(axis=0)
X_std_data = (X - X_mean) / X_std
print(f"Data shape: {X_std_data.shape}  (1000 customers × 5 features)")

# ── Step 2: Compute covariance matrix ────────────────────────────────
# Covariance matrix tells us: how do features vary together?
# Shape: (5, 5) — every feature pair has a covariance value
cov_matrix = np.cov(X_std_data.T)   # .T because np.cov expects features as rows
print(f"Covariance matrix shape: {cov_matrix.shape}")

# ── Step 3: Compute eigenvalues and eigenvectors ────────────────────
eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)
# np.linalg.eigh for symmetric matrices (covariance is always symmetric)
# Returns eigenvalues sorted ASCENDING — we want DESCENDING
idx          = np.argsort(eigenvalues)[::-1]
eigenvalues  = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]

# ── Step 4: Explained variance ────────────────────────────────────────
total_variance     = eigenvalues.sum()
explained_variance = eigenvalues / total_variance

feature_names = ['order_freq', 'avg_value', 'dist', 'night_pct', 'weekend_pct']

print(f"\nEigenvalues and explained variance:")
print(f"{'PC':<5} {'Eigenvalue':<14} {'Explained':<12} {'Cumulative':<12}")
print("─" * 46)
cumulative = 0
for i, (ev, pct) in enumerate(zip(eigenvalues, explained_variance)):
    cumulative += pct
    bar = '█' * int(pct * 40)
    print(f"  PC{i+1}  {ev:<14.4f} {pct*100:<10.1f}%  {cumulative*100:.1f}%  {bar}")

print(f"\nConclusion:")
print(f"  PC1 alone explains {explained_variance[0]*100:.1f}% of all variation")
print(f"  PC1 + PC2 explain  {(explained_variance[:2].sum())*100:.1f}% — good enough for most tasks")
print(f"  To preserve 95% variance, keep {np.searchsorted(np.cumsum(explained_variance), 0.95)+1} components")

# ── What does the first principal component 'mean'? ──────────────────
print(f"\nPC1 loadings (how much each feature contributes):")
pc1 = eigenvectors[:, 0]
for name, loading in zip(feature_names, pc1):
    bar = '█' * int(abs(loading) * 20)
    sign = '+' if loading > 0 else '-'
    print(f"  {name:<18}: {sign}{bar} {loading:.4f}")

import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

np.random.seed(42)

# Same DoorDash customer dataset
n = 1000
order_freq    = np.random.normal(15, 5, n)
avg_value     = 150 + 8 * order_freq + np.random.normal(0, 30, n)
delivery_dist = np.random.normal(4, 1.5, n)
night_pct     = np.random.normal(0.3, 0.1, n).clip(0, 1)
weekend_pct   = np.random.normal(0.4, 0.1, n).clip(0, 1)
X = np.column_stack([order_freq, avg_value, delivery_dist, night_pct, weekend_pct])
feature_names = ['order_freq', 'avg_value', 'dist', 'night_pct', 'weekend_pct']

# ── PCA from scratch — using eigenvalues ──────────────────────────────
def pca_from_scratch(X, n_components):
    """
    Full PCA implementation using eigenvalue decomposition.
    Every step is explicit so you can see exactly what is happening.
    """
    # Step 1: Standardise
    mean  = X.mean(axis=0)
    std   = X.std(axis=0)
    X_std = (X - mean) / std

    # Step 2: Covariance matrix
    cov = np.cov(X_std.T)                       # shape (n_features, n_features)

    # Step 3: Eigendecomposition
    eigenvalues, eigenvectors = np.linalg.eigh(cov)
    # Sort descending
    idx          = np.argsort(eigenvalues)[::-1]
    eigenvalues  = eigenvalues[idx]
    eigenvectors = eigenvectors[:, idx]

    # Step 4: Project data onto top n_components eigenvectors
    components   = eigenvectors[:, :n_components]     # shape (n_features, n_components)
    X_projected  = X_std @ components                  # shape (n_samples, n_components)

    explained_var = eigenvalues[:n_components] / eigenvalues.sum()
    return X_projected, explained_var, components

X_pca, exp_var, components = pca_from_scratch(X, n_components=2)

print(f"Original data shape:   {X.shape}")
print(f"Projected data shape:  {X_pca.shape}")
print(f"Variance explained:    PC1={exp_var[0]*100:.1f}%  PC2={exp_var[1]*100:.1f}%")
print(f"Total captured:        {exp_var.sum()*100:.1f}% of original information")

# ── sklearn PCA — one line, same result ───────────────────────────────
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

pca = PCA(n_components=2)
X_pca_sk = pca.fit_transform(X_scaled)

print(f"\nsklearn PCA result:")
print(f"  Explained variance ratio: {pca.explained_variance_ratio_.round(3)}")
print(f"  Cumulative: {pca.explained_variance_ratio_.cumsum().round(3)}")

# ── How many components to keep? ──────────────────────────────────────
pca_full = PCA()
pca_full.fit(X_scaled)
cumulative = np.cumsum(pca_full.explained_variance_ratio_)

print(f"\nHow many components to keep (95% threshold):")
for k, cum in enumerate(cumulative, 1):
    bar = '█' * int(cum * 30)
    flag = ' ← keep this many' if cum >= 0.95 and cumulative[k-2] < 0.95 else ''
    print(f"  k={k}: {bar} {cum*100:.1f}%{flag}")

# ── Real use: compress and reconstruct data ───────────────────────────
k = 2   # use only 2 dimensions
pca2 = PCA(n_components=k)
X_compressed   = pca2.fit_transform(X_scaled)   # 1000×5 → 1000×2
X_reconstructed = pca2.inverse_transform(X_compressed)  # 1000×2 → 1000×5

# Reconstruction error — how much did we lose?
reconstruction_error = np.mean((X_scaled - X_reconstructed) ** 2)
print(f"\nCompression: 5 features → {k} components")
print(f"Reconstruction MSE: {reconstruction_error:.4f}")
print(f"Information preserved: {pca2.explained_variance_ratio_.sum()*100:.1f}%")

import numpy as np
from sklearn.cluster import SpectralClustering

# ── Simplified PageRank demo using eigenvalues ─────────────────────────
# Imagine a tiny web with 5 pages: [home, products, about, blog, cart]
# Entry [i][j] = 1 means page i links to page j

link_matrix = np.array([
    # home  prod  about  blog  cart
    [  0,    1,    1,    1,    0  ],   # home links to products, about, blog
    [  1,    0,    0,    1,    1  ],   # products links to home, blog, cart
    [  1,    0,    0,    0,    0  ],   # about links only to home
    [  1,    1,    0,    0,    0  ],   # blog links to home, products
    [  1,    1,    0,    0,    0  ],   # cart links to home, products
], dtype=float)

page_names = ['home', 'products', 'about', 'blog', 'cart']

# Build transition matrix: normalise each row so probabilities sum to 1
row_sums           = link_matrix.sum(axis=1, keepdims=True)
transition_matrix  = link_matrix / row_sums

# PageRank is the principal eigenvector of the TRANSPOSE of transition matrix
eigenvalues, eigenvectors = np.linalg.eig(transition_matrix.T)

# Find the eigenvector with eigenvalue closest to 1 (the stationary distribution)
idx        = np.argmin(np.abs(eigenvalues - 1.0))
page_rank  = np.real(eigenvectors[:, idx])
page_rank  = np.abs(page_rank) / np.abs(page_rank).sum()  # normalise to sum=1

print("PageRank scores (importance of each page):")
for name, score in sorted(zip(page_names, page_rank), key=lambda x: x[1], reverse=True):
    bar = '█' * int(score * 100)
    print(f"  {name:<12}: {bar} {score:.4f}")
print("\n'home' and 'products' score highest because many pages link to them.")

# ── Spectral clustering demo ───────────────────────────────────────────
np.random.seed(42)
# Create two concentric ring clusters — K-Means would fail on these
from sklearn.datasets import make_circles
X_circles, y_circles = make_circles(n_samples=300, noise=0.05, factor=0.4)

# K-Means fails on ring-shaped clusters
from sklearn.cluster import KMeans
km_labels = KMeans(n_clusters=2, random_state=42, n_init=10).fit_predict(X_circles)
km_accuracy = max(
    (km_labels == y_circles).mean(),
    (km_labels != y_circles).mean()
)

# Spectral clustering handles them correctly using eigenvectors
sc_labels = SpectralClustering(
    n_clusters=2, affinity='rbf', random_state=42
).fit_predict(X_circles)
sc_accuracy = max(
    (sc_labels == y_circles).mean(),
    (sc_labels != y_circles).mean()
)

print(f"\nRing-shaped clusters:")
print(f"  K-Means accuracy:         {km_accuracy:.1%}  ← fails, thinks linearly")
print(f"  Spectral Clustering:       {sc_accuracy:.1%}  ← succeeds via eigenvectors")

An RNN or a deep network trains fine for hours, then the loss spikes to NaN with no obvious cause in the data. One of the first things an experienced engineer checks: the eigenvalues of the recurrent weight matrix (or, for very deep feedforward nets, the spectral norm of consecutive weight matrices). If the largest eigenvalue's magnitude is consistently above 1, the network is multiplying activations by a factor greater than 1 at every time step or layer — the exact mechanism behind exploding gradients. If it's consistently below 1, you get vanishing gradients instead.

This is the same A·e = λ·e relationship from Section 2, just applied to a weight matrix instead of a covariance matrix — the eigenvalues tell you how much repeated multiplication by that matrix stretches or shrinks a signal, which is exactly what happens to gradients flowing backward through many time steps or many layers.

import numpy as np

np.random.seed(42)

# ── Diagnosing RNN instability via the eigenvalues of the recurrent weight ──
# This is the check an ML engineer runs when a recurrent model trains fine
# for a while, then the loss suddenly explodes to NaN

def check_recurrent_stability(W_hidden):
    """
    W_hidden: the hidden-to-hidden weight matrix in a vanilla RNN cell.
    Its eigenvalues determine whether repeated multiplication (one per
    time step) makes signals grow, shrink, or stay roughly stable.
    """
    eigenvalues = np.linalg.eigvals(W_hidden)
    magnitudes  = np.abs(eigenvalues)
    spectral_radius = magnitudes.max()   # largest |eigenvalue|

    print(f"Spectral radius (max |λ|): {spectral_radius:.4f}")
    if spectral_radius > 1.05:
        print("  ⚠ > 1 — repeated multiplication over time steps will EXPLODE")
        print("     fix: gradient clipping, or switch to LSTM/GRU gating")
    elif spectral_radius < 0.95:
        print("  ⚠ < 1 — repeated multiplication over time steps will VANISH")
        print("     fix: orthogonal initialisation, or switch to LSTM/GRU")
    else:
        print("  ✓ ≈ 1 — signal magnitude roughly preserved across time steps")
    return spectral_radius

# A poorly initialised recurrent weight — this is what causes the 2am page
W_unstable = np.random.randn(64, 64) * 0.5   # random scale, no care taken
check_recurrent_stability(W_unstable)

# A properly initialised recurrent weight — orthogonal initialisation
# guarantees every eigenvalue has magnitude exactly 1
from scipy.stats import ortho_group
W_orthogonal = ortho_group.rvs(64, random_state=42)
print()
check_recurrent_stability(W_orthogonal)

# This is exactly why orthogonal initialisation (nn.init.orthogonal_ in
# PyTorch) is a standard recommendation for vanilla RNNs — it directly
# controls the eigenvalues of the recurrent weight matrix at the start
# of training, before anything has had a chance to explode or vanish.

PCA is the most commonly taught example, but it's one of at least half a dozen production uses: spectral clustering (graph Laplacian eigenvectors), PageRank (principal eigenvector of a link matrix), recommender systems (matrix factorisation is built on eigen/singular decomposition), and stability analysis of recurrent and very deep networks, where the eigenvalues of a weight matrix directly predict exploding or vanishing gradients. Treating this as "the PCA math" makes every one of those other appearances look unrelated when they're the exact same underlying tool.

Some matrices — called defective matrices — do not have enough independent eigenvectors to form a full basis for their dimension, which breaks algorithms that assume diagonalisability. The reason covariance matrices are so convenient is the spectral theorem: every real symmetric matrix is guaranteed to have a full set of real eigenvalues and orthogonal eigenvectors. That guarantee is precisely why np.linalg.eigh (for symmetric matrices) is both safer and faster than the general-purpose np.linalg.eig — it can rely on a property that doesn't hold for arbitrary matrices.

Explained variance is a useful heuristic, not a universal rule. The right number of components depends on the downstream task: for visualisation you often want exactly 2 or 3 regardless of variance captured, and for anomaly detection the small-eigenvalue directions that a 95%-variance cutoff would discard are sometimes exactly where anomalies show up, since anomalies are often defined as points that violate the low-variance structure of normal data. Picking a threshold without checking it against actual downstream performance is a common way PCA preprocessing quietly hurts a model.

Eigenvectors are only defined up to sign and scale — if e is an eigenvector, so is −e and so is 5e, all satisfying A·e = λ·e equally well. Worse, when a matrix has a repeated eigenvalue, any vector in the entire subspace spanned by the eigenvectors for that eigenvalue is a valid eigenvector, so there is no single "correct" choice at all — just a convention. This is exactly why PCA results can flip sign between library versions or random seeds without being wrong; it's a documented consequence of the math, not a bug to chase down.

A full eigendecomposition of an n×n matrix costs roughly O(n³) — fine for a 50-feature covariance matrix, prohibitive for a 50,000-dimensional one. This is exactly why production systems almost never call np.linalg.eig directly on large matrices: they use randomized SVD, truncated SVD, or iterative methods like the power method or Lanczos algorithm, which approximate only the top few eigenvalues/eigenvectors needed without ever forming or decomposing the full matrix. sklearn's PCA switches to a randomized solver automatically once the data gets large enough.

Most vectors change direction when you multiply them by a matrix — the matrix rotates and stretches them at the same time. An eigenvector is one of the special directions a specific matrix doesn't rotate at all; it only stretches or shrinks it, by a factor called the eigenvalue. Every square matrix has its own small set of these special directions, and they describe the "natural axes" the matrix operates along — which is why finding them tells you something fundamental about what the matrix does to any input.

The covariance matrix encodes how every pair of features varies together. Its eigenvectors point in the directions along which the data spreads out the most — the directions of maximum variance — and its eigenvalues tell you exactly how much variance each of those directions captures. Sorting by eigenvalue and keeping the top few eigenvectors gives you a lower-dimensional representation that preserves as much of the original variation as possible for the number of dimensions you're willing to keep.

eigh assumes the input is symmetric (or Hermitian) and exploits that guarantee to return real eigenvalues, use a numerically more stable algorithm, and run faster. eig makes no such assumption and can return complex eigenvalues even from a matrix that's supposed to be symmetric but has tiny floating-point asymmetries. In practice: use eigh for covariance matrices, graph Laplacians, and anything you know is symmetric by construction; use eig only for genuinely asymmetric matrices, like a Markov transition matrix in PageRank.

PageRank models the web as a transition matrix — each entry is the probability of moving from one page to another by following a link. A page's long-run importance is the stationary distribution of a random walk over this graph, and that stationary distribution is exactly the eigenvector of the transition matrix with eigenvalue 1. Pages that many important pages link to accumulate more of that eigenvector's weight, which is why PageRank captures something closer to "importance" than simply counting incoming links.

An RNN applies the same recurrent weight matrix repeatedly, once per time step, so a signal (or a gradient flowing backward) gets multiplied by that matrix over and over. The eigenvalues of that matrix determine what happens under repeated multiplication: if the largest eigenvalue's magnitude is above 1, the signal grows exponentially with sequence length and eventually overflows to NaN; if it's below 1, the signal vanishes instead. I'd check the spectral radius of the recurrent weight directly, and consider orthogonal initialisation, gradient clipping, or switching to an LSTM/GRU, whose gating mechanisms are specifically designed to control this.