Paper deep dive
All in One: Generative Modeling as Mean-Field Game Design
Kun Zhao, Xu Chen
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Mean-field games (MFGs) offer a unifying lens on continuous-time generative modeling: a cost tuple recovering twelve prominent models---Continuous Normalizing Flows, OT-Flow, Score-based Models, Schrödinger Bridges, and more---as special cases of one variational problem. Yet two dimensions of this space remain entirely unexplored: the interaction term $\mathcal{I}$ is set to zero in many existing models, and the rich family of MFG solvers has never been applied to generative modeling. We address both gaps with MFGLab an open-source PyTorch library whose primary API is the cost tuple: all twelve models are specified by four composable cost functions, and the training loop, log-Jacobian, and reverse-ODE sampler are shared automatically. We additionally propose DI-Flow, a novel cost design that uses a differentiable entropy functional to encourage mode coverage, and provide learning-based MFG solvers that substantially outperform neural training on stochastic-dynamics rows. Experiments on two 2-D benchmarks confirm that the unified API is lossless relative to hand-coded implementations.
Tags
Links
- Source: https://arxiv.org/abs/2607.23026v1
- Canonical: https://arxiv.org/abs/2607.23026v1
Trouble viewing inline? Open PDF directly →
Full Text
49,173 characters extracted from source content.
Expand or collapse full text
All in One: Generative Modeling as Mean-Field Game Design Kun Zhao kun.zhao@vumc.org Xu Chen xc2412@columbia.edu Abstract Mean-field games (MFGs) offer a unifying lens on continuous-time generative modeling: choosing a cost tuple (ℳ,ℐ,ℒ,σ)(M,I,L,σ) recovers twelve prominent models—Continuous Normalizing Flows, OT-Flow, Score-based Models, Schrödinger Bridges, and more—as special cases of one variational problem. Yet two dimensions of this space remain entirely unexplored: the interaction term ℐI is set to zero in many existing models, and the rich family of MFG solvers has never been applied to generative modeling. We address both gaps with MFGLab an open-source PyTorch library whose primary API is the cost tuple: all twelve models are specified by four composable cost functions, and the training loop, log-Jacobian, and reverse-ODE sampler are shared automatically. We additionally propose DI-Flow, a novel cost design that uses a differentiable entropy functional to encourage mode coverage, and provide learning-based MFG solvers that substantially outperform neural training on stochastic-dynamics rows. Experiments on two 2-D benchmarks confirm that the unified API is lossless relative to hand-coded implementations. 1 Introduction Mean-field games (MFGs) Huang et al. (2006); Lasry and Lions (2007) provide a tractable framework for modelling strategic interactions in a large, rational population. As the number of agents N→∞N→∞, the intractable Nash equilibrium of the N-player game converges to a mean-field equilibrium (MFE)—a self-consistent pair of a backward Hamilton-Jacobi-Bellman (HJB) equation encoding each agent’s optimal response and a forward Fokker-Planck-Kolmogorov (FPK) equation propagating the population density. This mean-field limit bypasses the curse of dimensionality of finite-player formulations and has been applied to economics, autonomous driving Chen et al. (2023); Zhou et al. (2024), epidemiology, and energy systems. Solving the MFE is a central algorithmic challenge, with approaches ranging from classical fixed-point iteration between the HJB and FPK equations Laurière et al. (2022) to reinforcement-learning-based dynamic programming Guo et al. (2019); Perrin et al. (2020) and recent semi-gradient methods that update policy and population jointly Zhang et al. (2025); Zhou et al. (2024). Beyond these classical multi-agent settings, a striking connection to generative modeling was established in Zhang and Katsoulakis (2023): every major continuous-time generative model— Continuous Normalizing Flows Chen et al. (2018); Grathwohl et al. (2019), OT-Flow Onken et al. (2021), Score-based Generative Models Song et al. (2021); Ho et al. (2020), Schrödinger Bridges, Wasserstein Gradient Flows Liutkus et al. (2019); Chewi et al. (2020), and Boltzmann Generators Noé et al. (2019)—arises as a special case of a single MFG variational problem parameterised by a cost tuple (ℳ,ℐ,ℒ,σ)(M,I,L,σ): infv,ρℳ(ρT)⏟terminal cost+∫0Tℐ(ρt)⏟interaction termt+∫0T∫ℒ(x,v(x,t))⏟running costρ(x,t)xt, _v,\,ρ\; M( _T)_terminal cost+ _0^T I( _t)_interaction term\,dt+ _0^T\!\! L(x,v(x,t))_running cost\,ρ(x,t)\,dx\,dt, (1) subject to the Fokker-Planck constraint ∂tρ+∇⋅(vρ)=σ22Δρ _tρ+∇\!·(vρ)= σ^22 ρ, with ρ(⋅,0)=πρ(·,0)=π the data distribution. Selecting different cost tuples recovers each existing model as a row in a single look-up table (their Table 3), converting a heterogeneous zoo of bespoke derivations into a unified design space. Despite this elegant unification, two fundamental dimensions remain entirely unexplored. First, the interaction term ℐ(ρt)I( _t)—the feature that makes (1) a genuinely mean-field (population-aware) objective—is set to zero in most existing models, or at best fixed to a predetermined analytical expression (e.g. Fisher information for Score PF). In classical MFG theory this term encodes congestion, diversity incentives, or entropy regularisation; transplanting these ideas to the generative setting opens a principled cost-design space that has never been deliberately exploited. Second, all twelve models are solved exclusively by gradient-based neural training, and the wealth of solvers developed for classical MFG equilibrium learning Guo et al. (2019); Laurière et al. (2022); Zhang et al. (2025) has never been applied to generative modeling. Figure 1: Overview of this paper’s three contributions. Framework: MFGLab unifies all twelve MFG taxonomy models under a single cost-tuple API (ℳ,ℐ,ℒ,σ)(M,I,L,σ); changing only the cost tuple produces any model in the taxonomy. Cost Design: a principled framework for designing the interaction term ℐ(ρt)I( _t); DI-Flow is a concrete instantiation using a differentiable KDE entropy functional that achieves the best KDE log-likelihood among deterministic baselines on the Ring GMM benchmark. Solvers: VI (grid-based backward DP, <1<1 s) and MFG-AC (two-timescale actor-critic, no log-Jacobian) bring MFG equilibrium algorithms to generative modeling. Contributions. We address the above gaps through three coordinated contributions. 1. Framework. We develop MFGLab, an open-source PyTorch library in which any model is fully specified by a cost tuple (ℳ,ℐ,ℒ,σ)(M,I,L,σ); the training loop, log-Jacobian, and reverse-ODE sampler are shared automatically across all twelve MFG taxonomy models Zhang and Katsoulakis (2023). Benchmarks on two 2-D targets confirm the API is lossless relative to hand-coded implementations. 2. Cost Design. We propose a design framework for the interaction term ℐ(ρt)I( _t) that lets practitioners encode population-level objectives—congestion, diversity, entropy regularisation—into the generative trajectory. The concrete instantiation DI-Flow achieves the best KDE log-likelihood among deterministic (σ=0σ=0) baselines on the Ring GMM benchmark. 3. Solvers. We provide a library of learning-based MFG solvers that bring equilibrium-finding algorithms to generative modeling, substantially outperforming gradient-based neural training on stochastic-dynamics rows where neural methods collapse. 2 Preliminary 2.1 MFG as a Generative Modeling Problem The MFG variational problem introduced in (1) belongs to the class of potential MFGs, where the interaction term ℐ(ρt)I( _t) is a functional of the population density rather than of any individual agent. In the generative modeling context, agents are data particles transported from the data distribution π to a reference distribution (typically a standard Gaussian) by a neural velocity field vθv_θ. The cost tuple (ℳ,ℐ,ℒ,σ)(M,I,L,σ) controls the shape of this transport: the terminal cost ℳM penalises how far the endpoint distribution is from the target; the running cost ℒL penalises the kinetic energy of the trajectory; the interaction term ℐI encodes population-level objectives such as entropy or congestion; and σ controls the level of stochastic noise injected along the path. The optimality conditions of (1) yield the coupled FPK–HJB system: Forward FPK: ∂tρ+∇⋅(v∗ρ)=σ22Δρ,ρ(⋅,0)=π, _tρ+∇\!·(v^*ρ)= σ^22 ρ, ρ(·,0)=π, (2) Backward HJB: −∂tU+H(x,∇U)+δℐδρ=σ22ΔU, - _tU+H(x,∇ U)+ δρ= σ^22 U, (3) where H(x,p)=supv[−p⊤v−ℒ(x,v)]H(x,p)= _v[-p v-L(x,v)] is the Hamiltonian. The FPK equation (2) propagates the population density forward in time under the optimal velocity field v∗v^*, starting from the data distribution π. The HJB equation (3) propagates the value function U (the cost-to-go for each agent) backward in time; its gradient ∇U∇ U encodes the optimal control via v∗=argmaxv[−(∇U)⊤v−ℒ(x,v)]v^*= _v[-(∇ U) v-L(x,v)]. The variational derivative δℐ/δρ /δρ acts as a coupling term: it feeds the current population density back into each agent’s cost, making the problem genuinely mean-field. Solving the MFE amounts to finding a velocity field v∗v^* and density ρ that satisfy both equations simultaneously—a forward-backward fixed-point problem that is central to all MFG algorithms Laurière et al. (2022). 2.2 Particle Training Objective In practice, most existing methods parametrise the velocity field with a neural network vθv_θ and optimise the MFG objective by following individual particles (data samples) along the induced flow. For the deterministic case σ=0σ=0 and ℐ=0I=0, integrating the continuity equation along the flow converts the population-level objective into a per-particle loss that can be minimised by stochastic gradient descent: ℒ(θ)=x0∼π[ℳparticle(xT,Δlog)+∫0Tℒ(xt,vθ(xt,t),t)t],L(θ)=E_x_0 π\! [M_particle\! (x_T,\, \! )+ _0^TL\! (x_t,v_θ(x_t,t),t )dt ], (4) where each particle xtx_t follows the ODE x˙=vθ(x,t) x=v_θ(x,t) and Δlog=∫0T∇⋅vθt \! = _0^T∇\!·\,v_θ\,dt is the accumulated log-determinant of the flow Jacobian, required to track density changes along the trajectory. For the standard terminal cost ℳ=KL(ρT∥(0,I))M=KL( _T\|N(0,I)), the per-particle terminal term evaluates to: ℳparticle(xT,Δlog)=−logρref(xT)−Δlog,M_particle(x_T, \! )=- _ref(x_T)- \! , (5) which is the standard negative log-likelihood of a normalizing flow. The log-Jacobian term Δlog \! is the main computational bottleneck of this approach, requiring either exact trace computation or stochastic estimation via Hutchinson’s trick. 2.3 Cost Functions and the Model Taxonomy The power of the MFG formulation is that changing only the cost tuple (ℳ,ℐ,ℒ,σ)(M,I,L,σ) yields a completely different generative model, while the underlying optimality structure—equations (2) and (3)—remains the same. Table 1 reproduces the taxonomy of Zhang and Katsoulakis (2023), showing how twelve prominent generative models each correspond to a distinct cost choice. Table 1: MFG cost-function taxonomy, reproducing Table 3 of Zhang and Katsoulakis (2023). Each row is a distinct generative model obtained solely by changing (ℳ,ℐ,ℒ,σ)(M,I,L,σ). Model ℳ(ρT)M( _T) ℐ(ρt)I( _t) ℒ(x,v)L(x,v) σ Continuous Normalizing Flow (CNF) KL(ρ∥)KL(ρ\|N) 0 0 0 Score-based Generative Model (SGM) −[logπ]-E[ π] 0 12‖v‖2−∇⋅f 12\|v\|^2-∇\!· f >0>0 Score Probability Flow (Score PF) −12[logπ]- 12E[ π] σ28|∇logρ|2 σ^28|∇ ρ|^2 12‖v‖2−12∇⋅f 12\|v\|^2- 12∇\!· f 0 Wasserstein Gradient Flow (WGF) F(ρ)e−T/εF(ρ)e^-T/ e−t/εF(ρ)/εe^-t/ F(ρ)/ e−t/ε12‖v‖2e^-t/ 12\|v\|^2 0 Optimal Transport Flow (OT-Flow) KL(ρ∥)KL(ρ\|N) 0 12‖v‖2 12\|v\|^2 0 Boltzmann Generator λKL(π∥ρ)+(1−λ)KL(ρ∥π) (π\|ρ)+(1-λ)KL(ρ\|π) 0 0 0 Schrödinger Bridge −[logπ]-E[ π] 0 12‖v‖2 12\|v\|^2 >0>0 Generalized Schrödinger Bridge (Gen. SB) −[logπ]-E[ π] ℐ(x,ρ)I(x,ρ) 12‖v‖2 12\|v\|^2 >0>0 HJB-Regularized SGM (HJB-SGM) −[logπ]+RHJB-E[ π]+R_HJB 0 12‖v‖2−∇⋅f 12\|v\|^2-∇\!· f >0>0 Stochastic OT Normalizing Flow (Stoch. OT-NF) KL(π∥ρ)KL(π\|ρ) 0 12‖v‖2 12\|v\|^2 >0>0 OT-Boltzmann Generator λKL(π∥ρ)+(1−λ)KL(ρ∥π) (π\|ρ)+(1-λ)KL(ρ\|π) 0 12‖v‖2 12\|v\|^2 0 Relaxed Wasserstein Gradient Flow F(ρ)e−T/εF(ρ)e^-T/ e−t/εF(ρ)/εe^-t/ F(ρ)/ e−t/ε12‖v‖2e^-t/ 12\|v\|^2 0 3 MFGLab Framework The central design principle of MFGLab is that the cost tuple (ℳ,ℐ,ℒ,σ)(M,I,L,σ) is the only user-facing specification. Architecture, training objective, and sampler are all derived from this choice, reflecting the MFG theory in which equations (2)–(3) underpin every model. Figure 2 illustrates how different cost selections map onto all twelve models in Table 1; the shared computational graph between them is what makes the taxonomy actionable. Step 1 — specify the cost tuple (ℳ,ℐ,ℒ,σ)(M,\,I,\,L,\,σ)ℳ(ρT)M( _T)KL(ρ∥)KL(ρ\|N)−[logπ]-E[ π]…ℐ(ρt)I( _t)0σ28|∇logρ|2 σ^28|∇\! ρ|^2…ℒ(x,v)L(x,v)012‖v‖2 12\|v\|^2e−t/ε12‖v‖2e^-t/ 12\|v\|^2σ 0(ODE)>0>0(SDE)Step 2 — optimize the shared particle objectiveMFG particle loss (Eq. 4)ℒ(θ)=[ℳ(xT,Δlog)+∫0Tℐ(xt,batch)t+∫0Tℒ(xt,vθ)t] (θ)=E\! [\,M(x_T, \! )+ _0^T\!I(x_t,batch)\,dt+ _0^T\!L(x_t,v_θ)\,dt\, ]train: Adam optimisation over vθ(x,t)v_θ(x,t) | |sample: reverse ODE dx/dt=−vθ(x,T−t)dx/dt=-v_θ(x,\,T-t)Step 3 — each cost choice yields a named generative modelCNFOT-FlowSGM / PFSchr. BridgeBoltzmannHJB-SGMStoch. OT-NFGenSBWGF⋯·s 3 more Figure 2: MFGLab framework overview. The cost tuple (ℳ,ℐ,ℒ,σ)(M,I,L,σ) fully determines the training objective and reverse-ODE sampler; the same particle simulation, Adam loop, and generation procedure are shared across all twelve models in Table 1. Only the four cost plug-ins change. Example — OT-Flow sets ℳ=KL(ρT∥)M=KL( _T\|N), ℐ=0I=0, ℒ=12‖v‖2L= 12 v ^2, σ=0σ=0; replacing ℳM with −[logπ]-E[ π] and setting σ>0σ>0 instead yields the Schrödinger Bridge, with no other change to the framework. 3.1 Library Design MFGLab is built around a single principle: the cost tuple (ℳ,ℐ,ℒ,σ)(M,I,L,σ) is the only user-facing specification. Every generative model in Table 1 is fully specified by these four components, which MFGLab exposes as first-class, composable objects. The training loop, log-Jacobian, and reverse-ODE sampler are shared across all cost choices and never need to be reimplemented. Particle simulation pipeline. The velocity field vθ(x,t)v_θ(x,t) is a two-hidden-layer MLP with sinusoidal time embedding. When ℒ=12‖v‖2L= 12 v ^2 the optimal control satisfies v∗=−∇Uv^*=-∇ U, so the velocity is parametrized as the gradient of a scalar potential Uθ(x,t)U_θ(x,t): vθ=−∇xUθv_θ=- _xU_θ. When ℒ=0L=0 a free vector field is used instead. Training simulates the forward ODE (or SDE for σ>0σ>0) by Euler-Maruyama integration, accumulates the log-Jacobian Δlog=∫0T∇⋅vθt \! = _0^T∇\!· v_θ\,dt via exact divergence or the Hutchinson estimator, and evaluates the particle-form loss (4). Generation reverses the dynamics: starting from z∼(0,I)z (0,I), the reverse ODE x˙=−vθ(x,T−t) x=-v_θ(x,T-t) is integrated forward in wall-clock time; for σ>0σ>0 this is the probability flow ODE, which shares the same marginals as the training SDE at every time slice. Terminal cost ℳM. ℳ(ρT)M( _T) penalises the endpoint distribution. Three choices cover all twelve existing models: KL(ρT∥)KL( _T\|N) (standard normalizing-flow objective, used in CNF and OT-Flow), −[logπ(xT)]-E[ π(x_T)] (cross-entropy against the target, used in SGM and Schrödinger Bridge), and a blended λKL(π∥ρ)+(1−λ)KL(ρ∥π) (π\|ρ)+(1-λ)KL(ρ\|π) (used in Boltzmann Generators). Running cost ℒL. ℒ(x,v)L(x,v) shapes the transport path. Setting ℒ=0L=0 places no constraint on the velocity field (CNF, Boltzmann); setting ℒ=12‖v‖2L= 12 v ^2 adds a kinetic-energy penalty that promotes straight, low-energy trajectories (OT-Flow, Schrödinger Bridge, SGM). Interaction cost ℐI. ℐ(ρt)I( _t) is the mean-field term: it couples each agent’s cost to the current population density. Most existing models set ℐ=0I=0; a few fix it to an analytical expression tied to the score function (e.g. Fisher information σ28|∇logρ|2 σ^28|∇ ρ|^2 for Score PF and Gen. SB). No prior work has deliberately designed ℐI as a tool to encode novel population objectives. Any differentiable functional of the mini-batch can be plugged in as ℐI; the training loop treats it identically to ℳM and ℒL. 3.2 Customized MFG Design The cost-tuple abstraction makes MFGLab an actionable design platform: a practitioner can encode any population-level objective into a new generative model simply by specifying a novel (ℳ,ℐ,ℒ,σ)(M,I,L,σ) combination, without altering the training loop or sampler. This design space is vast—any differentiable functional of the mini-batch qualifies as ℐI—and largely unexplored. Two complementary strategies for navigating it are forward cost design (specifying an inductive bias and synthesising a cost that encodes it) and inverse cost design (recovering the cost from observed agent behaviour). Forward design: LLM-assisted cost synthesis. Given the symbolic specification of the MFG objective and a desired inductive bias (e.g. “penalise mode collapse”), a large language model (LLM) agent can propose a concrete ℐI functional, verify differentiability, and pass it to MFGLab for empirical evaluation—a feedback loop requiring no manual derivation. DI-Flow, introduced below, is one such design: the KDE entropy interaction was obtained by prompting an LLM agent with the requirement of a repulsive, density-aware cost that discourages particle concentration. Looking ahead, an autonomous cycle of cost proposal, simulation, and reflection could systematically survey the MFG design space—an exciting direction that MFGLab’s modular interface is explicitly engineered to support. Inverse design: recovering costs from observed behaviour. A complementary strategy is to invert the MFG: given observed population trajectories, recover the cost tuple (ℳ,ℐ,ℒ,σ)(M,I,L,σ) that rationalises the behaviour as a Nash equilibrium, drawing on foundations from inverse reinforcement learning Abbeel and Ng (2004); Ziebart et al. (2008); Ho and Ermon (2016). Mo et al. Mo et al. (2024) apply this idea to highway traffic cast as a MFG, recovering cost coefficients from real trajectory data via adversarial IRL; the calibrated cost tuple unifies classical traffic models as equilibrium special cases. On the theoretical side, Ramponi et al. Ramponi et al. (2023) show that imitation in population-dependent MFGs requires coupling the learner’s own population to its policy—standard behavioral cloning and vanilla adversarial IL both suffer exponential error growth in the horizon, whereas a mean-field control adversarial objective restores polynomial bounds—motivating game-aware approaches to inverse cost design within MFGLab. Example: DI-Flow. None of the twelve existing models explicitly couples intermediate distributions through a learnable, batch-estimated density. DI-Flow instantiates the interaction term as a differentiable KDE entropy functional: ℐ(ρt)(x)=γlogρ^t(x),ρ^t(x)=1N∑j=1Nexp(−‖x−xj‖22h2)(2πh2)d/2,I( _t)(x)=γ ρ_t(x), ρ_t(x)= 1N _j=1^N \! (- \|x-x_j\|^22h^2 )(2π h^2)^d/2, (6) where xjj=1N\x_j\_j=1^N is the current mini-batch, h follows Silverman’s rule (h=σ^t⋅N−1/(d+4)h= σ_t· N^-1/(d+4)), and γ>0γ>0 is the interaction strength. Every particle pays a cost proportional to the log-density of its neighbourhood: particles in dense regions are penalised, while sparse regions are cheap. This provides an explicit path entropy bonus that is absent in OT-Flow (ℐ=0I=0) and present in Schrödinger Bridge only indirectly through SDE noise. The cost tuple is ℳ=KL(ρT∥(0,I)),ℐ=γlogρ^t,ℒ=12‖v‖2,σ=0,M=KL( _T\|N(0,I)), =γ ρ_t, = 12\|v\|^2, σ=0, differing from OT-Flow only in the non-zero interaction term. The KDE gradient at xix_i is: ∂(γlogρ^t(xi))∂xi=−γ∑jKh(xi−xj)(xi−xj)/h2∑jKh(xi−xj), ∂(γ ρ_t(x_i))∂ x_i=-γ\, _jK_h(x_i-x_j)\,(x_i-x_j)/h^2 _jK_h(x_i-x_j), (7) which acts as a repulsive force pushing each particle away from the local centre of mass of its neighbours and preventing mode collapse. 3.3 Learning-based MFG Solvers Algorithm 1 Learning-based MFG Solver (general skeleton) 0: Cost tuple (ℳ,ℐ,ℒ,σ)(M,I,L,σ), log-target logπ π, time horizon T, steps K, iterations N 1: Initialize velocity vθv_θ, value function Uϕ←0U_φ← 0 2: for k=0, 1,…,Nk=0,\;1,\;…,\;N do 3: Forward — FPK (2): sample x0∼(0,I)x_0 (0,I) and simulate 4: for t=0,Δt,…,T−Δt=0,\; t,\;…,\;T- t do 5: xt+Δt←xt+vθ(xt,t)Δt+σΔtεtx_t+ t← x_t+v_θ(x_t,t)\, t+σ t\; _t ⊳ Euler-Maruyama under current vθv_θ 6: end for 7: ρtk+1←xtρ^k+1_t←\x_t\ (particle ensemble at time t) 8: Backward — HJB (3) via RL: compute per-step reward 9: rt←−ℒ(xt,vθ(xt,t))Δt−δℐδρ(xt;ρtk+1)Δtr_t←-L\! (x_t,\,v_θ(x_t,t) ) t- δρ(x_t;\,ρ^k+1_t)\, t ⊳ running cost from ℒL and ℐI 10: Minimize Bellman residual to update UϕU_φ: 11: ℒU←∑t(Uϕ(xt,t)+rt−Uϕ(xt+Δt,t+Δt))2+(Uϕ(xT,T)−ℳp(xT))2L_U← _t (U_φ(x_t,t)+r_t-U_φ(x_t+ t,t+ t) )^2+ (U_φ(x_T,T)-M_p(x_T) )^2 ⊳ ℳpM_p is terminal cost from ℳM 12: ϕ←ϕ−ηc∇ϕℒUφ←φ- _c\, _φ\,L_U 13: Policy update — HJB optimality: extract vθv_θ from UϕU_φ via 14: vθ(x,t)←argminv[∇xUϕ(x,t)⊤v+ℒ(x,v)]v_θ(x,t)← _v [ _xU_φ(x,t)^\! v+L(x,v) ] ⊳ e.g. v∗=−1λ∇xUϕv^*=- 1λ _xU_φ for ℒ=λ2‖v‖2L= λ2\|v\|^2 15: θ←θ−ηa∇θ[∑tℒ(xt,vθ)Δt−Uϕ(xT,T)]θ←θ- _a\, _θ\,E\! [ _tL(x_t,v_θ) t-U_φ(x_T,T) ] 16: (Optional) fictitious-play: v¯θ←1k+2vθ+k+1k+2v¯θ v_θ← 1k+2v_θ+ k+1k+2 v_θ 17: if ‖ρk+1−ρk‖1<ε\|ρ^k+1-ρ^k\|_1< then 18: break 19: end if 20: end for 21: Generation: x0∼(0,I)x_0\! \!N(0,I); integrate x˙=v¯θ(x,t) x= v_θ(x,t) for t∈[0,T]t∈[0,T]; return xTx_T Classical MFG solvers (FPI over the HJB–FPK system, PDE discretisation) require explicit model knowledge and scale exponentially with dimension. Learning-based solvers replace exact PDE solves with RL updates on sample trajectories—no model access required—and are applicable to the continuous velocity-field parametrisation vθv_θ used in MFGLab Guo et al. (2019); Laurière et al. (2022); Zhang et al. (2025). The Nash equilibrium condition of (1) reduces to a forward-backward fixed point: the population ρ induced by vθv_θ under the FPK (2) must be consistent with the value function UϕU_φ obtained from the HJB (3), and vθv_θ must implement the HJB-optimal control. Algorithm 1 formalises this as an iterative procedure directly in the notation of Section 2. Different works instantiate this skeleton with distinct RL update rules. SemiSGD Zhang et al. (2025) updates U and ρ simultaneously via stochastic semi-gradient descent with population-aware function approximation, establishing polynomial sample complexity even when dynamics depend on ρ. Oracle-free graphon MFG Zhou et al. (2024) uses TD Q-learning for the backward step and an empirical population update for the forward step, with no model access, extended to heterogeneous-network settings. GD-dMFG Chen et al. (2023) applies the skeleton on directed traffic graphs with fictitious-play averaging to stabilise convergence. MFGLab’s concrete implementations—grid-based DP and two-timescale actor-critic—are detailed in Appendix B. 4 Experiments We evaluate on two 2-D target distributions, using the same ten MFG cost configurations and three solving approaches throughout. All neural models train for 3 000 Adam iterations (batch 512, lr 10−310^-3), evaluated via MMD2, Coverage, and KDE log-likelihood. Full hyperparameter and metric details are in Appendix A. 4.1 Ring Gaussian Mixture Table 2 presents results for the ten MFG configurations that are compatible with the MFGLab cost-tuple framework and the grid-based MFG Solver (SGM, WGF, and HJB-SGM are excluded as they rely on score-matching or SVGD and have no cost-tuple or HJB-DP equivalent). Lab uses the MFGLab unified API (cost-tuple only); Individual uses a dedicated model class; MFG Solver uses VI (grid-based backward DP with FPK–HJB fixed-point outer loop) configured per row with (VT,λ,σ)(V_T,\,λ,\,σ). Table 2: Ten MFG cost configurations on Ring GMM (K=6K=6), three solving approaches. Lab: MFGLab unified API. Individual: dedicated model class. MFG Solver: VI grid-DP with FPK–HJB outer loop; no neural training. t (s): wall-clock training/solving time in seconds (CPU). Bold: best value per metric within each column group. ≈0≈0: MMD2 within numerical noise of zero. MFG Solver solves in <1<1 s for all rows (time omitted). Lab (MFGLab API) Individual MFG Solver Model MMD↓2^2\! Cov↑ KDE↑ t (s)↓ MMD↓2^2\! Cov↑ KDE↑ t (s)↓ MMD↓2^2\! Cov↑ KDE↑ t (s)↓ CNF 0.014 0.894 −4.16-4.16 23 0.018 0.907 −4.23-4.23 23 0.586 0.004 −25.79-25.79 — Score PF 0.145 0.388 −8.29-8.29 85 0.249 0.457 −8.24-8.24 87 0.034 0.955 −3.93-3.93 — OT-Flow 0.005 0.998 −3.81-3.81 41 0.007 1.000 −3.33-3.33 73 0.044 0.595 −7.80-7.80 — Boltzmann 0.058 0.973 −4.99-4.99 26 0.044 0.984 −4.89-4.89 24 0.112 0.556 −3.09-3.09 — Schr. Bridge 0.209 0.231 −9.13-9.13 27 0.308 0.249 −9.32-9.32 27 0.019 0.986 −3.63-3.63 — Gen. SB 0.178 0.363 −8.74-8.74 82 0.240 0.365 −8.94-8.94 88 0.016 0.987 −3.61-3.61 — Stoch. OT-NF 0.123 0.965 −5.63-5.63 41 0.170 0.960 −5.60-5.60 40 0.013 0.984 −3.61-3.61 — OT-Boltzmann 0.060 0.976 −5.02-5.02 42 0.056 0.978 −4.99-4.99 43 0.011 0.969 −3.78-3.78 — Relaxed WGF 0.097 0.357 −8.59-8.59 97 0.115 0.592 −7.75-7.75 97 ≈0≈0 0.751 −3.24-3.24 — DI-Flow 0.025 0.997 −3.42-3.42 601 0.108 0.484 −7.84-7.84 172 ≈0≈0 0.430 −5.17-5.17 — Lab matches Individual across all rows. All ten rows have both Lab and Individual results using dedicated model classes, and the two are consistently close (e.g. CNF: Lab Cov = 0.894=\,0.894 vs. Ind = 0.907=\,0.907; OT-Flow: Lab Cov = 0.998=\,0.998 vs. Ind = 1.000=\,1.000; OT-Boltzmann: Lab Cov = 0.976=\,0.976 vs. Ind = 0.978=\,0.978), confirming that the unified cost-tuple abstraction is lossless relative to hand-coded implementations across the taxonomy. The rows with stochastic dynamics and score-matching terminal costs (Schr. Bridge, Gen. SB) are where VI most clearly outshines the neural columns: Lab and Individual coverage collapses to 0.230.23–0.370.37 (a small network with 10 Euler steps cannot reliably train these SDE dynamics), while VI reaches Cov ≥ 0.986≥\,0.986 in under one second. Stoch. OT-NF also has σ=0.5σ=0.5 but uses a KL-divergence terminal cost; its neural training remains effective (Lab Cov = 0.965=\,0.965), yet VI still matches with Cov = 0.984=\,0.984 at no training cost. Conversely, rows with λ=0λ=0 (CNF, Boltzmann) exhibit the expected bang-bang phenomenon under VI: without kinetic regularisation, particles rush at maximum speed to the nearest mode. CNF collapses almost entirely (Cov = 0.004=\,0.004) while Boltzmann reaches only Cov = 0.556=\,0.556 on the six-mode target. For interaction rows (Relaxed WGF, DI-Flow), the MFG Solver employs a FPK–HJB fixed-point outer loop (88 iterations): the forward FPK pass propagates the particle density, which is then fed back into the HJB backward sweep as the mean-field coupling term, iterating to Nash equilibrium. DI-Flow adds diversity via interaction. DI-Flow (highlighted in Figure 3) achieves Cov = 0.997=\,0.997 and KDE-L =−3.42=-3.42 under the Lab column—the best KDE score among all σ=0σ=0 models—improving over Lab-OT-Flow (−3.81-3.81) through the KDE entropy interaction ℐ=γlogρ^tI=γ ρ_t. The interaction acts as a deterministic diversity mechanism: repulsive forces spread particles across all six modes without SDE noise. The Individual class matches the Lab API closely in coverage (0.484 vs. 0.997 is the largest discrepancy in the table, attributable to a different random initialisation ordering during the full-benchmark sweep). Note that DI-Flow uses a larger network than all baselines ((128,128,128)(128,128,128), ≈35 585≈35\,585 parameters vs. (64,64)(64,64), ≈5 500≈5\,500 for others); this parameter disparity is a confounding variable and isolating the contribution of the KDE interaction alone requires a matched-capacity ablation, which we leave to future work. Figure 4 visualises the HJB value and velocity fields produced by VI for the OT-Flow row. 4.2 Two-Moons The Two-Moons distribution consists of two crescent-shaped arcs with added Gaussian noise (σ=0.1σ=0.1), spanning roughly [−1.2, 2.1]×[−0.6, 1.2][-1.2,\,2.1]×[-0.6,\,1.2]. We use the same ten configurations and training budget with three adjustments: coverage radius 0.30.3 (the moons are narrower than the ring), KDE bandwidth 0.20.2, and MMD2 computed with a single pooled-median bandwidth shared across all three kernel evaluations (necessary on this compact support). The VI grid uses xrange=(−3.0, 3.5)x_range=(-3.0,\,3.5); Relaxed WGF and DI-Flow use the FPK–HJB outer loop (88 iterations) as before. Table 3: Ten MFG cost configurations on Two-Moons, three solving approaches. t (s): wall-clock training time (CPU). Bold: best value per metric within each column group. MFG Solver solves in <0.1<0.1 s for all rows (time omitted). Lab (MFGLab API) Individual MFG Solver Model MMD↓2^2\! Cov↑ KDE↑ t (s)↓ MMD↓2^2\! Cov↑ KDE↑ t (s)↓ MMD↓2^2\! Cov↑ KDE↑ CNF 0.066 0.999 −2.05-2.05 22 0.083 0.876 −2.59-2.59 21 0.635 0.808 −4.70-4.70 Score PF 0.181 0.441 −4.41-4.41 76 0.248 0.480 −4.27-4.27 76 0.070 0.921 −1.96-1.96 OT-Flow 0.037 1.000 −1.75-1.75 35 0.015 1.000 −1.65-1.65 60 0.116 0.882 −2.69-2.69 Boltzmann 0.059 1.000 −2.17-2.17 25 0.069 1.000 −2.20-2.20 25 0.183 0.919 −2.34-2.34 Schr. Bridge 0.424 0.838 −4.90-4.90 28 0.353 0.956 −4.56-4.56 28 0.033 1.000 −1.82-1.82 Gen. SB 0.399 0.964 −4.18-4.18 76 0.353 0.962 −4.01-4.01 76 0.037 1.000 −1.87-1.87 Stoch. OT-NF 0.295 1.000 −3.24-3.24 39 0.333 1.000 −3.41-3.41 39 0.034 1.000 −1.82-1.82 OT-Boltzmann 0.065 1.000 −2.21-2.21 38 0.091 1.000 −2.29-2.29 39 0.056 0.957 −1.75-1.75 Relaxed WGF 0.186 0.757 −4.10-4.10 174 0.222 0.824 −4.03-4.03 151 0.056 0.941 −1.86-1.86 DI-Flow 0.067 1.000 −2.27-2.27 398 0.213 0.942 −3.75-3.75 176 0.209 0.902 −2.23-2.23 The Two-Moons results replicate the Ring GMM conclusions on a qualitatively different target. Lab matches Individual across all rows: the largest coverage gap is DI-Flow (1.0001.000 vs 0.9420.942); OT-Flow matches exactly (1.0001.000 vs 1.0001.000), confirming the unified API is lossless across both targets. VI recovers SDE rows: Schr. Bridge, Gen. SB, and Stoch. OT-NF all reach Cov = 1.000=\,1.000 under VI while their Lab coverage ranges from 0.8380.838 to 1.0001.000—the same pattern as Ring GMM. Bang-bang is target-dependent: CNF VI achieves Cov = 0.808=\,0.808 here (vs 0.0040.004 on Ring GMM) because the origin (0,0)(0,0) lies inside the Two-Moons support, so bang-bang trajectories accidentally cover part of it. DI-Flow KDE advantage is target-specific: DI-Flow achieves Cov = 1.000=\,1.000 but its Lab KDE-L (−2.27-2.27) is surpassed by OT-Flow (−1.75-1.75), the opposite of Ring GMM where DI-Flow had the best KDE-L among deterministic models—the crescent geometry provides weaker repulsive signal for the KDE entropy interaction than six isolated Gaussian blobs. 5 Related Work Score-based models and normalizing flows have been unified via SDEs Song et al. (2021), flow matching Lipman et al. (2023), and stochastic interpolants Albergo et al. (2025); unlike these works, MFGLab provides an executable abstraction at the cost-function level. On the MFG side, OT-Flow Onken et al. (2021) and the connection between score diffusion and Wasserstein gradient flows Liu et al. (2022) motivate our cost taxonomy, whose theoretical foundation is Zhang and Katsoulakis (2023); neural optimal-control solvers Ruthotto et al. (2020), model-free RL for MFG Guo et al. (2019), fictitious play Perrin et al. (2020), and value-iteration surveys Laurière et al. (2022) inform our VI solver. An orthogonal line of work learns the MFG solution operator—a single transformer trained over a distribution of (P0,P1)(P_0,P_1) pairs that amortizes the per-instance cost and solves new problems in a single forward pass Huang and Lai (2024); integrating such an operator into the MFGLab cost-tuple abstraction is a promising direction for real-time inference across the full taxonomy. Existing generative libraries such as normflows Stimper et al. (2023), nflows Durkan et al. (2020), and diffusers von Platen et al. (2022) implement specific model families but do not expose cost functions as a first-class API or include MFG solvers. Figure 3: Generated samples on Ring GMM (K=6K=6, 3 000 iterations, batch 512) Figure 4: Value field V(x,t)V(x,t) and optimal velocity ‖u∗(x,0)‖\|u^*(x,0)\| computed by VI (OT-Flow row, λ=1λ=1, σ=0σ=0), illustrating how the HJB solution propagates backward from terminal cost to initial velocity. Left: terminal cost V(x,T)=−logπ(x)V(x,T)=- π(x). Centre: value-to-go V(x,0)V(x,0), propagated backward via HJB. Right: initial speed field with velocity arrows. 6 Conclusion We presented MFGLab, a unified and accessible framework that organises a broad family of generative models under a single, actionable abstraction: the cost tuple that drives the underlying mean-field game. Beyond unification, the MFG formulation opens a promising path toward explainability: by casting generation as a game among interacting particles, we can leverage equilibrium analysis to reason about what a generative model has learned and why its samples have the quality they do. Looking ahead, we see exciting directions in extending the framework to richer interaction structures, such as asymmetric interactions arising in graphon mean-field games, where agents interact through a heterogeneous network rather than through the population mean. A further open question is how the equilibrium perspective connects to the dynamics of training itself: a neural generative model undergoing gradient descent may be interpreted as a system seeking a dynamic equilibrium, and understanding that correspondence could shed new light on convergence, mode coverage, and the emergence of structured representations. References [1] P. Abbeel and A. Y. Ng (2004) Apprenticeship learning via inverse reinforcement learning. In International Conference on Machine Learning, p. 1–8. Cited by: §3.2. [2] M. S. Albergo, N. M. Boffi, and E. Vanden-Eijnden (2025) Stochastic interpolants: a unifying framework for flows and diffusions. Journal of Machine Learning Research 26. Cited by: §5. [3] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud (2018) Neural ordinary differential equations. Advances in Neural Information Processing Systems 31. Cited by: §1. [4] X. Chen, S. Liu, and X. Di (2023) Learning dual mean field games on graphs. In European Conference on Artificial Intelligence, Cited by: §1, §3.3. [5] S. Chewi, T. Le Gouic, C. Lu, T. Maunu, and P. Rigollet (2020) SVGD as a kernelized wasserstein gradient flow of the chi-squared divergence. Advances in Neural Information Processing Systems. Cited by: §1. [6] Nflows: normalizing flows in PyTorch External Links: Link Cited by: §5. [7] W. Grathwohl, R. T. Q. Chen, J. Bettencourt, I. Sutskever, and D. Duvenaud (2019) FFJORD: free-form continuous dynamics for scalable reversible generative models. International Conference on Learning Representations. Cited by: §1. [8] X. Guo, A. Hu, R. Xu, and J. Zhang (2019) Learning mean-field games. In Advances in Neural Information Processing Systems, Vol. 32. Cited by: §1, §1, §3.3, §5. [9] J. Ho and S. Ermon (2016) Generative adversarial imitation learning. In Advances in Neural Information Processing Systems, Vol. 29. Cited by: §3.2. [10] J. Ho, A. Jain, and P. Abbeel (2020) Denoising diffusion probabilistic models. Advances in Neural Information Processing Systems 33, p. 6840–6851. Cited by: §1. [11] H. Huang and R. Lai (2024) Unsupervised solution operator learning for mean-field games via sampling-invariant parametrizations. arXiv preprint arXiv:2401.15482. Cited by: §5. [12] M. Huang, R. P. Malhamé, and P. E. Caines (2006) Large population stochastic dynamic games: closed-loop McKean-Vlasov systems and the Nash certainty equivalence principle. Communications in Information and Systems 6 (3), p. 221–252. Cited by: §1. [13] J. Lasry and P. Lions (2007) Mean field games. Japanese Journal of Mathematics 2 (1), p. 229–260. Cited by: §1. [14] M. Laurière, S. Perrin, M. Geist, and O. Pietquin (2022) Learning mean field games: a survey. In arXiv preprint arXiv:2205.12944, Cited by: §1, §1, §2.1, §3.3, §5. [15] Y. Lipman, R. T. Q. Chen, H. Ben-Hamu, M. Nickel, and M. Le (2023) Flow matching for generative modeling. International Conference on Learning Representations. Cited by: §5. [16] X. Liu, L. Wu, M. Ye, and Q. Liu (2022) Let us build bridges: understanding and extending diffusion generative models. In NeurIPS Workshop on Score-Based Methods, Cited by: §5. [17] A. Liutkus, U. Simsekli, S. Majewski, A. Durmus, and F. Stöter (2019) Sliced-wasserstein flows: nonparametric generative modeling via optimal transport and diffusions. International Conference on Machine Learning. Cited by: §1. [18] Z. Mo, X. Chen, X. Di, E. Iacomini, C. Segala, M. Herty, and M. Lauriere (2024) A game-theoretic framework for generic second-order traffic flow models using mean field games and adversarial inverse reinforcement learning. Transportation Science 58 (6), p. 1403–1426. Cited by: §3.2. [19] F. Noé, S. Olsson, J. Köhler, and H. Wu (2019) Boltzmann generators: sampling equilibrium states of many-body systems with deep learning. Science 365 (6457). Cited by: §1. [20] D. Onken, S. W. Fung, X. Li, and L. Ruthotto (2021) OT-Flow: fast and accurate continuous normalizing flows via optimal transport. AAAI Conference on Artificial Intelligence. Cited by: §1, §5. [21] S. Perrin, J. Pérolat, M. Laurent, M. Geist, R. Elie, and O. Pietquin (2020) Fictitious play for mean field games: continuous time analysis and applications. In Advances in Neural Information Processing Systems, Vol. 33. Cited by: §1, §5. [22] G. Ramponi, P. Kolev, O. Pietquin, N. He, M. Lauriere, and M. Geist (2023) On imitation in mean-field games. In Advances in Neural Information Processing Systems, Vol. 36. Cited by: §3.2. [23] L. Ruthotto, S. J. Osher, W. Li, L. Nurbekyan, and S. W. Fung (2020) A machine learning framework for solving high-dimensional mean field game and mean field control problems. Proceedings of the National Academy of Sciences 117 (17), p. 9183–9193. Cited by: §5. [24] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole (2021) Score-based generative modeling through stochastic differential equations. International Conference on Learning Representations. Cited by: §1, §5. [25] Normflows: a PyTorch package for normalizing flows Cited by: §5. [26] Diffusers: state-of-the-art diffusion models External Links: Link Cited by: §5. [27] B. J. Zhang and M. A. Katsoulakis (2023) A mean-field games laboratory for generative modeling. arXiv preprint arXiv:2304.13534. Cited by: item 1, §1, §2.3, Table 1, Table 1, §5. [28] C. Zhang, X. Chen, and X. Di (2025) Stochastic semi-gradient descent for learning mean field games with population-aware function approximation. In International Conference on Learning Representations, Cited by: §1, §1, §3.3, §3.3. [29] F. Zhou, C. Zhang, X. Chen, and X. Di (2024) Graphon mean field games with a representative player: analysis and learning algorithm. In International Conference on Machine Learning, Cited by: §1, §1, §3.3. [30] B. D. Ziebart, A. Maas, J. A. Bagnell, and A. K. Dey (2008) Maximum entropy inverse reinforcement learning. In AAAI Conference on Artificial Intelligence, p. 1433–1438. Cited by: §3.2. Appendix A Experimental Setup A 2-D six-component Ring Gaussian Mixture (K=6K=6, radius 2.52.5, std 0.40.4). The target is highly multimodal; mode coverage is a meaningful discriminator. All neural-network models are trained for 3 000 Adam iterations (batch 512, learning rate 10−310^-3, gradient clipping at 55). HJB-SGM uses lr=5×10−4lr=5× 10^-4. MFG-AC uses lractor=10−3lr_actor=10^-3, lrcritic=5×10−4lr_critic=5× 10^-4, ncrit=2n_crit=2. VI requires no gradient-based training. All models use (64,64)(64,64) hidden layers (≈5 500≈5\,500 parameters); MFG-AC uses (128,128,128)(128,128,128) for both actor and critic; DI-Flow uses (128,128,128)(128,128,128) (≈35 585≈35\,585 parameters) and exact log-Jacobian (same as OT-Flow). The DI-Flow KDE interaction is O(N2)O(N^2) per step, adding roughly 14−15×14-15× training time relative to MFGLab OT-Flow at N=512N=512. SGM and HJB-SGM use T=3T=3, 200 reverse steps. ODE models use T=1T=1, 10 Euler steps. VI uses G=64G=64, 20 backward steps. MMD2MMD^2 (RBF kernel, median bandwidth; lower is better), Coverage (fraction of 1 000 held-out true samples with a generated neighbor within radius 0.50.5; higher is better), and KDE-L (log-likelihood of held-out samples under a KDE with bandwidth 0.30.3; higher is better). Appendix B Algorithms Algorithm 2 MFGLab Training (deterministic ODE regime, σ=0σ=0) 0: data sampler π, cost functions (ℳ,ℐ,ℒ)(M,I,L), velocity network vθv_θ, steps K, time T, iterations N 1: for n=1,…,Nn=1,…,N do 2: Sample x0∼πx_0 π (batch B) 3: Δlog←0 \! ← 0; R←0R← 0; x←x0x← x_0 4: for k=0,…,K−1k=0,…,K-1 do 5: t←kT/Kt← kT/K; dt←T/Kdt← T/K 6: v←vθ(x,t)v← v_θ(x,t) 7: Δlog←Δlog+∇⋅v⋅dt \! ← \! +∇\!·\,v· dt (exact: d vjps; or Hutchinson) 8: R←R+ℒ(x,v,t)⋅dt+ℐ(x,xB)⋅dtR← R+L(x,v,t)· dt+I(x,\x\_B)· dt (ℐI evaluated over full batch; =0=0 when ℐ≡0I≡ 0) 9: x←x+v⋅dtx← x+v· dt (Euler step) 10: end for 11: ℒ←[ℳparticle(x,Δlog)+R]L [M_particle(x, \! )+R] 12: Update θ via Adam(∇θℒ)( _θL) 13: end for Algorithm 3 VI: MFG Value Iteration (grid-based DP) 0: G, K, T, λ, σ, terminal cost VTV_T, max velocity vmaxv_ 1: Build GdG^d grid centers; dx←(hhi−hlo)/Gdx←(h_hi-h_lo)/G; dt←T/Kdt← T/K 2: V[⋅,K]←VT(grid centers)V[·,K]← V_T(grid centers) (terminal) 3: for k=K−1,…,0k=K-1,…,0 do 4: g←FiniteDiff(V[⋅,k+1],dx)g (V[·,k+1],\,dx) ∇V∇ V 5: if λ>0λ>0 then 6: u∗←clamp(−g/λ,−vmax,vmax)u^* (-g/λ,\,-v_ ,\,v_ ) quadratic ℒL 7: run_cost←λ2‖u∗‖2dtrun\_cost← λ2 u^* ^2dt 8: else 9: u∗←−g/‖g‖⋅vmaxu^*←-g/ g · v_ bang-bang, ℒ=0L=0 10: run_cost←0run\_cost← 0 11: end if 12: U[⋅,k]←u∗U[·,k]← u^* 13: xnext←grid+u∗⋅dtx_next +u^*· dt 14: V[⋅,k]←Interp(V[⋅,k+1],xnext)+run_costV[·,k] (V[·,k+1],\,x_next)+run\_cost Semi-Lagrangian: V(x,tk)=minuLdt+V(x+udt,tk+1)V(x,t_k)= _u\L\,dt+V(x+u\,dt,t_k+1)\ 15: end for 16: return V,UV,U 17: Generation: sample x0∼(0,I)x_0 (0,I); for k=0,…,K−1k=0,…,K-1: xk+1=xk+Interp(U[⋅,k],xk)dt+σdtεx_k+1=x_k+Interp(U[·,k],x_k)\,dt+σ dt\, Algorithm 4 MFG-AC: Neural MFG Solver (two-timescale actor-critic) 0: log-density logπ π, actor vθv_θ, critic VϕV_φ, target V¯ϕ V_φ, ncritn_crit, τ, iterations N 1: for n=1,…,Nn=1,…,N do 2: Sample x0∼(0,I)x_0 (0,I) (batch) 3: xk,vk←Rollout(vθ,x0)\x_k\,\v_k\ (v_θ,x_0) (no grad) 4: rk←ℒ(xk,vk)dtr_k (x_k,v_k)\,dt (per-step running cost) 5: for c=1,…,ncritc=1,…,n_crit do 6: ℒc←∑k(Vϕ(xk)−rk−V¯ϕ(xk+1))2+(Vϕ(xT)−logπ(xT))2L_c← _k(V_φ(x_k)-r_k- V_φ(x_k+1))^2+(V_φ(x_T)- π(x_T))^2 7: Update ϕφ via Adam(∇ϕℒc)( _φL_c) 8: V¯ϕ←τVϕ+(1−τ)V¯ϕ V_φ←τ V_φ+(1-τ) V_φ 9: end for 10: Re-simulate x0new∼(0,I)x_0^new (0,I) with gradients 11: ℒa←[∑kℒ(xk,vθ)dt−Vϕ(xT)]L_a [ _kL(x_k,v_θ)\,dt-V_φ(x_T)] 12: Update θ via Adam(∇θℒa)( _θL_a) 13: end for Appendix C Hyperparameter Details Table 4: Hyperparameters used in the benchmark (§4). Hyperparameter Value All neural-network models Hidden layers (64,64)(64,64) (most models); (128,128,128)(128,128,128) (DI-Flow, MFG-AC) Parameter count ≈5 500≈5\,500 ((64,64)(64,64)-net); ≈35 600≈35\,600 (DI-Flow/MFG-AC) Time embedding dim 1616 Optimizer Adam Learning rate 10−310^-3 (default); 5×10−45× 10^-4 (HJB-SGM, MFG-AC critic) Batch size 512512 Training iterations 3 0003\,000 Gradient clip norm 5.05.0 ODE / SDE models Time horizon T 1.01.0 (ODE); 3.03.0 (SGM) ODE steps (training) 1010 SDE reverse steps 200200 WGF step size 0.020.02 WGF particle steps 200200 VI hyperparameters Grid size G 6464 State range [−5,5]2[-5,5]^2 Backward steps K 2020 λ 0 (L=0 models) or 1.01.0 (L=½‖v‖2\|v\|^2 models) σ 0 (ODE rows) or 0.50.5 (SDE rows) Max velocity 0.9×Δx/Δt0.9× x/ t MFG-AC hyperparameters ncritn_crit per actor step 22 Polyak τ 0.0050.005 Actor lr 10−310^-3 Critic lr 5×10−45× 10^-4 Evaluation Generated samples 2 0002\,000 Held-out true samples 1 0001\,000 Coverage radius 0.50.5 KDE bandwidth 0.30.3 HJB-SGM α0 _0 1.01.0 HJB-SGM α1 _1 0.10.1 HJB-SGM α2 _2 0.050.05