Introduction#

Compare PSD tests: eigenvalues vs. Cholesky; show failure on non-PSD and success on jittered fix.

Purpose#

Provide practical PSD diagnostics and repair.

Importance#

Prevents crashes/NaNs in kernel methods and GP solvers.

What this example demonstrates#

• Symmetrization, eigenvalue check, and Cholesky factorization.

• Adding jitter ($\epsilon I$) can restore PD for near-PSD matrices.

Background#

Cholesky is preferred for SPD solves; fails fast if not SPD.

Historical context#

Cholesky (early 1900s) for efficient SPD factorization.

Prevalence in ML#

Common in GP, KRR, and covariance handling.

Notes#

• Always symmetrize numerically; use small jitter (e.g., $1e-6$) when needed.

Connection to ML#

Stable training/inference for kernel models requires PSD kernels; jitter is standard.

Connection to Linear Algebra Theory#

Eigenvalue PSD characterization; $A=LL^\top$ iff $A$ is SPD.

Pedagogical Significance#

Hands-on PSD validation and fix.

References#

1. Golub & Van Loan (2013). Matrix Computations.

2. Rasmussen & Williams (2006). Gaussian Processes for ML.

Solution (Python)#

import numpy as np

np.random.seed(0)
A = np.random.randn(5, 5)
A = 0.5 * (A + A.T)  # symmetrize
A_bad = A.copy()
A_bad[0, 0] = -5.0  # make indefinite

def is_psd_eig(M, tol=1e-10):
    w = np.linalg.eigvalsh(M)
    return np.all(w >= -tol), w

for name, M in [("good", A), ("bad", A_bad)]:
    ok, w = is_psd_eig(M)
    print(f"{name}: min eig={w.min():.4f} PSD? {ok}")
    try:
        L = np.linalg.cholesky(M)
        print(f"{name}: Cholesky succeeded, diag min {np.min(np.diag(L)):.4f}")
    except np.linalg.LinAlgError:
        print(f"{name}: Cholesky failed")

# Jitter fix
eps = 1e-3
A_fix = A_bad + eps * np.eye(A_bad.shape[0])
ok, w = is_psd_eig(A_fix)
print("jittered: min eig=", round(w.min(), 4), "PSD?", ok)
L = np.linalg.cholesky(A_fix)
print("jittered Cholesky diag min:", round(np.min(np.diag(L)), 4))

Introduction#

Demonstrate PSD of covariance and perform whitening using eigenvalues; verify variance becomes identity.

Purpose#

Show practical PSD use in preprocessing and stability.

Importance#

Whitening decorrelates features; PSD guarantees nonnegative variances.

What this example demonstrates#

• $\Sigma=\tfrac{1}{n}X_c^\top X_c$ is PSD.

• Whitening transform $W=\Lambda^{-1/2} Q^\top$ yields covariance close to identity.

Background#

PCA/whitening are PSD-based transforms; eigenvalues must be nonnegative.

Historical context#

Classical in signal processing; ubiquitous in ML preprocessing.

Prevalence in ML#

Used in vision, speech, and as a step in ICA and some deep pipelines.

Notes#

• Add small floor to tiny eigenvalues for numerical stability.

Connection to ML#

Stabilizes downstream models; aligns scales across dimensions.

Connection to Linear Algebra Theory#

PSD eigen-structure; square roots and inverse square roots well-defined for PD.

Pedagogical Significance#

Concrete PSD-to-transform pipeline.

References#

1. Jolliffe (2002). PCA.

2. Shlens (2014). PCA tutorial.

Solution (Python)#

import numpy as np

np.random.seed(1)
n, d = 200, 5
X = np.random.randn(n, d) @ np.diag([3.0, 2.0, 1.0, 0.5, 0.2])
Xc = X - X.mean(axis=0, keepdims=True)

Sigma = (Xc.T @ Xc) / n
evals, evecs = np.linalg.eigh(Sigma)

# Whitening
floor = 1e-6
Lambda_inv_sqrt = np.diag(1.0 / np.sqrt(evals + floor))
W = Lambda_inv_sqrt @ evecs.T
Xw = (W @ Xc.T).T
Sigma_w = (Xw.T @ Xw) / n

