Stochastic Differential Equations

Dive into Deep Learning · §28.2

Brownian motion, Itô’s calculus, and the SDE
the forward process of every diffusion model.

Why add noise?

Motivation

No fixed map can forget every dataset: a data-independent bijection sends distinct laws to distinct endpoints, so no single deterministic flow carries every p_{\text{data}} to the same \mathcal N(\mathbf 0, \sigma^2 I). Noise does it before any learning:

d\mathbf X = \underbrace{\mathbf f(\mathbf X,t)\,dt}_{\text{drift}} + \underbrace{g(t)\,d\mathbf W}_{\text{diffusion}}.

(A learned flow sends one given p_{\text{data}} there: that is exactly flow matching, the score-matching-diffusion-and-flow-matching section.)

image/svg+xml Matplotlib v3.10.8, https://matplotlib.org/

01

Brownian motion

axioms, covariance, simulation

The Wiener process

Brownian motion

The \pm\sqrt{\Delta t} random walk has a non-degenerate limit W_t: W_0=0, independent increments, W_t-W_s\sim\mathcal N(0,t-s), continuous paths. Simulate with \Delta W = \sqrt{\Delta t}\,\xi, \xi\sim\mathcal N(0,1).

Covariance and non-differentiability

Brownian motion

\operatorname{Cov}(W_s,W_t) = \min(s,t).

Proof. Split W_t = W_s + (W_t-W_s); the cross term factorizes by independence and both factors are mean-zero. \blacksquare

But \Delta W/\Delta t = \xi/\sqrt{\Delta t}\to\pm\infty: paths are nowhere differentiable. We must integrate, never differentiate.

rng = np.random.default_rng(42)
T, n, n_paths = 1.0, 500, 20000
dt = T / n
t = np.linspace(0, T, n + 1)
dW = np.sqrt(dt) * rng.standard_normal((n_paths, n))
W = np.concatenate([np.zeros((n_paths, 1)), np.cumsum(dW, axis=1)], axis=1)
for s in [0.25, 0.5, 1.0]:
    k = round(s / dt)
    print(f'Var(W_t) at t={s:.2f}: {W[:, k].var():.4f}   (theory: {s:.2f})')
d2l.plot(t, [W[i] for i in range(8)] + [2 * np.sqrt(t), -2 * np.sqrt(t)],
         't', '$W_t$', fmts=['-'] * 8 + ['k--', 'k--'], figsize=(5, 3))

Var(W_t) at t=0.25: 0.2486   (theory: 0.25)
Var(W_t) at t=0.50: 0.4950   (theory: 0.50)
Var(W_t) at t=1.00: 1.0014   (theory: 1.00)

02

Itô calculus

quadratic variation, the integral, the chain rule

Quadratic variation

Why ordinary calculus fails

For a smooth path \sum(\Delta x_i)^2\to 0. For Brownian motion it converges to the elapsed time:

\sum_i (\Delta W_i)^2 \longrightarrow t.

The Itô multiplication table

Why ordinary calculus fails

Proof. \mathbb E[\sum(\Delta W_i)^2]=\sum\Delta t_i=t exactly, and the variance is \sum 2\Delta t_i^2 \le 2\delta t\to 0. \blacksquare

(dW)^2 = dt, \qquad dW\,dt = 0, \qquad (dt)^2 = 0. A squared increment is no longer negligible: it is dt.

The Itô integral locks the left endpoint

The integral

Evaluate the integrand before the increment, \int_0^t G_s\,dW_s = \lim\sum G_{t_i}(W_{t_{i+1}}-W_{t_i}). Then G_{t_i}\perp\Delta W_i, so

\mathbb E\Bigl[\textstyle\int_0^t G\,dW\Bigr] = 0, \qquad \mathbb E\Bigl[\bigl(\textstyle\int_0^t G\,dW\bigr)^2\Bigr] = \int_0^t \mathbb E[G_s^2]\,ds.

The Itô isometry turns a stochastic variance into an ordinary integral.

