A Machine Learning Framework for Solving High-Dimensional Mean Field Game and Mean Field Control Problems

• MFGs and MFCs model interactions within large agent populations across economics, finance, crowd motion, industrial engineering, data science, and material dynamics.Source: Introduction · p. 1
• Both problem classes couple a Hamilton-Jacobi-Bellman equation with a continuity equation, whose solution determines agents’ optimal actions without individual optimization.Source: Introduction · p. 1
• The framework targets potential MFGs and MFCs formulated as infinite-dimensional optimal control problems in density space.Source: Introduction · p. 1
• Grid-based numerical methods face a curse of dimensionality because computational complexity grows exponentially with spatial dimension.Source: Introduction · p. 1
• The method uses Lagrangian coordinates to solve the continuity equation and compute characteristics and Jacobian information from the potential.Source: Introduction · p. 1
• A specialized space-time neural network parameterizes the potential and enables direct penalization of HJB violations, transforming a broad class of MFGs into machine-learning problems.Source: Introduction · p. 1
• Experiments solve optimal transport and crowd-motion instances in up to 100 dimensions, validate against a two-dimensional Eulerian solver, and report improved accuracy and reduced cost from HJB penalization.Source: Introduction · p. 1
• The prototype implementation is released as open-source Julia software under a permissible license.Source: Introduction · p. 2

• Neural networks have been used to solve high-dimensional PDEs and control problems, while prior ML approaches addressed MFGs in spatial dimensions up to four.Source: Related work · p. 2, Mean Field Games and ML · p. 2
• High-dimensional optimal transport remains numerically challenging because Euclidean-coordinate methods require spatial discretization and are vulnerable to the curse of dimensionality.Source: Optimal Transport and ML · p. 2
• Normalizing flows provide a related ML approach that transforms samples from an unknown distribution toward an approximately normal distribution and can generate samples by reversing the flow.Source: Optimal Transport and ML · p. 2
• Earlier MFG methods were limited when formulations involved population density, whereas this framework estimates density along agents’ trajectories and supports nonlinear congestion terms.Source: Mean Field Games and ML · p. 2
• The neural architecture respects optimality-induced control structure and is used to approximately solve MFGs in 100 dimensions on a standard workstation.Source: Mean Field Games and ML · p. 2
• Prior Lagrangian approaches avoided some density-estimation difficulties through Fourier coordinates, while particle-in-cell density computation did not scale to high dimensions.Source: Lagrangian Methods in MFG · p. 2

• Mean field games model large populations of rational agents whose optimal behavior leads to a Nash equilibrium.Source: Mathematical Modeling · p. 2
• The population density evolves under agents’ strategies according to a continuity equation, with initial density ρ0.Source: Mean Field Games · p. 3
• An MFG equilibrium occurs when the forecast population distribution equals the actual distribution generated by optimal behavior.Source: Mean Field Games · p. 3
• Solving the coupled HJB and continuity PDEs yields optimal strategies for all agents simultaneously rather than optimizing each strategy separately.Source: Mean Field Games · p. 3
• Potential MFGs arise when the MFG system coincides with first-order optimality conditions of an infinite-dimensional constrained optimization problem.Source: Potential Mean Field Games · p. 3
• The paper re-parameterizes control through the potential and adds penalties for HJB violations without changing the minimizer.Source: Potential Mean Field Games · p. 4
• In MFC, a central planner selects one distributed strategy followed by all agents to minimize overall costs.Source: Mean Field Control · p. 4

• The method discretizes the optimal control problem by using Lagrangian coordinates to eliminate the continuity equation and other density-dependent terms.Source: Lagrangian Method · p. 4
• Characteristics are the trajectories of individual agents, making Lagrangian coordinates closely connected to the microscopic model and enabling stable mesh-free parallel computation.Source: Lagrangian Method · p. 4
• The framework focuses concretely on L2 transport costs while remaining adaptable to other transport-cost choices.Source: Lagrangian Method · p. 4
• Given a potential, characteristics are computed by solving an ODE with initial condition z(x,0)=x, and the continuity equation is evaluated along those curves.Source: Lagrangian Method · p. 4
• The method uses Jacobi’s identity to evolve the logarithm of the Jacobian determinant instead of directly computing the determinant, whose direct evaluation costs O(d^3) FLOPS.Source: Lagrangian Method · p. 4
• Running, terminal, and HJB-penalty contributions are accumulated along characteristics, with the HJB penalty enforcing equation violations during optimization.Source: Lagrangian Method · p. 5
• Quadrature converts the formulation into unconstrained optimization over the potential, and Monte Carlo integration is used for high-dimensional cases.Source: Optimization Problem · p. 5
• Independent sample trajectories make the computation trivially parallel, while stochastic approximation and sample average approximation provide optimization options.Source: Optimization Problem · p. 5

