Stochastic Calculus in Machine Learning: Sampling and Generative Modelling

Contents

• Introduction and overview
• 0. Review: conditional expectation, Gaussian vectors and discrete martingales
• 1. Stochastic processes and Brownian motion
• 2. Continuous-time martingales
• 3. The Itô integral
• 4. Itô’s formula, stochastic differential equations and their PDEs
• 5. Girsanov’s theorem and changes of measure
• 6. Martingale problems, weak and strong solutions
• 7. Langevin dynamics and the unadjusted Langevin algorithm
• 8. Diffusion models: generative modelling by time reversal
• 9. Research directions and open problems
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: an M2 student with a solid M1 in probability who uses Langevin samplers or diffusion models and wants to know why they work, when they fail, and how to prove it.

Empirical first: the unadjusted Langevin algorithm run on a standard Gaussian with step 0.5 is stable and mixes in a few dozen iterations, yet its samples have variance 1.33 instead of 1, whatever the number of iterations (ch. 7). On a two-mode target the same algorithm never leaves its starting mode, while a diffusion model driven by the exact scores recovers both modes with the right weights (ch. 8). The course explains both observations with one toolbox: Itô calculus.

The course builds stochastic calculus from Brownian motion up to stochastic differential equations, their partial differential equations and their changes of measure, then uses it to analyse two families of modern algorithms: diffusion-based Markov chain Monte Carlo and score-based generative models. It shares its objects with Computational Statistics (Monte Carlo and MCMC), Introduction to Probabilistic Graphical Models and Deep Generative Models and Optimal Transport for Machine Learning (Wasserstein distances), and it prepares Bayesian Machine Learning, Generative Models for Images and Diffusion and Multiscale Models in Artificial Intelligence.

What you will be able to do

• Construct Brownian motion, prove its main invariances, and compute its quadratic variation.
• Use continuous-time martingales: optional stopping, Doob’s inequalities, convergence theorems and the quadratic variation of a local martingale.
• Define the Itô integral through the Itô isometry, and compute stochastic integrals with Itô’s formula.
• Prove existence and uniqueness of strong solutions of Lipschitz stochastic differential equations, and simulate them with the Euler–Maruyama scheme.
• Derive the Kolmogorov backward and Fokker–Planck equations of a diffusion, and find its invariant distribution.
• Change the drift of a diffusion with Girsanov’s theorem, write path-space likelihood ratios, and bound the Kullback–Leibler divergence between two diffusions.
• Distinguish weak and strong solutions, and formulate a diffusion as a solution of a martingale problem.
• Prove non-asymptotic convergence bounds for the overdamped Langevin diffusion and for the unadjusted Langevin algorithm, and explain the bias of the latter.
• Derive the reverse-time SDE of a diffusion model, train a score by denoising score matching, and bound the error of the generated distribution.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com presents stochastic calculus as an essential tool for understanding and justifying modern machine learning techniques, such as the unadjusted Langevin algorithm and diffusion models. The course is an introduction to stochastic calculus and its applications in machine learning, with a particular focus on sampling and generative modelling.

It has two parts. The first covers the fundamentals of Itô calculus, including diffusion processes and their connections to partial differential equations. The second applies stochastic calculus to the analysis of diffusion-based MCMC methods and generative models.

The description states three objectives: present the foundational concepts of stochastic calculus and Itô calculus; develop familiarity with stochastic differential equations and diffusion processes; and apply these tools to the design and analysis of modern methods such as diffusion-based MCMC and generative models. The announced learning outcomes repeat these three points, adding the relevance of stochastic calculus in modern statistics, machine learning and applied contexts.

The description then announces eight lectures, matched item by item below. The Université Paris-Saclay course detail gives the volume (24 h of lectures, 5 ECTS) and no programme items.

Beyond the syllabus, the page adds the reflection principle (ch. 1.3), the strong order of Euler–Maruyama (ch. 4.4), and, in self-study, current research on sampling and diffusion models (ch. 9).

Where this course fits in the programme

Official prerequisites (master-mva.com): “first year level knowledge in statistics and probability”. The material assumes M1 measure-theoretic probability (conditional expectation, Gaussian vectors, modes of convergence), discrete-time martingales, and basic multivariate calculus; ch. 0 restates what is used.

• Preparation — the refreshers Probability and Statistics for MVA, Linear Algebra, Analysis and Optimization for MVA (gradients, Hessians, convexity, used in ch. 7) and Scientific Python, NumPy and PyTorch for MVA (the code of ch. 1–8).
• Related courses in semester 1 — the course description names Introduction to Probabilistic Graphical Models and Deep Generative Models and Computational Statistics as complements; Optimal Transport for Machine Learning (Wasserstein distances, used in ch. 7), Convex Optimization and Applications in Machine Learning (ULA as a noisy gradient descent), Theoretical Foundations of Deep Learning and Mathematical Methods for Neuroscience (stochastic models of neurons) share tools with it.
• Courses that build on it in semester 2 — the description also names Bayesian Machine Learning and Generative Models for Images (“Modèles génératifs pour l’image”). Further uses: Diffusion and Multiscale Models in Artificial Intelligence, Inverse Problems and Imaging: Statistical and Stochastic Approaches (Langevin posterior sampling), Deep Learning for Image Restoration and Synthesis, Interactions: Complexity, Stochasticity, Extrema and Rare Events, Stopping Times and Online Algorithms and Learning Physical Phenomena: Metamodelling and Uncertainty Quantification.

Course outline

• 0. Review: conditional expectation, Gaussian vectors and discrete martingales — the M1 tools in the form used later.
• 1. Stochastic processes and Brownian motion — construction, invariances, reflection principle, quadratic variation and the roughness of paths.
• 2. Continuous-time martingales — stopping, Doob’s inequality, convergence, local martingales and the quadratic variation process.
• 3. The Itô integral — the isometry, its extension, integrals against continuous martingales and semimartingales.
• 4. Itô’s formula, stochastic differential equations and their PDEs — the chain rule of stochastic calculus, strong solutions, Euler–Maruyama, Kolmogorov and Fokker–Planck equations, models from physics and biology.
• 5. Girsanov’s theorem and changes of measure — changing the drift, likelihoods on path space, the KL divergence between two diffusions, importance sampling.
• 6. Martingale problems, weak and strong solutions — Tanaka’s equation, the Stroock–Varadhan formulation and the Yamada–Watanabe theorems.
• 7. Langevin dynamics and the unadjusted Langevin algorithm — exponential convergence of the diffusion, non-asymptotic bounds for ULA, its limitations and extensions.
• 8. Diffusion models: generative modelling by time reversal — forward noising, the reverse-time SDE, score matching and KL guarantees.
• 9. Research directions and open problems — sampling complexity, Wasserstein gradient flows, theory of diffusion models (self-study).

Session plan

24 h in class = 8 sessions of 3 h.

The volume is the Université Paris-Saclay programme table (24 h of lectures). The master-mva.com description announces eight lectures, so each session below is one official lecture; this split is a proposal.

The order follows the official one with a single swap: lecture 8 (Langevin dynamics and ULA) is taught before lecture 7 (diffusion models). The forward noising process of a diffusion model is a Langevin diffusion whose convergence is proved in ch. 7, and the comparison of the two samplers on a multimodal target needs ULA first. Level 1 = foundations, 2 = core results, 3 = advanced and research-level material; the level never decreases.

Self-study, outside the 24 h (about 3 h in total): ch. 9 Research directions and open problems — sampling as optimisation, kinetic Langevin and beyond log-concavity, theory of diffusion models, bridges and interpolants, learned samplers (≈ 3 h, with the open-ended parts of Ex. 16–17).

How to study this module

• Read ch. 0 before the first session: the later proofs cite its numbered results without restating them.
• After each session, redo the proofs of the boxed theorems without notes; the written-exam proposal asks for exactly these arguments.
• Run every code block of a chapter, then change one parameter (step size, horizon, drift) and predict the effect before running again.
• Do the exercises of the Practice column right after the matching session; ★★★ exercises prepare the projects.
• Python environment: NumPy and SciPy for ch. 1–8, PyTorch for one block of ch. 8; every block fixes its seed and runs in under 30 seconds on one CPU thread.

0. Review: conditional expectation, Gaussian vectors and discrete martingales

Stochastic calculus rests on three M1 tools: conditional expectation, Gaussian vectors and discrete-time martingales. This chapter states them in the form used later; background is in the refresher Probability and Statistics for MVA. All random variables live on a probability space (\Omega, \mathcal F, \mathbb P).

0.1 Conditional expectation

The rules used in every proof:

• Tower property: if \mathcal H \subset \mathcal G, then \mathbb E[\mathbb E[X \mid \mathcal G] \mid \mathcal H] = \mathbb E[X \mid \mathcal H].
• Taking out what is known: if Z is \mathcal G-measurable and ZX is integrable, then \mathbb E[ZX \mid \mathcal G] = Z\, \mathbb E[X \mid \mathcal G].
• Independence: if X is independent of \mathcal G, then \mathbb E[X \mid \mathcal G] = \mathbb E[X].
• Conditional Jensen: for convex \varphi with \varphi(X) integrable, \varphi(\mathbb E[X \mid \mathcal G]) \le \mathbb E[\varphi(X) \mid \mathcal G].

0.2 Gaussian vectors

A Gaussian process (X_t)_{t \ge 0} has Gaussian finite-dimensional marginals; its law is fixed by its mean and covariance functions. Fact 2 makes the Wiener integrals of ch. 3 Gaussian.

0.3 Discrete-time martingales

1. Stochastic processes and Brownian motion

Brownian motion drives every model of the course. This chapter constructs it, lists its symmetries, and proves the property that forces a new calculus: finite quadratic variation but infinite total variation.

1.1 Processes, filtrations and continuity

The finite-dimensional distributions of X are the laws of the vectors (X_{t_1}, \dots, X_{t_n}). Two theorems of Kolmogorov turn them into a process, then into a continuous one.

Proof idea. On the dyadic grid of mesh 2^{-n}, Markov’s inequality bounds the probability that some increment exceeds 2^{-n\gamma} by C\, 2^{-n(\beta - \alpha\gamma)}, which is summable. By Borel–Cantelli, the Hölder bound holds on the grid for large n; chaining extends it to [0,T] (Revuz and Yor, Theorem I.2.1).

