Optimal Transport for Machine Learning

Contents

• Introduction and overview
• 0. Review: measures, convexity and linear programming
• 1. The Monge and Kantorovich problems
• 2. Kantorovich duality
• 3. Brenier’s theorem and Gaussian transport
• 4. Discrete optimal transport: linear programming and entropic regularization
• 5. Semi-discrete optimal transport
• 6. Dynamic formulation and Wasserstein gradient flows
• 7. Optimal transport in machine learning
• 8. Complements for self-study
• 9. Research directions and open problems
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: a student with a good M1 in mathematics, computer science or physics who wants a rigorous way to compare probability distributions, and to understand why that question sits at the centre of generative modelling.

Empirical first: two parallel segments in the plane, one shifted by \theta, are at Kullback–Leibler divergence +\infty and Jensen–Shannon divergence \log 2 for every \theta \neq 0, so a generator trained with these divergences gets no gradient, while their Wasserstein distance is |\theta| (Example 2.6). The Fréchet inception distance that scores image generators is itself a Wasserstein distance between two Gaussians (ch. 3).

This course builds optimal transport (OT) from its two classical formulations, Monge’s maps and Kantorovich’s couplings, to its duality theory, Brenier’s theorem and the dynamic view of Benamou and Brenier. It then develops the three numerical approaches announced by the course description: linear programming, semi-discrete solvers and entropic regularization.

The last in-class chapter uses Wasserstein gradient flows and duality to read the training of wide neural networks, the token dynamics of transformers, GANs and diffusion models.

It shares tools with Convex Optimization (duality), Stochastic Calculus in Machine Learning (Langevin dynamics, diffusion models), Computational Statistics and Kernel Methods (the maximum mean discrepancy).

What you will be able to do

• Formulate the Monge and Kantorovich problems, prove that optimal couplings exist, and compute Wasserstein distances on the line by sorting.
• Derive the dual of a transport problem, use c-transforms, and apply the Kantorovich–Rubinstein formula for W_1.
• State and use Brenier’s theorem, and recognize an optimal map as the gradient of a convex function.
• Compute the optimal map and the Wasserstein distance between two Gaussians with the Bures metric.
• Solve discrete transport problems exactly by linear programming, and approximately with Sinkhorn’s algorithm in PyTorch, with error bounds.
• Solve semi-discrete problems by gradient ascent on the concave dual, deterministic or stochastic.
• Explain the Benamou–Brenier formula and write the Wasserstein gradient flow of a functional, including the Fokker–Planck equation.
• Read the training of a two-layer network, the dynamics of self-attention, GANs and diffusion models as transport problems.
• Use entropic transport and Sinkhorn divergences as differentiable losses, and know their statistical behaviour in high dimension.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com presents optimal transport as a foundational mathematical theory that connects optimization, partial differential equations and probability. It offers a powerful framework for comparing probability distributions, and the description notes that it has become an important tool in machine learning, especially for designing and evaluating generative models.

The announced mathematical content is the Monge and Kantorovich formulations, Brenier’s theorem, the dual and dynamic formulations, the Bures metric on Gaussian distributions, and gradient flows. The announced numerical methods are linear programming, semi-discrete solvers and entropic regularization.

The announced applications are the training of neural networks via gradient flows, token dynamics in transformers, and the structure of GANs and diffusion models. The description states that the course focuses primarily on mathematical content rather than deep learning techniques.

It includes hands-on numerical programming sessions using PyTorch, along with exercise sessions, and the written exam reflects both components. The organisation announced is 8 courses of 3 h plus the exam. The Université Paris-Saclay course detail gives only the assessment type, with no programme of its own.

Every official topic has a section in these notes. The extensions beyond the syllabus are the review of ch. 0, flow matching with OT couplings in ch. 7.4, and two self-study chapters: the complements of ch. 8 (convergence of Sinkhorn, Sinkhorn divergences, statistical OT and sliced distances, barycenters, WGAN practice, numerical gradient flows) and the research directions of ch. 9.

Where this course fits in the programme

Official prerequisites: not displayed. In practice the course uses measure-theoretic probability, convex analysis, linear algebra, a little differential calculus and Python; ch. 0 reviews the essentials.

• Refreshers — Probability and Statistics (measures, laws, convergence), Linear Algebra, Analysis and Optimization (convexity, duality, matrix square roots) and Scientific Python, NumPy and PyTorch (the code of every chapter).
• Related courses of semester 1 — Convex Optimization (linear programming and Lagrange duality, ch. 2 and 4), Stochastic Calculus in Machine Learning (Langevin diffusions and diffusion models, ch. 6–7), Computational Statistics (Langevin sampling), Probabilistic Graphical Models and Deep Generative Models (GANs, VAEs), Theoretical Foundations of Deep Learning (training dynamics, ch. 7.1), Foundations of Large Language Models (transformers, ch. 7.2) and Geometric Data Analysis.
• Semester 2 courses that build on it — Generative Models for Images and Diffusion and Multiscale Models in AI (GANs, diffusion and flow models), Kernel Methods (the MMD at the other end of the Sinkhorn divergence), Geometry and Shape Spaces and Deformable Models and Geodesic Methods (geodesics and shape matching), Point Clouds and 3D Modelling (matching point clouds), and Deep Learning for Image Restoration and Synthesis.

Course outline

• 0. Review — push-forwards, weak convergence, convex conjugates, linear programming duality.
• 1. The Monge and Kantorovich problems — maps and couplings, existence, Wasserstein distances, the one-dimensional case.
• 2. Kantorovich duality — the discrete and continuous duals, c-transforms, the Kantorovich–Rubinstein formula.
• 3. Brenier’s theorem and Gaussian transport — optimal maps as gradients of convex functions; the Bures metric.
• 4. Discrete optimal transport — exact linear programming, entropic regularization, Sinkhorn’s algorithm in PyTorch.
• 5. Semi-discrete optimal transport — Laguerre cells, the concave dual, deterministic and stochastic solvers.
• 6. Dynamic formulation and Wasserstein gradient flows — displacement interpolation, Benamou–Brenier, JKO, Fokker–Planck.
• 7. Optimal transport in machine learning — mean-field training, token dynamics, GANs, diffusion and flow matching.
• 8. Complements for self-study — the material of the two announced sessions beyond the 18 h: proofs and practical tools.
• 9. Research directions and open problems — a self-study tour of current work.

Session plan

18 h in class = 6 sessions of 3 h. The volume of 18 h of lectures comes from the Université Paris-Saclay programme table, while master-mva.com announces 8 courses of 3 h plus the exam. The plan follows the official 18 h; the split into six sessions is a proposal, and the material of the two other announced sessions is gathered in the self-study chapter 8.

Level 1 = foundations, 2 = core results and algorithms, 3 = advanced and research-level material. The level never goes down from one session to the next.

Self-study, outside the 18 h (about 8 h in total): ch. 8 Complements: convergence of Sinkhorn, Sinkhorn divergences, statistical OT and sliced distances, barycenters, WGAN practice, numerical gradient flows (≈ 6 h, with Ex. 16) · ch. 9 Research directions and open problems (≈ 2 h).

How to study this module

• Read ch. 0 before session 1 if push-forwards or convex conjugates feel distant; ch. 1–3 use them without warning.
• After each session, do the exercises of its row, then run the code of the chapter and change one parameter: \varepsilon, the number of points, the step size, the temperature \beta.
• Implement Sinkhorn’s algorithm and the semi-discrete dual ascent yourself once before using a library; each is a dozen lines.
• The code uses NumPy, SciPy and PyTorch on the CPU, with synthetic data; every block runs in under half a minute on a laptop.
• The written exam allows no documents: learn the statements of Theorems 2.4, 3.2, 3.3 and 6.3, and the proof of Proposition 4.5, well enough to reproduce them.

0. Review: measures, convexity and linear programming

Optimal transport moves probability measures, so the course computes with push-forwards, convex conjugates and linear programs from the first session. This chapter collects the M1-level facts used later.

The Probability and Statistics and Linear Algebra, Analysis and Optimization refreshers cover the same ground at a slower pace.

0.1 Measures and push-forwards

Two cases are used constantly.

• For a discrete measure, the push-forward moves each atom: T_\sharp \Big( \sum_i a_i \delta_{x_i} \Big) = \sum_i a_i \delta_{T(x_i)} .
• For densities \rho_0, \rho_1 and a diffeomorphism T with T_\sharp \rho_0 = \rho_1, the change of variables formula gives the Jacobian equation \rho_0(x) = \rho_1\big(T(x)\big)\, |\det DT(x)| .

0.2 Convex conjugates

Two more facts are quoted without proof. A proper, convex, lower semicontinuous function satisfies \varphi^{**} = \varphi. A convex function that is finite on an open set is differentiable Lebesgue-almost everywhere there (Rockafellar, 1970, Theorem 25.5), and at such points \partial\varphi(x) = \{\nabla\varphi(x)\}.

0.3 Linear programming duality

1. The Monge and Kantorovich problems

How much does it cost to move a pile of sand into a hole of the same volume? This chapter states the question in Monge’s form, with a map, and in Kantorovich’s form, with a coupling. It proves that the second always has a solution and solves it explicitly on the real line.

1.1 Monge’s problem

Throughout, \mu and \nu are probability measures on spaces \mathcal{X} and \mathcal{Y}.

The cost of moving a unit of mass from x to y is c(x, y) \geq 0.

The running example is the cost \|x - y\|^p on \mathbb{R}^d, where \|\cdot\| is the Euclidean norm.

1.2 Kantorovich’s relaxation

Kantorovich (1942) allowed the mass at x to be split between several destinations. A transport plan is then a joint distribution of departures and arrivals.

The relaxation is consistent with Monge’s problem. The product \mu \otimes \nu is always a coupling, so \Pi(\mu,\nu) is never empty. A map T with T_\sharp\mu = \nu defines the coupling (\mathrm{Id}, T)_\sharp \mu with the same cost. Hence the Kantorovich value is at most the Monge value.

1.3 Wasserstein distances

Unlike the Kullback–Leibler divergence, W_p sees the geometry of the ground space. Two Dirac masses are at distance W_p(\delta_x, \delta_y) = \|x - y\| , which varies continuously with the points, while their Kullback–Leibler divergence is +\infty as soon as x \neq y. Ch. 2 turns this remark into the argument behind Wasserstein GANs.

1.4 Dimension one: sorting

On the real line, optimal transport has a closed form. It uses the cumulative distribution function F_\mu(t) = \mu((-\infty, t]) and the quantile function F_\mu^{-1}(s) = \inf\{t : F_\mu(t) \geq s\}.

The code checks Example 1.9, then solves random instances in three ways: sorting (Theorem 1.8), the assignment solver of SciPy (Proposition 1.5) and a generic linear-programming solver applied to Definition 1.3.

import numpy as np
from scipy.optimize import linear_sum_assignment, linprog

