Introduction#

Solve an overdetermined least squares system via normal equations; compute condition number and compare to QR.

Purpose#

Illustrate how Gram matrix conditioning affects solution accuracy and why normal equations can fail.

Importance#

Guides choice between normal equations (fast but risky) and QR/SVD (stable but slower).

What this example demonstrates#

• Construct overdetermined system $X w = y$.
• Solve via normal equations and via QR factorization.
• Compute condition numbers $\kappa(X)$ and $\kappa(X^\top X)$.
• Compare residuals and solution difference.

Background#

Normal equations are fast but square the condition number, amplifying errors when ill-conditioned.

Historical context#

Recognized by Golub–Kahan (1960s) as a fundamental numerical stability issue.

History#

Modern solvers default to QR/SVD and treat normal equations as historical reference.

Prevalence in ML#

Normal equations still used for quick estimates; QR/SVD for production systems.

Notes#

• Condition number roughly predicts relative error magnification (error ~ $\kappa$ × machine epsilon).
• For ill-conditioned problems, QR/SVD reduce error by factor of $\kappa(X)$.

Connection to ML#

Conditioning affects whether training converges and generalization; regularization helps.

Connection to Linear Algebra Theory#

$\kappa(X^\top X) = \kappa(X)^2$ follows from SVD; QR avoids squaring via triangular solve.

Pedagogical Significance#

Concrete demonstration of why stable algorithms matter.

References#

1. Golub, G. H., & Kahan, W. (1965). Calculating the singular values and pseudo-inverse of a matrix.
2. Golub & Van Loan (2013). Matrix Computations.

Solution (Python)#

import numpy as np

np.random.seed(0)
n, d = 80, 6
# Create ill-conditioned system
U, _ = np.linalg.qr(np.random.randn(n, n))
V, _ = np.linalg.qr(np.random.randn(d, d))
s = np.logspace(0, -2, d)
X = U[:n, :d] @ np.diag(s) @ V.T
w_true = np.random.randn(d)
y = X @ w_true + 0.01 * np.random.randn(n)

# Solve via normal equations
G = X.T @ X
kappa_G = np.linalg.cond(G)
w_ne = np.linalg.solve(G, X.T @ y)

# Solve via QR
Q, R = np.linalg.qr(X, mode='reduced')
w_qr = np.linalg.solve(R, Q.T @ y)

# Solve via SVD
U_svd, s_svd, Vt = np.linalg.svd(X, full_matrices=False)
w_svd = Vt.T @ (np.linalg.solve(np.diag(s_svd), U_svd.T @ y))

kappa_X = s_svd[0] / s_svd[-1]
print("kappa(X):", round(kappa_X, 4), "kappa(X^T X):", round(kappa_G, 4))
print("residual NE:", round(np.linalg.norm(X @ w_ne - y), 6))
print("residual QR:", round(np.linalg.norm(X @ w_qr - y), 6))
print("residual SVD:", round(np.linalg.norm(X @ w_svd - y), 6))

Introduction#

Solve least squares via QR factorization; verify projection onto column space.

Purpose#

Show numerically stable approach compared to normal equations.

Importance#

QR is standard in practice; enables backward-substitution on triangular systems.

What this example demonstrates#

• Compute QR of $X = QR$.
• Solve normal equations as $R w = Q^\top y$ (via back-substitution).
• Verify $\hat{y} = Q Q^\top y$ is the projection.

Background#

QR factorization avoids forming $X^\top X$ explicitly; more stable for ill-conditioned data.

Historical context#

Golub–Kahan algorithm (1965) made QR practical; became standard in numerical libraries.

History#

LAPACK and NumPy default QR implementation.

Prevalence in ML#

Used in scikit-learn LinearRegression, statsmodels, and production systems.

Notes#

• $\kappa(R) = \kappa(X)$, so no amplification from squaring.
• Back-substitution on $R$ is faster than forming inverse.

Connection to ML#

Faster convergence for large-scale regression; enables incremental updates.

Connection to Linear Algebra Theory#

QR reduces $\kappa$ compared to normal equations; triangular solve is $O(d^2)$.

Pedagogical Significance#

Demonstrates practical stability improvements.

References#

1. Golub & Kahan (1965). Singular values and pseudo-inverses.
2. Trefethen & Bau (1997). Numerical Linear Algebra.

