15.0pt+\int_{0}^{T}\Bigl{\\{}\Bigl{(}QX_{\tau_{p}(t)}+f\bigl{(}\overline{m}^{n}_{\tau_{p}(t)}\bigr{)}\Bigr{)}\cdot
Q\partial
X_{\tau_{t}(p)}+\alpha_{\tau_{p}(t)}\cdot\beta_{\tau_{p}(t)}\Bigr{\\}}\mathrm{d}t\biggr{]}.\end{split}$

We now solve the BSDE (under ${\mathbb{P}}^{{\bm{h}}^{n}}$)

$\mathrm{d}Y_{t}=-\bigl{(}Q^{\dagger}QX_{\tau_{p}(t)}+Q^{\dagger}f(\overline{m}^{n}_{\tau_{p}(t)})\bigr{)}\mathrm{d}t+Z_{t}^{\bm{B}}\mathrm{d}B_{t}+Z_{t}^{\bm{W}}\mathrm{d}W_{t}^{{\bm{h}}^{n}},\quad
Y_{T}=R^{\dagger}RX_{T}+R^{\dagger}g\bigl{(}\overline{m}^{n}_{T}\bigr{)}.$

By discrete integration by parts, we have

$\begin{split}&{\mathbb{E}}^{{\bm{h}}^{n}}\bigl{[}Y_{T}\cdot\partial
X_{T}\bigr{]}\\\
&=\sum_{\ell=0}^{N-1}{\mathbb{E}}^{{\bm{h}}^{n}}\Bigl{[}\Bigl{(}Y_{(\ell+1)T/p}-Y_{\ell
T/p}\Bigr{)}\cdot\partial X_{\ell T/p}+Y_{(\ell+1)T/p}\cdot\Bigl{(}\partial
X_{(\ell+1)T/p}-\partial X_{\ell T/p}\Bigr{)}\Bigr{]}\\\
&=\frac{T}{p}\sum_{\ell=0}^{N-1}{\mathbb{E}}^{{\bm{h}}^{n}}\Bigl{[}-\Bigl{(}Q^{\dagger}QX_{\ell
T/p}+Q^{\dagger}f(\overline{m}^{n}_{\ell T/p})\Bigr{)}\cdot\partial X_{\ell
T/p}+Y_{(\ell+1)T/p}\cdot\beta_{\ell T/p}\Bigr{]}.\end{split}$

Rewriting the above sum as an integral, we obtain

$\frac{\mathrm{d}}{\mathrm{d}\delta}\widetilde{J}^{p}({\bm{\alpha}}+\delta{\bm{\beta}})_{|\delta=0}={\mathbb{E}}^{{\bm{h}}^{n}}\biggl{[}\int_{0}^{T}\Bigl{(}Y_{\tau_{p}(t)+T/p}+\alpha_{\tau_{p}(t)}\Bigr{)}\cdot\beta_{\tau_{p}(t)}\mathrm{d}t\biggr{]},$

which means that the optimizer is given by

$\alpha_{\ell
T/p}^{n+1}=-{\mathbb{E}}^{{\bm{h}}^{n}}\Bigl{[}Y_{(\ell+1)T/p}^{n+1}\,|\,{\mathcal{F}}_{\ell
T/p}^{p,X_{0},{\bm{B}},{\bm{W}}}\Bigr{]}=-Y_{\ell
T/p}+\frac{T}{p}\Bigl{(}Q^{\dagger}QX_{\ell
T/p}^{n+1}+Q^{\dagger}f(\overline{m}^{n}_{\ell T/p})\Bigr{)},$ (2.41)

for $\ell\in\\{0,\cdots,p-1\\}$, where

$\begin{split}&\mathrm{d}X_{t}^{n+1}=\alpha_{t}^{n+1}\mathrm{d}t+\sigma\mathrm{d}B_{t}+\varepsilon\mathrm{d}W_{t}^{p,{\bm{h}}^{n}},\\\
&\mathrm{d}Y_{t}^{n+1}=-\bigl{(}Q^{\dagger}QX_{\tau_{p}(t)}^{n+1}+Q^{\dagger}f(\overline{m}_{\tau_{p}(t)}^{n})\bigr{)}\mathrm{d}t+Z_{t}^{n+1,{\bm{B}}}\mathrm{d}B_{t}+Z_{t}^{n+1,{\bm{W}}}\mathrm{d}W_{t}^{{\bm{h}}^{n}},\quad
t\in[0,T],\\\
&Y_{T}^{n+1}=R^{\dagger}RX_{T}^{n+1}+R^{\dagger}g(\overline{m}^{n}_{T}).\end{split}$

Importantly, we stress that ${\bm{X}}_{\tau_{p}(\cdot)}^{n+1}$ is
$({\mathbb{F}}^{p,X_{0},{\bm{B}},{\bm{W}}})_{0\leq t\leq T}$ adapted. We then
let, for the same solution ${\bm{\eta}}^{(p)}$ as in (2.24)–(2.26),

$\begin{split}\widetilde{h}_{\ell
T/p}^{n+1}&={\mathbb{E}}^{{\bm{h}}^{n}}\Bigl{[}Y_{(\ell+1)T/p}^{n+1}\,|\,{\mathcal{F}}_{\ell
T/p}^{X_{0},{\bm{B}},{\bm{W}}}\Bigr{]}-\eta_{\ell T/p}^{(p)}X^{n+1}_{\ell
T/p}\\\ &=Y_{\ell T/p}^{n+1}-\frac{T}{p}\Bigl{(}Q^{\dagger}QX_{\ell
T/p}^{n+1}+Q^{\dagger}f(\overline{m}^{n}_{\ell T/p})\Bigr{)}-\eta_{\ell
T/p}^{(p)}X^{n+1}_{\ell T/p},\end{split}$ (2.42)

for $\ell\in\\{0,\cdots,p-1\\}$. We obtain, for $\ell\in\\{0,\cdots,p-2\\}$,

$\begin{split}&\widetilde{h}_{(\ell+1)T/p}^{n+1}-\widetilde{h}_{\ell
T/p}^{n+1}=-\frac{T}{p}\Bigl{(}Q^{\dagger}QX_{(\ell+1)T/p}^{n+1}+Q^{\dagger}f\bigl{(}\overline{m}_{(\ell+1)T/p}^{n}\bigr{)}\Bigr{)}-\frac{T}{p}\eta_{\ell
T/p}^{(p)}\alpha^{n+1}_{\ell T/p}\\\ &\hskip
15.0pt-\Bigl{(}\eta_{(\ell+1)T/p}^{(p)}-\eta_{\ell
T/p}^{(p)}\Bigr{)}X_{(\ell+1)T/p}^{n+1}+\int_{\ell
T/p}^{(\ell+1)T/p}k_{s}^{n+1,{\bm{B}}}\mathrm{d}B_{s}+\int_{\ell
T/p}^{(\ell+1)T/p}k_{s}^{n+1,{\bm{W}}}\mathrm{d}W_{s}^{{\bm{h}}^{n}},\end{split}$

for some square integrable and
${\mathbb{F}}^{X_{0},{\bm{B}},{\bm{W}}}$-progressively measurable processes
${\bm{k}}^{n+1,{\bm{B}}}$ and ${\bm{k}}^{n+1,{\bm{W}}}$. Above, we used the
fact that
${\bm{W}}_{\tau_{p}(t)+T/p}^{p,{\bm{h}}^{n}}-{\bm{W}}_{\tau_{p}(t)}^{p,{\bm{h}}^{n}}={\bm{W}}_{\tau_{p}(t)+T/p}^{{\bm{h}}^{n}}-{\bm{W}}_{\tau_{p}(t)}^{{\bm{h}}^{n}}$.
Now, using (2.41) and conditioning on
${\mathcal{F}}_{\tau_{p}(t)}^{p,X_{0},{\bm{B}},{\bm{W}}}$ in (2.42), we obtain

$\begin{split}\alpha_{\ell T/p}^{n+1}=-\eta_{\ell T/p}^{(p)}X^{n+1}_{\ell
T/p}-{\mathbb{E}}^{{\bm{h}}^{n}}\Bigl{[}\widetilde{h}_{\ell
T/p}^{n+1}\,|\,{\mathcal{F}}_{\ell
T/p}^{p,X_{0},{\bm{B}},{\bm{W}}}\Bigr{]},\quad\ell\in\\{0,\cdots,p-1\\},\end{split}$

which permits to write

$\begin{split}X^{n+1}_{(\ell+1)T/p}&=X^{n+1}_{\ell
T/p}+\frac{T}{p}\alpha_{\ell T/p}^{n+1}+\sigma\Bigl{(}B_{(\ell+1)T/p}-B_{\ell
T/p}\Bigr{)}+\varepsilon\Bigl{(}W^{{\bm{h}}^{n}}_{(\ell+1)T/p}-W^{{\bm{h}}^{n}}_{\ell
T}\Bigr{)}\\\ &=\Bigl{(}I_{d}-\frac{T}{p}\eta_{\ell
T/p}^{(p)}\Bigr{)}X^{n+1}_{\ell
T/p}-\frac{T}{p}{\mathbb{E}}^{{\bm{h}}^{n}}\Bigl{[}\widetilde{h}_{\ell
T/p}^{n+1}\,|\,{\mathcal{F}}_{\ell T/p}^{p,X_{0},{\bm{B}},{\bm{W}}}\Bigr{]}\\\
&\hskip 15.0pt+\sigma\Bigl{(}B_{(\ell+1)T/p}-B_{\ell
T/p}\Bigr{)}+\varepsilon\Bigl{(}W^{{\bm{h}}^{n}}_{(\ell+1)T/p}-W^{{\bm{h}}^{n}}_{\ell
T}\Bigr{)}.\end{split}$

Modifying the processes ${\bm{k}}^{n+1,{\bm{B}}}$ and
${\bm{h}}^{n+1,{\bm{W}}}$ and using in addition the first equation for
${\bm{\eta}}^{(p)}$ in (2.24), we end-up with

$\begin{split}\widetilde{h}_{(\ell+1)T/p}^{n+1}-\widetilde{h}_{\ell
T/p}^{n+1}&=-\frac{T}{p}Q^{\dagger}f(\overline{m}_{(\ell+1)T/p}^{n})+\frac{T}{p}\Bigl{(}\eta_{\ell
T/p}^{(p)}+\rho_{\ell
T/p}^{(p)}\Bigr{)}{\mathbb{E}}^{{\bm{h}}^{n}}\Bigl{[}\widetilde{h}_{\ell
T/p}^{n+1}\,|\,{\mathcal{F}}_{\ell T/p}^{p,X_{0},{\bm{B}},{\bm{W}}}\Bigr{]}\\\
&\hskip 15.0pt+\int_{\ell
T/p}^{(\ell+1)T/p}k_{s}^{n+1,{\bm{B}}}\mathrm{d}B_{s}+\int_{\ell
T/p}^{(\ell+1)T/p}k_{s}^{n+1,{\bm{W}}}\mathrm{d}W_{s}^{{\bm{h}}^{n}},\end{split}$
(2.43)

for $\ell\in\\{0,\cdots,p-2\\}$, where

$\rho_{\ell
T/p}^{n+1}=\frac{T}{p}Q^{\dagger}Q+\Bigl{(}\eta_{(\ell+1)T/p}^{(p)}-\eta_{\ell
T/p}^{(p)}\Bigr{)}.$

By recalling that $\overline{m}_{(\ell+1)T/p}^{n}$ is
${\mathcal{F}}_{(\ell+1)T/p}^{p,{\bm{W}}}$-measurable and by modifying the
process ${\bm{k}}^{n+1,{\bm{W}}}$ accordingly, we can rewrite (2.43) in the
form

$\begin{split}\widetilde{h}_{(\ell+1)T/p}^{n+1}-\widetilde{h}_{\ell
T/p}^{n+1}&=-\frac{T}{p}{\mathbb{E}}^{{\bm{h}}^{n}}\Bigl{[}Q^{\dagger}f(\overline{m}_{(\ell+1)T/p}^{n})\,|\,{\mathcal{F}}_{\ell
T/p}^{p,{\bm{W}}}\Bigr{]}\\\ &\hskip 15.0pt+\frac{T}{p}\Bigl{(}\eta_{\ell
T/p}^{(p)}+\rho_{\ell
T/p}^{(p)}\Bigr{)}{\mathbb{E}}^{{\bm{h}}^{n}}\Bigl{[}\widetilde{h}_{\ell
T/p}^{n+1}|{\mathcal{F}}_{\ell T/p}^{p,X_{0},{\bm{B}},{\bm{W}}}\Bigr{]}\\\
&\hskip 15.0pt+\int_{\ell
T/p}^{(\ell+1)T/p}k_{s}^{n+1,{\bm{B}}}\mathrm{d}B_{s}+\int_{\ell
T/p}^{(\ell+1)T/p}k_{s}^{n+1,{\bm{W}}}\mathrm{d}W_{s}^{{\bm{h}}^{n}},\end{split}$
(2.44)