def ot_lp(a, b, C):
    """Kantorovich problem as a linear program (Definition 1.3)."""
    n, m = C.shape
    A_eq = np.vstack([np.kron(np.eye(n), np.ones(m)),     # row sums = a
                      np.kron(np.ones(n), np.eye(m))])    # column sums = b
    res = linprog(C.ravel(), A_eq=A_eq, b_eq=np.concatenate([a, b]),
                  bounds=(0, None), method="highs")
    return res.fun, res.x.reshape(n, m)

# Example 1.9: three points against three points, quadratic cost
x, y = np.array([0.0, 1.0, 3.0]), np.array([6.0, 2.0, 5.0])
print("sorted pairing:", np.mean((np.sort(x) - np.sort(y)) ** 2))   # 29/3

rng = np.random.default_rng(0)
n = 40
for d in (1, 2):
    X, Y = rng.normal(size=(n, d)), rng.normal(size=(n, d)) + 1.0
    C = ((X[:, None, :] - Y[None, :, :]) ** 2).sum(-1)
    r, c = linear_sum_assignment(C)                      # Monge: a permutation
    w_lp, P = ot_lp(np.full(n, 1 / n), np.full(n, 1 / n), C)
    line = f"d={d}: assignment {C[r, c].mean():.6f}  LP {w_lp:.6f}"
    if d == 1:
        line += f"  sorting {np.mean((np.sort(X[:, 0]) - np.sort(Y[:, 0])) ** 2):.6f}"
    print(line, " nonzeros in P:", int((P > 1e-12).sum()))

The three values agree (1.980791 in dimension 1, 2.369378 in dimension 2), and the linear program returns a plan with exactly 40 nonzero entries: a permutation, as Proposition 1.5 predicts.

2. Kantorovich duality

The Kantorovich problem is a linear program, so it has a dual. The dual variables are prices: they certify optimality, and for the distance cost they turn W_1 into a supremum over Lipschitz functions, the formula behind Wasserstein GANs.

2.1 The discrete dual

The dual has an economic reading, the shipper’s problem (Villani, 2009, ch. 5). A shipping company offers to do the transport. It charges f_i to pick up a unit of mass at x_i and g_j to deliver it at y_j.

The constraint f_i + g_j \leq C_{ij} keeps the offer competitive with doing the transport yourself. Duality says that the company, maximizing its revenue, earns exactly the optimal transport cost.

2.2 c-transforms and the general dual

2.3 The case of W_1: Kantorovich–Rubinstein

3. Brenier’s theorem and Gaussian transport

For the quadratic cost, the dual of ch. 2 becomes a problem over convex functions. Fenchel–Young equality then forces the optimal plan to be a map, the gradient of a convex function. Between Gaussian measures this map is affine, and the distance has a closed form built on the Bures metric.

3.1 Brenier’s theorem

Three remarks place the theorem.

• The density assumption cannot be dropped: Example 1.2 has no map at all.
• On the line, gradients of convex functions are the non-decreasing maps, and Theorem 1.8 is recovered.
• With densities \rho_0 and \rho_1 and a potential \varphi of class C^2, the Jacobian equation of ch. 0.1 for T = \nabla\varphi becomes the Monge–Ampère equation \rho_0(x) = \rho_1\big(\nabla\varphi(x)\big) \det D^2\varphi(x) .

3.2 Gaussian measures and the Bures metric

The Fréchet inception distance (FID) used to evaluate image generators (Heusel et al., 2017) is exactly this formula. It fits a Gaussian to the features of real images and another to the features of generated images, then reports W_2^2 between the two fits.

The code computes the closed form of Example 3.5, checks that A is symmetric with A \Sigma_0 A = \Sigma_1, and compares with the optimal assignment between samples.

import numpy as np
from scipy.linalg import sqrtm
from scipy.optimize import linear_sum_assignment

def gaussian_w2(m0, S0, m1, S1):
    """Squared W2 between Gaussians and the Brenier map matrix (Theorem 3.3)."""
    r0 = np.real(sqrtm(S0))
    M = np.real(sqrtm(r0 @ S1 @ r0))
    bures2 = np.trace(S0) + np.trace(S1) - 2 * np.trace(M)
    r0inv = np.linalg.inv(r0)
    A = r0inv @ M @ r0inv
    return np.sum((m0 - m1) ** 2) + bures2, A

m0, S0 = np.zeros(2), np.array([[1.0, 0.0], [0.0, 4.0]])
m1, S1 = np.array([1.0, 0.0]), np.array([[2.0, 1.0], [1.0, 2.0]])
w2, A = gaussian_w2(m0, S0, m1, S1)
print("closed form W2^2:", round(w2, 4))
print("A symmetric:", np.allclose(A, A.T), " A S0 A = S1:", np.allclose(A @ S0 @ A, S1))

rng = np.random.default_rng(0)
for n in (200, 800, 2000):                  # empirical OT between samples
    X = rng.multivariate_normal(m0, S0, n)
    Y = rng.multivariate_normal(m1, S1, n)
    C = ((X[:, None] - Y[None]) ** 2).sum(-1)
    r, c = linear_sum_assignment(C)
    print(f"n={n:5d}  empirical W2^2 = {C[r, c].mean():.4f}")

The closed form gives 1.7712. The empirical values (1.9386, 1.9208, 1.8143) approach it slowly from above. The plug-in estimate is biased upwards, because W_2^2 is jointly convex; ch. 8 shows how its convergence degrades with the dimension.

4. Discrete optimal transport: linear programming and entropic regularization

Numerical optimal transport starts with two point clouds. This chapter first solves the linear program exactly, then regularizes it with entropy. The regularized problem is solved by Sinkhorn’s algorithm: a few lines of PyTorch, parallel on a GPU and differentiable.

4.1 Exact solvers

Exact solvers exploit this structure. The network simplex moves from vertex to vertex of U(a, b); the Hungarian method (Kuhn, 1955) solves the assignment problem in polynomial time, O(n^3) in its later implementations; the auction algorithm (Bertsekas, 1988) is a parallel alternative. The code of ch. 1 used the HiGHS solver through scipy.optimize.linprog. All of them cost roughly cubic time in practice, which is prohibitive beyond a few thousand points.

4.2 Entropic regularization

The hard constraint f_i + g_j \leq C_{ij} of Theorem 2.1 has become a soft exponential penalty. Maximizing D_\varepsilon in f with g fixed has a closed form, and so has the maximization in g.

Each iteration costs O(nm) operations, all of them matrix–vector products or log-sum-exp reductions, which run in parallel on a GPU. The log domain avoids the underflow of e^{-C/\varepsilon} when \varepsilon is small.

Since \eta(K) grows like e^{2\max C / \varepsilon}, the rate degrades as \varepsilon \to 0: accuracy and speed trade off.

The PyTorch code compares Sinkhorn with the exact assignment for three values of \varepsilon, then differentiates the dual value through the iterations with autograd and checks Proposition 4.7.

import math, numpy as np, torch
from scipy.optimize import linear_sum_assignment
torch.set_num_threads(1)
torch.manual_seed(0)

def sinkhorn_log(a, b, C, eps, iters=2000):
    """Log-domain Sinkhorn (Algorithm 4.4): block ascent on the dual (f, g)."""
    f, g = torch.zeros_like(a), torch.zeros_like(b)
    for _ in range(iters):
        f = -eps * torch.logsumexp((g[None, :] - C) / eps + torch.log(b)[None, :], dim=1)
        g = -eps * torch.logsumexp((f[:, None] - C) / eps + torch.log(a)[:, None], dim=0)
    P = torch.exp((f[:, None] + g[None, :] - C) / eps) * a[:, None] * b[None, :]
    return P, f, g

n = 100
x = torch.randn(n, 2, dtype=torch.float64)
y = torch.randn(n, 2, dtype=torch.float64) * 0.5 + torch.tensor([2.0, 0.0], dtype=torch.float64)
a = b = torch.full((n,), 1.0 / n, dtype=torch.float64)
C = torch.cdist(x, y) ** 2
r, c = linear_sum_assignment(C.numpy())
exact = C.numpy()[r, c].mean()
for eps in (1.0, 0.1, 0.01):
    P, f, g = sinkhorn_log(a, b, C, eps)
    gap = (P * C).sum().item() - exact
    print(f"eps={eps:5.2f}  <C,P>-OT = {gap:.4f}   bound eps*log(n) = {eps*math.log(n):.4f}"
          f"   marginal error = {(P.sum(1)-a).abs().sum().item():.1e}")

# Gradient of the entropic cost with respect to the source positions (Proposition 4.7)
xs = x.clone().requires_grad_(True)
P, f, g = sinkhorn_log(a, b, torch.cdist(xs, y) ** 2, eps=0.1, iters=300)
dual_value = (f * a).sum() + (g * b).sum()        # = OT_eps at convergence
dual_value.backward()                              # autograd through 300 iterations
bary = (P @ y) / P.sum(1, keepdim=True)            # barycentric projection
print("grad = 2 a_i (x_i - T(x_i)):", torch.allclose(xs.grad, 2 * a[:, None] * (x - bary), atol=1e-6))

The gaps are 0.6729, 0.0670 and 0.0029 for \varepsilon = 1, 0.1, 0.01, well inside the bounds 4.61, 0.46 and 0.046 of Proposition 4.5. For \varepsilon = 0.01 the marginal error after 2000 iterations is still 9 \times 10^{-5}, the slowdown predicted by Theorem 4.6. The autograd gradient matches the barycentric formula.

5. Semi-discrete optimal transport

Between a density and a finite set of points, optimal transport partitions space into cells, one per point. The dual becomes a concave maximization over m numbers, solved by gradient ascent, Newton’s method or stochastic gradient. Applications include quantization, matching a continuous model to a finite data set, and the semi-discrete view of generative models.

Throughout, \mu has a density on a compact convex set \mathcal{X} \subset \mathbb{R}^d, \nu = \sum_{j=1}^m b_j \delta_{y_j} with distinct y_j, and c(x, y) = \|x - y\|^2.

Three algorithms follow.

• Gradient ascent: g \leftarrow g + \tau\,(b - \mu(L(g))). Each step computes a power diagram and the masses of its cells, exactly in 2D and 3D with computational geometry (Aurenhammer, Hoffmann and Aronov, 1998; Lévy, 2015).
• Damped Newton: the Hessian is minus a weighted graph Laplacian. Its entry (j, k), for j \neq k, is the integral of the density of \mu over the common facet of L_j and L_k, divided by 2\|y_j - y_k\|. Kitagawa, Mérigot and Thibert (2019) prove global convergence, quadratic near the solution.
• Stochastic gradient (Genevay, Cuturi, Peyré and Bach, 2016): sample x \sim \mu, find the index j^\star of its cell, and use the unbiased gradient b - e_{j^\star}. Only samples of \mu are needed, as with a generative model.