print("Sigma PSD? min eig=", round(evals.min(), 6))
print("Whitened covariance diag:", np.round(np.diag(Sigma_w), 4))

Introduction#

Build an RBF kernel matrix, verify PSD via eigenvalues/Cholesky, and show how jitter fixes near-singularity.

Purpose#

Connect kernel validity to PSD checks used in SVM/GP/KRR implementations.

Importance#

Kernel methods rely on PSD to remain convex and numerically stable.

What this example demonstrates#

• RBF kernel is PSD; small datasets can be nearly singular when points are close.

• Jitter improves conditioning and enables Cholesky.

Background#

Mercer kernels generate PSD Gram matrices; RBF is a classic example.

Historical context#

Kernel trick popularized SVMs/GPs; jitter common in GP practice.

Prevalence in ML#

Standard in GPs, KRR, and kernel SVMs.

Notes#

• Condition number matters for solves; inspect spectrum.

Connection to ML#

Stable training/inference for kernel models.

Connection to Linear Algebra Theory#

Eigenvalue PSD test; Cholesky existence for SPD.

Pedagogical Significance#

Shows PSD verification and repair on a kernel matrix.

References#

1. Schölkopf & Smola (2002). Learning with Kernels.

2. Rasmussen & Williams (2006). Gaussian Processes for ML.

Solution (Python)#

import numpy as np

np.random.seed(2)
n, d = 12, 3
X = np.random.randn(n, d) * 0.2  # close points -> near-singular

def rbf_kernel(A, sigma=0.5):
    A2 = (A**2).sum(1)[:, None]
    D2 = A2 + A2.T - 2 * A @ A.T
    return np.exp(-D2 / (2 * sigma**2))

K = rbf_kernel(X)
evals = np.linalg.eigvalsh(K)
print("min eig:", round(evals.min(), 8), "cond:", float(evals.max() / max(evals.min(), 1e-12)))

try:
    L = np.linalg.cholesky(K)
    print("Cholesky ok, min diag:", round(np.min(np.diag(L)), 6))
except np.linalg.LinAlgError:
    print("Cholesky failed, adding jitter")
    eps = 1e-4
    K = K + eps * np.eye(n)
    L = np.linalg.cholesky(K)
    print("Cholesky ok after jitter, min diag:", round(np.min(np.diag(L)), 6))

Introduction#

Construct an SPD metric, compute Mahalanobis distances, and relate to whitening.

Purpose#

Show how SPD matrices define learned distances and why PSD/PD matters.

Importance#

Metric learning, anomaly detection, and clustering rely on valid metrics.

What this example demonstrates#

• $d_M(x,y)^2 = (x-y)^\top M (x-y)$ with $M\succ 0$; relate to Euclidean distance after whitening.

Background#

SPD metrics generalize Euclidean geometry; learned metrics often constrain $M\succeq 0$.

Historical context#

Mahalanobis (1936); modern metric learning enforces PSD.

Prevalence in ML#

Used in k-NN with learned metrics, anomaly scores, and Gaussian modeling.

Notes#

• Ensure $M$ is SPD (eigenvalues > 0); use Cholesky or eig clip.

Connection to ML#

Improves retrieval and clustering by reweighting feature space.

Connection to Linear Algebra Theory#

SPD defines inner products; whitening via $M^{1/2}$ links to Euclidean distance.

Pedagogical Significance#

Concrete example of PSD/PD defining geometry.

References#

1. De Maesschalck et al. (2000). The Mahalanobis distance.

2. Murphy (2022). Probabilistic Machine Learning.

Solution (Python)#

import numpy as np

np.random.seed(3)
d = 4
A = np.random.randn(d, d)
M = A.T @ A + 0.5 * np.eye(d)  # SPD

x = np.random.randn(d)
y = np.random.randn(d)
diff = x - y
mah2 = diff.T @ M @ diff
eucl2 = np.dot(diff, diff)

evals = np.linalg.eigvalsh(M)
print("M min eig:", round(evals.min(), 6))
print("Mahalanobis^2:", round(mah2, 6), " Euclidean^2:", round(eucl2, 6))

Introduction#

Solve KRR with an SPD kernel matrix using Cholesky; show stability vs. explicit inverse.

Purpose#

Demonstrate practical SPD solve in a common ML method.

Importance#

Avoids numerical issues and is standard in GP/KRR implementations.