1.2 Brownian motion

Donsker’s theorem is a constructive alternative: if (\xi_k) are i.i.d. with mean 0 and variance 1, the interpolated walk n^{-1/2} \sum_{k \le nt} \xi_k converges in law to B on C([0,1]).

1.3 Invariances, Markov property and reflection

A stopping time is a random time \tau \in [0, \infty] such that \{\tau \le t\} \in \mathcal F_t for every t: at each time, one knows whether it has already occurred. The hitting time \tau_a = \inf\{t : B_t = a\} is one.

The events known at time \tau form the σ-algebra \mathcal F_\tau of the events A such that A \cap \{\tau \le t\} \in \mathcal F_t for every t (ch. 2.1 studies both notions).

The strong Markov property extends the third item to stopping times: if \tau is an almost surely finite stopping time, then B_{\tau + t} - B_\tau is a Brownian motion independent of \mathcal F_\tau (Revuz and Yor, ch. III.3).

1.4 Path regularity: the quadratic variation

Brownian paths are continuous but rough. One number measures the roughness and governs Itô calculus.

Hence \int H_s\, dB_s cannot be a pathwise Stieltjes integral, and ch. 3 needs a probabilistic construction.

Proof sketch for nowhere differentiability. If the path had a derivative bounded by K at some s \in [0,1], then for large n three consecutive increments over intervals of length 1/n next to s would be at most 7K/n in absolute value.

Each increment has density at most \sqrt{n/2\pi}, so a fixed position has probability at most (14K/\sqrt{2\pi n})^3. A union bound over the n positions gives O(n^{-1/2}) \to 0.

import numpy as np