for $\ell\in\\{0,\cdots,p-2\\}$. When $\ell=p-1$, we deduce from (2.42) that

$\begin{split}\widetilde{h}^{n+1}_{(p-1)T/p}&=R^{\dagger}R\,{\mathbb{E}}^{{\bm{h}}^{n}}\bigl{[}X_{T}^{n+1}|{\mathcal{F}}_{(p-1)T/p}^{X_{0},{\bm{B}},{\bm{W}}}\bigr{]}+R^{\dagger}g\bigl{(}\overline{m}^{n}_{T}\bigr{)}-\eta_{(p-1)T/p}^{(p)}X^{n+1}_{(p-1)T/p}\\\
&=\Bigl{[}R^{\dagger}R\Bigl{(}I_{d}-\frac{T}{p}\eta_{(p-1)T/p}^{(p)}\Bigr{)}-\eta_{(p-1)T/p}^{(p)}\Bigr{]}X^{n+1}_{(p-1)T/p}\\\
&\hskip
15.0pt+R^{\dagger}g\bigl{(}\overline{m}^{n}_{T}\bigr{)}-\frac{T}{p}R^{\dagger}R\,{\mathbb{E}}^{{\bm{h}}^{n}}\Bigl{[}\widetilde{h}_{(p-1)T/p}^{n+1}|{\mathcal{F}}_{(p-1)T/p}^{p,X_{0},{\bm{B}},{\bm{W}}}\Bigr{]},\end{split}$

which yields, as a boundary condition,

$\widetilde{h}^{n+1}_{(p-1)T/p}=R^{\dagger}g\bigl{(}\overline{m}^{n}_{T}\bigr{)}-\frac{T}{p}R^{\dagger}R\,{\mathbb{E}}^{{\bm{h}}^{n}}\Bigl{[}\widetilde{h}_{(p-1)T/p}^{n+1}|{\mathcal{F}}_{(p-1)T/p}^{p,X_{0},{\bm{B}},{\bm{W}}}\Bigr{]}.$
(2.45)

We now recall that, by construction, $\widetilde{h}_{\ell T/p}^{n+1}$ is
${\mathcal{F}}_{\ell T/p}^{X_{0},{\bm{B}},{\bm{W}}}$-measurable. Arguing as in
Remark 2.13, we can prove inductively (over $\ell$) that $\widetilde{h}_{\ell
T/p}^{n+1}$ is in fact ${\mathcal{F}}_{\ell T/p}^{p,{\bm{W}}}$-measurable. As
a consequence, ${\bm{k}}^{n+1,{\bm{B}}}$ has to be zero. By combining (2.44)
and (2.45), we hence recover BSDE (2.28).

∎

###  2.3. Exploitation versus exploration

We now address the regret, when we play the equilibrium strategy given by the
fictitious play.

####  2.3.1. Approximate Nash equilibrium formed by the solution of the MFG
with common noise

We first prove (which is easier) that the solution of the MFG with common
noise forms an approximate Nash equilibrium of the game without common noise.
In order to do so, we call ${\bm{\alpha}}^{\star}$ the equilibrium strategy of
the MFG with an $\varepsilon$ common noise, as given by Remark 2.4. Then

${\mathbb{E}}^{\bm{h}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}}^{\star};{\bm{m}};\varepsilon{\bm{W}}^{\bm{h}}\bigr{)}\Bigr{]},$

with ${\bm{m}}$ and ${\bm{h}}$ as in the statement of Remark 2.4, is the value
of the game with the $\varepsilon$-common noise, the intensity $\varepsilon$
being kept fixed in $(0,1]$ throughout the subsection. In particular, it must
be clear that both ${\bm{m}}$ and ${\bm{h}}$ depend on $\varepsilon$. We
claim:

######  Proposition 2.20.

There exists a constant $C$, only depending on the parameters $d$, $T$,
$\|f\|_{1,\infty}$, $\|g\|_{1,\infty}$, $|Q|$ and $|R|$, such that

$\inf_{\bm{\alpha}}{\mathbb{E}}^{\bm{h}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}};{\bm{m}};0\bigr{)}\Bigr{]}\geq{\mathbb{E}}^{\bm{h}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}}^{\star};{\bm{m}};0\bigr{)}\Bigr{]}-C\varepsilon,$

the infimum in the left-hand side being taken over
${\mathbb{F}}^{X_{0},{\bm{B}},{\bm{W}}}$-progressively measurable and square-
integrable processes $({\bm{\alpha}}_{t})_{0\leq t\leq T}$.

We start with the following lemma:

######  Lemma 2.21.

There exists a constant $C$, only depending on the parameters $d$, $T$,
$\|f\|_{1,\infty}$, $\|g\|_{1,\infty}$, $|Q|$ and $|R|$, such that, for any
${\mathbb{F}}^{X_{0},{\bm{B}},{\bm{W}}}$-progressively measurable process
$({\bm{\alpha}}_{t})_{0\leq t\leq T}$,

$\begin{split}&\biggl{|}{\mathbb{E}}^{\bm{h}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}};{\bm{m}};\varepsilon{\bm{W}}^{\bm{h}}\bigr{)}\Bigr{]}-{\mathbb{E}}^{\bm{h}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}};{\bm{m}};0\bigr{)}\Bigr{]}\biggr{|}\leq
C\biggl{(}1+{\mathbb{E}}\bigl{[}|X_{0}|^{2}\bigr{]}^{1/2}+{\mathbb{E}}^{\bm{h}}\biggl{[}\int_{0}^{T}|\alpha_{t}|^{2}\mathrm{d}t\biggr{]}^{1/2}\biggr{)}\varepsilon.\end{split}$

The symbol ${\mathbb{E}}$ instead of ${\mathbb{E}}^{{\bm{h}}}$ for the $L^{2}$
moment of the initial condition is fully justified by the fact that $X_{0}$ is
${\mathcal{F}}_{0}$-measurable.

###### Proof.

For the same initial condition $X_{0}$, we let

$\begin{split}&\mathrm{d}X_{t}^{\varepsilon}=\alpha_{t}\mathrm{d}t+\sigma\mathrm{d}B_{t}+\varepsilon\mathrm{d}W_{t}^{\bm{h}},\quad
t\in[0,T].\end{split}$

In particular, we write ${\bm{X}}^{0}=(X_{t}^{0})_{0\leq t\leq T}$ when there
is no common noise. Then, we obviously have

$\sup_{0\leq t\leq T}|X_{t}-X_{t}^{0}|\leq C\varepsilon\sup_{0\leq t\leq
T}|W_{t}^{\bm{h}}|.$

Meanwhile, we observe that

${\mathbb{E}}^{\bm{h}}\Bigl{[}\sup_{0\leq t\leq
T}\bigl{(}|X_{t}|^{2}+|X_{t}^{0}|^{2}\bigr{)}\Bigr{]}\leq
C\biggl{(}1+{\mathbb{E}}\bigl{[}|X_{0}|^{2}\bigr{]}+{\mathbb{E}}^{\bm{h}}\biggl{[}\int_{0}^{T}|\alpha_{t}|^{2}\mathrm{d}t\biggr{]}\biggr{)}.$

The result easily follows. ∎

The next lemma says that we can easily restrict ourselves to controls with a
finite energy.

######  Lemma 2.22.

There exists a constant $A$, only depending on the parameters $d$, $Q$, $R$,
$T$, $\|f\|_{1,\infty}$ and $\|g\|_{1,\infty}$ such that, for any
$\varepsilon\in(0,1]$ and any ${\mathbb{F}}^{{\bm{B}},{\bm{W}}}$-progressively
measurable process ${\bm{\alpha}}=({\bm{\alpha}}_{t})_{0\leq t\leq T}$,

$\begin{split}&{\mathbb{E}}^{\bm{h}}\biggl{[}\int_{0}^{T}|\alpha_{t}|^{2}\mathrm{d}t\biggr{]}\geq
A\Rightarrow{\mathbb{E}}^{\bm{h}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}};{\bm{m}};0\bigr{)}\Bigr{]}\geq{\mathbb{E}}^{\bm{h}}\Bigl{[}{\mathcal{R}}\bigl{(}0;{\bm{m}};\varepsilon{\bm{W}}\bigr{)}\Bigr{]}+1.\end{split}$

###### Proof.

It suffices to observe that

${\mathbb{E}}^{\bm{h}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}};{\bm{m}};0\bigr{)}\Bigr{]}\geq\frac{1}{2}{\mathbb{E}}^{\bm{h}}\biggl{[}\int_{0}^{T}|\alpha_{t}|^{2}\mathrm{d}t\biggr{]}.$

Meanwhile,

${\mathbb{E}}^{\bm{h}}\Bigl{[}{\mathcal{R}}\bigl{(}0;{\bm{m}};\varepsilon{\bm{W}}\bigr{)}\Bigr{]}\leq
C.$

The proof is easily completed. ∎

We now complete the

###### Proof of Proposition 2.20..

By Lemma 2.22, we can reduce ourselves to the set ${\mathcal{E}}_{A}$ of
processes ${\bm{\alpha}}$ such that

${\mathbb{E}}^{\bm{h}}\biggl{[}\int_{0}^{T}|\alpha_{t}|^{2}\mathrm{d}t\biggr{]}\leq
A.$

We then apply Lemma 2.21, which says that

$\sup_{{\bm{\alpha}}\in{\mathcal{E}}_{A}}\Bigl{|}{\mathbb{E}}^{\bm{h}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}};{\bm{m}};\varepsilon{\bm{W}}^{\bm{h}}\bigr{)}\Bigr{]}-{\mathbb{E}}^{\bm{h}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}};{\bm{m}};0\bigr{)}\Bigr{]}\Bigr{|}\leq
C\varepsilon.$

The proof is easily completed, using the fact that, for any
${\bm{\alpha}}\in{\mathcal{E}}_{A}$,

$\begin{split}{\mathbb{E}}^{\bm{h}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}}^{\star};{\bm{m}};\varepsilon{\bm{W}}^{\bm{h}}\bigr{)}\Bigr{]}&\leq{\mathbb{E}}^{\bm{h}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}};{\bm{m}};\varepsilon{\bm{W}}^{\bm{h}}\bigr{)}\Bigr{]}\leq{\mathbb{E}}^{\bm{h}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}};{\bm{m}};0\bigr{)}\Bigr{]}+C\varepsilon.\end{split}$

∎

####  2.3.2. Approximate Nash equilibrium formed by the return of the
fictitious play

We now state a similar result, but for the approximation returned by the
fictitious play subjected to the $\varepsilon$-randomization.

The point is hence to consider, under the environment ${\bm{m}}^{n+1}$, the
strategy ${\bm{\alpha}}^{n+1}$ (returned by the fictitious play, see (2.23))
as an approximate Nash equilibrium under the probability
${\mathbb{P}}^{{\bm{h}}^{n}}$. Notice that this is fully legitimated by the
fact that

$dX^{n+1}_{t}=\alpha_{t}^{n+1}\mathrm{d}t+\sigma\mathrm{d}B_{t}+\varepsilon\mathrm{d}W_{t}^{p,{\bm{h}}^{n}},\quad
t\in[0,T],$

together with

$m^{n+1}_{t}={\mathbb{E}}^{{\bm{h}}^{n}}\Bigl{[}X_{t}^{n+1}\,|\,\sigma({\bm{W}})\Bigr{]},\quad
t\in[0,T].$

In other words, the conditional state under the candidate strategy is the
environment itself, which is a pre-requisite in mean field games. Notice that,
here as well, ${\bm{\alpha}}^{n+1}$, ${\bm{h}}^{n}$ and ${\bm{m}}^{n+1}$ do
depend on $\varepsilon$.

######  Theorem 2.23.

There exists a constant $C$, only depending on the parameters $d$, $T$,
$\|f\|_{1,\infty}$, $\|g\|_{1,\infty}$, $|Q|$ and $|R|$, such that

$\begin{split}\inf_{\bm{\alpha}}{\mathbb{E}}^{{\bm{h}}_{n}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}};{\bm{m}}^{n+1};0\bigr{)}\Bigr{]}&\geq{\mathbb{E}}^{{\bm{h}}_{n}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}}^{n+1};{\bm{m}}^{n+1};0\bigr{)}\Bigr{]}-C\varepsilon\\\
&\hskip
15.0pt-\frac{C}{\varepsilon}\biggl{(}\frac{1}{n}+\frac{(1+\ln(p))^{1/2}}{p^{1/2}}\biggr{)}\Bigl{[}1+\bigl{(}\ln(np\varepsilon)\bigr{)}_{+}\Bigr{]},\end{split}$