What this example demonstrates#

• Solve $(K+\lambda I)\alpha = y$ via Cholesky; predictions $\hat{y}=K\alpha$.

Background#

Ridge regularization makes $K+\lambda I$ SPD.

Historical context#

Kernel methods mainstream since 1990s; Cholesky is the standard linear solve.

Prevalence in ML#

Regression/classification with kernels; GP regression uses same linear system.

Notes#

• Regularization size affects conditioning; add jitter if needed.

Connection to ML#

Core kernel regression solve; identical algebra to GP posterior mean.

Connection to Linear Algebra Theory#

SPD guarantees unique solution and stable Cholesky factorization.

Pedagogical Significance#

Shows end-to-end use of SPD property in a solver.

References#

1. Schölkopf & Smola (2002). Learning with Kernels.

2. Rasmussen & Williams (2006). Gaussian Processes for ML.

3. Murphy (2022). Probabilistic Machine Learning.

Solution (Python)#

import numpy as np

np.random.seed(4)
n, d = 40, 3
X = np.random.randn(n, d)

def rbf_kernel(A, B=None, sigma=1.0):
    if B is None:
        B = A
    A2 = (A**2).sum(1)[:, None]
    B2 = (B**2).sum(1)[None, :]
    D2 = A2 + B2 - 2 * A @ B.T
    return np.exp(-D2 / (2 * sigma**2))

K = rbf_kernel(X, X, sigma=1.2)
w_true = np.random.randn(n)
y = K @ w_true + 0.05 * np.random.randn(n)

lam = 1e-2
K_reg = K + lam * np.eye(n)
L = np.linalg.cholesky(K_reg)

# Solve via Cholesky: L L^T alpha = y
z = np.linalg.solve(L, y)
alpha = np.linalg.solve(L.T, z)
y_hat = K @ alpha

rmse = np.sqrt(np.mean((y_hat - y)**2))
print("Cholesky solve RMSE:", round(rmse, 6))
print("SPD? min eig(K_reg)=", round(np.linalg.eigvalsh(K_reg).min(), 6))

Introduction#

Compute the largest eigenvalue/vector of an SPD matrix via power iteration and compare to numpy’s eigvals.

Purpose#

Illustrate convergence rate and normalization; give a practical recipe.

Importance#

Many large-scale methods rely on top eigenpairs (PCA, spectral norm estimation).

What this example demonstrates#

• Convergence rate depends on eigengap $|\lambda_1/\lambda_2|$.

• Rayleigh quotient of iterates approximates $\lambda_1$.

Background#

Power methods date back a century and remain relevant for scalable eigensolvers.

Historical context#

von Mises/Householder refinements; modern variants include Lanczos/Arnoldi.

Prevalence in ML#

Used implicitly in randomized PCA, spectral norm estimation, and operator norm regularization.

Notes#

• Normalize every step; monitor residual $\lVert Ax-\lambda x\rVert$.

Connection to ML#

Top eigenvalue relates to Lipschitz constants; top eigenvector for PCA directions.

Connection to Linear Algebra Theory#

Convergence proof via eigen-expansion; Rayleigh quotient bounds.

Pedagogical Significance#

Shows a simple iterative algorithm achieving an eigenpair without full decomposition.

References#

1. Trefethen & Bau (1997). Numerical Linear Algebra.

2. Golub & Van Loan (2013). Matrix Computations.

Solution (Python)#

import numpy as np

np.random.seed(0)
d = 6
A = np.random.randn(d, d)
A = A.T @ A + 0.5 * np.eye(d)  # SPD

def power_iteration(A, iters=50):
		x = np.random.randn(A.shape[0])
		x /= np.linalg.norm(x)
		for _ in range(iters):
				x = A @ x
				x /= np.linalg.norm(x)
		lam = x @ (A @ x)
		return lam, x

lam_pi, v_pi = power_iteration(A, iters=40)
eigvals, eigvecs = np.linalg.eigh(A)
lam_true = eigvals[-1]
print("power iteration λ≈", round(lam_pi, 6), " true λ=", round(lam_true, 6))
print("angle to true v (deg):", round(np.degrees(np.arccos(np.clip(abs(v_pi @ eigvecs[:, -1]), -1, 1))), 6))