The code discretizes the uniform measure on [0, 1]^2 by a 200 \times 200 grid, runs gradient ascent on F for five weighted points, then Algorithm 5.4 with fresh uniform samples.

import numpy as np
rng = np.random.default_rng(1)

# mu = uniform density on [0,1]^2, discretised on a fine grid; nu = 5 weighted Diracs
s = (np.arange(200) + 0.5) / 200
X = np.stack(np.meshgrid(s, s), -1).reshape(-1, 2)          # 40 000 grid points
Y = rng.uniform(0.2, 0.8, size=(5, 2))
b = np.array([0.1, 0.15, 0.2, 0.25, 0.3])
C = ((X[:, None, :] - Y[None, :, :]) ** 2).sum(-1)          # c(x, y_j)

def dual_and_masses(g):
    """F(g) and the masses mu(L_j(g)) of the Laguerre cells (Proposition 5.2)."""
    vals = C - g[None, :]
    F = g @ b + vals.min(1).mean()
    return F, np.bincount(vals.argmin(1), minlength=len(b)) / len(X)

g = np.zeros(5)
for k in range(301):                                        # gradient ascent
    F, mass = dual_and_masses(g)
    if k in (0, 10, 50, 300):
        print(f"iter {k:3d}  F(g) = {F:.5f}  max |mu(L_j) - b_j| = {np.abs(mass - b).max():.4f}")
    g += 0.2 * (b - mass)

# Stochastic version (Algorithm 5.4): one fresh sample x ~ mu per step, averaged iterates
h, h_avg = np.zeros(5), np.zeros(5)
for k in range(1, 20001):
    x = rng.uniform(size=2)
    j = np.argmin(((x - Y) ** 2).sum(1) - h)
    grad = b.copy(); grad[j] -= 1.0                         # unbiased gradient of F
    h += 0.5 / np.sqrt(k) * grad
    h_avg += (h - h_avg) / k
F_sgd, mass_sgd = dual_and_masses(h_avg)
print(f"SGD:      F(h) = {F_sgd:.5f}  max |mu(L_j) - b_j| = {np.abs(mass_sgd - b).max():.4f}")

Gradient ascent raises F from 0.05908 to 0.06502 and balances the masses to grid precision within 50 iterations. The stochastic version, which never sees the grid, reaches F = 0.06501 with a largest mass error of 0.006 after 20 000 samples.

6. Dynamic formulation and Wasserstein gradient flows

Static optimal transport says where the mass goes; the dynamic view says how it travels. This makes \mathcal{P}_2(\mathbb{R}^d) a geodesic space with a formal Riemannian structure, in which partial differential equations such as the Fokker–Planck equation are gradient flows.

6.1 Displacement interpolation and the Benamou–Brenier formula

A family (\rho_t, v_t) solves the continuity equation \partial_t \rho_t + \mathrm{div}(\rho_t v_t) = 0 when, for every smooth compactly supported \varphi, \frac{d}{dt} \int \varphi \, d\rho_t = \int \langle \nabla\varphi, v_t \rangle \, d\rho_t . If particles follow \dot X_t = v_t(X_t), the law of X_t solves it: differentiate \mathbb{E}\,\varphi(X_t) in time.

With the momentum m = \rho v, the action \int\!\!\int \|m\|^2 / \rho is jointly convex and the constraint is linear, so the formula is a convex problem, which Benamou and Brenier (2000) solve numerically. Otto (2001) read it as a Riemannian metric: the squared norm of a tangent vector \dot\rho is the least kinetic energy \int \|v\|^2 \, d\rho of a velocity that produces it.

6.2 Wasserstein gradient flows

In \mathbb{R}^d, the gradient flow \dot x = -\nabla F(x) is the limit of the implicit Euler scheme x_{k+1} \in \arg\min_x F(x) + \|x - x_k\|^2 / (2\tau). Replacing the Euclidean distance by W_2 gives a scheme on measures.

The formula agrees with Proposition 6.6. The first variation is V + \log\rho + 1, and its gradient times \rho is \rho\nabla V + \nabla\rho, whose divergence is the right-hand side.

For potential and interaction energies, particles give exact solutions. If \dot x_i = -\nabla V(x_i) - \frac1N \sum_j \nabla W(x_i - x_j) , the empirical measure of the x_i solves the equation of Proposition 6.6. The entropy has no particle version, because it is infinite on discrete measures; the Brownian noise of the Langevin diffusion plays its role.

The code runs 100 000 Langevin particles from \mathcal{N}(3, 0.25) and measures W_2 to the standard Gaussian by sorting (Theorem 1.8).

import numpy as np
from scipy.stats import norm
rng = np.random.default_rng(0)

def w2_1d(samples, q):
    """W2 between an empirical measure and a law with quantile function q (Theorem 1.8)."""
    n = len(samples)
    return np.sqrt(np.mean((np.sort(samples) - q((np.arange(n) + 0.5) / n)) ** 2))

# Langevin dX = -X dt + sqrt(2) dB: the law follows the Fokker-Planck equation (Theorem 6.7)
n, dt, m, s = 100_000, 0.01, 3.0, 0.5
X = m + s * rng.standard_normal(n)                        # rho_0 = N(3, 0.25)
for k in range(301):
    t = k * dt
    if k % 100 == 0:
        mt, st = m * np.exp(-t), np.sqrt(s**2 * np.exp(-2*t) + 1 - np.exp(-2*t))
        exact = np.sqrt(mt**2 + (st - 1) ** 2)            # W2(rho_t, N(0,1)), Gaussian formula
        print(f"t={t:.1f}  W2(particles, N(0,1)) = {w2_1d(X, norm.ppf):.3f}"
              f"  exact = {exact:.3f}  e^-t W2(rho_0, N(0,1)) = {np.exp(-t) * np.sqrt(m**2 + (s-1)**2):.3f}")
    X += -X * dt + np.sqrt(2 * dt) * rng.standard_normal(n)

At t = 1, 2, 3 the particles give 1.096, 0.404 and 0.150, against the exact 1.105, 0.406 and 0.149, all below the contraction bound of Proposition 6.8 (1.119, 0.412, 0.151).

7. Optimal transport in machine learning

Four mechanisms of modern machine learning are optimal transport in disguise. Training a wide network is a Wasserstein gradient flow; the tokens of a transformer form an interacting particle system; a GAN minimizes a dual formula; diffusion and flow-matching models build transport paths.

7.1 Training two-layer networks as a Wasserstein gradient flow

A two-layer network with N neurons is an average over particles \theta_i = (a_i, w_i): f_\mu(x) = \int \Phi(\theta, x) \, d\mu(\theta), \qquad \Phi(\theta, x) = a\, \sigma(\langle w, x \rangle), \qquad \mu_N = \frac1N \sum_{i=1}^N \delta_{\theta_i} . The risk R(\mu) = \frac12 \mathbb{E}\, (f_\mu(X) - Y)^2 is convex in \mu, since f_\mu is linear in \mu, even though it is not convex in the \theta_i.

The code trains networks of increasing width, with the mean-field scaling, on a target made of six ReLU units.

import torch
torch.set_num_threads(1)
torch.manual_seed(0)

d, n_train = 2, 300
X = torch.randn(n_train, d)
teacher_w = torch.randn(6, d)
teacher_a = torch.tensor([1.0, -1.0, 1.0, -1.0, 1.0, -1.0])
Y = torch.relu(X @ teacher_w.T) @ teacher_a                   # target: 6 ReLU units

def train(N, steps=1500, lr=0.1, seed=0):
    """Two-layer net f(x) = (1/N) sum_i a_i relu(<w_i, x>): particles theta_i = (a_i, w_i)."""
    g = torch.Generator().manual_seed(seed)
    theta = torch.randn(N, d + 1, generator=g).requires_grad_(True)
    for k in range(steps):
        a, w = theta[:, 0], theta[:, 1:]
        risk = 0.5 * ((torch.relu(X @ w.T) @ a / N - Y) ** 2).mean()
        grad, = torch.autograd.grad(risk, theta)
        with torch.no_grad():
            theta -= lr * N * grad          # each particle follows -grad of dR/dmu (Prop. 7.1)
    return risk.item()

for N in (2, 6, 50, 500):
    risks = [train(N, seed=s) for s in range(3)]
    print(f"N = {N:3d} particles   final risk over 3 seeds: " + "  ".join(f"{r:.4f}" for r in risks))

With 2 particles the final risk depends on the seed (0.18, 3.81, 0.18); with 500 particles all three seeds reach 0.0002, as the mean-field picture predicts.

7.2 Token dynamics in transformers

A self-attention layer with a residual connection updates n tokens x_i \in \mathbb{R}^d by x_i \leftarrow x_i + \sum_{j} A_{ij} \, V x_j, \qquad A_{ij} = \frac{e^{\langle Q x_i, K x_j \rangle}}{\sum_k e^{\langle Q x_i, K x_k \rangle}} . With weights shared across layers and depth read as time, this is an Euler step of an ODE. The velocity of token x depends on the others only through their empirical measure \mu: \dot x = \Gamma_\mu(x), \qquad \Gamma_\mu(x) = \frac{\int e^{\langle Qx, Ky\rangle} \, V y \, d\mu(y)}{\int e^{\langle Qx, Ky \rangle} \, d\mu(y)} .

Geshkovski, Letrouit, Polyanskiy and Rigollet (2023) study Q = K = V = I, with tokens kept on the unit sphere by layer normalization and an inverse temperature \beta in the exponent. They prove that, in dimension d \geq 3 and in several regimes of \beta, almost every initial configuration converges to a single cluster. They also describe how, for large \beta, the convergence passes through long metastable states with several clusters.

Sander, Ablin, Blondel and Peyré (2022) replace the row normalization by Sinkhorn’s algorithm (ch. 4), which makes A doubly stochastic. They show that, for a symmetric attention kernel and unlike the softmax, this normalization lets the attention iterations be read as a discretized Wasserstein gradient flow.

import torch
torch.set_num_threads(1)
torch.manual_seed(0)

def n_clusters(x, tol=0.99):
    """Greedy count of groups of tokens with pairwise inner product above tol."""
    left, count = torch.ones(len(x), dtype=torch.bool), 0
    while left.any():
        i = int(torch.nonzero(left)[0])
        left &= (x @ x[i]) < tol
        count += 1
    return count

def attention_flow(x, beta, dt=0.1, steps=2000, report=(0, 50, 200, 2000)):
    """Tokens on the sphere: x_i' = P_{x_i}^perp( sum_j softmax_j(beta <x_i,x_j>) x_j )."""
    for k in range(steps + 1):
        if k in report:
            print(f"  time {k * dt:6.1f}  clusters = {n_clusters(x)}")
        A = torch.softmax(beta * x @ x.T, dim=1)           # row-stochastic attention
        v = A @ x
        v = v - (v * x).sum(1, keepdim=True) * x          # project on the tangent space
        x = x + dt * v
        x = x / x.norm(dim=1, keepdim=True)               # layer normalisation
    return x

