Applied and computational mathematics

Mathematics of machine learning

The linear algebra, optimization and probability that power algorithms learning patterns from data.

IntuitionLearning as descending a landscape

Suppose you want to predict something — a house price, a spam label, the next word in a sentence — using a model with adjustable knobs $\theta$ . Training the model means choosing $\theta$ so that a chosen error measure $L(\theta)$ , the loss, is as small as possible. If you picture $L$ as the height of a landscape over the space of all possible knob settings, learning becomes physical: start somewhere, and repeatedly step downhill.

z = \sin x \cos y

: a loss landscape with many hills and valleys. Unlike a single bowl, gradient descent starting near different points can slide into different valleys — different local minima.

Which shape you get depends on the model. A linear model trained with squared error has a loss that is exactly a paraboloid — a single bowl like the one used to illustrate convex optimization — so gradient descent always finds the best fit. A deep neural network, with many layers of nonlinear functions composed together, typically has a loss landscape closer to the bumpy surface above: full of local dips, flat plateaus, and saddle points.

UndergraduateEmpirical risk minimization

Definition: Loss function

A loss function $\ell(\hat y, y)$ measures how bad it is to predict $\hat y$ when the true value is $y$ : it is $0$ (or small) when $\hat y$ is close to $y$ , and grows as the prediction gets worse. For regression, a common choice is squared error $\ell(\hat y, y) = (\hat y - y)^2$ ; for classification, cross-entropy is standard.

Definition: Empirical risk

Given a training set of $n$ examples $\{(x_i, y_i)\}_{i=1}^n$ and a model $f_\theta$ , the empirical risk $\hat R(\theta)$ is the average loss over the training set. Training a model — including linear regression, logistic regression, and deep networks alike — means finding $\theta$ that minimizes $\hat R(\theta)$ .

\hat R(\theta) = \frac{1}{n}\sum_{i=1}^n \ell(f_\theta(x_i), y_i)

In reality, we care about performance on new data, not just the training set: the population (true) risk $R(\theta)$ averages the loss over the entire underlying distribution $\mathcal{D}$ that examples come from, not just the $n$ we happened to sample. Since $\mathcal{D}$ is unknown, $R(\theta)$ cannot be computed directly — $\hat R(\theta)$ is only a proxy, and the gap $R(\theta) - \hat R(\theta)$ is called the generalization gap.

R(\theta) = \mathbb{E}_{(x,y)\sim \mathcal{D}}[\ell(f_\theta(x), y)]

UndergraduateGradient descent and its convergence rate

Gradient descent minimizes $\hat R(\theta)$ (written $L(\theta)$ below for brevity) by repeatedly moving in the direction that decreases $L$ fastest: the negative gradient. At each step, we compute $\nabla L(\theta_k)$ and move a small distance controlled by the learning rate $\eta > 0$ .

\theta_{k+1} = \theta_k - \eta \nabla L(\theta_k)

Common training algorithms and their update rules
Method	Update rule	Per-step cost	Typical use
Batch gradient descent	$\theta_{k+1} = \theta_k - \eta \nabla L(\theta_k)$	$O(n)$ per step	Small datasets, exact gradient
Stochastic gradient descent	$\theta_{k+1} = \theta_k - \eta \nabla \ell_i(\theta_k)$	$O(1)$ per step	Huge datasets, noisy but fast updates
Mini-batch SGD	$\theta_{k+1} = \theta_k - \eta \dfrac{1}{b}\sum_{i \in B} \nabla \ell_i(\theta_k)$	$O(b)$ per step	Standard for deep learning; balances speed and stability
Momentum	$v_{k+1} = \beta v_k + \nabla L(\theta_k), \quad \theta_{k+1} = \theta_k - \eta v_{k+1}$	$O(n)$ or $O(b)$ per step	Accelerates through narrow valleys and plateaus

Theorem: Convergence rate of gradient descent