Neural-network parameterization

• A neural network parameterizes the potential, allowing HJB violations to be penalized during training and yielding a mesh-free Lagrangian scheme.The network is designed to approximate characteristics accurately and efficiently.Source: Machine Learning Framework for MFGs · p. 5, Neural Network Models for the Potential · p. 6
• The model maps s = (x, t) ∈ R^(d+1) to a scalar potential Φ(s), with trainable network parameters θ.Its architecture is modular, although architecture choices affect expressibility and training difficulty.Source: Neural Network Models for the Potential · p. 6

Derivative computation

• The residual network computes the potential gradient by back-propagation, with the first d components representing spatial derivatives and the final component representing the time derivative.These derivatives provide the inputs needed for characteristic computation.Source: Characteristics Computations · p. 6
• The Laplacian can be computed exactly through layerwise operations, with first-layer cost O(m · d) FLOPS.The computation squares and sums selected weight columns before one inner product.Source: Characteristics Computations · pp. 6–7

Computational complexity

• O(m^2 · d · M) FLOPS characterize each time-step computation, so cost is linear in input dimension and layers but quadratic in network width.This scaling motivates deep rather than wide architectures.Source: Characteristics Computations · pp. 6–7

Experimental design

• The experiments cover optimal transport and crowd motion, including dimensions d = 2, 10, 50, and 100, with a similar network architecture across cases.The study also evaluates scalability, the HJB penalty CHJB, optimization strategies, and validation against an Eulerian method.Source: Example 1: Dynamical Optimal Transport · p. 7, Numerical Experiments · p. 7, Example 2: Crowd Motion · p. 9

Optimal transport

• The two-dimensional optimal transport solution matches the initial and target densities, produces approximately straight characteristics, and partitions the target into eight approximately equal slices.The visualization uses push-forward and pull-back densities to assess the learned transformation.Source: Example 1: Dynamical Optimal Transport · p. 7, Figure 2 · p. 8
• Higher-dimensional optimal transport instances achieve quantitatively similar objective values for d = 10, 50, and 100, while runtime per iteration grows slower than predicted.The d = 50 case has a slight increase in terminal costs, while projected image quality remains similar across cases.Source: Example 1: Dynamical Optimal Transport · pp. 7–8

Training choices

• Removing the HJB penalty can prevent density matching and produce curved characteristics, whereas the penalty yields more accurate solutions at lower computational cost.Increasing characteristic time steps improves the penalty-free result but remains inferior and makes training four times more expensive.Source: Example 1: Dynamical Optimal Transport · p. 8
• BFGS with the SAA approach is more effective than ADAM with SA at reducing the HJB penalty at comparable computational cost, although both reach similar final results.The comparison is reported for the two-dimensional optimal transport instance.Source: Example 1: Dynamical Optimal Transport · p. 8

Validation

• Against a provably convergent Eulerian solver, the proposed method is competitive in d = 2 and scales to dimensions beyond the reach of Eulerian schemes.The comparison evaluates controls using an explicit Finite-Volume continuity-equation solver that was not used during optimization.Source: Example 1: Dynamical Optimal Transport · p. 8

Crowd motion

• The crowd motion experiment produces trajectories that avoid crowded regions and spread horizontally, with visually similar characteristics across dimensions and results comparable to a convergent Eulerian solver.For d = 2, the trajectories are approximately symmetric and bend around the central obstacle or congested region.Source: Example 2: Crowd Motion · pp. 8–9

• The framework combines Lagrangian PDE solvers with neural networks for mesh-free numerical solution of potential MFG and MFC problems.Source: Discussion and Outlook · p. 9
• The authors view the framework as a promising direction for large-scale MFG applications because it is mesh-free and suited to parallel computation.Source: Discussion and Outlook · p. 9
• The transformed MFG learning problems have unique traits that require further machine-learning attention.Source: Discussion and Outlook · p. 9
• Practical convergence results for neural networks in optimal control and firm theoretical results for regularization effects remain open mathematical issues.Source: Discussion and Outlook · p. 9
• More thorough parallel implementation may provide speed-up for more realistic MFG problems, while improved deep-learning optimization could increase efficiency.Source: Discussion and Outlook · p. 9
• Network architecture design lacks theoretical guidance and must support accurate evaluation of the gradients and Laplacians required by the Lagrangian scheme.Source: Discussion and Outlook · p. 9
• Applying the framework to more realistic problems is identified as an open application question, supported by the release of open-source code.Source: Discussion and Outlook · pp. 9–10