Itô’s lemma: a chain rule with a correction

The chain rule

Expand \phi(X) to second order; since (dX)^2=g^2\,dt survives,

d\phi = \bigl(\phi_t + f\,\phi_x + \tfrac12 g^2\,\phi_{xx}\bigr)dt + g\,\phi_x\,dW.

The \tfrac12 g^2\phi_{xx} term is the seed of the Fokker–Planck diffusion. Check \phi=x^2, X=W: d(W^2)=2W\,dW+dt, so \int_0^t W\,dW=\tfrac12(W_t^2-t):

rng = np.random.default_rng(3)
T, n = 1.0, 1000
dt = T / n
t = np.linspace(0, T, n + 1)
dW = np.sqrt(dt) * rng.standard_normal(n)
W = np.concatenate([[0.0], np.cumsum(dW)])
ito = np.concatenate([[0.0], np.cumsum(2 * W[:-1] * dW)])    # left endpoint
print(f'W_T^2 - int 2W dW = {W[-1]**2 - ito[-1]:.4f}   (Ito correction: T = {T})')
print(f'max |gap_t - t| along the path = {np.abs(W**2 - ito - t).max():.4f}')
d2l.plot(t, [W**2, ito, ito + t], 't', 'value',
         legend=['$W_t^2$', 'naive $\\int_0^t 2W\\,dW$', '$\\int_0^t 2W\\,dW + t$'],
         figsize=(5, 3))

W_T^2 - int 2W dW = 1.0150   (Ito correction: T = 1.0)
max |gap_t - t| along the path = 0.0235

03

SDEs and Euler–Maruyama

the precise meaning, simulation, convergence

An SDE is a distribution over paths

The objects

dX = f\,dt + g\,dW means X_t = X_0 + \int_0^t f\,ds + \int_0^t g\,dW_s. Here f is the mean velocity (\mathbb E[dX]=f\,dt) and g^2 the variance rate (\operatorname{Var}(dX)=g^2\,dt); a solution is an ensemble of jittery paths, and g=0 recovers the ODE.

Euler–Maruyama: Euler plus a √Δt kick

Simulation

X_{n+1} = X_n + f(X_n,t_n)\,\Delta t + g(X_n,t_n)\sqrt{\Delta t}\,\xi_n.

The noise scales as \sqrt{\Delta t}, not \Delta t; g=0 is exactly forward Euler:

def euler_maruyama(f, g, x0, t, dW):
    """Simulate dX = f dt + g dW on the grid t for a batch of paths."""
    X = np.empty((len(x0), len(t)))
    X[:, 0] = x0
    for k in range(len(t) - 1):
        h = t[k + 1] - t[k]
        X[:, k + 1] = X[:, k] + f(X[:, k], t[k]) * h + g(X[:, k], t[k]) * dW[:, k]
    return X

t = np.linspace(0, 2.0, 201)
X = euler_maruyama(lambda x, s: -x, lambda x, s: 0.0,
                   np.ones(1), t, np.zeros((1, 200)))
print('g = 0 recovers Euler on dx/dt = -x: '
      f'max |X_t - e^(-t)| = {np.abs(X[0] - np.exp(-t)).max():.5f}')
g = 0 recovers Euler on dx/dt = -x: max |X_t - e^(-t)| = 0.00185

Paths stray; laws agree

Convergence

Strong error tracks one path against a fine reference driven by the same increments; weak error compares only the terminal laws. The coarse path visibly strays while the histograms already coincide, and marginals are all a diffusion model needs:

The measured orders

Convergence

Strong order is \tfrac12 in general but 1 for additive noise: the omitted Milstein term carries \partial_x g, which vanishes. Weak order is 1 regardless:

strong-order slopes:  OU (additive) 1.02,  GBM (multiplicative) 0.45

Measured slopes: additive OU \approx 1, multiplicative GBM \approx 0.5.

04

The Ornstein–Uhlenbeck process

mean reversion, exact solution, the VP marginal

Mean reversion vs noise

OU