Let $L$ be convex and $L$ -smooth (its gradient is Lipschitz with constant $L$ : $\|\nabla L(x) - \nabla L(y)\| \le L\|x-y\|$ ), and let $\theta^\star$ minimize $L$ . Gradient descent with step size $\eta = 1/L$ satisfies $L(\theta_K) - L(\theta^\star) \le \dfrac{L\|\theta_0 - \theta^\star\|^2}{2K}$ after $K$ steps.

Why is it true?

Smoothness guarantees each gradient step decreases $L$ by an amount proportional to the squared gradient size, so $L(\theta_k)$ never goes up. Convexity lets us compare that per-step decrease to the still-remaining gap $L(\theta_k) - L(\theta^\star)$ . Adding up the guaranteed decreases over $K$ steps — a telescoping sum — shows the total decrease is bounded, and since the best point in a non-increasing sequence is at least as good as the average, the final gap must shrink like $1/K$ .

Proof

By $L$ -smoothness, for any $x, y$ : $L(y) \le L(x) + \nabla L(x)^\top (y-x) + \dfrac{L}{2}\|y-x\|^2$ . Setting $y = \theta_{k+1} = \theta_k - \dfrac{1}{L}\nabla L(\theta_k)$ and $x = \theta_k$ gives the descent lemma: $L(\theta_{k+1}) \le L(\theta_k) - \dfrac{1}{2L}\|\nabla L(\theta_k)\|^2$ .

By convexity of $L$ : $L(\theta_k) \le L(\theta^\star) + \nabla L(\theta_k)^\top(\theta_k - \theta^\star)$ . Adding this to the descent lemma: $L(\theta_{k+1}) - L(\theta^\star) \le \nabla L(\theta_k)^\top(\theta_k - \theta^\star) - \dfrac{1}{2L}\|\nabla L(\theta_k)\|^2$ .

Complete the square on the right-hand side: $\nabla L(\theta_k)^\top(\theta_k - \theta^\star) - \dfrac{1}{2L}\|\nabla L(\theta_k)\|^2 = \dfrac{L}{2}\left(\|\theta_k - \theta^\star\|^2 - \left\|\theta_k - \theta^\star - \dfrac{1}{L}\nabla L(\theta_k)\right\|^2\right) = \dfrac{L}{2}\left(\|\theta_k - \theta^\star\|^2 - \|\theta_{k+1} - \theta^\star\|^2\right)$ , since $\theta_{k+1} - \theta^\star = \theta_k - \theta^\star - \dfrac{1}{L}\nabla L(\theta_k)$ .

So $L(\theta_{k+1}) - L(\theta^\star) \le \dfrac{L}{2}\left(\|\theta_k - \theta^\star\|^2 - \|\theta_{k+1} - \theta^\star\|^2\right)$ . Summing this for $k = 0, \dots, K-1$ , the right-hand side telescopes to $\dfrac{L}{2}\left(\|\theta_0 - \theta^\star\|^2 - \|\theta_K - \theta^\star\|^2\right) \le \dfrac{L}{2}\|\theta_0 - \theta^\star\|^2$ .

The descent lemma also shows $L(\theta_k)$ is non-increasing, so $L(\theta_K)$ is at most the average of $L(\theta_1), \dots, L(\theta_K)$ : $L(\theta_K) - L(\theta^\star) \le \dfrac{1}{K}\sum_{k=1}^K \left(L(\theta_k) - L(\theta^\star)\right) \le \dfrac{L\|\theta_0 - \theta^\star\|^2}{2K}$ , which is exactly the claimed bound.

UndergraduateBackpropagation: the chain rule at scale

A network with $D$ layers computes $z^{(l)} = W^{(l)} a^{(l-1)} + b^{(l)}$ , then applies a nonlinearity $a^{(l)} = \sigma(z^{(l)})$ , for $l = 1, \dots, D$ , with $a^{(0)}$ the input. Naively applying the chain rule to compute $\partial L / \partial W^{(l)}$ for every layer separately would repeat the same sub-computations over and over, costing time that grows quadratically with depth.

