20. Linear Algebra for AI/ML (SVD, PCA, Tensors)

Level: Expert You'll use this for: dimensionality reduction, understanding embeddings, model compression (LoRA), and reading ML papers without glazing over.

Tensors and broadcasting

A tensor is the generalization of scalar (0D) → vector (1D) → matrix (2D) → n-D array. A batch of color images is a 4D tensor: [batch, height, width, channels]. Broadcasting is the rule set that lets numpy/PyTorch operate on tensors of different (but compatible) shapes without explicit loops.

import numpy as np
batch = np.random.rand(32, 128)      # 32 samples, 128 features each
bias = np.random.rand(128)           # one bias per feature
result = batch + bias                # broadcasts bias across all 32 rows

Misunderstanding broadcasting rules is one of the most common silent-bug sources in ML code — shapes that are "compatible" by broadcasting rules but not what you intended will run without error and produce wrong numbers.

Matrix rank and low-rank structure

The rank of a matrix is the number of linearly independent rows/columns — intuitively, the "true dimensionality" of the information it contains. A matrix can be 1000×1000 but have rank 10, meaning it's "secretly" much simpler than its size suggests.

Singular Value Decomposition (SVD)

Any matrix A (even non-square) can be decomposed:

A = U · Σ · Vᵀ

where Σ is a diagonal matrix of singular values, sorted largest to smallest, representing how much "information" each component carries. Keeping only the top k singular values gives the best possible rank-k approximation of A — the mathematical foundation of compression, denoising, and recommendation systems.

A = np.random.rand(100, 50)
U, S, Vt = np.linalg.svd(A, full_matrices=False)
k = 10
A_approx = U[:, :k] @ np.diag(S[:k]) @ Vt[:k, :]
print("compression: kept", k, "of", len(S), "singular values")

PCA — Principal Component Analysis, derived (not memorized)

PCA finds the directions (principal components) along which data varies the most. Mechanically: PCA = eigendecomposition of the covariance matrix, and its top components are the top singular vectors of the (centered) data matrix — i.e., PCA is SVD applied to centered data.

A 4-point PCA, by hand, to see each step concretely. Points: (0,0), (2,2), (4,4), (6,6) — all on a line, so PCA should find "the diagonal direction" as the one that captures all the variance.

1. Center the data (subtract the mean):
   mean = (3, 3)
   centered points: (-3,-3), (-1,-1), (1,1), (3,3)

2. Build the covariance matrix. For 2D centered data with columns x, y:
     cov[0][0] = mean(x·x) = (9+1+1+9)/4 = 5
     cov[1][1] = mean(y·y) = (9+1+1+9)/4 = 5
     cov[0][1] = cov[1][0] = mean(x·y) = (9+1+1+9)/4 = 5
   Cov = [5 5]
         [5 5]

3. Eigendecompose Cov (same characteristic-equation method as Ch. 11):
   det(Cov - λI) = (5-λ)² - 25 = 0  =>  λ² - 10λ = 0  =>  λ = 0 or λ = 10

   λ=10's eigenvector is the direction (1,1) (normalized: 1/√2, 1/√2) —
   this is the "diagonal" principal component, carrying ALL the variance
   (the other eigenvalue, 0, means zero variance perpendicular to it,
   which matches the data lying exactly on a line).

That's the entire mechanism np.linalg.eigh runs internally, just at whatever dimensionality your real dataset has — center, build the covariance matrix, and find its eigenvectors.

from numpy.linalg import eigh

def pca(X, k):
    X_centered = X - X.mean(axis=0)
    cov = np.cov(X_centered, rowvar=False)
    eigvals, eigvecs = eigh(cov)              # eigh: symmetric matrix, real eigenvalues
    order = np.argsort(eigvals)[::-1]
    top_k = eigvecs[:, order[:k]]
    return X_centered @ top_k                  # projected data

X = np.random.rand(200, 20)
reduced = pca(X, k=2)   # 20 dimensions -> 2, keeping max variance

LoRA (Low-Rank Adaptation) — this math in production

Fine-tuning a large model normally updates a full weight matrix W (potentially billions of parameters). LoRA observes that the update ΔW needed for fine-tuning often has low intrinsic rank, so instead of learning a full ΔW, you learn two small matrices B (n×r) and A (r×m) where r << n,m, and use ΔW = B·A. This is a direct, production-scale application of the low-rank idea from SVD above — it's why LoRA fine-tuning uses a tiny fraction of the memory/compute of full fine-tuning.

The arithmetic: one attention projection in a 7B-class model is 4096 × 4096.

full ΔW:       4096 · 4096          = 16,777,216 trainable numbers
LoRA, r = 8:   8 · (4096 + 4096)    =     65,536 trainable numbers   (0.39%)

Because B starts at zero, ΔW = 0 at step 0 and fine-tuning begins exactly at the pretrained model. That's a small detail of initialization that makes LoRA stable.

SVD, symbol by symbol, and by hand

A = U · Σ · Vᵀ          (A is m × n)

Read right to left, every matrix does three things: rotate (Vᵀ), stretch along the axes (Σ), rotate again (U). Any linear map, however messy, is a rotation, a stretch and another rotation.

By hand for A = [[3, 0], [4, 5]]:

AᵀA = [3 4]·[3 0] = [25 20]
      [0 5] [4 5]   [20 25]