x0 = torch.nn.functional.normalize(torch.randn(64, 3), dim=1)
for beta in (1.0, 9.0):
    print(f"beta = {beta}")
    attention_flow(x0.clone(), beta)

With 64 tokens on the sphere of \mathbb{R}^3: for \beta = 1 there is a single cluster from time 20 on; for \beta = 9, three clusters are still present at time 200.

7.3 GANs: generators, critics and dual formulas

A generator g_\theta pushes a latent law \zeta, for instance \mathcal{N}(0, I), to the model \mu_\theta = (g_\theta)_\sharp \zeta. The model has no density when the latent dimension is smaller than the data dimension, so its likelihood is not available. One minimizes instead a divergence D(\mu_\theta, \nu) to the data law \nu that has a dual formula, estimated with samples:

• Goodfellow et al. (2014): a discriminator \delta with values in (0, 1) reaches \max_\delta \; \mathbb{E}_\nu \log \delta + \mathbb{E}_{\mu_\theta} \log(1 - \delta) = 2\,\mathrm{JS}(\mu_\theta, \nu) - \log 4 .
• Wasserstein GAN (Arjovsky, Chintala and Bottou, 2017): by Proposition 2.5, \min_\theta W_1(\mu_\theta, \nu) = \min_\theta \; \max_{\mathrm{Lip}(f) \leq 1} \; \mathbb{E}_{z \sim \zeta} f\big(g_\theta(z)\big) - \mathbb{E}_{y \sim \nu} f(y) .

The critic f is the Kantorovich potential, and GAN training alternates ascent on the dual variable with descent on the generator.

Example 2.6 explains the gain. When the supports are disjoint, the Jensen–Shannon divergence is constant and gives no gradient, while W_1 gives a gradient of norm one. In practice the Lipschitz constraint is imposed by weight clipping or by a gradient penalty (ch. 8). The FID used to evaluate the samples is the Gaussian formula of Theorem 3.3.

7.4 Diffusion models and flow matching

Diffusion models. The forward noising process of a diffusion model is the Ornstein–Uhlenbeck process of Example 6.9. By Theorem 6.7, the law \rho_t of the noised data follows the Wasserstein gradient flow of \mathrm{KL}(\cdot \,|\, \mathcal{N}(0, I)) and converges exponentially fast.

The Fokker–Planck equation can be rewritten as a continuity equation: \partial_t \rho_t + \mathrm{div}(\rho_t v_t) = 0, \qquad v_t(x) = -x - \nabla \log \rho_t(x) . So the deterministic probability-flow ODE \dot x = v_t(x) has the same marginals. Integrated backwards from \mathcal{N}(0, I), it generates data. The only unknown is the score \nabla\log\rho_t, learned by denoising (Ho, Jain and Abbeel, 2020; Song et al., 2021).

Flow matching. Lipman et al. (2023) and Liu, Gong and Liu (2023) learn the velocity directly. Take a coupling of X_0 \sim \mathcal{N}(0, I) and data X_1, set X_t = (1 - t) X_0 + t X_1, and regress X_1 - X_0 on (X_t, t).

Optimal transport enters through the coupling. With the Brenier coupling, X_t is McCann’s interpolation (Definition 6.1). Paths do not cross, X_1 - X_0 is a function of X_t, and the learned flow is straight, with action W_2^2, the minimum of Theorem 6.3; one Euler step is then exact.

With independent pairs, paths cross, the velocity averages directions, and trajectories bend. Mini-batch OT couplings (Pooladian et al., 2023; Tong et al., 2024) straighten them; Exercise 14 measures the effect.

8. Complements for self-study

master-mva.com announces eight sessions of 3 h, while the official volume used for the session plan is 18 h. This chapter holds the material that the two extra sessions would cover: the proofs left out in class, and the tools that turn the core theory into working machine-learning losses. Each section can be read after the chapter it extends.

8.1 Why Sinkhorn converges: Hilbert’s projective metric

For K = e^{-C/\varepsilon}, \log \eta(K) = \frac1\varepsilon \max_{i, j, k, l} \big( C_{jk} + C_{il} - C_{ik} - C_{jl} \big) , so \lambda(K) \to 1 as \varepsilon \to 0. Two practical remedies follow.

• \varepsilon-scaling: solve with a large \varepsilon, then decrease it, warm-starting the potentials each time (Schmitzer, 2019).
• Accuracy targets: a plan with transport cost within \delta of the optimum is obtained in a number of operations that is n^2 times a polynomial in \|C\|_\infty / \delta, up to logarithmic factors (Altschuler, Weed and Rigollet, 2017).

Sinkhorn is also a special case of iterative Bregman projections: each half-step is the KL projection of the current plan onto one marginal constraint (Benamou et al., 2015). Exercise 9 checks it.

8.2 Sinkhorn divergences: a debiased loss between OT and MMD

The entropic cost of Definition 4.2 is biased: \mathrm{OT}_\varepsilon(\alpha, \alpha) > 0 unless \alpha is a Dirac mass. Minimizing \beta \mapsto \mathrm{OT}_\varepsilon(\alpha, \beta) therefore does not return \alpha; it returns a shrunken version of it.

The last point (Ramdas, García Trillos and Cuturi, 2017; Genevay, Peyré and Cuturi, 2018) explains the role of \varepsilon. Small values keep the geometry of optimal transport; large values give the statistical behaviour of kernel methods.

Estimated from n samples, entropic costs converge at the parametric rate n^{-1/2}, with a constant that blows up as \varepsilon \to 0 and grows with the dimension (Genevay et al., 2019; Mena and Niles-Weed, 2019). Unregularized W_p does much worse, as section 8.3 shows.

Genevay, Peyré and Cuturi (2018) train generators by minimizing a Sinkhorn divergence on mini-batches, with no critic to learn. The code fits an affine generator g(z) = Az + m to a correlated Gaussian. It runs Sinkhorn without building the autograd graph, then performs one differentiable update: by the envelope theorem (Proposition 4.7), this gives the exact gradient at convergence.

import math, torch
torch.set_num_threads(1)
torch.manual_seed(0)

def ot_eps(x, y, eps, iters=50):
    """Entropic OT between uniform point clouds (log-domain Sinkhorn, dual value)."""
    la, lb = -math.log(len(x)), -math.log(len(y))
    softmin = lambda C, h, lw: -eps * torch.logsumexp((h[None, :] - C) / eps + lw, dim=1)
    with torch.no_grad():                                  # iterate without the graph
        C = torch.cdist(x, y) ** 2
        f, g = torch.zeros(len(x)), torch.zeros(len(y))
        for _ in range(iters):
            f = softmin(C, g, lb)
            g = softmin(C.T, f, la)
    # one differentiable update: the envelope theorem gives the exact gradient at convergence
    f = softmin(torch.cdist(x, y) ** 2, g, lb)
    return f.mean() + g.mean()

def sinkhorn_div(x, y, eps):
    """Debiased Sinkhorn divergence S_eps (Definition 8.3)."""
    return ot_eps(x, y, eps) - 0.5 * ot_eps(x, x, eps) - 0.5 * ot_eps(y, y, eps)

m_true = torch.tensor([2.0, -1.0])
L_true = torch.tensor([[1.0, 0.0], [0.8, 0.5]])           # target N(m, L L^T)
A = torch.eye(2, requires_grad=True)                     # generator g(z) = A z + m
m = torch.zeros(2, requires_grad=True)
opt = torch.optim.Adam([A, m], lr=0.05)
sched = torch.optim.lr_scheduler.StepLR(opt, step_size=150, gamma=0.2)
for step in range(301):
    z, data = torch.randn(128, 2), torch.randn(128, 2) @ L_true.T + m_true
    loss = sinkhorn_div(z @ A.T + m, data, eps=0.1)
    opt.zero_grad(); loss.backward(); opt.step(); sched.step()
    if step % 100 == 0:
        print(f"step {step:3d}  S_eps = {loss.item():.4f}")
print("learned mean:", m.detach().numpy().round(2), " target:", m_true.numpy())
print("learned covariance:\n", (A @ A.T).detach().numpy().round(2))
print("target covariance:\n", (L_true @ L_true.T).numpy().round(2))

The learned mean is (1.97, -1.03) and the learned covariance has rows (0.97, 0.78) and (0.78, 0.84), against the target (2, -1) with rows (1, 0.8) and (0.8, 0.89). The remaining gap comes from the mini-batches of 128 points and the finite number of steps. The divergence values printed during training are noisy for the same reason.

8.3 Statistical optimal transport and sliced distances

How fast does the empirical measure \mu_n = \frac1n \sum_i \delta_{X_i} approach \mu in Wasserstein distance? In high dimension, very slowly.

In dimension 2, as in Example 3.5, the rate is still about n^{-1/2}, up to a logarithmic factor; the gap seen there is the upward bias of the plug-in estimate, plus sampling noise. In dimension 10 the rate is n^{-1/10}, and a loss built on sample W_p needs a number of samples exponential in d to reach a fixed accuracy. Weed and Bach (2019) show that the rate depends on the intrinsic dimension of the data rather than the ambient one. Entropic regularization (section 8.2) and slicing are the two standard remedies.

Each slice is a one-dimensional problem, solved by sorting (Theorem 1.8). A Monte Carlo average over L directions costs O(L\, n \log n). Since P_\theta is 1-Lipschitz, SW_2 \leq W_2.

SW_2 is a distance, because a measure is determined by its one-dimensional projections (Cramér–Wold). It was introduced for texture mixing by Rabin, Peyré, Delon and Bernot (2011) and studied for barycenters by Bonneel, Rabin, Peyré and Pfister (2015). In dimension one, empirical measures converge at the rate n^{-1/2} under mild conditions on \mu (Bobkov and Ledoux, 2019), which is what slicing exploits.

8.4 Barycenters and Gaussian mixtures

For instance, the equal-weight barycenter of \mathcal{N}(-2, 0.25) and \mathcal{N}(2, 1) is \mathcal{N}(0, 0.75^2): the means and the standard deviations are averaged. The mixture of the two, with equal weights, is bimodal instead.

In any dimension, the barycenter of Gaussians is Gaussian. Its covariance solves the fixed-point equation \Sigma = \sum_k \lambda_k \big( \Sigma^{1/2} \Sigma_k \Sigma^{1/2} \big)^{1/2} , which the iteration of Álvarez-Esteban, del Barrio, Cuesta-Albertos and Matrán (2016) solves.

For Gaussian mixture models, Delon and Desolneux (2020) restrict the couplings to Gaussian mixtures. The resulting distance MW_2 is a discrete optimal transport problem between the mixture weights, with the ground cost W_2^2 between components given by Theorem 3.3. It is an upper bound on W_2 between the mixtures, cheap to compute, and it combines ch. 3 and ch. 4.

8.5 Wasserstein GANs in practice