dX = -\theta X\,dt + \sigma\,dW, \qquad \theta > 0.

A restoring drift pulls toward 0 while constant diffusion pushes out:

Solved by an integrating factor

OU

Apply Itô to e^{\theta t}X: the correction vanishes (\phi_{xx}=0) and the drift cancels:

X_t = X_0 e^{-\theta t} + \sigma\!\int_0^t e^{-\theta(t-s)}\,dW_s.

The isometry gives the variance, so the transition kernel is

X_t\mid X_0 \sim \mathcal N\!\Bigl(X_0 e^{-\theta t},\; \tfrac{\sigma^2}{2\theta}\bigl(1-e^{-2\theta t}\bigr)\Bigr).

The stationary distribution

OU

As t\to\infty the mean dies and the variance saturates, independent of the start:

X_\infty \sim \mathcal N\!\bigl(0,\ \sigma^2/(2\theta)\bigr),

with relaxation time 1/\theta. The path cloud and its endpoint histogram confirm it:

t=0.5:  mean  1.2237 vs  1.2131,   std 0.4997 vs 0.5060
t=1.0:  mean  0.7336 vs  0.7358,   std 0.5912 vs 0.5918
t=2.0:  mean  0.2712 vs  0.2707,   std 0.6304 vs 0.6305
t=4.0:  mean  0.0406 vs  0.0366,   std 0.6399 vs 0.6363
fraction of paths inside the +-2 sigma band at t=4: 0.952   (Gaussian: 0.954)

Variance-preserving = DDPM’s forward marginal

VP

Choose \sigma^2 = 2\theta and write \bar\alpha_t = e^{-2\theta t}:

X_t = \sqrt{\bar\alpha_t}\,X_0 + \sqrt{1-\bar\alpha_t}\,\epsilon, \qquad \operatorname{Var}(X_t) = 1\ \ \forall t.

rng = np.random.default_rng(2)
theta = 1.0
sigma = np.sqrt(2 * theta)                   # the VP choice: sigma^2 = 2 theta
T, n, n_paths = 3.0, 300, 50000
t = np.linspace(0, T, n + 1)
x0 = rng.choice([-1.0, 1.0], size=n_paths)   # unit-variance, far from Gaussian
dW = np.sqrt(T / n) * rng.standard_normal((n_paths, n))
X = euler_maruyama(lambda x, s: -theta * x, lambda x, s: sigma, x0, t, dW)
print(f'{"t":>5} {"Var(X_t)":>9} {"E[X_t^4]":>9} {"3 - 2 alpha_t^2":>16}')
for s in [0.0, 0.25, 1.0, 3.0]:
    k = round(s / (T / n))
    alpha = np.exp(-2 * theta * s)
    print(f'{s:5.2f} {X[:, k].var():9.4f} {np.mean(X[:, k]**4):9.3f} '
          f'{3 - 2 * alpha**2:16.3f}')
    t  Var(X_t)  E[X_t^4]  3 - 2 alpha_t^2
 0.00    1.0000     1.000            1.000
 0.25    1.0030     2.286            2.264
 1.00    1.0073     2.999            2.963
 3.00    1.0139     3.062            3.000

This is exactly DDPM’s forward marginal in continuous time; the VP-SDE is dX=-\tfrac12\beta(t)X\,dt+\sqrt{\beta(t)}\,dW.

Recap

Wrap-up

  • Brownian motion: \mathcal N(0,t) increments, \operatorname{Cov}=\min(s,t), nowhere differentiable.
  • Quadratic variation (dW)^2=dt forces a new calculus.
  • Itô integral: zero mean, isometry; Itô’s lemma adds \tfrac12 g^2\phi_{xx}.
  • An SDE is a law over paths; Euler–Maruyama adds a \sqrt{\Delta t} kick; weak order 1.
  • OU mean-reverts to \mathcal N(0,\sigma^2/2\theta).
  • VP normalization (\sigma^2=2\theta) is DDPM’s forward marginal.

Next: track the whole density, the Fokker–Planck equation.