the infimum in the left-hand side being taken over
${\mathbb{F}}^{X_{0},{\bm{B}},{\bm{W}}}$-progressively measurable and square-
integrable processes $({\bm{\alpha}}_{t})_{0\leq t\leq T}$.

The reader must realize that, differently from what is done in (2.23), the
environment used in the cost functional is the conditional state of the
reference particle when using the strategy ${\bm{\alpha}}^{n+1}$. This looks
subtle, but this makes a big difference in the analysis.

Notice that Proposition 2.19 just gives a bound on the difference between
$\overline{\bm{m}}^{n+1}$ and ${\bm{m}}$, while we here need a bound on the
difference between ${\bm{m}}^{n+1}$ and ${\bm{m}}$. This prompts us to start
with the following refinement of Proposition 2.19.

######  Proposition 2.24.

With the same notation as in the statement of Theorem 2.14, there exists a
constant $C$, only depending on $d$, $T$, $\|f\|_{1,\infty}$ and
$\|g\|_{1,\infty}$ such that

${\mathbb{E}}^{{\bm{h}}}\Bigl{[}\sup_{0\leq t\leq
T}\Bigl{(}|X_{t}^{n}-X_{t}|^{2}+|m_{t}^{n}-m_{t}|^{2}\Bigr{)}\Bigr{]}+{\mathbb{E}}^{{\bm{h}}}\biggl{[}\int_{0}^{T}|\alpha_{t}^{n}-\alpha_{t}^{\star}|^{2}\mathrm{d}t\biggr{]}\leq
C\Bigl{(}\frac{1}{n^{2}}+\frac{1+\ln(p)}{p}\Bigr{)}.$

###### Proof of Proposition 2.24.

We recall from Lemma 2.11 and (2.31) that

$\mathrm{d}X^{n+1}_{t}=-\Bigl{(}\eta_{\tau_{p}(t)}^{(p)}X_{\tau_{p}(t)}^{n+1}+h_{t}^{n+1}\Bigr{)}\mathrm{d}t+\sigma\mathrm{d}B_{t}+\varepsilon\mathrm{d}W_{t}^{p,{\bm{h}}^{n}},$

which we rewrite

$\begin{split}\mathrm{d}X^{n+1}_{t}&=-\Bigl{(}\eta_{\tau_{p}(t)}^{(p)}X_{\tau_{p}(t)}^{n+1}+h^{n+1}_{t}\Bigr{)}\mathrm{d}t+\bigl{(}h_{t}^{n}-h_{t}\bigr{)}\mathrm{d}t\\\
&\hskip
15.0pt+\sigma\mathrm{d}B_{t}+\varepsilon\mathrm{d}\Bigl{(}W_{t}^{p,{\bm{h}}^{n}}-W_{t}^{{\bm{h}}^{n}}\Bigr{)}+\varepsilon\mathrm{d}W_{t}^{{\bm{h}}},\quad
t\in[0,T].\end{split}$

Meanwhile, we write

$\mathrm{d}X_{t}=-\bigl{(}\eta_{t}X_{t}+h_{t}\bigr{)}\mathrm{d}t+\sigma\mathrm{d}B_{t}+\varepsilon\mathrm{d}W_{t}^{{\bm{h}}},\quad
t\in[0,T].$

Following the proof of Lemma 2.18,

${\mathbb{E}}^{{\bm{h}}}\Bigl{[}\sup_{|t-s|\leq
1/p}|X_{t}-X_{s}|^{2}\Bigr{]}\leq C\frac{1+\ln(p)}{p}.$

We can rewrite the dynamics for ${\bm{X}}$ in the form

$\begin{split}\mathrm{d}X_{t}&=-\Bigl{(}\eta_{t}^{(p)}X_{\tau_{p}(t)}+h_{t}\Bigr{)}\mathrm{d}t+\Bigl{(}\eta_{t}^{(p)}X_{\tau_{p}(t)}-\eta_{t}X_{t}\Bigr{)}\mathrm{d}t+\sigma\mathrm{d}B_{t}+\varepsilon\mathrm{d}W_{t}^{{\bm{h}}},\quad
t\in[0,T].\end{split}$

And then, by Remark 2.12 and Lemma 2.18,

$\begin{split}{\mathbb{E}}^{{\bm{h}}}\Bigl{[}\sup_{0\leq r\leq
t}|X_{r}^{n+1}-X_{r}|^{2}\Bigr{]}&\leq
C\int_{0}^{t}{\mathbb{E}}^{{\bm{h}}}\Bigl{[}\sup_{0\leq r\leq
s}|X_{r}^{n+1}-X_{r}|^{2}\Bigr{]}\mathrm{d}r+C\frac{1+\ln(p)}{p}\\\ &\hskip
15.0pt+{\mathbb{E}}^{{\bm{h}}}\int_{0}^{T}|h_{t}^{n}-h_{t}|^{2}\mathrm{d}t+{\mathbb{E}}^{{\bm{h}}}\int_{0}^{T}|h_{t}^{n+1}-h_{t}|^{2}\mathrm{d}t.\end{split}$

By Proposition 2.19, the two terms on the last line are less than
$C/n^{2}+C(1+\ln(p))/p$. The result for ${\bm{X}}^{n+1}-{\bm{X}}$ then follows
from Gronwall’s lemma. The result for ${\bm{m}}^{n+1}-{\bm{m}}$ can be proved
in the same way. As for the bound for
${\bm{\alpha}}^{n+1}-{\bm{\alpha}}^{\star}$, it follows from the feedback form
of the control. ∎

We now turn to

###### Proof of Theorem 2.23.

First Step. Recalling the shape of the optimizer ${\bm{\alpha}}^{n+1}$ in
Lemma 2.11, together with the fact that $\widetilde{\bm{h}}^{n+1}$ is bounded,
uniformly in $n\geq 1$, see Lemma 2.16, we deduce that

$\Bigl{|}{\mathbb{E}}^{{\bm{h}}_{n}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}}^{n+1};{\bm{m}}^{n+1};0\bigr{)}\Bigr{]}\Bigr{|}\leq
C.$

Here, we must invoke the fact that the problem

$\inf_{\bm{\alpha}}{\mathbb{E}}^{{\bm{h}}_{n}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}};{\bm{m}}^{n+1};0\bigr{)}\Bigr{]}$

is a linear quadratic-problem in a random environment. Similar to the
minimization problem studied in Remark 2.4, its solution is given in feedback
form. If we call $\widehat{\bm{\alpha}}^{n+1}$ the minimizer, it reads

$\widehat{\alpha}_{t}^{n+1}=-\bigl{(}\eta_{t}\widehat{X}_{t}^{n+1}+\widehat{h}_{t}^{n+1}\bigr{)},\quad
t\in[0,T],$

with

$\mathrm{d}\widehat{X}_{t}^{n+1}=-\bigl{(}\eta_{t}\widehat{X}_{t}^{n+1}+\widehat{h}_{t}^{n+1}\bigr{)}\mathrm{d}t+\sigma\mathrm{d}B_{t},$

and where, ${\mathbb{P}}^{{\bm{h}}^{n}}$ almost surely,

$\forall t\in[0,T],\quad|\widehat{h}_{t}^{n+1}|\leq C.$

As a result, with probability 1 under ${\mathbb{P}}^{{\bm{h}}^{n}}$,

$\begin{split}&\sup_{0\leq t\leq T}|\widehat{X}_{t}^{n+1}|\leq
C\Bigl{(}1+\sup_{0\leq t\leq T}|B_{t}|\Bigr{)},\\\ &\sup_{0\leq t\leq
T}|\widehat{\alpha}_{t}^{n+1}|\leq C\Bigl{(}1+\sup_{0\leq t\leq
T}|B_{t}|\Bigr{)}.\end{split}$ (2.46)

In particular,

${\mathcal{R}}\bigl{(}\widehat{\bm{\alpha}}^{n+1};{\bm{m}}^{n+1};0\bigr{)}\leq
C\Bigl{(}1+\sup_{0\leq t\leq T}|B_{t}|^{2}\Bigr{)}.$ (2.47)

Second Step. As a result, for any $\varrho>0$,

$\begin{split}&\inf_{\bm{\alpha}}{\mathbb{E}}^{{\bm{h}}_{n}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}};{\bm{m}}^{n+1};0\bigr{)}\Bigr{]}={\mathbb{E}}^{{\bm{h}}_{n}}\Bigl{[}{\mathcal{R}}\bigl{(}\widehat{\bm{\alpha}}^{n+1};{\bm{m}}^{n+1};0\bigr{)}\Bigr{]}\\\
&\geq{\mathbb{E}}^{{\bm{h}}_{n}}\Bigl{[}{\mathcal{R}}\bigl{(}\widehat{\bm{\alpha}}^{n+1};{\bm{m}}^{n+1};0\bigr{)}{\mathbf{1}}_{\\{{\mathcal{R}}(\widehat{\bm{\alpha}}^{n+1};{\bm{m}}^{n+1};0)\leq\varrho\\}}\Bigr{]}.\end{split}$

Using the upper bound on $d_{\textrm{\rm
TV}}\bigl{(}{\mathbb{P}}^{\bm{h}},{\mathbb{P}}^{{\bm{h}}_{n}}\bigr{)}$ (see
Remark 2.15), we deduce that

$\begin{split}\inf_{\bm{\alpha}}{\mathbb{E}}^{{\bm{h}}_{n}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}};{\bm{m}}^{n+1};0\bigr{)}\Bigr{]}&\geq{\mathbb{E}}^{{\bm{h}}}\Bigl{[}{\mathcal{R}}\bigl{(}\widehat{\bm{\alpha}}^{n+1};{\bm{m}}^{n+1};0\bigr{)}{\mathbf{1}}_{\\{{\mathcal{R}}(\widehat{\bm{\alpha}}^{n+1};{\bm{m}}^{n+1};0)\leq\varrho\\}}\Bigr{]}\\\
&\hskip
15.0pt-C\frac{\varrho}{\varepsilon}\Bigl{(}\frac{1}{n}+\frac{(1+\ln(p))^{1/2}}{p^{1/2}}\Bigr{)}.\end{split}$

And then, by (2.47), we obtain

$\begin{split}\inf_{\bm{\alpha}}{\mathbb{E}}^{{\bm{h}}_{n}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}};{\bm{m}}^{n+1};0\bigr{)}\Bigr{]}&\geq{\mathbb{E}}^{{\bm{h}}}\Bigl{[}{\mathcal{R}}\bigl{(}\widehat{\bm{\alpha}}^{n+1};{\bm{m}}^{n+1};0\bigr{)}\Bigr{]}\\\
&\hskip
15.0pt-C\exp\bigl{(}-\eta\varrho\bigr{)}-C\frac{\varrho}{\varepsilon}\Bigl{(}\frac{1}{n}+\frac{(1+\ln(p))^{1/2}}{p^{1/2}}\Bigr{)},\end{split}$
(2.48)

for some $\eta>0$. Now, we notice from (2.46) that

${\mathbb{E}}^{{\bm{h}}}\Bigl{[}\sup_{0\leq t\leq
T}\bigl{|}\widehat{X}_{t}^{n+1}\bigr{|}^{2}\Bigr{]}\leq C.$

Therefore, by Proposition 2.24, we get

${\mathbb{E}}^{{\bm{h}}}\Bigl{[}\bigl{|}{\mathcal{R}}\bigl{(}\widehat{\bm{\alpha}}^{n+1};{\bm{m}}^{n+1};0\bigr{)}-{\mathcal{R}}\bigl{(}\widehat{\bm{\alpha}}^{n+1};{\bm{m}};0\bigr{)}\bigr{|}\Bigr{]}\leq
C\Bigl{(}\frac{1}{n}+\frac{(1+\ln(p))^{1/2}}{p^{1/2}}\Bigr{)}.$

Plugging the latter into (2.48), we obtain

$\begin{split}\inf_{\bm{\alpha}}{\mathbb{E}}^{{\bm{h}}_{n}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}};{\bm{m}}^{n+1};0\bigr{)}\Bigr{]}&\geq{\mathbb{E}}^{{\bm{h}}}\Bigl{[}{\mathcal{R}}\bigl{(}\widehat{\bm{\alpha}}^{n+1};{\bm{m}};0\bigr{)}\Bigr{]}\\\
&\hskip
15.0pt-C\exp\bigl{(}-\eta\varrho\bigr{)}-C\frac{\varrho}{\varepsilon}\Bigl{(}\frac{1}{n}+\frac{(1+\ln(p))^{1/2}}{p^{1/2}}\Bigr{)},\end{split}$

and then, by Proposition 2.20,

$\begin{split}&\inf_{\bm{\alpha}}{\mathbb{E}}^{{\bm{h}}_{n}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}};{\bm{m}}^{n+1};0\bigr{)}\Bigr{]}\\\
&\hskip
15.0pt\geq{\mathbb{E}}^{{\bm{h}}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}}^{\star};{\bm{m}};0\bigr{)}\Bigr{]}-C\varepsilon-C\exp\bigl{(}-\eta\varrho\bigr{)}-C\frac{\varrho}{\varepsilon}\Bigl{(}\frac{1}{n}+\frac{(1+\ln(p))^{1/2}}{p^{1/2}}\Bigr{)}.\end{split}$
(2.49)