Introduction#

Compute covariance, eigenpairs, and project data onto top-$k$ components; verify variance captured.

Purpose#

Connect PCA steps to eigenvalues/eigenvectors and retained variance.

Importance#

Core dimensionality reduction tool in ML.

What this example demonstrates#

• $\Sigma = \tfrac{1}{n} X_c^\top X_c$ is PSD; eigenvalues = variances along eigenvectors.

• Retained variance ratio from top-$k$ eigenvalues.

Background#

Eigen-decomposition of covariance equals SVD-based PCA.

Historical context#

PCA roots in Pearson/Hotelling; ubiquitous in data analysis.

Prevalence in ML#

Widely used in preprocessing, visualization, and compression.

Notes#

• Center data; consider scaling features.

Connection to ML#

Variance retention guides choice of $k$; whitening uses inverse sqrt of eigenvalues.

Connection to Linear Algebra Theory#

PSD eigen-structure; orthogonal projections onto principal subspaces.

Pedagogical Significance#

Demonstrates PSD eigen-decomposition in a practical workflow.

References#

1. Jolliffe (2002). PCA.

2. Shlens (2014). PCA tutorial.

Solution (Python)#

import numpy as np

np.random.seed(1)
n, d, k = 120, 8, 3
X = np.random.randn(n, d) @ np.diag(np.linspace(3, 0.5, d))
Xc = X - X.mean(axis=0, keepdims=True)

Sigma = (Xc.T @ Xc) / n
evals, evecs = np.linalg.eigh(Sigma)
idx = np.argsort(evals)[::-1]
evals, evecs = evals[idx], evecs[:, idx]

Vk = evecs[:, :k]
X_proj = Xc @ Vk
variance_retained = evals[:k].sum() / evals.sum()
print("Top-k variance retained:", round(variance_retained, 4))

Introduction#

Compute PageRank on a small directed graph using power iteration on a stochastic matrix with damping.

Purpose#

Show Perron–Frobenius in action and convergence to the dominant eigenvector.

Importance#

Seminal large-scale eigen-application; template for Markov-chain ranking.

What this example demonstrates#

• Transition matrix $P$ has eigenvalue 1; power iteration converges to stationary distribution.

• Damping ensures irreducibility/aperiodicity.

Background#

Web graph ranking; damping factor (e.g., 0.85) handles dead ends/spiders.

Historical context#

Brin–Page (1998) launched web search revolution.

Prevalence in ML#

Graph ranking, recommendation propagation, random-walk-based features.

Notes#

• Normalize columns to sum to 1; add teleportation for damping.

Connection to ML#

Graph-based semi-supervised learning often reuses random-walk ideas.

Connection to Linear Algebra Theory#

Perron–Frobenius guarantees positive dominant eigenvector; spectral gap drives convergence.

Pedagogical Significance#

Concrete, small-scale power iteration on a stochastic matrix.

References#

1. Brin, S., & Page, L. (1998). The anatomy of a large-scale hypertextual web search engine.

2. Langville, A., & Meyer, C. (2006). Google’s PageRank and Beyond.

Solution (Python)#

import numpy as np

# Small directed graph adjacency
A = np.array([[0,1,1,0],
							[1,0,0,1],
							[1,1,0,0],
							[0,1,1,0]], dtype=float)

# Column-stochastic transition (out-links per column)
col_sums = A.sum(axis=0, keepdims=True)
P = A / np.where(col_sums == 0, 1, col_sums)

alpha = 0.85
n = P.shape[0]
J = np.ones((n, n)) / n
M = alpha * P + (1 - alpha) * J  # damping

v = np.ones(n) / n
for _ in range(50):
		v = M @ v
		v = v / v.sum()

eigvals, eigvecs = np.linalg.eig(M)
idx = np.argmax(np.real(eigvals))
v_true = np.real(eigvecs[:, idx])
v_true = v_true / v_true.sum()

print("Power iteration PageRank:", np.round(v, 4))
print("Eigenvector PageRank:", np.round(v_true, 4))

Introduction#

Compute Gershgorin disks for a matrix and compare to actual eigenvalues to illustrate spectrum localization.

Purpose#

Provide quick sanity bounds on eigenvalues without full eigendecomposition.

Importance#