The Lipschitz constraint of Proposition 2.5 is the hard part of a Wasserstein GAN. Arjovsky, Chintala and Bottou (2017) clip the critic’s weights; Gulrajani et al. (2017) penalize its gradient instead, guided by the following property.

The WGAN-GP critic loss adds \lambda\, \mathbb{E}\, (\|\nabla f(\hat x)\| - 1)^2, where \hat x is drawn uniformly on segments between generated and real samples. It replaces the coupling of Proposition 8.9, which is unknown, by random pairs.

8.6 Computing Wasserstein gradient flows

• Entropic JKO: replacing W_2^2 by its entropic version turns each JKO step into a KL-proximal problem, solved by Sinkhorn-like scalings on a grid (Peyré, 2015).
• Particles: potential and interaction energies are simulated exactly by particles (ch. 6), and the entropy by Brownian noise (Langevin).
• Mean-field Langevin: adding noise to the training of a two-layer network (ch. 7) gives the gradient flow of the risk plus an entropy. This regularized flow converges to the unique minimizer of the regularized risk (Mei, Montanari and Nguyen, 2018), exponentially fast under log-Sobolev-type conditions (Chizat, 2022; Nitanda, Wu and Suzuki, 2022).

9. Research directions and open problems

This chapter is a self-study tour of current research that builds on the course. It reports results with their authors and years; it is not a complete survey.

9.1 Estimating transport maps in high dimension

The curse of dimensionality of Theorem 8.5 concerns distances. Estimating the Brenier map itself is harder still. Hütter and Rigollet (2021) give minimax rates for maps of given smoothness, which improve with the smoothness and degrade with the dimension.

Entropic estimators, which apply the barycentric projection of Proposition 4.7 to samples, are computationally simple. Pooladian and Niles-Weed (2021) give rates for them. The lecture notes of Chewi, Niles-Weed and Rigollet (2024) survey the field.

A recurring open question is which estimators adapt, at once, to the smoothness of the map and to a low intrinsic dimension of the data (Weed and Bach, 2019).

9.2 Neural optimal transport solvers

Brenier’s theorem suggests parametrizing the potential \varphi by an input convex neural network (Amos, Xu and Kolter, 2017) and solving the dual of Proposition 3.1 by stochastic optimization (Makkuva, Taghvaei, Oh and Lee, 2020).

The benchmark of Korotin et al. (2021), built on pairs of measures with known optimal maps, shows that such solvers can generate good samples while recovering the true map poorly. Training neural OT solvers reliably remains open.

9.3 Transport paths for generative models

• Flow matching with OT couplings. Mini-batch couplings (Pooladian et al., 2023; Tong et al., 2024) straighten the paths of ch. 7.4. Their bias, and its cost at large scale, are active questions.
• Schrödinger bridges. The Schrödinger problem is entropic optimal transport on path space (Léonard, 2014). De Bortoli, Thornton, Heng and Doucet (2021) solve it with score-based models, which gives diffusion models with a prescribed endpoint.
• Is the diffusion map an OT map? Khrulkov et al. (2023) observed that the probability-flow map of a diffusion model is close to the optimal transport map. Lavenant and Santambrogio (2022) proved that the flow map of the Fokker–Planck equation is not optimal in general. For Gaussian data it is (Exercise 15); how close it is for real data is open.

9.4 Mean-field theory of neural networks

The two-layer results of ch. 7.1 raise several questions.

• Which regime? Chizat, Oyallon and Bach (2019) separate the “lazy” regime, where parameters barely move and the network behaves like a kernel method, from the mean-field regime, where features are learned.
• Quantitative rates. Entropic regularization gives exponential convergence (section 8.6). Without it, rates for the noiseless mean-field flow and finite-N corrections are largely open.
• Depth. Extending the measure-valued description beyond two layers is an open research programme.

9.5 Transformers as interacting particles

The clustering results of Geshkovski, Letrouit, Polyanskiy and Rigollet (2023) concern simplified dynamics: shared weights, identity matrices, no MLP block.

• How do learned Q, K, V, causal masks and the MLP layers change the picture?
• What are the time scales of the metastable states, and do they explain how trained transformers behave at finite depth?

Doubly stochastic attention (Sander et al., 2022) offers a variational structure that softmax attention lacks; whether it helps at scale is open.

9.6 Beyond balanced Euclidean transport

• Unbalanced OT relaxes the marginal constraints with divergence penalties, to compare measures of different masses. Its scaling algorithms extend Sinkhorn (Chizat, Peyré, Schmitzer and Vialard, 2018), and its geometry is studied by Liero, Mielke and Savaré (2018).
• Gromov–Wasserstein distances compare metric measure spaces that need not share an ambient space, by matching distances rather than points (Mémoli, 2011). Peyré, Cuturi and Solomon (2016) compute them with entropic regularization.
• Gaussian mixture models get their own transport distance (Delon and Desolneux, 2020; section 8.4).

9.7 A short list of open problems

1. Regularity of the Brenier map (Caffarelli, 1992) under assumptions that real data satisfy, and its use for estimation.
2. Global convergence of the noiseless mean-field dynamics of two-layer networks, with rates and finite-width corrections.
3. Convergence of GAN training, a min–max problem in which the dual variable is itself learned.
4. The best computational complexity of exact and approximate discrete optimal transport.
5. The time scales of metastable clusters in attention dynamics, and their role in trained models.
6. Straight generative flows at scale: unbiased OT couplings that remain cheap with mini-batches.

Exercises

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

1. ★ Monotone transport by hand. Let \mu = \frac12 \delta_0 + \frac12 \delta_2 and \nu = \frac13 (\delta_1 + \delta_3 + \delta_4). Compute the optimal plan, W_1(\mu, \nu) and W_2^2(\mu, \nu) with Theorem 1.8.

   Solution: the quantile function of \mu is 0 on (0, \frac12] and 2 after; that of \nu is 1, 3, 4 on the thirds of (0, 1].

   The points \frac13, \frac12, \frac23 cut the unit interval into four pieces. On them, the pairs are (0, 1), (0, 3), (2, 3), (2, 4), with masses \frac13, \frac16, \frac16, \frac13: this is the optimal plan.

   Integrating the distances and their squares against these masses gives W_1 = \tfrac13 + \tfrac36 + \tfrac16 + \tfrac23 = \tfrac53, \qquad W_2^2 = \tfrac13 + \tfrac96 + \tfrac16 + \tfrac43 = \tfrac{10}{3} .

2. ★ Book shifting. Let \mu be uniform on [0, 1] and \nu uniform on [1, 2].

   Compare the translation x \mapsto x + 1 and the flip x \mapsto 2 - x for the costs |x - y| and |x - y|^2.

   Solution: both maps push \mu to \nu. For |x - y|, the translation costs 1 and the flip costs \int_0^1 (2 - 2x) \, dx = 1. Since W_1 = 1 by Theorem 1.8, both are optimal: optimal maps need not be unique for W_1.

   For the quadratic cost, the translation costs 1 and the flip \int_0^1 (2 - 2x)^2 \, dx = \frac43. Only the translation, which is non-decreasing, is optimal, in line with Theorem 3.2.

3. ★ A Brenier map on the line. Let U be uniform on [0, 1] and let \nu be the law of U^2.

   Find the density of \nu. Show that T(x) = x^2 is the optimal map from the uniform law to \nu, and compute W_2^2.

   Solution: the Jacobian equation (ch. 0) gives the density \frac{1}{2\sqrt{y}} on (0, 1]. On [0, 1], T is the derivative of the convex function x^3/3, extended by 0 on the negative half-line; by the converse of Theorem 3.2 it is optimal. Then W_2^2 = \int_0^1 (x - x^2)^2 \, dx = \tfrac13 - \tfrac12 + \tfrac15 = \tfrac{1}{30} .

4. ★ Bures distance by hand. (a) Compute the optimal map and W_2^2 between \mathcal{N}(0, \mathrm{diag}(1, 4)) and \mathcal{N}((1, 0), \mathrm{diag}(4, 1)).

   2. Compute W_2^2 between \mathcal{N}(0, I_2) and \mathcal{N}(0, \Sigma), where \Sigma has rows (2, 1) and (1, 2).

   Solution: (a) the covariances commute, so A = \mathrm{diag}(2, \frac12) and T(x) = (1, 0) + A x. Definition 3.4 gives W_2^2 = 1 + (1 - 2)^2 + (2 - 1)^2 = 3 .

   2. I commutes with \Sigma, whose eigenvalues are 3 and 1. Hence W_2^2 = (1 - \sqrt3)^2 + (1 - 1)^2 = 4 - 2\sqrt3 \approx 0.536.

5. ★ Sinkhorn on a 2 × 2 problem. Let a = b = (\frac12, \frac12) and let C have rows (0, 1) and (1, 0). Compute P^\varepsilon, its transport cost, and the number of Sinkhorn iterations needed from v = \mathbf{1}.

   Solution: by symmetry P^\varepsilon has p on the diagonal and \frac12 - p off it. The Gibbs form (Proposition 4.3) with constant potentials gives the ratio p / (\frac12 - p) = e^{1/\varepsilon}, so p = \frac{1}{2 (1 + e^{-1/\varepsilon})}, \qquad \langle C, P^\varepsilon \rangle = 1 - 2p = \frac{1}{1 + e^{1/\varepsilon}} . This is 0.269 for \varepsilon = 1, below the bound \log 2 \approx 0.693 of Proposition 4.5, and exponentially small as \varepsilon \to 0.

   One iteration suffices: from v = \mathbf{1}, u is a constant vector, and the next v is again \mathbf{1}.

6. ★ The price of regularization (Python). With the code of ch. 4, plot \langle C, P^\varepsilon \rangle - \mathrm{OT} against \varepsilon on a log–log scale, together with the bound \varepsilon \log n. Then count the iterations needed to reach a marginal error of 10^{-6} for \varepsilon \in \{1, 0.3, 0.1, 0.03, 0.01\}.

   Hints and expected result: the gap stays well below the bound and decreases roughly linearly in \varepsilon. With the point clouds of the code, the iteration counts are about 18, 66 and 280 for \varepsilon = 1, 0.3, 0.1, then tens of thousands for \varepsilon = 0.03 and 0.01: Theorem 4.6 in action. \varepsilon-scaling (ch. 8.1) reduces them.

7. ★★ W_1 on the line. Show that for measures with bounded support on \mathbb{R}, W_1(\mu, \nu) = \int_{\mathbb{R}} |F_\mu(t) - F_\nu(t)| \, dt . Check it on Exercise 1.

   Solution: for a 1-Lipschitz f, integration by parts gives \int f \, d(\mu - \nu) = -\int f'(t)\, (F_\mu - F_\nu)(t) \, dt. Since |f'| \leq 1, this is at most \int |F_\mu - F_\nu|, with equality for f' = -\mathrm{sign}(F_\mu - F_\nu). Proposition 2.5 concludes.

   In Exercise 1, |F_\mu - F_\nu| equals \frac12, \frac16, \frac23, \frac13 on the unit intervals between 0 and 4. The sum is \frac53, as found there.

