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 Swiggy 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 rupees)
# 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 Swiggy 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")