key: cord-0126633-n5uakzhz authors: Chow, Yat Tin; Fung, Samy Wu; Liu, Siting; Nurbekyan, Levon; Osher, Stanley title: A numerical algorithm for inverse problem from partial boundary measurement arising from mean field game problem date: 2022-04-11 journal: nan DOI: nan sha: d939a969c57930e67d2a9b50b79ef5c6e936d7a3 doc_id: 126633 cord_uid: n5uakzhz In this work, we consider a novel inverse problem in mean-field games (MFG). We aim to recover the MFG model parameters that govern the underlying interactions among the population based on a limited set of noisy partial observations of the population dynamics under the limited aperture. Due to its severe ill-posedness, obtaining a good quality reconstruction is very difficult. Nonetheless, it is vital to recover the model parameters stably and efficiently in order to uncover the underlying causes for population dynamics for practical needs. Our work focuses on the simultaneous recovery of running cost and interaction energy in the MFG equations from a emph{finite number of boundary measurements} of population profile and boundary movement. To achieve this goal, we formalize the inverse problem as a constrained optimization problem of a least squares residual functional under suitable norms. We then develop a fast and robust operator splitting algorithm to solve the optimization using techniques including harmonic extensions, three-operator splitting scheme, and primal-dual hybrid gradient method. Numerical experiments illustrate the effectiveness and robustness of the algorithm. The basis for the MFG framework is the concept of Nash equilibrium, where agents cannot unilaterally improve their objectives. Under suitable regularity assumptions, a common MFG model reduces to the following system of partial differential equations (PDE): −∂ t φ(x, t) − ν∆φ(x, t) + H(x, ∇ x φ(x, t)) = F (x, ρ(·, t)), in Ω × (0, 1), ∂ t ρ(x, t) − ν∆ρ(x, t) − ∇ x · (ρ(x, t) ∇ p H(x, ∇ x φ(x, t))) = 0, in Ω × (0, 1), ρ(x, 0) = ρ 0 (x), φ(x, 1) = g(x), in Ω . Introduced in [19, 17, 18] and [14, 13] , MFG is an actively growing field significantly advancing the understanding of social cooperation and economics [1, 5, 10] , biological systems [30] , election dynamics [31] , population games [20] , robotic control [25] , machine learning [28, 22] , dynamics of multiple populations [7] . Recently, they are utilized to understand pandemic modeling and control such as COVID-19 [21] . With the significant descriptive power of MFGs, it is vital to consider inverse problems arising in MFGs. We aim to reconstruct MFG parameters for a class of nonlocal problems, including the geometry of the underlying space and the interactions between large crowds, based on partial population observations. More specifically, we are interested in the following problem. for finitely many examples of ρ 0 and terminal cost g, can we numerically recover the speed field κ(x) and the interaction kernel K(x, y) from observations? Such a model-recovery algorithm can help understand the underlying population dynamics in numerous problems, such as migration flow or contagious rate of COVID-19. We further envision applications to include rescue and exploration team management, policymaking, diplomacy, election modeling, catastrophe management, and evacuation planning. Note that m(x, s) = −ρ (x, s) ∇ p H(x, ∇ x φ(x, s)) represents the flux of the agents through the state x at time s as a result of optimal actions. The interpretation of the flux is straightforward for crowd-dynamics models and can be measured by counting people crossing checkpoints or parts of the border. For such models, the value function φ(x, s) could represent the travel cost for a traveller who is at location x at time s. Hence, one could also consider an inverse problem where one observes the value function, instead of the flux, by looking at travel companies' expenses or consumer ticket prices (discounted for the companies' profit margins). For economic and finance models [11, 1, 5] the state variables typically represent asset (wealth, income, inventory) levels instead of a physical location. Hence, the value function represents maximal utility for agents with a given asset level, and the flux represents the total amount of transactions performed by them. Interestingly, in economic models one often has implicit mean-field interactions reflected in market-clearing type conditions instead of an explicit interaction functional F (x, ρ). Hence, a related inverse problem is to find an appropriate market-clearing condition or tune its parameters for a given economy. This manuscript addresses explicit models with flux data leaving the implicit ones with other data types for future work. 1.1. Our contributions. We propose a new MFG inverse problem with non-invasive partial boundary measurements. Based on insights from [26, 23, 24] we postulate a feature expansion representation for the interaction kernel K and formulate the forward problem (1.1) as a convex-concave saddle point problem. Furthermore, we design a three-operator splitting scheme [8] for the resulting inverse problem with a saddle-point constraint. The algorithm reduces to a forward-backward splitting for the parameter updates, a primaldual hybrid gradient for the forward saddle point problem update, and a proximal-point algorithm for the adjoint problem update. Intriguingly, our algorithm applies to inverse problems whose forward problem has a saddle point structure beyond MFG. 1.2. Related work. Despite of the large body of work on theory, numerical methods, and applications [2] , inverse problems arisen from MFG is still quite an unexplored terrain. To the best of our knowledge, only [9, 16, 4 ] study such problems. The work in [9] is the closest to our objective but considers the case with a full space-time measurement of data in the sampling domain. However, most inverse problems in practice only have partial boundary measurements available, either obtained via non-invasive measurement methods or because of the limited access to the sampling domain. Compared with the case with full space-time measurement in the domain, inverse problems with only partial boundary measurements are generally known to be more severely ill-posed. In this work, we focus on the recovery problem with only boundary measurements coming from several measurement events. The rest of the paper is organized as follows. In Section 2, we introduce an abstract inverse problem with a saddle point constraint and a generic algorithm to solve it. In Section 3, we present the inverse MFG formulation. Next, in Section 4 we discuss the implementation of the algorithm in Section 2 for the inverse MFG in Section 3. Section 5 contains three numerical examples to demonstrate the robustness and effectiveness of our algorithm. Finally, Section 6 contains a discussion and concluding remarks. 2. An inverse problem with a saddle point forward model In this section, we formulate an abstract inverse problem with a saddle point forward model. We discuss suitable Karush-Kuhn-Tucker (KKT) conditions and a generic algorithm to solve such inverse problems. A forward saddle point problem. Consider a saddle point problem Here, x is the primal variable, and y is the dual variable in the forward problem. Commonly, y is either used to handle constraints in the forward problem or linearize nonlinear components via some splitting scheme. The variable c represents model parameters associated with the functional F , while u represents boundary and initialterminal conditions. Given model parameters c, we define a boundary measurement map Λ c as follows: where Π B,x , Π B,y denote a projection operator that represent the partial boundary measurements of x, y available. We note that u corresponds to boundary conditions of the forward problem, whereas B is the subset of the domain where the partial measurements are collected. The inverse problem and a generic algorithm. Assume that are noisy measurements for a given Our goal is to recover c ∈ D. We formulate this problem as a constrained optimization problem where R is a suitable regularizer and · are suitable choices of (semi)-norms. Introducing Lagrange multipliers (dual variables) (λ x i , λ y i ), (2.2) reduces to Thus, the KKT condition for (2.2), (2.3) are as follows: for i = 1, . . . , N . Here, Π * B,x , Π * B,y are the adjoints of Π B,x , Π B,y , respectively. Finally, we formulate these KKT conditions as an inclusion problem Note that A is monotone but C is not known to be monotone in general. Here, we outline an iterative algorithm for solving (2.3) . At (n+1)-th iteration, we first update the adjoint variables {(λ x i , λ y i )} N i=1 using the Chambolle-Pock method [6] ; then we update c for the inverse problem by taking a proximal gradient step; next we use the Chambolle-Pock method again to compute forward problems Summarizing, a high level description of the (n + 1)-th iteration is as follows: In what follows, we specify the MFG inverse problem and the implementation of the algorithm above for it. Here, we explain the saddle point problem formulation of nonlocal MFG [26, 23, 24] and formulate the inverse MFG problem of our interest. 3.1. Saddle point formulation of MFG via feature-space expansions. Consider the following MFG system with nonlocal couplings: in Ω . We assume that K is positive definite and translation invariant, which yields that the mean-field interaction satisfies the Lasry-Lions monotonicity condition [19] and agents are crowd averse. Moreover, (3.1) admits a saddle point formulation Here, χ Z (z) is the indicator function over the set Z defined by Modeling the interaction term 1 2 Ω ×Ω K(x, y)ρ(x, t)ρ(y, t)dxdy directly is costly for both forward model and the inverse problem. Moreover, based on the works from [26, 23, 24, 3] , we model and approximate this term using feature-space expansions. More specifically, based on Bochner's theorem [27] , we postulate that we obtain K(x, y) = ζ(x; µ, ω) · ζ(y; µ, ω), µ, ω ∈ C odd=even . Using this representation, we obtain where a(t) = (a 1,1 (t) , a 1,2 (t) , .., a r,1 (t) , a r,2 (t)) are auxiliary dual variables. The last equality is a result from [26] . Hence, For more details on representation of nonlocal MFG interactions via a basis and computational methods, see [26, 24, 23, 3] . We also attach an example Algorithm 3 for calculating the nonlocal mean-field game problem in the appendix. An inverse mean-field game problem. Denoting by u = (ρ 0 , g), x = (ρ, m), y = (φ, a), c = (κ, µ), we place the MFG forward model in the abstract framework (2.1). Next,we assume that Ω ⊂ Ω and κ(x) is known in the domain Ω \Ω. We refer to Ω and Ω as sampling and computational domains, respectively. An example is shown in Figure 1 , where the Ω is the large square domain, while Ω is the inner square with its boundary highlighted in red. Denote ρ as the solution to the mean-field game system. From left to right, the pictures display the density distribution ρ at time t = 0.1, 0.5, 0.9. The solid red line represents the boundary of domain Ω. In this mean-filed game, the density travels from the right towards the left, crossing the boundary ∂Ω twice. Next, we take for the partial boundary measurement along the boundary ∂Ω. Here, ∂Ω + means that the normal vector n is pointing outward. Measuring the density and flux through ∂Ω is reasonable based on physical meaning of the variables. We cannot measure a directly because it is a non-physical auxiliary variable introduced specifically for an efficient representation of nonlocal interactions. We assume that the ground truth parameters (κ, µ) represent a disturbance of background parameters (κ 0 , µ 0 ). Therefore, given an additional parameter ε ≥ 0, we would also like to have a regularization term in the form of R β = χ β≥ε (β). We also write R(κ, µ) = It is also possible to have other choices of regularization for (κ, µ), such as T V, H 1 , Wavelet norms. Now, we can formulate the inverse MFG as follows: (κ, µ) ) . We propose an inverse algorithm adapted from the three-operator splitting method [8] , which has also been shown to predict Nash equilibria in traffic flows [12] . We also discuss stabilizing techniques that are essential in practice. 4.1. The three-operator splitting scheme. Denoting by λ (ρ,m) := (λ ρ , λ m ) and λ (a,φ) := (λ a , λ φ ) and applying the framework in Section 2 to (3.4) we obtain the following inclusion formulation of the inverse MFG problem: The three-operator splitting scheme in [8] applies to optimization problems of the form find z ∈ H such that 0 ∈ Az + Bz + Cz, (4.2) where A, B, C are maximal monotone operators defined on a Hilbert space H, and C is cocoercive. Denote by I H the identity map in H, and J S := (I + S) −1 the resolvent of a monotone operator S. The splitting scheme for solving (4.2) can be summarized as follows where γ is a scalar. If an operator S is of the sub-differential forms; that is, S = ∂f S for some functional f S , the resolvent J S reduces to the proximal map x → arg min Overall, the algorithm for (4.1) follows three components of the generic framework in Section 2.3, upon some modification. In what follows, we discuss each component separately. where the right hand side denotes the standard H 1 (Ω) semi-norm. Assuming appropriate regularity of (ρ, m), we recall that the operator Π B,(ρ,m) is the restriction/trace operator onto the appropriate Sobolev space on the boundary L 2 [0, T ], With the aforementioned choice of the semi-norms, we naturally have the (formal) adjoint of Π B,(ρ,m) , Π * B,(ρ,m) , as the Dirichlet and Neumann harmonic extension operators by definition; that is, Here we use pr 1 , pr 2 to denote the projection from the noisy data. The harmonic extension is taken at each time t ∈ [0, 1] independently. In the implementation, we use a standard finite difference scheme to compute the harmonic extension on spatial grids for each time grid point. Note that if we assume κ = κ 0 to be known outside of domain Ω, the measurements of m and ∇φ are equivalent, as −κ 0 (x)ρ(x)∂ n φ(x) = m(x) · n on ∂Ω + . We remark that the techniques of harmonic extension have been applied to various other problems, e.g. over point clouds and in machine learning [29] . It is clear to see that λ φ is redundant and λ a = 0 whenever 0 ∈ A + B + C. Hence, we can consider only C (κ, µ), ((ρ, m), (a, φ)), (λ (ρ,m) , (0, 0)) . In this case, we preform a primal-dual hybrid gradient method for updating λ (ρ,m) : In this part, we focus on the update for the inverse problem variables (κ, µ). where S α (r) is the shrinkage operator given as S α (r) = sign(r) max{|r| − α, 0}, and Since we have the µ, ω ∈ C odd=even := {(x 1,1 , x 1,2 , x 2,1 , x 2,2 , ..., x r,1 , x r,2 ) : x i,1 = x i,2 ∀i = 1, ..., r}, we write the projector ∂χ C odd=even (where we identify C odd=even with R r ) as Π µ : R 2r → C odd=even ∼ = R r (x 1,1 , x 1,2 , x 2,1 , x 2,2 , ..., x r,1 , x r,2 ) → ( and its adjoint as Update of the forward problem. As for the forward problem, we use primal-dual hybrid gradient method (PDHG) [6] to update ((ρ i , m i ) , (φ i , a i )) for each event i, for 1 ≤ i ≤ N . The iterative updates contains three parts: firstly a proximal gradient descent step for (ρ i , m i ) with stepsizes (α n ρ i , α n m i ); then a proximal gradient ascent step for (φ i , a i ) of stepsizes (α n φ i , α n a i ); lastly an extrapolating step for (φ i , a i ). Note that we make the choice of norm x,t for φ, based on the General-proximal Primal-Dual Hybrid Gradient (G-prox PDHG) method [15] that can be interpreted as a preconditioning step for obtaining a mesh-size-free convergence rate for the algorithm. Overall, the computation for the forward model follows the computational method proposed in [23, 24] . = argmin (ρ,m) L (ρ 0,i , g i ), (ρ, m), (a n i , φ n i ), (κ n+1 , µ n+1 ) Assembling all three components described above, we arrive at the following algorithm for solving (4.1). Input: (ρ 0,i , g i ,r B,i ) for i = 1, ..., N , (κ 0 , µ 0 ) Output: (κ n , µ n ) for n = 1, ..., N max while iteration n < N maximal do 1.Update for the adjoint problem: Here, we discuss key numerical strategies for stabilizing Algorithm 1. We refer to Appendix 2 (in particular, Algorithm 2) for more implementation details. While the change of κ n+1 is made from the accumulation of all measurement events (through (λ n+1 ρ i , λ n+1 m i )), there is sometimes unexpected change of κ n+1 (x) that makes the algorithm highly unstable. For instance, there may be a large κ n+1 (x) at a single grid point. Moreover, we are using harmonic expansion method to update (λ n+1 ρ i , λ n+1 m i ), which causes large variances of κ(x) along the boundary ∂Ω. Therefore, we add a cut-off function and a convolution kernel to the step (4.6) to have a smoother change in κ n+1 (x) in space. Specifically, we have where T mask is a cut-off function that truncates the change of κ near ∂Ω given by for a function ξ(x) vanishing near ∂Ω. As for the convolution where the convolution kernel ψ(x) satisfies Ω ψ(x)dx = 1. On the other hand, after the inverse problem parameters (κ n+1 , µ n+1 ) are updated, we get a new pair of parameters for a set of mean-field game problems. It is unclear whether starting from (ρ n i , m n i , φ n i , a n i ) and taking the update rule (4.7) once produces physical solutions for the new mean-field game system due to highly nonlinear dependence of the solution on the system parameters. Therefore, instead of preforming one iteration for the forward problem, we apply the PDHG algorithm for the forward problem until its error reaches a preset tolerance. More specifically, at every iteration n, with new system parameters (κ n+1 , µ n+1 ), we use (ρ n i , m n i , φ n i , a n i ) as an initial guess and calculate the mean-field game solution accurately so that the primal-dual gap is smaller than residual the preset tolerance. This section demonstrates the efficiency and robustness of the inverse mean-field game algorithm with three examples. We also discuss details on the rule we used to choose the best reconstruction parameters. In order to collect our observed data of the forward problem, we solve a set of meanfield game problem (3.3) with given (ρ 0,i , g i ) and (κ, µ) by finite difference method with a mesh of size (0.05, 0.04) in space-time. Each problem is solved via primal-dual optimization approach with primal-dual gap e tol < 2e − 3. The initial density function ρ 0,i is the average of two Gaussian functions with centers x G ∈ Ω \Ω. The final cost function g(x) is smooth and has a smaller value around a single point x g ∈ Ω \Ω such that densities are concentrated in the neighborhood of x g at the final time. We want to point out that there is room to improve the initial density function and final cost function choices. We choose this set of (ρ 0,i , g i ) to ensure that the density's movement covers the domain Ω as completely as possible. We also expect the nonlocal interaction among agents to be better reflected at the partial boundary measurements by setting the initial density as two Gaussians rather than one. We only take 16 forward measurement events for each of the following numerical examples. The partial boundary measurement means that we only collect the ρ, m along the boundary ∂Ω in each event. Therefore, the resulting inverse problem is severely ill-posed. To test the robustness of our reconstruction algorithm, we add some random noise to the measurements as follows: From the noisy observed data {(ρ, m · n) δ (t i , x j )} i=1,..,I,j=1,...,J on the sampling points of the measurement surface, we then use the algorithm to reconstruct the forward problem parameters (κ, µ). Recall that we paramatrized the running cost L(x, v) := 1 2κ(x) |v| 2 by κ and nonlocal kernel K(x, y) := ζ(x; µ, ω) · ζ(y; µ, ω) by µ. Since we aim at recovering the model on a given domain with fixed grid points, we fix the choice of ω, and only seek sparse recovery of µ. In the following examples, we use a set of parameters uniformly, without tuning. γ c = 0.2, γ µ = 0.1, α c = 0.1, α µ = 0.1, α λ = 10. We set the lower-bound projection parameter ε 1 = κ c , this is based on the additional assumption of the model parameters that κ(x) ≥ κ c for x ∈ Ω . The projection parameter for kernel coefficient is ε 2 = 1. To account for unknown ground truth of the model parameters, we introduce Res to quantify the quality of the reconstructed parameters. where ρ n i , m n i are the solution of the forward mean-field game problem with i-th choice of initial density and final cost function with the reconstructed parameter (κ n , µ n ) at the n-th iteration of the algorithm. The boundary residual Res measures how much the new boundary measurements of the mean-field game model with the recovered parameters deviate from the given partial measurements. If (κ, µ) = (κ true , µ true ), we would expect that Res is close to 0. Therefore, we pick the reconstructed parameters at n opt -th iteration by taking n opt = arg min n Res n , (κ opt , µ opt ) = (κ nopt , µ nopt ). When we implemented the algorithm, we observed that the quantity Res n first decreased then increased with respect to the iteration. We also observed that with large enough number of iterations, (for example, 1500), the inverse problem is contaminated and the reconstruction of mean-field game coefficients are very bad. In the following examples, we take fixed number of iterations N = 1500 for the inverse algorithm, and pick the reconstructed model parameters accordingly. 5.2. Example 1. This example tests a running cost κ(x) with a bump at point (0.25, 0.25), which means the density that travel crossing near this point has a lower cost than other routes. The density are also expected to accelerate when they travel across this point. The nonlocal kernel K(x, y) is constructed via a Gaussian function plus some sparse terms in forms of µ 2 k cos(ω k · x). The nonlocal kernel, in general, penalizes being too concentrated. The amplify of certain Fourier frequencies determines the agents' particular interaction preferences. Specifically, we have the following: K s (x, y) = 0.2094 2 (cos (π(x 1 − y 1 )) + cos (π(x 2 − y 2 ))) + 0.2613 2 (cos (π(x 1 − y 1 ) + π(x 2 − y 2 )) + cos (−π(x 1 − y 1 ) + π(x 2 − y 2 ))) , We have µ 0 , which represents K 0 (x, y) via the expansions form µ 2 k cos(ω k · x), known. The variable k 0 is a given constant value that makes the kernel integration K(x, y)dxdy = 1. Varying this constant corresponding to changing the coefficient of the zero Fourier mode (0, 0). This constant k 0 does not change the intensity of repulsion effect among the agents, since k 0 ρ(y)dy = k 0 ρ 0 (y)dy is uniform over the domain Ω . With K(x, y) = ζ(x; µ, ω) · ζ(y; µ, ω) for µ, ω ∈ C odd=even , we omit the even entries (eg. µ k,2 , ω k,2 ) and express the kernel K s as follows: µ s = (0.2094, 0.2094, 0.2613, 0.2613), ω s = ((π, 0), (0, π), (π, π), (−π, π)). Here, we also assume that κ c = 2 and κ(x) = 2 for x ∈ Ω \Ω is known. Given (κ 0 , µ 0 ) and the noisy partial boundary measurements with corresponding event parameters (ρ 0,i , g i ,r B,i ), we apply our inverse algorithm.The results are shown in Figure 2 ,3. In Figure 2 , we plot the residual Res n and the max x κ n (x) along the iteration. We see that the residual oscillates and decreases first, then bounces back and increases. In Figure 3 , we show the reconstruction of model parameters by taking n opt = arg min n Res n , (κ opt , µ opt ) = (κ nopt , µ nopt ). We see that the reconstructed κ opt (x) has a single bump sits near (0.25, 0.25). The shape of the bump is not as sharp as the ground truth κ. The maximal value of running cost max x κ true (x) = 6; while max x κ opt (x) = 4. As for the non-local kernel, we have µ opt nicely reconstructed, where µ true = µ 0 + µ s ≈ µ opt . This example shows that our inverse algorithm is robust to noise and can recover the model parameters (κ, µ) simultaneously. Figure 3 . From left to right: the true running cost κ(x); the reconstructed running cost κ opt (x) at iteration n opt ; the coefficient representation of nonlocal kernel K(x, y) in vector form, where x-axis represents different Fourier mode ω and the y-axis corresponds to the coefficients µ. In this example, we make the κ(x) more complicated by having two bumps sitting diagonally. We except that if the density travels across these two bumps, it will accelerate twice. The model set-up is as follows: Figure 4 . From left to right: the true running cost κ(x); the reconstructed running cost κ opt (x) at iteration n opt ; the coefficient representation of nonlocal kernel K(x, y) in vector form. We can see from the Figure 4 that recovered bumps are well separated, and their locations are accurately captured. Reconstructed bumps are more spread compared to the ground truth, and there is some noise on upper left and bottom right corners of the domain Ω . The nonlocal kernel is reconstructed nicely as shown in Figure 4 (right). A precise sparse representation of K s (x, y) is recovered. Considering the severe ill-posedness of the inverse problem with 10% multiplicative noise added to the boundary measurements, the reconstruction quality is quite satisfactory. 5.4. Example 3. In this example, we modify the κ(x) by having two bumps sitting in parallel. Similar to the Example 2, the density would prefer to move crossing these bumps. We set the nonlocal kernel with K s containing Fourier modes with higher frequency. ω s = ((2π, −π), (2π, π), (π, 2π), (π, −2π)). Figure 5 . From left to right: the ground true running cost κ(x); the reconstructed running cost κ opt (x) at iteration n opt ; the coefficient representation of nonlocal kernel K(x, y) in vector form. We have the reconstruction result shown in Figure 5 . The two parallel sitting bumps are well separated and located with reasonable accuracy. Again, the bumps are diffused with some noise near the upper boundary of Ω. The nonlocal kernel is recovered very nicely. In this paper, we formulate a new class of inverse mean-field game problems given only partial boundary measurements. A novel model recovery algorithm is proposed based on the saddle point formulation of MFGs. We demonstrate the robustness and effectiveness of the numerical inverse algorithm with several examples, where the MFG model parameters are reconstructed accurately. Our algorithm can be further generalized to other inverse problems with saddle point structure in the forward problem. Partial differential equation models in macroeconomics Mean field games Random features for highdimensional nonlocal mean-field games Data assimilation in price formation Mean field game of controls and an application to trade crowding A first-order primal-dual algorithm for convex problems with applications to imaging Multi-population mean field games systems with neumann boundary conditions A three-operator splitting scheme and its optimization applications. Set-valued and variational analysis A mean field game inverse problem Economic models and mean-field games theory Regularity theory for mean-field game systems Learn to predict equilibria via fixed point networks Large-population cost-coupled LQG problems with nonuniform agents: individual-mass behavior and decentralized -Nash equilibria Large population stochastic dynamic games: closed-loop McKean-Vlasov systems and the Nash certainty equivalence principle Solving large-scale optimization problems with a convergence rate independent of grid size Inverse problem for non-viscous mean field control: Example from traffic Jeuxà champ moyen. I. Le cas stationnaire Jeuxà champ moyen. II. Horizon fini et contrôle optimal Mean field games Convergence of large population games to mean field games with interaction through the controls Controlling propagation of epidemics via mean-field control Alternating the population and control neural networks to solve high-dimensional stochastic mean-field games Computational methods for first-order nonlocal mean field games with applications Splitting methods for a class of non-potential mean field games A mean field game approach to swarming robots control Fourier approximation methods for first-order nonlocal mean-field games A bochner theorem and applications. Duke mathematical journal A machine learning framework for solving high-dimensional mean field game and mean field control problems Harmonic extension on the point cloud Mean-field games for bio-inspired collective decision-making in dynamical networks Learning deep mean field games for modeling large population behavior We denote by L(·) := L (ρ 0,i , g i ), (ρ i , m i ), (a i , φ i ), (κ, µ) , the function L defined in Equation 3 .3 for simplicity. The KKT conditions for the inverse mean-field game problem are thenFurthermore, the derivatives of L are given by Input: (ρ 0,i , g i ,r B,i ) for i = 1, ..., N , (κ 0 , µ 0 ) Output: (κ n , µ n ) for n = 1, ..., N max while iteration n < N maximal do 1.Update for the adjoint problem by computing (λ n+1 Algorithm 3 Iterative algorithm for the nonlocal mean-field game systemInput: (ρ 0 , g), (κ, µ), a set of initial guess (ρ 0 , m 0 , φ 0 , a 0 ), e tol , a set of stepsizes (α j ρ , α j m , α j φ , α j a ) Output: (ρ * , m * , φ * , a * ) while iteration j < J max and primal-dual gap PD(ρ j , m j , φ j , a j ) ≥ e tol dox,t (φ j+1,temp , a j+1,temp ) = argmin (a,φ) − L (ρ 0 , g), (ρ n+1 , m n+1 ), (a, φ), (κ, µ)(φ j+1 , a j+1 ) = 2(φ j+1,temp , a j+1,temp ) − (φ j , a j ) j ← j + 1 end while return (ρ j , m j , φ j , a j )