eigenvalues of AᵀA:  (25-λ)² - 400 = 0  ->  λ = 45 or 5
singular values:     σ₁ = √45 ≈ 6.708,  σ₂ = √5 ≈ 2.236

(Same [[a,b],[b,a]] pattern as Ch. 11, with eigenvalues a ± b.)

Theorem (Eckart–Young, 1936). The best rank-k approximation of A, measured by the sum of squared errors (Frobenius norm), is the truncated SVD Aₖ, and its error is exactly the energy you dropped:

‖A - Aₖ‖²_F = σ²ₖ₊₁ + σ²ₖ₊₂ + ...

For our matrix, ‖A‖²_F = 9 + 0 + 16 + 25 = 50 = 45 + 5. Keeping only σ₁ (rank 1) leaves an error of exactly σ₂² = 5, so one number pair keeps 90% of the energy. Real data is far more extreme. The singular values of natural images, user–item ratings and LLM weight updates fall off a cliff, which is why compression, recommendation and LoRA work at all.

Attention, the equation behind every transformer

Attention(Q, K, V) = softmax( Q·Kᵀ / √d_k ) · V

Attention is a soft dictionary lookup. A Python dict returns the one value whose key matches exactly. Attention returns a blend of all values, weighted by how well each key matches.

A tiny example by hand. One query q = (1, 0), three tokens with keys (1, 0), (0, 1), (1, 1) and scalar values 10, 20, 90, with d_k = 2:

scores  = q·kᵢ        = (1, 0, 1)
scaled  = / √2        = (0.707, 0, 0.707)
weights = softmax     = (0.401, 0.198, 0.401)
output  = 0.401·10 + 0.198·20 + 0.401·90  ≈  44.1

The two keys that match the query get twice the weight of the one that doesn't.

Why divide by √d_k? The variance argument (Ch. 8 and 14 earning their keep). If the entries of q and k are independent with mean 0 and variance 1, then q·k = Σᵢ qᵢkᵢ is a sum of d_k terms each with variance 1, so Var(q·k) = d_k and its standard deviation is √d_k. With d_k = 64, raw scores have a standard deviation of about 8. Softmax of scores like (8, 0, -8) is (0.9997, 0.0003, 0.0000), which is nearly one-hot. In that regime the softmax gradient p - y (Ch. 13) is almost zero and learning stalls. Dividing by √d_k = 8 restores a standard deviation of 1, and the same scores become (0.67, 0.24, 0.09). The √d_k that appears in every transformer paper is a probability fact used to fix a calculus problem.

High-dimensional geometry: why embeddings work

Intuition from 2D and 3D fails badly in 768 dimensions. Draw two random unit vectors in ℝᵈ. Their cosine similarity has mean 0 and standard deviation about 1/√d:

d = 3:    std ≈ 0.58     random vectors often look "similar"
d = 768:  std ≈ 0.036    random vectors are almost exactly perpendicular

In high dimensions almost every pair of random directions is nearly orthogonal, so a space can hold an exponential number of nearly independent concepts. That's why a cosine similarity of 0.3 between two embeddings is a meaningful signal and not noise. It's also why approximate nearest-neighbour indexes (HNSW, IVF-PQ in FAISS and vector databases) can prune most of the space quickly.

Networking in the wild: PCA finds anomalies across a whole network

Stack a backbone network's link loads over time into a matrix: rows are 5-minute intervals, columns are links. Normal traffic is driven by a few shared forces (time of day, weekly cycle, a handful of big customers), so the matrix is close to low rank, and a few principal components explain most of the variance. Lakhina, Crovella and Diot (SIGCOMM 2004) split each time step into the part explained by the top components (the "normal subspace") and the residual. A DDoS, a routing change or an outage on a few links shows up as a residual spike, even when no single link crosses its own alert threshold. Network-wide anomaly detection is Eckart–Young applied to SNMP counters.

Real-life engineering ties

• AI: PCA for feature reduction/visualization, SVD for recommendation systems (matrix factorization: Users × Items ≈ low-rank approximation), LoRA/QLoRA for efficient fine-tuning, embedding compression.
• Security: PCA-based anomaly detection (project traffic/log features onto principal components; large residual = anomalous, doesn't fit the "normal" subspace).
• Network: less direct, but spectral graph analysis (eigenvectors of a network's Laplacian matrix) detects community structure/partitions.

More real-life examples:

• The 2006-2009 Netflix Prize competition was won using matrix-factorization techniques directly descended from SVD, applied to the giant sparse users x movies ratings matrix — arguably the event that popularized this math for recommendation systems industry-wide.
• Visualizing high-dimensional word or sentence embeddings (e.g., a 768-dimensional BERT embedding space) is typically done by reducing to 2-3 dimensions with PCA or the related t-SNE/UMAP techniques, so humans can look at a scatter plot.
• JPEG image compression uses the Discrete Cosine Transform, a close mathematical relative of SVD, to concentrate an image's information into a few significant coefficients and discard the rest — the same "keep the top-k components" idea as the SVD compression example above.
• Hugging Face's widely used PEFT library implements LoRA and QLoRA in production, letting engineers fine-tune multi-billion-parameter models on a single consumer GPU by training only the small low-rank matrices described in this chapter.

Common mistake

Running PCA without centering (subtracting the mean) first — PCA is mathematically defined on centered data; skipping this step gives a result dominated by the mean offset, not the actual variance structure.