# A mesh-independent method for second-order potential mean field games

Kang Liu, Laurent Pfeiffer

► To cite this version:

Kang Liu, Laurent Pfeiffer. A mesh-independent method for second-order potential mean field games. IMA Journal of Numerical Analysis, 2025, 45 (2), pp.1226-1266. 10.1093/imanum/drae061. hal-04269396

HAL Id: hal-04269396

https://inria.hal.science/hal-04269396v1

Submitted on 3 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A MESH-INDEPENDENT METHOD FOR SECOND-ORDER POTENTIAL MEAN FIELD GAMES

KANG LIU¹,² AND LAURENT PFEIFFER¹

**ABSTRACT.** This article investigates the convergence of the Generalized Frank-Wolfe (GFW) algorithm for the resolution of potential and convex second-order mean field games. More specifically, the impact of the discretization of the mean-field-game system on the effectiveness of the GFW algorithm is analyzed. The article focuses on the theta-scheme introduced by the authors in a previous study. A sublinear and a linear rate of convergence are obtained, for two different choices of stepsizes. These rates have the mesh-independence property: the underlying convergence constants are independent of the discretization parameters.

## 1. INTRODUCTION

1.1. **Context and main contributions.** This article is concerned with the numerical resolution of second-order mean field games (MFGs). These models describe the asymptotic behavior of Nash equilibria in stochastic differential games, as the number of players goes to infinity. They were introduced independently in [29] and [24] and have applications in various domains, such as economics, biology, finance and social networks, see [20]. In this article, we consider the following standard second-order MFG on the space $Q := [0, 1] \times \mathbb{T}^d$,

