Eigenvalue Decomposition for Statistics

Learning Objectives

By the end of this section, you will be able to:

Where You'll Apply This: Dimensionality reduction, feature extraction, data whitening, spectral clustering, recommender systems, image compression, and understanding the conditioning of optimization problems.

The Big Picture

Every linear transformation can be understood through its eigenvalues and eigenvectors. These special values reveal the intrinsic structure of a transformation — the directions along which it simply stretches or compresses, without rotating.

Historical Context

The eigenvalue problem has a rich history spanning centuries of mathematical development.

Eigenvalue Definition

Let's formally define eigenvalues and eigenvectors, then build intuition for what they mean.

This equation says: when we apply matrix A to eigenvector $\mathbf{v}$, the result is simply $\mathbf{v}$ scaled by $\lambda$. The direction is preserved (or flipped if $\lambda < 0$).

Geometric Interpretation

Think of a matrix as a transformation that stretches, rotates, and shears space. Most vectors get "scrambled" — they change both direction and length. But eigenvectors are invariant directions: they only stretch or shrink.

Interactive: 2D Eigenvalue Explorer

Explore how different 2×2 matrices transform the unit circle. The eigenvectors show the directions that only get scaled, and the eigenvalues tell you by how much.

Eigendecomposition

If a matrix has n linearly independent eigenvectors, we can write it in a special factorized form called the eigendecomposition.

This decomposition is powerful because it lets us understand any power of A easily:

The Spectral Theorem

For symmetric matrices (where $A = A^T$), something beautiful happens. The Spectral Theorem guarantees:

1. All eigenvalues are real — no complex numbers.
2. Eigenvectors are orthogonal — they form a perpendicular coordinate system.
3. Decomposition simplifies to $A = Q \Lambda Q^T$ where Q is orthogonal ($Q^{-1} = Q^T$).

Trace and Determinant Relationships

Two fundamental matrix properties are directly linked to eigenvalues:

Covariance Matrix Eigenanalysis

The covariance matrix is perhaps the most important matrix in statistics. Its eigendecomposition reveals the fundamental structure of multivariate data.

Key Properties of Covariance Matrices:

• Symmetric: $\Sigma = \Sigma^T$ (by construction)
• Positive Semi-Definite: All eigenvalues $\lambda_i \geq 0$
• Diagonal = Variances: $\Sigma_{ii} = \text{Var}(X_i)$
• Off-diagonal = Covariances: $\Sigma_{ij} = \text{Cov}(X_i, X_j)$

Interactive: Covariance Eigendecomposition

Generate correlated 2D data and watch how the eigendecomposition reveals the principal directions of variance. This is PCA in action!

The PCA Connection

Principal Component Analysis is nothing more than eigendecomposition of the covariance matrix! Each principal component is an eigenvector, and the variance it explains is the corresponding eigenvalue.

Computing Eigenvalues: Power Iteration

How do we actually compute eigenvalues? For large matrices, we can't solve the characteristic polynomial directly. The power iteration method is an elegant, iterative approach.

Why does this work? Express the initial vector in the eigenbasis:

After k iterations:

If $|\lambda_1| > |\lambda_2|$, the first term dominates as k → ∞.

Interactive: Power Iteration

Watch the power iteration algorithm converge to the dominant eigenvector. Notice how convergence speed depends on the ratio of eigenvalues.

Real-World Applications

AI/ML Connections

Eigenvalue decomposition appears throughout deep learning, often in subtle but crucial ways.

Python Implementation

Let's implement eigenvalue decomposition in Python. We'll cover both the NumPy approach and a from-scratch power iteration.

Covariance matrices are always symmetric, so always use eigh for PCA!

1import numpy as np
2
3# Robust eigenvalue decomposition
4def eigendecompose(matrix, descending=True):
5    """
6    Compute eigenvalues and eigenvectors.
7
8    For symmetric matrices, uses eigh (more stable).
9    For general matrices, uses eig.
10
11    Returns eigenvalues sorted in descending order.
12    """
13    # Convert to numpy array
14    A = np.array(matrix, dtype=np.float64)
15
16    # Check if symmetric
17    if np.allclose(A, A.T):
18        eigenvalues, eigenvectors = np.linalg.eigh(A)
19    else:
20        eigenvalues, eigenvectors = np.linalg.eig(A)
21
22    # Sort by eigenvalue (descending for PCA convention)
23    idx = np.argsort(eigenvalues)[::-1] if descending else np.argsort(eigenvalues)
24    eigenvalues = eigenvalues[idx]
25    eigenvectors = eigenvectors[:, idx]
26
27    return eigenvalues, eigenvectors