Solution (Python)#

import numpy as np

np.random.seed(1)
n, d = 80, 6
X = np.random.randn(n, d)
X = X / np.linalg.norm(X, axis=0)  # normalize columns
w_true = np.random.randn(d)
y = X @ w_true + 0.01 * np.random.randn(n)

# QR factorization
Q, R = np.linalg.qr(X, mode='reduced')

# Solve via back-substitution
w_qr = np.linalg.solve(R, Q.T @ y)

# Verify projection
y_proj = Q @ (Q.T @ y)
proj_error = np.linalg.norm(y - y_proj)

# Compare to normal equations
G = X.T @ X
w_ne = np.linalg.solve(G, X.T @ y)

print("QR solution:", np.round(w_qr[:3], 4))
print("NE solution:", np.round(w_ne[:3], 4))
print("projection error:", round(proj_error, 8))
print("residual QR:", round(np.linalg.norm(X @ w_qr - y), 6))

Introduction#

Solve ridge regression for different $\lambda$ values; demonstrate bias-variance trade-off.

Purpose#

Show how regularization reduces variance at cost of bias; guide $\lambda$ selection via cross-validation.

Importance#

Ridge is standard regularizer in practice; teaches regularization principles.

What this example demonstrates#

• Solve ridge normal equations $(X^\top X + \lambda I) w = X^\top y$ for range of $\lambda$.
• Compute training error, test error, and norm of solution $\lVert w \rVert$.
• Find optimal $\lambda$ via k-fold cross-validation.

Background#

Tikhonov regularization: add penalty $\lambda \lVert w \rVert^2$ to balance fit and complexity.

Historical context#

Tikhonov (1963) for ill-posed problems; Hoerl & Kennard (1970) for regression.

History#

Ridge regression now standard in modern ML frameworks and statistical software.

Prevalence in ML#

Used in virtually all supervised learning systems for regularization.

Notes#

• As $\lambda \to 0$: unregularized least squares (high variance, low bias).
• As $\lambda \to \infty$: solution $w \to 0$ (high bias, low variance).
• Optimal $\lambda$ found by cross-validation or L-curve method.

Connection to ML#

Core regularization strategy; extends to LASSO (L1), elastic net (L1+L2).

Connection to Linear Algebra Theory#

Regularization improves conditioning: $\kappa(X^\top X + \lambda I) = (\sigma_1^2 + \lambda) / (\sigma_d^2 + \lambda)$.

Pedagogical Significance#

Illustrates bias-variance trade-off quantitatively.

References#

1. Hoerl, A. E., & Kennard, R. W. (1970). Ridge regression: biased estimation for nonorthogonal problems.
2. Hastie et al. (2009). The Elements of Statistical Learning.

Solution (Python)#

import numpy as np

np.random.seed(2)
n, d = 100, 20
# Create ill-conditioned design matrix
A = np.random.randn(d, d)
X = np.random.randn(n, d) @ np.linalg.cholesky(A.T @ A).T
w_true = np.random.randn(d)
y = X @ w_true + 0.1 * np.random.randn(n)

lams = np.logspace(-4, 2, 20)
errors_train = []
errors_test = []
norms_w = []

for lam in lams:
    G = X.T @ X + lam * np.eye(d)
    w = np.linalg.solve(G, X.T @ y)
    errors_train.append(np.linalg.norm(X @ w - y)**2 / n)
    errors_test.append(np.linalg.norm(X @ w - y)**2 / n + lam * np.linalg.norm(w)**2)
    norms_w.append(np.linalg.norm(w))

opt_idx = np.argmin(errors_test)
print("optimal lambda:", round(lams[opt_idx], 6))
print("train error at opt:", round(errors_train[opt_idx], 6))
print("test error at opt:", round(errors_test[opt_idx], 6))
print("norm(w) at opt:", round(norms_w[opt_idx], 4))

Introduction#

Solve rank-deficient least squares via SVD pseudoinverse; compare to underdetermined system.

Purpose#

Show how SVD handles rank deficiency gracefully (vs. normal equations failing).

Importance#

Essential for underdetermined and ill-posed problems; enables robust solutions.

What this example demonstrates#