8. ★★ The Gelbrich bound. Let \mu, \nu \in \mathcal{P}_2(\mathbb{R}^d) have means m, m' and covariances \Sigma, \Sigma'. Show that W_2^2(\mu, \nu) \geq \|m - m'\|^2 + \mathcal{B}(\Sigma, \Sigma')^2 , with equality for Gaussians.

   Solution: for a coupling (X, Y) with cross-covariance K, \mathbb{E}\,\|X - Y\|^2 = \|m - m'\|^2 + \mathrm{tr}\,\Sigma + \mathrm{tr}\,\Sigma' - 2\, \mathrm{tr}\, K . The joint covariance is positive semidefinite, which forces K = \Sigma^{1/2} R\, \Sigma'^{1/2} with \|R\|_{\mathrm{op}} \leq 1.

   Hence \mathrm{tr}\, K \leq \|\Sigma'^{1/2} \Sigma^{1/2}\|_*, the nuclear norm, which equals \mathrm{tr} (\Sigma^{1/2} \Sigma' \Sigma^{1/2})^{1/2}. Theorem 3.3 gives equality for Gaussians (Gelbrich, 1990).

9. ★★ Sinkhorn as alternating KL projections. For a positive matrix Q, show that the solution of \min \{ \mathrm{KL}(P \,|\, Q) : P \mathbf{1} = a \} is \mathrm{diag}(a \oslash Q\mathbf{1})\, Q. Deduce that each half-step of Algorithm 4.4 is such a projection.

   Solution: the Lagrangian condition \log(P_{ij}/Q_{ij}) + 1 = \lambda_i gives P_{ij} = Q_{ij}\, e^{\lambda_i - 1}, a row rescaling. The constraint fixes the factor of row i to a_i / (Q\mathbf{1})_i.

   With Q = \mathrm{diag}(u) K \mathrm{diag}(v), the new row factors give u \leftarrow a \oslash (Kv), the first Sinkhorn update. The column update is the projection onto \{P^\top \mathbf{1} = b\} (Benamou et al., 2015).

10. ★★ Semi-discrete transport on a segment. Let \mu be uniform on [0, 1] and \nu = \sum_j b_j \delta_{y_j} with y_1 < \dots < y_m.

    Show that the optimal cells are the intervals [s_{j-1}, s_j], where s_j = b_1 + \dots + b_j.

    Find g for y = (0, \frac12, 1) and b = (\frac12, \frac14, \frac14).

    Solution: the optimal map is non-decreasing (Theorem 1.8), so the cells are consecutive intervals of lengths b_j. As in Example 5.5, the boundary between cells j and j + 1 is s_j = \frac{y_j + y_{j+1}}{2} + \frac{g_j - g_{j+1}}{2\, (y_{j+1} - y_j)} . Here s_1 = \frac12 and s_2 = \frac34 give g_1 - g_2 = \frac14 and g_2 - g_3 = 0, up to an additive constant.

11. ★★ Interpolating Gaussians on the line. Compute the displacement interpolation (Definition 6.1) between \mathcal{N}(m_0, \sigma_0^2) and \mathcal{N}(m_1, \sigma_1^2), and check Proposition 6.2. Compare its variance at t = \frac12 with that of the mixture, for m_0 = -2, m_1 = 2 and \sigma_0 = \sigma_1 = \frac12.

    Solution: the Brenier map is T(x) = m_1 + \frac{\sigma_1}{\sigma_0}(x - m_0), so \mu_t is Gaussian with mean (1 - t) m_0 + t m_1 and standard deviation (1 - t)\sigma_0 + t \sigma_1.

    The means of \mu_s and \mu_t differ by (t - s)(m_1 - m_0), and their standard deviations by (t - s)(\sigma_1 - \sigma_0). Definition 3.4 gives Proposition 6.2: W_2(\mu_s, \mu_t) = |t - s| \sqrt{(m_1 - m_0)^2 + (\sigma_1 - \sigma_0)^2} . At t = \frac12, \mu_t = \mathcal{N}(0, \frac14), while the mixture has variance \frac14 + 4 = 4.25 and two modes.

12. ★★ One JKO step for a potential energy. Let F(\rho) = \int V \, d\rho. Show that the JKO step (Definition 6.4) is \rho_{k+1} = (S_\tau)_\sharp \rho_k, where S_\tau(x) minimizes V(y) + \frac{1}{2\tau}\|y - x\|^2. Make it explicit for V(x) = \frac{\lambda}{2}\|x\|^2 and let \tau \to 0.

    Solution: take any \rho and an optimal coupling \pi of \rho_k and \rho. Since the second marginal of \pi is \rho, the JKO objective satisfies \int \Big[ V(y) + \frac{1}{2\tau}\|x - y\|^2 \Big] \, d\pi(x, y) \;\geq\; \int \min_{y} \Big[ V(y) + \frac{1}{2\tau}\|x - y\|^2 \Big] \, d\rho_k(x) . The measure (S_\tau)_\sharp\rho_k, with the coupling (\mathrm{Id}, S_\tau)_\sharp\rho_k, reaches this lower bound.

    For the quadratic potential, S_\tau(x) = x / (1 + \tau\lambda). After t / \tau steps, each point is multiplied by (1 + \tau\lambda)^{-t/\tau} \to e^{-\lambda t}: particles follow \dot x = -\nabla V(x).

13. ★★ Doubly stochastic attention (Python). Modify the code of ch. 7.2 to normalize the attention matrix with 50 log-domain Sinkhorn iterations instead of the softmax. Compare the number of clusters over time for \beta = 1 and \beta = 9.

    Hints and expected result: replace torch.softmax by alternate row and column normalizations of beta * x @ x.T with torch.logsumexp. Both normalizations cluster the tokens. With the seed of the code, the softmax gives 12 then 3 clusters at times 5 and 20 for \beta = 9; the Sinkhorn version gives 18 then 4, and both give 3 clusters at time 200. Relate the slower merging to the doubly stochastic constraint (Sander et al., 2022).

14. ★★ Flow matching with and without OT couplings (Python). Train a small MLP velocity field v(x, t) by flow matching (Proposition 7.5), from \mathcal{N}(0, I_2) to a mixture of 8 Gaussians of standard deviation 0.3 on the circle of radius 4. Use either independent pairs or pairs coupled by an exact assignment within each mini-batch of 256. Generate with K = 1, 2, 4, 16 Euler steps and measure W_2 to 1000 test samples.

    Hints and expected result: use linear_sum_assignment for the coupling and for the evaluation. With an MLP of three hidden layers of width 256 and 3000 Adam steps at learning rate 10^{-3}, W_2 is about 3.5, 1.5–1.9, 0.7–0.8 and 0.5 with independent pairs, and about 0.75, 0.55–0.65, 0.55 and 0.5 with OT pairs; the exact values depend on the network and the seed.

    OT couplings make one-step generation possible. The common floor near 0.5 is the error of the evaluation itself: two independent 1000-point samples of the target are already at W_2 \approx 0.45–0.7, mostly because the eight modes receive slightly different numbers of points. Compute this baseline before reading the other numbers.

15. ★★★ Is the diffusion map optimal? Let \rho_0 = \mathcal{N}(0, \Sigma_0) evolve by the Ornstein–Uhlenbeck process, and let X_t follow the probability-flow ODE of ch. 7.4 from X_0 \sim \rho_0. Show that X_0 \mapsto X_t is the optimal map from \rho_0 to \rho_t. Then discuss non-Gaussian data.

    Hints and expected result: \rho_t = \mathcal{N}(0, \Sigma_t) with \Sigma_t = e^{-2t} \Sigma_0 + (1 - e^{-2t}) I, and the velocity is v_t(x) = (\Sigma_t^{-1} - I)\, x. All the \Sigma_t are functions of \Sigma_0, so they commute, and the flow map is M_t = \exp \int_0^t (\Sigma_s^{-1} - I) \, ds, symmetric positive definite.

    By the converse of Theorem 3.2, x \mapsto M_t x is optimal. In dimension one every flow map is increasing, hence optimal, so counterexamples need d \geq 2. Lavenant and Santambrogio (2022) construct one for the Fokker–Planck flow; measuring the gap on real data is open (ch. 9.3).

16. ★★★ The curse of dimensionality, measured (Python). For d \in \{1, 2, 5, 10\} and n \in \{64, \dots, 1024\}, estimate \mathbb{E}\, W_2^2 between two independent n-samples of the uniform law on [0, 1]^d with linear_sum_assignment, and fit the slope in log–log scale. Repeat with SW_2^2 (50 random directions) and with the Sinkhorn divergence; explain the results with Theorem 8.5.

    Hints and expected result: with 3 repetitions per point, the fitted slopes of W_2^2 are about −0.9, −0.45 and −0.25 for d = 2, 5, 10, close to -2/d. The sliced slopes stay near -1 in every dimension. Discuss what this implies for losses built on sample Wasserstein distances.

Projects and assessment

The official assessment, on master-mva.com, asks students to attend all courses and do the numerical parts, which are not graded. The grade comes from a final written exam, with no supporting documents allowed, which reflects both the numerical and the exercise components. The Université Paris-Saclay course detail describes the assessment as “Evaluation Continue non Intégrale”.

No project is graded in the official format. The subjects below extend the numerical sessions for students who want to go further, or prepare the research internship. The written exam is a model of the final exam, and the rubric applies to the projects.

Solo projects

1. [Solo · ≈ 20 h] Sinkhorn divergences between optimal transport and MMD. Does the debiased Sinkhorn divergence really interpolate between W_1 and the energy distance, and does the bias of \mathrm{OT}_\varepsilon matter for gradient flows? A paper study of Feydy, Séjourné, Vialard, Amari, Trouvé and Peyré (2019).
   • Implement: log-domain Sinkhorn with \varepsilon-scaling (ch. 8.1); S_\varepsilon for c = \|x - y\|; the energy distance, which is the MMD with kernel -\|x - y\|.
   • Data: synthetic point clouds in 2D (a Gaussian, a ring, two moons from sklearn.datasets.make_moons).
   • Experiments: S_\varepsilon against \varepsilon from 10^{-2} to 10^{2}; particle gradient flows of S_\varepsilon and of \mathrm{OT}_\varepsilon toward a fixed target.
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters: ch. 4, 8.1–8.2.
   • Success: S_\varepsilon within 5 % of the exact W_1 (linear program of ch. 1) at the smallest \varepsilon and of half the energy distance at the largest; the flow of S_\varepsilon matches the target variance within 5 %, while \mathrm{OT}_\varepsilon shrinks it.