Useful for diagnosing stability (e.g., Jacobians) and conditioning.

What this example demonstrates#

• All eigenvalues lie within the union of Gershgorin disks centered at $a_{ii}$ with radius $\sum_{j\neq i} |a_{ij}|$.

Background#

Classic theorem (1911) for eigenvalue localization.

Historical context#

Still a staple for quick qualitative checks.

Prevalence in ML#

Less direct, but useful for reasoning about spectrum without expensive computation.

Notes#

• Tightness varies; row/column scaling can sharpen disks.

Connection to ML#

Stability analysis of iterative methods; rough spectral norm estimates.

Connection to Linear Algebra Theory#

Eigenvalue inclusion sets; diagonal dominance implications.

Pedagogical Significance#

Gives a geometric picture of eigenvalue bounds.

References#

1. Gershgorin, S. (1911). Über die Abgrenzung der Eigenwerte einer Matrix.

2. Horn & Johnson (2013). Matrix Analysis.

Solution (Python)#

import numpy as np

np.random.seed(2)
A = np.random.randn(4, 4)

centers = np.diag(A)
radii = np.sum(np.abs(A), axis=1) - np.abs(centers)
eigvals = np.linalg.eigvals(A)

print("Gershgorin centers:", np.round(centers, 3))
print("Gershgorin radii:", np.round(radii, 3))
print("Eigenvalues:", np.round(eigvals, 3))

Introduction#

Perform unnormalized spectral clustering on a simple graph with two clusters; use second-smallest eigenvector (Fiedler) for separation.

Purpose#

Show how Laplacian eigenvectors reveal cluster structure.

Importance#

Nonlinear/non-convex cluster discovery.

What this example demonstrates#

• $L=D-W$; smallest eigenvalue 0 with eigenvector $\mathbf{1}$; Fiedler vector splits the graph.

Background#

Graph cuts and Laplacian spectra; normalized variants common in practice.

Historical context#

Spectral clustering surged in the 2000s for manifold data.

Prevalence in ML#

Image segmentation, manifold learning, community detection.

Notes#

• For normalized Laplacian, use $L_{\text{sym}}$; results similar on this toy example.

Connection to ML#

Embedding nodes using a few eigenvectors before k-means is standard pipeline.

Connection to Linear Algebra Theory#

Properties of Laplacian eigenvalues (nonnegative; multiplicity of 0 equals number of components).

Pedagogical Significance#

Concrete end-to-end spectral clustering demonstration.

References#

1. Ng, A., Jordan, M., & Weiss, Y. (2002). On Spectral Clustering.

2. von Luxburg, U. (2007). A tutorial on spectral clustering.

Solution (Python)#

import numpy as np

# Two clusters with strong intra-cluster edges
W = np.array([
		[0,1,1,0,0,0],
		[1,0,1,0,0,0],
		[1,1,0,0,0,0],
		[0,0,0,0,1,1],
		[0,0,0,1,0,1],
		[0,0,0,1,1,0],
], dtype=float)

D = np.diag(W.sum(axis=1))
L = D - W
evals, evecs = np.linalg.eigh(L)

fiedler = evecs[:, 1]
clusters = (fiedler > 0).astype(int)
print("Eigenvalues:", np.round(evals, 4))
print("Fiedler vector:", np.round(fiedler, 4))
print("Cluster assignment via sign:", clusters)

Introduction#

Show how null space reveals linear dependencies among features and why OLS becomes non-unique when $\operatorname{rank}(X)<d$.

Purpose#

Compute null space vectors and connect them to redundant directions; use pseudoinverse for a minimal-norm solution.

Importance#

Avoids unstable fits and clarifies identifiability in models.

What this example demonstrates#

• If $v\in\text{null}(X)$, $X(w+\alpha v) = Xw$ for all $\alpha$; infinitely many OLS solutions.

• Pseudoinverse $w^* = X^+ y$ yields the minimal-norm solution.

Background#

Rank deficiency arises from duplicate/derived features or insufficient data.

Historical context#

Gauss/Legendre least squares; Penrose pseudoinverse enables solutions in singular cases.

Prevalence in ML#

High-dimensional regression, feature engineering pipelines, polynomial expansions.

Notes#

• Use SVD to diagnose numerical rank; add Ridge to regularize.

Connection to ML#