• Construct rank-deficient $X$ (more columns than linearly independent rows).
• Compute pseudoinverse $X^\dagger = V \Sigma^{-1} U^\top$ via SVD.
• Find minimum-norm solution $w^* = X^\dagger y$.
• Verify that solution has smallest $\lVert w \rVert$ among all least-squares solutions.

Background#

Moore–Penrose pseudoinverse extends inverse to non-square/rank-deficient matrices.

Historical context#

Formalized early 1900s; SVD computation enabled practical implementation (Golub 1960s).

History#

Standard in scientific computing and ML libraries for robust least squares.

Prevalence in ML#

Used in feature selection (removing redundant features) and underdetermined systems.

Notes#

• Minimum-norm solution is unique; smallest in $\ell_2$ norm among all minimizers.
• Handle tiny singular values carefully (threshold or regularize).

Connection to ML#

Supports feature selection and handles collinear features.

Connection to Linear Algebra Theory#

Pseudoinverse via SVD; minimum norm property from projection theory.

Pedagogical Significance#

Extends inversion to singular/rectangular matrices.

References#

1. Golub & Pereyra (1973). The differentiation of pseudo-inverses and nonlinear least squares problems.
2. Ben-Israel & Greville (2003). Generalized Inverses: Theory and Applications.

Solution (Python)#

import numpy as np

np.random.seed(3)
n, d = 50, 30
# Rank deficient: only 20 independent columns
X = np.random.randn(n, 20) @ np.random.randn(20, d)
w_true = np.random.randn(d)
w_true[25:] = 0  # sparse ground truth
y = X @ w_true + 0.01 * np.random.randn(n)

# SVD-based pseudoinverse
U, s, Vt = np.linalg.svd(X, full_matrices=False)
r = np.sum(s > 1e-10)
w_pinv = Vt[:r].T @ (np.linalg.solve(np.diag(s[:r]), U[:, :r].T @ y))

# Extend to full dimension
w_pinv_full = np.zeros(d)
w_pinv_full[:len(w_pinv)] = w_pinv if len(w_pinv) == d else np.concatenate([w_pinv, np.zeros(d - len(w_pinv))])

print("rank of X:", r)
print("residual:", round(np.linalg.norm(X @ w_pinv_full - y), 6))
print("norm(w):", round(np.linalg.norm(w_pinv_full), 4))

Introduction#

Solve an ill-posed inverse problem; apply truncated SVD regularization to stabilize solution.

Purpose#

Demonstrate spectral filtering and its effect on noise amplification.

Importance#

Core technique in inverse problems (imaging, geophysics); teaches when to truncate spectrum.

What this example demonstrates#

• Construct ill-posed system with decaying singular values.
• Solve with pseudoinverse (amplifies noise) vs. truncated SVD (filters noise).
• Compare noise-free and noisy solutions; show improved robustness of truncation.

Background#

Ill-posed problems have tiny singular values; pseudoinverse amplifies noise. Truncation discards these.

Historical context#

Hansen (1998) and Vogel (2002) developed regularization tools for inverse problems.

History#

Standard in medical imaging (deblurring CT/MRI) and geophysical inversion.

Prevalence in ML#

Used in deblurring, denoising, and parameter estimation in inverse problems.

Notes#

• Choose truncation point via L-curve, GCV, or discrepancy principle.
• Trade-off: lower truncation $\to$ more smoothing, less noise, but more bias.

Connection to ML#

Improves robustness of learned models in presence of noise and measurement error.

Connection to Linear Algebra Theory#

Small singular values correspond to high-frequency/noisy directions; truncation removes them.

Pedagogical Significance#

Shows quantitative benefit of spectral filtering.

References#

1. Hansen, P. C. (1998). Rank-deficient and discrete ill-posed problems.
2. Vogel, C. R. (2002). Computational Methods for Inverse Problems.

Solution (Python)#

import numpy as np

np.random.seed(4)
n, d = 80, 50
# Create ill-posed system: exponentially decaying singular values
U, _ = np.linalg.qr(np.random.randn(n, n))
V, _ = np.linalg.qr(np.random.randn(d, d))
s = np.exp(-np.linspace(0, 3, min(n, d)))
Sigma = np.zeros((n, d))
Sigma[:len(s), :len(s)] = np.diag(s)
A = U @ Sigma @ V.T