2. [Solo · ≈ 20 h] A transport distance between Gaussian mixtures. How far is the cheap distance MW_2 above the true W_2 between two mixtures? A paper study of Delon and Desolneux (2020).
   • Implement: MW_2 as the discrete OT problem between mixture weights with the ground cost of Theorem 3.3; the Gelbrich lower bound of Exercise 8; W_2 between samples of 2000 points by exact assignment.
   • Data: mixtures fitted with sklearn.mixture.GaussianMixture to synthetic data in dimensions 2 and 5, with 2 to 6 components.
   • Deliverables: code, a 4–6 page report with the geodesic between two mixtures drawn at five times, a 10-minute talk. Chapters: ch. 3, 4, 6.1, 8.4.
   • Success: the chain Gelbrich ≤ W_2 ≤ MW_2 checked on 20 random pairs, up to the sampling error measured by repetition; the ratio MW_2 / W_2 plotted against the separation of the components.
3. [Solo · ≈ 15 h] The unadjusted Langevin algorithm in Wasserstein distance. How much does discretizing the Langevin diffusion of Theorem 6.7 cost? A proof-and-experiment project after Dalalyan (2017) and Durmus and Moulines (2017).
   • Prove: for V with \lambda I \preceq D^2 V \preceq L I, the chain below contracts W_2 by the factor \max(|1 - h\lambda|, |1 - hL|) per step, by synchronous coupling as in Proposition 6.8: x_{k+1} = x_k - h \nabla V(x_k) + \sqrt{2h}\, \xi_k .
   • Compute: for V(x) = \lambda \|x\|^2 / 2, the stationary law \mathcal{N}(0, \frac{2}{\lambda(2 - h\lambda)} I) and its exact W_2 bias (Definition 3.4).
   • Experiment: the bias against h in 1D by sorting (Theorem 1.8), for the Gaussian and for the double well x^4/4 - x^2/2, with a quadrature of e^{-V} as reference.
   • Deliverables: a written proof, code, a 4-page report. Chapters: ch. 1.4, 3.2, 6.2. Success: the Gaussian bias matches the closed form within Monte Carlo error, with log–log slope 1 in h.

Group projects

1. [Group of 3 · ≈ 30 h each] Straight generative flows with mini-batch OT couplings. How much do OT couplings straighten flow-matching paths, and what does the mini-batch bias cost? A paper study of Tong et al. (2024), with Pooladian et al. (2023).
   • Implement: flow matching (Proposition 7.5) with independent pairs, exact mini-batch assignment and mini-batch Sinkhorn couplings; Euler sampling with K steps.
   • Data: 8 Gaussians on a circle, two moons and a 2D Swiss roll (sklearn.datasets); a Gaussian pair, whose Brenier map is known by Theorem 3.3, as sanity check.
   • Split: member A writes the couplings and the Gaussian check; member B the velocity networks and samplers; member C the metrics (exact W_2 to test samples, SW_2, path straightness against the action of Theorem 6.3).
   • Deliverables: repository, an 8-page report, a 15-minute defence. Chapters: ch. 3, 4, 6, 7.4.
   • Success: one-step generation with OT couplings at least twice as close in W_2 as with independent pairs, as in Exercise 14; the bias plotted against the batch size.
2. [Group of 3 · ≈ 30 h each] Critics and losses for Wasserstein GANs. Does the critic of a WGAN estimate W_1, and which way of imposing the Lipschitz constraint works best? A paper study of Arjovsky, Chintala and Bottou (2017) and Gulrajani et al. (2017), with the critic-free loss of Genevay, Peyré and Cuturi (2018).
   • Implement: weight clipping, gradient penalty (ch. 8.5) and a Sinkhorn-divergence generator (ch. 8.2), all with small MLPs.
   • Data: 8 and 25 Gaussians in 2D, and the 2D Swiss roll.
   • Split: one member per method; each one also owns one shared metric: the exact W_1 by the LP of ch. 1 on 500 samples, the exact W_2 by assignment, and the critic’s gradient norm along transport segments (Proposition 8.9).
   • Deliverables: repository, an 8-page report, a 15-minute defence. Chapters: ch. 2, 4, 7.3, 8.2, 8.5.
   • Success: on a fixed pair of measures, the gradient-penalty critic estimates W_1 within 10 % of the LP value; the three generators compared by W_2 over 3 seeds.
3. [Group of 2 · ≈ 25 h each] Clusters and metastability in attention dynamics. When do tokens driven by self-attention form one cluster, and how long do multi-cluster states last? A paper study of Geshkovski, Letrouit, Polyanskiy and Rigollet (2023), with Sander, Ablin, Blondel and Peyré (2022).
   • Implement: the sphere dynamics of ch. 7.2 for 32 to 256 tokens in dimensions 2, 3 and 8, with softmax and with Sinkhorn normalization (Exercise 13).
   • Split: member A writes the simulations and the cluster detection; member B proves that for \beta = 0 the dynamics is a projected gradient ascent of \|\bar x\|^2, where \bar x is the mean token, and analyses the plateaus.
   • Deliverables: code, a 6-page report, a 15-minute joint talk. Chapters: ch. 4, 6.2, 7.2.
   • Success: the proof for \beta = 0; the number of clusters against \log t for five values of \beta; the plateau length plotted against \beta.

Written exam (proposal)

No documents allowed, as announced on master-mva.com. Proposed format: three hours, no electronic device, four independent questions, 20 points. The last question covers the numerical sessions, as the course description announces.

1. Question (4 points). State the Monge and Kantorovich problems between two discrete measures. Give an example where Monge’s problem has no admissible map. Write the dual of the discrete Kantorovich problem and the complementary slackness conditions.

   Expected:

   • Definitions 1.1 and 1.3, with the transport polytope U(a, b) (1 pt).
   • One Dirac mass against two Dirac masses: a map sends an atom to a single point (Example 1.2) (1 pt).
   • The dual of Theorem 2.1, with constraints f_i + g_j \leq C_{ij} (1 pt).
   • Slackness: P_{ij} > 0 implies f_i + g_j = C_{ij}, with the shipping-company reading (1 pt).

2. Question (5 points). Let \mu_0 = \mathcal{N}(0, \mathrm{diag}(4, 1)) and \mu_1 = \mathcal{N}((3, 0), \mathrm{diag}(1, 9)).

   1. Give the optimal map from \mu_0 to \mu_1 and justify its optimality. (b) Compute W_2^2(\mu_0, \mu_1). (c) Describe the displacement interpolation at t = \frac12 and compute W_2(\mu_0, \mu_{1/2}).

   Expected:

   • 1. The covariances commute, so A = \mathrm{diag}(\frac12, 3) and T(x) = (3, 0) + A x. T is the gradient of a convex quadratic function, hence optimal by the converse of Theorem 3.2 (2 pts).
   • 2. By Definition 3.4 (1.5 pts), W_2^2 = 9 + (2 - 1)^2 + (1 - 3)^2 = 14 .
   • 3. \mu_{1/2} is Gaussian with mean (1.5, 0) and standard deviations 1.5 and 2. By Proposition 6.2, W_2(\mu_0, \mu_{1/2}) = \sqrt{14} / 2 \approx 1.87 (1.5 pts).

3. Question (5 points). For the entropic problem of Definition 4.2: (a) show that, for fixed g, the maximizer of the dual D_\varepsilon in f is the Sinkhorn update of Algorithm 4.4; (b) prove that the transport cost of the entropic plan exceeds the optimal cost by at most \varepsilon \log \min(n, m).

   Expected:

   • 1. D_\varepsilon is concave and separable in the f_i. Setting the partial derivative in f_i to zero gives a_i = a_i \, e^{f_i / \varepsilon} \sum_j b_j \, e^{(g_j - C_{ij})/\varepsilon} , hence the log-sum-exp formula (2 pts).
   • 2. Compare the entropic objectives at P^\varepsilon and at an unregularized optimum P^0; KL is non-negative, and \mathrm{KL}(P^0 | a \otimes b) is a mutual information, at most \min(H(a), H(b)) \leq \log \min(n, m). This is Proposition 4.5 (3 pts).

4. Question (6 points). A student writes the following solver.

   import numpy as np
   
   def solve(a, b, C, eps, iters):
       K = np.exp(-C / eps)
       v = np.ones(len(b))
       for _ in range(iters):
           u = a / (K @ v)
           v = b / (K.T @ u)
       return u[:, None] * K * v[None, :]

   1. Which problem does it solve, and which marginal constraint holds exactly when it returns? (b) Why does it return nan for costs of order 100 and \varepsilon = 0.1, and how is it fixed? (c) The student needs the gradient of the entropic cost with respect to the source positions x_i, for C_{ij} = \|x_i - y_j\|^2. How can it be computed without differentiating through the loop?

   Expected:

   • 1. Entropic OT (Definition 4.2) by Sinkhorn’s matrix scaling. The last update is in v, so the column sums equal b exactly; the row sums converge linearly (Theorem 4.6) (2 pts).
   • 2. e^{-1000} underflows to 0, so K v vanishes and a / 0 appears. The fix is the log-domain iteration of Algorithm 4.4 with log-sum-exp, and \varepsilon-scaling for speed (2 pts).
   • 3. Envelope theorem (Proposition 4.7): the gradient is 2 a_i (x_i - T_\varepsilon(x_i)), with the barycentric projection of the returned plan (2 pts).

Grading rubric (proposal)

References

Official references (bibliography of the master-mva.com course page).

• G. Peyré, Optimal Transport for Machine Learners, course notes, 2025 — PDF.
• G. Peyré, M. Cuturi, Computational Optimal Transport, Foundations and Trends in Machine Learning, Now Publishers, 2019 — arXiv:1803.00567.
• F. Santambrogio, Optimal Transport for Applied Mathematicians: Calculus of Variations, PDEs, and Modeling, Birkhäuser, 2015 — PDF.

Textbooks and surveys

• C. Villani, Topics in Optimal Transportation, Graduate Studies in Mathematics 58, American Mathematical Society, 2003.
• C. Villani, Optimal Transport: Old and New, Grundlehren der mathematischen Wissenschaften 338, Springer, 2009.
• L. Ambrosio, N. Gigli, G. Savaré, Gradient Flows in Metric Spaces and in the Space of Probability Measures, Birkhäuser, 2005.
• Q. Mérigot, B. Thibert, “Optimal transport: discretization and algorithms”, in Handbook of Numerical Analysis, vol. 22, Elsevier, 2021.
• P. Billingsley, Convergence of Probability Measures, 2nd ed., Wiley, 1999 — ch. 0.
• R. T. Rockafellar, Convex Analysis, Princeton University Press, 1970 — ch. 0.
• D. Bertsimas, J. N. Tsitsiklis, Introduction to Linear Optimization, Athena Scientific, 1997 — ch. 0, 4.

Foundations (ch. 1–3, 6)