Feature selection and regularization strategies hinge on rank awareness.

Connection to Linear Algebra Theory#

FTLA: residuals in $\text{null}(X^\top)$; solution set $w_0 + \text{null}(X)$.

Pedagogical Significance#

Makes the geometry of “non-unique solutions” tangible.

References#

1. Golub & Van Loan (2013). Matrix Computations.

2. Penrose (1955). Moore–Penrose pseudoinverse.

3. Hoerl & Kennard (1970). Ridge regression.

Solution (Python)#

import numpy as np

np.random.seed(0)
n, d = 20, 6
X = np.random.randn(n, d)
X[:, 5] = X[:, 0] + X[:, 1]  # make a perfectly colinear feature
w_true = np.array([1.0, -0.5, 0.3, 0.0, 2.0, 0.3])
y = X @ w_true + 0.1 * np.random.randn(n)

U, S, Vt = np.linalg.svd(X, full_matrices=False)
rank = np.sum(S > 1e-8)
nullspace_basis = Vt[rank:].T  # columns spanning null(X)

print("rank(X)=", rank, " d=", d, " nullity=", d - rank)
print("Nullspace basis shape:", nullspace_basis.shape)

# Minimal-norm solution via pseudoinverse
w_min = Vt.T @ (np.where(S > 1e-12, (U.T @ y) / S, 0.0))
print("||w_min||2=", np.linalg.norm(w_min))
print("OLS residual norm:", np.linalg.norm(y - X @ w_min))

Introduction#

Verify empirically that centered covariance has rank at most $n-1$ regardless of feature dimension.

Purpose#

Explain why PCA cannot produce more than $n-1$ nonzero eigenvalues and how this affects high-dimensional settings.

Importance#

Shapes expectations for PCA on small data; prevents overinterpretation.

What this example demonstrates#

• With $X_c\in\mathbb{R}^{n\times d}$ centered, $\operatorname{rank}(X_c) \le \min(n-1, d)$; hence $\operatorname{rank}(\Sigma) \le n-1$.

Background#

Centering imposes a linear constraint across rows, reducing rank by at least one when $n>1$.

Historical context#

PCA theory and practice emphasize centering for correct variance structure.

Prevalence in ML#

Common in text, genomics, and other $d\gg n$ problems.

Notes#

• Always center before PCA; whitening depends on accurate rank.

Connection to ML#

Model selection of $k$ principal components must respect $n-1$ limit.

Connection to Linear Algebra Theory#

Row-sum constraint places $\mathbf{1}$ in $\text{null}(X_c^\top)$.

Pedagogical Significance#

Reinforces how constraints reduce rank.

References#

1. Jolliffe (2002). PCA.

2. Shlens (2014). PCA tutorial.

Solution (Python)#

import numpy as np

np.random.seed(1)
n, d = 30, 200
X = np.random.randn(n, d)
Xc = X - X.mean(axis=0, keepdims=True)

U, S, Vt = np.linalg.svd(Xc, full_matrices=False)
rank = np.sum(S > 1e-10)
print("rank(Xc)=", rank, " <= min(n-1,d)=", min(n-1, d))

Introduction#

Construct a synthetic user–item rating matrix with known low rank and recover it via truncated SVD.

Purpose#

Demonstrate latent-factor modeling and show reconstruction error scales with tail singular values.

Importance#

Illustrates the power of rank reduction in recommender systems.

What this example demonstrates#

• $R\approx U_k \Sigma_k V_k^\top$ captures most variance when spectrum decays.

Background#

Matrix factorization underlies collaborative filtering; ALS/SGD optimize latent vectors.

Historical context#

Post-Netflix Prize, low-rank methods became industry standard.

Prevalence in ML#

Ubiquitous in recommendation and implicit feedback modeling.

Notes#

• For missing data, completion requires specialized optimization (not shown here).

Connection to ML#

Latent dimensions reflect user/item factors; rank controls capacity.

Connection to Linear Algebra Theory#

Eckart–Young guarantees best rank-$k$ approximation.

Pedagogical Significance#

Shows direct link from SVD to practical factor models.

References#

1. Koren, Bell, Volinsky (2009). Matrix factorization techniques for recommender systems.

2. Candès & Recht (2009). Exact matrix completion via convex optimization.