Third Step. We now revert the computations achieved in the second step. To do
so, we observe that the bounds (2.46) and (2.47) remain true if we replace
$\widehat{\bm{X}}^{n+1}$ by ${\bm{X}}^{n+1}$ and similarly
$\widehat{\bm{\alpha}}^{n+1}$ by ${\bm{\alpha}}^{n+1}$, with ${\bm{X}}^{n+1}$
and ${\bm{\alpha}}^{n+1}$ being defined as in Lemma 2.11. Then, by Proposition
2.24, (2.49) yields

$\begin{split}&\inf_{\bm{\alpha}}{\mathbb{E}}^{{\bm{h}}_{n}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}};{\bm{m}}^{n+1};0\bigr{)}\Bigr{]}\\\
&\hskip
15.0pt\geq{\mathbb{E}}^{{\bm{h}}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}}^{n+1};{\bm{m}}^{n+1};0\bigr{)}\Bigr{]}-C\varepsilon-C\exp\bigl{(}-\eta\varrho^{2}\bigr{)}-C\frac{\varrho^{2}}{\varepsilon}\Bigl{(}\frac{1}{n}+\frac{(1+\ln(p))^{1/2}}{p^{1/2}}\Bigr{)}\\\
&\hskip
15.0pt\geq{\mathbb{E}}^{{\bm{h}}^{n}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}}^{n+1};{\bm{m}}^{n+1};0\bigr{)}{\mathbf{1}}_{\\{{\mathcal{R}}({\bm{\alpha}}^{n+1};{\bm{m}}^{n+1};0)\leq\varrho\\}}\Bigr{]}\\\
&\hskip
30.0pt-C\varepsilon-C\exp\bigl{(}-\eta\varrho^{2}\bigr{)}-C\frac{\varrho^{2}}{\varepsilon}\Bigl{(}\frac{1}{n}+\frac{(1+\ln(p))^{1/2}}{p^{1/2}}\Bigr{)},\end{split}$

and then,

$\begin{split}&\inf_{\bm{\alpha}}{\mathbb{E}}^{{\bm{h}}_{n}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}};{\bm{m}}^{n+1};0\bigr{)}\Bigr{]}\\\
&\hskip
15.0pt\geq{\mathbb{E}}^{{\bm{h}}^{n}}\Bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}}^{n+1};{\bm{m}}^{n+1};0\bigr{)}\Bigr{]}-C\varepsilon-C\exp\bigl{(}-\eta\varrho^{2}\bigr{)}-C\frac{\varrho}{\varepsilon}\Bigl{(}\frac{1}{n}+\frac{(1+\ln(p))^{1/2}}{p^{1/2}}\Bigr{)}.\end{split}$

Fourth Step. Choosing

$\eta\varrho=\biggl{(}-\ln\biggl{[}\frac{1}{\varepsilon}\biggl{(}\frac{1}{n}+\frac{(1+\ln(p))^{1/2}}{p^{1/2}}\biggr{)}\biggr{]}\biggr{)}_{+},$

and noticing that

$\eta\varrho\leq\biggl{(}-\ln\biggl{[}\frac{1}{\varepsilon}\biggl{(}\frac{1}{n}+\frac{1}{p^{1/2}}\biggr{)}\biggr{]}\biggr{)}_{+}\leq\biggl{(}\ln(\varepsilon)+\ln(n)+\ln(p)\biggr{)}_{+},$

we complete the proof. ∎

##  3\. Numerical experiments

This section is devoted to several numerical experiments based on our learning
algorithm. We first provide in Subsection 3.1 a theoretical analysis of a
benchmark example that we use throughout the section. In Subsection 3.2, we
present several variants of the implemented version of the algorithm. We also
explain how to compute a reference solution. Then, we study in Subsection 3.3
the numerical behavior of the algorithm for a fixed intensity $\varepsilon>0$
of the common noise. Finally, in Subsection 3.4, we discuss a numerical
vanishing viscosity method based on the learning method.

###  3.1. A benchmark example

As a particular instance of the cost functional $J$ in (1.2), we focus here on

$J\bigl{(}{\bm{\alpha}};{\bm{m}}\bigr{)}=\frac{1}{2}{\mathbb{E}}\biggl{[}\bigl{|}X_{T}+g(m_{T})\bigr{|}^{2}+\int_{0}^{T}|\alpha_{t}|^{2}\mathrm{d}t\biggr{]},$
(3.1)

which implicitly means that $f=0$. The examples that are addressed below are
in dimension $d=1$ or $d=2$, with $T=1$ in any cases. For simplicity, we also
work with $X_{0}=0$.

In dimension $d=1$, we choose

$g(x)=\cos(\kappa x),\quad x\in{\mathbb{R}},$ (3.2)

for a free positive parameter $\kappa$ that we tune from smaller to larger
values. Obviously, the Lipschitz constant of $g$ becomes larger with $\kappa$.
Accordingly, the coupling between the two forward and backward equations in
(1.7) becomes stronger as $\kappa$ gets larger, which makes it more difficult
to solve, especially when $\varepsilon=0$. We illustrate this in Lemma 3.1
below: it says that, when $\varepsilon=0$, (1.7) is uniquely solvable when
$|\kappa|<2$. Even more, the analysis performed in [24] shows that the more
standard (and obviously simpler) Picard scheme (1.9)–(1.10) would converge
when $\kappa$ is small, whether there is a common noise or not. Our result is
thus especially relevant when $\kappa$ gets larger.

In dimension $2$, we work with a similar terminal cost $g=(g_{1},g_{2})$:

$\begin{split}&g_{1}(x_{1},x_{2})=\cos(\kappa x_{1})\cos(\kappa x_{2}),\quad
g_{2}(x_{1},x_{2})=\sin(\kappa x_{1})\sin(\kappa x_{2}),\quad x_{1},\
x_{2}\in{\mathbb{R}},\end{split}$ (3.3)

where, as in dimension 1, $\kappa$ is a free positive parameter that we tune
form smaller to larger values.

Regardless of the choice of $g$, the Riccati equation (2.2) associated with
(3.1) writes

$\dot{\eta}_{t}-\eta_{t}^{2}=0,\quad t\in[0,1]\ ;\quad\eta_{1}=1,$ (3.4)

the solution of which is given by

$\eta_{t}=\frac{1}{2-t},\quad t\in[0,T].$ (3.5)

####  3.1.1. Solutions without common noise

Even though this is not needed for all the examples we handle below, we feel
useful to have a preliminary discussion about the shape of the solutions when
there is no common noise, keeping in mind that the MFG without common noise is
our original model. Recalling (1.7), we indeed have that, whenever there is no
common noise (i.e., $\varepsilon=0$), the equilibria are given as the
solutions of the deterministic system

$\dot{m}_{t}=-\bigl{(}\eta_{t}m_{t}+h_{t}\bigr{)},\quad\dot{h}_{t}=\eta_{t}h_{t},\quad
t\in[0,1]\ ;\quad h_{1}=g(m_{1}),$ (3.6)

with $m_{0}=0$, which prompts us to perform the changes of variable:

$\widetilde{m}_{t}=\frac{m_{t}}{2-t},\quad\widetilde{h}_{t}=(2-t)h_{t},\quad
t\in[0,1].$ (3.7)

We then have that $(m_{t},h_{t})_{0\leq t\leq 1}$ solves (3.6) if and only if

$\begin{split}&\dot{\widetilde{m}}_{t}=-\frac{1}{2-t}h_{t}=-\frac{1}{(2-t)^{2}}\widetilde{h}_{t},\quad\dot{\widetilde{h}}_{t}=0,\quad
t\in[0,1]\
;\quad\widetilde{h}_{1}=g\bigl{(}\widetilde{m}_{1}\bigr{)},\end{split}$

from which we deduce the following simpler characterization (recalling that
$m_{0}=0$)

$\widetilde{m}_{1}=-g\bigl{(}\widetilde{m}_{1}\bigr{)}\int_{0}^{1}\frac{\mathrm{d}t}{(2-t)^{2}}=-\frac{g(\widetilde{m}_{1})}{2}.$
(3.8)

There are as many equilibria as solutions to the equation $x=-g(x)/2$. We thus
have the following obvious lemma:

######  Lemma 3.1.

When $\varepsilon=0$ and regardless of the dimension, the solutions of the MFG
associated with the cost functional (3.1) are given by the roots
$\widetilde{m}_{1}$ of the equation $2x+g(x)=0$ and then by the changes of
variable (3.7). In particular, if the Lipschitz constant of the function $g$
is strictly less than 2, then the MFG has one and only one solution.

####  3.1.2. Potential structure when $\varepsilon=0$

Following (1.13) and (1.14), we may associate a mean field control problem
with the MFG (with $\varepsilon=0$) if the function $g$ derives from a
potential $G$. The cost functional is given by

${\mathcal{J}}({\bm{\alpha}})={\mathbb{E}}\biggl{[}\frac{1}{2}|X_{1}|^{2}+G\bigl{(}{\mathbb{E}}(X_{1})\bigr{)}+\frac{1}{2}\int_{0}^{1}|\alpha_{t}|^{2}\mathrm{d}t\biggr{]},$

where, in the right-hand side,

$X_{1}=\int_{0}^{1}\alpha_{t}\mathrm{d}t.$

The analysis of the minimization problem
$\inf_{{\bm{\alpha}}}{\mathcal{J}}({\bm{\alpha}})$ is straightforward. Indeed,
by an obvious convexity argument, we observe that

${\mathcal{J}}({\bm{\alpha}})\geq{\mathcal{J}}\biggl{(}{\mathbb{E}}\int_{0}^{1}\alpha_{t}\mathrm{d}t\biggr{)},$

where, in the right-hand side, it is obviously understood that the argument of
${\mathcal{J}}$ is a constant control. This shows that the minimizers are
deterministic and constant in time. Using the generic notation $\beta$ for
$\int_{0}^{1}\alpha_{t}\mathrm{d}t$, the problem is thus to minimize

${\mathscr{J}}(\beta)=\beta^{2}+G(\beta),\quad\beta\in{\mathbb{R}}.$

Of course, we observe that any minimizer of ${\mathscr{J}}$ is also a solution
of the equation $\beta=-g(\beta)/2$: We recover the fact that any solution to
the mean field control problem is an MFG solution. Obviously, the converse may
not be true. Accordingly, the set of solutions to the mean field control
problem may be strictly included in the set of equilibria of the corresponding
mean field game. This principle is illustrated by Figure 5 below with $g$ as
in (3.2): For all the values of $\kappa\in\\{2,3,\cdots,10\\}$, the potential
${\mathscr{J}}$ has a unique minimizer; but, for $\kappa\in\\{7,8,9,10\\}$,
the derivative of ${\mathscr{J}}$ has several zeros, hence proving that the
MFG has several solutions.

(a) Evolution of ${\mathscr{J}}$ with $\kappa$.

(b) Evolution of $\frac{\mathrm{d}}{\mathrm{d}\beta}{\mathscr{J}}$ with
$\kappa$.

Figure 5. Minimizers of the mean field control problem and equlibria of the
mean field game for $\kappa\in\\{2,\cdots,10\\}$ when $d=1$ and
$g(x)=\cos(\kappa x)$: Solutions of the mean field control problem are the
minimizers of ${\mathcal{J}}$; Solutions of the mean field game are the roots
of $\frac{\mathrm{d}}{\mathrm{d}\beta}{\mathscr{J}}$.

Interestingly, the recent contribution [23] says that, when the MFG is
potential, the equilibria that are not associated with a minimizer of the
potential should be ruled out. Equivalently, only the MFG solutions that
minimize (globally) the potential should be selected. Even more, the guess
(which has been rigorously established in [23, 30] but in different settings
that the one addressed here) is that the minimizers of the potential are
precisely those that appear by forcing uniqueness with a common noise and then
by tuning down the intensity of the common noise. We check this prediction
numerically in Subsection 3.4, by using our learning method.

###  3.2. Algorithms for a fixed intensity of the common noise

We here explain how to implement the fictitious play in the two benchmark
examples (3.2) and (3.3) along the lines of Figure 3. Throughout the
subsection, the intensity of the common noise is fixed. For simplicity, the
intensity of the common noise is chosen as $\varepsilon=1$ and the intensity
of the independent noise is chosen as $\sigma=0$ or $1$. Numerical experiments
are discussed in the next subsection.

####  3.2.1. Numerical reference solution

We first explain how to compute a reference solution by solving numerically
the decoupled forward-backward equation (2.6). Observe that there is no
contradiction in using a numerical method for testing our learning algorithm:
the numerical method relies on the explicit knowledge of the coefficient $g$,
whilst the learning algorithm is intended to work on the sole observations of
the costs.