$$ \left\{ \begin{array}{ll} (i) & -\partial_t u - \sigma \Delta u + H^c(t, x, \nabla u(t, x)) = f^c(t, x, m(t)) & (t, x) \in Q, \\ (ii) & v(t, x) = -H_p^c(t, x, \nabla u(t, x)) & (t, x) \in Q, \\ (iii) & \partial_t m - \sigma \Delta m + \text{div}(vm) = 0 & (t, x) \in Q, \\ (iv) & m(0, x) = m_0^c(x), \quad u(1, x) = g^c(x) & x \in \mathbb{T}^d, \end{array} \right. \qquad (\text{MFG}) $$

where the Hamiltonian $H^c$ is related to the Fenchel conjugate of a running cost $\ell^c: Q \times \mathbb{R}^d \to \mathbb{R}$:

$$ H^c(t, x, p) := \sup_{v \in \mathbb{R}^d} \langle -p, v \rangle - \ell^c(t, x, v). $$

The existence and uniqueness of the classical solution of (MFG) is proved in [29] under assumptions on the coupling function $f^c$.

Several discretization schemes have been proposed and analyzed for the resolution of (MFG). They consist of a backward discrete Hamilton-Jacobi-Bellman (HJB) equation and a forward discrete Fokker-Planck (FP) equation: They preserve the nature of the problem as a coupled system of two equations. An implicit finite difference scheme has been introduced in [2] and convergence results for this scheme have been obtained in [1] and [4]. Other schemes, based on semi-Lagrangian discretizations have been investigated in [16, 17, 23] for first-order and (possibly) degenerate second-order MFGs. This work will focus on a scheme called theta-scheme, recently introduced by the authors in [11]. In short, this scheme involves a Crank-Nicolson discretization of the diffusion term and an explicit discretization of the first-order non-linear term.

The resolution of the discretized coupled system is in general a difficult task. We restrict our attention to the case of potential MFGs, for which optimization methods can be leveraged. The system (MFG) is said to be potential (or variational) if there exists a function $F^c: [0, 1] \times \mathcal{D} \to \mathbb{R}$

<footnote>Date: November 1, 2023.</footnote>
such that for any $t \in [0, 1]$ and $m_1, m_2 \in \mathcal{D}$ (see the definition of $\mathcal{D}$ in (4.1)),

$$F^c(t, m_1) - F^c(t, m_2) = \int_0^1 \int_{x \in \mathbb{T}^d} f^c(t, x, m_1 + s(m_2 - m_1))(m_2(x) - m_1(x))dxds. \quad (1.1)$$

We further restrict our attention to the case of a convex potential function $F^c$. In the presence of such $F^c$, the system (MFG) can be interpreted as the first-order optimality condition of an optimal control problem driven by the FP equation,

$$\left\{
\begin{array}{l}
\displaystyle \inf_{(m,v)} J^c(m,v) := \int_Q \ell^c(t,x,v)m(t,x)dtdx + \int_0^1 F^c(t,m(t))dt + \int_{\mathbb{T}^d} g^c(x)m(T,x)dx, \\
\text{such that} \quad \left\{
\begin{array}{ll}
\partial_t m - \sigma \Delta m + \text{div}(vm) = 0, & \forall (t,x) \in Q, \\
m(0,x) = m_0^c(x), & \forall x \in \mathbb{T}^d.
\end{array}
\right.
\end{array}
\right.
\quad (\text{OC})$$

Problem (OC) is equivalent to a convex optimal control problem, obtained through the classical Benamou-Brenier transform [8]. Then, some numerical algorithms can be applied to find a solution of (OC), such as ADMM [7, 6], the Chambolle-Pock algorithm [3], the fictitious play [23] and the generalized Frank-Wolfe (GFW) algorithm [30]. Some articles propose to discretize the optimal control problem (OC), see for example [27, 6]. In this context, it is very desirable that the potential structure of the continuous MFG is preserved at the level of the discretized coupled system, so that one can apply in a direct fashion suitable optimization methods to the discrete system. This is in particular the case for the implicit scheme proposed in [1] and solved in [3] with the Chambolle-Pock algorithm. As we establish in this article, the theta-scheme of [11] also preserves the potential structure of the MFG system.

We focus in this article on the resolution of the discrete MFG system with the Generalized Frank-Wolfe (GFW) algorithm see [12]. This algorithm is an iterative method, consisting in solving at each iteration a partially linearized version of the potential problem (OC). The linearized problem to be solved is equivalent to a stochastic optimal control problem that can be solved by dynamic programming. As we will explain more in detail, this allows to interpret the GFW method as a best-response procedure. For a specific choice of stepsize, it coincides with the fictitious play method of [14]. Others works have investigated the fictitious play method for MFGs: [32] proves the convergence of the continuous method, in a discrete setting with common noise; [23] proves a general result for fully discrete MFGs, which can be applied to discretized first-order MFGs. The article [21] shows the connexion between fictitious play and the Frank-Wolfe algorithm for potential discrete MFGs.

The general objective of the article is to show that the performance of the GFW algorithm is not impacted by a refinement of the discretization grid. The main results of our article are two mesh-independence properties for the resolution of (MFG) with the theta-scheme and the GFW algorithm. The terminology mesh-independence was coined in the article [5]. It is said that an algorithm satisfy a mesh-independence property when approximately the same number of iterations is required to satisfy a stopping criterion, when comparing an infinite-dimensional problem and its discrete counterpart. In a more precise fashion, we will say that the GFW algorithm has a mesh-independent sublinear rate of convergence if there exists a constant $C > 0$, independent of the discretization parameters, such that

$$J_h(m_h^k, v_h^k) - J_h^* \le \frac{C}{k}, \quad \forall k \ge 1.$$

In the above estimate, $J_h$ denotes the discretized counterpart of $J$, $J_h^*$ denotes the value of the discretized optimal control problem, and $(m_h^k, v_h^k)$ denotes the candidate to optimality obtained at iteration $k$. Similarly, we will say that the GFW algorithm has a mesh-independent linear rate of
convergence if there exist two constants $C > 0$ and $\delta \in (0, 1)$ such that

$$J_h(m_h^k, v_h^k) - J_h^* \leq C\delta^k, \quad \forall k \geq 1.$$

We establish that the GFW algorithm has a mesh-independent sublinear (resp. a linear) rate for two different choices of stepsize. Our analysis is close to the one performed in [30], in which the sublinear and the linear convergence of the GFW algorithm is demonstrated for the continuous model and for the same choices of stepsizes as in the present study. While the sublinear convergence of the GFW method (in a general setting) is classical, the linear rate of convergence relies on recent techniques from [26]. To the best of our knowledge, in the context of mean field games, the mesh-independence property has never been established so far for any other method. Though it seems a natural property, it may not hold in general. In particular, it might not hold for primal-dual methods, whose application relies on a saddle-point formulation of the convex counterpart of (OC) of the form, in which the Fokker-Planck equation is “dualized”. This saddle-point formulation involves a linear operator, encoding the (discrete) Fokker-Planck equation (see for example [3, Sec. 3.2]). As the discretization parameters decrease, the operator norm of these operators (for the Euclidean norm) increases, which has an impact on the convergence properties of methods such as the Chambolle-Pock algorithm. In contrast, the discrete Fokker-Planck equation remains satisfied at each iteration of the GFW equation.

This article is organized as follows. In Section 2, we introduce some preliminary results and notations. In Section 3, we introduce a general class of potential discrete MFGs, containing the theta-scheme. We establish the sublinear and the linear convergence of the GFW method in this discrete setting. We give explicit formulas for the convergence constants. These constants essentially depend on the Lipschitz-modulus of the coupling function of the MFG and on two bounds, for different norms, of the solution of the discretized Fokker-Planck equation, denoted $C_1$ and $C_2$. We recall the theta-scheme in Section 4, we show that it preserves the potential structure of the MFG, and we prove that the GFW algorithm has a mesh-independent sublinear and linear rates of convergence. The technical analysis relies on precise estimates of the constants $C_1$ and $C_2$, obtained thanks to a general energy estimate and an $L^\infty$-estimate for the discrete Fokker-Planck equation.

## 1.2. Notation.
We discretize the interval $[0, 1]$ with a time step $\Delta t = 1/T$, where $T \in \mathbb{N}_+$. The time set is denoted by $\mathcal{T}$ ($\tilde{\mathcal{T}}$ if the final time step $T$ is included). Given a finite subset $S$ of $\mathbb{R}^d$, we denote by $\mathbb{R}(S)$ (resp. $\mathbb{R}^d(S)$) the set of functions from $S$ to $\mathbb{R}$ (resp. $\mathbb{R}^d$). We also denote by $\mathcal{P}(S)$ the set of probability measures over $S$. We call curve of probability measures any function $m: \tilde{\mathcal{T}} \times S \to \mathbb{R}$ such that $m(t, \cdot) \in \mathcal{P}(S)$, for any $t \in \mathcal{T}$. The set of probability curves is denoted by $\mathcal{P}(\tilde{\mathcal{T}}, S)$. In mathematical terms,

$$\mathcal{T} = \{0, 1, \dots, T-1\}, \quad \tilde{\mathcal{T}} = \{0, 1, \dots, T\};$$

$$\mathbb{R}(S) = \{m: S \to \mathbb{R}\}, \quad \mathbb{R}^d(S) = \{m: S \to \mathbb{R}^d\};$$

$$\mathcal{P}(S) = \{m \in \mathbb{R}(S) \mid \forall x \in S, m(x) \ge 0, \sum_{y \in S} m(y) = 1\};$$

$$\mathcal{P}(\tilde{\mathcal{T}}, S) = \{m \in \mathbb{R}(\tilde{\mathcal{T}} \times S) \mid \forall t \in \tilde{\mathcal{T}}, m(t, \cdot) \in \mathcal{P}(S)\}.$$

We denote by $\| \cdot \|$ and $\langle \cdot, \cdot \rangle$ the Euclidean norm and the scalar product in $\mathbb{R}^d$. Let $S_1$ and $S_2$ be two finite sets. Let $\mu \in \mathbb{R}^d(S_1 \times S_2)$. For any $x \in S_1$ and for any $p_1$ and $p_2 \in [1, \infty]$, we denote by $\|\mu(x, \cdot)\|_{p_2}$ the $L^{p_2}$-norm of the function $y \mapsto \mu(x, y)$, defined as follows:

$$\|\mu(x, \cdot)\|_{p_2} = \begin{cases} \left( \sum_{y \in S_2} \|\mu(x, y)\|^{p_2} \right)^{1/p_2}, & \text{if } p_2 \in [1, \infty), \\ \max_{y \in S_2} \|\mu(x, y)\|, & \text{if } p_2 = \infty. \end{cases}$$
We next define

$$ \|\mu\|_{p_1, p_2} = \begin{cases} \left( \sum_{x \in S_1} \|\mu(x, \cdot)\|_{p_2}^{p_1} \right)^{1/p_1}, & \text{if } p_1 \in [1, \infty), \\ \max_{x \in S_1} \|\mu(x, \cdot)\|_{p_2}, & \text{if } p_1 = \infty. \end{cases} $$

**Lemma 1.1** (Hölder's inequality). Let $S_1$ and $S_2$ be two finite sets. Let $\mu$ and $\nu \in \mathbb{R}^n(S_1 \times S_2)$ and let $p_1$ and $p_2 \in [1, \infty]$. It holds that

$$ \sum_{x_1 \in S_1} \sum_{x_2 \in S_2} \left| \langle \mu(x_1, x_2), \nu(x_1, x_2) \rangle \right| \le \|\mu\|_{p_1, p_2} \|\nu\|_{p_1^*, p_2^*}, $$

where $1/p_i + 1/p_i^* = 1$, for $i = 1, 2$.

**Definition 1.2** (Nemytskii operators). Let $\zeta$ be a function from $\mathcal{X} \times \mathcal{Y}$ to $\mathcal{Z}$ and let $u$ be a function from $\mathcal{X}$ to $\mathcal{Y}$. Then the Nemytskii operator is the mapping $\zeta[u]$ from $\mathcal{X}$ to $\mathcal{Z}$ defined by

$$ \zeta[u](x) = \zeta(x, u(x)), \quad \forall x \in \mathcal{X}. $$

## 2. POTENTIAL DISCRETE MEAN FIELD GAMES

In this section we introduce a general class of discrete MFGs containing the $\theta$-scheme of [11]. We provide a first potential formulation of the discrete MFGs and show their equivalence with a convex optimization problem, using the classical Benamou-Brenier transformation. The analysis being rather standard (see for example [9]), we mostly give succinct proofs.

**2.1. Problem formulation.** We fix $T \in \mathbb{N}_+$ and a finite subset $S$ of $\mathbb{R}^d$. For the description of the MFG model, we fix a running cost $\ell$, a coupling cost $f$, an initial condition $m_0$, and a terminal cost $g$, where

$$ \ell: \mathcal{T} \times S \times \mathbb{R}^d \to \mathbb{R} \cup \{\infty\}, \quad f: \mathcal{T} \times S \times \mathbb{R}(S) \to \mathbb{R}, \quad m_0 \in \mathcal{P}(S), \quad g \in \mathbb{R}(S). $$

For a given $m \in \mathcal{P}(\tilde{\mathcal{T}}, S)$, we denote by $\bar{\ell}_m$ the map defined by

$$ \bar{\ell}_m: (t, x, \omega) \in \mathcal{T} \times S \mapsto \ell(t, x, \omega) + f(t, x, m(t)), $$

where $m(t) = (m(t, x))_{x \in S}$. We require the following assumptions for $\ell, m, f$, and $m_0$.

**Assumption 1.** (1) *Bounded domain.* For all $(t, x) \in \mathcal{T} \times S$, $\ell(t, x, \cdot)$ is lower semi-continuous with a non-empty domain. There exists $D > 0$ such that for any $t \in \mathcal{T}$, for any $x \in S$, and for any $\|v\| > D$, we have $\ell(t, x, v) = +\infty$.

(2) *Regularity.* There exists $L_f$ such that for any $(t, x)$ and for any $m_1$, and $m_2$ in $\mathcal{P}(S)$, we have

$$ |f(t, x, m_1) - f(t, x, m_2)| \le L_f \|m_1 - m_2\|_2. $$

(3) *Strong convexity.* There exists $\alpha > 0$ such that for any $t \in \mathcal{T}$ and for any $x \in S$, the function $\ell(t, x, \cdot)$ is $\alpha$-strongly convex, i.e.,

$$ \ell(t, x, v_2) \ge \ell(t, x, v_1) + \langle p, v_2 - v_1 \rangle + \frac{\alpha}{2} \|v_2 - v_1\|^2, $$

for all $v_1$ and $v_2$ in $\mathbb{R}^n$ and for all $p \in \partial_p \ell(t, x, v_1)$.

(4) *Monotonicity.* For any $t \in \mathcal{T}$, for any $m_1$ and $m_2$ in $\mathcal{P}(S)$,

$$ \sum_{x \in S} (f(t, x, m_1) - f(t, x, m_2))(m_1(x) - m_2(x)) \ge 0. $$

We now fix two elements $\pi_0 \in \mathbb{R}(\mathcal{T} \times S^2)$ and $\pi_1 \in \mathbb{R}^d(\mathcal{T} \times S^2)$ and define the map $\pi$ by

$$ \pi: (t, x, y, \omega) \in \mathcal{T} \times S \times S \times \mathbb{R}^d \mapsto \pi_0(t, x, y) + \Delta t \langle \pi_1(t, x, y), \omega \rangle. $$
The map $\pi$ describes the probability of an agent located at time $t$ in state $x$, using the control $\omega$, to reach state $y$ at time $t+1$. Our analysis will exploit the fact that $\pi$ is affine with respect to $\omega$. Recalling the constant $D$ introduced in Assumption 1, we consider the following assumption.

**Assumption 2.** The elements $\pi_0$ and $\pi_1$ satisfy the following:

$$ 
\left\{ 
\begin{array}{ll}
\pi_0(t, x, \cdot) \in \mathcal{P}(S), & \forall (t, x) \in \mathcal{T} \times S, \\
\sum_{y \in S} \pi_1(t, x, y) = 0, & \forall (t, x) \in \mathcal{T} \times S, \\
\pi_0(t, x, y) \ge \Delta t D \|\pi_1(t, x, y)\|, & \forall (t, x, y) \in \mathcal{T} \times S \times S.
\end{array}
\right.
\qquad (2.1) $$

An immediate consequence of Assumption 2 is the following: For all $(t, x) \in \mathcal{T} \times S$, for all $\omega \in \mathbb{R}^d$, if $\|\omega\| \le D$, then $\pi(t, x, \cdot, \omega) \in \mathcal{P}(S)$.

**Assumption 3.** There exists a function $F: \mathcal{T} \times \mathcal{P}(S) \to \mathbb{R}^d$ such that for any $t \in \mathcal{T}$ and for any $m_1$ and $m_2$ in $\mathcal{P}(S)$, it holds

$$ F(t, m_1) - F(t, m_2) = \int_0^1 \sum_{x \in S} f(t, x, m_1 + s(m_2 - m_1))(m_2(x) - m_1(x))ds. $$

We have the following convexity property for the potential function $F$.

**Lemma 2.1.** For any $t \in \mathcal{T}$ and for any $m_1$ and $m_2$ in $\mathcal{P}(S)$, it holds that

$$ F(t, m_2) \ge F(t, m_1) + \sum_{x \in S} f(t, x, m_1)(m_2(x) - m_1(x)). \qquad (2.2) $$

*Proof.* Inequality (2.2) follows from the definition of $F$ and the monotonicity in Assumption 1. ■

Until the end of Section 3, we assume that Assumptions 1, 2, and 3 are satisfied. Following [11, Sec. 2.3], we consider the following discrete mean field game system, involving the variables $u \in \mathbb{R}(\tilde{\mathcal{T}} \times S)$, $v \in \mathbb{R}^d(\mathcal{T} \times S)$, and $m \in \mathcal{P}(\tilde{\mathcal{T}}, S)$:

$$ 
\left\{ 
\begin{array}{ll}
(i) & u = \mathbf{HJB}(m), \\
(ii) & v = \mathbf{V}(u), \\
(iii) & m = \mathbf{FP}(v),
\end{array}
\right.
\qquad (\text{DMFG}) $$

where the Hamilton-Jacobi-Bellman mapping **HJB**, the optimal control mapping **V**, and the Fokker-Planck mapping **FP** are defined as follows:

* Given $m \in \mathcal{P}(\tilde{\mathcal{T}}, S)$, $u = \mathbf{HJB}(m) \in \mathbb{R}(\tilde{\mathcal{T}} \times S)$ is the solution to

$$ 
\left\{ 
\begin{array}{l}
u(t, x) = \inf_{\omega \in \mathbb{R}^d} \bar{\ell}_m(t, x, \omega) \Delta t + \sum_{y \in S} \pi(t, x, y, \omega) u(t+1, y), \quad \forall (t, x) \in \mathcal{T} \times S, \\
u(T, x) = g(x), \quad \forall x \in S.
\end{array}
\right.
\qquad (2.3) $$

* Given $u \in \mathbb{R}(\mathcal{T} \times S)$, $v = \mathbf{V}(u) \in \mathbb{R}^d(\mathcal{T} \times S)$ is defined by

$$ v(t, x) = \underset{\omega \in \mathbb{R}^d}{\operatorname{argmin}} \ell(t, x, \omega) \Delta t + \sum_{y \in S} \pi(t, x, y, \omega) u(t+1, y), \quad \forall (t, x) \in \mathcal{T} \times S. \qquad (2.4) $$

* Given $v \in \mathbb{R}^d(\mathcal{T} \times S)$, $m = \mathbf{FP}(v) \in \mathbb{R}(\tilde{\mathcal{T}} \times S)$ is defined as the solution to

$$ 
\left\{ 
\begin{array}{ll}
m(t+1, y) = \sum_{x \in S} \pi(t, x, y, v(t, x))m(t, x), & \forall (t, y) \in \mathcal{T} \times S, \\
m(0, x) = m_0(x), & \forall x \in S.
\end{array}
\right.
\qquad (2.5) $$

**Lemma 2.2** (Continuity of **HJB**). For any $m_1$ and $m_2$ in $\mathcal{P}(\tilde{\mathcal{T}}, S)$, we have

$$ \|\mathbf{HJB}(m_1) - \mathbf{HJB}(m_2)\|_{\infty, \infty} \le L_f \|m_1 - m_2\|_{\infty, 2}. \qquad (2.6) $$
Proof. See [11, Eq. A.6].

<figure><img src="image_2.png" /></figure>

**Lemma 2.3** (Continuity of FP). Let $v_1$ and $v_2$ in $\mathbb{R}^d(\mathcal{T} \times S)$ be such that $\|v_1\|_{\infty,\infty} \le D$ and $\|v_2\|_{\infty,\infty} \le D$. There exists a constant $C$, independent of $v_1$ and $v_2$, such that

$$ \|\mathbf{FP}(v_1) - \mathbf{FP}(v_2)\|_{\infty,2}^2 \le C\Delta t \sum_{t \in \mathcal{T}} \sum_{x \in S} \|(v_1 - v_2)m_1(t, x)\|^2. \quad (2.7) $$

As an immediate consequence, the mapping **FP** is continuous.

Proof. Let $m_1 = \mathbf{FP}(v_1)$ and $m_2 = \mathbf{FP}(v_2)$. Let $\mu = m_1 - m_2$. It is easy to verify that

$$ \left\{ \begin{array}{l} \mu(t+1, y) = \sum_{x \in S} \pi(t, x, y, v_1(t, x))\mu(t, x) + \Delta t \sum_{x \in S} \langle \pi_1(t, x, y), (v_2 - v_1)m_2(t, x) \rangle, \\ \mu(0, y) = 0. \end{array} \right. \qquad (2.8) $$

Inequality (2.7) immediately follows from Gronwall's inequality, keeping in mind that all norms are equivalent on the finite-dimensional vector space $\mathbb{R}^d(\mathcal{T} \times S)$. ■

<figure><img src="image_3.png" /></figure>

**Theorem 2.4** (Existence). Under Assumptions 1 and 2, system (DMFG) has a solution.

Proof. We follow the proof of [11, Thm. 3.6]. We note first that by Assumption 2, the composed mapping **FP** ◦ **V** ◦ **HJB** is valued in $\mathcal{P}(\tilde{\mathcal{T}}, S)$. By Brouwer's fixed point theorem, it suffices to show that **FP** ◦ **V** ◦ **HJB** is continuous. The continuity of **HJB** and **FP** was established in Lemmas 2.2-2.3. The continuity of **V** is deduced from the strong convexity of $\ell$, see step 2 of the proof of [11, Thm. 3.6]. ■

**2.2. Potential formulation.** Similarly to the case of continuous MFGs (see [29, 8] for example), the system (DMFG) has a potential formulation. Consider the following optimal control problem:

$$ \inf_{\substack{m \in \mathcal{P}(\tilde{\mathcal{T}}, S) \\ v \in \mathbb{R}^d(\mathcal{T} \times S)}} J(m, v), \quad \text{subject to: } (m, v) \in \mathcal{A}, \qquad (P) $$

where the cost function $J$ and the set $\mathcal{A}$ are defined by

$$ J(m, v) = \Delta t \sum_{t \in \mathcal{T}} \sum_{x \in S} \ell[v](t, x)m(t, x) + \Delta t \sum_{t \in \mathcal{T}} F(t, m(t)) + \sum_{x \in S} g(x)m(T, x); \\ \mathcal{A} = \{ (m, v) \in \mathcal{P}(\tilde{\mathcal{T}}, S) \times \mathbb{R}^d(\mathcal{T} \times S) \mid m = \mathbf{FP}(v), \|v\|_{\infty,\infty} \le D \}. $$

Problem (P) is a non-convex problem which can be made convex with the classical Benamou-Brenier transform (see [8] for example) defined by

$$ \chi: (m, v) \in \mathcal{A} \mapsto (m, mv) \in \mathcal{P}(\tilde{\mathcal{T}}, S) \times \mathbb{R}^d(\mathcal{T} \times S). \qquad (2.9) $$

Here $mv$ is the pointwise product of $m$ and $v$: $mv(t, x) := m(t, x)v(t, x)$, for all $t \in \mathcal{T}$ and for all $x \in S$. We set $\tilde{\mathcal{A}} = \chi(\mathcal{A})$ and consider the cost function $\tilde{J}$, defined by

$$ \tilde{J}(m, w) = \Delta t \sum_{t \in \mathcal{T}} \sum_{x \in S} \tilde{\ell}[m, w](t, x) + \Delta t \sum_{t \in \mathcal{T}} F(t, m(t)) + \sum_{x \in S} g(x)m(T, x), $$

where the function $\tilde{\ell}[m, w]: \mathcal{T} \times S \to \mathbb{R}$ is defined by

$$ \tilde{\ell}[m, w](t, x) = \begin{cases} \ell(t, x, \frac{w(t, x)}{m(t, x)})m(t, x), & \text{if } m(t, x) \neq 0, \\ 0, & \text{if } m(t, x) = w(t, x) = 0, \\ +\infty, & \text{otherwise.} \end{cases} $$

The new problem of interest is

$$ \inf_{\substack{m \in \mathcal{P}(\tilde{\mathcal{T}}, S) \\ w \in \mathbb{R}^d(\mathcal{T} \times S)}} \tilde{J}(m, w), \quad \text{subject to: } (m, w) \in \tilde{\mathcal{A}}. \qquad (\tilde{P}) $$
In the rest of the section, we investigate some properties of Problems (P) and (P̃), as well as their relationship with (DMFG)

**Lemma 2.5.** It holds that $\text{val}(P) = \text{val}(\tilde{P})$.

*Proof.* It follows from the definitions of $\mathcal{A}$, $J$, $\tilde{J}$, and $\chi$ that for any $(m, v) \in \mathcal{A}$, we have $\tilde{J}(\chi(m, v)) = J(m, v)$. The lemma follows immediately. ■

We next discuss the convexity of Problem (P̃).

**Lemma 2.6 (Convexity).** The set $\tilde{\mathcal{A}}$ is convex. The cost function $\tilde{J}(m, w)$ is convex.

*Proof.* Take any $(m_1, w_1)$, $(m_2, w_2)$ in $\tilde{\mathcal{A}}$ and any $\lambda \in (0, 1)$. By the definition of $\tilde{\mathcal{A}}$, there exist $v_1, v_2 \in \mathbb{R}^d(\mathcal{T} \times S)$, such that $(m_i, v_i) \in \mathcal{A}$ and $w_i = m_i v_i$ for $i = 1, 2$. Let

$$m = \lambda m_1 + (1 - \lambda) m_2,$$

$$w = \lambda w_1 + (1 - \lambda) w_2,$$

$$v(t, x) = \begin{cases} 0, & \text{if } m_1(t, x) = m_2(t, x) = 0, \\ \frac{\lambda m_1 v_1 + (1-\lambda)m_2 v_2}{\lambda m_1 + (1-\lambda)m_2}(t, x), & \text{otherwise.} \end{cases}$$

We can check that $(m, v) \in \mathcal{A}$ and that $(m, w) = \chi(m, v)$. The convexity of $\tilde{\mathcal{A}}$ follows.

The function $\tilde{J}$ is defined as the sum of three terms. The last one is linear, thus convex. The second one is also convex, by Lemma 2.1. Finally, for any $(t, x)$, the map $(m, w) \mapsto \tilde{\ell}[m, w](t, x)$ is convex (see [19, Proposition 2.3]). The convexity of $\tilde{J}$ follow. ■

Given $m' \in \mathcal{P}(\tilde{T}, S)$, we consider the cost function $\tilde{J}_{m'}$, defined by

$$\tilde{J}_{m'}(m, w) = \Delta t \sum_{t \in \mathcal{T}} \sum_{x \in S} \tilde{\ell}[m, w](t, x) + f(t, x, m'(t))m(t, x) + \sum_{x \in S} g(x)m(T, x),$$

for $(m, w) \in \mathcal{P}(\tilde{T}, S) \times \mathbb{R}^d(\mathcal{T} \times S)$. We will regard $\tilde{J}_{m'}$ as a partial linearization of $J$ around $m'$. We define the corresponding optimal control problem:

$$\inf_{(m,w) \in \tilde{\mathcal{A}}} \tilde{J}_{m'}(m,w). \qquad (P_{m'})$$

**Lemma 2.7.** Let $(m_1, w_1)$ and $(m_2, w_2)$ be in $\tilde{\mathcal{A}}$. Then

$$\tilde{J}(m_2, w_2) - \tilde{J}(m_1, w_1) \geq \tilde{J}_{m_1}(m_2, w_2) - \tilde{J}_{m_1}(m_1, w_1).$$

*Proof.* This is an immediate consequence of the definitions of $\tilde{J}$, $\tilde{J}_{m'}$, and Lemma 2.1. ■

The next lemma provides us with a solution to Problem ($P_{m'}$).

**Lemma 2.8.** Let $m' \in \mathcal{P}(\tilde{T}, S)$. Let us set $\tilde{v} = \mathbf{V} \circ \mathbf{HJB}(m')$, $\tilde{m} = \mathbf{FP} \circ \mathbf{V} \circ \mathbf{HJB}(m')$, and $\tilde{w} = \tilde{m}\tilde{v}$. Then $(\tilde{m}, \tilde{w})$ is the unique solution to ($P_{m'}$). Moreover, for any $(m, w) \in \tilde{\mathcal{A}}$ and for any $v$ such that $w = mv$, it holds that

$$\tilde{J}_{m'}(m, w) - \tilde{J}_{m'}(\tilde{m}, \tilde{w}) \geq \frac{\alpha}{2} \Delta t \sum_{t \in \mathcal{T}} \sum_{x \in S} \|(v - \tilde{v})(t, x)\|^2 m(t, x). \qquad (2.10)$$

*Proof.* The first inequality is proved in [11, Sec. 3.5, page 14]. It is of similar nature to the fundamental equality of [1]. As a consequence of (2.10), the pair $(\tilde{m}, \tilde{w})$ is a solution to ($P_{m'}$). It remains to prove uniqueness. Let $(\hat{m}, \hat{w}) \in \tilde{\mathcal{A}}$ be a solution to ($P_{m'}$). Let $\hat{v}$ be such that $\hat{w} = \hat{m}\hat{v}$. Then, by inequality (2.10), we have $\sum_{t \in \mathcal{T}} \sum_{x \in S} \|\hat{v}(t, x) - \tilde{v}(t, x)\|^2 \hat{m}(t, x) = 0$. It follows next that $(\hat{v} - \tilde{v})(t, x)\hat{m}(t, x) = 0$ for all $(t, x)$. Applying Lemma 2.3, we immediately obtain that $\hat{m} = \tilde{m}$. Finally, we have $\hat{w} - \tilde{w} = \hat{m}\hat{v} - \tilde{m}\tilde{v} = \hat{m}(\hat{v} - \tilde{v}) = 0$, which concludes the proof of uniqueness. ■
**Lemma 2.9.** System (DMFG) has a unique solution $(\bar{m}, \bar{u}, \bar{v})$. Moreover, $(\bar{m}, \bar{w}) := \chi(\bar{m}, \bar{v})$ is the unique solution to $(\tilde{P})$.

*Proof.* Let $(m', u', v')$ be a solution to (DMFG). Let $w' = m'v'$. Combining Lemma 2.7 and Lemma 2.8, we deduce that for any $(m, w) \in \tilde{\mathcal{A}}$,

$$ \tilde{J}(m, w) - \tilde{J}(m', w') \geq \tilde{J}_{m'}(m, w) - \tilde{J}_{m'}(m', w') \geq 0. $$

Thus $(m', w')$ is a solution to $(\tilde{P})$.

We next prove that $(m', w')$ is the unique solution to $(\tilde{P})$. Let $(m, w)$ be a solution to $(\tilde{P})$. The above inequality shows that $(m, w)$ is also a solution to $(P_{m'})$. Thus by Lemma 2.8, $(m, w) = (m', w')$.

It remains to prove the uniqueness of the solution to (DMFG). Let $(m, u, v)$ be a solution to (DMFG). As was proved above, $(m, mv)$ is a solution to $(\tilde{P})$ and therefore $m = m'$. It follows that $u = \mathbf{HJB}(m) = \mathbf{HJB}(m') = u'$ and that $v = \mathbf{V}(u) = \mathbf{V}(u') = v'$, which concludes the proof. ■

### 3. GENERALIZED FRANK-WOLFE ALGORITHM: THE DISCRETE CASE

We investigate in this section the convergence of the GFW algorithm, applied to the (convex) potential problem $(\tilde{P})$. In this section, Assumptions 1-3 are supposed to be satisfied. We recall that (DMFG) has a unique solution $(\bar{u}, \bar{v}, \bar{m})$ and that by Lemma 2.9, $(\bar{m}, \bar{w}) = \chi(\bar{m}, \bar{v})$ is the unique solution of problem $(\tilde{P})$.

**3.1. Algorithm and convergence results.** We first define the mapping $\mathbf{BR}: \mathcal{P}(\tilde{T}, S) \to \tilde{\mathcal{A}}$. Given $m' \in \mathcal{P}(\tilde{T}, S)$, we obtain $(\tilde{m}, \tilde{w}) = \mathbf{BR}(m')$ by successively computing

$$ \tilde{v} = \mathbf{V} \circ \mathbf{HJB}(m'), \quad \tilde{m} = \mathbf{FP}(\tilde{v}), \quad \text{and} \quad \tilde{w} = \tilde{m}\tilde{v}. $$

We refer to $\mathbf{BR}$ as the best-response mapping: given a prediction $m'$ of the equilibrium distribution of the agents, $\tilde{v}$ (as defined above) is the optimal feedback for the underlying optimal control problem and $\tilde{m}$ the resulting distribution. As was demonstrated in Lemma 2.8, $\mathbf{BR}(m')$ is also the unique solution to the linearized problem $(P_{m'})$. This allows us to write the GFW algorithm for the resolution of $(\tilde{P})$ in the form of a best-response algorithm.

<caption>**Algorithm 1:** Generalized Frank-Wolfe Algorithm</caption>

Initialization: $(m^0, w^0) \in \tilde{\mathcal{A}}$;
First iteration: $(m^1, w^1) = (\bar{m}^0, \bar{w}^0) = \mathbf{BR}(m^0)$;
for $k = 1, 2, \dots$ do
| **Step 1:** Resolution of the partial linearized problem.
| Set $(\bar{m}^k, \bar{w}^k) = \mathbf{BR}(m^k)$;
| **Step 2:** Update.
| Choose $\lambda_k \in [0, 1]$;
| Set $(m^{k+1}, w^{k+1}) = (1 - \lambda_k)(m^k, w^k) + \lambda_k(\bar{m}^k, \bar{w}^k)$;
end

Note that the choice of the stepsize $\lambda_k$ will be discussed in Proposition 3.1. We introduce now three constants, $C_1$, $C_2$, and $C_3$, that will be used for the convergence analysis of Algorithm 1. The constants $C_1$ and $C_2$ are defined by

$$ C_1 = \sup_{\|v\|_{\infty,\infty} \le D} \|\mathbf{FP}(v)\|_{\infty,2}^2 \quad \text{and} \quad C_2 = \sup_{m \in \mathcal{P}(\tilde{T}, S)} \|\mathbf{FP} \circ \mathbf{V} \circ \mathbf{HJB}(m)\|_{\infty,\infty}. \quad (3.1) $$
The finiteness of $C_1$ and $C_2$ follows from the compactness of $\mathcal{P}(\tilde{T}, S)$ and from the continuity of the three mappings **HJB**, **V**, and **FP**. Note that

$$C_1 = \sup_{(m,v) \in \mathcal{A}} \|m\|_{\infty,2}^2 = \sup_{(m,w) \in \tilde{\mathcal{A}}} \|m\|_{\infty,2}^2. \quad (3.2)$$

By Lemma 2.3, there exists a constant $C_3 > 0$ such that for any $(m_1, v_1)$ and $(m_2, v_2)$ in $\mathcal{A}$,

$$\|m_1 - m_2\|_{\infty,2}^2 \le C_3 \Delta t \sum_{t \in \mathcal{T}} \sum_{x \in S} \|(v_1 - v_2)m_1(t, x)\|^2. \quad (3.3)$$

We next introduce three other constants, $D_1$, $D_2$, and $c$, defined by:

$$D_1 = C_1 L_f |S|^{1/2}, \quad D_2 = (L_f + |S|^{1/2}) \sqrt{\frac{2C_2C_3}{\alpha}}, \quad c = \max \left\{ 1 - \frac{\alpha}{4C_2C_3L_f|S|^{1/2}}, \frac{1}{2} \right\}. \quad (3.4)$$

**Proposition 3.1.** We consider the sequence $(m_k, w_k)_{k \ge 1}$ generated by Algorithm (1).

(1) Sublinear rate. Assume that $\lambda_k = 2/(k+2)$, for all $k \ge 1$. Then,

$$\tilde{J}(m^k, w^k) - \tilde{J}(\bar{m}, \bar{w}) \le \frac{8D_1}{k}, \quad \forall k \ge 1. \quad (3.5)$$

(2) Linear rate. Assume that

$$\lambda_k = \min \left\{ \frac{\tilde{J}_{m^k}(m^k, w^k) - \tilde{J}_{m^k}(\bar{m}^k, \bar{w}^k)}{L_f |S|^{1/2} \|m^k - \bar{m}^k\|_{\infty,2}^2}, 1 \right\}, \quad (3.6)$$

for all $k \ge 1$. Then,

$$\tilde{J}(m^k, w^k) - \tilde{J}(\bar{m}, \bar{w}) \le 4D_1 c^k, \quad \forall k \ge 1. \quad (3.7)$$

The following subsection is dedicated to the proof of Proposition 3.1. The used stepsize rule (3.6) is motivated in the proof of convergence (see (3.13)).

## 3.2. Convergence analysis.

**Lemma 3.2 (A priori bounds).** For any $k \ge 1$, we have $(m^k, w^k) \in \tilde{\mathcal{A}}$. As a consequence,

$$\|m^k\|_{\infty,2}^2 \le C_1 \quad \text{and} \quad \|m^k\|_{\infty,\infty} \le C_2. \quad (3.8)$$

*Proof.* Using the definitions of $C_1$ and $C_2$ (given in (3.1)) and using (3.2), we have that

$$\|\bar{m}^k\|_{\infty,2}^2 \le C_1, \quad \|\bar{m}^k\|_{\infty,\infty} \le C_2, \quad \text{and} \quad (\bar{m}^k, \bar{w}^k) \in \tilde{\mathcal{A}}, \quad (3.9)$$

for any $k \ge 1$. To conclude the proof of the lemma, it suffices to observe that for any $k \ge 1$, $(m^k, w^k)$ is a convex combination of $(\bar{m}^k, \bar{w}^k) \in \tilde{\mathcal{A}}$ for $0 \le \kappa \le k-1$. Then we have $(m^k, w^k) \in \tilde{\mathcal{A}}$, since $\tilde{\mathcal{A}}$ is convex, by Lemma 2.6. Inequality (3.8) follows from (3.9) and from the triangle inequality. ■

**Lemma 3.3.** For any $t \in \mathcal{T}$ and for any $m_1$ and $m_2$ in $\mathcal{P}(S)$, we have

$$F(t, m_2) \le F(t, m_1) + \sum_{x \in S} f(t, x, m_1)(m_2(x) - m_1(x)) + \frac{L_f |S|^{1/2}}{2} \|m_2 - m_1\|_2^2. \quad (3.10)$$

*Proof.* Combining Assumption 3 and the Lipschitz-continuity of $f$ (Assumption 1), we deduce that

$$F(t, m_2) \le F(t, m_1) + \sum_{x \in S} f(t, x, m_1)(m_2(x) - m_1(x)) + \frac{L_f}{2} \|m_2 - m_1\|_2 \|m_2 - m_1\|_1.$$

By Hölder's inequality, we have that $\|m_2 - m_1\|_1 \le |S|^{1/2} \|m_2 - m_1\|_2$. Inequality (3.10) follows. ■
Algorithm 1 generates two sequences $(m^k, w^k)_{k \ge 0}$ and $(\bar{m}^k, \bar{w}^k)_{k \ge 0}$ in $\tilde{\mathcal{A}}$. For the analysis, we need to fix two sequences $(v^k)_{k \ge 0}$ and $(\bar{v}^k)_{k \ge 0}$ such that $w^k = m^k v^k$ and $\bar{w}^k = \bar{m}^k \bar{v}^k$. We introduce the following six sequences of positive numbers:

$$
\begin{align*}
\delta_k &= \|m^k - \bar{m}\|_{\infty,2}^2, & \bar{\delta}_k &= \|m^k - \bar{m}^k\|_{\infty,2}^2; \\
\gamma_k &= \tilde{J}(m^k, w^k) - \tilde{J}(\bar{m}, \bar{w}), & \bar{\gamma}_k &= \tilde{J}_{mk}(m^k, w^k) - \tilde{J}_{mk}(\bar{m}^k, \bar{w}^k); \\
\epsilon_k &= \Delta t \sum_{t \in \mathcal{T}} \sum_{x \in S} \|v^k - \bar{v}\|^2 m^k(t, x), & \bar{\epsilon}_k &= \Delta t \sum_{t \in \mathcal{T}} \sum_{x \in S} \|v^k - \bar{v}^k\|^2 m^k(t, x).
\end{align*}
$$

The following lemma establishes various relationships between these sequences, independent of the choice of stepsize $(\lambda_k)_{k \ge 1}$. For convenience, we fix $\lambda_0 = 1$ and make use of the fact that $(m^1, w^1) = (1 - \lambda_0)(m^0, w^0) + \lambda_0(\bar{m}^0, \bar{w}^0)$.

**Lemma 3.4.** For any choice of stepsizes $(\lambda_k)_{k \ge 1}$, we have $\bar{\delta}_k \le 4C_1$, for any $k \ge 1$. Moreover we have that

$$
\gamma_k \le \bar{\gamma}_k \quad \text{and} \quad \gamma_{k+1} \le \gamma_k - \lambda_k \bar{\gamma}_k + \lambda_k^2 \frac{L_f |S|^{1/2}}{2} \bar{\delta}_k. \qquad (3.11)
$$

We also have the following estimates:

$$
\frac{\alpha}{2C_2C_3}\delta_k \le \frac{\alpha}{2}\epsilon_k \le \gamma_k \quad \text{and} \quad \frac{\alpha}{2C_2C_3}\bar{\delta}_k \le \frac{\alpha}{2}\bar{\epsilon}_k \le \bar{\gamma}_k. \qquad (3.12)
$$

Proof. Step 1. The inequality $\bar{\delta}_k = \|m^k - \bar{m}^k\|_{\infty,2}^2 \le 4C_1$ follows from the bounds $\|m^k\|_{\infty,2}^2 \le C_1$ and $\|\bar{m}^k\|_{\infty,2}^2 \le C_1$ obtained in Lemma 3.2.

**Step 2.** We next prove that $\gamma_k \le \bar{\gamma}_k$. Recalling that $(\bar{m}^k, \bar{w}^k)$ minimizes $\tilde{J}_{mk}(\cdot)$ over $\tilde{\mathcal{A}}$, we obtain that

$$
\bar{\gamma}_k = \tilde{J}_{mk}(m^k, w^k) - \tilde{J}_{mk}(\bar{m}^k, \bar{w}^k) \ge \tilde{J}_{mk}(m^k, w^k) - \tilde{J}_{mk}(\bar{m}^k, \bar{w}^k).
$$

Using next Lemma 2.7, we deduce that $\bar{\gamma}_k \ge J(m^k, w^k) - J(\bar{m}, \bar{w}) = \gamma_k$. Let us prove the upper bound of $\gamma_{k+1}$. Since $\tilde{\ell}$ is convex (see the proof of Lemma 2.6), we have

$$
\tilde{\ell}[m^{k+1}, w^{k+1}](t, x) \le (1 - \lambda_k)\tilde{\ell}[m^k, w^k](t, x) + \lambda_k \tilde{\ell}[\bar{m}^k, \bar{w}^k](t, x).
$$

Moreover, by Lemma 3.3, we have

$$
F(t, m^{k+1}(t)) \le F(t, m^k(t)) + \lambda_k \sum_{x \in S} f(t, x, m^k(t)) (\bar{m}^k(t, x) - m^k(t, x)) + \frac{\lambda_k^2 L_f |S|^{1/2}}{2} \bar{\delta}_k.
$$

Then (3.11) follows from the above two inequalities and the definitions of $\tilde{J}$, $\tilde{J}_{mk}$, $\gamma_k$, and $\bar{\gamma}_k$.

**Step 3.** Let us prove (3.12). From the definition of $C_3$ and Lemma 3.2 we have that

$$
\begin{align*}
\delta_k &\le C_3 \delta t \sum_{t \in \mathcal{T}} \sum_{x \in S} \|(v_k - \bar{v})(t, x) m_k(t, x)\|^2 \\
&\le C_2 C_3 \delta t \sum_{t \in \mathcal{T}} \sum_{x \in S} \|(v_k - \bar{v})(t, x)\|^2 m_k(t, x) = C_2 C_3 \epsilon_k.
\end{align*}
$$

Moreover, using Lemmas 2.7 and 2.8, we obtain that

$$
\gamma_k \ge \tilde{J}_{\bar{m}}(m^k) - \tilde{J}_{\bar{m}}(\bar{m}) \ge \frac{\alpha}{2} \epsilon_k.
$$

We prove in a similar fashion that $\bar{\delta}_k \le C_2 C_3 \bar{\epsilon}_k$ and that $\bar{\gamma}_k \ge \frac{\alpha}{2} \bar{\epsilon}_k$.

**Lemma 3.5.** In Algorithm 1, let $u^k = \mathbf{HJB}(m^k)$. Then, whatever the choice of stepsizes $(\lambda_k)_{k \ge 1}$, we have

$$
\|u^k - \bar{u}\|_{\infty,\infty} + \|m^k - \bar{m}\|_{\infty,1} \le D_2 \sqrt{\gamma_k}.
$$
Proof. By Cauchy-Schwarz inequality and Lemma 3.4, we have

$$ \|m^k - \bar{m}\|_{\infty,1} \le \sqrt{|S|} \sqrt{\delta_k} \le \sqrt{|S|} \sqrt{\frac{2C_2C_3}{\alpha}} \sqrt{\gamma_k}. $$

By Lemma 2.2, we have $\|u^k - \bar{u}\|_{\infty,\infty} \le L_f \sqrt{\delta_k}$. Recalling the definition of $D_2$, we obtain the announced result. ■

**Proof of Proposition 3.1. Sublinear case.** Combining the two inequalities in (3.11) and using the bound $\delta_k \le 4C_1$, we obtain that

$$ \gamma_{k+1} \le (1 - \lambda_k)\gamma_k + 2\lambda_k^2 C_1 L_f |S|^{1/2} = (1 - \lambda_k)\gamma_k + 2\lambda_k^2 D_1. $$

In particular, for $k = 0$, we have $\lambda_0 = 1$, therefore $\gamma_1 \le 2D_1$. It is then easy to prove by induction the inequality $\gamma_k \le 8D_1/k$ (see [25, Thm. 1] for example).

**Linear case.** Note first that the stepsize rule (3.6) writes

$$ \lambda_k = \min \left\{ \frac{\bar{\gamma}_k}{L_f |S|^{1/2} \bar{\delta}_k}, 1 \right\}. $$

It is easy to verify that $\lambda_k$ minimizes the upper-bound (3.11), that is to say:

$$ \lambda_k = \underset{\lambda \in [0,1]}{\operatorname{argmin}} - \lambda \bar{\gamma}_k + \lambda^2 \frac{L_f |S|^{1/2}}{2 \bar{\delta}_k}. \quad (3.13) $$

Let $k \ge 1$. We consider the following two cases:

(1) If $\bar{\gamma}_k \ge L_f |S|^{1/2} \bar{\delta}_k$, then $\lambda_k = 1$. We deduce from (3.11) that

$$ \gamma_{k+1} \le \gamma_k - \bar{\gamma}_k + \frac{L_f |S|^{1/2}}{2} \bar{\delta}_k \le \gamma_k - \frac{\bar{\gamma}_k}{2} \le \frac{\gamma_k}{2} \le c\gamma_k, $$

where the second inequality follows from the assumption $\bar{\gamma}_k \ge L_f |S|^{1/2} \bar{\delta}_k$, the third one from the inequality $\gamma_k \le \bar{\gamma}_k$ and the last one from the fact that $c \ge 1/2$.

(2) If $\bar{\gamma}_k \le L_f |S|^{1/2} \bar{\delta}_k$, then $\lambda_k = \frac{\bar{\gamma}_k}{L_f |S|^{1/2} \bar{\delta}_k}$. We deduce from (3.11) and from the inequality $\gamma_k \le \bar{\gamma}_k$ that

$$ \gamma_{k+1} \le \gamma_k - \frac{\bar{\gamma}_k}{2L_f |S|^{1/2} \bar{\delta}_k} \bar{\gamma}_k \le \left(1 - \frac{\bar{\gamma}_k}{2L_f |S|^{1/2} \bar{\delta}_k}\right) \gamma_k. $$

By (3.12), we know that $\bar{\gamma}_k/\bar{\delta}_k \ge \frac{\alpha}{2C_2C_3}$. It follows that

$$ \gamma_{k+1} \le \left(1 - \frac{\alpha}{4C_2C_3L_f|S|^{1/2}}\right) \gamma_k \le c\gamma_k. $$

It follows that $\gamma_k \le \gamma_1 c^{k-1}$ for all $k \ge 1$. Since $c^{-1} \le 2$ and $\gamma_1 \le 2D_1$, estimate (3.7) holds true. ■

## 4. MESH-INDEPENDENT CONVERGENCE OF THE GFW ALGORITHM FOR SECOND-ORDER MFGs

We establish in this section two mesh-independence principles for the GFW algorithm, applied to a discretization of (MFG) with the $\theta$-scheme proposed in [11]. First, we recall the theta-scheme and the main convergence result of [11]. Then, we show that the potential structure of the continuous (MFG) is preserved by the theta-scheme at the discrete level. Therefore, Algorithm 1 can be applied to the theta-scheme. The mesh-independence principles are stated in Theorems 4.5 and 4.9 and demonstrated in Subsections 4.3 and 4.4.
4.1. **The theta-scheme and error estimates.** In this subsection, we recall the theta-scheme of the continuous system (MFG) investigated in [11] and its main result. Let us define

$$ \mathcal{D} = \{ \mu \in \mathbb{L}^2(\mathbb{T}^d) \mid \mu \ge 0, \int_{\mathbb{T}^d} \mu(x)dx = 1 \}. \quad (4.1) $$

We make the following assumptions on the data function $\ell^c: Q \times \mathbb{R}^d \to \mathbb{R}$, $g^c: \mathbb{T}^d \to \mathbb{R}$, $f^c: Q \times \mathcal{D} \to \mathbb{R}$, and $m^c: \mathbb{T}^d \to \mathbb{R}$ of the continuous model (MFG).

**Assumption A.** (1) *Regularity*. The function $\ell^c$ is continuously differentiable with respect to $v$. There exist positive constants $L_\ell^c$, $L_g^c$, and $L_f^c$ such that for any $(t, x) \in Q$, for any $v \in \mathbb{R}^d$, and for any $m \in \mathcal{D}$,

* $\ell^c(\cdot, x, v)$, $\ell^c(t, \cdot, v)$, and $\ell_v^c(\cdot, x, v)$ are $L_\ell^c$-Lipschitz continuous

* $g^c$ is $L_g^c$-Lipschitz continuous

* $f^c(\cdot, x, m)$, $f^c(t, \cdot, m)$, and $f^c(t, x, \cdot)$ are $L_f^c$-Lipschitz continuous (with respect to the $\|\cdot\|_{\mathbb{L}^2}$-norm for the third variable).

(2) *Strong convexity*. There exists $\alpha^c > 0$ such that for any $(t, x) \in Q$, $\ell^c(t, x, \cdot)$ is strongly convex with modulus $\alpha^c$, i.e.

$$ \ell^c(t, x, v_2) \ge \ell^c(t, x, v_1) + \langle \ell_v^c(t, x, v_1), v_2 - v_1 \rangle + \frac{\alpha^c}{2} \|v_2 - v_1\|^2, \quad \forall v_1, v_2 \in \mathbb{R}^d. $$

(3) *Monotonicity*. The function $f^c$ is monotone, i.e., for any $t \in [0, T]$, for any $m_1$ and $m_2 \in \mathcal{D}$,

$$ \int_{\mathbb{T}^d} (f^c(t, x', m_1) - f^c(t, x', m_2))(m_1(x') - m_2(x'))dx' \ge 0. $$

**Assumption B.** The continuous mean field game (MFG) admits a unique solution $(u^*, v^*, m^*)$, with $u^*, m^* \in C^{1+r/2,2+r}(Q)$ and $v^* \in C^r(Q) \cap \mathbb{L}^\infty([0, 1]; C^{1+r}(\mathbb{T}^d))$, where $r \in (0, 1)$.

Note that more explicit assumptions on the problem data allows to verify Assumption B, see [11, Appendix B].

The time and space discretization parameters are denoted $\Delta t > 0$ and $h > 0$. As in the discrete model, we suppose that $\Delta t = 1/T$ and that $h = 1/N$, for two integers $T$ and $N$. The discrete state space is defined as

$$ S = \mathbb{T}^d/h\mathbb{Z}^d = \{(i_1, i_2, \dots, i_d)h \mid i_1, \dots, i_d \in \mathbb{Z}/N\mathbb{Z}\}. $$

Given $x \in \mathbb{T}^d$, let us set $B_h(x) = \prod_{i=1}^n [x - he_i/2, x + he_i/2)$. Given $y \in \mathbb{T}^d$, we denote by $x_h[y]$ the unique point $x$ in $S$ such that $y \in B_h(x)$. We consider the mappings $\mathcal{I}_h: \mathbb{R}(\mathbb{T}^d) \to \mathbb{R}(S)$ and $\mathcal{R}_h: \mathbb{R}(S) \to \mathbb{R}(\mathbb{T}^d)$, defined as follows: For any $m^c \in \mathbb{R}(\mathbb{T}^d)$ and for any $m \in \mathbb{R}(S)$,

$$ \mathcal{I}_h(m^c)(x) = \int_{B_h(x)} m^c(y)dy, \quad \forall x \in S; $$

$$ \mathcal{R}_h(m)(y) = \frac{m(x_h[y])}{h^d}, \quad \forall y \in \mathbb{T}^d. $$

We consider the constant $M > 0$, defined by

$$ M = \frac{1}{\alpha^c} \left( 2 \max_{(t,x) \in Q} \|\ell_v^c(t, x, 0)\| + \sqrt{d}(L_\ell^c + L_f^c + L_g^c) \right). \quad (4.2) $$

Note that $M$ is independent of $\Delta t$ and $h$. We define the truncated running cost $\hat{\ell}^c$ as follows:

$$ \hat{\ell}^c(t, x, v) = \begin{cases} \ell^c(t, x, v), & \text{if } \|v\| \le M, \\ +\infty, & \text{otherwise.} \end{cases} $$
The discrete counterparts of $\ell^c$, $g^c$, $f^c$, and $m_0^c$ are defined as follows: For any $t \in \mathcal{T}$, for any $x \in S$, for any $v \in \mathbb{R}^d$, for any $m \in \mathcal{P}(S)$,

$$
\begin{align}
\ell(t, x, p) &= \hat{\ell}^c(t\Delta t, x, v) \nonumber \\
g(x) &= g^c(x), \nonumber \\
f(t, x, m) &= \frac{1}{h^d} \int_{B_h(x)} f^c(t\Delta t, y, \mathcal{R}_h(m)) dy, \tag{4.3} \\
m_0(x) &= \mathcal{I}_h(m_0^c)(x). \nonumber
\end{align}
$$

We denote by $H: \mathcal{T} \times S \times \mathbb{R}^d$ the associated Hamiltonian, defined by

$$
H(t, x, p) = \sup_{v \in \mathbb{R}^d} -\langle p, v \rangle - \ell(t, x, v). \quad (4.4)
$$

Let us now discretize the differential operators appearing in the continuous MFG system. Let $(e_i)_{i=1,...,d}$ denote the canonical basis of $\mathbb{R}^d$. The discrete Laplace, gradient and divergence operators for the centered finite-difference scheme are defined as follows:

$$
\begin{align*}
\Delta_h \mu(x) &= \sum_{i=1}^{d} \frac{\mu(x + he_i) + \mu(x - he_i) - 2\mu(x)}{h^2}, && \forall \mu \in \mathbb{R}(S), \forall x \in S, \\
\nabla_h \mu(x) &= \left( \frac{\mu(x + he_i) - \mu(x - he_i)}{2h} \right)_{i=1}^{d}, && \forall \mu \in \mathbb{R}(S), \forall x \in S, \\
\mathrm{div}_h \mu(x) &= \sum_{i=1}^{d} \frac{\mu_i(x + he_i) - \mu_i(x - he_i)}{2h}, && \forall \mu \in \mathbb{R}^d(S), \forall x \in S,
\end{align*}
$$

where $\mu_i$ is the $i$-th coordinate of $\mu$.

Finally, we fix $\theta \in (1/2, 1)$. The theta-scheme for (MFG), introduced in [11, Sec. 2.3], writes:

$$
\left\{
\begin{array}{ll}
(i) & u = \mathbf{HJB}_{\theta}(m), \\
(ii) & v = \mathbf{V}_{\theta}(u), \\
(iii) & m = \mathbf{FP}_{\theta}(v),
\end{array}
\right.
\qquad (\text{Theta-mfg})
$$

where the Hamilton-Jacobi-Bellman mapping $\mathbf{HJB}_\theta: \mathcal{P}(\tilde{\mathcal{T}}, S) \to \mathbb{R}(\tilde{\mathcal{T}} \times S)$, $m \mapsto u$, is defined by

$$
\left\{
\begin{array}{l}
-\displaystyle\frac{u(t+1,x)-u(t+1/2,x)}{\Delta t} - \theta \sigma \Delta_h u(t+1/2,x) = 0, \\
\\
-\displaystyle\frac{u(t+1/2,x)-u(t,x)}{\Delta t} - (1-\theta) \sigma \Delta_h u(t+1/2,x) + H[\nabla_h u(\cdot+1/2, \cdot)](t,x) = f(t,x,m(t)), \\
\\
u(T,x) = g(x),
\end{array}
\right.
$$

the optimal control mapping $\mathbf{V}_{\theta}: \mathbb{R}(\tilde{\mathcal{T}} \times S) \to \mathbb{R}^d(\mathcal{T} \times S)$, $u \mapsto v$, is defined by

$$
\left\{
\begin{array}{l}
-\displaystyle\frac{u(t+1,x)-u(t+1/2,x)}{\Delta t} - \theta\sigma\Delta_h u(t+1/2,x) = 0, \\
v(t,x) = -H_p(t,x, \nabla_h u(t+1/2,x)),
\end{array}
\right.
$$

and the Fokker-Planck mapping $\mathbf{FP}_{\theta}: \mathbb{R}^d(\mathcal{T} \times S) \to \mathbb{R}(\tilde{\mathcal{T}} \times S)$, $v \mapsto m$, is defined by

$$
\left\{
\begin{array}{l}
\displaystyle\frac{m(t+1,x)-m(t,x)}{\Delta t} - \theta\sigma\Delta_h m(t+1,x) - (1-\theta)\sigma\Delta_h m(t,x) + \mathrm{div}_h vm(t,x) = 0; \\
m(0,x) = m_0(x).
\end{array}
\right.
$$

Let us briefly motivate the theta-scheme. For the Fokker-Planck equation, a Crank-Nicolson discretization of the Laplace operator and an explicit discretization of the first-order term are utilized. An adjoint scheme is used for the HJB equation.
Given $u(t+1, \cdot)$, one first needs to solve an implicit scheme for the heat equation, corresponding to a diffusion equal to $\theta\sigma$. This yields the intermediate function $u(t+1/2, \cdot)$, which is not “saved” in the output $\mathbf{HJB}_{\theta}$. Then $u(t, \cdot)$ is obtained by solving an explicit scheme, containing the first order term and a diffusion equal to $(1-\theta)\sigma$.

**Remark 4.1.** The evaluation of the mapping $\mathbf{HJB}_{\theta}$, which is a crucial step in the GFW algorithm, requires to solve successively linear equations (for the implicit part) and explicit equations. It is therefore easier to implement than fully implicit schemes, which would require to solve general non-linear equations (with a policy iteration algorithm, for example).

Let us consider the following CFL condition:

$$ \Delta t \le \frac{h^2}{2d(1-\theta)\sigma}, \quad h \le \frac{2(1-\theta)\sigma}{M}. \qquad (\text{CFL}) $$

**Lemma 4.2.** Let Assumption A and condition (CFL) hold true. Then the system (Theta-mfg) is a particular case of (DMFG), with $\ell, f, m_0$ and $g$ defined by (4.3). Furthermore, Assumptions 1-2 hold with the constants

$$ D = M, \quad L_f = L_f^c h^{-d/2}, \quad \text{and} \quad \alpha = \alpha^c. $$

*Proof.* See [11, Lem. 4.1, Lem. 4.2, Thm. 4.4].

<figure><img src="image_4.png" /></figure>

Let us emphasize that the well-posedness of the mappings $\mathbf{HJB}_{\theta}$, $\mathbf{V}_{\theta}$, and $\mathbf{FP}_{\theta}$ is a corollary of the fact that (Theta-mfg) is a discrete MFG. The explicit formulas for $\pi_0$ and $\pi_1$ are provided in [11, Lemma 4.1]. We define the traces $u_h^* \in \mathbb{R}(\tilde{T} \times S)$ and $m_h^* \in \mathcal{P}(\tilde{T}, S)$ of the continuous solution $(u^*, m^*)$ as follows:

$$ u_h^*(t, x) = u^*(t\Delta, x) \quad \text{and} \quad m_h^*(t, \cdot) = \mathcal{I}_h(m^*(t\Delta t, \cdot)), \quad \forall t \in \tilde{T}, \forall x \in S. \qquad (4.5) $$

The main result of [11] is the following error estimate.

**Theorem 4.3.** Let $1/2 < \theta < 1$. Let Assumptions A-B and condition (CFL) hold true. Then (Theta-mfg) admits a unique solution $(u_h, v_h, m_h)$. There exists $C^*$ independent of $\Delta t$ and $h$ such that

$$ \|u_h - u_h^*\|_{\infty,\infty} + \|m_h - m_h^*\|_{\infty,1} \le C^* h^r. $$

*Proof.* See [11, Thm. 2.10].

<figure><img src="image_5.png" /></figure>

**4.2. GFW algorithm for the theta-scheme and main results.** Let us first give the assumption on the potential structure of $f$.

**Assumption C.** There exists a function $F^c: [0, 1] \times \mathcal{D} \to \mathbb{R}^d$ such that (1.1) holds.

The following lemma shows that the discretization of $f^c$ proposed in (4.3) preserves the potential structure of the MFG system at the discrete level.

**Lemma 4.4.** Let Assumption C hold true. Then, $f$ defined by (4.3) satisfies Assumption 3, with the primitive function $F$ defined by

$$ F(t, m) = F^c(t\Delta t, \mathcal{R}_h(m)). $$
Proof. Taking any $(t, m_1, m_2) \in \mathcal{T} \times \mathcal{P}(S)^2$, we have that

$$
\begin{align*}
F(t, m_1) - F(t, m_2) &= F^c(t\Delta t, \mathcal{R}_h(m_1)) - F^c(t\Delta t, \mathcal{R}_h(m_2)) \\
&= \int_0^1 \int_{x \in \mathbb{T}^d} f^c(t\Delta t, x, \mathcal{R}_h(m_1 + s(m_2 - m_1)))(\mathcal{R}_h(m_2)(x) - \mathcal{R}_h(m_1)(x))dxds \\
&= \int_0^1 \sum_{x \in S} \frac{1}{h^d} \int_{y \in B_h(x)} f^c(t\Delta t, y, \mathcal{R}_h(m_1 + s(m_2 - m_1)))dy(m_2(x) - m_1(x))ds \\
&= \int_0^1 \sum_{x \in S} f(t, x, m_1 + s(m_2 - m_1))(m_2(x) - m_1(x))ds,
\end{align*}
$$

where the second equality follows from Assumption C and the linearity of operator $\mathcal{R}_h$, the third equality comes from the fact that $\mathcal{R}_h(m)$ is constant in $B_h(x)$, and the last equality derives from the definition of $f$.

Following the definition of **BR** in Sec. 3, we define the best response mapping for (Theta-mfg):
for any $m' \in \mathcal{P}(\tilde{T}, S)$,

$$ \mathbf{BR}_{\theta}(m') = \chi(\tilde{m}, \tilde{v}), \quad \text{where } \tilde{v} = \mathbf{V}_{\theta} \circ \mathbf{HJB}_{\theta}(m') \text{ and } \tilde{m} = \mathbf{FP}_{\theta}(\tilde{v}). $$

The GFW algorithm writes as follows.

<caption>**Algorithm 2:** Generalized Frank-Wolfe Algorithm for (Theta-mfg)</caption>

Initialization: $m_h^0 \in \mathcal{P}(\tilde{T}, S)$;
First iteration: $(m_h^1, w_h^1) = \mathbf{BR}_\theta(m_h^0)$;
**for** $k = 1, 2, \dots$ **do**
| **Step 1:** Resolution of the partial linearized problem.
| Set $(\bar{m}_h^k, \bar{w}_h^k) = \mathbf{BR}_\theta(m_h^k)$;
| **Step 2:** Update.
| Choose $\lambda_k \in [0, 1]$;
| Set $(m_h^{k+1}, w_h^{k+1}) = (1 - \lambda_k)(m_h^k, w_h^k) + \lambda_k(\bar{m}_h^k, \bar{w}_h^k)$;
| **end**

From now on, all notations introduced in Sections 2 and 3 for general discrete MFGs will be restricted to (Theta-mfg), without the adjunction of the subscript $\theta$. For example, we will denote by $\gamma_k$ the optimality gap in Algorithm 2. The following result is our first mesh-independence principle. Recall that $(u_h, m_h)$ and $(u_h^*, m_h^*)$ have been introduced in Theorem 4.3.

**Theorem 4.5 (Sublinear rate).** Let Assumptions A-C and condition (CFL) hold true. In Algorithm 2, take $\lambda_k = 2/(k+2)$, for any $k \ge 1$. Then there exist two constants $C_\theta$ and $\bar{C}_\theta$, both independent of $\Delta t$ and $h$, such that

$$ \gamma_k \le \frac{C_\theta}{k}, \quad \forall k \ge 1, \tag{4.6} $$

and such that

$$
\begin{align}
& \|u_h^k - u_h\|_{\infty, \infty} + \|m_h^k - m_h\|_{\infty, 1} \le \frac{\bar{C}_\theta}{\sqrt{kh^{d/2}}}, \tag{4.7} \\
& \|u_h^k - u_h^*\|_{\infty, \infty} + \|m_h^k - m_h^*\|_{\infty, 1} \le \bar{C}_\theta \left( h^r + \frac{1}{\sqrt{kh^{d/2}}} \right), \nonumber
\end{align}
$$

for any $k \ge 1$. In particular, for $k \ge h^{-(2r+d/2)}$, we have

$$ \|u_h^k - u_h^*\|_{\infty, \infty} + \|m_h^k - m_h^*\|_{\infty, 1} \le 2\bar{C}_\theta h^r. \tag{4.8} $$
Remark 4.6. Let us emphasize that only the estimate (4.6) is mesh-independent. The estimates provided in (4.7) get worse as $h \to 0$.

The proof of Theorem 4.5 is given in Subsec. 4.3. It is based on Proposition 3.1, more precisely on estimates of the three fundamental constants $C_1$, $C_2$, and $C_3$. We will derive from these estimates some estimates of the constants $D_1$, $D_2$, and $c$ (as defined in (3.4)). A key point is that the estimate of the constant $D_1$ is mesh-independent. Our analysis mainly relies on a stability result for the discrete Fokker-Planck equation, for the $\ell^2$-norm.

In order to establish mesh-independent estimates of $D_2$ and $c$, so as to have a mesh-independent linear rate of convergence, we need to establish an $\ell^\infty$-estimate for the discrete Fokker-Planck, assuming that the involved vector field $v$ derives from a value function $u$. In the continuous case, such an estimate is relatively easy to deduce from the semi-concavity of $u$. The transposition of this analysis to a discrete setting is quite delicate; our analysis is restricted to the case of a separable running cost (with respect to the control variables).

A function $l: \mathbb{R}^d \to \mathbb{R}$ is said to be semi-concave if $l(x) - L\|x\|^2/2$ is concave, for some $L \ge 0$. This definition makes no sense when $l$ is a function defined on torus. In fact, one can check that a periodical function is concave if and only if it is constant. Besides the previous definition, a second one based on a quadratic inequality is introduced in [13, Def. 1.1.1] for functions on an open set. We use the latter one for our definition of semi-concavity on a torus and its discretization.

**Definition 4.7.** [Semi-concave functions on the torus] Let $L$ be a positive constant. A function $l^c: \mathbb{T}^d \to \mathbb{R}$ is said to be $L$-semi-concave if

$$l^c(x) \ge \frac{l^c(x+y) + l^c(x-y)}{2} - L\|y\|^2, \quad \forall x, y \in \mathbb{T}^d. \quad (4.9)$$

Remark 4.8. If a function $l^c: \mathbb{T}^d \to \mathbb{R}$ is $C^2$, then it is semi-concave, as can be easily verified.

We consider the following assumption.

**Assumption D.** (1) Semi-concavity. There exists $L^c > 0$, such that for any $(t, v, m) \in [0, 1] \times \mathbb{R}^d \times \mathcal{D}$ with $\|v\| \le M$, the functions $l^c(t, x, v)$, $f^c(t, x, m)$, and $g^c(x)$ are $L^c$-semi-concave with respect to $x$.

(2) Separability. There exist functions $l_i^c: Q \times \mathbb{R} \to \mathbb{R}$ for $i = 1, 2, \dots, d$, such that for any $(t, x, v) \in Q \times \mathbb{R}^d$,

$$l^c(t, x, v) = \sum_{i=1}^{d} l_i^c(t, x, v_i),$$

where $v_i$ is the $i$-th coordinate of $v$.

In the next theorem we establish a linear and mesh-independent rate of convergence for the GFW algorithm with the line-search rule (3.6).

**Theorem 4.9 (Linear rate).** Let Assumptions A-D and condition (CFL) hold true. In Algorithm 2, set $\lambda_k$ with the rule (3.6), for all $k \ge 1$. Then there exist three constants $c_\theta \in (0, 1)$, $C_\theta > 0$ and $\bar{C}_\theta$, independent of $\Delta t$ and $h$, such that

$$\gamma_k \le C_\theta c_\theta^k, \quad \forall k \ge 1, \quad (4.10)$$

and such that

$$\begin{align}
& \|u_h^k - u_h\|_{\infty, \infty} + \|m_h^k - m_h\|_{\infty, 1} \le \bar{C}_\theta c_\theta^{k/2}, \nonumber \\
& \|u_h^k - u_h^*\|_{\infty, \infty} + \|m_h^k - m_h^*\|_{\infty, 1} \le \bar{C}_\theta (h^r + c_\theta^{k/2}), \tag{4.11}
\end{align}$$

for any $k \ge 1$. In particular, taking $k \ge 2r \log(h)/\log(c_\theta)$, we have

$$\|u_h^k - u_h^*\|_{\infty, \infty} + \|m_h^k - m_h^*\|_{\infty, 1} \le 2\bar{C}_\theta h^r. \quad (4.12)$$
Note that both estimates (4.10) and (4.11) are now mesh-independent. The proof of Theorem 4.9 is given in Subsec. 4.4.

### 4.3. Proof of the sublinear rate of convergence.
We prove in this section Theorem 4.5. Assumptions A-C are supposed to be satisfied all along the subsection, as well as the condition (CFL). Our analysis relies on an energy estimate, obtained in Lemma 4.10, which allows us to find first estimates of the convergence constants $C_1, C_2, C_3, D_1$, and $D_2$ (in Lemma 4.11).

We define the forward discrete gradient as follows:

$$ \nabla_h^+ \nu(x) = \left( \frac{\nu(x + he_i) - \nu(x)}{h} \right)_{i=1}^d, \quad \forall \nu \in \mathbb{R}(S), \forall x \in S. $$

**Lemma 4.10 (Energy estimate).** Let $\mu \in \mathbb{R}(\tilde{T} \times S)$ satisfy the following equation:

$$ \left\{ \begin{array}{l} (Id - \theta \sigma \Delta t \Delta_h) \mu(t+1) = (Id + (1-\theta)\sigma \Delta t \Delta_h) \mu(t) - \Delta t \operatorname{div}_h (v(t)\mu(t)) - \Delta t \operatorname{div}_h (\delta_v(t)), \\ \mu(0) = \mu_0, \end{array} \right. $$

where $\delta_v(t) \in \mathbb{R}^d(S)$ and $\|v\|_{\infty,\infty} \le M$. Then,

$$ \max_{t \in \tilde{T}} \|\mu(t)\|_2^2 \le c(\sigma, \theta, M) \left( \|\mu_0\|_2^2 + (1-\theta)\sigma \Delta t \|\nabla_h^+ \mu_0\|_2^2 + \frac{1}{\sigma(2\theta-1)} \sum_{\tau \in \mathcal{T}} \Delta t \|\delta_v(\tau)\|_2^2 \right), \quad (4.13) $$

where

$$ c(\sigma, \theta, M) := 1 + \frac{M^2}{\sigma(2\theta - 1)} \exp \left( \frac{M^2}{\sigma(2\theta - 1)} \right). \quad (4.14) $$

*Proof.* We deduce from the proof of [11, Prop. 4.5] that

$$ \begin{aligned} & \frac{1}{2} \left( \|\mu(t+1)\|_2^2 - \|\mu(t)\|_2^2 \right) + \theta \sigma \Delta t \|\nabla_h^+ \mu(t+1)\|_2^2 \\ & \le -(1-\theta) \sigma \Delta t \langle \nabla_h^+ \mu(t+1), \nabla_h^+ \mu(t) \rangle + \Delta t (\gamma_1 + \gamma_2), \end{aligned} \quad (4.15) $$

where

$$ \gamma_1 = \sum_{x \in S} \langle \nabla_h \mu(t+1, x), \mu v(t, x) \rangle, \quad \gamma_2 = \sum_{x \in S} \langle \nabla_h \mu(t+1, x), \delta_v(t, x) \rangle. $$

Let $\alpha_1 = \sigma(2\theta - 1) > 0$. Applying Young's inequality to the right-hand side of (4.15), we have the following inequalities:

$$ \begin{aligned} -\langle \nabla_h^+ \mu(t+1), \nabla_h^+ \mu(t) \rangle &\le \frac{1}{2} \|\nabla_h^+ \mu(t+1)\|_2^2 + \frac{1}{2} \|\nabla_h^+ \mu(t)\|_2^2; \\ \gamma_1 &\le \frac{\alpha_1}{2} \|\nabla_h \mu(t+1)\|_2^2 + \frac{1}{2\alpha_1} \|\mu v(t)\|_2^2 \le \frac{\alpha_1}{2} \|\nabla_h^+ \mu(t+1)\|_2^2 + \frac{M^2}{2\alpha_1} \|\mu(t)\|_2^2; \\ \gamma_2 &\le \frac{\alpha_1}{2} \|\nabla_h^+ \mu(t+1)\|_2^2 + \frac{1}{2\alpha_1} \|\delta_v(t)\|_2^2. \end{aligned} $$

Combining the above inequalities and (4.15), it follows that

$$ \|\mu(t+1)\|_2^2 - \|\mu(t)\|_2^2 + (1-\theta)\sigma \Delta t \left( \|\nabla_h^+ \mu(t+1)\|_2^2 - \|\nabla_h^+ \mu(t)\|_2^2 \right) \le \Delta t \left( \frac{M^2}{\alpha_1} \|\mu(t)\|_2^2 + \frac{1}{\alpha_1} \|\delta_v(t)\|_2^2 \right). $$

Summing the previous inequality over $t$, it follows that for any $t \in \mathcal{T}$,

$$ \|\mu(t+1)\|_2^2 \le \gamma + \frac{M^2 \Delta t}{\alpha_1} \sum_{\tau=0}^{t} \|\mu(\tau)\|_2^2, $$

where

$$ \gamma = \|\mu_0\|_2^2 + (1 - \theta) \sigma \Delta t \|\nabla_h^+ \mu_0\|_2^2 + \frac{1}{\alpha_1} \sum_{\tau \in \mathcal{T}} \Delta t \|\delta_v(\tau)\|_2^2. $$
We deduce from the discrete Gronwall inequality [18] that

$$ \max_{t \in \bar{T}} \|\mu(t)\|_2^2 \le \gamma + T\gamma \frac{M^2 \Delta t}{\alpha_1} \exp \left( \sum_{\tau \in \bar{T}} \frac{M^2 \Delta t}{\alpha_1} \right) = \gamma \left( 1 + \frac{M^2}{\alpha_1} \exp \left( \frac{M^2}{\alpha_1} \right) \right). $$

The conclusion follows.

We define next two constants $E_1$ and $E_3$, both independent of $\Delta t$ and $h$:

$$ E_1 = c(\sigma, \theta, M) \left( \|m_0^c\|_{L^\infty}^2 + \frac{1}{2} \|\nabla m_0^c\|_{L^\infty}^2 \right) \quad \text{and} \quad E_3 = \frac{c(\sigma, \theta, M)}{\sigma(2\theta - 1)}. $$

A constant $E_2$ will be introduced later on.

**Lemma 4.11.** The constants $C_1$, $C_2$, and $C_3$ (as defined in (3.1)-(3.3)) satisfy the following inequalities:

$$ C_1 \le E_1 h^d, \quad C_2 \le \sqrt{C_1} \le E_1^{1/2} h^{d/2}, \quad C_3 \le E_3. \quad (4.16) $$

As a consequence, for the constants $D_1$ and $D_2$ defined in (3.4), we have

$$ D_1 \le E_1 L_f^c, \quad D_2 \le (1 + L_f^c) \sqrt{\frac{2E_1^{3/2} E_3 L_f^c}{\alpha^c}} h^{-d/4}. \quad (4.17) $$

*Proof.* The condition (CFL) implies that

$$ d(1 - \theta)\sigma\Delta t \le \frac{h^2}{2} \le \frac{1}{2}. \quad (4.18) $$

We let the reader verify that

$$ \|m_0\|_2^2 + (1 - \theta)\sigma\Delta t \|\nabla_h^+ m_0\|_2^2 \le \left( \|m_0^c\|_{L^\infty}^2 + d(1 - \theta)\sigma\Delta t \|\nabla m_0^c\|_{L^\infty}^2 \right) h^d. $$

Combining the above two estimates, we deduce that

$$ \|m_0\|_2^2 + (1 - \theta)\sigma\Delta t \|\nabla_h^+ m_0\|_2^2 \le \left( \|m_0^c\|_{L^\infty}^2 + \frac{1}{2} \|\nabla m_0^c\|_{L^\infty}^2 \right) h^d. $$

Using this inequality in Lemma 4.10, applied with $\delta_v(t) = 0$ and $\mu_0 = m_0$, we obtain the estimate of $C_1$. Then, the estimate of $C_2$ is deduced from Hölder's inequality. The estimate of $C_3$ follows from equality (2.8) and Lemma 4.10 by taking $\mu_0 = 0$ and $\delta_v = (v_1 - v_2)m_1$. Finally, (4.17) follows from (3.4) and the previous estimates.

**Proof of Theorem 4.5.** Inequality (4.6) directly follows from Proposition 3.1 and from the estimate $D_1 \le E_1 L_f^c$, with $C_\theta = 8E_1 L_f^c$. Then, using Lemma 3.5 and the estimate of $D_2$ obtained in (4.17), we deduce that

$$ \|u^k - u_h\|_{\infty, \infty} + \|m^k - m_h\|_{\infty, 2} \le D_2 \sqrt{\gamma_k} \le 4E_1 L_f^c (1 + L_f^c) \sqrt{\frac{E_1 E_3}{\alpha^c}} \frac{1}{\sqrt{kh^{d/2}}}. $$

Recall that $C^*$ was introduced in Theorem 4.3. The inequalities in (4.7) hold true with

$$ \bar{C}_\theta = \max \left\{ 4E_1 L_f^c (1 + L_f^c) \sqrt{\frac{E_1 E_3}{\alpha^c}}, C^* \right\}. $$

The theorem is proved.
**4.4. Proof of the linear rate of convergence.** We prove in this subsection Theorem 4.9. Assumptions A-D are supposed to be satisfied all along the subsection, as well as the condition (CFL). At a technical level, we look for a more precise estimate of the constant $C_2$ (see Lemma 4.22), using an $l^\infty$-stability result for the map $FP_\theta \circ V_\theta \circ HJB_\theta$ (see Lemma 4.21), whose proof is inspired from the continuous case (see for example [15, Lem. 5.3]). Our analysis begins with a series of technical lemmas which will allow us to establish the semi-concavity of the value function.

**Definition 4.12.** [Semi-concave functions on $S$] A function $l: S \to \mathbb{R}$ is said to be $L$-semi-concave if

$$l(x) \geq \frac{l(x+y) + l(x-y)}{2} - L\|y\|^2, \quad \forall x, y \in S. \qquad (4.19)$$

**Lemma 4.13.** Let $L$, $L_1$, and $L_2$ be positive constants. The following statements hold:

(1) Let $l_1: S \to \mathbb{R}$ be $L_1$-semi-concave and $l_2: S \to \mathbb{R}$ be $L_2$-semi-concave. For any $\lambda_1, \lambda_2 \geq 0$, the function $\lambda_1 l_1 + \lambda_2 l_2$ is $(\lambda_1 L_1 + \lambda_2 L_2)$-semi-concave.

(2) Let $\{l_\omega: S \to \mathbb{R}\}_{\omega \in \Omega}$ be a family of $L$-semi-concave functions. Let $l: S \to \mathbb{R}$, $l(x) = \inf_{\omega \in \Omega} l_\omega(x)$. Suppose that for all $x \in S$, it holds that $l(x_0) > -\infty$. Then $l$ is $L$-semi-concave.

**Proof.** The first point is obtained by the definition (4.19) and the non-negativity of $\lambda_1$ and $\lambda_2$. Let us prove the second point. Let $x_0 \in S$. We have $l(x_0) > -\infty$. Then for any $\epsilon > 0$, there exists $\omega_\epsilon \in \Omega$ such that

$$l(x_0) \geq l_{\omega_{\epsilon}}(x_0) - \epsilon \geq \frac{l_{\omega_{\epsilon}}(x_0 + y) + l_{\omega_{\epsilon}}(x_0 - y)}{2} - L\|y\|^2 - \epsilon \geq \frac{l(x_0 + y) + l(x_0 - y)}{2} - L\|y\|^2 - \epsilon, \\ \text{where the second inequality follows from the semi-concavity of } l_{\omega_{\epsilon}}. \text{ We deduce the } L\text{-semi-concavity at the point } x_0 \text{ by the arbitrariness of } \epsilon. \blacksquare$$

**Lemma 4.14.** The functions $l$, $f$, and $g$, defined in (4.3), are $L^c$-semi-concave with respect to $x$.

**Proof.** The semi-concavity of $l$ and $g$ is a direct consequence of their definitions in (4.3). Let us prove the semi-concavity of $f$. By the definition of $f$ in (4.3), taking any $t \in T$, $m \in P(T, S)$, and $x, y \in S$, we have

$$\begin{align*}
f(t, x + y, m) + f(t, x - y, m) \\
&= \frac{1}{h^d} \int_{z \in B_h(x+y)} f^c(t, z, \mathcal{R}_h(m)) dz + \frac{1}{h^d} \int_{z \in B_h(x-y)} f^c(t, z, \mathcal{R}_h(m)) dz \\
&= \frac{1}{h^d} \int_{z \in B_h(0)} f^c(t, x + y + z, \mathcal{R}_h(m)) + f^c(t, x - y + z, \mathcal{R}_h(m)) dz \\
&\leq \frac{1}{h^d} \int_{z \in B_h(0)} 2f^c(t, x + z, \mathcal{R}_h(m)) + 2L^c\|y\|^2 dz \\
&= \frac{1}{h^d} \int_{z \in B_h(x)} 2f^c(t, z, \mathcal{R}_h(m)) dz + 2L^c\|y\|^2 = 2f(t, x, m) + 2L^c\|y\|^2.
\end{align*}$$

The conclusion follows.

**Definition 4.15.** Let $A$ be a function from $\mathbb{R}(S)$ to $\mathbb{R}$. We say that $A$ is translation invariant if for any $X \in \mathbb{R}(S)$ and any $y \in S$,

$$A(X(\cdot)) = A(X(\cdot - y)).$$

**Lemma 4.16.** Suppose that $X, Y \in \mathbb{R}(S)$ satisfy the following equation for some $c > 0$:

$$(Id - c\Delta t\Delta_h)Y = X. \qquad (4.20)$$

Let $A: \mathbb{R}(S) \to \mathbb{R}$ be a translation invariant function. We have the following statements:
* If A is convex and l.s.c, then $A(Y) \le A(X)$.

* If A is concave and u.s.c, then $A(Y) \ge A(X)$.

*Proof.* For any $Z \in \mathbb{R}(S)$ and $y \in S$, we define $\tau^y Z \in \mathbb{R}(S)$ by

$$ \tau^y Z(\cdot) = Z(\cdot - y). $$

Let $\gamma = c\Delta t/h^2$. We define $\mathbb{S}_X: \mathbb{R}(S) \to \mathbb{R}(S)$,

$$ \mathbb{S}_X(\mu) = \frac{1}{1 + 2d\gamma} \left\{ X + \gamma \sum_{j=1}^{d} \left( \tau^{he_j} \mu + \tau^{-he_j} \mu \right) \right\}. $$

By the proof of [11, Lemma 2.6], $\mathbb{S}_X$ is a contraction mapping and $Y$ is the fixed point of $\mathbb{S}_X$.

Suppose that A is l.s.c. and convex. Suppose that $\mu \in \mathbb{R}(S)$ satisfy $A(\mu) \le A(X)$. By the convexity of A, we have

$$ \begin{aligned} A(\mathbb{S}_X(\mu)) &\le \frac{1}{1 + 2d\gamma} \left\{ A(X) + \gamma \sum_{j=1}^{d} \left( A(\tau^{he_j} \mu) + A(\tau^{-he_j} \mu) \right) \right\} \\ &= \frac{1}{1 + 2d\gamma} \left\{ A(X) + \gamma \sum_{j=1}^{d} (A(\mu) + A(\mu)) \right\} \le A(X), \end{aligned} $$

where the second line follows from the translation invariance of A. Therefore, $A(\mathbb{S}_X^k(X)) \le A(X)$ for any $k \ge 1$. Since $Y = \lim_{k \to \infty} \mathbb{S}_X^k(X)$, by the lower-semi-continuity of A, we have

$$ A(Y) \le \liminf_{k \to \infty} A(\mathbb{S}_X^k(X)) \le A(X). $$

For the case where A is u.s.c and concave, it suffices to apply the previous result to $-A$.

**Lemma 4.17.** Let $X, Y \in \mathbb{R}(S)$ satisfy (4.20). Then, the following statements hold.

(1) Maximum/minimum principle:

$$ \min_{x \in S} X(x) \le \min_{x \in S} Y(x) \le \max_{x \in S} Y(x) \le \max_{x \in S} X(s). $$

(2) Conservation of the mass:

$$ \sum_{x \in S} Y(s) = \sum_{x \in S} X(s). $$

(3) Conservation of the Lipschitz constant: If $X(x)$ is L-Lipschitz, then $Y(x)$ is L-Lipschitz.

(4) Conservation of the semi-concavity constant: If $X(x)$ is L-semi-concave, then $Y(x)$ is L-semi-concave.

*Proof.* We use Lemma 4.16 for the proof. The key point is the choice of the translation invariant function A in Lemma 4.16. Keep in mind that the maximum (resp. minimum) of a family of linear functions is l.s.c. and convex (resp. u.s.c. and concave) in finite dimensions.

For point (1), it suffices to take $A(X) = \min_{x \in S}\{X(x)\}$ and $A(X) = \max_{x \in S}\{X(x)\}$. For point (2), we take $A(X) = \sum_{x \in S} X(x)$. For point (3), we take

$$ A(X) = \max_{x \in S} \max_{y \in S, y \ne 0} \frac{X(x+y) - X(x)}{\|y\|}. $$

Finally, for point (4), we take

$$ A(X) = \max_{x \in S} \max_{y \in S, y \ne 0} \frac{X(x+y) + X(x-y) - 2X(x)}{\|y\|^2}. $$

The conclusion follows.
**Lemma 4.18** (Semi-concavity of the value function). Let Assumptions A, D(1) and condition (CFL) hold true. Then for any $m \in \mathcal{P}(\tilde{T}, S)$, the function $u = \mathbf{HJB}_{\theta}(m)$ is $3L^c$-semi-concave with respect to $x$.

*Proof.* Observe that $\mathbf{HJB}_{\theta}(m)$ is equivalent to the formulation below:

$$ 
\left\{ 
\begin{array}{l} 
(Id - \theta \sigma \Delta t \Delta_h) u(t + 1/2) = u(t + 1); \\ 
u(t, x) = \Delta t \inf_{\|\omega\| \le M} \left\{ f(t, x, m(t)) + \ell(t, x, \omega) + \langle \omega, \nabla_h u(t + 1/2, x) \rangle \right\} \\ 
\qquad + (Id + (1 - \theta) \sigma \Delta t \Delta_h) u(t + 1/2)(x), \quad \forall x \in S; \\ 
u(T, x) = g(x), \quad \forall x \in S. 
\end{array} 
\right. 
\qquad (4.21) 
$$

We prove the lemma by induction. For $t = T$, by the terminal condition, it is obvious that

$$ u(T, x) \geq \frac{u(T, x + y) + u(T, x - y)}{2} - L^c \|y\|^2, \quad \forall x, y \in S. $$

Suppose that for some $t \in \mathcal{T}$, we have

$$ u(t + 1, x) \geq \frac{u(t + 1, x + y) + u(t + 1, x - y)}{2} - (2(T - 1 - t)\Delta t + 1)L^c \|y\|^2, \quad \forall x, y \in S. $$

Since $u(t+1)$ and $u(t+1/2)$ satisfy the implicit scheme (4.20), by Lemma 4.17(4), we have

$$ u(t + 1/2, x) \geq \frac{u(t + 1/2, x + y) + u(t + 1/2, x - y)}{2} - (2(T - 1 - t)\Delta t + 1)L^c \|y\|^2, \quad \forall x, y \in S. $$

Let $r' = (1 - \theta)\sigma\Delta t/h^2$. The second equation in (4.21) can be written as follows:

$$ u(t, x) = \inf_{\|\omega\| \le M} l_{\omega}(t, x), $$

where

$$ 
\begin{aligned} 
l_{\omega}(t, x) := & (1 - 2dr')u(t + 1/2, x) + \sum_{i=1}^{d} \left(r' + \frac{\omega_i}{2h}\right)u(t + 1/2, x + he_i) \\ 
& + \sum_{i=1}^{d} \left(r' - \frac{\omega_i}{2h}\right)u(t + 1/2, x - he_i) + \Delta t(f(t, x, m(t)) + \ell(t, x, \omega)). 
\end{aligned} 
$$

By condition (CFL), the coefficients of the above equation are positive for any $\|w\| \le M$. Then by Lemma 4.13(1), the semi-concavity of $u(t+1/2, \cdot)$, $f$, and $\ell$, we have that $l_{\omega}$ is $(2(T-t)\Delta t + 1)L^c$-semi-concave. Since $u(t,x) > -\infty$ for any $x \in S$, we deduce from Lemma 4.13(2) that

$$ u(t, x) \geq \frac{u(t, x + y) + u(t, x - y)}{2} - (2(T - t)\Delta t + 1)L^c \|y\|^2, \quad \forall x, y \in S. $$

The conclusion follows by induction.

<figure><img src="image_6.png" /></figure>

We have the following regularity result for the discrete Hamiltonian $H$ (defined in (4.4)).

**Lemma 4.19.** Let Assumptions A and D(2) hold true. Then, for any $(t, x) \in \mathcal{T} \times S$ and $\|p\| \le \sqrt{d}(L_{\ell}^c + L_{f}^c + L_{g}^c)$,

$$ H(t, x, p) = \sum_{i=1}^{d} H_{i}^{c}(t \Delta t, x, p_{i}), \quad H_{p}(t, x, p) = \left( \frac{\partial H_{i}^{c}}{\partial p_{i}}(t \Delta t, x, p_{i}) \right)_{i=1}^{d}, \qquad (4.22) $$

where $p_i$ is the i-th coordinate of p and

$$ H_{i}^{c}(t, x, p_{i}) = \sup_{v_{i}} -v_{i}p_{i} - \ell_{i}^{c}(t, x, v_{i}). \qquad (4.23) $$
Moreover, $\frac{\partial H_c^c}{\partial p_i}(t\Delta t, x, p_i)$ is $L_c^c/\alpha^c$-Lipschitz with respect to $x$.

**Proof.** Equality (4.22) is from [11, Lemma 5.1] and the separable form of $\ell^c$. The Lipschitz continuity of $\frac{\partial H_c^c}{\partial p_i}$ is proved with the same argument as the one the proof of [11, Lemma 2.7]. ■

**Lemma 4.20** (Lipschitz continuity of the value function). Let Assumption A and condition (CFL) hold true. For any $m \in \mathcal{P}(\tilde{T}, S)$, let $u = \mathbf{HJB}_\theta(m)$ and $v = \mathbf{V}_\theta(u)$. Then $u$ is $(L_c^c + L_f^c + L_g^c)$-Lipschitz with respect to $x$ and $\|v\|_{\infty,\infty} \le M$.

**Proof.** See [11, Lemma 4.3]. ■

**Lemma 4.21** ($\ell^\infty$-stability). Let Assumptions A, D and condition (CFL) hold true. Then,

$$ \sup_{\mu \in \mathcal{P}(\tilde{T}, S)} \| \mathbf{FP}_\theta \circ \mathbf{V}_\theta \circ \mathbf{HJB}_\theta(\mu) \|_{\infty, \infty} \le \exp \left( \frac{d(L_c^c + 6L^c)}{\alpha^c} \right) \|m_0\|_{\infty}. $$

**Proof.** Let $\mu \in \mathcal{P}(\tilde{T}, S)$, let $u = \mathbf{HJB}_\theta(\mu)$, let $v = \mathbf{V}_\theta(u)$, and let $m = \mathbf{FP}_\theta(v)$. Observe that $m = \mathbf{FP}_\theta(v)$ is equivalent to the formulation below:

$$ \left\{ \begin{array}{l} m(t + 1/2) = (Id + (1 - \theta)\sigma\Delta t\Delta_h)m(t) - \Delta t \mathrm{div}_h(v(t)m(t)); \\ (Id - \theta\sigma\Delta t\Delta_h)m(t + 1) = m(t + 1/2); \\ m(0) = m_0. \end{array} \right. \qquad (4.24) $$

Let us first compare $\|m(t + 1/2, \cdot)\|_{\infty}$ and $\|m(t + 1, \cdot)\|_{\infty}$. Since $m(t + 1)$ and $m(t + 1/2)$ satisfy the implicit scheme (4.20), by Lemma 4.17(1), we have

$$ \|m(t + 1, \cdot)\|_{\infty} \le \|m(t + 1/2, \cdot)\|_{\infty}. \qquad (4.25) $$

Then, we compare $\|m(t, \cdot)\|_{\infty}$ and $\|m(t + 1/2, \cdot)\|_{\infty}$. Let $r' = (1 - \theta)\sigma\Delta t/h^2$. The first equation in (4.24) shows that for any $(t, x) \in \mathcal{T} \times S$,

$$ m(t + 1/2, x) = (1 - 2d\gamma')m(t, x) + \sum_{i=1}^{d} \left( \gamma' - \Delta t \frac{v_i(t, x + he_i)}{2h} \right) m(t, x + he_i) \\ + \sum_{i=1}^{d} \left( \gamma' + \Delta t \frac{v_i(t, x - he_i)}{2h} \right) m(t, x - he_i). $$

Condition (CFL) implies that all the coefficients in the above equation are positive. Therefore, for any $(t, x) \in \mathcal{T} \times S$,

$$ m(t + 1/2, x) \le \left( 1 - \Delta t \sum_{i=1}^{d} \frac{v_i(t, x + he_i) - v_i(t, x - he_i)}{2h} \right) \|m(t, \cdot)\|_{\infty} = (1 - \Delta t \mathrm{div}_h v) \|m(t, \cdot)\|_{\infty}. \qquad (4.26) $$
By Lemma 4.20 and Lemma 4.17(3), we have $||\nabla_h u(t + 1/2, x)|| \le \sqrt{d}(L_{\ell}^c + L_{f}^c + L_{g}^c)$ for any $(t, x) \in \mathcal{T} \times S$. Then, formula (4.22) implies that

$$ 
\begin{aligned}
-\mathrm{div}_h v(t, x) &= \frac{1}{2h} \sum_{i=1}^{d} \frac{\partial H_i^c}{\partial p_i} (t\Delta t, x + he_i, (\nabla_h u(t + 1/2, x + he_i))_i) \\
&\quad - \frac{\partial H_i^c}{\partial p_i} (t\Delta t, x - he_i, (\nabla_h u(t + 1/2, x - he_i))_i) \\
&\le \frac{dL_{\ell}^c}{\alpha^c} + \frac{1}{2h} \sum_{i=1}^{d} \frac{\partial H_i^c}{\partial p_i} (t\Delta t, x, (\nabla_h u(t + 1/2, x + he_i))_i) \\
&\quad - \frac{\partial H_i^c}{\partial p_i} (t\Delta t, x, (\nabla_h u(t + 1/2, x - he_i))_i),
\end{aligned} 
\quad (4.27) $$

where the last inequality follows from the Lipschitz-continuity of $\frac{\partial H_i^c}{\partial p_i}$ with respect to $x$ (established in Lemma 4.19). Since $H_i^c$ is convex on $p_i$, the derivative $\frac{\partial H_i^c}{\partial p_i}$ is non-decreasing with respect to $p_i$. Furthermore, we know that $\frac{\partial H_i^c}{\partial p_i}$ is $1/\alpha^c$-Lipschitz on $p_i$ by the strong convexity of $\ell^c$. It follows that for any $(t, x) \in Q$ and $p_i^1, p_i^2 \in \mathbb{R}$,

$$ \frac{\partial H_i^c}{\partial p_i}(t, x, p_i^1) - \frac{\partial H_i^c}{\partial p_i}(t, x, p_i^2) \le \max \left\{ 0, \frac{1}{\alpha^c}(p_i^1 - p_i^2) \right\}. $$

Applying the above inequality to (4.27), we have

$$ 
\begin{aligned}
- \mathrm{div}_h v(t, x) &\le \frac{dL_{\ell}^c}{\alpha^c} + \frac{1}{2\alpha^c h} \sum_{i=1}^d \max \{0, (\nabla_h u(t + 1/2, x + he_i))_i - (\nabla_h u(t + 1/2, x - he_i))_i\} \\
&= \frac{dL_{\ell}^c}{\alpha^c} + \frac{1}{\alpha^c} \sum_{i=1}^d \max \left\{0, \frac{u(t + 1/2, x + 2he_i) + u(t + 1/2, x - 2he_i) - 2u(t + 1/2, x)}{4h^2}\right\}.
\end{aligned} $$

By Lemma 4.18, for any $(t, x, y) \in \mathcal{T} \times S^2$ and $y \ne 0$, we have

$$ \frac{u(t + 1/2, x + y) + u(t + 1/2, x - y) - 2u(t + 1/2, x)}{\|y\|^2} \le 6L^c. $$

Taking $y = 2he_i$, it follows that

$$ -\mathrm{div}_h v(t, x) \le \frac{d(L_{\ell}^c + 6L^c)}{\alpha^c}. \quad (4.28) $$

Combining (4.25), (4.26), and (4.28), we have

$$ \|m(t + 1, \cdot)\|_{\infty} \le \left(1 + \Delta t \frac{d(L_{\ell}^c + 6L^c)}{\alpha^c}\right) \|m(t, \cdot)\|_{\infty}. $$

Since $\Delta t = 1/T$, the conclusion follows.

<figure><img src="image_7.png" /></figure>

We are now ready to derive an improved estimate of $C_2$ (in comparison with the one in (4.16)). We define the constant $E_2$ as follows:

$$ E_2 = \exp \left( \frac{d(L_{\ell}^c + 6L^c)}{\alpha^c} \right) \|m_0^c\|_{\mathbb{L}^{\infty}}. $$

It is independent of $\Delta t$ and $h$.
**Lemma 4.22.** The constants $C_1$, $C_2$, and $C_3$ defined in (3.1)-(3.3) satisfy the following inequalities:

$$C_1 \le E_1 h^d, \quad C_2 \le E_2 h^d, \quad C_3 \le E_3. \qquad (4.29)$$

As a consequence, for the constants defined in (3.4), we have

$$D_1 \le E_1 L_f^c, \quad D_2 \le (1 + L_f^c) \sqrt{\frac{2E_1 E_2 E_3 L_f^c}{\alpha^c}}, \quad c \le c_\theta := \max \left\{ 1 - \frac{\alpha^c}{4E_2 E_3 L_f^c}, \frac{1}{2} \right\}. \qquad (4.30)$$

*Proof.* The estimates of $C_1$ and $C_3$ are the same as in (4.16), and the estimate of $C_2$ is a direct consequence of Lemma 4.21 and the regularity of $m_0^c$. Then (4.30) is deduced from (4.29) and (3.4). ■

*Proof of Theorem 4.9.* Inequality (4.10) holds true with $C_\theta = 4E_1 L_f^c$, as a direct consequence of Proposition 3.1. Inequality (4.11) is established in similar fashion to inequality (4.7). It holds true with

$$\bar{C}_\theta = \max \left\{ 2E_1 L_f^c (1 + L_f^c) \sqrt{\frac{2E_2 E_3}{\alpha^c}}, C^* \right\}.$$

The theorem is proved. ■

**4.5. Discussion on convergence constants.** In this subsection, we study the dependence of the convergence constants $C_\theta$ and $c_\theta$ (appearing in Theorem 4.9) with respect to the viscosity parameter $\sigma$ and the Lipschitz constant $L_f^c$ of the coupling term $f^c$. First, let us recall the constant $c(\sigma, \theta, M)$, introduced in (4.14),

$$c(\sigma, \theta, M) = 1 + \frac{M^2}{\sigma(2\theta - 1)} \exp \left( \frac{M^2}{\sigma(2\theta - 1)} \right).$$

It is not difficult to see that $c(\sigma, \theta, M)$ decreases and converges to 1 as $\sigma$ goes to $+\infty$. The constant $E_2$ is independent of $\sigma$ and $L_f^c$ by its definition (assuming that the change of $L_f^c$ has no impact on the semi-concavity constant of $f^c$).

By the proofs in the previous subsection, we can give the following explicit formulas of $C_\theta$ and $c_\theta$ in Theorem 4.9 (without using $E_1$ and $E_3$):

$$C_{\theta} = 4c(\sigma, \theta, M)L_{f}^{c}\left(\|m_{0}^{c}\|_{L^{\infty}}^{2} + \frac{1}{2}\|\nabla m_{0}^{c}\|_{L^{\infty}}^{2}\right), \quad c_{\theta} = \max\left\{1 - \frac{\sigma\alpha^{c}(2\theta - 1)}{4c(\sigma, \theta, M)L_{f}^{c}E_{2}}, \frac{1}{2}\right\}. \quad (4.31)$$

**Lemma 4.23.** For the constants in (4.31), we have the following.

(1) Fix $L_f^c$ and $\theta$. There exists $\sigma^* > 0$ and $C_1^* > 0$, such that for any $\sigma \ge \sigma^*$, we have

$$C_{\theta} \le C_{1}^{*} \left( \|m_{0}^{c}\|_{L^{\infty}}^{2} + \frac{1}{2} \|\nabla m_{0}^{c}\|_{L^{\infty}}^{2} \right), \quad c_{\theta} = \frac{1}{2}.$$

(2) Fix $\sigma$ and $\theta$. There exists $L^* > 0$ and $C_2^* > 0$, such that for any $L_f^c \le L^*$, we have

$$C_{\theta} \le C_{2}^{*} \left( \|m_{0}^{c}\|_{L^{\infty}}^{2} + \frac{1}{2} \|\nabla m_{0}^{c}\|_{L^{\infty}}^{2} \right), \quad c_{\theta} = \frac{1}{2}.$$

*Proof.* Point (1) follows from the monotonicity of $c(\sigma, \theta, M)$ w.r.t. $\sigma$ and the fact that $M$ is independent of $\sigma$. We can take $C_1^* = 5L_f^c$ for example. To prove (2), we first notice that if $L_f^c \le 1$, then $M \le M^*$, where $M^*$ is defined by (4.2), replacing $L_f^c$ with 1. The monotonicity of $c(\sigma, \theta, M)$ w.r.t. $M$ shows that $c(\sigma, \theta, M) \le c(\sigma, \theta, M^*)$. Since $c(\sigma, \theta, M^*)L_f^c$ goes to 0 as $L_f^c$ goes to 0, we prove the existence of $L^*$, and $C_2^* = 4c(\sigma, \theta, M^*)L^*$. ■
From the proof of Proposition 3.1, we know that $c_\theta = 1/2$ implies that $\lambda_k = 1$ for any $k \ge 0$. In other words, Algorithm 2 is equivalent to the so-called best-response iteration, i.e., for any $k \ge 0$,

$$ (m^{k+1}, w^{k+1}) = \mathbf{BR}_{\theta}(m^k). $$

Combined with Lemma 4.23, we have the following observations.

(1) *High-viscosity case:* let $L_f^c$ be fixed, if $\sigma$ is large enough, then the best response iteration has a linear convergence rate with a factor $1/2$.

(2) *Weak-coupling case:* let $\sigma$ be fixed, if $L_f^c$ is small enough, then the best response iteration has a linear convergence rate with a factor $1/2$.

## 5. NUMERICAL TESTS

### 5.1. Problem formulation.
In this section, we consider an example of (MFG) in dimension one. We identify the torus with the segment $[0, 1]$. The initial distribution is concentrated around the point 0.5, the running cost is a quadratic function of the control, and the terminal cost $g(x)$ decreases to 0 as $x$ goes to zero. Additionally, we consider a non-local congestion term which penalizes the density of the agents within the intervals $[0.2, 0.3]$ and $[0.7, 0.8]$. We will refer to $[0.2, 0.3] \cup [0.7, 0.8]$ as the *congestion-sensitive zone*.

To model this situation, let us introduce the functions $\varphi_{A,k} \in C^\infty(\mathbb{R})$ and $\phi_{A,k,l_1,l_2} \in C^\infty(\mathbb{R})$, parameterized by $A > 0, k > 0, 0 < l_1 < l_2 < 1$ and defined by

$$ \varphi_{A,k}(x) = \begin{cases} Ae^{-\frac{1}{1-(kx)^2}}, & \text{if } |x| < \frac{1}{k}, \\ 0 & \text{otherwise,} \end{cases} \quad \phi_{A,k,l_1,l_2}(x) = \begin{cases} \varphi_{A,k}(x-l_1), & \text{if } x < l_1, \\ Ae^{-1} & \text{if } l_1 \le x \le l_2, \\ \varphi_{A,k}(x-l_2), & \text{otherwise.} \end{cases} $$

The function $\varphi_{A,k}$ is a smooth approximation of the piecewise constant function equal to $Ae^{-1}$ on $[x - 1/k, x + 1/k]$ and zero elsewhere. The function $\phi_{A,k,l_1,l_2}$ is a smooth approximation of the piecewise constant function equal to $Ae^{-1}$ on $[l_1 - 1/k, l_2 + 1/k]$ and zero elsewhere.

The data of our one-dimensional MFG is parameterized by five positive numbers $a_1, a_2, k_0, k_1$, and $k_2$ and defined by: For any $(t, x) \in Q$, any $v \in \mathbb{R}$, and any $m \in D(\mathbb{T})$,

* $\ell^c(t, x, v) = \frac{1}{2}v^2$;

* $m_0^c(x) = \phi_{1,k_0,0.49,0.51}(x)/\|\phi_{1,k_0,0.49,0.51}\|_{\mathbb{L}^1}$;

* $g^c(x) = \phi_{a_1,k_1,1/k_1,1-1/k_1}(x)$;

* $f^c(t, x, m) = h^c(x) \int_0^1 h^c(y)m(y)dy$, where $h^c(x) = \phi_{a_2,k_2,0.24,0.25}(x) + \phi_{a_2,k_2,0.75,0.76}(x)$.

We take $a_1 = 2, a_2 = 20, k_0 = 10, k_1 = 3$, and $k_2 = 20$. The functions $m_0^c, g^c$, and $h^c$ are shown in Figure 1. Moreover, we fix the viscosity coefficient $\sigma = 0.02$. In this context, the agents have their initial condition around 0.5. The terminal cost $g^c$ is an incentive to move to the point 0 (which is the same as the point 1, since we identify the torus with $[0, 1]$). If there was no congestion cost and no diffusion coefficient, the agents would move at a constant speed because of the quadratic running cost. The congestion cost $f^c$ is an incentive to spend less time in the congestion-sensitive area.

We can verify that this one-dimensional MFG satisfies Assumptions A, C, D, with the constants in Assumption A satisfying

$$ \alpha^c = 1, \quad L_{\ell}^c = 0, \quad L_g^c \le a_1 k_1, \quad L_f^c \le \frac{a_2^2 k_2}{e}. $$

Furthermore, [11, Assumption C, Appx. B] holds true for this example, which implies Assumption B for any $r < 1$ by [11, Thm. B.2].
<figure><img src="image_8.png" /><figcaption>FIGURE 1. Data of the one-dimensional MFG with $a_1 = 2$, $a_2 = 20$, $k_0 = 10$, $k_1 = 3$, and $k_2 = 20$.</figcaption></figure>

5.2. **Results.** For the discretization of the system, we first choose the parameters

$$ \theta = 0.8, \quad h = 1/300, \quad \text{and} \quad \Delta t = \frac{h^2}{2(1-\theta)\sigma} = 1/720, $$

and we present the outcome of Algorithm (Theta-mfg) after 1000 iterations of the GFW algorithm, for step-sizes determined by line-search. For a better interpretation of the result, we also present the solution of the problem obtained by removing the congestion term $f^c$, which is a simple stochastic optimal control problem that can be solved in one iteration of the GFW method.

We present the equilibrium distribution of the agents in Figure 2a (without congestion term on the left, with congestion on the right). Note that the vertical axis corresponds to the time variable and is oriented downwards. We also present the restriction of the equilibrium distribution to the time interval [0.2, 0.8] in Figure 2b, with another color scale. As the time progresses, the agents are transported towards the target points 0 and 1. The congestion term leads to a reduced density in the congestion-sensitive zone: We see two dark blue vertical areas corresponding to this zone. We also see that at time $t \approx 0.35$, a significant part of the agents is still located around 0.5 and has not crossed yet the sensitive zone, in comparison with the case without $f^c$. Similarly, we present the optimal control $v$ in Figure 3 (without congestion term on the left, with congestion term on the right). Unsurprisingly, the agents must have a high velocity (in absolute value) in the sensitive zone. It is interesting to see that for $t$ close to zero and for the agents not that close to 0.5, there is an incentive to “rush” to the sensitive zone. Finally, we display the value functions for the two problems in Figure 3. In the present setting, note that the optimal control is the discrete gradient of the value function.

We next investigate the convergence of Algorithm 2 (for the same discretization parameters as above). We execute Algorithm 2 with 1000 iterations, utilizing the open-loop choice $\lambda_k = 2/(k+2)$ and the closed-loop choice (3.6) (referred to as the line-search method). We present the convergence results in Figure 4. Evaluating $\gamma_k$, equal to $\mathcal{J}(m_k, w_k) - \mathcal{J}(\bar{m}, \bar{w})$ by definition, is difficult since the exact solution $(\bar{m}, \bar{w})$ is not known. On the other hand, the quantity $\bar{\gamma}_k$, which serves as an upper bound of $\gamma_k$ by (3.11) can directly computed in view of its definition, based on $(m_k, w_k)$ and $(\bar{m}_k, \bar{w}_k)$. Therefore, instead of evaluating $\gamma_k$, we display the evolution of $\bar{\gamma}_k$, see Figure 4. The two figures of Figure 4 are the same, with different scales for the horizontal axis. In the left part of Figure 4, we see that Algorithm 2 exhibits a convergence rate of order $\mathcal{O}(1/k^4)$ for the choice $\lambda_k = 2/(k+2)$, which is better than the theoretical convergence rate $\mathcal{O}(1/k)$ obtained from (3.5).
<figure><img src="image_9.png" /><figcaption>(A) Comparison of distributions in the time horizon [0, 1]: the case without $f^c$ (left), the case with $f^c$ (right).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(B) Comparison of distributions in the time horizon [0.2, 0.8]: the case without $f^c$ (left), the case with $f^c$ (right).</figcaption></figure>

FIGURE 2. Distributions

In the right part of Figure 4, a linear convergence rate can be observed for the line-search case, as predicted in (3.7).

Finally, we present numerical results concerning the mesh-independence of Algorithm 2 applied to (Theta-mfg). To see this, we discretize the state space with steps sizes: $h = 1/250$, $h = 1/500$, and $h = 1/1000$. The corresponding step sizes for the time space are: $\Delta t = 1/500$, $\Delta t = 1/2000$, and $\Delta t = 1/8000$. The convergence results associated with these discretization steps are displayed in Figure 5. From the left part of Figure 5, it can be observed that the convergence rate of Algorithm 2 remains unaffected by the choice of $h$ when $\lambda_k = 2/(k+2)$. The right part of Figure 5 shows that the convergence rate of Algorithm 2 can even benefit from a refinement of the discretization
<figure><img src="image_11.png" /><figcaption>(A) Comparison of optimal controls: the case without $f^c$ (left), the case with $f^c$ (right).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(B) Comparison of value functions: the case without $f^c$ (left), the case with $f^c$ (right).</figcaption></figure>

<caption>FIGURE 3. Optimal controls and value functions</caption>

parameters in the line-search case. These results are consistent with mesh-independence properties outlined in Theorems 4.5 and 4.9.

## 6. CONCLUSION

We have established mesh-independent convergence results for the resolution of potential MFGs with the generalized Frank-Wolfe algorithm. This robustness property makes the GFW algorithm a method of choice for the resolution of potential MFGs. Our analysis has benefited from the intrinsic simplicity of the convergence proof of the GFW algorithm, which at the discrete level only
<figure><img src="image_13.png" /><figcaption>FIGURE 4. Convergence results of Algorithm 2.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 5. Mesh-independence property of Algorithm 2.</figcaption></figure>

required us to prove natural $L^2$ energy estimates for the Fokker-Planck equation (in the sublinear case) and an estimate of the semi-concavity modulus of the HJB equation (in the linear case). We expect that our analysis can be extended to the combination of the GFW algorithm with other discretization schemes, such as the implicit scheme of [1]. Let us stress, however, that the implementation of the GFW algorithm is made difficult for such schemes, since they would imply to solve fully implicit discrete HJB equations, involving nonlinear implicit relations, while for the theta-scheme, the implicit relations are linear, thus much easier to handle.

Let us underline again that the application of the GFW algorithm requires us to interpret the MFG system as the first-order necessary optimality condition of an optimization problem, it is therefore restricted to the case of potential problems. The convexity of the potential problem is crucial in the current analysis, however, we mention that the Frank-Wolfe algorithm was extended to non-convex problems (see [28]). In a non-convex setting, the convergence to a global solution cannot be ensured, yet the convergence of some stationarity criterion can be demonstrated. Future work will aim at proving mesh-independent principles for those convergence results, in the context of non-convex mean-field-type optimal control problems.

Another line of research could focus on the extension of the current framework to the case of degenerate potential mean-field games. It would be natural to investigate the combination of the
Frank-Wolfe algorithm with semi-Lagrangian schemes, which can handle models with a possibly degenerate diffusion, see [22] and the references therein.

Finally, we mention the case of first-order MFGs in their Lagrangian formulation, in which the equilibrium configuration is described by a probability measure on some trajectory set. It has a prescribed marginal $m_0$, describing the initial conditions of the agents. In our preprint [31], we propose a general class of optimization problems, containing the potential Lagrangian MFGs. We propose a tractable variant of the Frank-Wolfe algorithm, based on our work [10], that is combined with a discretization of $m_0$ as an empirical distribution associated with a set of $N$ points. Our method exhibits a sublinear rate of convergence which can be qualified as mesh-independent, since it improves as $N$ increases.

## REFERENCES

[1] Y. Achdou, F. Camilli, and I. Capuzzo-Dolcetta. Mean field games: Convergence of a finite difference method. *SIAM Journal on Numerical Analysis*, 51(5):2585-2612, 2013.

[2] Y. Achdou and I. Capuzzo-Dolcetta. Mean field games: Numerical methods. *SIAM Journal on Numerical Analysis*, 48(3):1136-1162, 2010.

[3] Y. Achdou and M. Laurière. Mean field games and applications: Numerical aspects. *Mean Field Games: Cetraro, Italy 2019*, pages 249-307, 2020.

[4] Y. Achdou and A. Porretta. Convergence of a finite difference scheme to weak solutions of the system of partial differential equations arising in mean field games. *SIAM Journal on Numerical Analysis*, 54(1):161-186, 2016.

[5] E.L. Allgower, K. Böhmer, F.A. Potra, and W.C. Rheinboldt. A mesh-independence principle for operator equations and their discretizations. *SIAM Journal on Numerical Analysis*, 23(1):160-169, 1986.

[6] R. Andreev. Preconditioning the augmented lagrangian method for instationary mean field games with diffusion. *SIAM Journal on Scientific Computing*, 39(6):A2763-A2783, 2017.

[7] J.-D. Benamou and G. Carlier. Augmented lagrangian methods for transport optimization, mean field games and degenerate elliptic equations. *Journal of Optimization Theory and Applications*, 167(1):1-26, 2015.

[8] J.-D. Benamou, G. Carlier, and F. Santambrogio. Variational mean field games. In *Active Particles, Volume 1*, pages 141-171. Springer, 2017.

[9] J.F. Bonnans, P. Lavigne, and L. Pfeiffer. Discrete potential mean field games: duality and numerical resolution. *Mathematical Programming*, 202:241-278, 2023.

[10] J.F. Bonnans, K. Liu, N. Oudjane, L. Pfeiffer, and C. Wan. Large-scale nonconvex optimization: randomization, gap estimation, and numerical resolution. *arXiv preprint arXiv:2204.02366*, 2022.

[11] J.F. Bonnans, K. Liu, and L. Pfeiffer. Error estimates of a theta-scheme for second-order mean field games. *ESAIM: Mathematical Modelling and Numerical Analysis*, 57(4):2493-2528, 2023.

[12] K. Bredies, D. A Lorenz, and P. Maass. A generalized conditional gradient method and its connection to an iterative shrinkage method. *Computational Optimization and applications*, 42(2):173-193, 2009.

[13] P. Cannarsa and C. Sinestrari. *Semiconcave functions, Hamilton-Jacobi equations, and optimal control*, volume 58. Springer Science & Business Media, 2004.

[14] P. Cardaliaguet and S. Hadikhanloo. Learning in mean field games: the fictitious play. *ESAIM: Control, Optimization and Calculus of Variations*, 23(2):569-591, 2017.

[15] P. Cardaliaguet and C.-A. Lehalle. Mean field game of controls and an application to trade crowding. *Mathematics and Financial Economics*, 12(3):335-363, 2018.

[16] E. Carlini and F.J. Silva. A fully discrete semi-lagrangian scheme for a first order mean field game problem. *SIAM Journal on Numerical Analysis*, 52(1):45-67, 2014.

[17] E. Carlini and F.J. Silva. A semi-lagrangian scheme for a degenerate second order mean field game system. *Discrete & Continuous Dynamical Systems*, 35(9):4269, 2015.

[18] D.S. Clark. Short proof of a discrete Gronwall inequality. *Discrete applied mathematics*, 16(3):279-281, 1987.

[19] P.L. Combettes. Perspective functions: Properties, constructions, and examples. *Set-Valued and Variational Analysis*, 26(2):247-264, 2018.

[20] B. Djehiche, A. Tcheukam, and H. Tembine. Mean-field-type games in engineering. *AIMS Electronics and Electrical Engineering*, 1(1):18-37, 2017.

[21] M. Geist, J. Pérolat, M. Laurière, R. Elie, S. Perrin, O. Bachem, R. Munos, and O. Pietquin. Concave utility reinforcement learning: The mean-field game viewpoint. In *Proceedings of the 21st International Conference on Autonomous Agents and Multiagent Systems*, pages 489-497, 2022.

[22] J. Gianatti and F.J. Silva. Approximation of deterministic mean field games with control-affine dynamics. *Foundations of Computational Mathematics*, pages 1-45, 2023.
[23] S. Hadikhanloo and F.J. Silva. Finite mean field games: fictitious play and convergence to a first order continuous mean field game. *Journal de Mathématiques Pures et Appliquées*, 132:369–397, 2019.

[24] M. Huang, R.P. Malhamé, and P.E. Caines. Large population stochastic dynamic games: closed-loop Mckean-Vlasov systems and the Nash certainty equivalence principle. *Communications in Information & Systems*, 6(3):221–252, 2006.

[25] M. Jaggi. Revisiting Frank-Wolfe: Projection-free sparse convex optimization. In *International Conference on Machine Learning*, pages 427–435. PMLR, 2013.

[26] K. Kunisch and D. Walter. On fast convergence rates for generalized conditional gradient methods with backtracking stepsize. *Numerical Algebra, Control and Optimization*, pages 0–0, 2022.

[27] A. Lachapelle, J. Salomon, and G. Turinici. Computation of mean field equilibria in economics. *Mathematical Models and Methods in Applied Sciences*, 20(04):567–588, 2010.

[28] S. Lacoste-Julien. Convergence rate of Frank-Wolfe for non-convex objectives. *arXiv preprint arXiv:1607.00345*, 2016.

[29] P.-L. Lasry, J.-M. and Lions. Mean field games. *Japanese journal of mathematics*, 2(1):229–260, 2007.

[30] P. Lavigne and L. Pfeiffer. Generalized conditional gradient and learning in potential mean field games. *Applied Mathematics and Optimization*, 88, article 89, 2023.

[31] K. Liu and L. Pfeiffer. Mean field optimization problems: stability results and lagrangian discretization. *ArXiv preprint*, 2023.

[32] S. Perrin, J. Pérolat, M. Laurière, M. Geist, R. Elie, and O. Pietquin. Fictitious play for mean field games: Continuous time analysis and applications. *Advances in Neural Information Processing Systems*, 33:13199–13213, 2020.

<footnote n="1">¹UNIVERSITÉ PARIS-SACLAY, CNRS, CENTRALESUPÉLEC, INRIA, LABORATOIRE DES SIGNAUX ET SYSTÈMES, 91190, GIF-SUR-YVETTE, FRANCE.</footnote>

<footnote n="2">²INSTITUT POLYTECHNIQUE DE PARIS, CNRS, ECOLE POLYTECHNIQUE, CMAP, 91120 PALAISEAU, FRANCE.
Email address: kang.liu@polytechnique.edu, laurent.pfeiffer@inria.fr</footnote>