Solution (Python)#

import numpy as np

np.random.seed(2)
u, i, k = 80, 60, 5
U_true = np.random.randn(u, k)
V_true = np.random.randn(i, k)
R = U_true @ V_true.T + 0.1 * np.random.randn(u, i)

U, S, Vt = np.linalg.svd(R, full_matrices=False)
Rk = (U[:, :k] * S[:k]) @ Vt[:k]
err = np.linalg.norm(R - Rk, 'fro')**2
tail = (S[k:]**2).sum()
print("Fro error:", round(err, 6), " Tail sum:", round(tail, 6), " Close?", np.allclose(err, tail, atol=1e-5))

Introduction#

Solve $Ax=b$ when $A$ is rectangular or rank-deficient; verify $AA^+A=A$ and minimal norm among all solutions.

Purpose#

Provide a robust recipe for under-/overdetermined systems.

Importance#

Avoids fragile inverses and clarifies the solution geometry.

What this example demonstrates#

• $x^*=A^+ b$ minimizes $\lVert x\rVert_2$ subject to $Ax=b$ for consistent systems.

• Penrose conditions hold numerically.

Background#

Pseudoinverse defined via SVD; used in control, signal processing, ML.

Historical context#

Penrose (1955) established the four defining equations.

Prevalence in ML#

Closed-form layers, analytic baselines, and data-fitting routines.

Notes#

• Use SVD-backed implementations; threshold small singular values.

Connection to ML#

Stable baselines and analytic steps inside pipelines.

Connection to Linear Algebra Theory#

Projects onto $\text{row}(A)$/$\text{col}(A)$; minimal-norm in $\text{null}(A)$ components.

Pedagogical Significance#

Bridges algebraic definition to numerical practice.

References#

1. Penrose, R. (1955). A generalized inverse for matrices.

2. Golub & Van Loan (2013). Matrix Computations.

Solution (Python)#

import numpy as np

np.random.seed(3)
m, d = 10, 12
A = np.random.randn(m, d)
# Make A rank-deficient by zeroing a singular value via colinearity
A[:, 0] = A[:, 1] + A[:, 2]
b = np.random.randn(m)

U, S, Vt = np.linalg.svd(A, full_matrices=False)
S_inv = np.where(S > 1e-10, 1.0 / S, 0.0)
A_plus = Vt.T @ np.diag(S_inv) @ U.T

x_star = A_plus @ b
print("Penrose A A^+ A ~ A?", np.allclose(A @ (A_plus @ A), A, atol=1e-8))
print("Residual ||Ax-b||:", np.linalg.norm(A @ x_star - b))
print("||x_star||2:", np.linalg.norm(x_star))

Introduction#

Show that the attention score matrix $S=QK^\top$ has rank at most $\min(n, d_k)$ and explore implications for head dimension.

Purpose#

Connect feature dimension to expressivity through rank bounds.

Importance#

Head size choices affect the diversity of attention patterns.

What this example demonstrates#

• For $Q\in\mathbb{R}^{n\times d_k}$ and $K\in\mathbb{R}^{n\times d_k}$, $\operatorname{rank}(QK^\top) \le \min(n, d_k)$.

Background#

Rank of product bounded by inner dimension; scaled dot-products preserve rank.

Historical context#

Transformers leverage multiple heads to increase effective rank/expressivity.

Prevalence in ML#

All transformer models; multi-head concatenation increases representational capacity.

Notes#

• Multi-head attention can be seen as block structures that raise overall rank after concatenation.

Connection to ML#

Guides architecture design (choosing $d_k$ and number of heads).

Connection to Linear Algebra Theory#

Rank bounds and product properties.

Pedagogical Significance#

Links a practical hyperparameter to a crisp linear algebra bound.

References#

1. Vaswani, A. et al. (2017). Attention Is All You Need.

2. Devlin, J. et al. (2019). BERT.

Solution (Python)#

import numpy as np

np.random.seed(4)
for n, dk in [(32, 8), (32, 32), (64, 16)]:
	 Q = np.random.randn(n, dk)
	 K = np.random.randn(n, dk)
	 S = Q @ K.T
	 r = np.linalg.matrix_rank(S)
	 print(f"n={n}, d_k={dk}, rank(S)={r}, bound={min(n, dk)}")