rng = np.random.default_rng(1)
T, n_max = 1.0, 2 ** 20
dB = rng.standard_normal(n_max) * np.sqrt(T / n_max)    # finest increments
B = np.concatenate([[0.0], np.cumsum(dB)])                # one Brownian path on [0, 1]
for k in [4, 8, 12, 16, 20]:
    incr = np.diff(B[:: n_max // 2 ** k])                 # increments on a dyadic grid of 2^k steps
    qv, tv = np.sum(incr ** 2), np.sum(np.abs(incr))
    print(f"2^{k:<2d} steps: quadratic variation {qv:.4f}   total variation {tv:7.1f}")

2. Continuous-time martingales

Stochastic integrals are martingales, controlled by martingale inequalities. This chapter transfers ch. 0.3 to continuous time and introduces the quadratic variation process.

From now on, (\mathcal F_t) is a filtration satisfying the usual conditions (it is right-continuous and contains the null sets), and B is an (\mathcal F_t)-Brownian motion.

2.1 Martingales and stopping times

import numpy as np

rng = np.random.default_rng(2)
a, b, dt, n_paths = 1.0, 2.0, 1e-3, 20_000
x, tau = np.zeros(n_paths), np.zeros(n_paths)
alive, t = np.ones(n_paths, dtype=bool), 0.0
while alive.any():                                        # Euler steps until every path has left (-a, b)
    t += dt
    x[alive] += np.sqrt(dt) * rng.standard_normal(alive.sum())
    out = alive & ((x <= -a) | (x >= b))
    tau[out], alive = t, alive & ~out
print(f"P(exit at b) = {np.mean(x >= b):.4f}   (exact a/(a+b) = {a / (a + b):.4f})")
print(f"E[tau]       = {tau.mean():.4f}   (exact ab = {a * b:.4f})")

2.2 Doob’s inequality and convergence

Proof idea. Doob’s upcrossing inequality along dyadic times, uniform in the grid, extends to the continuous path, which therefore cannot oscillate forever. Part 2 adds Theorem 2.6 to the L^2 argument of Theorem 0.5 (Revuz and Yor, ch. II.2).

2.3 Local martingales

Stochastic integrals of merely locally square-integrable integrands need not be martingales; localisation names them.

2.4 Quadratic variation and the Doob–Meyer decomposition

For Brownian motion, B^2 - t is a martingale and t is the quadratic variation. This holds for every continuous local martingale.

The bracket of two continuous local martingales is defined by polarisation: \langle M, N \rangle = \tfrac14 \big( \langle M + N \rangle - \langle M - N \rangle \big). It is the unique continuous finite-variation process with \langle M, N \rangle_0 = 0 such that MN - \langle M, N \rangle is a local martingale. For a d-dimensional Brownian motion, \langle B^i, B^j \rangle_t = \delta_{ij}\, t.

3. The Itô integral

We want to define \int_0^t H_s\, dB_s, the limit of \sum_i H_{t_i} (B_{t_{i+1}} - B_{t_i}), for random integrands H. Corollary 1.9 rules out a pathwise definition. Itô defines the integral as a limit in L^2(\mathbb P), using one exact identity: the Itô isometry.

3.1 Simple integrands and the isometry

3.2 Extension to square-integrable integrands

Localisation. Let M be a continuous local martingale and H a progressive process such that, almost surely for every t, \int_0^t H_s^2\, d\langle M \rangle_s < \infty . Stopping when this integral or |M| exceeds n reduces the construction to Theorem 3.4.

Then H \cdot M is a continuous local martingale with the same bracket; its mean need not be 0.

For deterministic f \in L^2([0,t]), the Wiener integral \int_0^t f(s)\, dB_s is an L^2 limit of Gaussian variables, hence (Proposition 0.3) \int_0^t f(s)\, dB_s \sim \mathcal N\Big(0, \int_0^t f(s)^2\, ds\Big).

3.3 Riemann sums and the role of the left end point

import numpy as np

rng = np.random.default_rng(3)
n, n_paths = 500, 5_000
dB = rng.standard_normal((n_paths, n)) * np.sqrt(1.0 / n)
B = np.hstack([np.zeros((n_paths, 1)), np.cumsum(dB, axis=1)])
ito = np.sum(B[:, :-1] * dB, axis=1)                      # left end points
strat = np.sum(0.5 * (B[:, :-1] + B[:, 1:]) * dB, axis=1) # average of the two end points
B1 = B[:, -1]
print(f"mean |Ito sum - (B_1^2 - 1)/2|   = {np.mean(np.abs(ito - (B1**2 - 1) / 2)):.4f}")
print(f"mean |Strat. sum - B_1^2 / 2|     = {np.mean(np.abs(strat - B1**2 / 2)):.1e}")
print(f"E[Ito integral] = {ito.mean():+.4f}   E[(Ito integral)^2] = {np.mean(ito**2):.4f}  (isometry: 0.5)")
print(f"E[Strat. integral] = {strat.mean():+.4f}")

3.4 Continuous semimartingales

Solutions of SDEs, X_t = X_0 + \int_0^t b_s\, ds + \int_0^t \sigma_s\, dB_s, are continuous semimartingales with \langle X \rangle_t = \int_0^t \sigma_s^2\, ds.

4. Itô’s formula, stochastic differential equations and their PDEs

This chapter gives the chain rule of stochastic calculus and uses it to solve SDEs, to link their laws to PDEs, and to analyse two classical models of physics and biology.

4.1 Itô’s formula

In differential notation, the formula is computed with the multiplication table dt \cdot dt = 0, \qquad dt \cdot dB^i_t = 0, \qquad dB^i_t \cdot dB^j_t = \delta_{ij}\, dt . For instance d(B_t^2) = 2 B_t\, dB_t + dt, which is Example 3.6 again.

4.2 Two applications

Lévy’s theorem recognises Brownian motions in ch. 5 and ch. 6. For moments, the Burkholder–Davis–Gundy inequalities extend Doob’s: for p > 0 and M_0 = 0, \mathbb E \sup_{s \le t} |M_s|^p is comparable to \mathbb E\langle M \rangle_t^{p/2} (Revuz and Yor, Theorem IV.4.1).

4.3 Stochastic differential equations: strong solutions

4.4 The Euler–Maruyama scheme

Simulation replaces the SDE by the recursion, for a step h and k \ge 0, X^h_{k+1} = X^h_k + h\, b(X^h_k) + \sigma(X^h_k) \big( B_{(k+1)h} - B_{kh} \big), \qquad X^h_0 = \xi .

The weak error is of order h: |\mathbb E f(X_T) - \mathbb E f(X^h_{T/h})| \le C h for smooth f (Talay and Tubaro, 1990). Sampling cares about the weak error; ch. 7 measures it on the invariant law.

import numpy as np

rng = np.random.default_rng(4)
mu, sigma, x0, T = 1.0, 0.8, 1.0, 1.0
n_paths, n_fine = 5_000, 2 ** 10
dW = rng.standard_normal((n_paths, n_fine)) * np.sqrt(T / n_fine)
exact = x0 * np.exp((mu - sigma**2 / 2) * T + sigma * dW.sum(axis=1))   # closed-form GBM at time T
hs, errs = [], []
for k in [2, 4, 6, 8, 10]:
    m = 2 ** k
    incr = dW.reshape(n_paths, m, -1).sum(axis=2)       # the same Brownian path on a coarser grid
    x = np.full(n_paths, x0)
    for j in range(m):                                  # Euler-Maruyama
        x = x + mu * x * (T / m) + sigma * x * incr[:, j]
    hs.append(T / m)
    errs.append(np.sqrt(np.mean((x - exact) ** 2)))
    print(f"h = 2^-{k:<2d}  strong L2 error at T = {errs[-1]:.4f}")
print(f"fitted order: {np.polyfit(np.log(hs), np.log(errs), 1)[0]:.2f}")

4.5 Generator, Kolmogorov and Fokker–Planck equations

Here X solves dX_t = b(X_t)\, dt + \sigma(X_t)\, dB_t with Lipschitz coefficients, a = \sigma \sigma^\top, and X^x starts at x.

With a killing rate V \ge 0, the same proof gives the Feynman–Kac formula: u(t,x) = \mathbb E\big[g(X^x_t)\, e^{-\int_0^t V(X^x_s)\, ds}\big] solves \partial_t u = \mathcal L u - V u.

Ch. 7 starts here: to sample from \pi \propto e^{-U}, simulate this diffusion.

4.6 Diffusions in physics and biology

The course description stresses models of physics and biology; two classical ones use exactly the tools above.

5. Girsanov’s theorem and changes of measure

Changing the drift of a diffusion amounts to reweighting its paths, and Girsanov’s theorem gives the weight. It yields likelihoods for diffusions, importance samplers, and the Kullback–Leibler formula between two diffusions behind the guarantees of ch. 7 and ch. 8.

5.1 Warm-up: shifting a Gaussian

Let Z \sim \mathcal N(0,1) under \mathbb P and define \mathbb Q by d\mathbb Q / d\mathbb P = e^{\theta Z - \theta^2/2}. For every u, \mathbb E_{\mathbb Q}\, e^{iuZ} = e^{-\theta^2/2}\, \mathbb E\, e^{(\theta + iu) Z} = e^{iu\theta - u^2/2}, so Z \sim \mathcal N(\theta, 1) under \mathbb Q. Girsanov’s theorem does this for every increment of a Brownian path, with a shift that may depend on the past.

5.2 Girsanov’s theorem

A condition is needed: with \theta \equiv 1 on [0, \infty), Z_t \to 0 almost surely (Example 2.8), and Z_\infty = 0 is no density.

5.3 Radon–Nikodym derivatives between diffusion laws

Let \mathcal C = C([0,T], \mathbb R^d) with its canonical process \omega_t. For i = 1, 2, let \mathbb P^i be the law on \mathcal C of the solution of dX_t = b_i(t, X)\, dt + dB_t, \qquad X_0 \sim \nu , where the drift b_i(t, \omega) may depend on the path up to time t (it is progressive). Let \mathbb W^\nu be the law of \xi + B with \xi \sim \nu independent of B.

Path-dependent drifts matter: the interpolation of an Euler scheme has drift b(\omega_{kh}) on [kh, (k+1)h), so Theorem 5.3 compares an algorithm with the diffusion it discretises (ch. 7 and ch. 8).

The inequality on the left is data processing: mapping a path to its value at time T cannot increase the divergence. A drift error along the path thus bounds the error of the final sample.

5.4 Likelihood ratios and importance sampling

The diffusion coefficient is different. For dX = \sigma\, dB, one observed path on [0,T] has quadratic variation exactly \sigma^2 T (Theorem 1.8). The laws for two values of \sigma are mutually singular: no likelihood ratio exists, and \sigma is read off a single path without error.

import numpy as np
from scipy.stats import norm

rng = np.random.default_rng(5)
a, n, n_paths = 4.0, 200, 20_000
dt = 1.0 / n
# Crude Monte Carlo of P(B_1 > a): about 0.6 successes expected
crude = np.mean(rng.standard_normal(n_paths) > a)
# Importance sampling: simulate under Q, where B_t = W_t + a t, and reweight by dP/dQ
dB = a * dt + np.sqrt(dt) * rng.standard_normal((n_paths, n))
B = np.hstack([np.zeros((n_paths, 1)), np.cumsum(dB, axis=1)])
w = np.exp(-a * B[:, -1] + a**2 / 2)                       # Girsanov weight on F_1
x0, x1 = B[:, :-1], B[:, 1:]
p_cross = np.where((x0 < a) & (x1 < a), np.exp(-2 * (a - x0) * (a - x1) / dt), 1.0)
hit = 1 - np.prod(1 - p_cross, axis=1)                      # P(the bridge reaches a | grid values)
for name, est, exact in [("P(B_1 > a)      ", w * (B[:, -1] > a), norm.sf(a)),
                         ("P(max B_t >= a) ", w * hit, 2 * norm.sf(a))]:
    m, se = est.mean(), est.std() / np.sqrt(n_paths)
    print(f"{name} IS = {m:.4e} +/- {1.96 * se:.1e}   exact = {exact:.4e}")
print(f"crude Monte Carlo of P(B_1 > a) with {n_paths} draws: {crude:.1e}")

The zero-variance change of measure is the conditional law given the event: a diffusion whose extra drift is the gradient of the log-probability of the event given the current state (a Doob h-transform). The reverse drift of diffusion models (ch. 8) has the same structure.

6. Martingale problems, weak and strong solutions

Theorem 4.5 builds a solution from a given Brownian motion, but some equations have no such solution. This chapter separates what a solution is from how it is built: weak solutions, the martingale problem of Stroock and Varadhan, and the Yamada–Watanabe theorems.

6.1 Weak solutions and Tanaka’s equation

Weak existence and uniqueness in law can thus hold without strong existence; for sampling, only the law matters.

6.2 The martingale problem

By Dynkin’s formula (Proposition 4.10), f(X_t) - \int_0^t \mathcal L f(X_s)\, ds is a martingale. Stroock and Varadhan made this the definition of a diffusion with generator \mathcal L.

The formulation pays off in three ways.

• Existence by compactness. Laws of approximations (Euler schemes) are tight on path space, and limit points solve the martingale problem.
• Uniqueness gives the Markov property. \mathrm{MP}(b, a) is well posed when a is bounded, continuous and uniformly elliptic and b is bounded measurable (Stroock and Varadhan, 1979); its solutions form a strong Markov family.
• Diffusion limits of algorithms. If the generators of rescaled Markov chains converge to \mathcal L, their laws converge to the solution of \mathrm{MP}(b,a); random-walk Metropolis in high dimension thus behaves like a Langevin diffusion (Roberts, Gelman and Gilks, 1997).

6.3 From weak to strong: the Yamada–Watanabe theorems

Proof idea. Put two weak solutions on one space with a common Brownian motion and initial value, conditionally independent given these. Pathwise uniqueness makes them equal, and a process conditionally independent of itself given (X_0, B) is a function of (X_0, B) (Revuz and Yor, ch. IX.1).

Noise can also restore uniqueness. For dX = b(X)\, dt + dB with b bounded measurable, \dot x = b(x) may have many solutions, yet the SDE has a unique strong solution (Zvonkin, 1974, in dimension one; Veretennikov, 1981). Weak existence already follows from Girsanov’s theorem with \theta = b(X).

7. Langevin dynamics and the unadjusted Langevin algorithm

We want samples from a density \pi(x) \propto e^{-U(x)} on \mathbb R^d known up to its normalising constant, such as a Bayesian posterior. This chapter proves that the Langevin diffusion of Corollary 4.13 converges to \pi, discretises it into the unadjusted Langevin algorithm (ULA), and prices the discretisation.

7.1 The sampling problem and the Langevin diffusion

In Bayesian logistic regression with a Gaussian prior of precision \lambda, U(x) = \lambda |x|^2 / 2 + \sum_i \log(1 + e^{-y_i \langle a_i, x \rangle}): the normalising constant is a d-dimensional integral, but \nabla U is cheap.

The overdamped Langevin diffusion is dX_t = -\nabla U(X_t)\, dt + \sqrt 2\, dB_t . It leaves \pi invariant (Corollary 4.13) and, for Lipschitz \nabla U, has a unique strong solution (Theorem 4.5). Write \mu P_t for the law of X_t when X_0 \sim \mu.

Assumption (A). U \in C^2(\mathbb R^d) and, for some 0 < m \le L, its Hessian satisfies m I \preceq \nabla^2 U(x) \preceq L I for all x.

The condition number is \kappa = L / m, and x^\star denotes the minimiser of U.

In the logistic example, m \ge \lambda and L \le \lambda + \frac14 \|A\|^2, with A the matrix of rows a_i. Laws are compared in Wasserstein distance, W_2(\mu, \nu) = \inf \big\{ (\mathbb E|X - Y|^2)^{1/2} : X \sim \mu,\ Y \sim \nu \big\}, an infimum over couplings, which is attained.

7.2 Exponential convergence of the diffusion

Convergence in entropy needs a functional inequality. With the relative Fisher information I(\nu \| \pi) = \int |\nabla \log (d\nu / d\pi)|^2\, d\nu, \pi satisfies a logarithmic Sobolev inequality with constant C_{\mathrm{LS}} if, for all smooth \nu \ll \pi, \mathrm{KL}(\nu \| \pi) \le \frac{C_{\mathrm{LS}}}{2}\, I(\nu \| \pi). By the Bakry–Émery criterion, (A) gives C_{\mathrm{LS}} = 1/m; for the standard Gaussian, C_{\mathrm{LS}} = 1 (Bakry, Gentil and Ledoux, 2014, ch. 5).

Log-Sobolev inequalities survive bounded perturbations of U (Holley and Stroock, 1987), so Theorem 7.2 covers some non-convex potentials. But with two modes separated by a barrier of height \Delta U, C_{\mathrm{LS}} grows like e^{\Delta U}; ch. 8 shows the consequence.

7.3 The unadjusted Langevin algorithm

import numpy as np

rng = np.random.default_rng(6)
U = lambda x: 0.5 * x**2                                    # target N(0, 1)
grad_U = lambda x: x
n_chains, n_iter = 20_000, 400
for gamma in [0.1, 0.5, 1.0]:
    x_ula = np.zeros(n_chains)
    x_mala = np.zeros(n_chains)
    for _ in range(n_iter):
        x_ula = x_ula - gamma * grad_U(x_ula) + np.sqrt(2 * gamma) * rng.standard_normal(n_chains)
        y = x_mala - gamma * grad_U(x_mala) + np.sqrt(2 * gamma) * rng.standard_normal(n_chains)
        log_q_fwd = -(y - x_mala + gamma * grad_U(x_mala)) ** 2 / (4 * gamma)
        log_q_bwd = -(x_mala - y + gamma * grad_U(y)) ** 2 / (4 * gamma)
        log_r = -U(y) + U(x_mala) + log_q_bwd - log_q_fwd   # Metropolis-Hastings correction
        accept = np.log(rng.uniform(size=n_chains)) < log_r
        x_mala = np.where(accept, y, x_mala)
    print(f"gamma = {gamma:.1f}: Var ULA = {x_ula.var():.3f} (theory {1 / (1 - gamma / 2):.3f})"
          f"   Var MALA = {x_mala.var():.3f}")

Such bounds are due to Dalalyan (2017) and Durmus and Moulines (2017, 2019); Dalalyan and Karagulyan (2019) get the constant 1.65 for \gamma \le 2/(m + L). The bias term is of order \gamma^{1/2}, while Example 7.4 shows order \gamma; with a Lipschitz Hessian, a finer analysis recovers order \gamma (Durmus and Moulines, 2019).

Girsanov’s theorem gives a second, information-theoretic route, reused in ch. 8.

With Pinsker’s inequality and Theorem 7.2, this bound of order L^2 T d\gamma gives total-variation guarantees for ULA, as in Dalalyan (2017).

7.4 Limitations

The bias. At fixed step, ULA targets \pi_\gamma \ne \pi. Remedies: decreasing steps (Durmus and Moulines, 2017) or a Metropolis correction.

MALA is exact, but its step must shrink like d^{-1/3} from a stationary start (Roberts and Rosenthal, 1998).

Instability. The globally Lipschitz gradient of (A) is not a technicality.

Tamed ULA replaces \gamma \nabla U by \gamma \nabla U / (1 + \gamma |\nabla U|) and is stable (Brosse, Durmus, Moulines and Sabanis, 2019); for superlinear coefficients, Euler moments diverge even as h \to 0 (Hutzenthaler, Jentzen and Kloeden, 2011).

Multimodality. With separated modes, the log-Sobolev constant and the mixing times are exponentially large (ch. 7.2); ch. 8 sidesteps this by learning.

7.5 Extensions of overdamped Langevin dynamics

• Kinetic Langevin: dX = V\, dt, dV = -\nabla U(X)\, dt - \lambda V\, dt + \sqrt{2\lambda}\, dB leaves \pi \otimes \mathcal N(0, I_d) invariant (Exercise 17; ch. 9.2).
• Stochastic gradients: a mini-batch estimate of \nabla U gives SGLD (Welling and Teh, 2011).
• Non-smooth potentials: a Moreau–Yosida envelope gives MYULA (Durmus, Moulines and Pereyra, 2018).
• Preconditioning: dX = -M \nabla U(X)\, dt + \sqrt{2M}\, dB still leaves \pi invariant for a constant matrix M \succ 0.

8. Diffusion models: generative modelling by time reversal

A diffusion model learns to sample from a data distribution \mu known through samples. It noises the data with a simple diffusion, then runs it backwards in time, which needs the gradient of the log-density of the noisy data: the score. This chapter derives the reverse equation, learns the score, and bounds the error with Girsanov’s theorem.

8.1 The forward noising process

Let \mu be a probability measure on \mathbb R^d with a finite second moment. The forward process is the Ornstein–Uhlenbeck diffusion dX_t = -X_t\, dt + \sqrt 2\, dB_t, \qquad X_0 \sim \mu, \qquad t \in [0, T].

By Example 4.6, with \alpha_t = e^{-t} and \sigma_t^2 = 1 - e^{-2t}, X_t = \alpha_t X_0 + \sigma_t Z, \qquad Z \sim \mathcal N(0, I_d) \text{ independent of } X_0 . For t > 0, the law p_t of X_t has a smooth positive density (a Gaussian convolution). The forward process is the Langevin diffusion with U(x) = |x|^2/2 and invariant law \gamma_d = \mathcal N(0, I_d), whose log-Sobolev constant is 1. Theorem 7.2 gives, for 0 < t_0 < T, \mathrm{KL}(p_T \| \gamma_d) \le e^{-2(T - t_0)}\, \mathrm{KL}(p_{t_0} \| \gamma_d) . The noising forgets the data exponentially fast, whatever \mu. The variance-preserving SDE of Song et al. (2021), dX = -\frac12 \beta(t) X\, dt + \sqrt{\beta(t)}\, dB, is a time change of this process.

8.2 Time reversal

8.3 Learning the score

The reverse SDE needs the unknown \nabla \log p_t. Two identities make it learnable from samples.

Part 2 replaces an unknown regression target by a known one: a network s_\theta(t, x) minimises \int_{t_0}^T \lambda(t)\, \mathbb E \Big| s_\theta(t, X_t) + \frac{Z}{\sigma_t} \Big|^2 dt over samples of X_0, Z and t. With \lambda(t) = \sigma_t^2 and s_\theta = -\varepsilon_\theta / \sigma_t, this is the noise-prediction loss \mathbb E|\varepsilon_\theta(t, X_t) - Z|^2 of Ho, Jain and Abbeel (2020).

Memorisation. With the empirical measure of n training points and its exact score, p_t is a mixture of n Gaussians and the reverse process ends on the training points. A diffusion model generalises only because its learned score is not the exact empirical score (ch. 9).

8.4 Sampling guarantees

The algorithm replaces p_T by \gamma_d and the score by an estimate \hat s, and stops at a small t_0 > 0, since the score may blow up as t \to 0 for data near a low-dimensional set: d\hat Y_t = \big( \hat Y_t + 2\, \hat s(T - t, \hat Y_t) \big)\, dt + \sqrt 2\, dW_t, \qquad \hat Y_0 \sim \gamma_d, \qquad t \le T - t_0 . Let \hat q be the law of \hat Y_{T - t_0}, the output.

Three remarks put the theorem in context.

• The score term is the denoising score-matching loss with weight \lambda \equiv 1, up to C_t: training minimises a bound on the KL error of the output.
• Nothing is assumed on the shape of \mu: no log-concavity, no smoothness. Compared with ch. 7, the hardness has moved into learning the score.
• Discretisation adds a third term, comparing the score along the path with its values at grid times. Chen et al. (2023) bound it for Lipschitz scores; Benton, De Bortoli, Doucet and Deligiannidis (2024) need a number of steps nearly linear in d; Conforti, Durmus and Gentiloni Silveri (2025) prove KL guarantees under minimal data assumptions.

import numpy as np

rng = np.random.default_rng(7)
w, m, s = np.array([0.2, 0.8]), np.array([-3.0, 3.0]), np.array([0.5, 0.5])   # data: Gaussian mixture

def score(t, x):
    """Exact score of p_t, the law of e^{-t} X_0 + sqrt(1 - e^{-2t}) Z (forward OU)."""
    al = np.exp(-t)
    v = al**2 * s**2 + 1 - al**2                              # component variances at time t
    logc = np.log(w) - 0.5 * np.log(v) - (x[:, None] - al * m) ** 2 / (2 * v)
    r = np.exp(logc - logc.max(axis=1, keepdims=True))
    r /= r.sum(axis=1, keepdims=True)                         # responsibilities
    return np.sum(r * (al * m - x[:, None]) / v, axis=1)

n, T, N = 10_000, 6.0, 400
h = T / N
y_sde = rng.standard_normal(n)                                # start from N(0, 1), not p_T
y_ode = y_sde.copy()
for k in range(N):
    t = T - k * h                                             # forward time seen by the reverse step
    y_sde = y_sde + h * (y_sde + 2 * score(t, y_sde)) + np.sqrt(2 * h) * rng.standard_normal(n)
    y_ode = y_ode + h * (y_ode + score(t, y_ode))             # probability-flow ODE, Euler step
x_ula = np.full(n, -3.0)                                      # ULA on U = -log p_0, started in the left mode
for _ in range(2_000):
    x_ula = x_ula + 0.01 * score(0.0, x_ula) + np.sqrt(0.02) * rng.standard_normal(n)
for name, y in [("reverse SDE", y_sde), ("prob. flow ODE", y_ode), ("ULA", x_ula)]:
    print(f"{name:15s} weight of the right mode = {np.mean(y > 0):.3f}   "
          f"spread of |y| around 3 = {np.abs(y).std():.3f}")

import numpy as np
import torch

torch.manual_seed(0)
torch.set_num_threads(1)
means = torch.tensor([[0.0, 2.0], [-1.7, -1.0], [1.7, -1.0]])
weights = torch.tensor([0.5, 0.3, 0.2])

def sample_data(n):                                       # three-component Gaussian mixture in 2D
    idx = torch.multinomial(weights, n, replacement=True)
    return means[idx] + 0.3 * torch.randn(n, 2)

net = torch.nn.Sequential(torch.nn.Linear(3, 128), torch.nn.SiLU(),
                          torch.nn.Linear(128, 128), torch.nn.SiLU(), torch.nn.Linear(128, 2))
opt = torch.optim.Adam(net.parameters(), lr=2e-3)
T, delta = 5.0, 1e-3
for step in range(3_000):                                 # denoising score matching, noise-prediction form
    x0 = sample_data(512)
    t = delta + (T - delta) * torch.rand(512, 1)
    z = torch.randn_like(x0)
    xt = torch.exp(-t) * x0 + torch.sqrt(1 - torch.exp(-2 * t)) * z
    loss = ((net(torch.cat([xt, t], dim=1)) - z) ** 2).sum(dim=1).mean()
    opt.zero_grad()
    loss.backward()
    opt.step()

@torch.no_grad()
def score(t, x):                                          # s_theta(t, x) = -eps_theta(t, x) / sigma_t
    tt = torch.full((x.shape[0], 1), t)
    return -net(torch.cat([x, tt], dim=1)) / np.sqrt(1 - np.exp(-2 * t))

y, N = torch.randn(5_000, 2), 500
ts = np.linspace(T, delta, N + 1)
for k in range(N):                                        # Euler-Maruyama on the reverse SDE
    h = ts[k] - ts[k + 1]
    y = y + h * (y + 2 * score(ts[k], y)) + np.sqrt(2 * h) * torch.randn_like(y)
label = torch.cdist(y, means).argmin(dim=1)
print("final training loss:", round(loss.item(), 3))
print("mode weights:", np.round(torch.bincount(label, minlength=3).numpy() / len(y), 3), "target [0.5 0.3 0.2]")
print("within-mode std:", round((y - means[label]).std().item(), 3), "target 0.3")

9. Research directions and open problems

This chapter is for self-study. It surveys active questions around the two applications of the course, with entry points into the literature. The works cited are listed in the References.

9.1 Sampling as optimisation in the space of measures

Jordan, Kinderlehrer and Otto (1998) showed that the Fokker–Planck equation of the Langevin diffusion is the gradient flow of the functional \nu \mapsto \mathrm{KL}(\nu \| \pi) for the Wasserstein distance W_2. Theorem 7.2 is then the statement that this functional decreases along its own gradient flow.

This viewpoint imports convex optimisation into sampling. Wibisono (2018) reads ULA as a forward–flow splitting of this gradient flow, which explains its bias. Durmus, Majewski and Miasojedow (2019) analyse Langevin Monte Carlo as a stochastic convex optimisation algorithm.

Open questions. What is the optimal number of gradient queries to sample a strongly log-concave target to accuracy \varepsilon, and how does it depend on d and \kappa? Upper bounds have improved quickly, while matching lower bounds are known only in restricted settings.

9.2 Beyond overdamped Langevin

• Kinetic Langevin dynamics. The velocity variable makes the generator degenerate, and convergence relies on hypocoercivity (Villani, 2009). Eberle, Guillin and Zimmer (2019) give quantitative contraction rates by couplings. Cheng, Chatterji, Bartlett and Jordan (2018) show that a discretisation reaches W_2 accuracy \varepsilon in O(\sqrt d / \varepsilon) steps under (A), against O(d / \varepsilon^2) for ULA.
• Better discretisations. The randomised midpoint method of Shen and Lee (2019) improves the dependence on d further for kinetic Langevin dynamics.
• Beyond log-concavity. Vempala and Wibisono (2019) prove KL convergence of ULA under a log-Sobolev inequality and a Lipschitz gradient, with no convexity, which covers bounded perturbations of convex potentials (Holley and Stroock, 1987). For targets with well-separated modes, all local samplers remain exponentially slow; tempering and annealing are the classical answers.
• Stochastic and inexact gradients. Stochastic gradient Langevin dynamics (Welling and Teh, 2011) uses mini-batches; Dalalyan and Karagulyan (2019) quantify the effect of inexact gradients on ULA. Proximal versions (Pereyra, 2016; Durmus, Moulines and Pereyra, 2018) handle the non-smooth priors of imaging.

9.3 Theory of diffusion models

• Convergence. After Lee, Lu and Tan (2022), the analyses of Chen et al. (2023) and Chen, Lee and Lu (2023) showed that, given an L^2-accurate score, sampling needs no log-concavity of the data. Benton, De Bortoli, Doucet and Deligiannidis (2024) obtain a nearly linear dependence on the dimension; Conforti, Durmus and Gentiloni Silveri (2025) prove KL guarantees under minimal assumptions on the data.
• Low-dimensional data. De Bortoli (2022) studies convergence under the manifold hypothesis, where the score blows up as t \to 0 and early stopping is necessary.
• Statistics. Oko, Akiyama and Suzuki (2023) prove that diffusion models with neural-network score estimators attain minimax rates for distribution estimation in Besov classes.
• Deterministic samplers and design. The probability-flow ODE underlies DDIM (Song, Meng and Ermon, 2021); Karras, Aittala, Aila and Laine (2022) study the design space of noise schedules and solvers.

Open questions. Why do trained models generalise instead of memorising their training set (ch. 8.3)? How should the score be estimated when d is large but the data have low intrinsic dimension? Which guarantees hold for deterministic samplers, which have no noise to correct their errors?

9.4 Bridges, flows and interpolants

Diffusion models transport a Gaussian to the data along a fixed noising path. Other paths are possible.

• Schrödinger bridges find the diffusion closest in KL to a reference process with prescribed initial and final laws (Léonard, 2014). De Bortoli, Thornton, Heng and Doucet (2021) solve them by iterated diffusion models.
• Flow matching and stochastic interpolants learn the velocity field of a chosen interpolation between noise and data by regression, as in denoising score matching (Lipman et al., 2023; Albergo and Vanden-Eijnden, 2023; Liu, Gong and Liu, 2023).

Open question. Which interpolation minimises the total error, including the learning error, for a given computational budget?

9.5 Sampling with learned diffusions

• Unnormalised targets. When only e^{-U} is available, no data exist to learn a score. The path integral sampler (Zhang and Chen, 2022) and denoising diffusion samplers (Vargas, Grathwohl and Doucet, 2023) learn the drift by minimising a path-space KL divergence, computed with Girsanov’s theorem (ch. 5).
• Posterior sampling in inverse problems. A diffusion model trained on images is a prior; sampling the posterior given measurements requires the score of a conditional law. Approximations such as diffusion posterior sampling (Chung et al., 2023) work well in practice, but exact posterior sampling with a diffusion prior has no general efficient algorithm. See Inverse Problems and Imaging and Diffusion and Multiscale Models in Artificial Intelligence.

Exercises

The exercises are sorted by level: ★ applies the course, ★★ combines results or proves a short statement, ★★★ is open-ended or research-flavoured. Exercises 6, 9, 14, 16 and 17 are computational.

1. ★ The Brownian bridge. Let \beta_t = B_t - t B_1 for t \in [0,1]. Show that \beta is a centred Gaussian process with \mathrm{Cov}(\beta_s, \beta_t) = s(1-t) for s \le t, and that it is independent of B_1.

   Solution. For s \le t, \mathrm{Cov}(\beta_s, \beta_t) = s - ts - st + st = s(1 - t). Also \mathrm{Cov}(\beta_t, B_1) = t - t = 0. The vector (\beta_{t_1}, \dots, \beta_{t_n}, B_1) is Gaussian, so zero covariance means independence (Proposition 0.3). Hence \beta is Brownian motion conditioned on B_1 = 0, the bridge used in Example 5.6.

2. ★ Exit from an asymmetric interval. Brownian motion starts at 0, and \tau is its exit time from (-1, 3). Compute \mathbb P(B_\tau = 3) and \mathbb E\tau. Generalise to a start at x \in (-a, b).

   Solution. Example 2.5 with a = 1, b = 3 gives \mathbb P(B_\tau = 3) = \tfrac14, \qquad \mathbb E\tau = 3 . From x, the same two martingales give \mathbb P(B_\tau = b) = (x + a)/(a + b) and \mathbb E\tau = (x + a)(b - x).

3. ★ Itô’s formula in practice. (a) Compute d(B_t^3), d(t B_t) and d(e^{B_t}). (b) Deduce \mathbb E B_t^4 = 3t^2 from Itô’s formula. (c) Show that \int_0^t s\, dB_s = t B_t - \int_0^t B_s\, ds, and give its law.

   Solution. (a) d(B^3) = 3B^2 dB + 3B\, dt, d(tB) = B\, dt + t\, dB, and d(e^B) = e^B dB + \frac12 e^B dt.

   (b) d(B^4) = 4B^3 dB + 6B^2 dt. The dB integral is a true martingale, since \mathbb E \int_0^t B_s^6\, ds < \infty. Hence \mathbb E B_t^4 = 6 \int_0^t s\, ds = 3t^2.

   (c) Integrate d(tB_t). The left side is a Wiener integral, so \int_0^t s\, dB_s \sim \mathcal N(0, t^3/3).

4. ★ Covariance of the Ornstein–Uhlenbeck process. Let dX = -\theta X\, dt + \sigma\, dB with X_0 \sim \mathcal N(0, \sigma^2/(2\theta)) independent of B. Show that \mathrm{Cov}(X_s, X_t) = \frac{\sigma^2}{2\theta} e^{-\theta|t-s|}. Deduce the variance of the time average \frac1T \int_0^T X_t\, dt for large T.

   Solution. For s \le t, X_t = e^{-\theta(t-s)} X_s + \sigma \int_s^t e^{-\theta(t-r)} dB_r, and the integral is independent of \mathcal F_s. So \mathrm{Cov}(X_s, X_t) = e^{-\theta(t-s)} \mathrm{Var}(X_s). Then \mathrm{Var}\Big(\frac1T \int_0^T X_t\, dt\Big) = \frac{2}{T^2} \int_0^T (T - u)\, \frac{\sigma^2}{2\theta} e^{-\theta u}\, du \sim \frac{\sigma^2}{\theta^2 T}.

5. ★ Invariant densities in dimension one. Let dX = b(X)\, dt + \sigma(X)\, dB with \sigma > 0. (a) Show that a density with zero probability flux, -b\pi + \frac12 (\sigma^2 \pi)' = 0, satisfies \mathcal L^* \pi = 0, and solve for \pi.

   (b) Apply this to the Ornstein–Uhlenbeck process and to the Cox–Ingersoll–Ross process dX = (\alpha - \beta X)\, dt + c\sqrt X\, dB on (0, \infty), with \alpha, \beta, c > 0.

   Solution. (a) \mathcal L^* \pi is the derivative of the flux. With g = \sigma^2 \pi, the equation reads g' = (2b/\sigma^2)\, g, so \pi(x) \propto \frac{1}{\sigma(x)^2} \exp\Big( \int_0^x \frac{2 b(y)}{\sigma(y)^2}\, dy \Big). (b) For the Ornstein–Uhlenbeck process this gives e^{-\theta x^2/\sigma^2}, the law \mathcal N(0, \sigma^2/(2\theta)). For the CIR process it gives x^{2\alpha/c^2 - 1} e^{-2\beta x / c^2}: a Gamma law with shape 2\alpha/c^2 and rate 2\beta/c^2.

6. ★ ULA on a high-dimensional Gaussian (Python). Let \pi = \mathcal N(0, I_d) with d = 100. Run ULA from 0 with \gamma \in \{0.05, 0.1, 0.2\} for 300 iterations on 2 000 parallel chains. Estimate the stationary variance v per coordinate, compute W_2 = \sqrt d\, |\sqrt v - 1|, and compare with Example 7.4 and Theorem 7.5.

   Hints and expected result. Here m = L = 1 and v = 1/(1 - \gamma/2) = 1.026, 1.053, 1.111. The distances are W_2 = 0.127,\ 0.260,\ 0.541, while the bound 1.5 \sqrt{\gamma d} of Theorem 7.5 gives 3.35, 4.74 and 6.71. On Gaussians the bound is pessimistic by an order of magnitude, because the true bias is of order \gamma, not \gamma^{1/2}.

7. ★★ A strict local martingale. Let B be a Brownian motion in \mathbb R^3 started at x \ne 0; admit that it never hits 0.

   (a) Show that M_t = 1/|B_t| is a continuous local martingale. (b) Show that \mathbb E M_t \to 0, so that M is not a martingale.

   Solution. (a) f(y) = 1/|y| is harmonic on \mathbb R^3 \setminus \{0\}. Itô’s formula, stopped at the hitting times of small spheres around 0 (which tend to infinity), gives M_t = M_0 + \int_0^t \nabla f(B_s) \cdot dB_s.

   (b) For R > 0, split according to |B_t| \ge R or not. The density of B_t is at most (2\pi t)^{-3/2} and \int_{|y| < R} |y|^{-1} dy = 2\pi R^2, so \mathbb E M_t \le 1/R + 2\pi R^2 (2\pi t)^{-3/2}.

   With R = \sqrt t, this is O(t^{-1/2}), while M_0 = 1/|x|. As in Definition 2.9, M is a positive supermartingale, but its mean decreases strictly.

8. ★★ Kullback–Leibler divergence between two Ornstein–Uhlenbeck processes. Let \mathbb P^i be the law on [0, T] of dX = -\theta_i X\, dt + dB, both started from \mathcal N(0, 1/(2\theta_1)). Show that \mathrm{KL}(\mathbb P^1 \| \mathbb P^2) = \frac{(\theta_1 - \theta_2)^2\, T}{4\theta_1}. For \theta_1 = 1 and \theta_2 = 1.1, how long must the path be for the divergence to reach 1?

   Solution. By Corollary 5.4, the divergence is \frac12 \mathbb E_{\mathbb P^1} \int_0^T (\theta_1 - \theta_2)^2 X_t^2\, dt, and \mathbb E X_t^2 = 1/(2\theta_1) by stationarity under \mathbb P^1. Here the divergence is 0.0025\, T, which reaches 1 at T = 400. Distinguishing two drifts at distance \delta takes a time of order 1/\delta^2, in line with Exercise 9.

9. ★★ Maximum likelihood for the drift of an Ornstein–Uhlenbeck process (Python). Simulate dX = -\theta X\, dt + dB with \theta = 2, from stationarity, with Euler step 0.01 on [0, T] for T \in \{10, 100, 1000\}, with 200 replications.

   Compute \hat\theta_T of Example 5.5 with sums in place of integrals, compare its spread with \sqrt{2\theta/T}, and prove consistency.

   Hints and expected result. \hat\theta_T - \theta = -M_T / \langle M \rangle_T with M = \int X\, dB. The ergodic theorem gives \langle M \rangle_T / T \to 1/(2\theta), and M_T / \langle M \rangle_T \to 0 almost surely when \langle M \rangle_\infty = \infty (strong law of large numbers for martingales).

   A run gives empirical means 2.17, 2.02 and 2.01, and standard deviations 0.68, 0.197 and 0.063, against \sqrt{2\theta/T} = 0.632, 0.200, 0.063. The excess at T = 10 is a finite-sample bias of the ratio estimator.

10. ★★ Cameron–Martin and small balls. With h(t) = t on [0, 1], use Girsanov’s theorem to show that \frac{\mathbb P(\sup_{t \le 1} |B_t - t| < \varepsilon)}{\mathbb P(\sup_{t \le 1} |B_t| < \varepsilon)} \longrightarrow e^{-1/2} \qquad (\varepsilon \to 0).

    Solution. Under d\mathbb Q / d\mathbb P = e^{B_1 - 1/2}, W = B - h is a Brownian motion (Theorem 5.1 with \theta \equiv 1).

    Writing \mathbb P in terms of \mathbb Q and then renaming W as B gives \mathbb P\Big(\sup_t |B_t - t| < \varepsilon\Big) = \mathbb E\Big[ \mathbf 1_{\{\sup_t |B_t| < \varepsilon\}}\, e^{-B_1 - 1/2} \Big]. On the event, |B_1| < \varepsilon, so the ratio lies between e^{-1/2 - \varepsilon} and e^{-1/2 + \varepsilon}. The factor e^{-\frac12 \int h'^2} is the energy cost of following the path h.

11. ★★ Tanaka’s equation in detail. Let X be a Brownian motion and B = \int \mathrm{sgn}(X)\, dX. (a) Prove that \int_0^t \mathbf 1_{\{X_s = 0\}} ds = 0 almost surely, and deduce that \int_0^t \mathbf 1_{\{X_s = 0\}} dX_s = 0.

    (b) Conclude that X and -X both solve dY = \mathrm{sgn}(Y)\, dB.

    Solution. (a) By Fubini, \mathbb E \int_0^t \mathbf 1_{\{X_s = 0\}} ds = \int_0^t \mathbb P(X_s = 0)\, ds = 0. The stochastic integral has bracket \int \mathbf 1_{\{X = 0\}} ds = 0, so it vanishes (Proposition 2.13).

    (b) With the convention \mathrm{sgn}(0) = -1, one has \mathrm{sgn}(-x) = -\mathrm{sgn}(x) - 2 \cdot \mathbf 1_{\{x = 0\}}. Hence \int \mathrm{sgn}(-X)\, dB = -\int \mathrm{sgn}(X)\, dB - 2 \int \mathbf 1_{\{X = 0\}}\, dB = -X - 0.

12. ★★ Tweedie’s formula for Gaussian data. Let \mu = \mathcal N(m, s^2 I_d) and v_t = \alpha_t^2 s^2 + \sigma_t^2.

    (a) Compute p_t, \nabla \log p_t and \mathbb E[X_0 \mid X_t], and check Proposition 8.3. (b) Compute C_t. (c) What does denoising score matching return over affine functions s(x) = Ax + c?

    Solution. (a) p_t = \mathcal N(\alpha_t m, v_t I_d), \nabla \log p_t(x) = -(x - \alpha_t m)/v_t, and \mathbb E[X_0 \mid X_t = x] = m + \frac{\alpha_t s^2}{v_t} (x - \alpha_t m). Substituting into Tweedie’s formula, and using 1 - \alpha_t^2 s^2 / v_t = \sigma_t^2 / v_t, returns the score.

    (b) By the Pythagorean identity, C_t = \mathbb E|Z|^2/\sigma_t^2 - \mathbb E|\nabla \log p_t(X_t)|^2 = d/\sigma_t^2 - d/v_t.

    (c) The affine class contains the exact score, so the minimiser is A = -I/v_t, c = \alpha_t m / v_t.

13. ★★ Fixation in the Wright–Fisher diffusion. (a) Check that u(x) = -2(x \log x + (1-x) \log(1-x)) solves \frac12 x(1-x) u'' = -1 with u(0) = u(1) = 0. (b) Compute u(0.1) and interpret it. (c) Prove \mathbb P(X_\tau = 1) = x with Dynkin’s formula.

    Solution. (a) u'(x) = -2 \log\frac{x}{1-x} and u''(x) = -\frac{2}{x(1-x)}. (b) u(0.1) = 0.650: an allele at frequency 10 % disappears or fixes after 0.65\, N generations on average, and it is lost with probability 0.9.

    (c) For f(x) = x, \mathcal L f = 0, so X is a bounded martingale, and optional stopping gives \mathbb E X_\tau = x.

14. ★★★ ULA beyond Lipschitz gradients (Python). Let U(x) = x^2/2 + x^4/4 on \mathbb R.

    (a) Explain why Theorem 7.5 does not apply and what Proposition 7.9 predicts. (b) Run ULA from 0 with \gamma \in \{0.01, 0.02, 0.05, 0.1\} on 10^4 chains for 10^4 steps. Record the fraction of chains that diverge and the bias of \mathbb E X^2, using quadrature for the exact value. (c) Estimate the order of the bias in \gamma. (d) Repeat with tamed ULA, and discuss.

    Hints and expected result. (a) \nabla U is not globally Lipschitz; from |x| \ge 2/\sqrt\gamma the chain diverges with positive probability. (b) From 0, no chain diverges in practice: the noise almost never reaches 2/\sqrt\gamma \ge 6.3.

    The exact value is \mathbb E_\pi X^2 = 0.4679, and the ULA biases are about +0.003, +0.005, +0.014 and +0.031. Long runs are needed to see them: compute the stationary law of each chain on a grid, or average over chains and iterations.

    (c) The bias is linear in \gamma, about 0.3\gamma, as the weak order 1 of Euler–Maruyama predicts. (d) Taming gives biases of about +0.011, +0.023, +0.061 and +0.134: stability everywhere costs a fourfold larger constant. Open question: design a scheme that is both stable and as accurate as ULA where ULA is stable.

15. ★★★ Diffusion model on Gaussian data: exact analysis. Let \mu = \mathcal N(0, s^2) in dimension one, with the exact score, the reverse SDE started from \mathcal N(0,1) in continuous time, and t_0 = 0. (a) Show that the output is Gaussian with variance s^2 + \frac{(1 - s^2)\, s^4}{(e^{2T} + s^2 - 1)^2}. (b) For s^2 = 4 and T \in \{1, 2, 3\}, compare \mathrm{KL}(\mu \| \hat q), \mathrm{KL}(p_T \| \gamma_1) and e^{-2T} \mathrm{KL}(\mu \| \gamma_1). (c) Which step of the proof of Theorem 8.4 is loose here?

    Solution. (a) With v_u = 1 + (s^2 - 1)e^{-2u}, the score is -y/v_{T-t}, and the reverse variance solves V' = 2(1 - 2/v_{T-t})V + 2. The deviation from the exact solution is multiplied by \exp\big(2\int_0^T (1 - 2/v_u)\, du\big) = e^{2T} s^4 / (e^{2T} + s^2 - 1)^2, and it starts at 1 - v_T = (1 - s^2) e^{-2T}.

    (b) All laws are centred Gaussians, so each divergence is explicit. For T = 1, 2, 3:

    Quantity	T = 1	T = 2	T = 3
    \mathrm{KL}(\mu \Vert \hat q)	3.6 \times 10^{-3}	3.3 \times 10^{-6}	1.3 \times 10^{-9}
    \mathrm{KL}(p_T \Vert \gamma_1)	3.3 \times 10^{-2}	7.3 \times 10^{-4}	1.4 \times 10^{-5}
    e^{-2T}\, \mathrm{KL}(\mu \Vert \gamma_1)	0.109	0.015	0.0020

    (c) Two steps are loose. The log-Sobolev rate e^{-2T} is attained only by mean shifts, so a variance error decays like e^{-4T}. Data processing ignores that the reverse dynamics itself contracts the initial error, which gives the observed e^{-8T} decay in KL.

16. ★★★ Memorisation (Python, open-ended). Replace \mu in Example 8.5 by the empirical measure of n = 20 samples of it, and use the exact score of the resulting Gaussian mixture. (a) Show that, as t_0 \to 0, the outputs concentrate on the training points: measure the distance to the nearest training point. (b) Propose a remedy, such as a larger t_0, kernel smoothing of the empirical measure, or a smooth parametric score. Study the trade-off between fidelity (a Kolmogorov–Smirnov distance to \mu) and novelty (distance to the training set).

    Hints and expected result. With the exact empirical score, the output at time t_0 is a training point plus noise of standard deviation \sigma_{t_0} \approx \sqrt{2 t_0}, so the median distance to the training set shrinks like \sqrt{t_0}. Larger t_0 trades memorisation for blur; this trade-off is an open question for trained networks (ch. 9.3).

17. ★★★ Kinetic Langevin dynamics. Consider dX = V\, dt, dV = -\nabla U(X)\, dt - \lambda V\, dt + \sqrt{2\lambda}\, dB. (a) Write its generator and show that \rho(x, v) \propto e^{-U(x) - |v|^2/2} is invariant. (b) Explain why the proof of Theorem 7.2 fails. (c) Implement a splitting scheme and compare it with ULA on \mathcal N(0, \mathrm{diag}(1, 10^{-2})) at an equal number of gradient evaluations (Python, open-ended).

    Solution of (a)–(b). The generator is \mathcal L = v \cdot \nabla_x - \nabla U \cdot \nabla_v - \lambda v \cdot \nabla_v + \lambda \Delta_v. The Hamiltonian part gives -\nabla_x \cdot (v \rho) + \nabla_v \cdot (\nabla U \rho) = v \cdot \nabla U\, \rho - \nabla U \cdot v\, \rho = 0, and the friction–noise part gives \lambda \nabla_v \cdot (v\rho + \nabla_v \rho) = 0. The entropy now decreases at rate \lambda \int |\nabla_v \log(\mu_t/\rho)|^2 d\mu_t, which vanishes for laws of the form f(x)\, e^{-|v|^2/2} that are not at equilibrium. So no log-Sobolev argument applies directly, and convergence needs hypocoercivity (Villani, 2009).

Projects and assessment

The official assessment, on master-mva.com, is a written exam; no weighting, homework or project is displayed, and the Université Paris-Saclay course detail gives the assessment at the level of the constituent element only. The projects below are therefore proposals for practice and for self-assessment, sized like a typical MVA paper study. The written-exam proposal mirrors the official format.

Solo projects

1. [Solo · ≈ 20 h] Step size and bias of ULA for Bayesian logistic regression. How does the bias of ULA scale with the step, and when does MALA pay off? A paper study of Durmus and Moulines (2019).
   • Implement: ULA for several steps \gamma \le 1/L, MALA, and a long MALA run as reference; compute m and L of the posterior.
   • Data: the breast-cancer data set bundled with scikit-learn (569 samples, 30 standardised features), Gaussian prior of precision 1.
   • Deliverables: code, 4–6 page report, 10-minute talk. Chapters: ch. 4, 7.
   • Success criterion: the bias of posterior means and variances, against \gamma, has a log–log slope between 0.8 and 1.2; the report compares it with Theorem 7.5 and finds the cheapest method for a fixed accuracy.
2. [Solo · ≈ 20 h] Strong and weak orders of SDE schemes. Do Euler–Maruyama and Milstein reach their theoretical orders, and what happens for non-Lipschitz coefficients? Based on Kloeden and Platen (1992) and Talay and Tubaro (1990).
   • Implement: both schemes for geometric Brownian motion (exact solution known) and for the Cox–Ingersoll–Ross process, with a positivity-preserving truncation for the latter.
   • Deliverables: code, 4–6 page report with error-versus-step plots, 10-minute talk. Chapters: ch. 4, 6.
   • Success criterion: fitted strong orders 0.5 \pm 0.1 and 1.0 \pm 0.1 on geometric Brownian motion, weak order 1.0 \pm 0.15, and a documented failure mode for CIR.
3. [Solo · ≈ 20 h] Girsanov importance sampling for a rare exit. Estimate the probability that the Ornstein–Uhlenbeck process dX = -X\, dt + \sqrt 2\, dB, started at 0, exceeds a level a before time T = 1, with a chosen so that this probability is about 10^{-6}.
   • Implement: crude Monte Carlo, constant-drift importance sampling, and a state-dependent drift change guided by large deviations, with the Brownian-bridge correction of Example 5.6.
   • Reference value: finite differences for the Kolmogorov backward equation of the hitting probability.
   • Deliverables: code, 4–6 page report, 10-minute talk. Chapters: ch. 1, 4, 5.
   • Success criterion: relative error below 5 % (95 % level) with 10^4 paths, matching the PDE value.

Group projects

1. [Group of 2–3 · ≈ 30 h each] Score-based generative modelling through SDEs. A paper study of Song, Sohl-Dickstein, Kingma, Kumar, Ermon and Poole (2021).
   • Implement: denoising score matching for the variance-preserving SDE, the reverse-SDE sampler, the probability-flow ODE, and a predictor–corrector sampler with a Langevin corrector.
   • Data: two-moons and Swiss-roll samples from scikit-learn, and the 8×8 digits data set bundled with scikit-learn (1 797 images).
   • Split: member A trains the score networks; member B writes the samplers; member C writes the evaluation (sliced Wasserstein distance, likelihood through the ODE).
   • Deliverables: code, 8-page joint report, 20-minute talk. Chapters: ch. 4, 5, 8.
   • Success criterion: sliced Wasserstein distance to held-out data within twice the distance between two halves of the held-out data.
2. [Group of 2–3 · ≈ 30 h each] Checking the KL guarantees of diffusion models. A paper study of Conforti, Durmus and Gentiloni Silveri (2025), with Chen et al. (2023).
   • Setting: Gaussian mixtures in dimensions 1, 2 and 10, where the scores are exact; controlled score errors are added.
   • Implement: Euler and exponential-integrator samplers; KL estimates by quadrature in low dimension and by Monte Carlo with exact densities otherwise.
   • Split: member A derives each term of Theorem 8.4 for mixtures; member B runs the samplers; member C estimates the divergences and their uncertainty.
   • Deliverables: code, 8-page joint report, 20-minute talk. Chapters: ch. 5, 7, 8.
   • Success criterion: the initialisation, score and discretisation terms reproduced with the predicted dependence on T, on the score error and on the step, over at least one decade each.
3. [Group of 2–3 · ≈ 30 h each] Kinetic against overdamped Langevin on ill-conditioned targets. A paper study of Cheng, Chatterji, Bartlett and Jordan (2018).
   • Implement: ULA, MALA and a kinetic Langevin discretisation.
   • Targets: Gaussians with condition number 100 in dimensions 10 to 1 000 (exact W_2), and the Bayesian logistic regression of solo project 1.
   • Split: member A writes the samplers; member B the Gaussian experiments; member C the logistic experiments and effective sample sizes.
   • Deliverables: code, 8-page joint report, 20-minute talk. Chapters: ch. 7, 9.
   • Success criterion: for each method, the exponents of d and of 1/\varepsilon in the number of gradient evaluations needed for W_2 \le \varepsilon, fitted over at least one decade of each.
   • On the Gaussian targets, the exact ULA bias of Example 7.4 predicts \sqrt d / \varepsilon evaluations up to logarithms, and the fitted exponents must match it within 0.15.
   • The report compares every fit with the worst-case bounds d / \varepsilon^2 for ULA (Corollary 7.6) and \sqrt d / \varepsilon for kinetic Langevin (Cheng et al.), and explains the gaps.

Written exam (proposal)

Three hours, no documents, four independent questions.

1. Question (5 points). State and prove the convergence in L^2 of the quadratic variation of Brownian motion along partitions of [0, t]. Deduce that Brownian paths have infinite variation. State the Itô isometry for simple integrands and prove it.

   Expected: the variance computation \mathbb E(Q_n - t)^2 = 2 \sum \delta_i^2 \le 2t |\pi_n| (Theorem 1.8); the bound \sum \Delta_i^2 \le \max_i |\Delta_i| \sum |\Delta_i| (Corollary 1.9); vanishing cross terms by conditioning and diagonal terms \mathbb E[\xi_i^2] \delta_i (Proposition 3.2).

2. Question (5 points). Let dX = -\theta X\, dt + \sigma\, dB. Solve the equation, give the law of X_t and the stationary law, and check the stationary law with the Fokker–Planck equation. For \sigma = \sqrt 2, write ULA with step \gamma and compute its stationary variance.

   Expected: X_t = e^{-\theta t} X_0 + \sigma \int_0^t e^{-\theta(t-s)} dB_s, Gaussian with variance \sigma^2(1 - e^{-2\theta t})/(2\theta); the zero-flux equation (Example 4.6, Theorem 4.12); ULA X_{k+1} = (1 - \gamma\theta) X_k + \sqrt{2\gamma} Z has stationary variance 2/(\theta(2 - \gamma\theta)), a bias of order \gamma.

3. Question (5 points). State Girsanov’s theorem. Prove that the KL divergence between the path laws of dX = b_i(X)\, dt + dB, i = 1, 2, with the same initial law, equals \frac12 \mathbb E_1 \int_0^T |b_1 - b_2|^2(X_t)\, dt. Apply it to two Ornstein–Uhlenbeck drifts.

   Expected: Theorem 5.1 with Novikov’s condition; the substitution d\omega = b_1 dt + dW^1 and the identity (b_1 - b_2) \cdot b_1 - \frac12 |b_1|^2 + \frac12 |b_2|^2 = \frac12 |b_1 - b_2|^2 (Corollary 5.4); the value (\theta_1 - \theta_2)^2 T / (4\theta_1) (Exercise 8).

4. Question (5 points). For the forward process dX = -X\, dt + \sqrt 2\, dB, show that the marginals of the SDE with drift y + 2\nabla \log p_{T-t}(y) are p_{T-t}. Then prove that denoising score matching and score matching have the same minimisers.

   Expected: the Fokker–Planck computation of Theorem 8.1, using p \nabla \log p = \nabla p; Tweedie’s formula and the Pythagorean identity of Proposition 8.3.

Total: 5 + 5 + 5 + 5 = 20 points.

Grading rubric (proposal)

References

Official references (master-mva.com, 2026–2027)

• D. Revuz, M. Yor, Continuous Martingales and Brownian Motion, 3rd ed., Grundlehren der Mathematischen Wissenschaften 293, Springer-Verlag, Berlin, 1999 — the reference for ch. 1–6.
• G. Conforti, A. Durmus, M. Gentiloni Silveri, “KL convergence guarantees for score diffusion models under minimal data assumptions”, SIAM Journal on Mathematics of Data Science 7(1), 86–109, 2025 — ch. 8.4.
• A. Durmus, É. Moulines, “High-dimensional Bayesian inference via the unadjusted Langevin algorithm”, Bernoulli 25(4A), 2854–2882, 2019 — ch. 7.

Textbooks

• I. Karatzas, S. E. Shreve, Brownian Motion and Stochastic Calculus, 2nd ed., Graduate Texts in Mathematics 113, Springer, 1991.
• J.-F. Le Gall, Brownian Motion, Martingales, and Stochastic Calculus, Graduate Texts in Mathematics 274, Springer, 2016.
• B. Øksendal, Stochastic Differential Equations: An Introduction with Applications, 6th ed., Universitext, Springer, 2003.
• D. W. Stroock, S. R. S. Varadhan, Multidimensional Diffusion Processes, Grundlehren der Mathematischen Wissenschaften 233, Springer, 1979.
• P. E. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, Applications of Mathematics 23, Springer, 1992.
• G. A. Pavliotis, Stochastic Processes and Applications: Diffusion Processes, the Fokker–Planck and Langevin Equations, Texts in Applied Mathematics 60, Springer, 2014.
• D. Bakry, I. Gentil, M. Ledoux, Analysis and Geometry of Markov Diffusion Operators, Grundlehren der Mathematischen Wissenschaften 348, Springer, 2014.

Stochastic calculus and diffusions

• R. E. A. C. Paley, N. Wiener, A. Zygmund, “Notes on random functions”, Mathematische Zeitschrift 37, 647–668, 1933.
• K. Itô, “Stochastic integral”, Proceedings of the Imperial Academy 20(8), 519–524, 1944.
• I. V. Girsanov, “On transforming a certain class of stochastic processes by absolutely continuous substitution of measures”, Theory of Probability and Its Applications 5(3), 285–301, 1960.
• T. Yamada, S. Watanabe, “On the uniqueness of solutions of stochastic differential equations”, Journal of Mathematics of Kyoto University 11(1), 155–167, 1971.
• A. K. Zvonkin, “A transformation of the phase space of a diffusion process that removes the drift”, Mathematics of the USSR-Sbornik 22(1), 129–149, 1974.
• A. Yu. Veretennikov, “On strong solutions and explicit formulas for solutions of stochastic integral equations”, Mathematics of the USSR-Sbornik 39(3), 387–403, 1981.
• D. Talay, L. Tubaro, “Expansion of the global error for numerical schemes solving stochastic differential equations”, Stochastic Analysis and Applications 8(4), 483–509, 1990.
• M. Hutzenthaler, A. Jentzen, P. E. Kloeden, “Strong and weak divergence in finite time of Euler’s method for stochastic differential equations with non-globally Lipschitz continuous coefficients”, Proceedings of the Royal Society A 467(2130), 1563–1576, 2011.
• A. Einstein, “Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen”, Annalen der Physik 17, 549–560, 1905.
• P. Langevin, “Sur la théorie du mouvement brownien”, Comptes Rendus de l’Académie des Sciences 146, 530–533, 1908.
• R. Holley, D. Stroock, “Logarithmic Sobolev inequalities and stochastic Ising models”, Journal of Statistical Physics 46(5–6), 1159–1194, 1987.
• R. Jordan, D. Kinderlehrer, F. Otto, “The variational formulation of the Fokker–Planck equation”, SIAM Journal on Mathematical Analysis 29(1), 1–17, 1998.
• C. Villani, Hypocoercivity, Memoirs of the American Mathematical Society 202(950), 2009.

Langevin sampling

• G. O. Roberts, R. L. Tweedie, “Exponential convergence of Langevin distributions and their discrete approximations”, Bernoulli 2(4), 341–363, 1996.
• G. O. Roberts, A. Gelman, W. R. Gilks, “Weak convergence and optimal scaling of random walk Metropolis algorithms”, Annals of Applied Probability 7(1), 110–120, 1997.
• G. O. Roberts, J. S. Rosenthal, “Optimal scaling of discrete approximations to Langevin diffusions”, Journal of the Royal Statistical Society B 60(1), 255–268, 1998.
• M. Welling, Y. W. Teh, “Bayesian learning via stochastic gradient Langevin dynamics”, ICML, 2011.
• M. Pereyra, “Proximal Markov chain Monte Carlo algorithms”, Statistics and Computing 26(4), 745–760, 2016.
• A. S. Dalalyan, “Theoretical guarantees for approximate sampling from smooth and log-concave densities”, Journal of the Royal Statistical Society B 79(3), 651–676, 2017.
• A. Durmus, É. Moulines, “Nonasymptotic convergence analysis for the unadjusted Langevin algorithm”, Annals of Applied Probability 27(3), 1551–1587, 2017.
• A. Durmus, É. Moulines, M. Pereyra, “Efficient Bayesian computation by proximal Markov chain Monte Carlo: when Langevin meets Moreau”, SIAM Journal on Imaging Sciences 11(1), 473–506, 2018.
• X. Cheng, N. S. Chatterji, P. L. Bartlett, M. I. Jordan, “Underdamped Langevin MCMC: a non-asymptotic analysis”, COLT, 2018.
• A. Wibisono, “Sampling as optimization in the space of measures: the Langevin dynamics as a composite optimization problem”, COLT, 2018.
• A. S. Dalalyan, A. Karagulyan, “User-friendly guarantees for the Langevin Monte Carlo with inaccurate gradient”, Stochastic Processes and their Applications 129(12), 5278–5311, 2019.
• N. Brosse, A. Durmus, É. Moulines, S. Sabanis, “The tamed unadjusted Langevin algorithm”, Stochastic Processes and their Applications 129(10), 2019.
• A. Durmus, S. Majewski, B. Miasojedow, “Analysis of Langevin Monte Carlo via convex optimization”, Journal of Machine Learning Research 20(73), 1–46, 2019.
• A. Eberle, A. Guillin, R. Zimmer, “Couplings and quantitative contraction rates for Langevin dynamics”, Annals of Probability 47(4), 1982–2010, 2019.
• S. Vempala, A. Wibisono, “Rapid convergence of the unadjusted Langevin algorithm: isoperimetry suffices”, NeurIPS, 2019.
• R. Shen, Y. T. Lee, “The randomized midpoint method for log-concave sampling”, NeurIPS, 2019.

Diffusion models

• B. D. O. Anderson, “Reverse-time diffusion equation models”, Stochastic Processes and their Applications 12(3), 313–326, 1982.
• U. G. Haussmann, É. Pardoux, “Time reversal of diffusions”, Annals of Probability 14(4), 1188–1205, 1986.
• A. Hyvärinen, “Estimation of non-normalized statistical models by score matching”, Journal of Machine Learning Research 6, 695–709, 2005.
• P. Vincent, “A connection between score matching and denoising autoencoders”, Neural Computation 23(7), 1661–1674, 2011.
• C. Léonard, “A survey of the Schrödinger problem and some of its connections with optimal transport”, Discrete and Continuous Dynamical Systems A 34(4), 1533–1574, 2014.
• J. Sohl-Dickstein, E. Weiss, N. Maheswaranathan, S. Ganguli, “Deep unsupervised learning using nonequilibrium thermodynamics”, ICML, 2015.
• Y. Song, S. Ermon, “Generative modeling by estimating gradients of the data distribution”, NeurIPS, 2019.
• J. Ho, A. Jain, P. Abbeel, “Denoising diffusion probabilistic models”, NeurIPS, 2020.
• Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, B. Poole, “Score-based generative modeling through stochastic differential equations”, ICLR, 2021.
• J. Song, C. Meng, S. Ermon, “Denoising diffusion implicit models”, ICLR, 2021.
• V. De Bortoli, J. Thornton, J. Heng, A. Doucet, “Diffusion Schrödinger bridge with applications to score-based generative modeling”, NeurIPS, 2021.
• V. De Bortoli, “Convergence of denoising diffusion models under the manifold hypothesis”, Transactions on Machine Learning Research, 2022.
• H. Lee, J. Lu, Y. Tan, “Convergence for score-based generative modeling with polynomial complexity”, NeurIPS, 2022.
• T. Karras, M. Aittala, T. Aila, S. Laine, “Elucidating the design space of diffusion-based generative models”, NeurIPS, 2022.
• Q. Zhang, Y. Chen, “Path integral sampler: a stochastic control approach for sampling”, ICLR, 2022.
• S. Chen, S. Chewi, J. Li, Y. Li, A. Salim, A. R. Zhang, “Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions”, ICLR, 2023.
• H. Chen, H. Lee, J. Lu, “Improved analysis of score-based generative modeling: user-friendly bounds under minimal smoothness assumptions”, ICML, 2023.
• K. Oko, S. Akiyama, T. Suzuki, “Diffusion models are minimax optimal distribution estimators”, ICML, 2023.
• F. Vargas, W. Grathwohl, A. Doucet, “Denoising diffusion samplers”, ICLR, 2023.
• H. Chung, J. Kim, M. T. McCann, M. L. Klasky, J. C. Ye, “Diffusion posterior sampling for general noisy inverse problems”, ICLR, 2023.
• Y. Lipman, R. T. Q. Chen, H. Ben-Hamu, M. Nickel, M. Le, “Flow matching for generative modeling”, ICLR, 2023.
• M. S. Albergo, E. Vanden-Eijnden, “Building normalizing flows with stochastic interpolants”, ICLR, 2023.
• X. Liu, C. Gong, Q. Liu, “Flow straight and fast: learning to generate and transfer data with rectified flow”, ICLR, 2023.
• J. Benton, V. De Bortoli, A. Doucet, G. Deligiannidis, “Nearly d-linear convergence bounds for diffusion models via stochastic localization”, ICLR, 2024.