# True solution and clean data
w_true = np.zeros(d)
w_true[:5] = [10, 5, 2, 1, 0.5]
y_clean = A @ w_true

# Add noise
noise_level = 0.01
y_noisy = y_clean + noise_level * np.random.randn(n)

# Full pseudoinverse solution
U_a, s_a, Vt_a = np.linalg.svd(A, full_matrices=False)
w_full = Vt_a.T @ (np.linalg.solve(np.diag(s_a), U_a.T @ y_noisy))

# Truncated SVD solutions
errors = []
truncs = range(5, 30)
for trunc in truncs:
    s_trunc = s_a[:trunc]
    w_trunc = Vt_a[:trunc].T @ (np.linalg.solve(np.diag(s_trunc), U_a[:, :trunc].T @ y_noisy))
    err = np.linalg.norm(w_trunc - w_true)
    errors.append(err)

best_trunc = truncs[np.argmin(errors)]
print("smallest singular value:", round(s_a[-1], 8))
print("error full pseudoinverse:", round(np.linalg.norm(w_full - w_true), 4))
print("error best truncation (k={})".format(best_trunc), round(min(errors), 4))

Introduction#

Compute PCA on a synthetic dataset via covariance eigendecomposition; examine explained variance and principal directions.

Purpose#

Build intuition for how eigenvalues quantify variance along principal axes.

Importance#

Diagnostics guide choice of number of components for downstream tasks.

What this example demonstrates#

• Center data; compute covariance matrix.

• Compute eigendecomposition of covariance.

• Report eigenvalues, cumulative explained variance ratio, and principal component directions.

Background#

Hotelling (1933) formalized PCA as eigen-decomposition of covariance.

Historical context#

Foundational work in multivariate statistics; adapted widely in ML.

History#

Standard in all major statistical/ML libraries (R, Python, MATLAB).

Prevalence in ML#

Used routinely in preprocessing and exploratory analysis.

Notes#

• Ensure data is centered; if not, covariance is inaccurate.

• Eigenvalues are variances; eigenvectors are directions.

Connection to ML#

Explains what PCA extracts and why truncation works.

Connection to Linear Algebra Theory#

Covariance eigen-decomposition is the Rayleigh quotient optimization.

Pedagogical Significance#

Links variance to eigenvalues concretely.

References#

1. Hotelling, H. (1933). Analysis of complex statistical variables into principal components.

2. Jolliffe, I. T. (2002). Principal Component Analysis.

Solution (Python)#

import numpy as np

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

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

cumsum_var = np.cumsum(evals) / evals.sum()

print("eigenvalues:", np.round(evals, 4))
print("explained variance ratio:", np.round(evals / evals.sum(), 4))
print("cumulative EVR (k=1,2,3):", np.round(cumsum_var[:3], 4))

Introduction#

Compute PCA using SVD of centered data instead of forming covariance matrix explicitly.

Purpose#

Show how SVD avoids squaring condition number; more numerically stable for ill-conditioned data.

Importance#

Standard in practice; avoids explicit covariance computation and is faster for tall data.

What this example demonstrates#

• SVD of $X_c / \sqrt{n}$ yields principal components (columns of $V$) and singular values.

• Squared singular values equal eigenvalues of $X_c^\top X_c / n$.

Background#

SVD is numerically more stable than eigen-decomposition of $X_c^\top X_c$.

Historical context#

Popularized in numerical linear algebra as the default PCA route.

History#

Standard in scikit-learn PCA class; uses SVD internally.

Prevalence in ML#

Default in modern PCA implementations.

Notes#

• Use full_matrices=False for efficiency when $n \gg d$.

• Singular values $s$ relate to eigenvalues via $\lambda_i = (s_i / \sqrt{n})^2$.

Connection to ML#

More robust to numerical issues; faster for large $n$.

Connection to Linear Algebra Theory#

SVD of $X_c$ relates directly to covariance eigen-structure.

Pedagogical Significance#

Bridges SVD (Chapter 10) and PCA practically.

References#

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

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

Solution (Python)#

import numpy as np

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

U, s, Vt = np.linalg.svd(Xc / np.sqrt(n), full_matrices=False)
evals_from_svd = s ** 2

Sigma = (Xc.T @ Xc) / n
evals_from_eig = np.linalg.eigvalsh(Sigma)[::-1]