The numerical method we use is tailor-made to our problem. Indeed, we observe
that the solution to the backward equation in (2.6) writes

$h_{t}={\mathbb{E}}^{{\bm{h}}}\biggl{[}\exp\biggl{(}-\int_{t}^{1}\eta_{s}\mathrm{d}s\biggr{)}g\bigl{(}m_{1}\bigr{)}\,\big{|}\,{\mathcal{F}}^{{\bm{W}}}_{t}\biggr{]},\quad
t\in[0,1],$ (3.9)

where $(m_{t})_{0\leq t\leq 1}$ solves is an Ornstein-Uhlenbeck process (with
independent coordinates)

$\mathrm{d}m_{t}=-\eta_{t}m_{t}\mathrm{d}t+\mathrm{d}W_{t},\quad t\in[0,1]\
;\quad m_{0}=0.$ (3.10)

We then employ a Picard scheme for solving the fixed point equation (3.11), by
iterating

$h_{t}^{\textrm{\rm ref},n+1}={\mathbb{E}}^{{\bm{h}}^{\textrm{\rm
ref},n}}\biggl{[}\exp\biggl{(}-\int_{t}^{1}\eta_{s}\mathrm{d}s\biggr{)}g\bigl{(}m_{1}\bigr{)}\,\big{|}\,{\mathcal{F}}^{{\bm{W}}}_{t}\biggr{]},\quad
t\in[0,1],$ (3.11)

Numerically, we use a time grid $\\{t_{k}=k/p\\}_{k}$ with $p$ uniform steps
(for the same $p$ as the one used in the learning method, which makes the
comparison easier). Following the seminal work of [40], we find, at each
iteration $n$ of the Picard sequence and at each node $t_{j}>0$ of the time
mesh, an approximation of the conditional expectation in the right-hand side
of (3.11) in the form of a (deterministic) function of $m_{t_{j}}$ (the
forward component in (3.11)). This function is taken in the linear span of
$\\{H_{\ell}^{d}(\cdot/\sigma_{t_{j}})\\}_{|\ell|\leq D}$, where
$\sigma_{t_{j}}$ is the common standard deviation of the coordinates of
$m_{t_{j}}$ (which is explicitly computable) and
$\\{H_{\ell}^{d}(\cdot)\\}_{|\ell|\leq D}$ is the collection of Hermite
polynomials of dimension $d$ ($d=1,2$ in our case) and of degree less than
$D$, for a given integer $D$ (here $\ell$ is a $d$-tuple of integers and
$|\ell|$ is the sum of its entries). For completeness, we recall that, when
$d=1$,

$H_{\ell}^{1}(x)=\frac{(-1)^{\ell}}{\sqrt{2^{\ell}\ell!}}e^{x^{2}}\frac{\mathrm{d}}{\mathrm{d}x^{\ell}}[e^{-x^{2}}],\quad
x\in{\mathbb{R}}.$

When $d=2$,

$H_{(\ell_{1},\ell_{2})}^{2}(x_{1},x_{2})=H_{\ell_{1}}^{1}(x_{1})H_{\ell_{2}}^{1}(x_{2}),\quad
x_{1},x_{2}\in{\mathbb{R}}.$

The weights in the Hermite decomposition are computed by a regression method:

Algorithm 1. [Reference solution]

  1. Input:

Introduce $\Delta_{t_{k}}w^{(j)}$ with $j\in\\{1,\cdots,N\\}$ and
$k\in\\{1,\cdots,p\\}$, a collection of $N$ realizations of independent
Gaussian variables ${\mathcal{N}}(0,I_{d})$.

  2. Task:

For each $i$, compute $(m^{(j)}_{t_{k}})_{j=0,\cdots,p}$ the realizations of
the Euler scheme associated with (3.10) and driven by the simulations

$\biggl{(}w^{(j)}_{t_{k}}=\frac{1}{\sqrt{p}}\Bigl{[}\Delta_{t_{1}}w^{(j)}+\cdots+\Delta_{t_{k}}w^{(j)}\Bigr{]}\biggr{)}_{k=1,\cdots,p}.$

  3. Loop:

One iteration consists of one Picard iteration. The input at iteration number
$n$ is encoded in the form of an array $(h_{t_{k}}^{n,(j)})_{j,k}$ with
entries in ${\mathbb{R}}^{d}$. The following procedure returns the next input
for iteration number $n+1$.

    1. (a)

Look for $(h_{t_{k}}^{n+1,(j)})_{j,k}$ in terms of $(h_{t_{k}}^{n,(j)})_{j,k}$
by solving, for each $k=0,\cdots,p-1$, the minimization problem

$\min_{{\bm{c}}=(c_{\ell})_{|\ell|\leq
D}}\frac{1}{N}\sum_{j=1}^{N}{\mathcal{E}}^{n,(j)}_{k}\biggl{|}\exp\biggl{(}-\sum_{s=k}^{p-1}\eta_{t_{s}}\biggr{)}g(m_{1}^{(j)})-\sum_{|\ell|\leq
D}c_{\ell}H^{d}_{\ell}\bigl{(}m_{t_{k}}^{(j)}\bigr{)}\biggr{|}^{2},$

where each $c_{\ell}$ in ${\bm{c}}=(c_{\ell})_{|\ell|\leq D}$ is in
${\mathbb{R}}^{d}$ and with

${\mathcal{E}}^{n,(j)}_{k}=\exp\biggl{(}\sum_{s=k}^{p-1}h_{t_{s}}^{n,(j)}\cdot\Delta_{t_{k+1}}w^{(j)}-\frac{1}{2p}\sum_{s=k}^{p-1}|h_{t_{s}}^{n,(j)}|^{2}\biggr{)}.$

(When $k=0$, only $c_{\bm{0}}$, where $|{\bm{0}}|=0$, matters and it is given
by a mere empirical mean.)

    2. (b)

To enforce some stability, update each coordinate of
$(h_{t_{k}}^{n+1,(j)})_{j,k}$ by projecting it onto $[-1,1]$.

######  Remark 3.2.

