Yves Achdou and Mathieu Laurière · "Mean Field Games and Applications: Numerical Aspects" · lecture notes posted on arXiv in March 2020 · my read of the math in them
My friend Johnathan sent me this one with a question: could a paper like this be made visual? It's a survey, written for a summer school, of how you actually compute a mean field game. A mean field game is what you get when a huge number of identical, selfish agents each optimize their own path through a crowd, and the crowd is nothing but the sum of their paths. The math comes out as two partial differential equations pointing in opposite directions in time, which is a nightmare for the usual "march forward from the initial condition" way of solving things. Most of the paper is about getting around that nightmare: a finite-difference scheme that respects the structure, a convex-optimization view that turns the two equations into one minimization, the linear-algebra machinery to make that fast, and Newton's method for the cases where the convex view doesn't apply. Then three applications with real simulations: two crowds sharing a hallway, a crowd evacuating a room, and an economy of households saving against a borrowing limit.
So yes, it can be made visual, and the way to do it is to run everything. Every figure below is a small made-up game solved for real, in your browser, by the methods the paper describes: the paper's upwind scheme, its fixed-point sweeps, Newton on its map $\Xi$, the Chambolle–Pock iteration on its discrete variational problem, multigrid with both coarsenings, its two-population hallway, its evacuation comparison, its Huggett economy. None of the numbers on the curves are the paper's; its own numbers appear only in the tables, where I say so. Where the paper's model has a hall of fifty meters and 3300 people, mine has a unit interval and a density. The shapes are the point.
Start with one agent. Its position $X_t$ lives on the circle $\mathbb{T} = \mathbb{R}/\mathbb{Z}$ (the paper uses the torus $\mathbb{T}^d$ to dodge boundary conditions; I'll stay in one dimension), it picks a velocity $v$, and noise shakes it:
$$dX_t = v(t, X_t)\, dt + \sqrt{2\nu}\, dW_t .$$Moving costs something, and so does being somewhere the crowd is dense. The agent pays a running cost $f(x, m, v)$ that depends on where it is, on the density $m$ of everyone else there, and on how hard it is pushing, plus a terminal cost $\phi(x, m)$ at time $T$. Here's the whole game in two sentences.
A mean field game is a pair $(\hat m, \hat v)$ of a flow of probability densities and a feedback control such that (1) $\hat v$ minimizes the expected cost $\mathbb{E}\!\left[\int_0^T f(X_t, \hat m(t, X_t), v(t, X_t))\, dt + \phi(X_T, \hat m(T, X_T))\right]$ among all controls $v$, with the crowd $\hat m$ treated as given, and (2) for every $t$, $\hat m(t, \cdot)$ is the law of $X_t$ when every agent uses $\hat v$.
Condition (1) is one person optimizing against a fixed crowd. Condition (2) says the crowd is made of people doing exactly that. A pair that satisfies both is a Nash equilibrium: nobody can do better by changing their own rule, because the crowd they'd be changing it against is already the crowd their rule produces. Before any of that can be solved, the crowd has to become a function. That's what "mean field" means, and it's the first picture.
The purple curve in that figure comes from the equation that moves a density when every particle in it follows the same drift. For a given velocity field $v$, the density $m^v$ of the law of $X_t$ satisfies the Kolmogorov–Fokker–Planck equation
$$\partial_t m - \nu \Delta m + \operatorname{div}(m\, v) = 0, \qquad m(0, \cdot) = m_0 .$$Diffusion spreads it, the divergence term carries it along $v$. Now the other half: how does one agent choose $v$? Through the Hamiltonian, which is a Legendre transform of the running cost.
With the drift $b(x, m, \gamma)$ and running cost $f(x, m, \gamma)$ of a control $\gamma$, the Lagrangian is $L(x, m, \gamma, p) = f(x, m, \gamma) + \langle b(x, m, \gamma), p \rangle$ and the Hamiltonian is
$$H(x, m, p) = \max_{\gamma \in \mathbb{R}^d} \big(-L(x, m, \gamma, p)\big).$$The paper assumes $H$ is $C^1$ in $(x, p)$ and strictly convex in $p$.
Suppose the drift is the control, $b = \gamma$, and the cost splits as $f(x, m, \gamma) = L_0(x, \gamma) + f_0(x, m)$: an effort term and a crowd term. Find $H$, and specialize to $L_0 = \tfrac{1}{2}|\gamma|^2$.
$H(x, m, p) = \max_\gamma \{-L_0(x, \gamma) - \langle \gamma, p \rangle\} - f_0(x, m) = H_0(x, p) - f_0(x, m)$, where $H_0$ is the Legendre transform of $L_0$. The crowd only enters through $f_0$, subtracted. With $L_0 = \tfrac{1}{2}|\gamma|^2$, the maximum of $-\tfrac{1}{2}|\gamma|^2 - \langle \gamma, p\rangle$ is at $\hat\gamma = -p$ and equals $\tfrac{1}{2}|p|^2$, so $H = \tfrac{1}{2}|p|^2 - f_0(x, m)$. The best velocity is minus the gradient of the value function: go downhill, at a speed equal to the slope.
Put the two halves together. The value function $u(t, x)$, the least cost an agent can still expect from position $x$ at time $t$, satisfies a Hamilton–Jacobi–Bellman equation, which runs backward from the terminal cost. The density satisfies the KFP equation with the optimal velocity $-H_p$, which runs forward from $m_0$.
The equilibrium best response is $\hat v(t, x) = \operatorname{argmax}_a \{-f(x, m, a) - \langle b(x, m, a), \nabla u \rangle\}$, where $(u, m)$ solves
$$\begin{cases} -\partial_t u - \nu \Delta u + H(x, m, \nabla u) = 0 & \text{in } [0, T) \times \mathbb{T}^d, \\[2pt] \partial_t m - \nu \Delta m - \operatorname{div}\big(m\, H_p(x, m, \nabla u)\big) = 0 & \text{in } (0, T] \times \mathbb{T}^d, \\[2pt] u(T, x) = \phi(x, m(T, x)), \qquad m(0, x) = m_0(x). \end{cases}$$Freeze the crowd $m$. Dynamic programming for the single agent says $u(t, x) = \min_v \mathbb{E}[\text{cost on } [t, t + dt] + u(t + dt, X_{t + dt})]$; expanding with Itô's formula gives $-\partial_t u - \nu \Delta u = \min_v \{f + \langle b, \nabla u\rangle\} = -\max_v\{-L\} = -H$, which is the first equation, and the minimizing $v$ is the stated $\hat v$. For the second, the optimal drift is $b(x, m, \hat v)$, and by the envelope theorem that equals $-H_p(x, m, \nabla u)$ (differentiate $H = -L(\hat\gamma)$ with respect to $p$; the $\partial \hat\gamma / \partial p$ term vanishes at the optimum). Put that drift into the KFP equation. Condition (2) of the definition is what lets the same $m$ appear in both. $\blacksquare$
If $m$ solves the KFP equation on the torus, then $\int_{\mathbb{T}^d} m(t, x)\, dx = 1$ for all $t$.
Integrate the equation over $\mathbb{T}^d$. The Laplacian integrates to zero, since $\int \Delta m = \int \operatorname{div}(\nabla m)$ and the divergence of any periodic field has zero integral (the flux out of one side comes back in the other). The same for $\operatorname{div}(m H_p)$. So $\frac{d}{dt}\int m = 0$, and $\int m_0 = 1$. $\blacksquare$
Now the structural trouble. The first equation needs $m$ on all of $[0, T]$ to be solved backward from $T$; the second needs $u$ on all of $[0, T]$ to be solved forward from $0$. Neither can go first. The simplest thing to try is to guess $m$, solve for $u$, solve for a new $m$, and repeat, which is what the paper's own sweeps do once the equations are on a grid. Here's that loop on my made-up game: a crowd starting at $x = 0.25$, drawn toward $x = 0.75$ by a cost $V(x) = 2\cos 2\pi(x - 0.25)$, charged $m$ for standing in a crowd, with $\nu = 0.1$ and $T = 1$.
The paper's scheme, from Achdou and Capuzzo-Dolcetta (2010), is a finite-difference method with one idea at its center: discretize the Hamiltonian so that the monotonicity of the continuous problem survives, and then the discrete KFP equation is forced on you as the adjoint of the discrete HJB equation. Here are the pieces, in one dimension with $N_h$ cells of width $h$ and $N_T$ steps of length $\Delta t$, writing $U^n_i \approx u(t_n, x_i)$ and $M^n_i \approx m(t_n, x_i)$.
The difference operators are $(D_t W)^n = (W^{n+1} - W^n)/\Delta t$, $(DW)_i = (W_{i+1} - W_i)/h$, $(\Delta_h W)_i = (W_{i+1} - 2W_i + W_{i-1})/h^2$, and the two-sided gradient $[\nabla_h W]_i = \big((DW)_i, (DW)_{i-1}\big)$: the forward and backward differences at $i$, kept separate. A discrete Hamiltonian $\tilde H(x, p_1, p_2)$ takes both, and must satisfy
Let $K = \mathbb{R}_- \times \mathbb{R}_+$ and $P_K$ the projection onto it. Then $\tilde H(p_1, p_2) = \tfrac{1}{2}|P_K(p_1, p_2)|^2 = \tfrac{1}{2}\big((p_1^-)^2 + (p_2^+)^2\big)$ satisfies ($\tilde H_1$) to ($\tilde H_4$) with $H_0(p) = \tfrac{1}{2}p^2$.
Write $a = p_1^- = \max(-p_1, 0)$ and $b = p_2^+ = \max(p_2, 0)$, so $\tilde H = \tfrac{1}{2}(a^2 + b^2)$. As $p_1$ increases, $a$ decreases or stays, so $\tilde H$ is nonincreasing in $p_1$; as $p_2$ increases, $b$ increases or stays, so $\tilde H$ is nondecreasing in $p_2$: ($\tilde H_1$). For ($\tilde H_2$), $p^- = 0$ or $p^+ = 0$ according to the sign of $p$, and $(p^-)^2 + (p^+)^2 = p^2$ either way. For ($\tilde H_3$), $\partial \tilde H / \partial p_1 = -p_1^-$ and $\partial \tilde H/\partial p_2 = p_2^+$, both continuous (zero on one side, linear on the other, meeting at $0$). For ($\tilde H_4$), $p \mapsto (p^-)^2$ is convex (its derivative $-2p^-$ is nondecreasing) and so is $p \mapsto (p^+)^2$, and a sum of convex functions of separate variables is convex. $\blacksquare$
"Upwind" because the forward difference only counts when it's negative (the function drops ahead of you) and the backward one only when it's positive (it rises behind you): each side listens to the direction the information flows from. The discrete system is then
$$-(D_t U_i)^n - \nu(\Delta_h U^n)_i + \tilde H(x_i, [\nabla_h U^n]_i) = f_0(x_i, M^{n+1}_i), \qquad U^{N_T}_i = \phi(M^{N_T}_i),$$ $$(D_t M_i)^n - \nu(\Delta_h M^{n+1})_i - \mathcal{T}_i(U^n, M^{n+1}) = 0, \qquad M^0_i = \bar m_0(x_i),$$with the discrete transport operator
$$\mathcal{T}_i(U, M) = \frac{1}{h}\Big( M_i \tilde H_{p_1}(x_i, [\nabla_h U]_i) - M_{i-1} \tilde H_{p_1}(x_{i-1}, [\nabla_h U]_{i-1}) + M_{i+1} \tilde H_{p_2}(x_{i+1}, [\nabla_h U]_{i+1}) - M_i \tilde H_{p_2}(x_i, [\nabla_h U]_i)\Big).$$Both equations are implicit: the HJB in $U^n$ given $U^{n+1}$, the KFP in $M^{n+1}$ given $M^n$. That transport operator looks like a mess until you see where it comes from.
For every grid function $W$, $\displaystyle\sum_i \mathcal{T}_i(U, M)\, W_i = -\sum_i M_i \big\langle \tilde H_p(x_i, [\nabla_h U]_i),\, [\nabla_h W]_i \big\rangle$. Consequently the operator $M \mapsto -\nu \Delta_h M - \mathcal{T}(U, M)$ is the adjoint of the linearization of $U \mapsto -\nu \Delta_h U + \tilde H(x, \nabla_h U)$, and the scheme conserves the total mass $h \sum_i M^n_i$.
Abbreviate $A_i = \tilde H_{p_1}(x_i, [\nabla_h U]_i)$ and $B_i = \tilde H_{p_2}(x_i, [\nabla_h U]_i)$. Then $h \sum_i \mathcal{T}_i W_i = \sum_i (M_i A_i - M_{i-1} A_{i-1}) W_i + \sum_i (M_{i+1} B_{i+1} - M_i B_i) W_i$. Shift the index in the subtracted terms (periodicity lets every sum run over all $i$): the first sum is $\sum_i M_i A_i (W_i - W_{i+1})$ and the second is $\sum_i M_i B_i (W_{i-1} - W_i)$. Dividing by $h$ gives $-\sum_i M_i \big(A_i (DW)_i + B_i (DW)_{i-1}\big)$, which is the claimed inner product with $[\nabla_h W]_i = ((DW)_i, (DW)_{i-1})$. The right side is exactly the Gateaux derivative of $\sum_i M_i \tilde H(x_i, [\nabla_h U]_i)$ in the direction $W$, with a minus sign, which is the adjoint statement. Taking $W \equiv 1$ makes every difference vanish, so $\sum_i \mathcal{T}_i(U, M) = 0$; summing the discrete KFP equation over $i$ (the Laplacian sums to zero by periodicity) gives $\sum_i M^{n+1}_i = \sum_i M^n_i$. $\blacksquare$
That is the design principle of the whole paper: the discrete KFP is not a separate approximation of the continuous KFP, it is the transpose of the discrete HJB's linearization. The continuous equations have the same relationship, and keeping it on the grid is what makes everything downstream (existence, positivity, the variational structure) work. Monotonicity ($\tilde H_1$) is what keeps the density nonnegative, and here is what happens without it.
Both equations are implicit in time, so there is no stability restriction tying $\Delta t$ to $h$. The paper points out that this is what lets the same scheme run over long horizons and adapt to the ergodic (infinite-horizon, stationary) problem, where there's no time to march in at all. The cost is a nonlinear solve per time step in the HJB, which is small: Newton on $N_h$ unknowns with a tridiagonal-plus-corners Jacobian, converging in three or four steps.
The convergence theorem the paper puts its weight on is the weak one, from Achdou and Porretta (2016), because it needs so little.
Assume $\nu > 0$, a continuous terminal cost $u_T(x)$, a bounded nonnegative $m_0$ of mass one, a continuous coupling cost $f_0$ bounded below, a continuous convex $C^1$ Hamiltonian $H_0$, and a discrete Hamiltonian with ($\tilde H_1$) to ($\tilde H_4$) plus two growth conditions. Then, along a subsequence of grids $h, \Delta t \to 0$, the piecewise-constant interpolations of the discrete solutions $(U, M)$ converge in $L^\beta$ to a weak solution $(\tilde u, \tilde m)$ of the MFG system. If weak solutions are unique (which Lasry–Lions monotonicity guarantees), the whole sequence converges.
Two things to notice. The theorem does not assume a solution exists; it produces one, so the scheme doubles as an existence proof. And it does not need the monotonicity of Lasry and Lions, only for uniqueness. What it does need is $\nu > 0$. The next figure is why that matters in practice: a monotone scheme keeps working as $\nu \to 0$, where the solution stops being smooth.
The paper's version of that picture is a two-dimensional ergodic game with $\nu = 0.001$, a potential with two wells, and the smoothing coupling $F[m] = (1 - \Delta)^{-2} m$: the value function comes out only Lipschitz and the density is two spikes at its minima. Mine is the one-dimensional cousin with the same coupling, solved down to $\nu = 0.005$ by Newton on the stationary system (with the kinks of $\tilde H$ smoothed by $\varepsilon = 0.01$, the regularization the paper mentions before it uses Newton). Same story: corners in $u$, spikes in $m$, and a scheme that doesn't blink.
Here's the move that I find most beautiful in the paper. When the cost splits as in the example above, with $f_0$ nondecreasing in $m$, the two coupled PDEs are the optimality conditions of a single convex minimization problem. The trick is to change variables from the velocity $\gamma$ to the momentum $w = m\gamma$, the amount of crowd times how fast it's going, and to rewrite the effort cost in terms of $(m, w)$.
The perspective of $L_0$ is $L(x, m, w) = m\, L_0\!\left(x, \dfrac{w}{m}\right)$ for $m > 0$, with $L(x, 0, 0) = 0$ and $L = +\infty$ otherwise. With $F(x, m) = \int_0^m f_0(x, s)\, ds$ and $\Phi(x, m) = \int_0^m \phi(x, s)\, ds$ (both $+\infty$ for $m < 0$), the variational problem is to minimize
$$\mathcal{J}(m, w) = \int_0^T\!\!\int_{\mathbb{T}^d} \big[L(x, m, w) + F(x, m)\big]\, dx\, dt + \int_{\mathbb{T}^d} \Phi(x, m(T, x))\, dx$$over pairs $(m, w)$ with $m \geq 0$ satisfying the continuity equation $\partial_t m - \nu \Delta m + \operatorname{div} w = 0$, $m(0) = m_0$.
Why that's a good idea: the constraint is now linear in $(m, w)$, and the cost is convex, so the whole thing is a convex program with a linear constraint, and the entire toolbox of convex optimization opens up. The convexity rests on one fact about the perspective.
If $L_0(x, \cdot)$ is convex, then $(m, w) \mapsto m L_0(x, w/m)$ is convex on $\{m > 0\}$.
Take $(m, w)$ and $(m', w')$ with $m, m' > 0$ and $\lambda \in [0, 1]$. Set $\bar m = \lambda m + (1 - \lambda) m'$ and $\bar w = \lambda w + (1 - \lambda) w'$. Then
$$\bar m\, L_0\!\left(\frac{\bar w}{\bar m}\right) = \bar m\, L_0\!\left(\theta \frac{w}{m} + (1 - \theta)\frac{w'}{m'}\right) \quad\text{with } \theta = \frac{\lambda m}{\bar m} \in [0, 1],$$because $\theta \cdot \frac{w}{m} + (1 - \theta)\frac{w'}{m'} = \frac{\lambda w + (1 - \lambda) w'}{\bar m}$. Convexity of $L_0$ gives $\bar m L_0(\bar w/\bar m) \leq \bar m\big[\theta L_0(w/m) + (1 - \theta) L_0(w'/m')\big] = \lambda m L_0(w/m) + (1 - \lambda) m' L_0(w'/m')$. $\blacksquare$
On the grid, the paper does the same thing with a discrete cost $\Theta(M, W)$ (the perspective of $\tilde H^*$, the convex conjugate of the discrete Hamiltonian, summed over the grid, plus $F$ and $\Phi$) and a linear map $\Sigma(M, W) = (AM + BW, M^0)$ where $A$ is the discrete heat operator and $B$ the discrete divergence, so that the constraint reads $\Sigma(M, W) = (0, \bar M^0)$. Then it proves the thing that justifies the whole section.
The discrete problem $\min \{\Theta(M, W) : \Sigma(M, W) = (0, \bar M^0)\}$ has a minimizer, and if $\nu > 0$ then $M^n_i > 0$ for $n \geq 1$. Its Fenchel–Rockafellar dual is a minimization over $U$, and the primal and dual solutions satisfy: $(U, M)$ solves the finite difference scheme above, and $W^n_i = -M^{n+1}_i P_K([\nabla_h U^n]_i)$. If $f_0$ and $\phi$ are strictly increasing, the solution is unique.
Take $\tilde H = \tfrac{1}{2}|P_K p|^2$, whose conjugate is $\tilde H^*(\gamma) = \tfrac{1}{2}|\gamma|^2$ on $K$ and $+\infty$ off it, so the node cost is $\tilde L(m, w) = |w|^2 / 2m$ when $-w/m \in K$. Form the Lagrangian with a multiplier $-U$ for the constraint: $\Theta(M, W) - \langle U, AM + BW \rangle - \langle \varphi, M^0 - \bar M^0\rangle$. Stationarity in $W^n_i$: $W^n_i / M^{n+1}_i = (B^* U)^n_i = -[\nabla_h U^n]_i$, clamped to the cone, which is $W^n_i = -M^{n+1}_i P_K([\nabla_h U^n]_i)$. Stationarity in $M^n_i$ for $n \geq 1$: $-|W^{n-1}_i|^2 / 2(M^n_i)^2 + F'(M^n_i) = (A^* U)^n_i$. The first term is $-\tfrac{1}{2}|P_K \nabla_h U^{n-1}|^2 = -\tilde H(\nabla_h U^{n-1})$, $F' = f_0$, and $(A^* U)^n_i = (U^{n-1}_i - U^n_i)/\Delta t - \nu(\Delta_h U^{n-1})_i$ by transposing the heat operator. Rearranged: $(U^{n-1}_i - U^n_i)/\Delta t - \nu (\Delta_h U^{n-1})_i + \tilde H([\nabla_h U^{n-1}]_i) = f_0(x_i, M^n_i)$, the discrete HJB at level $n - 1$. The constraint itself is the discrete KFP, since $BW$ with $W = -M P_K(\nabla_h U)$ is exactly $-\mathcal{T}(U, M)$. $\blacksquare$
The paper's statement writes $W = M P_K(\nabla_h U)$ and requires $w \in K$ in its definition of $\tilde L$; the sign has to be the one above for $BW$ to match the transport operator, and I list that among the slips at the end. The figure below checks it numerically.
The scheme was built by hand, equation by equation, with the adjoint structure as a guiding principle. The duality theorem says the same scheme is the KKT system of a convex program. So two proofs of existence (the monotone-scheme one and the convex-analysis one) land on the same discrete object, and every convergence theorem ever proved for ADMM or Chambolle–Pock now applies to a pair of PDEs that point in opposite directions in time.
With a convex program in hand you don't need the fixed-point sweeps at all. The paper describes two primal-dual methods: ADMM on the dual (the Benamou–Brenier idea, imported into mean field games by Benamou and Carlier) and the Chambolle–Pock algorithm (Briceño-Arias, Kalise and Silva). Both alternate two proximal steps. One of them is a projection onto the linear constraint, which means solving a linear system with the operator $\Sigma\Sigma^*$: fourth order in space when $\nu > 0$ (the heat operator applied twice), second order in time. The other is pointwise: at each grid node, a little convex problem in three variables. Here is Chambolle–Pock running on my game, and the duality theorem passing its test.
Chambolle–Pock with primal step $r$, dual step $s$ and extrapolation $\tau \in [0, 1]$ converges to a primal-dual solution when $rs < 1$. ADMM with penalty $r > 0$ converges (Eckstein and Bertsekas, 1992) when $\Lambda\Lambda^*$ is symmetric positive definite and a primal-dual solution exists, to a triple $(\bar U, \Lambda^* \bar U, \bar\sigma)$ solving the dual and primal problems.
In ADMM, the multiplier $\sigma^{(k+1)} = (M^{(k+1)}, W^{(k+1)})$ satisfies $M^{(k+1)} \geq 0$ at every iteration, and $W^{(k+1)}$ vanishes wherever $M^{(k+1)} = 0$.
The second step's optimality condition says $\sigma^{(k+1)} \in \partial \Psi(Q^{(k+1)})$ for the dual functional $\Psi$, with $\Psi(Q^{(k+1)}) < \infty$. For a convex function, $\sigma \in \partial \Psi(Q)$ means $\Psi^*(\sigma) = \langle \sigma, Q \rangle - \Psi(Q)$, which is finite. But $\Psi^*$ is the primal cost $\Theta$ in disguise, and $\Theta(M, W) = +\infty$ unless $M \geq 0$ and $W = 0$ where $M = 0$. So finiteness forces both. $\blacksquare$
Try the step sizes yourself. The condition $rs < 1$ is sharp, and inside it the balance between $r$ and $s$ sets the speed.
The expensive step in both primal-dual methods is that linear solve. When $\nu > 0$ the operator is fourth order in space, so its condition number grows like $h^{-4}$: halve the grid spacing and the system gets sixteen times harder for an iterative solver. The paper's answer is a multigrid preconditioner, and the section explaining it is a small, self-contained lecture on why multigrid works. The core observation is that simple iterations like Jacobi or Gauss–Seidel are terrible solvers but excellent smoothers.
For the discrete Laplacian $-\Delta_h$ on $N$ interior points of $(0, 1)$ with $h = 1/(N+1)$, weighted Jacobi with weight $\omega$ multiplies the $k$-th eigenmode $\sin(k\pi x_i)$ by $1 - \omega(1 - \cos k\pi h)$ at every sweep. With $\omega = 2/3$, every mode with $k \geq (N+1)/2$ is multiplied by a factor of absolute value at most $1/3$.
The matrix is $A = \frac{1}{h^2}\operatorname{tridiag}(-1, 2, -1)$, with eigenvectors $v_k = (\sin(k\pi x_i))_i$ and eigenvalues $\lambda_k = \frac{2}{h^2}(1 - \cos k\pi h)$ (plug $v_k$ in and use $\sin(a - b) + \sin(a + b) = 2\sin a \cos b$). Weighted Jacobi is $x \leftarrow x + \omega D^{-1}(f - Ax)$ with $D = \frac{2}{h^2} I$, so the error $e = x - x^*$ is multiplied by $I - \omega D^{-1} A$, whose eigenvalue on $v_k$ is $1 - \omega \frac{h^2}{2}\lambda_k = 1 - \omega(1 - \cos k\pi h)$. For $k \geq (N + 1)/2$ we have $\cos k\pi h \leq 0$, so the factor lies between $1 - 2\omega = -1/3$ and $1 - \omega = 1/3$. $\blacksquare$
So three sweeps leave at most $4\%$ of any fast mode and barely touch a slow one. What's left after smoothing is smooth, and a smooth error can be represented on a grid with half the points, where it looks twice as fast and the smoother bites again. Repeat down to a grid small enough to solve directly, then interpolate the corrections back up. That is a V-cycle.
Now the part specific to mean field games. The linear system lives on a space-time grid, and a first attempt coarsens it like any other grid: every other point in $x$, every other point in $t$. The paper reproduces a table from Achdou and Perez (2012) where this full coarsening preconditioner is used inside BiCGStab on a linearized MFG system, and it falls apart as the grid is refined. Those are the paper's numbers.
| $\nu$ \ grid | $32 \times 32 \times 32$ | $64 \times 64 \times 64$ | $128 \times 128 \times 64$ |
|---|---|---|---|
| 0.6 | 40 | 92 | fail |
| 0.36 | 24 | 61 | fail |
| 0.2 | 21 | 45 | fail |
Full-coarsening multigrid preconditioner, 4 levels: average number of preconditioned BiCGStab iterations to cut the residual by a factor of 0.01. From the paper's Table 1.
The reason is anisotropy. The operator is fourth order in $x$ and second order in $t$, so the unknowns are coupled far more tightly along $x$ than along $t$, and the smoother only smooths along $x$. The error stays rough in $t$, the coarse grid in $t$ can't see it, and the correction misses. The fix is semi-coarsening: halve the grid in $x$ only. Here are both, run on a small version of the operator.
With semi-coarsening, full weighting in $x$ (the $\tfrac{1}{4}, \tfrac{1}{2}, \tfrac{1}{4}$ stencil) and linear interpolation, the paper notes a classical rule: the orders of the two transfer operators should add up to at least the order of the differential operator, here $2 + 2 = 4$, and then a single smoothing sweep per level is already enough for convergence. The paper's own numbers, from the Chambolle–Pock runs of Briceño-Arias and coauthors on a two-dimensional game with $H_0(p) = |p|^{q'}/q'$, a quadratic crowd cost and a potential with three cosines, show how little the viscosity matters once the preconditioner is right.
| $\nu$ | residual cut by $10^{-3}$ | solved to $10^{-8}$ | ||||
|---|---|---|---|---|---|---|
| $32^3$ | $64^3$ | $128^3$ | $32^3$ | $64^3$ | $128^3$ | |
| 0.6 | 1.65 | 1.86 | 2.33 | 3.33 | 3.40 | 3.38 |
| 0.36 | 1.62 | 1.90 | 2.43 | 3.10 | 3.21 | 3.83 |
| 0.2 | 1.68 | 1.93 | 2.59 | 3.07 | 3.31 | 4.20 |
| 0.12 | 1.84 | 2.25 | 2.65 | 3.25 | 3.73 | 4.64 |
| 0.046 | 1.68 | 2.05 | 2.63 | 2.88 | 3.59 | 4.67 |
Average BiCGStab iterations per linear solve with the semi-coarsened multigrid preconditioner, two smoothing sweeps per level, $T = 1$. From the paper's Table 4.
| $\nu$ | grid $64 \times 64 \times 64$ | grid $128 \times 128 \times 128$ | ||||
|---|---|---|---|---|---|---|
| total (s) | first prox (s) | CP iterations | total (s) | first prox (s) | CP iterations | |
| 0.6 | 116.3 | 11.50 | 20 | 921.1 | 107.2 | 20 |
| 0.36 | 120.4 | 11.40 | 21 | 952.3 | 118.0 | 21 |
| 0.2 | 119.0 | 11.26 | 22 | 1028.8 | 127.6 | 22 |
| 0.12 | 129.1 | 14.11 | 22 | 1036.4 | 135.5 | 23 |
| 0.046 | 225.0 | 23.28 | 39 | 1982.2 | 260.0 | 42 |
Chambolle–Pock with the multigrid preconditioner: time to converge, time spent in the first (multigrid) proximal step, and iterations, stopping when successive densities differ by less than $10^{-6}$. From the paper's Table 3.
Read the two tables together. The number of inner BiCGStab iterations barely moves as the grid grows or $\nu$ shrinks: the preconditioner has made the linear algebra grid-independent. And the first proximal step, the one that uses multigrid, takes about a tenth of the total time; the pointwise second step, which doesn't, takes the rest. The outer iteration count does climb at the smallest $\nu$, but that is a property of Chambolle–Pock on a sharper problem, not of the preconditioner. The same long-horizon example also shows the turnpike that the paper mentions: a time-dependent game over a long horizon drives onto the stationary game and stays there until the end approaches.
The variational structure needs the Hamiltonian to split as $H_0(x, p) - f_0(x, m)$. Congestion models don't: if moving through a crowd is harder than moving through empty space, the effort cost depends on $m$ and $p$ together, there is no convex program, and the primal-dual methods are off the table. The paper's general-purpose tool for that case is Newton's method on the full discrete system, with two pieces of craft. The first is which system. Rather than Newton on the unknowns $(U, M)$ directly, it introduces the coupling values $f^n_i$ and $\Phi_i$ and defines a map $\Xi$: solve the discrete HJB with those as data, then the discrete KFP with the resulting $U$, then read off $g^n_i = f_0(x_i, M^{n+1}_i)$ and $\Psi_i = \phi(M^{N_T}_i)$. The equation to solve is $(I - \Xi)(f, \Phi) = 0$. Every evaluation of $\Xi$ includes a KFP solve, so every iterate's density is nonnegative with the right mass; Newton applied directly to $(U, M)$ would not guarantee that. The second piece is the initial guess: Newton needs a good one, and a solution at a slightly larger $\nu$ is one, so you walk $\nu$ down in steps, a continuation method.
The paper never builds the Jacobian of $\Xi$. It's dense and huge. Each Newton step solves its linear system with a nested iterative method such as BiCGStab, which only needs the product of the Jacobian with a vector, and that is a finite difference of $\Xi$ along the vector. I do the same below, with GMRES. One more thing the paper says in passing and that cost me an afternoon: it assumes the discrete Hamiltonians are $C^2$, "after a regularization procedure." The upwind Hamiltonian has kinks where a difference quotient changes sign, and Newton through a kink converges linearly at best. Smooth the kinks (I replace $p^-$ by $\tfrac{1}{2}(\sqrt{p^2 + \varepsilon^2} - p)$ with $\varepsilon = 0.01$) and quadratic convergence comes back.
The paper is honest that there is no general convergence proof for this strategy; "having these methods work efficiently is part of the know-how of the numerical analyst." It is, in its experience, much faster than the primal-dual methods whenever both apply, and it handles congestion, several populations, and couplings through the control, none of which the variational methods can touch. For a method that does come with a proof, it describes the recursive solver of Chassagneux, Crisan and Delarue: over a short enough time interval the forward-backward map is a contraction (Banach's fixed-point theorem gives existence and uniqueness constructively), so cut the horizon into $K$ short slices and solve them in the only order the coupling allows, which is backward, recursively, and repeatedly.
The first application is the one the paper calls state of the art, because it stacks every difficulty at once: congestion (so no variational structure), two populations (so no Lasry–Lions uniqueness), walls, entrances and exits (so real boundary conditions), and a stationary equilibrium with no ergodicity to lean on (people keep arriving and leaving). Two crowds share a hall. Each wants to reach certain doors, each finds it harder to move where the other crowd is dense, and each dislikes being the local minority. The model is four equations, two per population, and here are the paper's exact costs.
For populations $i \in \{0, 1\}$ and $j = 1 - i$, the congestion Hamiltonian and the coupling cost are
$$H_i(x, p; m_i, m_j) = \frac{|p|^2}{1 + m_i + 5 m_j}, \qquad \Phi_i(x, m_i, m_j) = 0.5 + 0.5\left(\frac{m_i}{m_i + m_j + \epsilon} - 0.5\right)_- + (m_i + m_j - 4)_+,$$and the system is, for each $i$,
$$\partial_t u_i + \nu \Delta u_i - H_i(\nabla u_i; m_i, m_j) = -\Phi_i, \qquad \partial_t m_i - \nu \Delta m_i - \operatorname{div}\!\Big(m_i \frac{\partial H_i}{\partial p}(\nabla u_i; m_i, m_j)\Big) = 0 .$$Read the constants. The $5$ in front of $m_j$ says the other crowd slows you five times as much as your own: people with different goals move in different directions, and crossing them is worse than following your own kind. The $0.5$ is the cost per unit time of being in the hall at all, the next term kicks in only when your population is in the local minority (the negative part), and the last term is pure crowd aversion, charged only above a density of four. None of the three terms is differentiable where it switches on, and the paper notes it regularizes them with a second small parameter before differentiating. The boundary conditions are the physical content.
Walls: $\partial_n u_i = 0$ and $\partial_n m_i = 0$ (the agent's process reflects, and with the first condition the zero-flux condition on $m_i$ reduces to the second). Exits for population $i$: $u_i = $ exit cost and $m_i = 0$ (an agent who reaches the door is out of the game; a negative cost is a reward). Entrances: $u_i = $ exit cost (set high so nobody leaves that way) and $\nu\, \partial_n m_i + m_i\, n \cdot \partial_p H_i = $ entry flux.
| door | population 0 | population 1 |
|---|---|---|
| north-west | entrance, flux 1, cost 0 | exit, cost 0 |
| south-west | exit, cost −8.5 | exit, cost −7 |
| north-east | exit, cost −4 | exit, cost −4 |
| south-east | exit, cost 0 | entrance, flux 1, cost 0 |
| south | exit, cost −7 | exit, cost −4 |
The paper's doors: a hall 5 long and 2.5 wide, the south door 1 wide, the others 0.33.
How do you compute a stationary equilibrium of a system with no stationarity built in? The paper's trick: solve the finite-horizon problem, then throw away the ends. Take the solution at the midpoint $T/2$, use $u_i(T/2, \cdot)$ as the new terminal condition and $m_i(T/2, \cdot)$ as the new initial condition, and solve again. Repeat, and the influence of the artificial initial and terminal data fades until $(u_i, m_i)$ stop depending on time. Within each finite-horizon solve, the four discrete equations are attacked the Newton way, but with one more layer: the two HJB equations don't involve $M$ at all once the coupling values $Y$ and the congestion factors $Z = (1 + M_i + 5 M_j)^{-1}$ are treated as data, so they can be marched backward, then the two KFP equations marched forward, and Newton is applied to the fixed point of that round trip, with continuation in $\nu$ from about $1$ down to the target. My toy below skips the time-dependent layer and iterates the stationary equations directly, which is cheaper and, for a toy, enough.
What the paper sees in its hall at $\nu = 0.30$: population 0 enters north-west and mostly leaves south-west; population 1 enters south-east and splits between the south and south-west doors, and the two overlap in the crossing. At $\nu = 0.16$ the picture changes qualitatively: each population keeps to its own region, and the paper reports the same down to $\nu = 10^{-3}$, with the transition zones between the crowds getting sharper as $\nu$ shrinks. My toy's lanes are the one-dimensional-wall version of that, and the overlap integral tells the same story in one number.
Everything so far is a Nash equilibrium: each agent optimizes against the others. Now ask a different question with the same ingredients. Suppose all agents follow one shared rule, chosen by a planner to minimize the total cost, and the planner knows that changing the rule changes the crowd. That is mean field type control (or control of McKean–Vlasov dynamics), and it is a genuinely different problem, because the density is now inside the thing being optimized rather than frozen while each agent chooses.
A mean field type control problem seeks a feedback $v^*$ minimizing $J(v) = \mathbb{E}\!\left[\int_0^T f(X^v_t, m^v(t, X^v_t), v(t, X^v_t))\, dt + \phi(X^v_T, m^v(T, X^v_T))\right]$, where $m^v$ is the density of $X^v$ itself, which solves $dX^v_t = b(X^v_t, m^v(t, X^v_t), v(t, X^v_t))\, dt + \sqrt{2\nu}\, dW_t$.
If a smooth optimal feedback $v^*$ exists, then $v^*(t, x) = \operatorname{argmax}_a\{-f(x, m, a) - \langle b(x, m, a), \nabla u\rangle\}$ where $(u, m)$ solves the MFG system with one extra term in the HJB equation,
$$-\partial_t u - \nu \Delta u + H(x, m, \nabla u) + m\, \frac{\partial H}{\partial m}(x, m, \nabla u) = 0,$$and the matching extra term $m\, \partial_m \phi$ in the terminal condition. (For local dependence on $m$, $\partial/\partial m$ is an ordinary partial derivative. The paper writes this term as an integral over the torus; its own example that follows uses the local form above.)
Write the cost as $J = \int_0^T \int \ell(x, m, v)\, m\, dx\, dt$ with $m$ constrained by the KFP equation $\partial_t m - \nu \Delta m + \operatorname{div}(mv) = 0$ (take $b = v$ and $\phi = 0$ to keep it short). Introduce a multiplier $-u$ and the Lagrangian $\int\!\!\int \ell m - \int\!\!\int u\,[\partial_t m - \nu \Delta m + \operatorname{div}(mv)]$. Vary $v$: $\int\!\!\int \delta v \cdot [m\, \ell_v + m \nabla u] = 0$ after integrating the divergence by parts, so $\ell_v = -\nabla u$, which is the optimality condition in the definition of $H(x, m, p) = \max_\gamma\{-\ell - \langle \gamma, p\rangle\}$ at $p = \nabla u$. Vary $m$, integrating by parts in $t$ and $x$: $\ell + m\, \ell_m + \partial_t u + \nu \Delta u + v \cdot \nabla u = 0$. At the optimal $v$, $\ell + v \cdot \nabla u = -H(x, m, \nabla u)$, and by the envelope theorem $\partial_m H = -\ell_m$ (differentiate $H = -\ell(\hat\gamma) - \langle \hat\gamma, p\rangle$ in $m$; the $\partial\hat\gamma/\partial m$ term vanishes at the maximum). Substituting: $-\partial_t u - \nu \Delta u + H + m\, \partial_m H = 0$. The extra term is the planner accounting for the cost its rule imposes on everyone else; a selfish agent, holding $m$ fixed, never sees it. $\blacksquare$
The paper's evacuation Hamiltonian is $H(x, m, p) = \dfrac{8|p|^2}{(1 + m)^{3/4}} - \dfrac{1}{3200}$ (with the sign that makes it convex in $p$; see the slips). Find the planner's Hamiltonian $H + m\, \partial_m H$.
$\partial_m H = 8|p|^2 \cdot \left(-\tfrac{3}{4}\right)(1 + m)^{-7/4}$, so $m\, \partial_m H = -6m|p|^2 (1 + m)^{-7/4}$. Then
$$H + m\, \partial_m H = |p|^2 \frac{8(1 + m) - 6m}{(1 + m)^{7/4}} - \frac{1}{3200} = |p|^2 \frac{2(1 + m) + 6}{(1 + m)^{7/4}} - \frac{1}{3200} = \left(\frac{2}{(1 + m)^{3/4}} + \frac{6}{(1 + m)^{7/4}}\right)|p|^2 - \frac{1}{3200},$$which is the coefficient in the paper's equation (71). The planner's agents behave as if motion were cheaper in a crowd than the selfish agents think ($2 + 6/(1 + m) < 8$ for $m > 0$), because slowing down near the exits so the people behind can move is worth it to the group. The KFP equation and the boundary conditions are unchanged: same physics, different rule.
The paper's experiment is a square hall fifty meters on a side with five rectangular obstacles, two exit doors, $\nu = 0.05$, a horizon of fifty minutes, no terminal cost, and 3300 people at density four per square meter at $t = 0$, solved by the Newton strategy above for both models. The result: the planner's crowd has lower density peaks at the corridor junctions and leaves the room faster. The difference between the two total costs is the price of anarchy, and it is positive. Here is my one-dimensional version, with the paper's congestion shape and a room of length one.
Every pedestrian in this model solves a partial differential equation before each step, and the equation already contains the future positions of everyone else. I have been in crowds. Nobody was solving anything. The model doesn't claim they were; it claims that a crowd nobody can improve on by themselves moves as if they had. That is somehow worse.
The last application is the oldest. Economists had mean field games before the name existed: the heterogeneous-agent models of Bewley, Huggett and Aiyagari, in which a continuum of households each save and consume against their own luck, and prices are whatever makes their choices add up. The paper's contribution, from Achdou, Han, Lasry, Lions and Moll, is the continuous-time version, and the part it dwells on is a boundary condition: households cannot borrow more than a fixed amount, which is a constraint on the state, and the way that constraint shows up in the equations is the whole story.
In the Huggett model a household has wealth $x$ and income $y_t \in \{y_1, y_2\}$, $y_1 < y_2$, a two-state Markov chain switching at rates $\lambda_1$ (from $y_1$ up) and $\lambda_2$ (from $y_2$ down). It chooses consumption $c_t \geq 0$ and its wealth evolves as $dx_t = (y_t + r x_t - c_t)\, dt$ at interest rate $r$, subject to the borrowing limit $x_t \geq \underline{x}$. It maximizes $\mathbb{E}\int_0^\infty e^{-\rho t} u(c_t)\, dt$ for a discount rate $\rho$ and a strictly increasing, strictly concave utility $u$, for instance the CRRA utility $u(c) = c^{1 - \gamma}/(1 - \gamma)$. The value functions $v_1(x), v_2(x)$ satisfy, for $x > \underline{x}$,
$$-\rho v_i + H(x, y_i, \partial_x v_i) + \lambda_i (v_j - v_i) = 0, \qquad H(x, y, p) = \max_{c \geq 0}\big((y + rx - c)p + u(c)\big).$$The Hamiltonian here is a maximum over consumption: $p$ is the marginal value of wealth, $(y + rx - c)p$ is what saving is worth, $u(c)$ what consuming is worth. For CRRA utility it has a closed form, and a shape that explains the boundary condition.
For $u(c) = c^{1 - \gamma}/(1 - \gamma)$ with $\gamma > 0$, $\gamma \neq 1$, the maximizer is $c^* = p^{-1/\gamma}$ and
$$H(x, y, p) = (y + rx)\, p + \frac{\gamma}{1 - \gamma}\, p^{1 - 1/\gamma}.$$Its minimum over $p$ is at $p_{\min} = (y + rx)^{-\gamma}$, where $c^* = y + rx$ and the drift $y + rx - c^*$ is zero. The drift is negative for $p < p_{\min}$ and positive for $p > p_{\min}$. With $H^\uparrow(p) = \max_{0 \leq c \leq y + rx}(\cdots)$ and $H^\downarrow(p) = \max_{c \geq \max(0, y + rx)}(\cdots)$, the nondecreasing and nonincreasing envelopes, $H = H^\uparrow + H^\downarrow - \min_p H$.
The objective $(y + rx - c)p + u(c)$ is concave in $c$ with derivative $-p + c^{-\gamma}$, zero at $c^* = p^{-1/\gamma}$. Substituting, $u(c^*) = p^{(1-\gamma)(-1/\gamma)}/(1 - \gamma) = p^{1 - 1/\gamma}/(1 - \gamma)$ and $-c^* p = -p^{1 - 1/\gamma}$, so $H = (y + rx)p + p^{1 - 1/\gamma}\big(\tfrac{1}{1 - \gamma} - 1\big) = (y + rx)p + \tfrac{\gamma}{1 - \gamma} p^{1 - 1/\gamma}$. Differentiating in $p$: $H'(p) = (y + rx) + \tfrac{\gamma}{1-\gamma}\big(1 - \tfrac{1}{\gamma}\big) p^{-1/\gamma} = (y + rx) - p^{-1/\gamma} = y + rx - c^*(p)$, the drift, which vanishes at $p^{-1/\gamma} = y + rx$, that is $p_{\min} = (y + rx)^{-\gamma}$, is negative below (large $c^*$) and positive above. So $H$ decreases then increases. $H^\uparrow$ restricts $c \leq y + rx$, which forces a nonnegative drift: for $p \geq p_{\min}$ the unconstrained maximizer already satisfies this and $H^\uparrow = H$; for $p < p_{\min}$ the constraint binds at $c = y + rx$ and $H^\uparrow = u(y + rx) = \min H$. Symmetrically $H^\downarrow = H$ for $p \leq p_{\min}$ and $= \min H$ above. Adding the two double-counts the constant $\min H$ once. $\blacksquare$
At the floor $x = \underline{x}$ a household may not dissave, so the HJB equation there uses $H^\uparrow$ in place of $H$. That is the state-constraint boundary condition, in the viscosity-solution sense of Soner and Capuzzo-Dolcetta–Lions, and it has a consequence on the other side of the problem. The long-run distribution of households has to allow for people stuck at the floor.
The ergodic measures $m_1, m_2$ on $[\underline{x}, \infty)$ are the long-run distributions of wealth among households with income $y_1$ and $y_2$: $\lim_{t \to \infty}\mathbb{E}[\varphi(x_t) \mid y_t = y_i] = \langle m_i, \varphi\rangle$. Each is written $m_i(A) = \int_A g_i(x)\, dx + \mu_i\, \delta_{\underline{x}}(A)$: a density $g_i$ plus a possible Dirac mass $\mu_i$ at the borrowing limit. In the interior they satisfy the stationary KFP equations $-\partial_x\big(g_i H_p(\cdot, y_i, \partial_x v_i)\big) + \lambda_i g_i - \lambda_j g_j = 0$, and at the floor a weak formulation that includes $\mu_i$ with $H^\uparrow_p$ in place of $H_p$.
$\displaystyle\int dm_1 = \frac{\lambda_2}{\lambda_1 + \lambda_2}$ and $\displaystyle\int dm_2 = \frac{\lambda_1}{\lambda_1 + \lambda_2}$.
In the weak formulation, take the test functions $\varphi_1 \equiv 1$ and $\varphi_2 \equiv 0$. Every term with a derivative of a test function disappears, and what remains is $-\lambda_1 \int dm_1 + \lambda_2 \int dm_2 = 0$ (the density terms and the Dirac terms together make up the full measures). With $\int dm_1 + \int dm_2 = 1$, solve the two linear equations. In words: in the long run, the flow of households switching from low income to high, $\lambda_1 \int dm_1$, equals the flow switching back. $\blacksquare$
The coupling between households is thin. There is no crowd cost in the HJB equation; what ties everyone together is the interest rate $r$, pinned down by one scalar equation, that the economy neither creates nor destroys wealth: $\int x\, dm_1 + \int x\, dm_2 = 0$. The paper's theorem then says what the equilibrium looks like.
Let $c^*_i(x)$ be the optimal consumption at income $y_i$, and suppose $r < \rho$. (1) The low-income saving policy $rx + y_1 - c^*_1(x)$ is negative for every $x > \underline{x}$. (2) Under a mild condition on $u$ at the floor, there is an explicit $\nu_1 > 0$ with $rx + y_1 - c^*_1(x) \sim -\sqrt{2\nu_1}\sqrt{x - \underline{x}}$ near $\underline{x}$, so low-income households hit the floor in finite time; then $\mu_1 > 0$ and $g_1(x) \sim \gamma_1/\sqrt{x - \underline{x}}$. (3) If the relative risk aversion is bounded, there is a target wealth $\bar x \geq \underline{x}$ with $rx + y_2 - c^*_2(x) > 0$ below it and $< 0$ above, and $\mu_2 = 0$.
So the only way out of the floor is a lucky draw of the high income, and a household with the high income saves up to a target and no further. The scheme that computes this is the paper's upwind scheme again, with one difference at the first node: the Hamiltonian there is $H^\uparrow$ of the forward difference alone. And the discrete KFP equation is obtained exactly as in the first section, by requiring $\langle D\mathcal{F}(V) W, M\rangle = 0$ for every $W$, the transpose of the linearized HJB operator. The paper shows that this discrete equation is consistent with the weak formulation if $h\, m_{0, j}$ is read as the Dirac mass $\mu_j$: the first cell's density, times its width, is the mass at the floor. In my toy, that cell is the red bar.
Aiyagari's model replaces the two-state income by a diffusing productivity $y_t$ and closes the loop through production instead of a bond market.
In the Aiyagari model wealth evolves as $dx_t = (w y_t + r x_t - c_t)\, dt$ and productivity as $dy_t = \mu(y_t)\, dt + \sigma(y_t)\, dW_t$, with the same borrowing limit. Firms produce with a Cobb–Douglas technology $F(X, Y) = A X^\alpha Y^{1 - \alpha}$ from the aggregate capital $X = \int x\, dm$ and labor $Y = \int y\, dm$, and the prices are the marginal products, $r = \partial_X F - \delta$ and $w = \partial_Y F$, with $\delta$ the depreciation rate.
In two dimensions the floor is a line, $\{x = \underline{x}\}$, and the singular part of the ergodic measure lives on it as a measure $\eta(dy)$ in the productivity variable. The paper points out that there is no PDE for $\eta$ and no local boundary condition for the density $g$, because the boundary term $\int H^\uparrow_p\, \partial_x \varphi\, d\eta$ in the weak formulation can't be written as a distribution acting on $\varphi(\underline{x}, \cdot)$ alone. The weak formulation is all there is, and the finite difference scheme discretizes it directly. The paper's figure shows the result: a ridge of mass along the floor at low productivity, exactly where the optimal drift points down into it.
The paper's own conclusion is modest: finite differences for the forward-backward system, a variational reading where one exists, and the algorithms to match; two applications argued in detail; and a list of what was left out (semi-Lagrangian schemes, gradient flows for monotone games, impulse control, and the high-dimensional problems where finite differences stop at about four dimensions and neural networks take over). Its closing sentence is that the most interesting open problems are in the applications, in the contact with social science and economics.
What I take from it is a single discipline applied five times. The discrete KFP equation is the transpose of the discrete HJB's linearization, and that one choice gives you mass conservation, positivity, a convex dual when there is one, the right boundary rows at a wall or a floor, and a Newton map whose every iterate is a genuine density. Then the linear algebra follows the structure of the operator: smooth in the strongly coupled direction, coarsen in that direction only. And then the applications are the same two equations wearing different costs: a crowd that pays to cross another crowd, a planner who charges each agent for the slowdown it causes, a household that cannot fall below a floor. Johnathan asked whether it could be made visual. Everything on this page was computed by the methods it describes, in a browser, in a few seconds. I think the answer is yes.
Slips I noticed in the source while building this, which I've left as they are in the paper: the discrete problem defines $\tilde L$ with $w \in K$ but needs $-w/m \in K$, and the duality theorem's $W = M P_K(\nabla_h U)$ is off by a sign from the transport operator; the mean field type control condition is written with an integral over the torus where its own example uses the local term $m\, \partial_m H$; the evacuation Hamiltonian is stated with the opposite sign from the HJB written under it; ADMM's second step is written with the multiplier from the step after it; and a handful of index ranges, typos and figure labels that don't change any mathematics.