1import numpy as np
2from numpy.linalg import eigh, norm
3
4# ============================================
5# Example 1: Eigendecomposition of Covariance
6# ============================================
7np.random.seed(42)
8
9# Generate correlated 2D data
10n = 200
11mu = np.array([0, 0])
12cov_true = np.array([[2.0, 1.5],
13                      [1.5, 1.0]])
14
15# Sample from multivariate normal
16data = np.random.multivariate_normal(mu, cov_true, n)
17
18# Compute sample covariance
19cov_sample = np.cov(data.T)
20print("Sample Covariance Matrix:")
21print(cov_sample)
22
23# Eigendecomposition
24eigenvalues, eigenvectors = eigh(cov_sample)
25# Sort descending
26idx = np.argsort(eigenvalues)[::-1]
27eigenvalues = eigenvalues[idx]
28eigenvectors = eigenvectors[:, idx]
29
30print(f"\nEigenvalues: {eigenvalues}")
31print(f"Variance explained: {eigenvalues / eigenvalues.sum() * 100}%")
32print(f"\nPC1: {eigenvectors[:, 0]}")
33print(f"PC2: {eigenvectors[:, 1]}")
34
35# ============================================
36# Example 2: Power Iteration from Scratch
37# ============================================
38def power_iteration(A, num_iterations=100, tol=1e-10):
39    """
40    Compute dominant eigenvalue and eigenvector.
41
42    Parameters:
43        A: Square matrix
44        num_iterations: Maximum iterations
45        tol: Convergence tolerance
46
47    Returns:
48        eigenvalue, eigenvector
49    """
50    n = A.shape[0]
51    v = np.random.randn(n)
52    v = v / norm(v)
53
54    eigenvalue_old = 0
55
56    for i in range(num_iterations):
57        # Apply matrix
58        Av = A @ v
59
60        # Rayleigh quotient for eigenvalue estimate
61        eigenvalue = v @ Av
62
63        # Normalize
64        v = Av / norm(Av)
65
66        # Check convergence
67        if abs(eigenvalue - eigenvalue_old) < tol:
68            print(f"Converged after {i+1} iterations")
69            break
70
71        eigenvalue_old = eigenvalue
72
73    return eigenvalue, v
74
75# Test on our covariance matrix
76eigenvalue_power, eigenvector_power = power_iteration(cov_sample)
77print(f"\n=== Power Iteration Results ===")
78print(f"Dominant eigenvalue: {eigenvalue_power:.6f}")
79print(f"True eigenvalue:     {eigenvalues[0]:.6f}")
80print(f"Eigenvector: {eigenvector_power}")
81
82# ============================================
83# Example 3: PCA from Scratch
84# ============================================
85def pca(X, n_components=None):
86    """
87    Principal Component Analysis.
88
89    Parameters:
90        X: Data matrix (n_samples, n_features)
91        n_components: Number of components to keep
92
93    Returns:
94        projected_data, components, explained_variance_ratio
95    """
96    # Center the data
97    X_centered = X - X.mean(axis=0)
98
99    # Compute covariance matrix
100    cov = np.cov(X_centered.T)
101
102    # Eigendecomposition
103    eigenvalues, eigenvectors = eigh(cov)
104    idx = np.argsort(eigenvalues)[::-1]
105    eigenvalues = eigenvalues[idx]
106    eigenvectors = eigenvectors[:, idx]
107
108    # Select components
109    if n_components is None:
110        n_components = len(eigenvalues)
111
112    components = eigenvectors[:, :n_components]
113    explained_variance = eigenvalues[:n_components]
114    explained_variance_ratio = explained_variance / eigenvalues.sum()
115
116    # Project data
117    projected = X_centered @ components
118
119    return projected, components, explained_variance_ratio
120
121# Apply PCA
122projected, components, var_ratio = pca(data, n_components=2)
123print(f"\n=== PCA Results ===")
124print(f"Explained variance ratio: {var_ratio}")
125print(f"Total variance explained: {var_ratio.sum():.2%}")

Knowledge Check

Test your understanding of eigenvalue decomposition with this interactive quiz.

Summary

Looking Ahead: In the next section, we'll explore Singular Value Decomposition (SVD), which generalizes eigendecomposition to rectangular matrices and is even more widely used in machine learning applications.