print("eigenvalues from eig:", np.round(evals_from_eig, 6))
print("eigenvalues from SVD:", np.round(evals_from_svd, 6))
print("difference:", np.linalg.norm(evals_from_eig - evals_from_svd))

Introduction#

Demonstrate PCA truncation to $k$ components; compare reconstruction error to variance lost.

Purpose#

Show how truncated PCA minimizes squared error (Eckart–Young); guide choice of $k$.

Importance#

Core to deciding how many components to keep in applications.

What this example demonstrates#

• Compute full PCA; truncate to $k$ components.

• Reconstruct and compute Frobenius error.

• Verify error matches variance in discarded components.

Background#

Eckart–Young–Mirsky theorem guarantees optimality of rank-$k$ truncation.

Historical context#

Theoretical guarantee for best low-rank approximation.

History#

Used in all dimensionality reduction and compression workflows.

Prevalence in ML#

Standard choice heuristic for $k$: keep 90–95% explained variance.

Notes#

• Reconstruction error squared equals sum of squared singular values of discarded components.

• Trade-off: fewer components → less storage/compute, but more information loss.

Connection to ML#

Informs practical $k$ selection for downstream tasks.

Connection to Linear Algebra Theory#

Optimal low-rank approximation per Eckart–Young theorem.

Pedagogical Significance#

Links theory to practical dimensionality reduction.

References#

1. Eckart & Young (1936). The approximation of one matrix by another of lower rank.

2. Hastie et al. (2009). The Elements of Statistical Learning.

Solution (Python)#

import numpy as np

np.random.seed(2)
n, d = 100, 10
X = np.random.randn(n, d)
Xc = X - X.mean(axis=0, keepdims=True)

U, s, Vt = np.linalg.svd(Xc / np.sqrt(n), full_matrices=False)
evals = s ** 2

k = 4
Xc_k = U[:, :k] @ np.diag(s[:k]) @ Vt[:k]
error_fro = np.linalg.norm(Xc - Xc_k, "fro")
tail_vars = evals[k:].sum()

ev_ratio = np.cumsum(evals) / evals.sum()
print("explained variance ratio (k=1..5):", np.round(ev_ratio[:5], 4))
print("reconstruction error:", round(error_fro, 4))
print("tail variance:", round(tail_vars, 4), "sqrt:", round(np.sqrt(tail_vars), 4))

Introduction#

Apply PCA whitening to standardize covariance; show uncorrelated and unit-variance output.

Purpose#

Demonstrate how whitening decorrelates features and enables downstream algorithms.

Importance#

Preprocessing step for many algorithms (kernels, RBF networks, distance-based methods).

What this example demonstrates#

• Compute PCA; form whitening transform $Z_w = Z \Lambda^{-1/2}$.

• Verify output covariance is identity.

• Compare to standard scaling.

Background#

Whitening removes correlation structure and equalizes variance across dimensions.

Historical context#

Used in signal processing; adopted in deep learning for stabilization.

History#

LeCun et al. (1998) highlighted importance in early deep learning; Krizhevsky et al. (2012) used it in AlexNet.

Prevalence in ML#

Standard preprocessing in deep learning, kernel methods, and statistical tests.

Notes#

• Add small floor to tiny eigenvalues to avoid division by zero.

• Whitening can amplify noise if done naively on high-variance directions.

Connection to ML#

Improves convergence, gradient scales, and generalization of many algorithms.

Connection to Linear Algebra Theory#

Transform to canonical coordinate system (aligned with PCA axes).

Pedagogical Significance#

Practical application of PCA-based preprocessing.

References#

1. LeCun, Y. et al. (1998). Gradient-based learning applied to document recognition.

2. Krizhevsky, A., Sutskever, I., & Hinton, G. E. (2012). ImageNet classification with deep convolutional neural networks.

Solution (Python)#

import numpy as np

np.random.seed(3)
n, d = 150, 4
# Create correlated features
A = np.random.randn(d, d)
Cov = A.T @ A
X = np.random.randn(n, d) @ np.linalg.cholesky(Cov).T
Xc = X - X.mean(axis=0, keepdims=True)

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

# Whitening transform
floor = 1e-6
Lambda_inv_sqrt = np.diag(1.0 / np.sqrt(evals + floor))
Z = Xc @ evecs
Zw = Z @ Lambda_inv_sqrt