z^{(l)} = W^{(l)} a^{(l-1)} + b^{(l)}, \qquad a^{(l)} = \sigma(z^{(l)})

Backpropagation avoids the repeated work by computing one intermediate quantity per layer, the error term $\delta^{(l)} = \partial L / \partial z^{(l)}$ , starting from the output layer and working backward. Each $\delta^{(l)}$ is built directly from $\delta^{(l+1)}$ , reusing it instead of recomputing it, and the weight gradient falls out immediately from $\delta^{(l)}$ and the layer's input $a^{(l-1)}$ :

\delta^{(l)} = \left((W^{(l+1)})^\top \delta^{(l+1)}\right) \odot \sigma'(z^{(l)}), \qquad \frac{\partial L}{\partial W^{(l)}} = \delta^{(l)} (a^{(l-1)})^\top

Because each $\delta^{(l)}$ is computed once and reused by the layer before it, the total cost of one backward pass is proportional to $D$ , the same order as one forward pass — not $D^2$ . This is what makes training networks with dozens or hundreds of layers computationally feasible.

UndergraduateDimensionality reduction: SVD and PCA

Real datasets often have redundant, correlated features. The singular value decomposition (SVD) factors any matrix $A$ (for example, $n$ centered data points as rows) into $A = U \Sigma V^\top$ , where $U$ and $V$ have orthonormal columns and $\Sigma$ is diagonal with non-negative entries $\sigma_1 \ge \sigma_2 \ge \cdots \ge 0$ , the singular values.

A = U \Sigma V^\top

Principal component analysis (PCA) uses the right singular vectors $v_1, v_2, \dots$ (columns of $V$ ) as new coordinate axes: $v_1$ is the direction along which the (centered) data spreads out the most, $v_2$ the next most, and so on, each orthogonal to the previous ones. Keeping only the top $k$ directions and discarding the rest compresses the data while keeping as much of its variation as any $k$ -dimensional projection possibly can.

Theorem: Eckart–Young–Mirsky theorem

Let $A \in \mathbb{R}^{m \times n}$ have singular values $\sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_r > 0$ , and for $k < r$ let $A_k = \sum_{i=1}^k \sigma_i u_i v_i^\top$ keep only the top $k$ singular components. Then for every matrix $B$ with $\operatorname{rank}(B) \le k$ , $\|A - B\|_2 \ge \sigma_{k+1}$ , and this bound is attained: $\|A - A_k\|_2 = \sigma_{k+1}$ . So $A_k$ is a best rank- $k$ approximation of $A$ in the operator norm.

Why is it true?

The singular values measure how much $A$ stretches vectors along each orthogonal direction $v_i$ ; keeping the largest ones and dropping the smallest throws away the directions $A$ stretches least. Any other rank- $k$ matrix $B$ must fail to reconstruct some direction among the top $k+1$ singular directions (there are too many of them to fit in a $k$ -dimensional image), and that failure costs at least $\sigma_{k+1}$ .

Proof

Achievability. Since $U, V$ have orthonormal columns, $\|A - A_k\|_2 = \|U(\Sigma - \Sigma_k)V^\top\|_2 = \|\Sigma - \Sigma_k\|_2$ , where $\Sigma_k$ keeps the top $k$ singular values and zeroes the rest. $\Sigma - \Sigma_k$ is diagonal with entries $0, \dots, 0, \sigma_{k+1}, \dots, \sigma_r$ , so its operator norm is its largest entry, $\sigma_{k+1}$ .

Optimality. Let $B$ be any matrix with $\operatorname{rank}(B) \le k$ ; its null space (kernel) has dimension at least $n - k$ . Let $S = \operatorname{span}(v_1, \dots, v_{k+1})$ , a $(k+1)$ -dimensional subspace. Since $\dim(\ker B) + \dim(S) \ge (n-k) + (k+1) = n+1 > n$ , the two subspaces must intersect in more than just the origin: there exists a unit vector $z \in S$ with $Bz = 0$ .

