Linear Algebra &
Differential Equations

01Vectors, norms, and inner products

A vector $\mathbf{x}\in\mathbb{R}^n$ is an ordered tuple $\mathbf{x}=(x_1,\dots,x_n)^\top$. Vectors live in a vector space: a set closed under addition and scalar multiplication with the usual axioms. In ML, data points, gradients, and parameters are all vectors; in quant finance, asset returns, portfolio weights, and yield-curve shifts are vectors.

Inner product

The Euclidean inner product of $\mathbf{x},\mathbf{y}\in\mathbb{R}^n$ is $$\langle\mathbf{x},\mathbf{y}\rangle \;=\; \mathbf{x}^\top\mathbf{y} \;=\; \sum_{i=1}^n x_i y_i \;=\; \|\mathbf{x}\|\,\|\mathbf{y}\|\cos\theta.$$ It measures alignment: orthogonality $\iff \langle\mathbf{x},\mathbf{y}\rangle=0$. The Cauchy–Schwarz inequality states $|\langle\mathbf{x},\mathbf{y}\rangle|\le\|\mathbf{x}\|\,\|\mathbf{y}\|$.

Norms

Geometry you must internalize

• Angle: $\cos\theta = \dfrac{\mathbf{x}^\top\mathbf{y}}{\|\mathbf{x}\|\|\mathbf{y}\|}$ — this is cosine similarity.
• Projection of $\mathbf{x}$ onto $\mathbf{y}$: $\mathrm{proj}_\mathbf{y}\mathbf{x}=\dfrac{\mathbf{x}^\top\mathbf{y}}{\mathbf{y}^\top\mathbf{y}}\,\mathbf{y}$.
• Parallelogram law: $\|\mathbf{x}+\mathbf{y}\|^2+\|\mathbf{x}-\mathbf{y}\|^2=2\|\mathbf{x}\|^2+2\|\mathbf{y}\|^2$.

02Matrices and their algebra

A matrix $A\in\mathbb{R}^{m\times n}$ represents a linear map $\mathbb{R}^n\to\mathbb{R}^m$ via $\mathbf{x}\mapsto A\mathbf{x}$. Column $j$ of $A$ is the image of $\mathbf{e}_j$, the $j$th standard basis vector. Matrix multiplication $AB$ composes maps; it is associative, distributive, but not commutative.

Matrix norms

03Rank, span, basis, and the four fundamental subspaces

The span of vectors $\{\mathbf{v}_1,\dots,\mathbf{v}_k\}$ is the set of all their linear combinations. A set is linearly independent if $\sum c_i\mathbf{v}_i=\mathbf{0}\implies c_i=0$ for all $i$. A basis of a subspace is a linearly independent spanning set; its size is the dimension.

For $A\in\mathbb{R}^{m\times n}$, Strang's four subspaces are:

Rank is preserved under invertible maps: $\mathrm{rank}(PAQ)=\mathrm{rank}(A)$ for invertible $P,Q$. Also $\mathrm{rank}(AB)\le\min(\mathrm{rank}(A),\mathrm{rank}(B))$ and $\mathrm{rank}(A^\top A)=\mathrm{rank}(A)$.

04Determinants and trace

05Eigenvalues and eigenvectors

If $A$ has $n$ linearly independent eigenvectors, it is diagonalizable: $$A=V\Lambda V^{-1},\qquad \Lambda=\mathrm{diag}(\lambda_1,\dots,\lambda_n),\quad V=[\mathbf{v}_1\,\cdots\,\mathbf{v}_n].$$ Then $A^k=V\Lambda^k V^{-1}$ and $f(A)=V f(\Lambda) V^{-1}$ for analytic $f$. This is why eigendecomposition accelerates matrix powers (Markov chains, PageRank, recurrent neural network — RNN — dynamics) and matrix functions (matrix exponentials).

Facts worth memorizing