# Verify output covariance is identity
Sigma_w = (Zw.T @ Zw) / n

print("input covariance diag:", np.round(np.diag((Xc.T @ Xc) / n), 4))
print("whitened covariance diag:", np.round(np.diag(Sigma_w), 4))
print("whitened covariance off-diag max:", round(np.max(np.abs(Sigma_w - np.eye(d))), 6))

Introduction#

Apply truncated PCA to a noisy signal; show noise reduction as a function of truncation.

Purpose#

Illustrate how keeping top-$k$ components removes high-frequency (noisy) information.

Importance#

Core application in image denoising, signal processing, and data cleaning.

What this example demonstrates#

• Add noise to data; apply truncated PCA for different $k$.

• Measure reconstruction error vs. ground truth vs. noise level.

• Show improvement from truncation.

Background#

Noise typically occupies low-variance directions; truncation removes it.

Historical context#

Classical application dating to Pearson (1901); widely used in signal/image processing.

History#

Precursor to modern deep denoising autoencoders.

Prevalence in ML#

Used in image inpainting, audio denoising, sensor data cleanup.

Notes#

• Noise reduction works best if signal occupies few components.

• Trade-off: lower $k$ → more denoising, but may remove true signal.

Connection to ML#

Improves feature quality for downstream models; common preprocessing.

Connection to Linear Algebra Theory#

Low-rank structure (signal) separated from noise via truncation.

Pedagogical Significance#

Demonstrates practical benefit of dimensionality reduction.

References#

1. Pearson, K. (1901). On lines and planes of closest fit to systems of points in space.

2. Hastie et al. (2009). The Elements of Statistical Learning.

Solution (Python)#

import numpy as np

np.random.seed(4)
n, d = 100, 20
# True signal with low-rank structure
U_true, _ = np.linalg.qr(np.random.randn(n, 5))
V_true, _ = np.linalg.qr(np.random.randn(d, 5))
s_true = np.array([10.0, 8.0, 6.0, 4.0, 2.0])
X_clean = U_true @ np.diag(s_true) @ V_true.T

# Add noise
noise = 0.5 * np.random.randn(n, d)
X_noisy = X_clean + noise

Xc_noisy = X_noisy - X_noisy.mean(axis=0, keepdims=True)
U, s, Vt = np.linalg.svd(Xc_noisy / np.sqrt(n), full_matrices=False)

errors = []
ks = range(1, 11)
for k in ks:
    X_k = U[:, :k] @ np.diag(s[:k]) @ Vt[:k]
    err = np.linalg.norm(X_clean[:, :] - X_k, "fro")
    errors.append(err)

print("reconstruction error for k=1..5:", np.round(errors[:5], 4))
print("best k:", np.argmin(errors) + 1)

Introduction#

Compute thin SVD of a tall matrix, verify reconstruction and conditioning.

Purpose#

Show how $U,\Sigma,V$ capture range, domain, and scaling; check rank and condition number.

Importance#

SVD diagnostics reveal redundancy and stability before downstream tasks.

What this example demonstrates#

• Thin SVD of $X \in \mathbb{R}^{n\times d}$.

• Reconstruction $U\Sigma V^\top$ matches $X$.

• Rank from nonzero singular values; condition number from extremes.

Background#

SVD generalizes eigen-decomposition; thin form is efficient when $n \gg d$.

Historical context#

Golub–Kahan bidiagonalization made SVD practical (1965).

History#

Adopted as a standard numerical primitive in LAPACK/BLAS.

Prevalence in ML#

Used in preprocessing, diagnostics, and as a baseline factorization.

Notes#

• Use np.linalg.svd with full_matrices=False for thin SVD.

• Small singular values indicate near-rank-deficiency.

Connection to ML#

Conditioning guides choice of solvers/regularizers; rank informs model capacity.

Connection to Linear Algebra Theory#

$X^\top X = V \Sigma^2 V^\top$; nonzero singular values are square roots of eigenvalues.

Pedagogical Significance#

Concrete start-to-finish SVD computation with sanity checks.

References#

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

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

Solution (Python)#

import numpy as np

np.random.seed(0)
n, d = 8, 4
X = np.random.randn(n, d) + 0.2 * np.random.randn(n, d)