Here, the coordinates in (3.10) have the same marginal variances because
$(\eta_{t})_{0\leq t\leq T}$ is scalar-valued (or equivalent is a $d$-diagonal
matrix). This is no longer true when $R$ and $Q$ in (1.2) are general
matrices. In such a case (and more generally when the components are not
centered), the components of $(m_{t})_{0\leq t\leq T}$ are no longer
independent and the ‘normalized’ Hermite polynomials
$\\{H_{\ell}^{d}(\cdot/\sigma_{t_{j}})\\}_{|\ell|\leq D}$ used above should be
instead replaced by
$\\{H_{\ell}^{d}({\mathbb{K}}^{-1/2}(m_{t_{j}})(\cdot-{\mathbb{E}}(m_{t_{j}}))\\}_{|\ell|\leq
D}$, where ${\mathbb{K}}(m_{t_{j}})$ is the covariance matrix of $m_{t_{j}}$
and ${\mathbb{K}}^{1/2}(m_{t_{j}})$ is any root of it (e.g., as given by the
Cholesky decomposition).

####  3.2.2. Numerical approximation by fictitious play

We now address the implementation of the fictitious play according to the
principle stipulated by Figure 3. The black-box therein is discretized in a
suitable manner. It is used to solve numerically the optimization problem
(2.23), and then to implement the updating rule (2.32) for the environment.
Since Lemma 2.11 says that the corresponding optimal law has an affine Markov
feedback form, we thus perform the optimization in (2.23) over closed loop
controls that are affine in the state variable, with the linear coefficient
being a scalar only depending on time and with the intercept possibly
depending on the environment. From the practical point of view, this choice is
fully justified if the model is expected to be linear-quadratic in the
space/action variables (as it is in (1.1)–(1.2), with $Q$ and $R$ being
scalars), but this does not require the coefficients $Q$, $R$, $f$ and $g$ to
be known explicitly. Mathematically, this writes as follows. We restrict
ourselves to controls ${\bm{\alpha}}=(\alpha_{t_{k}})_{{k=1,\cdots,p}}$ of the
form

$\alpha_{t_{k}}=a_{t_{k-1}}X_{t_{k-1}}+C_{t_{k-1}},$

where $a_{t_{k-1}}$ is a scalar and $C_{t_{k-1}}$ is a $d$-dimensional random
vector. In comparison with the formula (2.27), $a_{t_{k-1}}$ should be
understood as a proxy for $-\eta_{t_{k-1}}$ and $C_{t_{k-1}}$ as a proxy for
the intercept therein. Part of the numerical difficulty is thus to have a
tractable regression method for capturing the randomness of $C_{t_{k-1}}$.
Recalling that the intercept ${\bm{h}}$ in the solution of the mean field game
with common noise is a function of the environment $(m_{t})_{0\leq t\leq 1}$,
see Remark 2.4, a key point in our numerical method is to look for
$C_{t_{k-1}}$ as a $\sigma(\overline{m}_{t_{k-1}})$-measurable random
variable, where $\overline{m}_{t_{k-1}}$ is the current proxy for the state of
the environment at time $t_{k-1}$. Numerically, we use a regression method
based on a $d$-dimensional Hermite polynomials of degree less than $D$, for a
fixed value of $D$.

We now make this construction explicit when $\sigma=\varepsilon=1$, noticing
that the algorithm must rely on simulations for the independent and common
noises. Below, we call $M$ the number of particles and $N$ the number of
simulations of the whole system, with $i$ denoting the generic label for a
particle and $j$ denoting the generic label for a realization. We are thus
given a collection

$\Delta_{t_{k}}b^{(i,j)},\ \Delta_{t_{k}}w^{(j)},\quad i\in\\{1,\cdots,M\\},\
j\in\\{1,\cdots,N\\},\ k\in\\{1,\cdots,p\\},$ (3.12)

of realizations of independent ${\mathcal{N}}(0,I_{d})$ Gaussian variables. At
iteration number $n$ of the fictitious play, the current proxy for the
environment is thus given in the form of a collection
$(\overline{m}^{n,(j)}_{t_{k}})_{j,k}$, with $j$ running from $1$ to $N$ and
$k$ running from $0$ to $p$. Similarly, the proxy for the intercept in the
optimal law (2.27) is given in the form of a collection
$(h^{n,(j)}_{t_{k}})_{j,k}$, also with $j$ running from $1$ to $N$ and $k$
running from $0$ to $p-1$. The fact that the two proxies are independent of
$i$ is fully consistent with the fact that $\overline{\bm{m}}^{n}$ and
${\bm{h}}^{n}$ in Remark 2.13 are adapted to the realization of the common
noise. Following our introductory discussion, we solve, as an approximation of
the stochastic control problem (2.23), the minimization problem:

$\min\biggl{\\{}{\mathcal{E}}^{n,(j)}\frac{1}{MN}\sum_{i=1}^{M}\sum_{j=1}^{N}\biggl{[}\frac{1}{2p}\sum_{k=1}^{p}|\alpha_{t_{k}}^{(i,j)}|^{2}+\frac{1}{2}\Bigl{|}x^{(i,j)}_{1}+g\Bigl{(}\overline{m}_{1}^{n,(j)}\Bigr{)}\Bigr{|}^{2}\biggr{]}\biggr{\\}},$
(3.13)

where

${\mathcal{E}}^{n,(j)}=\exp\Biggl{(}-\sqrt{\frac{1}{p}}\sum_{k=0}^{p-1}h^{n,(j)}_{t_{k}}\cdot\Delta_{t_{k+1}}w^{(j)}-\frac{1}{2p}\sum_{k=0}^{p-1}|h^{n,(j)}_{t_{k}}|^{2}\Biggr{)}.$

Moreover, in the minimization problem, $\alpha_{t_{k}}^{(i,j)}$,
$x_{t_{k}}^{(i,j)}$ are required to be of the form

$\begin{split}&x_{t_{k}}^{(i,j)}=x_{t_{k-1}}^{(i,j)}+\frac{1}{p}\alpha_{t_{k}}^{(i,j)}+{\frac{1}{\sqrt{p}}}\Delta_{t_{k}}b^{(i,j)},\quad\ell=1,\cdots,p\
;\quad x_{0}^{(i,j)}=x_{0},\\\
&\alpha_{t_{k}}^{(i,j)}=a_{t_{k-1}}x_{t_{k-1}}^{(i,j)}+C_{t_{k-1}}^{(j)}+h^{n,(j)}_{t_{k-1}}+\sqrt{p}\Delta_{t_{k}}w^{(j)},\quad
k=1,\cdots,p.\end{split}$ (3.14)

In particular, $(x_{t_{k}}^{(i,j)})_{k=0,\cdots,p}$ solves the following Euler
scheme:

$x_{t_{k}}^{(i,j)}=x_{t_{k-1}}^{(i,j)}+\frac{1}{p}\Bigl{(}a_{t_{k-1}}x_{t_{k-1}}^{(i,j)}+C_{t_{k-1}}^{(j)}+h^{n,(j)}_{t_{k-1}}\Bigr{)}+{\frac{1}{\sqrt{p}}}\Delta_{t_{k}}b^{(i,j)}+{\frac{1}{\sqrt{p}}}\Delta_{t_{k}}w^{(j)},\quad
k=1,\cdots,p.$

Moreover, $C^{(j)}_{t_{k}}$ is required to be of the form

$C^{(j)}_{t_{k}}=\sum_{|\ell|\leq
D}c_{t_{k}}(\ell)H^{d}_{\ell}\biggl{(}\bigl{(}U_{t_{k}}^{n}\bigr{)}^{-1}\Bigl{(}\overline{m}^{n,(j)}_{t_{k}}-\frac{1}{N}\sum_{r=1}^{N}\overline{m}^{n,(r)}_{k}\Bigr{)}\biggr{)},$
(3.15)

where $c_{k}(\ell)\in{\mathbb{R}}^{d}$ and $U_{t_{k}}^{n}$ is the diagonal
matrix obtained by taking the roots of the diagonal of the empirical
covariance matrix:

$\Sigma_{t_{k}}^{n}=\frac{1}{N}\sum_{j=1}^{N}\Bigl{(}\overline{m}^{n,(j)}_{t_{k}}-\frac{1}{N}\sum_{r=1}^{N}\overline{m}^{n,(r)}_{t_{k}}\Bigr{)}\otimes\Bigl{(}\overline{m}^{n,(j)}_{t_{k}}-\frac{1}{N}\sum_{r=1}^{N}\overline{m}^{n,(r)}_{t_{k}}\Bigr{)}.$

######  Remark 3.3.

Following Remark 3.2, we can define $U_{t_{k}}^{n}$ (in a more general
fashion) as the upper triangular matrix given by the Cholesky decomposition of
$\Sigma_{t_{k}}^{n}$ in the form
$\Sigma_{t_{k}}^{n}=(U_{t_{k}}^{n})^{\dagger}U_{t_{k}}^{n}$.

The minimization problem is thus defined over
${\bm{a}}=(a_{t_{k}})_{k}\in{\mathbb{R}}^{p}$ and
${\bm{c}}=(c_{k}(\ell))_{k,\ell}\in({\mathbb{R}}^{d})^{p\times L}$, where $L$
is the number of distinct $d$-tuples of (non-negative) integers whose sum is
less than $D$.

We can summarize the algorithm in the following form.

Algorithm 2. [Learning with two noises, $\sigma=\varepsilon=1$]

  1. Input:

Take as input the realizations (3.12), which are the same at any iteration of
the algorithm.

  2. Loop:

At rank $n+1$ of the fictitious play:

    1. (a)

Take as input the proxies $(\overline{m}^{n,j}_{t_{k}})_{j,k}$ and
$(h^{n,j}_{t_{k}})_{j,k}$ for (respectively) the environment and the
intercept.

    2. (b)

Solve the minimization problem (3.13) over ${\bm{a}}=(a_{t_{k}})_{k}$ in
(3.14) and ${\bm{c}}=(c_{t_{k}}(\ell))_{k,\ell}$ in (3.15). Call
${\bm{a}}^{n+1}=(a^{n+1}_{t_{k}})_{k}$ and
${\bm{c}}^{n+1}=(c^{n+1}_{t_{k}}(\ell))_{k,\ell}$ the optimal points (returned
by any optimization method).

    3. (c)

With ${\bm{a}}^{n+1}=(a_{t_{k}}^{n+1})_{k}$ and
${\bm{c}}^{n+1}=(c^{n+1}_{t_{k}}(\ell))_{k,\ell}$, associate
$(x^{n+1,(i,j)}_{t_{k}})_{i,j,k}$ as in (3.14) and
$(C_{t_{k}}^{n+1,(j)})_{k,j}$ as in (3.15).

    4. (d)

Update the proxies by letting

$\begin{split}&m^{n+1,(j)}_{t_{k}}=\frac{1}{M}\sum_{i=1}^{M}x^{n+1,(i,j)}_{t_{k}},\quad\overline{m}^{n+1,(j)}_{t_{k}}=\frac{1}{n+1}m^{n+1,(j)}_{t_{k}}+\frac{n}{n+1}\overline{m}^{n,(j)}_{t_{k}},\\\
&h^{n+1,(j)}_{t_{k}}=-C_{t_{k}}^{n+1,(j)}.\end{split}$

In some of the numerical examples below, we will use variants of Algorithm 2.
We first focus on the shape of the algorithm when $\sigma=0$, recalling that
our results do not require the presence of an idiosyncratic noise. Implicitly,
we then have a single particle only, i.e. $M=1$. Moreover,
$\Delta_{t_{k}}b^{(i,j)}$ no longer appears in (3.14). We then have the
following variant of Algorithm 2, which is of lower complexity:

Algorithm 3. [Learning with common noise only, $\sigma=0$, $\varepsilon=1$]

  1. Input:

Take as input the realizations $(\Delta_{t_{k}}w^{(j)})_{j,k}$ in (3.12),
which are the same at any iteration of the algorithm.

  2. Loop:

At rank $n+1$ of the fictitious play, apply the same loop as in Algorithm 2,
but with $M=1$ and with $\Delta_{t_{k}}b^{(1,j)}=0$ in (3.14).

We end up our presentation with the case when there is no common noise, i.e.
$\sigma=1$ and $\varepsilon=0$. In that case, our strategy no longer applies.
The point is thus to implement the standard version of the fictitious play.
Accordingly, there is no Girsanov density in (3.13) and
$\alpha_{t_{k}}^{(i,j)}$ in (3.14) merely writes

$\alpha_{t_{k}}^{(i,j)}=a_{t_{k-1}}x_{t_{k-1}}^{(i,j)}+C_{t_{k-1}},\quad
k=1,\cdots,p.$

with $C_{t_{k}}$ being a deterministic $d$-dimensional vector (which is
independent of $j$). In particular, there is no need to use the proxy
$(h^{n,j}_{t_{k}})_{j,k}$ for the intercept. Equivalently, we can assume that
$N=1$ and $D=0$. The rest of the algorithm is similar. It may be written as
follows.

Algorithm 4. [Learning with idiosyncratic noise only, $\sigma=1$,
$\varepsilon=0$]

  1. Input:

Take as input the realizations $(\Delta_{t_{k}}b^{(i,1)})_{k}$ in (3.12),
which are the same at any iteration of the algorithm.

  2. Loop:

At rank $n+1$ of the fictitious play:

    1. (a)

Take as input the proxy $(\overline{m}^{n}_{t_{k}})_{k}$ for the environment
and the intercept.

    2. (b)

Solve the minimization problem (3.13) over ${\bm{a}}=(a_{t_{k}})_{k}$ in
(3.14) and ${\bm{c}}=(c_{t_{k}})$ in (3.15), with ${\mathcal{E}}^{n,(j)}=1$ in
(3.13) and with the last line in (3.14) being replaced by

$\alpha_{t_{k}}^{(i,1)}=a_{t_{k-1}}x_{t_{k-1}}^{(i,1)}+c_{t_{k-1}}(0),\quad
k=1,\cdots,p.$

Call ${\bm{a}}^{n+1}=(a^{n+1}_{t_{k}})_{k}$ and
${\bm{c}}^{n+1}(0)=(c^{n+1}_{t_{k}}(0))_{k}$ the optimal points.

    3. (c)

With ${\bm{a}}^{n+1}=(a_{t_{k}}^{n+1})_{k}$ and
${\bm{c}}^{n+1}=(c^{n+1}_{t_{k}}(0))_{k}$, associate
$(x^{n+1,(i,1)}_{t_{k}})_{i,k}$ as in (3.14) (with the same prescription as
above).

    4. (d)

Update the proxy by letting

$\begin{split}&m^{n+1,(1)}_{t_{k}}=\frac{1}{M}\sum_{i=1}^{M}x^{n+1,(i,1)}_{t_{k}},\quad\overline{m}^{n+1,(1)}_{t_{k}}=\frac{1}{n+1}m^{n+1,(1)}_{t_{k}}+\frac{n}{n+1}\overline{m}_{t_{k}}^{n,(1)}.\end{split}$

In all the numerical experiments, we use ADAM optimizer (as implemented in
TensorFlow) in order to solve the optimization step in Algorithms 2, 3 and 4.
We recall from the discussion in Subsection 1.7 that the code in TensorFlow
relies on some automatic differentiation and thus makes an explicit use of the
linear-quadratic form of the coefficients. In this sense, our numerical
implementation requires part of the model to be known, but, as we already
explained in Subsection 1.7, this does not really affect the scope of our
conclusions: Provided the optimization method used in Algorithms 2, 3 and 4 is
sufficiently accurate, the whole works well and demonstrates the interest for
exploring the state space by means of the common noise.

###  3.3. Numerical experiments

The numerical examples that we provide below are here to illustrate several
features of our form of fictitious play. In all the examples treated below, we
choose $d=2$, $f\equiv 0$ and $g$ as in (3.3) with $\kappa=10$.

We emphasize that our aim is to demonstrate, numerically, the positive impact
of the common noise, not to compare the numerical rates of convergence
obtained in the numerical experiments with the bounds obtained in the
theoretical analysis. The reason is that the algorithms we present below
include additional features (like numerical optimization methods and
regressions) that are not addressed in the theoretical analysis. This makes
difficult any precise comparison. Also, it is fair to say that the time of
execution of the full algorithm becomes quite long (several hours for the
longest experiments), even for low values of $n$ and $p$ ($n$ less than 50 and
$p$ less than 100).

####  3.3.1. Evolution of the learnt cost with one type of noise only

We first test the influence of each of the two types of noises onto the
behavior of the algorithm. We thus compare the evolution of the cost learnt by
the fictitious play in two scenarios: In the first scenario, there is an
independent noise, but no common noise (Algorithm 4); in the second scenario,
there is a common noise but no independent noise (Algorithm 3). We will see
next what happens when there are the two noises.

(a) With independent but no common noise

(b) With common but no independent noise

Figure 6. Comparison of the learnt cost depending on the type of noise. In
$x$-axis: number of iterations; In $y$-axis; learnt cost.

The experiments are computed with: $n=20$ learning iterations, $p=30$ time
steps, $M=4\times 10^{5}$, $N=1$, $\sigma=1$ and $\varepsilon=0$ in case (A)
and $n=p=20$, $M=1$, $N=4\times 10^{5}$, $\sigma=0$, $\varepsilon=1$ and $D=4$
in case (B). In ADAM method, the learning rate is $0.01$, with 15 epochs.

In plot (B), the orange line is the theoretical equilibrium cost as computed
by the BSDE method explained in Subsection 3.2.1. In plot (A), there is no
computed reference cost. As we already explained, there might not be a unique
equilibrium and the notion of reference cost no longer makes sense. Notice by
the way that the equilibrium cost in case (B) is not an equilibrium cost in
case (A) since the two problems are different (as being set over different
forms of dynamics). Anyway, the conclusion is clear: the learnt cost does not
convergence in case (A), whilst it does in case (B). This is a clear evidence
of the numerical impact of the common noise onto the behavior of the
fictitious play.

####  3.3.2. Error in the optimal path/intercept with common noise and no
independent noise

We now compute the error achieved by the learning algorithm. Using the same
notation as in the statement of Proposition 2.19, we focus on the following
$L^{2}$ error:

$\biggl{[}{\mathbb{E}}\int_{0}^{1}\Bigl{(}|\overline{m}_{t}^{n}-m_{t}|^{2}+|h_{t}^{n}-h_{t}|^{2}\Bigr{)}\mathrm{d}t\biggr{]}^{1/2},$

which slightly differs from the error appearing in Proposition 2.19. Focusing
here on the time-averaged error is obviously more advantageous from the
numerical point of view, but it is sufficient to demonstrate the efficiency of
the algorithm and also to identify some of its limitations. Notice also that
the expectation in the above error is taken under the non-tilted expectation.
Since $\varepsilon$ is here equal to 1, the Girsanov density
${\mathcal{E}}({\bm{h}})$ and its inverse can indeed be easily estimated.
Also, since we have bounds for the moments of $\overline{\bm{m}}^{n}$ and
${\bm{m}}$, Proposition 2.19 shows that the above expectation tends to 0 as
$n$ tends to $\infty$. This is sufficient for our purpose since we do not try
to identify the numerical rate of convergence precisely.

Numerically, the error is approximated by

$\biggl{[}\frac{1}{Np}\sum_{j=1}^{N}\sum_{k=0}^{p-1}|\overline{m}_{t_{k}}^{n,(j)}-m_{t_{k}}^{(j)}|^{2}+|h_{t_{k}}^{n,(j)}-h_{t_{k}}^{(j)}|^{2}\biggr{]}^{1/2},$
(3.16)

where $(\overline{m}_{t_{k}}^{n,(j)})_{j,k}$ and $(h_{t_{k}}^{n,(j)})_{j,k}$
are the returns of Algorithm 2 or 3 (depending on the value of $\sigma$) and
$(\overline{m}_{t_{k}}^{(j)})_{j,k}$ and $(h_{t_{k}}^{(j)})_{j,k}$ are the
reference solutions as computed with Algorithm 1 (with a sufficiently high
number of Picard iterations: 10 in practice).

The first experiment is performed by running Algorithm 3 (no idiosyncratic
noise), but assuming that the solution $(\eta_{t})_{0\leq t\leq 1}$ to the
Riccati equation (3.4) is known, see (3.5): In Algorithm 3, $a_{t_{k}}$ is
replaced by $\eta_{t_{k}}$. This amounts to say that the coefficients $Q$ and
$R$ in (1.2) are explicitly known. The evolution of the error with the number
of iterations is plotted in (A), Figure 7. The second experiment is performed
by running Algorithm 3, but without any further knowledge of the solution of
the Riccati equation (3.4). The evolution of the error with the number of
iterations is plotted in (B), Figure 7.

(a) Riccati known

(b) Riccati unknown

Figure 7. Comparison of the error returned by Algorithm 3 (common noise only),
depending on whether the solution to the Riccati equation is known or not.

The experiments are computed with: $n=20$ learning iterations, $p=30$ time
steps, $M=1$, $N=5\times 10^{4}$, $\sigma=0$, $\varepsilon=1$ and $D=4$. On
Plot (A), the error is less than $0.01$ after iteration 3.

The conclusion is quite clear: When randomness only comes from the common
noise, it is very difficult for the optimization method to distinguish between
$a_{t_{k-1}}x_{t_{k-1}}^{(1,j)}$ and $C_{t_{k-1}}^{(j)}$ in (3.14).

####  3.3.3. Error in the optimal path/intercept with both noises

We now proceed with the same analysis but putting the two noises. As the total
number of simulated paths is $N\times M$, this may be however rather costly.
We proceed below by freezing the quantity $N\times M$.

(a) $M=4\times 10^{4}$, $N=1$, time (cpu): 4867

(b) $M=2\times 10^{4}$, $N=20$, time (cpu): 4096

(c) $M=10^{4}$, $N=40$, time (cpu): 8409

Figure 8. Comparison of the error returned by Algorithm 2, depending on the
number of particles/realizations.

The experiments are computed with: $n=10$ learning iterations, $p=30$ time
steps, $\sigma=1$ and $\varepsilon=1$; 4000 units in cpu time is around 1h.
The conclusion is that idiosyncratic noises demonstrate to be useful from the
numerical point of view, even though the results stated in Section 2 remain
the same whatever the value of $\sigma$. A careful inspection of the numerical
results show that, in fact, idiosyncratic noises provide a better fit of the
solution to the Riccati equation, which is fully consistent with the
conclusion of the previous paragraph. In turn, we guess that a model-based
method, in which $Q$ and $R$ first, and then the solution to the Riccati
equation, would be learnt separately, would make sense in this specific
setting.

####  3.3.4. Validation

The reader could worry about a possible overfitting in our numerical
experiments. Actually, in order to validate our results, we can learn the
coefficients ${\bm{a}}^{n}$ and ${\bm{c}}^{n}$ in Algorithm 2 for a first
series of data in (3.12). In brief, this first series of data is used to learn
the equilibrium feedback. Then, we can use a second series of data (say, of
the same size) in (3.12) in order to compute the error. In clear, given this
second series of data, we can implement Algorithm 1 and then compute the
resulting error (3.16) when $(\overline{m}_{t_{k}}^{n,(j)})_{j,k}$ and
$(h_{t_{k}}^{n,(j)})_{j,k}$ therein are obtained from the formulas (3.14) and
(3.15) by using the coefficients ${\bm{a}}^{n}$ and ${\bm{c}}^{n}$ returned by
the first series of data and then by implementing the Euler scheme with
respect to the second series of data. The resulting plots are very similar to
those represented in Figure 8. For this reason, we feel useless to insert them
here. Intuitively, what happens is that there are sufficiently many
realizations of the common noise in our experiments to guarantee, in the
regression (3.13), a convenient form of averaging with respect to all these
realizations.

###  3.4. Vanishing viscosity

The next question in our numerical experiments is to address the behavior of
the algorithm as the viscosity tends to $0$. As made clear by the analysis
performed in Section 2, the influence of the small viscosity may manifest in
several manners, depending on the error that is addressed: recall indeed that
the bound in the statement of Proposition 2.19 does not depend on
$\varepsilon$, whilst the bound provided by the main Theorem 2.14 blows up as
$\varepsilon$ tends to $0$.

For sure, this is a main issue from the numerical point of view since we may
expect the numerical method to become somewhat instable as $\varepsilon$ tends
to $0$. Independently of the bounds provided by the results, the numerical
difficulty is easy to understand: The Girsanov density underlying our
algorithm becomes highly singular and this may produce high variances in the
computations. Our strategy to overcome this obvious drawback is inspired from
simulated annealing method: The viscosity parameter is decreased step by step;
At each new step, we use the return of the algorithm at the previous value of
the viscosity as a new guess in our algorithm. As exemplified below (see
Figure 10), this turns out to be an efficient method in the benchmark example
introduced in Subsection 3.1.

The use of the previous return (at the previous value of the viscosity) may be
formalized as follows. Assume indeed that we have a descreasing sequence
$(\varepsilon_{q})_{q\geq 1}$ that converges to $0$, with each
$\varepsilon_{q}>0$ being understood as the value of the viscosity at step
$q$. With each $q$, we may associate a return ${\bm{h}}^{\infty,q}$. The
symbol $\infty$ in superscript should be understood as an ‘idealization’ in
the interpretation of ${\bm{h}}^{\infty,q}$: This is the return of our
algorithm, when run with viscosity $\varepsilon_{q}$, after sufficiently many
iterations $n$ of the fictitious play; Since we do not quantify here which
value of $n$ we use (this would be doable, but this would lead to additional
cumbersome notations and computations), we just think of $n$ as being large,
ideally equal to $\infty$, whence our notation. At the next step $q+1$, when
using viscosity $\varepsilon_{q+1}$, we then implement (1.19) (for simplicity,
we also take the time discretization parameter $1/p$ to be equal to $0$, but
the same idea can be obviously implemented when $p$ is finite as in (2.23))
with ${\bm{h}}^{n}-{\bm{h}}^{\infty,q}$ in lieu of ${\bm{h}}^{n}$. In words,
(1.19) becomes

${\bm{\alpha}}^{n+1}=\textrm{\rm
argmin}_{\bm{\alpha}}{\mathbb{E}}^{{\bm{h}}^{n}-{\bm{h}}^{\infty,q}}\bigl{[}{\mathcal{R}}\bigl{(}{\bm{\alpha}};\overline{\bm{m}}^{n};{\bm{W}}^{{\bm{h}}^{n}-{\bm{h}}^{\infty,q}}\bigr{)}\bigr{]}.$
(3.17)

Implicitly, the viscosity that is used in the above formula is
$\varepsilon_{q+1}$. Also, it may be useful to recall that $n$ here denotes
the iteration in the fictitious play: Whatever the value of $n$, we use the
same guess ${\bm{h}}^{\infty,q}$.

The reader must understand that this modification of the guess does not change
our theoretical analysis. In short, the analysis performed in Section 2 is
written with $0$ as guess. The fact that our analysis still holds can be
justified by changing the probability space. We may indeed equip
$(\Omega,{\mathcal{A}})$ with the new probability measure
$\check{\mathbb{P}}={\mathbb{P}}^{-{\bm{h}}^{\infty,q}}$ and with the new
common noise $\check{\bm{W}}={\bm{W}}^{-{\bm{h}}^{\infty,q}}$. Under
$\check{\mathbb{P}}$, the common noise $\check{\bm{W}}$ is a Brownian motion.
Since the latter is the same in (3.17) for any $n$ (which is a prescription in
our method), we can invoke the results of Section 2 in order to study (3.17).
For instance, Proposition 2.8 gives

$\check{\mathbb{E}}^{\bm{h}}\Bigl{[}\sup_{0\leq t\leq
T}\Bigl{(}|\overline{m}_{t}^{n}-m_{t}|^{2}+|h_{t}^{n}-h_{t}|^{2}\Bigr{)}\Bigr{]}\leq\frac{C}{n^{2}},$

where ${\bm{h}}$ is a sloppy notation for the backward component of the
solution to the MFG with common noise of intensity $\varepsilon_{q+1}$. If, as
we said above, ${\bm{h}}^{\infty,q+1}$ is indeed understood as the limit of
${\bm{h}}^{n}$ as $n$ tends to $\infty$, then we can identify ${\bm{h}}$ with
${\bm{h}}^{\infty,q+1}$ itself, hence leading to

$\check{\mathbb{E}}^{{\bm{h}}^{\infty,q+1}}\Bigl{[}\sup_{0\leq t\leq
T}\Bigl{(}|\overline{m}_{t}^{n}-m_{t}^{\infty,q+1}|^{2}+|h_{t}^{n}-h_{t}^{\infty,q+1}|^{2}\Bigr{)}\Bigr{]}\leq\frac{C}{n^{2}},$

where, for consistency, we wrote ${\bm{m}}^{\infty,q+1}$ for the forward
component of the solution to the MFG with common noise of intensity
$\varepsilon_{q+1}$. Recalling the definition of $\check{\bm{P}}$, the above
may be rewritten as

${\mathbb{E}}^{{\bm{h}}^{\infty,q+1}-{\bm{h}}^{\infty,q}}\Bigl{[}\sup_{0\leq
t\leq
T}\Bigl{(}|\overline{m}_{t}^{n}-m_{t}^{\infty,q+1}|^{2}+|h_{t}^{n}-h_{t}^{\infty,q+1}|^{2}\Bigr{)}\Bigr{]}\leq\frac{C}{n^{2}}.$

Now, if we tune the increment $\varepsilon_{q}-\varepsilon_{q+1}$ in such a
manner that ${\bm{h}}^{\infty,q+1}$ and ${\bm{h}}^{\infty,q}$ are sufficiently
close compared to $\varepsilon_{q}$, it becomes reasonable to expect that the
same expectation but under the historical probability measure (hence without
the Girsanov transformation) is small, namely

$\lim_{n\rightarrow 0}{\mathbb{E}}\Bigl{[}\sup_{0\leq t\leq
T}\Bigl{(}|\overline{m}_{t}^{n}-m_{t}^{\infty,q+1}|^{2}\Bigr{)}\Bigr{]}=0,$

with a precise decay that would depend on the bound for the Girsanov density.
We feel better not to address this question mathematically here. This will be
the subject of a separate note. However, we prove the guess to be satisfied in
a numerical experiment based on the aforementioned benchmark example. We
choose a variant of $g$ in (3.2), with $d=1$:

$g(x)=\cos\bigl{(}\kappa(x-x_{0})\bigr{)}-2x_{0},$ (3.18)

with $\kappa=10$ and where $x_{0}$ is a root of the equation

$\cos\bigl{(}\kappa x_{0}\bigr{)}=2x_{0}.$

Numerically, we find that a choice is $x_{0}\approx-0.384$. The motivation for
such an $x_{0}$ is that $0$ is a solution of (3.8). In other words, $0$ is an
equilibrium. Even more, we observe that, if the viscosity is zero, then the
iterative sequence defined through the two updating rules (1.9) and (1.10)
remains in $({\bm{m}}^{n},{\bm{h}}^{n})=(0,-2x_{0})$ for any $n\geq 1$ if
${\bm{m}}^{0}$ is chosen as $0$. In other words, the standard fictitious play
converges (as predicted by the theory since this 1-dimensional MFG is
potential, see §3.1.2) and chooses the $0$-equilibrium. Now, in order to
implement the same prediction method as in §3.1.2, we plot the corresponding
potential (which is given by a primitive $G$ of $g$). The plots are given in
Figure 9: We observe that $0$ is just a local minimizer of the potential and
that the global minimizer is around $-0.5$: Our selection method is hence to
rule the $0$-equilibrium and to retain the other one. In order to check this
numerically, we apply the vanishing viscosity method (3.17), with the
following values: We choose $\sigma=0$ (no idiosyncratic noise), $p=10$ time
steps, $M=10^{6}$ Monte-Carlo simulations, $n=10$ iterations in the fictitious
play. Our regression on the Hermite polynomials goes up to polynomials of
degree 6; We start from $\varepsilon=1$ as viscosity and we decrease by $0.1$
its value up to viscosity $0.1$; at new each step, we use the last return of
the fictitious play as guess in (3.17). The plot below (Figure 10) shows that
our method is able to capture the global minimizer even if initialized from
the local minimizer.

Figure 9. Equilibrium predicted by the potential rule: Potential $G$ on the
left pane; Zeros of the function $g$ in the right pane.

Figure 10. Selection of a solution by vanishing viscosity: Plots of the
histogram of the approximate mean state of the population at terminal time.

## Acknowledgement

The authors are grateful to Mathieu Laurière for very helpful discussions on
the topic.

## References

  * [1] Y. Achdou, F. Camilli, and I. Capuzzo-Dolcetta.  Mean field games: numerical methods for the planning problem.  SIAM J. Control Optim., 50(1):77–109, 2012.
  * [2] Y. Achdou, F. Camilli, and I. Capuzzo-Dolcetta.  Mean field games: convergence of a finite difference method.  SIAM J. Numer. Anal., 51(5):2585–2612, 2013.
  * [3] Y. Achdou and I. Capuzzo-Dolcetta.  Mean field games: numerical methods.  SIAM J. Numer. Anal., 48(3):1136–1162, 2010.
  * [4] Y. Achdou and M. Laurière.  Mean field type control with congestion (II): An augmented Lagrangian method.  Appl. Math. Optim., 74(3):535–578, 2016.
  * [5] Y. Achdou and M. Laurière.  Mean field games and applications: Numerical aspects.  In Mean Field Games, Cetraro, Italy 2019, Cardaliaguet, Pierre, Porretta, Alessio (Eds.), LNM 2281, pages 203–248. Springer, 2021.
  * [6] C. Alasseur, I. Ben Taher, and A. Matoussi.  An extended mean field game for storage in smart grids.  J. Optim. Theory Appl., 184(2):644–670, 2020.
  * [7] E. Bayraktar, A. Cecchin, A. Cohen, and F. Delarue.  Finite state mean field games with Wright-Fisher common noise.  J. Math. Pures Appl., 147:98–162, 2021.
  * [8] C. Beck, W. E, and A. Jentzen.  Machine learning approximation algorithms for high-dimensional fully nonlinear partial differential equations and second-order backward stochastic differential equations.  J. Nonlinear Sci., 29(4):1563–1619, 2019.
  * [9] J.-D. Benamou and G. Carlier.  Augmented Lagrangian methods for transport optimization, mean field games and degenerate elliptic equations.  J. Optim. Theory Appl., 167(1):1–26, 2015.
  * [10] C. Bender and J. Zhang.  Time discretization and Markovian iteration for coupled FBSDEs.  Ann. Appl. Probab., 18(1):143–177, 2008.
  * [11] L. Briceño Arias, D. Kalise, Z. Kobeissi, M. Laurière, A. Mateos González, and F. J. Silva.  On the implementation of a primal-dual algorithm for second order time-dependent mean field games with local couplings.  In CEMRACS 2017—numerical methods for stochastic models: control, uncertainty quantification, mean-field, volume 65 of ESAIM Proc. Surveys, pages 330–348. EDP Sci., Les Ulis, 2019.
  * [12] L. M. Briceño Arias, D. Kalise, and F. J. Silva.  Proximal methods for stationary mean field games with local couplings.  SIAM J. Control Optim., 56(2):801–836, 2018.
  * [13] P. Cardaliaguet and S. Hadikhanloo.  Learning in mean field games: the fictitious play.  ESAIM Control Optim. Calc. Var., 23(2):569–591, 2017.
  * [14] E. Carlini and F. J. Silva.  A fully discrete semi-Lagrangian scheme for a first order mean field game problem.  SIAM J. Numer. Anal., 52(1):45–67, 2014.
  * [15] E. Carlini and F. J. Silva.  A semi-Lagrangian scheme for a degenerate second order mean field game system.  Discrete Contin. Dyn. Syst., 35(9):4269–4292, 2015.
  * [16] R. Carmona and F. Delarue.  Probabilistic theory of mean field games with applications. I, volume 83 of Probability Theory and Stochastic Modelling.  Springer, Cham, 2018.  Mean field FBSDEs, control, and games.
  * [17] R. Carmona and F. Delarue.  Probabilistic theory of mean field games with applications. II, volume 84 of Probability Theory and Stochastic Modelling.  Springer, Cham, 2018.  Mean field games with common noise and master equations.
  * [18] R. Carmona and D. Lacker.  A probabilistic weak formulation of mean field games and applications.  Ann. Appl. Probab., 25(3):1189–1231, 2015.
  * [19] R. Carmona and M. Laurière.  Convergence analysis of machine learning algorithms for the numerical solution of mean-field control and games: I-the ergodic case.  arXiv e-prints, arXiv:1907.05980, 2019.
  * [20] R. Carmona and M. Laurière.  Convergence analysis of machine learning algorithms for the numerical solution of mean-field control and games: II-the finite horizon case.  arXiv e-prints, arXiv:1908.01613, 2019.
  * [21] R. Carmona, M. Laurière, and Z. Tan.  Linear-quadratic mean-field reinforcement learning: convergence of policy gradient methods.  arXiv e-prints, arXiv:1910.04295, 2019.
  * [22] R. Carmona, M. Laurière, and Z. Tan.  Model-free mean field reinforcement learning: mean field MDP and mean-field Q-learning.  arXiv e-prints, arXiv:1910.12802, 2019.
  * [23] A. Cecchin and F. Delarue.  Selection by vanishing common noise for potential finite state mean field games.  arXiv e-prints, 2020.
  * [24] J.-F. Chassagneux, D. Crisan, and F. Delarue.  Numerical method for FBSDEs of McKean-Vlasov type.  Ann. Appl. Probab., 29(3):1640–1684, 2019.
  * [25] P. E. Chaudru de Raynal and C. A. Garcia Trillos.  A cubature based algorithm to solve decoupled McKean-Vlasov forward-backward stochastic differential equations.  Stochastic Process. Appl., 125(6):2206–2255, 2015.
  * [26] D. Crisan and E. McMurray.  Cubature on Wiener space for McKean-Vlasov SDEs with smooth scalar interaction.  Ann. Appl. Probab., 29(1):130–177, 2019.
  * [27] J. Cvitanić and J. Zhang.  The steepest descent method for forward-backward SDEs.  Electron. J. Probab., 10:1468–1495, 2005.
  * [28] F. Delarue.  On the existence and uniqueness of solutions to FBSDEs in a non-degenerate case.  Stochastic Processes and their Applications, 99(2):209 – 286, 2002\.
  * [29] F. Delarue.  Restoring uniqueness to mean-field games by randomizing the equilibria.  Stochastics and Partial Differential Equations: Analysis and Computations, 7:598–678, 2019.
  * [30] F. Delarue and R. Foguen Tchuendom.  Selection of equilibria in a linear quadratic mean field game.  Stochastic Processes and their Applications, 130(2):1000–1040, 2020\.
  * [31] F. Delarue and S. Menozzi.  A forward-backward stochastic algorithm for quasi-linear PDEs.  Ann. Appl. Probab., 16(1):140–184, 2006.
  * [32] P. Dorato and A. Levis.  Optimal linear regulators: The discrete-time case.  IEEE Transactions on Automatic Control, 16:613–620, 1971.
  * [33] K. Doya.  Reinforcement learning in continuous time and space.  Neural Comput., 12(1):219–245, 2000.
  * [34] W. E, J. Han, and A. Jentzen.  Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations.  Commun. Math. Stat., 5(4):349–380, 2017.
  * [35] R. Elie, L. Pérolat, M. Laurière, M. Geist, and O. Pietquin.  On the convergence of model free learning in mean field games.  AAAI, 2020.
  * [36] D. Firoozi and S. Jaimungal.  Exploratory LQG Mean Field Games with Entropy Regularization.  arXiv:2011.12946, 2020.
  * [37] M. Fischer and G. Nappo.  On the moments of the modulus of continuity of Itô processes.  Stoch. Anal. Appl., 28(1):103–122, 2010.
  * [38] R. Foguen Tchuendom.  Uniqueness for linear-quadratic mean field games with common noise.  Dyn. Games Appl., 8(1):199–210, 2018.
  * [39] M. Geist, J. Pérolat, M. Laurière, R. Élie, S. Perrin, O. Bachem, R. Munos, and O. Pietquin.  Concave Utility Reinforcement Learning: the Mean-field Game viewpoint.  arXiv:2106.03787, 2021.
  * [40] E. Gobet, J.-P. Lemor, and X. Warin.  A regression-based Monte Carlo method to solve backward stochastic differential equations.  Ann. Appl. Probab., 15(3):2172–2202, 2005.
  * [41] X. Guo, A. Hu, R. Xu, and J. Zhang.  Learning mean-field games.  In Proceedings of NeurIPS. 2019.
  * [42] S. Hadikhanloo.  Learning in anonymous nonatomic games with applications to first-order mean field games.  arXiv e-prints, arXiv:1704.00378, 2017.
  * [43] S. Hadikhanloo and F. J. Silva.  Finite mean field games: fictitious play and convergence to a first order continuous mean field game.  J. Math. Pures Appl. (9), 132:369–397, 2019.
  * [44] K. Huang, X. Di, Q. Du, and X. Chen.  A game-theoretic framework for autonomous vehicles velocity control: bridging microscopic differential games and macroscopic mean field games.  Discrete Contin. Dyn. Syst. Ser. B, 25(12):4869–4903, 2020.
  * [45] M. Huang, P. Caines, and R. Malhamé.  Individual and mass behavior in large population stochastic wireless power control problems: centralized and Nash equilibrium solutions.  pages 98 – 103, 2003.
  * [46] M. Huang, P. E. Caines, and R. P. Malhamé.  The Nash certainty equivalence principle and McKean-Vlasov systems: An invariance principle and entry adaptation.  In Decision and Control, 2007 46th IEEE Conference on, pages 121–126. IEEE, 2007.
  * [47] 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.  Commun. Inf. Syst., 6(3):221–251, 2006.
  * [48] T. P. Huijskens, M. J. Ruijter, and C. W. Oosterlee.  Efficient numerical Fourier methods for coupled forward-backward SDEs.  J. Comput. Appl. Math., 296:593–612, 2016.
  * [49] J.-M. Lasry and P.-L. Lions.  Jeux à champ moyen. I. Le cas stationnaire.  C. R. Math. Acad. Sci. Paris, 343(9):619–625, 2006.
  * [50] J.-M. Lasry and P.-L. Lions.  Jeux à champ moyen. II. Horizon fini et contrôle optimal.  C. R. Math. Acad. Sci. Paris, 343(10):679–684, 2006.
  * [51] J.-M. Lasry and P.-L. Lions.  Mean field games.  Jpn. J. Math., 2(1):229–260, 2007.
  * [52] P.-L. Lions.  Cours au Collège de France, Equations aux dérivées partielles et applications.  https://www.college-de-france.fr/site/pierre-louis-lions/course.htm.
  * [53] J. Ma, P. Protter, and J. Yong.  Solving forward-backward stochastic differential equations explicitly – a four step scheme.  Probability Theory and Related Fields, 98(3):339–359, 1994.
  * [54] J. Ma, J. Shen, and Y. Zhao.  On numerical approximations of forward-backward stochastic differential equations.  SIAM J. Numer. Anal., 46(5):2636–2661, 2008.
  * [55] G. N. Milstein and M. V. Tretyakov.  Numerical algorithms for semilinear parabolic equations with small parameter based on approximation of stochastic equations.  Math. Comp., 69(229):237–267, 2000.
  * [56] R. Munos.  A study of reinforcement learning in the continuous case by the means of viscosity solutions.  Machine Learning, 40:265–299, 2000.
  * [57] R. Murray and M. Palladino.  A model for system uncertainty in reinforcement learning, 2018.
  * [58] S. Perrin, M. Laurière, J. Pérolat, M. Geist, R. Élie, and O. Pietquin.  Mean Field Games Flock! The Reinforcement Learning Way.  arXiv:2105.07933, 2021.
  * [59] 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.  Thirty-fourth Conference on Neural Information Processing Systems (NeurIPS 2020), December 2020.
  * [60] R. S. Sutton and A. G. Barto.  Reinforcement learning: an introduction.  Adaptive Computation and Machine Learning. MIT Press, Cambridge, MA, second edition, 2018.
  * [61] H. Wang, T. Zariphopoulou, and X. Y. Zhou.  Reinforcement learning in continuous time and space: A stochastic control approach.  Journal of Machine Learning Research, 21(198):1–34, 2020.
  * [62] T. Z. Xin Guo, Renyuan Xu.  Entropy Regularization for Mean Field Games with Learning.  arXiv:2010.00145, 2020.
  * [63] J. Yong and X. Y. Zhou.  Stochastic controls, volume 43 of Applications of Mathematics (New York).  Springer-Verlag, New York, 1999.  Hamiltonian systems and HJB equations.