• $\sum\lambda_i=\mathrm{tr}(A)$ and $\prod\lambda_i=\det A$.
• Similar matrices $A$ and $P^{-1}AP$ share eigenvalues.
• Real symmetric matrices have real eigenvalues and an orthonormal eigenbasis (Spectral Theorem).
• If $A\mathbf{v}=\lambda\mathbf{v}$ then $(A-\mu I)\mathbf{v}=(\lambda-\mu)\mathbf{v}$ — the shift trick used in deflation and iterative solvers.
• Gershgorin discs: each eigenvalue lies in $\bigcup_i \{z : |z-A_{ii}|\le \sum_{j\neq i}|A_{ij}|\}$.
• Power iteration: $\mathbf{v}_{k+1}=A\mathbf{v}_k/\|A\mathbf{v}_k\|$ converges to the dominant eigenvector when $|\lambda_1|>|\lambda_2|$.

06Symmetric, positive-definite, orthogonal

Tests for positive definiteness

1. All eigenvalues $> 0$.
2. All leading principal minors $> 0$ (Sylvester's criterion).
3. There exists non-singular $R$ with $A=R^\top R$ (e.g., Cholesky).
4. $\mathbf{x}^\top A \mathbf{x}>0$ for all nonzero $\mathbf{x}$.

07Singular value decomposition and the decomposition toolkit

The SVD is the single most important factorization in machine learning (ML). It encodes geometry (any linear map is rotate → scale axes → rotate), solves rank-deficient least squares, and powers principal component analysis (PCA), recommender systems, and model compression.

Consequences

• $\|A\|_2=\sigma_1$, $\|A\|_F=\sqrt{\sum\sigma_i^2}$, $\|A\|_*=\sum\sigma_i$.
• Condition number $\kappa_2(A)=\sigma_1/\sigma_r$.
• Eckart–Young–Mirsky: the best rank-$k$ approximation in Frobenius (or spectral) norm is $A_k=\sum_{i=1}^k \sigma_i\mathbf{u}_i\mathbf{v}_i^\top$ with error $\|A-A_k\|_F^2=\sum_{i>k}\sigma_i^2$.
• Relation to eigendecomposition: $A^\top A = V\Sigma^\top\Sigma V^\top$ (columns of $V$ are eigenvectors of $A^\top A$, $\sigma_i^2$ are eigenvalues).

The decomposition menu

08Projections, least squares, pseudoinverse

Given a subspace $\mathcal{S}=\mathcal{C}(A)$ with $A$ of full column rank, the orthogonal projection onto $\mathcal{S}$ is $$P=A(A^\top A)^{-1}A^\top,$$ which satisfies $P^2=P=P^\top$. Geometrically, $P\mathbf{b}$ is the closest point in $\mathcal{S}$ to $\mathbf{b}$.

Least squares

For an overdetermined system $A\mathbf{x}\approx\mathbf{b}$ the normal equations are $$A^\top A\mathbf{x}=A^\top \mathbf{b}\quad\Longleftrightarrow\quad \mathbf{x}^* = (A^\top A)^{-1}A^\top\mathbf{b}.$$ Numerically, solve via QR: if $A=QR$, then $\mathbf{x}^*=R^{-1}Q^\top\mathbf{b}$. Squaring $A^\top A$ squares the condition number, so avoid forming it in ill-conditioned problems.

Moore–Penrose pseudoinverse

For any $A$ with SVD $A=U\Sigma V^\top$, the pseudoinverse is $$A^+ = V\Sigma^+ U^\top,\qquad \Sigma^+_{ii}=\begin{cases}1/\sigma_i & \sigma_i>0\\ 0 & \text{else}\end{cases}.$$ The minimum-norm least-squares solution to $A\mathbf{x}=\mathbf{b}$ is $\mathbf{x}^*=A^+\mathbf{b}$. When $A$ has full column rank, $A^+=(A^\top A)^{-1}A^\top$.

Ridge regression

Adding $\ell_2$ regularization yields $$\mathbf{x}_\lambda = (A^\top A + \lambda I)^{-1}A^\top\mathbf{b},$$ which always has a unique solution (the matrix is PD for $\lambda>0$) and shrinks components along small singular directions: in SVD coordinates, each coefficient is scaled by $\sigma_i/(\sigma_i^2+\lambda)$.

09Matrix calculus for machine learning

Gradients drive training. Adopt numerator layout throughout: if $f:\mathbb{R}^n\to\mathbb{R}^m$, the Jacobian $\partial f/\partial \mathbf{x}\in\mathbb{R}^{m\times n}$; the gradient of a scalar is a row vector (or its transpose, at your preference — just be consistent).

The cookbook

In the table below, $\mathbf{x}\in\mathbb{R}^n$ is the variable, $\mathbf{a},\mathbf{b}\in\mathbb{R}^n$ are constant vectors, and $A,B$ are constant matrices of compatible size (so each expression is well-posed).

Chain rule and backprop

For $\mathbf{h}=g(\mathbf{x})$ and $\ell=f(\mathbf{h})$, the gradient $\nabla_\mathbf{x}\ell = J_g(\mathbf{x})^\top\,\nabla_\mathbf{h}\ell$. Backpropagation is just repeated application of this rule — with computation graph structure used so each Jacobian-vector product runs in cost proportional to the forward pass.

Hessian

$H=\nabla^2 f\in\mathbb{R}^{n\times n}$ is the Jacobian of the gradient. It is symmetric (for $C^2$ functions). Convexity $\iff H\succeq 0$ everywhere. Newton's method: $\mathbf{x}_{k+1}=\mathbf{x}_k - H^{-1}\nabla f$.

10ML and quant applications

The covariance matrix

The sample covariance from a data matrix $X\in\mathbb{R}^{n\times d}$ with centered rows is $\hat\Sigma = \tfrac{1}{n-1}X^\top X\in\mathbb{R}^{d\times d}$ — PSD always, and positive-definite generically when $n>d$. The same object powers PCA (§10), the Mahalanobis distance $\|\mathbf{x}-\boldsymbol\mu\|_\Sigma := \sqrt{(\mathbf{x}-\boldsymbol\mu)^\top\Sigma^{-1}(\mathbf{x}-\boldsymbol\mu)}$, whitening transforms $\mathbf{z} = \Sigma^{-1/2}(\mathbf{x}-\boldsymbol\mu)$, and the Gaussian log-likelihood. In finance, $\Sigma$ is the asset-return covariance; portfolio variance for weights $\mathbf{w}$ is $\mathbf{w}^\top\Sigma\mathbf{w}$.

Kernel functions

Three useful intuitions: (i) $k(\mathbf{x},\mathbf{x}')$ quantifies how alike $\mathbf{x}$ and $\mathbf{x}'$ are, with $k(\mathbf{x},\mathbf{x})=\|\phi(\mathbf{x})\|_\mathcal{H}^2$ the "self-similarity"; (ii) the PD condition says no weighted average of pairwise similarities can be negative — similarity has no "imaginary" mass; (iii) kernels let us lift non-linear problems into linear ones: running a linear method in feature space is the same as running a non-linear method in the original space, through $k$ alone.

New kernels are built from old by a small algebra — sums, positive-scalar multiples, products, and pointwise limits of PD kernels are PD; $k(\mathbf{x},\mathbf{x}') = k_1(\mathbf{x},\mathbf{x}')\cdot k_2(\mathbf{x},\mathbf{x}')$ corresponds to the tensor-product feature map. This closure is what makes compositional kernel design tractable.

Kernel Gram matrices

Common kernels:

• Linear: $k(\mathbf{x},\mathbf{x}') = \mathbf{x}^\top\mathbf{x}'$ — feature map is identity.
• Polynomial: $k(\mathbf{x},\mathbf{x}') = (\mathbf{x}^\top\mathbf{x}' + c)^p$ — finite-dim feature space of monomials up to degree $p$.
• Gaussian / RBF: $k(\mathbf{x},\mathbf{x}') = \exp(-\|\mathbf{x}-\mathbf{x}'\|^2/(2\ell^2))$ — infinite-dimensional feature space, length-scale $\ell$.
• Matérn family — RBF's roughness-controlled cousin, standard in Gaussian-process models.

The Gram matrix is the computational object that kernel methods (support vector machines, kernel ridge regression, Gaussian processes) manipulate instead of the possibly infinite-dimensional feature vectors — the "kernel trick". Training and prediction can often be written using only $K$ and its inverse, never touching $\phi(\mathbf{x})$.

PCA in one page

1. Center data: $\tilde X = X - \mathbf{1}\bar{\mathbf{x}}^\top$, where $\mathbf{1}\in\mathbb{R}^n$ denotes the all-ones column vector and $\bar{\mathbf{x}}\in\mathbb{R}^d$ is the column mean.
2. Compute covariance $\Sigma = \tfrac{1}{n-1}\tilde X^\top\tilde X$ (or use SVD of $\tilde X$ directly).
3. Eigendecompose $\Sigma=Q\Lambda Q^\top$. Columns of $Q$ are principal directions; $\lambda_i$ are variances.
4. Project: $Z = \tilde X Q_k$ where $Q_k$ holds the top-$k$ eigenvectors.

Equivalently, top-$k$ components come from truncated SVD of $\tilde X$. Variance explained by the first $k$ components is $\sum_{i=1}^k\lambda_i/\sum_i\lambda_i$.

Mean–variance portfolio (Markowitz)

Given expected returns $\boldsymbol\mu\in\mathbb{R}^d$, covariance $\Sigma\succ 0$, and target return $r$, minimize $\mathbf{w}^\top\Sigma\mathbf{w}$ subject to $\mathbf{1}^\top\mathbf{w}=1$ (weights sum to one) and $\boldsymbol\mu^\top\mathbf{w}=r$. Lagrangian yields a closed form where the optimal $\mathbf{w}^*$ is an affine function of $r$, with $\Sigma^{-1}\boldsymbol\mu$ and $\Sigma^{-1}\mathbf{1}$ as building blocks. With a risk-free asset of rate $r_f$, the tangency portfolio is $\propto \Sigma^{-1}(\boldsymbol\mu-r_f\mathbf{1})$.

Sampling from a Gaussian

To sample $\mathbf{x}\sim\mathcal{N}(\boldsymbol\mu,\Sigma)$ with $\Sigma=LL^\top$ (Cholesky): draw $\mathbf{z}\sim\mathcal{N}(\mathbf{0},I)$, set $\mathbf{x}=\boldsymbol\mu + L\mathbf{z}$. This underpins Monte Carlo simulation of correlated assets.

Condition number

$\kappa_2(A)=\sigma_1/\sigma_n$ controls how errors amplify when solving $A\mathbf{x}=\mathbf{b}$: $\|\Delta\mathbf{x}\|/\|\mathbf{x}\|\lesssim \kappa_2(A)\cdot\|\Delta\mathbf{b}\|/\|\mathbf{b}\|$. In deep nets, poorly conditioned Hessians explain why vanilla stochastic gradient descent (SGD) struggles without preconditioning / adaptive methods.

11Exercises · Linear algebra

Each exercise is tagged by difficulty: ★ warm-up, ★★ core, ★★★ stretch. Solutions are hidden behind a toggle — try before peeking.

12First-order ordinary differential equations

An ordinary differential equation (ODE) relates a function $y(t)$ to its derivatives. A first-order ODE has the form $y'=f(t,y)$. Under mild conditions (Picard–Lindelöf: $f$ Lipschitz in $y$), a unique local solution exists given $y(t_0)=y_0$.

Three tractable forms

Canonical examples

Exponential decay / growth: $y'=\lambda y\Rightarrow y(t)=y_0 e^{\lambda t}$.
Logistic: $y'=ry(1-y/K)\Rightarrow y(t)=K/(1+(K/y_0-1)e^{-rt})$. Models population and sigmoidal activations in continuous time.
Newton's cooling: $y'=-k(y-y_\infty)\Rightarrow y=y_\infty+(y_0-y_\infty)e^{-kt}$. The same first-order linear template appears everywhere: exponential moving averages, Ornstein–Uhlenbeck drift, resistor–capacitor (RC) circuits.

13Second-order linear ODEs

Consider $y''+py'+qy=g(t)$ with constants $p,q$. The general solution is $y=y_h+y_p$ where $y_h$ solves the homogeneous equation and $y_p$ is any particular solution.

Homogeneous case, constant coefficients

Try $y=e^{rt}$. Characteristic equation: $r^2+pr+q=0$, with roots $r_{1,2}=\tfrac{-p\pm\sqrt{p^2-4q}}{2}$.

Particular solutions

• Undetermined coefficients: if $g(t)$ is a polynomial, exponential, sine/cosine, or products/sums thereof, guess the same form with unknown coefficients and match.
• Variation of parameters: for $y_h=C_1 y_1+C_2 y_2$, take $y_p = u_1(t)y_1 + u_2(t)y_2$ with $u_1,u_2$ solving a 2×2 system involving the Wronskian $W=y_1 y_2' - y_2 y_1'$.

Damped harmonic oscillator

$my'' + cy' + ky = 0$. With $\omega_0^2 = k/m$, $\zeta = c/(2\sqrt{mk})$:

• $\zeta<1$ — underdamped: oscillation with decaying amplitude, frequency $\omega_0\sqrt{1-\zeta^2}$.
• $\zeta=1$ — critically damped: fastest non-oscillatory return to equilibrium.
• $\zeta>1$ — overdamped: slow exponential decay, no oscillation.

This framework classifies momentum-based optimizers (heavy-ball, Nesterov) viewed as discretizations of damped oscillators.

14Systems of linear ODEs and the matrix exponential

A linear system $\mathbf{x}'(t) = A\mathbf{x}(t)$ with $A\in\mathbb{R}^{n\times n}$ has the unique solution $\mathbf{x}(t) = e^{At}\mathbf{x}(0)$, where the matrix exponential is $$e^{At} = \sum_{k=0}^\infty \frac{(At)^k}{k!}.$$

When $A$ is diagonalizable, $A=V\Lambda V^{-1}$ and $e^{At}=V e^{\Lambda t} V^{-1}$ with $e^{\Lambda t}=\mathrm{diag}(e^{\lambda_i t})$. The solution decomposes into modes: writing $\mathbf{x}(0)=\sum c_i \mathbf{v}_i$, $$\mathbf{x}(t)=\sum_i c_i e^{\lambda_i t}\mathbf{v}_i.$$ Each eigenvalue is a mode: $\mathrm{Re}(\lambda_i)>0$ means growth, $<0$ decay, $\mathrm{Im}(\lambda_i)\neq 0$ oscillation.

Non-homogeneous systems

For $\mathbf{x}'=A\mathbf{x}+\mathbf{g}(t)$, the variation-of-parameters formula gives $$\mathbf{x}(t)=e^{At}\mathbf{x}(0)+\int_0^t e^{A(t-s)}\mathbf{g}(s)\,ds.$$

Jordan form for defective matrices

If $A$ is not diagonalizable, use Jordan form $A=PJP^{-1}$ where $J$ is block-diagonal with Jordan blocks $J_k = \lambda I + N$, $N$ nilpotent. Then $e^{J_k t}=e^{\lambda t}\sum_{j=0}^{k-1}\tfrac{(Nt)^j}{j!}$ — finite sum.

15Stability and phase portraits

An equilibrium $\mathbf{x}^*$ of $\mathbf{x}'=\mathbf{f}(\mathbf{x})$ satisfies $\mathbf{f}(\mathbf{x}^*)=\mathbf{0}$. Linearize: $\mathbf{y}=\mathbf{x}-\mathbf{x}^*$ evolves approximately as $\mathbf{y}'=J\mathbf{y}$ where $J=\partial\mathbf{f}/\partial\mathbf{x}\big|_{\mathbf{x}^*}$.

Lyapunov / spectral criterion

• Asymptotically stable if all eigenvalues of $J$ have $\mathrm{Re}(\lambda_i)<0$.
• Unstable if any $\mathrm{Re}(\lambda_i)>0$.
• Purely imaginary eigenvalues: linearization inconclusive (need higher-order analysis or Lyapunov functions).

2D classification

For a 2×2 Jacobian with trace $\tau$ and determinant $\Delta$, the equilibrium is:

16Partial differential equations — the essentials

Three archetypal linear second-order partial differential equations (PDEs) span most applications in physics, engineering, and quant.

Heat / diffusion equation

$u_t = \alpha u_{xx}$ where $\alpha>0$ is the diffusivity (thermal conductivity / thermal capacity). Fundamental solution on $\mathbb{R}$: $K(x,t)=\tfrac{1}{\sqrt{4\pi\alpha t}}e^{-x^2/(4\alpha t)}$ (the heat kernel). For general initial condition $u_0$, $u(x,t)=(K(\cdot,t)*u_0)(x)$ where $*$ denotes convolution. The Black–Scholes PDE reduces to the heat equation under a change of variables.

Wave equation

$u_{tt} = c^2 u_{xx}$ where $c>0$ is the propagation speed. d'Alembert: $u(x,t)=\tfrac{1}{2}[u_0(x-ct)+u_0(x+ct)] + \tfrac{1}{2c}\int_{x-ct}^{x+ct}v_0(s)\,ds$, given initial position $u_0$ and velocity $v_0$. Finite speed of propagation $c$.

Laplace / Poisson equation

$\Delta u = 0$ (Laplace) or $\Delta u = f$ (Poisson). Models steady-state diffusion, electrostatics, harmonic problems. Solutions inherit strong regularity (harmonic functions satisfy the mean-value property).

Separation of variables (template)

1. Assume $u(x,t)=X(x)T(t)$.
2. Substitute; divide to obtain $T'/T = \alpha X''/X = -\lambda$ (separation constant).
3. Solve the spatial eigenvalue problem with boundary conditions — this yields countably many $\lambda_n$ and modes $X_n$.
4. Superpose: $u(x,t)=\sum_n c_n T_n(t) X_n(x)$, with $c_n$ fixed by initial conditions (Fourier coefficients).

17Stochastic differential equations (quant)

A Brownian motion (or Wiener process) $W_t$ is a continuous stochastic process with $W_0=0$, independent increments, and $W_t - W_s \sim \mathcal{N}(0,t-s)$. Its paths are continuous but nowhere differentiable; informally, $dW_t \sim \mathcal{N}(0,dt)$.

An Itô SDE has the form $$dX_t = \mu(X_t,t)\,dt + \sigma(X_t,t)\,dW_t,$$ with drift $\mu$ and diffusion $\sigma$. Existence and uniqueness hold under Lipschitz and linear-growth conditions on $\mu,\sigma$.

Itô's lemma (the chain rule)

For a function $f(x,t)$ that is twice continuously differentiable in $x$ and once in $t$ (written $f\in C^{2,1}$), and $X_t$ as above, $$df(X_t,t) = \left(f_t + \mu f_x + \tfrac{1}{2}\sigma^2 f_{xx}\right) dt + \sigma f_x\, dW_t.$$ The extra $\tfrac{1}{2}\sigma^2 f_{xx}$ term is the hallmark of Itô calculus; it arises because $(dW_t)^2 = dt$ in the mean-square sense.

Geometric Brownian motion

Geometric Brownian motion (GBM): $dS_t = \mu S_t\,dt + \sigma S_t\,dW_t$. Apply Itô to $\log S_t$: $$d(\log S_t) = (\mu - \tfrac{1}{2}\sigma^2)\,dt + \sigma\, dW_t.$$ Hence $S_t = S_0 \exp\big((\mu-\tfrac{1}{2}\sigma^2)t + \sigma W_t\big)$, lognormally distributed. This is the Black–Scholes stock model.

Black–Scholes PDE

Under the risk-neutral measure (with constant risk-free rate $r$), the price $V(S,t)$ of a European option written on an asset $S_t$ following $dS=rS\,dt+\sigma S\,dW$ satisfies $$V_t + rS V_S + \tfrac{1}{2}\sigma^2 S^2 V_{SS} - rV = 0,$$ with terminal condition $V(S,T)=\mathrm{payoff}(S)$ at maturity $T$. For a European call, the closed-form solution is the celebrated Black–Scholes formula.

Ornstein–Uhlenbeck

Ornstein–Uhlenbeck (OU) process: $dX_t = \theta(\bar X - X_t)\,dt + \sigma\,dW_t$, where $\theta>0$ is the mean-reversion speed, $\bar X$ is the long-run mean, and $\sigma>0$ is the diffusion amplitude. Mean-reverting; stationary distribution $\mathcal{N}(\bar X,\sigma^2/(2\theta))$. Used for interest-rate models (Vasicek) and noise schedules in diffusion generative models.

18Numerical methods for ODEs and SDEs

ODE integrators

We approximate the solution of $\mathbf{x}'=\mathbf{f}(\mathbf{x},t)$ on a grid $t_0

Stability and stiffness

A method is $A$-stable if for the test equation $y'=\lambda y$ with $\mathrm{Re}(\lambda)<0$, $y_n\to 0$ for any step $h>0$. Explicit methods have bounded stability regions; stiff problems (wide spread of time scales) demand implicit methods or exponential integrators.

SDE schemes

• Euler–Maruyama: $X_{n+1}=X_n+\mu\,h + \sigma\,\Delta W_n$ with the Brownian increment $\Delta W_n = W_{t_{n+1}} - W_{t_n} \sim\mathcal{N}(0,h)$. Strong order $1/2$, weak order $1$.
• Milstein: adds an Itô correction term $\tfrac{1}{2}\sigma\sigma_x(\Delta W_n^2-h)$, where $\sigma_x = \partial\sigma/\partial x$, for strong order $1$.
• For GBM, simulate $\log S$ directly via the exact solution to avoid discretization bias.

19ML and quant applications of differential equations

Gradient flow

The continuous-time limit of gradient descent is $\dot{\mathbf{x}} = -\nabla f(\mathbf{x})$. For a quadratic $f(\mathbf{x})=\tfrac{1}{2}\mathbf{x}^\top A\mathbf{x}$ with $A\succ 0$, $\mathbf{x}(t) = e^{-At}\mathbf{x}(0)$. Decay along eigendirection $\mathbf{v}_i$ has rate $\lambda_i$: small eigenvalues give slow modes, which is why preconditioning (dividing by $A$) flattens the convergence landscape.

Momentum as a damped oscillator

Heavy-ball descent, in the continuous limit, becomes $\ddot{\mathbf{x}} + \gamma \dot{\mathbf{x}} + \nabla f(\mathbf{x}) = 0$, where $\gamma>0$ is the damping coefficient. Optimal tuning critically damps: $\gamma = 2\sqrt{\mu}$ for a $\mu$-strongly convex $f$ (i.e., $\nabla^2 f \succeq \mu I$ with strong-convexity constant $\mu>0$). This gives the accelerated $O(\sqrt{\kappa}\log(1/\varepsilon))$ rate of Nesterov's method, where $\kappa = L/\mu$ is the condition number ($L$ being the Lipschitz constant of $\nabla f$) and $\varepsilon$ is the target accuracy.

Neural ODEs

A residual block $\mathbf{h}_{n+1}=\mathbf{h}_n + f_\theta(\mathbf{h}_n)$ is one Euler step of $\dot{\mathbf{h}} = f_\theta(\mathbf{h},t)$. Neural ODEs replace finite depth with the flow of an ODE solved by an adaptive integrator, with gradients via the adjoint method (itself an ODE).

Diffusion generative models

The forward noising process is an SDE $d\mathbf{x}_t = -\tfrac{1}{2}\beta_t\mathbf{x}_t\,dt + \sqrt{\beta_t}\,d\mathbf{W}_t$ (an OU with time-varying coefficient), where $\beta_t>0$ is a prescribed noise schedule controlling how fast data is corrupted toward pure noise, and $\mathbf{W}_t$ is a standard multidimensional Brownian motion. Let $p_t(\mathbf{x})$ denote the marginal density of $\mathbf{x}_t$. The reverse-time SDE (Anderson, 1982) exists and has the form $d\mathbf{x}_t = [-\tfrac{1}{2}\beta_t\mathbf{x}_t - \beta_t\nabla_{\mathbf{x}}\log p_t(\mathbf{x})]\,dt + \sqrt{\beta_t}\,d\bar{\mathbf{W}}_t$, where $\bar{\mathbf{W}}_t$ is a Brownian motion running backward in time. The quantity $\nabla_{\mathbf{x}}\log p_t(\mathbf{x})$ is called the score. Training a score network to approximate it then lets us integrate the reverse SDE (or the deterministic probability-flow ODE with the same marginals) to generate samples.

Quant: term structure and risk

• Vasicek model: $dr_t = a(b-r_t)\,dt + \sigma\,dW_t$ (OU for short rate).
• CIR (Cox–Ingersoll–Ross): $dr_t = a(b-r_t)\,dt + \sigma\sqrt{r_t}\,dW_t$ (non-negative short rate).
• Heston: stochastic volatility model with a CIR variance process; yields a 2D PDE for option prices.
• Feynman–Kac links PDE pricing and expectations under risk-neutral measure: $V(x,t) = \mathbb{E}[\text{discounted payoff} \mid X_t=x]$ solves the corresponding backward PDE.

20Exercises · Differential equations