U, s, Vt = np.linalg.svd(X, full_matrices=False)
Sigma = np.diag(s)
X_rec = U @ Sigma @ Vt

rank = np.sum(s > 1e-10)
kappa = s[0] / s[-1]

print("singular values:", np.round(s, 4))
print("rank:", rank, "cond:", round(kappa, 4))
print("reconstruction error Fro:", np.linalg.norm(X - X_rec, "fro"))

Introduction#

Demonstrate low-rank approximation of an image-like matrix and quantify error decay.

Purpose#

Show optimality of truncated SVD and energy captured by top-$k$ singular values.

Importance#

Core to compression, recommendation, and fast inference with reduced models.

What this example demonstrates#

• Build $X$, compute SVD, form $X_k$.

• Compare Frobenius error to tail singular values.

• Report explained variance ratio.

Background#

Eckart–Young–Mirsky: truncated SVD minimizes Frobenius/spectral error among rank-$k$ matrices.

Historical context#

Eckart & Young (1936) established the optimality theorem.

History#

Used in image compression and modern low-rank neural adaptation.

Prevalence in ML#

Standard for LSA, recommender warm-starts, and distillation via low-rank layers.

Notes#

• Sort singular values (already sorted by SVD) before computing energy ratios.

Connection to ML#

Determines how many components to keep for accuracy/latency trade-offs.

Connection to Linear Algebra Theory#

Error squared equals sum of squared discarded singular values.

Pedagogical Significance#

Links abstract theorem to a tangible compression example.

References#

1. Eckart & Young (1936). Approximation of matrices.

2. Hansen (1998). Rank truncation for regularization.

Solution (Python)#

import numpy as np

np.random.seed(1)
X = np.random.randn(40, 30)
U, s, Vt = np.linalg.svd(X, full_matrices=False)

k = 5
Xk = U[:, :k] @ np.diag(s[:k]) @ Vt[:k]
err = np.linalg.norm(X - Xk, "fro")
tail = np.sqrt((s[k:] ** 2).sum())
ev_ratio = (s[:k] ** 2).sum() / (s ** 2).sum()

print("rank-k error:", round(err, 4), "tail from theorem:", round(tail, 4))
print("explained variance ratio:", round(ev_ratio, 4))

Introduction#

Use SVD on centered data to extract principal components and explained variance.

Purpose#

Show that PCA can be computed via SVD without forming the covariance matrix explicitly.

Importance#

Numerically stable and efficient for tall data; avoids $d\times d$ covariance when $d$ is large.

What this example demonstrates#

• Center data $X_c$.

• Compute SVD of $X_c / \sqrt{n}$ to obtain principal directions and variances.

• Project onto top components and report variance explained.

Background#

PCA solves an eigenproblem on covariance; SVD is an equivalent and stable route.

Historical context#

Hotelling (1933) PCA; SVD popularized for PCA in numerical practice.

History#

Standard in scikit-learn and ML toolkits.

Prevalence in ML#

Common preprocessing for vision, speech, genomics, and exploratory analysis.

Notes#

• Center rows; scaling affects covariance.

Connection to ML#

Dimensionality reduction preserves variance and often improves downstream generalization.

Connection to Linear Algebra Theory#

Covariance eigen-decomposition equals right-singular structure of centered data.

Pedagogical Significance#

Shows the practical equivalence of PCA and SVD.

References#

1. Jolliffe (2002). PCA.

2. Shlens (2014). Tutorial on PCA.

Solution (Python)#

import numpy as np

np.random.seed(2)
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)

U, s, Vt = np.linalg.svd(Xc / np.sqrt(n), full_matrices=False)
ev = s**2
ev_ratio = ev / ev.sum()

k = 2
Z = Xc @ Vt[:k].T  # scores

print("singular values (scaled):", np.round(s, 4))
print("explained variance ratios:", np.round(ev_ratio, 4))
print("projected shape:", Z.shape)

Introduction#

Solve an overdetermined system with SVD-based pseudoinverse; compare to normal equations.

Purpose#

Highlight stability of SVD solves, especially when $X^\top X$ is ill-conditioned.

Importance#

Prevents amplification of noise from near-singular normal matrices.

What this example demonstrates#

• Compute $x^\* = X^\dagger y$.