• The prototype extends Flux and is written in Julia rather than being fully described in the paper.Source: A. Supplemental Information for… · p. 12
• Unit tests are provided to verify the validity of the code.Source: A. Supplemental Information for… · p. 12
• Experiments ran on a shared Ubuntu workstation with 40 CPU cores and 1 TB of memory, so reported runtimes are rough averages.Source: A. Supplemental Information for… · p. 12

Dynamical Optimal Transport

• The default training setup uses SAA with resampling every 25 iterations and BFGS with a backtracked Armijo line search.Source: Dynamical Optimal Transport · p. 12
• The experiments use a Gaussian-mixture initial density and centered Gaussian target with covariance 0.3 · I across dimensions.Source: Dynamical Optimal Transport · p. 13
• Across dimensions, projected densities and characteristics are qualitatively similar, with trajectories moving toward the target along almost straight paths.Source: Dynamical Optimal Transport · p. 13
• Without the HJB penalty, the mean-field objective decreases substantially, but the learned densities and characteristics do not match the desired solution well.Source: Dynamical Optimal Transport · p. 13
• Using nt = 8 characteristic steps quadruples training cost but produces more meaningful results, likely because approximate transport costs become more accurate.Source: Dynamical Optimal Transport · p. 13
• BFGS converges in fewer iterations and reduces the HJB penalty more sharply early on, whereas ADAM has lower cost per iteration; overall costs are approximately comparable.Source: Dynamical Optimal Transport · p. 13

Crowd Motion

• The crowd-motion setup uses nt = 4 Runge-Kutta steps because characteristics are expected to bend around the obstacle.Source: Crowd Motion · p. 13
• The model uses a Kullback-Leibler terminal cost with λKL = 5 and HJB-penalty parameters α1 = 10 and α2 = 1.Source: Crowd Motion · p. 13
• Crowd-motion characteristics bend to avoid congestion and the domain center, effects represented by the running costs F.Source: Crowd Motion · p. 13
• Unlike optimal transport, where characteristics would be straight and parallel, the crowd-motion trajectories respond to obstacle and congestion costs.Source: Crowd Motion · pp. 13–14
• Across tested dimensions, the model learns similar dynamics, with agents avoiding the centrally placed high-travel-cost obstacle region.Source: Crowd Motion · p. 14

Eulerian reference

• The Eulerian reference solves a discretized convex optimization problem with staggered momentum and density variables on a regular space-time grid.Its implementation uses a conservative scheme and a primal-dual Newton method with multilevel initialization.Source: Comparison to Eulerian Methods · pp. 14–15

Method

• The comparison uses identical explicit finite-volume schemes to evaluate controls from the machine learning and Eulerian methods fairly.The trained neural network is evaluated on the Eulerian method’s space-time grid before solving the continuity equation.Source: Comparison to Eulerian Methods · p. 15

Eulerian comparison

• For the two-dimensional optimal-transport comparison, both methods are visually comparable, while the machine-learning objective is about 0.59% higher.The machine-learning result has higher image similarity, whereas the Eulerian method has lower transport costs.Source: Comparison to Eulerian Methods · p. 15, Table 2 · p. 17
• The Lagrangian method is less than 1% sub-optimal with the HJB penalty and two time steps, while scaling to dimensions beyond Eulerian schemes.This performance is obtained despite non-convex training and a vastly reduced parameter count.Source: Comparison to Eulerian Methods · p. 15

HJB penalty

• The HJB penalty improves density matching and characteristic geometry; without it, optimization can fail, while eight time steps remain inferior and cost four times more.With the penalty and two time steps, the pull-back and push-forward match their targets more closely and characteristics are nearly straight.Source: Figure 6 · p. 16

Optimization

• BFGS converges in fewer iterations than ADAM and produces a more substantial reduction of the HJB penalty.Figure 7 tracks the overall objective, mean field term, and HJB penalty using validation data.Source: Figure 7 · p. 17