Write $z = \sum_{i=1}^{k+1} c_i v_i$ with $\sum_{i=1}^{k+1} c_i^2 = 1$ (since $z$ has unit norm and the $v_i$ are orthonormal). Because $Bz = 0$ , $\|(A-B)z\| = \|Az\| = \left\|\sum_{i=1}^{k+1} c_i \sigma_i u_i\right\| = \sqrt{\sum_{i=1}^{k+1} c_i^2 \sigma_i^2}$ , using orthonormality of the $u_i$ .

Since $\sigma_i \ge \sigma_{k+1}$ for every $i \le k+1$ , $\sum_{i=1}^{k+1} c_i^2 \sigma_i^2 \ge \sigma_{k+1}^2 \sum_{i=1}^{k+1} c_i^2 = \sigma_{k+1}^2$ . Hence $\|A-B\|_2 \ge \|(A-B)z\| \ge \sigma_{k+1}$ for every rank- $\le k$ matrix $B$ , which together with achievability proves the theorem.

The same $A_k$ also minimizes the Frobenius-norm error among rank- $\le k$ matrices, with the discarded error equal to the discarded singular values: $\|A - A_k\|_F = \sqrt{\sum_{i=k+1}^r \sigma_i^2}$ . Since the total variance of the data equals $\sum_i \sigma_i^2$ , this is exactly the variance PCA leaves behind when it keeps only $k$ components — the smaller this quantity, the more faithfully $k$ dimensions summarize the original data.

\|A - A_k\|_F = \sqrt{\sum_{i=k+1}^r \sigma_i^2}

UndergraduateReal-World Applications and Worked Examples

These tools appear together in almost every practical machine learning system: gradient descent (and its variants) trains recommendation engines, image classifiers, and language models; backpropagation is what makes training deep networks for speech recognition and computer vision tractable; and PCA is used to compress genomic data, denoise sensor readings, and visualize high-dimensional financial or scientific datasets in two or three dimensions.

Example: One step of gradient descent for house-price prediction

A real-estate site models predicted price (in $100{,}000$ ) as $\hat y = \theta x$ , where $x$ is house size (in $1{,}000$ sq ft). Three training houses give $(x_i, y_i)$ : $(1, 2)$ , $(2, 3)$ , $(3, 5)$ . Using empirical squared-error risk $\hat R(\theta) = \tfrac{1}{3}\sum_{i=1}^3 (\theta x_i - y_i)^2$ and starting from $\theta_0 = 1$ , find $\theta_1$ after one gradient descent step with $\eta = 0.05$ .

Solution

First differentiate: $\nabla \hat R(\theta) = \tfrac{2}{3}\sum_{i=1}^3 x_i(\theta x_i - y_i) = \tfrac{2}{3}\big[(\theta - 2) + 2(2\theta - 3) + 3(3\theta - 5)\big] = \tfrac{2}{3}(14\theta - 23)$ .

Evaluate at $\theta_0 = 1$ : predictions are $\hat y_i = 1, 2, 3$ ; residuals $\hat y_i - y_i = -1, -1, -2$ ; so $\nabla \hat R(1) = \tfrac{2}{3}(14 \cdot 1 - 23) = \tfrac{2}{3}(-9) = -6$ .

Apply the update rule: $\theta_1 = \theta_0 - \eta \nabla \hat R(\theta_0) = 1 - 0.05 \times (-6) = 1 + 0.3 = 1.3$ .

The gradient was negative, so gradient descent correctly pushed $\theta$ upward: increasing the price-per-size slope reduces the squared error on this training set, moving $\theta$ toward the true least-squares optimum (which turns out to be $\theta^\star = 23/14 \approx 1.64$ , found by setting $\nabla \hat R(\theta) = 0$ ).