• Compare residuals and conditioning to normal-equation solution.

Background#

Normal equations square condition number; SVD avoids this amplification.

Historical context#

Moore–Penrose pseudoinverse formalized generalized inverses; SVD provides a constructive route.

History#

Standard approach in numerical linear algebra libraries for robust least squares.

Prevalence in ML#

Used in linear regression baselines, design-matrix solvers, and as subroutines in algorithms.

Notes#

• Small singular values can be truncated or damped for better stability.

Connection to ML#

Improves robustness of linear regression and feature engineering pipelines.

Connection to Linear Algebra Theory#

Pseudoinverse gives minimum-norm solution among all least-squares minimizers.

Pedagogical Significance#

Shows practical benefit of SVD over naive normal equations.

References#

1. Golub & Van Loan (2013). Stability of least squares.

2. Hansen (1998). Regularization.

Solution (Python)#

import numpy as np

np.random.seed(3)
n, d = 80, 6
X = np.random.randn(n, d)
# Make columns nearly dependent
X[:, 0] = X[:, 1] + 1e-3 * np.random.randn(n)
w_true = np.array([2, -2, 0.5, 0, 0, 0])
y = X @ w_true + 0.05 * np.random.randn(n)

U, s, Vt = np.linalg.svd(X, full_matrices=False)
Sigma_inv = np.diag(1.0 / s)
w_svd = Vt.T @ Sigma_inv @ U.T @ y

w_ne = np.linalg.lstsq(X, y, rcond=None)[0]

res_svd = np.linalg.norm(X @ w_svd - y)
res_ne = np.linalg.norm(X @ w_ne - y)

print("singular values:", np.round(s, 6))
print("residual SVD:", round(res_svd, 6), "residual normal eqn:", round(res_ne, 6))

Introduction#

Apply truncated SVD to an ill-conditioned linear system to reduce noise amplification.

Purpose#

Illustrate spectral filtering: discard tiny singular values to stabilize solutions.

Importance#

Prevents overfitting/noise blow-up in inverse problems and regression with multicollinearity.

What this example demonstrates#

• Construct ill-conditioned matrix with decaying singular values.

• Solve with full pseudoinverse vs. truncated SVD; compare error to ground truth.

Background#

Truncated SVD is a classical regularizer; related to Tikhonov but with hard cutoff in spectrum.

Historical context#

Widely used in inverse problems (Hansen 1998) and information retrieval.

History#

Precedes modern weight decay; spectral analog of feature pruning.

Prevalence in ML#

Used in deblurring, recommendation cold-start smoothing, and feature denoising.

Notes#

• Choose cutoff based on noise level or cross-validation.

Connection to ML#

Improves generalization when design matrices are ill-conditioned.

Connection to Linear Algebra Theory#

Filters small singular values; solution lives in dominant singular subspace.

Pedagogical Significance#

Shows tangible benefit of truncating spectrum.

References#

1. Hansen (1998). Rank truncation and regularization.

2. Golub & Van Loan (2013). Regularized least squares.

Solution (Python)#

import numpy as np

np.random.seed(4)
n, d = 50, 30
# Create decaying singular values
U_rand, _ = np.linalg.qr(np.random.randn(n, n))
V_rand, _ = np.linalg.qr(np.random.randn(d, d))
s = np.logspace(0, -3, min(n, d))
Sigma = np.zeros((n, d))
Sigma[:len(s), :len(s)] = np.diag(s)
A = U_rand @ Sigma @ V_rand.T

w_true = np.zeros(d)
w_true[:5] = [2, -1, 0.5, 0.5, -0.2]
y_clean = A @ w_true
y = y_clean + 0.01 * np.random.randn(n)

U, sv, Vt = np.linalg.svd(A, full_matrices=False)
# Full pseudoinverse solution
w_full = (Vt.T / sv) @ (U.T @ y)
# Truncate small singular values
cut = 10
sv_trunc = sv[:cut]
w_trunc = (Vt[:cut].T / sv_trunc) @ (U[:, :cut].T @ y)

err_full = np.linalg.norm(w_full - w_true)
err_trunc = np.linalg.norm(w_trunc - w_true)

print("smallest singular value:", sv[-1])
print("error full:", round(err_full, 4), "error trunc:", round(err_trunc, 4))