• G. Monge, “Mémoire sur la théorie des déblais et des remblais”, Histoire de l’Académie Royale des Sciences de Paris, 1781.
• L. Kantorovich, “On the translocation of masses”, Doklady Akademii Nauk SSSR, 1942.
• Y. Brenier, “Polar factorization and monotone rearrangement of vector-valued functions”, Communications on Pure and Applied Mathematics, 1991.
• M. Gelbrich, “On a formula for the L2 Wasserstein metric between measures on Euclidean and Hilbert spaces”, Mathematische Nachrichten, 1990.
• L. Caffarelli, “The regularity of mappings with a convex potential”, Journal of the American Mathematical Society, 1992.
• R. J. McCann, “A convexity principle for interacting gases”, Advances in Mathematics, 1997.
• J.-D. Benamou, Y. Brenier, “A computational fluid mechanics solution to the Monge–Kantorovich mass transfer problem”, Numerische Mathematik, 2000.
• R. Jordan, D. Kinderlehrer, F. Otto, “The variational formulation of the Fokker–Planck equation”, SIAM Journal on Mathematical Analysis, 1998.
• F. Otto, “The geometry of dissipative evolution equations: the porous medium equation”, Communications in Partial Differential Equations, 2001.
• H. Lavenant, F. Santambrogio, “The flow map of the Fokker–Planck equation does not provide optimal transport”, Applied Mathematics Letters, 2022.

Algorithms (ch. 4, 5, 8)

• H. W. Kuhn, “The Hungarian method for the assignment problem”, Naval Research Logistics Quarterly, 1955.
• D. P. Bertsekas, “The auction algorithm: a distributed relaxation method for the assignment problem”, Annals of Operations Research, 1988.
• G. Birkhoff, “Extensions of Jentzsch’s theorem”, Transactions of the American Mathematical Society, 1957.
• J. Franklin, J. Lorenz, “On the scaling of multidimensional matrices”, Linear Algebra and its Applications, 1989.
• M. Cuturi, “Sinkhorn distances: lightspeed computation of optimal transport”, Advances in Neural Information Processing Systems, 2013.
• J.-D. Benamou, G. Carlier, M. Cuturi, L. Nenna, G. Peyré, “Iterative Bregman projections for regularized transportation problems”, SIAM Journal on Scientific Computing, 2015.
• J. Altschuler, J. Weed, P. Rigollet, “Near-linear time approximation algorithms for optimal transport via Sinkhorn iteration”, Advances in Neural Information Processing Systems, 2017.
• B. Schmitzer, “Stabilized sparse scaling algorithms for entropy regularized transport problems”, SIAM Journal on Scientific Computing, 2019.
• F. Aurenhammer, F. Hoffmann, B. Aronov, “Minkowski-type theorems and least-squares clustering”, Algorithmica, 1998.
• B. Lévy, “A numerical algorithm for L2 semi-discrete optimal transport in 3D”, ESAIM: Mathematical Modelling and Numerical Analysis, 2015.
• J. Kitagawa, Q. Mérigot, B. Thibert, “Convergence of a Newton algorithm for semi-discrete optimal transport”, Journal of the European Mathematical Society, 2019.
• A. Genevay, M. Cuturi, G. Peyré, F. Bach, “Stochastic optimization for large-scale optimal transport”, Advances in Neural Information Processing Systems, 2016.
• G. Peyré, “Entropic approximation of Wasserstein gradient flows”, SIAM Journal on Imaging Sciences, 2015.

Statistics, barycenters and Gaussian mixtures (ch. 8)

• R. M. Dudley, “The speed of mean Glivenko–Cantelli convergence”, Annals of Mathematical Statistics, 1969.
• N. Fournier, A. Guillin, “On the rate of convergence in Wasserstein distance of the empirical measure”, Probability Theory and Related Fields, 2015.
• J. Weed, F. Bach, “Sharp asymptotic and finite-sample rates of convergence of empirical measures in Wasserstein distance”, Bernoulli, 2019.
• S. Bobkov, M. Ledoux, “One-dimensional empirical measures, order statistics, and Kantorovich transport distances”, Memoirs of the American Mathematical Society, 2019.
• A. Ramdas, N. García Trillos, M. Cuturi, “On Wasserstein two-sample testing and related families of nonparametric tests”, Entropy, 2017.
• A. Genevay, G. Peyré, M. Cuturi, “Learning generative models with Sinkhorn divergences”, AISTATS, 2018.
• J. Feydy, T. Séjourné, F.-X. Vialard, S. Amari, A. Trouvé, G. Peyré, “Interpolating between optimal transport and MMD using Sinkhorn divergences”, AISTATS, 2019.
• A. Genevay, L. Chizat, F. Bach, M. Cuturi, G. Peyré, “Sample complexity of Sinkhorn divergences”, AISTATS, 2019.
• G. Mena, J. Niles-Weed, “Statistical bounds for entropic optimal transport: sample complexity and the central limit theorem”, Advances in Neural Information Processing Systems, 2019.
• J. Rabin, G. Peyré, J. Delon, M. Bernot, “Wasserstein barycenter and its application to texture mixing”, Scale Space and Variational Methods in Computer Vision (SSVM 2011), Lecture Notes in Computer Science, Springer, 2012.
• N. Bonneel, J. Rabin, G. Peyré, H. Pfister, “Sliced and Radon Wasserstein barycenters of measures”, Journal of Mathematical Imaging and Vision, 2015.
• M. Agueh, G. Carlier, “Barycenters in the Wasserstein space”, SIAM Journal on Mathematical Analysis, 2011.
• P. C. Álvarez-Esteban, E. del Barrio, J. A. Cuesta-Albertos, C. Matrán, “A fixed-point approach to barycenters in Wasserstein space”, Journal of Mathematical Analysis and Applications, 2016.
• J. Delon, A. Desolneux, “A Wasserstein-type distance in the space of Gaussian mixture models”, SIAM Journal on Imaging Sciences, 2020.

Machine learning (ch. 7, 8)

• I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, Y. Bengio, “Generative adversarial nets”, Advances in Neural Information Processing Systems, 2014.
• M. Arjovsky, S. Chintala, L. Bottou, “Wasserstein generative adversarial networks”, ICML, 2017.
• I. Gulrajani, F. Ahmed, M. Arjovsky, V. Dumoulin, A. Courville, “Improved training of Wasserstein GANs”, Advances in Neural Information Processing Systems, 2017.
• M. Heusel, H. Ramsauer, T. Unterthiner, B. Nessler, S. Hochreiter, “GANs trained by a two time-scale update rule converge to a local Nash equilibrium”, Advances in Neural Information Processing Systems, 2017 — the FID.
• S. Mei, A. Montanari, P.-M. Nguyen, “A mean field view of the landscape of two-layer neural networks”, Proceedings of the National Academy of Sciences, 2018.
• L. Chizat, F. Bach, “On the global convergence of gradient descent for over-parameterized models using optimal transport”, Advances in Neural Information Processing Systems, 2018.
• L. Chizat, E. Oyallon, F. Bach, “On lazy training in differentiable programming”, Advances in Neural Information Processing Systems, 2019.
• J. Sirignano, K. Spiliopoulos, “Mean field analysis of neural networks: a law of large numbers”, SIAM Journal on Applied Mathematics, 2020.
• L. Chizat, “Mean-field Langevin dynamics: exponential convergence and annealing”, Transactions on Machine Learning Research, 2022.
• A. Nitanda, D. Wu, T. Suzuki, “Convex analysis of the mean field Langevin dynamics”, AISTATS, 2022.
• M. E. Sander, P. Ablin, M. Blondel, G. Peyré, “Sinkformers: Transformers with doubly stochastic attention”, AISTATS, 2022.
• B. Geshkovski, C. Letrouit, Y. Polyanskiy, P. Rigollet, “A mathematical perspective on Transformers”, arXiv preprint, 2023.
• J. Ho, A. Jain, P. Abbeel, “Denoising diffusion probabilistic models”, Advances in Neural Information Processing Systems, 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.
• Y. Lipman, R. T. Q. Chen, H. Ben-Hamu, M. Nickel, M. Le, “Flow matching for generative modeling”, ICLR, 2023.
• X. Liu, C. Gong, Q. Liu, “Flow straight and fast: learning to generate and transfer data with rectified flow”, ICLR, 2023.
• A.-A. Pooladian, H. Ben-Hamu, C. Domingo-Enrich, B. Amos, Y. Lipman, R. T. Q. Chen, “Multisample flow matching: straightening flows with minibatch couplings”, ICML, 2023.
• A. Tong, K. Fatras, N. Malkin, G. Huguet, Y. Zhang, J. Rector-Brooks, G. Wolf, Y. Bengio, “Improving and generalizing flow-based generative models with minibatch optimal transport”, Transactions on Machine Learning Research, 2024.
• A. S. Dalalyan, “Theoretical guarantees for approximate sampling from smooth and log-concave densities”, Journal of the Royal Statistical Society: Series B, 2017.
• A. Durmus, É. Moulines, “Nonasymptotic convergence analysis for the unadjusted Langevin algorithm”, Annals of Applied Probability, 2017.

Research directions (ch. 9)

• J.-C. Hütter, P. Rigollet, “Minimax estimation of smooth optimal transport maps”, Annals of Statistics, 2021.
• A.-A. Pooladian, J. Niles-Weed, “Entropic estimation of optimal transport maps”, arXiv preprint, 2021.
• S. Chewi, J. Niles-Weed, P. Rigollet, Statistical Optimal Transport, lecture notes, arXiv preprint, 2024.
• B. Amos, L. Xu, J. Z. Kolter, “Input convex neural networks”, ICML, 2017.
• A. Makkuva, A. Taghvaei, S. Oh, J. Lee, “Optimal transport mapping via input convex neural networks”, ICML, 2020.
• A. Korotin, L. Li, A. Genevay, J. Solomon, A. Filippov, E. Burnaev, “Do neural optimal transport solvers work? A continuous Wasserstein-2 benchmark”, Advances in Neural Information Processing Systems, 2021.
• C. Léonard, “A survey of the Schrödinger problem and some of its connections with optimal transport”, Discrete and Continuous Dynamical Systems A, 2014.
• V. De Bortoli, J. Thornton, J. Heng, A. Doucet, “Diffusion Schrödinger bridge with applications to score-based generative modeling”, Advances in Neural Information Processing Systems, 2021.
• V. Khrulkov, G. Ryzhakov, A. Chertkov, I. Oseledets, “Understanding DDPM latent codes through optimal transport”, ICLR, 2023.
• L. Chizat, G. Peyré, B. Schmitzer, F.-X. Vialard, “Scaling algorithms for unbalanced optimal transport problems”, Mathematics of Computation, 2018.
• M. Liero, A. Mielke, G. Savaré, “Optimal entropy-transport problems and a new Hellinger–Kantorovich distance between positive measures”, Inventiones Mathematicae, 2018.
• F. Mémoli, “Gromov–Wasserstein distances and the metric approach to object matching”, Foundations of Computational Mathematics, 2011.
• G. Peyré, M. Cuturi, J. Solomon, “Gromov–Wasserstein averaging of kernel and distance matrices”, ICML, 2016.