Example: Finding the principal component of a covariance matrix

Two centered features have covariance matrix $\Sigma = \begin{pmatrix} 4 & 2 \\ 2 & 1 \end{pmatrix}$ . Find the direction of the first principal component and the fraction of total variance it explains.

Solution

The principal directions are eigenvectors of $\Sigma$ . The characteristic equation is $\det(\Sigma - \lambda I) = (4-\lambda)(1-\lambda) - 2 \cdot 2 = \lambda^2 - 5\lambda = 0$ , giving eigenvalues $\lambda_1 = 5$ , $\lambda_2 = 0$ .

For $\lambda_1 = 5$ : solve $(\Sigma - 5I)v = 0$ , i.e. $\begin{pmatrix} -1 & 2 \\ 2 & -4 \end{pmatrix} v = 0$ , giving $v \propto (2, 1)$ . Normalizing, the first principal component direction is $u_1 = \tfrac{1}{\sqrt5}(2, 1)$ .

Total variance is $\operatorname{tr}(\Sigma) = 4 + 1 = 5$ , which equals $\lambda_1 + \lambda_2 = 5 + 0$ . The fraction of variance explained by the first component is $\lambda_1/(\lambda_1+\lambda_2) = 5/5 = 100\%$ .

Since $\lambda_2 = 0$ , this dataset is exactly one-dimensional: every point lies exactly along the direction $(2,1)$ , so projecting onto that single line (a rank-1 approximation) loses no information at all, matching the Eckart–Young formula $\|A-A_1\|_F = \sqrt{\lambda_2} = 0$ .

Research frontier as of 2026

As of 2026, a central open question is why heavily overparameterized deep networks generalize well despite being able to fit random labels perfectly, as Zhang et al. (2017) demonstrated empirically, undermining classical capacity-based explanations of generalization. Belkin et al. (2019) showed this fits into a broader double descent curve: as model size grows past the point where training error first reaches zero, test error can rise sharply and then fall again, contradicting the textbook U-shaped bias–variance trade-off. On the optimization side, the neural tangent kernel framework of Jacot, Gabriel and Hongler (2018) shows that in the infinite-width limit, gradient descent on a neural network behaves like a fixed kernel method, giving one rigorous (if idealized) lens on why training such non-convex landscapes so often finds good solutions in practice; extending these results to realistic finite-width, finite-time training regimes remains an active research direction.

The loss is $L(\theta) = (\theta - 4)^2$ . Starting from $\theta_0 = 0$ with learning rate $\eta = 0.1$ , what is $\theta_1$ after one gradient descent step?

Which formula correctly defines the empirical risk $\hat R(\theta)$ for a training set of $n$ examples $\{(x_i,y_i)\}$ and loss function $\ell$ ?

A matrix has singular values $\sigma_1 = 6$ , $\sigma_2 = 3$ , $\sigma_3 = 2$ . What is the Frobenius-norm error $\|A - A_1\|_F$ of the best rank-1 approximation $A_1$ ?

A self-driving car's vision network has 50 layers. Computing every weight's gradient with backpropagation costs about the same as one extra forward pass through the network, not 50 times more. Which property of backpropagation explains this?

References

Ian Goodfellow, Yoshua Bengio, Aaron Courville (2016). Deep Learning
Christopher M. Bishop (2006). Pattern Recognition and Machine Learning
Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht, Oriol Vinyals (2017). Understanding deep learning requires rethinking generalization · arXiv:1611.03530
Arthur Jacot, Franck Gabriel, Clément Hongler (2018). Neural Tangent Kernel: Convergence and Generalization in Neural Networks · arXiv:1806.07572
Mikhail Belkin, Daniel Hsu, Siyuan Ma, Soumik Mandal (2019). Reconciling modern machine learning practice and the classical bias-variance trade-off · arXiv:1812.11118