# Decoding Mean Field Games from Population and Environment Observations By
Gaussian Processes

Jinyan Guo1, Chenchen Mou2, Xianjin Yang3,∗, Chao Zhou1 1Department of
Mathematics and Risk Management Institute, National University of Singapore,
Singapore<EMAIL_ADDRESS>2Department of Mathematics, City University of
Hong Kong, Hong Kong SAR, China<EMAIL_ADDRESS>3Department of
Computing and Mathematical Sciences, California Institute of Technology,
Pasadena, CA 91125, USA<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

This paper presents a Gaussian Process (GP) framework, a non-parametric
technique widely acknowledged for regression and classification tasks, to
address inverse problems in mean field games (MFGs). By leveraging GPs, we aim
to recover agents’ strategic actions and the environment’s configurations from
partial and noisy observations of the population of agents and the setup of
the environment. Our method is a probabilistic tool to infer the behaviors of
agents in MFGs from data in scenarios where the comprehensive dataset is
either inaccessible or contaminated by noises.

∗Corresponding author.

##  1\. Introduction

In this paper, we present a Gaussian Process (GP) framework to infer
strategies and environment configurations in mean field games (MFGs) from
partially noisy observations of agent populations and environmental setups.
MFGs, as proposed in [33, 34, 35, 29, 27, 28, 26], study the actions of a
large number of indistinguishable rational agents. As the number of agents
approaches infinity, the Nash equilibrium of a standard MFG is described by a
coupled system of two distinct partial differential equations (PDEs): the
Hamilton–Jacobi–Bellman (HJB) and the Fokker–Plank (FP) equations. The HJB
determines the agent’s value function, while the FP governs the evolution of
the distribution of agents. We refer the readers to [24, 21, 38, 20, 22, 19,
37] for applications of MFGs in various fields.

To date, the consistency and completeness of MFGs have been explored
thoroughly across different contexts. The earliest findings can be traced back
to the foundational studies by Lasry and Lions, and further expanded upon in
Lions’ lectures at Collège de France (referenced by [40]). Notably, only a few
MFG models provide explicit solutions, highlighting the significance of
numerical calculations in MFGs.

A typical time-dependent MFG admits the following form

$\displaystyle\begin{cases}-\partial_{t}u(t,x)-\nu\Delta
u+H(x,Du)=F[m](x),&\forall(t,x)\in(0,T)\times\mathbb{R}^{d},\\\
\partial_{t}m(t,x)-\nu\Delta
m-\operatorname{div}(mD_{p}H(x,Du))=0,&\forall(t,x)\in(0,T)\times\mathbb{R}^{d},\\\
m(0,x)=\mu(x),u(T,x)=g(x),&\forall x\in\mathbb{R}^{d},\end{cases}$ (1.1)

where $T$ is the terminal time, $u$ stands for a representative agent’s value
function, $m$ is the probability density function representing the
distribution of agents, $\nu$ denotes the volatility of agents’ movements, $H$
is referred as the Hamiltonian, $F$ is the mean field term describing the
effect of agents’ interactions, $\mu$ is the initial distribution, and $g$
represents the terminal cost. If agents play a MFG in (1.1), the optimal
strategy of agents is given by $-D_{p}H(x,Du)$ [33, 34, 35].

In applications, MFGs have emerged as a prominent framework for modeling and
analyzing the strategic interactions of a large number of agents.
Traditionally, a thorough understanding of these interactions requires
comprehensive data on individual agent strategies, population dynamics, and
the intricacies of the environment in which these agents operate. However, in
many real-world scenarios, our observations are often partial, noisy, or even
sporadic. The challenge then is to distill meaningful information and insights
from such limited data.

In this research, we focus on deducing agents’ strategies and environmental
information by observing the agents’ distribution and leveraging limited data
about the surrounding environment. Specifically, we are concerned with the
following problem.

###### Problem 1.

Suppose that agents are engaged in a MFG within a competitive environment (For
instance, suppose that agents play a MFG specified in (1.1)). Based on limited
and noisy observations of the agents’ population ($m$ in (1.1)) and partial
observations on the environment ($\nu,H,F$, and $g$ in (1.1)), we infer the
complete population (values of $m$), the optimal strategy employed by the
agents ($-D_{p}H(x,Du)$ in (1.1)), and the setup of the environment ($\nu$,
$H$, $F$, and $g$ in (1.1)).

We further explore Problem 1 by delving into the realms of both general
stationary and time-dependent MFGs, as detailed in Sections 3 and 4
respectively. Stationary MFGs’ inverse problems are of interest due to their
representation as asymptotic limits of time-dependent MFGs [7, 16, 9]. Agents
typically reach a steady state rapidly under appropriate conditions [7, 16,
9]. Therefore, in practical scenarios, observed MFG data is likely in this
steady state, emphasizing the need to recover strategies and environmental
information at equilibrium.

Inverse problems in MFGs often suffer from non-uniqueness, implying multiple
PDE profiles can fit the same data. Particularly, observing only the
distribution of agents rarely yields a unique profile for quantities
$(m,u,\nu,H,F,\text{ and }g)$ as defined in Problem 1. To motivate the
solvability of Problem 1 and highlight the necessity of observations on the
environment, we consider a scenario where agents engage in a MFG described in
(1.1), with a Hamiltonian expressed as $H(x,Du)=\frac{|Du(x)|^{2}}{2}+V(x)$.
In this setting, the optimal strategy for a typical agent is represented by
$-Du$ [33, 34, 35]. First, we consider the case when we know the viscosity
$\nu=1$. If we observe the aggregate population, i.e., we have access to the
population density $m$, we can uniquely determine the velocity field $-Du$
across the regions where $m$ is supported. This is evident from the FP
equation in (1.1). To illustrate, let $u_{1}$ be a solution to (1.1). Assume
there exists another value function, $u_{2}$, fulfilling the condition
$\partial_{t}m(t,x)-\Delta m-\operatorname{div}(mDu_{2})=0$. By subtracting
this equation from the FP equation in (1.1), we obtain
$\operatorname{div}(m(Du_{1}-Du_{2}))=0$. Further, upon multiplying this
resultant equation by $u_{1}-u_{2}$ and applying integration by parts, it
follows that $\int m|Du_{1}-Du_{2}|^{2}\operatorname{d\\!}x=0$. Consequently,
we deduce that $Du_{1}=Du_{2}$ almost everywhere with respect to the measure
$m$. Therefore, even without explicit knowledge of the potential function $V$
in the Hamiltonian or the specific interactions among agents, as described by
$F$, it is possible to ascertain the strategies $-Du$ of agents solely by
observing the population’s distribution. However, when the viscosity $\nu$ is
not known, it is impossible to recover a unique pair $(\nu,-Du)$ just only
from the Fokker–Plank equation. Indeed, given $m$, if $(\nu,u)$ solve
$\partial_{t}m(t,x)-\nu\Delta m-\operatorname{div}(mDu)=0$, $(\nu+C,u-C\log
m)$ is also a solution for any constant $C$. Therefore, to uniquely determine
the optimal strategy of agents, in addition to the knowledge of the density
$m$, we need to know either the viscosity $\nu$ or some information about the
Hamiltonian and the coupling function $F$ in the HJB equation in (1.1). Thus,
without any information about the environment characterized by $(\nu,H,F,g)$,
it is unlikely to determine the optimal strategies of agents uniquely.

Problem 1 expands upon the scenarios outlined in the preceding paragraph. It
delves into the recovery of the strategy, population, and environment through
partial and noisy observations of both the population and the environment. We
tackle Problem 1 through Gaussian Processes (GPs). a powerful non-parametric
statistical method renowned for its ability to model and predict complex
phenomena with a degree of uncertainty. The beauty of GPs lies in their
inherent flexibility, allowing them to capture a wide range of behaviors and
patterns without being tied down by a strict parametric form. This makes them
an appealing choice for tackling the multifaceted challenges posed by MFGs. We
provide a general GP framework which determines the unknown variables using
Maximum A Posteriori (MAP) estimators. These estimators are conditioned on the
fulfillment of PDEs at designated sample points. Additionally, the process
incorporates regularization through the inclusion of observations that are
subject to noises.

In Section 5, we demonstrate the effectiveness of our methodology through a
variety of examples. Our tests focus on reconstructing comprehensive profiles
of populations, strategies of agents, and complete environmental setups from
partial and noisy data about populations and environments. In particular,
Subsection 5.1 details the calculation of the viscosity constant, an essential
parameter for gauging uncertainty in dynamic systems. Subsections 5.2 and 5.3
are devoted to identifying the coupling functions in MFG, which reflect the
interactions between agents. In Subsection 5.4, we address a misspecification
problem. This section demonstrates the robustness and effectiveness of our
methodology by recovering agents’ information using a misspecified model. The
data originates from a MFG characterized by local interactions, yet our model
employs a non-local coupling approach for this reconstruction. This contrast
underscores the adaptability and precision of our technique in handling
complex data scenarios. Finally, in Subsection 5.5, we apply our methods to a
time-dependent MFG, highlighting its flexibility in time-dependent
environments.

###  1.1. Related Works in MFGs

Although numerous numerical algorithms exist for solving MFGs (see [1, 3, 48,
43, 44, 6, 5, 10, 11, 52, 39, 36, 2, 23, 47, 45]), there is a notably limited
amount of research addressing the inverse problems in MFGs [41, 17, 15, 32,
31, 30, 42, 18]. To the best of our knowledge, only three numerical papers
have been published that specifically address the inverse problems in MFGs.
The work [17] introduces models to reconstruct ground metrics and interaction
kernels in MFG running costs. Efficiency and robustness of the models are
demonstrated through numerical simulations. The paper [15] presents a
numerical algorithm to solve an inverse problem in MFGs using partial boundary
measurements. Authors in [32] propose a globally convergent convexification
numerical method for recovering the global interaction term from the single
measurement data. In a recent study, [54] tackles mean field control inverse
problems through operator learning, introducing a novel algorithm designed for
broader applications in operator learning. Their approach, as detailed in
[54], involves training models on pairs of input and output data, all
calculated under the same parameters. Contrasting with this method, our paper
presents a distinct approach where we estimate the unknowns by analyzing the
data in a single observation.

In contrast to previous research on inverse problems in MFGs, the GP framework
introduced in this paper offers a comprehensive approach to simultaneously
infer all unknown variables using only noisy, partial observations of the
population and environment. Notably, to the best of our knowledge, this is the
first study to estimate the viscosity parameter $\nu$ in MFGs. This parameter
is vital as it reflects the uncertainty in agents’ dynamics.

###  1.2. Related Works in GPs

GP Regression, a Bayesian non-parametric technique for supervised learning,
stands out for its ability to quantify uncertainty. GPs have been successfully
applied in solving and learning ODEs [25, 55] and PDEs [13, 14, 56, 50, 51].
In the context of MFGs, GPs have also been deployed to solve MFG systems, see
[47, 45].

This paper extends the earlier work of learning unknowns in PDEs using GPs in
[13] for inverse problems in MFGs. The key difference between (1.1) and a
general PDE system lies in the fact that (1.1) consists of a forward and
backward PDE system, where $m$ is given the initial condition and $u$ admits
the terminal condition. Hence, the time discretization scheme presented in
[13] can not be simply applied in the MFG setting. Instead, in Section 4, we
adopt a forward-backward discretization which is inspired by the finite
difference method proposed in [1] for solving time-dependent inverse problems
of MFGs.

###  1.3. Outlines

The structure of this paper is as follows. In Section 2, we provide a review
of the fundamental concepts of GP regression for the sake of readers. In
Section 3, we delve into a comprehensive analysis of Problem 1, focusing on
its application in stationary MFGs. This is followed by formulating the
inverse problem of recovering unknown elements in stationary MFGs as an
optimal recovery problem. There, we establish the existence of minimizers.
Progressing to Section 4, we introduce a method of time-discretization to
address inverse problems in time-dependent MFGs. The efficacy of our proposed
techniques is validated through various numerical experiments, as presented in
Section 5. The paper concludes in Section 6, where we engage in further
discussion and discuss potential future work.

###### Notations A0.

In this paper, a real-valued vector $\boldsymbol{v}$ is consistently displayed
in boldface, except when it denotes a point in the physical domain. The
Euclidean norm of $\boldsymbol{v}$ is denoted by $|\boldsymbol{v}|$, while its
transpose is represented as $\boldsymbol{v}^{T}$. For a vector
$\boldsymbol{v}$, $v_{i}$ represents its $i^{\textit{th}}$ element. For a
matrix $M$, we denote by $|M|$ the Frobenius norm. Given a function $u$ and a
vector $\boldsymbol{v}$, $u(\boldsymbol{v})$ denotes the vector
$(u(v_{1}),\dots,u(v_{N}))$, with $N$ being the length of $\boldsymbol{v}$. We
denote by $\delta_{x}$ the Dirac delta function concentrated at $x$. Let
$\Omega$ be a subset of $\mathbb{R}^{d}$. The terms $\operatorname{int}\Omega$
and $\partial\Omega$ refer to the interior and boundary of $\Omega$,
respectively. We denote by $\mathcal{N}(0,\gamma^{2}I)$ the standard
multivariate normal distribution with the covariance matrix given by
$\gamma^{2}I$, where $\gamma>0$. For a normed vector space labeled $V$, its
norm is represented by $\|\cdot\|_{V}$.

Let’s consider a Banach Space $\mathcal{U}$ associated with a quadratic norm
$\|\cdot\|_{\mathcal{U}}$. The dual space of $\mathcal{U}$ is represented as
$\mathcal{U}^{*}$, and we use the notation $[\cdot,\cdot]$ for the duality
pairing. We propose the existence of a covariance operator,
$\mathcal{K}_{\mathcal{U}}\mathrel{\mathop{\mathchar
58\relax}}\mathcal{U}^{*}\rightarrow\mathcal{U}$, which is linear, bijective,
symmetric (meaning
$[\mathcal{K}_{\mathcal{U}}\phi,\psi]=[\mathcal{K}_{\mathcal{U}}\psi,\phi]$),
and maintains positivity ($[\mathcal{K}_{\mathcal{U}}\phi,\phi]>0$ whenever
$\phi\neq 0$). The norm in $\mathcal{U}$ is defined as
$\|u\|_{\mathcal{U}}^{2}=[\mathcal{K}^{-1}_{\mathcal{U}}u,u],\text{for every
}u\in\mathcal{U}$. Let
$\boldsymbol{\phi}=(\phi_{1},\dots,\phi_{P}),P\in\mathbb{N}$, which is
situated within the product space ${(\mathcal{U}^{*})}^{\bigotimes P}$. For
any element $u$ in $\mathcal{U}$, we express the pairing
$[\boldsymbol{\phi},u]$ as $[\boldsymbol{\phi},u]\mathrel{\mathop{\mathchar
58\relax}}=([{\phi}_{1},u],\dots,[{\phi}_{P},u])$.

Let $\\{\mathcal{U}_{i}\\}_{i=1}^{N}$ be a collection of RKHSs and define the
product space
$\boldsymbol{\mathcal{U}}=\mathcal{U}_{1}\times\dots\times\mathcal{U}_{N}$. We
define the dual $\boldsymbol{\mathcal{U}}^{*}$ of $\boldsymbol{\mathcal{U}}$
as the product space
$\mathcal{U}^{*}_{1}\times\dots\times\mathcal{U}_{N}^{*}$. Let
$\boldsymbol{\phi}=(\phi_{1},\dots,\phi_{N})\in\boldsymbol{\mathcal{U}}^{*}$
and $\boldsymbol{u}=(u_{1},\dots,u_{N})\in\boldsymbol{\mathcal{U}}$, we define
$\boldsymbol{\phi}\odot\boldsymbol{u}=([\phi_{1},u_{1}],\dots,[\phi_{N},u_{N}])$.
Let $\Phi$ be a matrix of linear operators with the $i^{\textit{th}}$ column,
denoted by $\boldsymbol{\phi}_{i}$, in $\boldsymbol{\mathcal{U}}^{*}$. For
$\boldsymbol{u}=(u_{1},\dots,u_{N})\in\mathbb{\mathcal{U}}$, we represent the
matrix $[\Phi,\boldsymbol{u}]$ such that the $i^{\textit{th}}$ column is
defined as $\boldsymbol{\phi}_{i}\odot\boldsymbol{u}$. For a matrix
$\boldsymbol{Y}$, let $\overset{\rightarrow}{\boldsymbol{Y}}$ be the vector
formed by concatenating ${\boldsymbol{Y}}$’s columns. The symbol $C$ stands
for a variable positive real number which can vary from one instance to
another.

##  2\. Prerequisites for GP Regression

In this section, we briefly revisit the fundamentals of GP Regression and its
connection to Reproducing Kernel Hilbert Spaces (RKHSs). We refer readers to
[53, 49] for an in-depth exploration of the subject.

###  2.1. Learning Scalar-valued Functions with GPs

Let $\Omega\subseteq\mathbb{R}^{d}$ be an open subset. A real-valued GP,
$f\mathrel{\mathop{\mathchar 58\relax}}\Omega\to\mathbb{R}$, is defined such
that for any finite set $\boldsymbol{x}$ of points in $\Omega$,
$f(\boldsymbol{x})$ follows a joint Gaussian distribution. A GP $f$ is usually
characterized by a mean function $\mu\mathrel{\mathop{\mathchar
58\relax}}\Omega\to\mathbb{R}$ and a covariance function
$K\mathrel{\mathop{\mathchar 58\relax}}\Omega\times\Omega\to\mathbb{R}$ such
that $E(f(x))=\mu(x)$ and
$\operatorname{Cov}(f(x),f(x^{\prime}))=K(x,x^{\prime}),\forall
x,x^{\prime}\in\Omega$. In this case, we denote $f\sim\mathcal{GP}(\mu,K)$. In
supervised learning, we aim to build a GP estimator $f^{\dagger}$ from a
training set $(x_{i},y_{i})_{i=1}^{N}$. To do that, we assume
$f\sim\mathcal{GP}(0,K)$ and define the estimator $f^{\dagger}$ as the mean of
the posterior of $f$ conditioned on the training set, i.e.,

$\displaystyle f^{\dagger}=E[f|f(x_{i})=y_{i},i=1,\dots,N].$ (2.1)

The estimator $f^{\dagger}$ admits an explicit formula given by

$\displaystyle
f^{\dagger}(x)=K(x,\boldsymbol{x})K(\boldsymbol{x},\boldsymbol{x})^{-1}\boldsymbol{y},$
(2.2)

where $\boldsymbol{x}=(x_{i})_{i=1}^{N},\boldsymbol{y}=(y_{i})_{i=1}^{N}$,
$K(x,\boldsymbol{x})=(K(x,x_{i}))_{i=1}^{N}$, and
$K(\boldsymbol{x},\boldsymbol{x})$ is a matrix with entries $K(x_{i},x_{j})$.
The estimator $f^{\dagger}$ can also be understood from an optimization
perspective. Let $K$ be a positive definite covariance function. Then,
according to the Moore–Aronszajn theorem, there is a RKHS $\mathcal{U}$
associated with $K$, such that $f(x)=\langle f,K(x,\cdot)\rangle$. Then, the
estimator $f^{\dagger}$ in (2.1) and (2.2) is a minimizer of the following
optimal recover problem

$\displaystyle\begin{cases}\min\limits_{f\in\mathcal{U}}\|f\|_{\mathcal{U}}^{2}\\\
\operatorname{s.t.}f(x_{i})=y_{i},\forall i\in\\{1,\dots,N\\}.\end{cases}$
(2.3)

We call (2.2) the representer formula for the optimization problem (2.3) and
it is the cornerstone of GP Regression. Let $\mathcal{U}^{*}$ be the dual of
$\mathcal{U}$ and denote by $[,]$ the duality pairing. Let $\delta_{x}$ bet
the Dirac function centered at $x$. Then, $\delta_{x}\in\mathcal{U}^{*}$ for
$x\in\Omega$. Define
$\boldsymbol{\delta}=(\delta_{x_{1}},\dots,\delta_{x_{N}})$. Then, for
simplicity, we write (2.3) as

$\displaystyle\begin{cases}\min\limits_{f\in\mathcal{U}}\|f\|_{\mathcal{U}}^{2}\\\
\operatorname{s.t.}[\boldsymbol{\delta},f]=\boldsymbol{y}.\end{cases}$

For a kernel $K$ and a collection $\boldsymbol{\phi}$ of linear operators in
the dual space of the RKHS associated with $K$, let $K(x,\boldsymbol{\phi})$
be the vector with entries
$\int_{\Omega}K(x,x^{\prime})\phi_{i}(x^{\prime})\operatorname{d\\!}x^{\prime}$
and let $K(\boldsymbol{\phi},\boldsymbol{\phi})$, called the Gram matrix, be
with entries
$\int_{\Omega}\int_{\Omega}K(x,x)\phi_{i}(x^{\prime})\phi_{j}(x^{\prime})\operatorname{d\\!}x\operatorname{d\\!}x^{\prime}$.
Then, we reformulate $f^{\dagger}$ in (2.2) as

$\displaystyle
f^{\dagger}(x)=K(x,\boldsymbol{\delta})K(\boldsymbol{\delta},\boldsymbol{\delta})^{-1}\boldsymbol{y}.$

Next, we rewrite the training set $(x_{i},y_{i})_{i=1}^{N}$ as
$(\delta_{x_{i}},y_{i})_{i=1}^{N}$. Then, if we replace the training set
$(\delta_{x_{i}},y_{i})_{i=1}^{N}$ with $(\phi_{i},y_{i})_{i=1}^{N}$, where
$\\{\phi_{i}\\}_{i=1}^{N}$ contains other independent linear operators in
$\mathcal{U}^{*}$, the mean of the posterior of the prior $f$,
$f\sim\mathcal{GP}(0,K)$, conditioned on the training set
$(\phi_{i},y_{i})_{i=1}^{N}$ is given by [49, Sec 17.8]

$\displaystyle
f^{\dagger}(x)=K(x,\boldsymbol{\phi})K(\boldsymbol{\phi},\boldsymbol{\phi})^{-1}\boldsymbol{y},$
(2.4)

which is the minimizer of the following optimal recovery problem

$\displaystyle\begin{cases}\min\limits_{f\in\mathcal{U}}\|f\|_{\mathcal{U}}^{2}\\\
\operatorname{s.t.}[\boldsymbol{\phi},f]=\boldsymbol{y}.\end{cases}$

###  2.2. Learning Vector-valued Functions with GPs

The above discussion can be extended to consider regressing vector-valued
functions with GPs. We refer readers to [4, 46] for further details.

Let $\Omega\subseteq\mathbb{R}^{d}$ be an open subset. A vector-valued GP,
$\boldsymbol{f}\mathrel{\mathop{\mathchar 58\relax}}\Omega\to\mathbb{R}^{m}$,
is defined such that for any $\boldsymbol{X}\in\Omega^{\bigotimes N}$,
$\boldsymbol{f}(\boldsymbol{X})\in\mathbb{R}^{N\times m}$ follows a joint
Gaussian distribution. A vector-valued GP $\boldsymbol{f}$ is characterized by
a mean function $\boldsymbol{\mu}\mathrel{\mathop{\mathchar
58\relax}}\Omega\to\mathbb{R}^{m}$ and a covariance function
$K\mathrel{\mathop{\mathchar
58\relax}}\Omega\times\Omega\to\mathbb{R}^{m\times m}$ such that
$E(\boldsymbol{f}(x))=\boldsymbol{\mu}(x)$ and
$\operatorname{Cov}(\boldsymbol{f}(x),\mathbf{f}(x^{\prime}))=K(x,x^{\prime})$,
for all $x,x^{\prime}\in\Omega$. We denote this as
$\boldsymbol{f}\sim\mathcal{GP}(\boldsymbol{\mu},K)$. In the setting of
learning vector-valued functions, our goal is to construct a GP estimator
$\boldsymbol{f}^{\dagger}$ from a training set
$(x_{i},\boldsymbol{Y}_{i})_{i=1}^{N}$, where
$\boldsymbol{Y}_{i}\in\mathbb{R}^{m}$. We assume
$\boldsymbol{f}\sim\mathcal{GP}(\boldsymbol{0},K)$ and define the estimator
$\boldsymbol{f}^{\dagger}$ as the mean of the posterior of $\boldsymbol{f}$
conditioned on the training set, i.e.,

$\displaystyle\boldsymbol{f}^{\dagger}=E[\boldsymbol{f}|\boldsymbol{f}(x_{i})=\boldsymbol{Y}_{i},i=1,\dots,N].$

Let $\boldsymbol{Y}$ be the matrix with the $i^{\textit{th}}$ column being
$\boldsymbol{Y}_{i}$, and let $\overset{\rightarrow}{\boldsymbol{Y}}$ bet the
vector obtained by concatenating the columns of $\boldsymbol{Y}$. The
estimator $\boldsymbol{f}^{\dagger}$ is explicitly given by

$\displaystyle\boldsymbol{f}^{\dagger}(x)=K(x,\boldsymbol{x})K(\boldsymbol{x},\boldsymbol{x})^{-1}\overset{\rightarrow}{\boldsymbol{Y}},$

where $K(x,\boldsymbol{x})$ is a matrix with $m$ rows and $N\times m$ columns
by concatenating matrices $K(x,x_{i})$ for $i=1,\dots,N$, and
$K(\boldsymbol{x},\boldsymbol{x})$ is defined as

$\displaystyle
K(\boldsymbol{x},\boldsymbol{x}^{\prime})=\begin{bmatrix}K(x_{1},x_{1})&\dots&K(x_{1},x_{N})\\\
\dots&\dots&\dots\\\ K(x_{N},x_{1})&\dots&K(x_{N},x_{N})\end{bmatrix}.$

The estimator $\boldsymbol{f}^{\dagger}$ also has an optimization
interpretation. Let $\boldsymbol{\mathcal{U}}$ be the vector-valued RKHS
associated with $K$, where
$\boldsymbol{f}(x)^{T}\boldsymbol{c}=\langle\boldsymbol{f},K(\cdot,x)\boldsymbol{c}\rangle$,
for any $\boldsymbol{c}\in\mathbb{R}^{m}$ and
$\boldsymbol{f}\in\boldsymbol{\mathcal{U}}$. The estimator
$\boldsymbol{f}^{\dagger}$ minimizes the following optimal recovery problem

$\displaystyle\begin{cases}\min\limits_{\boldsymbol{f}\in\mathcal{U}}\|\boldsymbol{f}\|_{\mathcal{U}}^{2}\\\
\operatorname{s.t.}\boldsymbol{f}(x_{i})=\boldsymbol{Y}_{i},\forall
i\in\\{1,\dots,N\\}.\end{cases}$

In this paper, we are interested in the special case where each output of
$\boldsymbol{f}$ can be modeled as an independent GP, i.e., $K(x,x^{\prime})$
is a diagonal matrix for all $x,x^{\prime}\in\Omega$. Denote
$K(x,x^{\prime})=\operatorname{diag}(K_{1}(x,x^{\prime}),\dots,K_{m}(x,x^{\prime}))$.
For simplicity, we represent $\boldsymbol{f}(x)=(f_{1}(x),\dots,f_{m}(x))$ and
$f_{i}\in\mathcal{U}_{i}$ with $\mathcal{U}_{i}$ being a RKHS associated with
the kernel $K_{i}$ for $i=1,\dots,m$. Hence, the RKHS
$\boldsymbol{\mathcal{U}}$ is equivalent to
$\mathcal{U}_{1}\times\dots\times\mathcal{U}_{m}$. Accordingly,
$\|\boldsymbol{f}\|_{\boldsymbol{\mathcal{U}}}^{2}=\sum_{i=1}^{N}\|f_{i}\|_{\mathcal{U}_{i}}^{2},\forall\boldsymbol{f}\in\boldsymbol{\mathcal{U}}$.

Let
$\boldsymbol{\mathcal{U}}^{*}=\mathcal{U}_{1}^{*}\times\dots\times\mathcal{U}_{m}^{*}$.
Consider the training set with
$(\boldsymbol{\phi}_{i},\boldsymbol{Y}_{i})_{i=1}^{N}$, where
$\boldsymbol{\phi}_{i}$ is a vector of linear operators with
$\phi_{i,j}\in\mathcal{U}_{j}^{*}$. Let $\Phi$ be the matrix with the
$i^{\textit{th}}$ column being $\boldsymbol{\phi}_{i}$. Let
$\boldsymbol{\phi}\mathrel{\mathop{\mathchar
58\relax}}=\overset{\rightarrow}{\Phi}$, the vector concatenating
$(\boldsymbol{\phi}_{i})_{i=1}^{N}$. Then, the posterior mean of the prior
$\boldsymbol{f}$, assumed as
$\boldsymbol{f}\sim\mathcal{GP}(\boldsymbol{0},K)$, is given by
$\boldsymbol{f}^{\dagger}=E[\boldsymbol{f}|[\Phi,\boldsymbol{f}]=\boldsymbol{Y}]$
admits the following representer formula

$\displaystyle\boldsymbol{f}^{\dagger}(x)=K(x,\boldsymbol{\phi})K(\boldsymbol{\phi},\boldsymbol{\phi})^{-1}\overset{\rightarrow}{\boldsymbol{Y}},$
(2.5)

where $\overset{\rightarrow}{\boldsymbol{Y}}$ is the vector concatenating
columns of $\boldsymbol{Y}$ and $\boldsymbol{f}^{\dagger}$ minimizes the
following optimal recovery problem:

$\displaystyle\begin{cases}\min\limits_{\boldsymbol{f}\in\boldsymbol{\mathcal{U}}}\|\boldsymbol{f}\|_{\boldsymbol{\mathcal{U}}}^{2}\\\
\operatorname{s.t.}[\Phi,\boldsymbol{f}]=\boldsymbol{Y}.\end{cases}$

##  3\. An Inverse Problem for Stationary MFGs

In this section, we elaborate Problem 1 in the context of stationary MFGs.
Inverse problems in stationary MFGs are intriguing as they represent the long-
term limits of time-dependent MFGs [7, 16, 9]. Under certain conditions,
agents’ states converge exponentially to a steady state [7, 16, 9]. Hence, MFG
data is likely to reflect this steady state, making it practical to focus on
recovering strategies and environmental information in this equilibrium. Then,
we introduce a GP framework to solve the inverse problem. This approach builds
upon the GP method from [13], adapting it to coupled PDE systems with integral
constraints. Subsection 3.1 redefines Problem 1 for stationary MFGs,
Subsection 3.2 formulates the inverse problem as an optimal recovery problem
and provides a general framework for solving inverse problems of stationary
MFGs. To illustrate the framework, we describe our method through a concrete
example in Subsection 3.3. Subsection 3.4 proves the existence of a minimizer,
confirming the solvability of our proposed inverse problem. Subsection 3.5
offers a probabilistic interpretation of our method, linking it to maximum
likelihood estimation.

###  3.1. An Inverse Problem for Stationary MFGs

Let $\Omega$ be a subset of $\mathbb{R}^{d}$. Suppose that the stationary MFGs
of our interests admit the following form

$\displaystyle\begin{cases}\boldsymbol{\mathcal{P}}(u^{*},m^{*},\overline{H}^{*},\boldsymbol{V}^{*})(x)=0,&\forall
x\in\operatorname{int}\Omega,\\\
\mathcal{B}(u^{*},m^{*},\boldsymbol{V}^{*})(x)=0,&\forall
x\in\partial\Omega,\\\
\int_{\Omega}u^{*}\operatorname{d\\!}x=0,\int_{\Omega}m^{*}\operatorname{d\\!}x=1,\end{cases}$
(3.1)

where $\boldsymbol{V}^{*}$ contains a set of functions, which represent
environment configurations. Given $\boldsymbol{V}^{*}$,
$(u^{*},m^{*},\overline{H}^{*})$ solves the stationary MFG (3.1) such that
$u^{*}$ stands for the value function, $m^{*}$ represents the density of
populations of agents, and $\overline{H}^{*}$ is real number ensuring the unit
mass of $m^{*}$. In this context, $\boldsymbol{\mathcal{P}}$ is a collection
of nonlinear differential operators, while $\mathcal{B}$ stands for a boundary
operator. We proceed under the assumption that equation (3.1) has a unique
classical solution $(u^{*},m^{*},\overline{H}^{*})$ when the underlying
environment configuration $\boldsymbol{V}^{*}$ is given. The equation (3.1)
represents a general form of the long-time limit of the time-dependent MFG
described in (1.1). For discussions on the long-time behavior of MFGs, see [7,
16, 9]. We detail Problem 1 for stationary MFGs of the form (3.1) in the
following problem.

###### Problem 2.

Let $\boldsymbol{V}^{*}=\\{V^{*}_{i}\\}_{i=1}^{N_{v}}$ be a collection of
functions on $\Omega$, where $N_{v}\in\mathbb{N}$. Assume that given
$\boldsymbol{V}^{*}$, (3.1) admits a unique classical solution
$(u^{*},m^{*},\overline{H}^{*})$. Suppose that we only have some partial noisy
observations on $m^{*}$ and on a subset $\boldsymbol{V}_{o}^{*}$ of
$\boldsymbol{V}^{*}$. The objective is to infer the values of $u^{*}$,
$m^{*}$, $\overline{H}^{*}$, and $\boldsymbol{V}^{*}$. To elaborate, we have

  1. 1.

Partial noisy observations on $m^{*}$. We have a set of linear operators
$\\{\phi_{l}^{o}\\}_{l=1}^{N_{m}}$, $N_{m}\in\mathbb{N}$ where some noisy
observations on $m^{*}$ are available. These observations are represented as
$\boldsymbol{m}^{o}$, i.e.
$\boldsymbol{m}^{o}=([\phi_{1}^{o},m^{*}],\dots,[\phi_{N_{m}}^{o},m^{*}])+\boldsymbol{\epsilon}$,
$\boldsymbol{\epsilon}\sim\mathcal{N}(0,\gamma^{2}I)$. We call
$\boldsymbol{\phi}^{o}=(\phi_{1}^{o},\dots,\phi_{N_{m}}^{o})$ the vector of
observation operators. For instance, if we only have observations of $m^{*}$
at a finite set of collocation points, $\boldsymbol{\phi}^{o}$ contains Diracs
centered at these points.

  2. 2.

Partial noisy observations of $\boldsymbol{V}^{*}$. Let
$\boldsymbol{V}_{o}^{*}=\\{V_{j}^{*}\\}_{j\in\mathcal{O}_{v}}$ with
$\mathcal{O}_{v}\subset\\{1,\dots,N_{v}\\}$ be a subset of
$\boldsymbol{V}^{*}$. For this subset, we possess noisy observations at
collections of linear operators
$\\{\\{\psi^{o,j}_{e}\\}_{e=1}^{N^{j}_{v}}\\}_{j\in\mathcal{O}_{v}}$. These
observations can be collected into a matrix $\boldsymbol{V}^{o}$ such that the
$j^{\textit{th}}$ row of $\boldsymbol{V}^{o}$, denoted by
$\boldsymbol{V}_{j}^{o}$ is given by

$\displaystyle\boldsymbol{V}^{o}_{j}=([\psi^{o,j}_{1},V_{j}],\dots,[\psi^{o,j}_{N_{v}^{j}},V_{j}])+\boldsymbol{\epsilon}_{j},\boldsymbol{\epsilon}_{j}\sim\mathcal{N}(0,\gamma_{j}^{2}I).$

Inverse Problem Statement
Suppose that agents are involved in the stationary MFG in (3.1), we infer
$u^{*}$, $m^{*}$, $\overline{H}^{*}$, and $\boldsymbol{V}^{*}$ based on
$\boldsymbol{m}^{o}$ and $\boldsymbol{V}^{o}$.

###  3.2. The Optimal Recovery Problem

To solve Problem 2, we approximate
$(u^{*},m^{*},\overline{H}^{*},\boldsymbol{V}^{*})$ by GPs conditioned on PDEs
at sampled collocation points in $\Omega$. Then, we compute the solution by
calculating the MAP points of such conditioned GPs. More precisely, we take a
set of samples $\\{x_{i}\\}_{i=1}^{M}$ in such a way that
$x_{1},\dots,x_{M_{\Omega}}\in\operatorname{int}\Omega$ and
$x_{M_{\Omega}+1},\dots,x_{M}\in\partial\Omega$ for $1\leqslant
M_{\Omega}\leqslant M$. We put GP priors on $u^{*},m^{*},\overline{H}^{*}$ and
$\boldsymbol{V}^{*}$ by assuming that $u^{*}\sim\mathcal{GP}(0,K_{u})$,
$m^{*}\sim\mathcal{GP}(0,K_{m})$, $\overline{H}^{*}\sim\mathcal{N}(0,1)$, and
$\boldsymbol{V}^{*}\sim\mathcal{GP}(\boldsymbol{0},K_{\boldsymbol{V}})$, where
$\boldsymbol{V}^{*}$ is a vector-valued GP with independent outputs. Let
$\mathcal{U}$, $\mathcal{M}$, and $\boldsymbol{\mathcal{V}}$ be RKHSs
associated with $K_{u},K_{m}$, and $K_{\boldsymbol{V}}$, separately. Denote
$\boldsymbol{\phi}^{o}=(\phi^{o}_{1},\dots,\phi_{{N_{m}}}^{o})$. Let
$\boldsymbol{m}^{o}$ be a noisy observations of
$[\boldsymbol{\phi}^{o},m^{*}]$ such that
$[\boldsymbol{\phi}^{o},m^{*}]-\boldsymbol{m}^{o}\in\mathcal{N}(0,\gamma^{2}I)$.
Let $\Psi$ be a matrix of linear functionals with each column in
$\boldsymbol{\mathcal{V}}^{*}$, the dual of $\boldsymbol{\mathcal{V}}$, such
that we have a noisy observation $\boldsymbol{V}^{o}$ of
$[\Psi,\boldsymbol{V}^{*}]$. Then, to solve Problem 2, we introduce a positive
penalization parameter $\beta$ and approximate
$(u^{*},m^{*},\overline{H}^{*},\boldsymbol{V}^{*})$ by a minimizer of the
following optimal recovery problem

$\displaystyle\begin{cases}\min\limits_{(u,m,\overline{H},\boldsymbol{V})\in\mathcal{U}\times\mathcal{M}\times\mathbb{R}\times\boldsymbol{\mathcal{V}}}\|u\|_{\mathcal{U}}^{2}+\|m\|_{\mathcal{M}}^{2}+|\overline{H}|^{2}+\|\boldsymbol{V}\|_{\boldsymbol{\mathcal{V}}}^{2}+\frac{1}{\gamma^{2}}|[\boldsymbol{\phi}^{o},m]-\boldsymbol{m}^{o}|^{2}+|\Sigma^{-1}([\Psi,\boldsymbol{V}]-\boldsymbol{V}^{o})|^{2}\\\
\quad\quad\quad\quad\quad\quad\quad+\beta\bigg{|}\frac{1}{M_{\Omega}}\sum\limits_{i=1}^{M_{\Omega}}u(x_{i})\bigg{|}^{2}+\beta\bigg{|}\frac{1}{M_{\Omega}}\sum\limits_{i=1}^{M_{\Omega}}m(x_{i})-1\bigg{|}^{2}\\\
\text{s.t.}\quad\boldsymbol{\mathcal{P}}(u,m,\overline{H},\boldsymbol{V})(x_{i})=0,\quad\text{for}\
i=1,\dots M_{\Omega},\\\ \quad\quad\
\mathcal{B}(u,m,\boldsymbol{V})(x_{j})=0,\quad\text{for}\
j=M_{\Omega+1},\dots,M,\end{cases}$ (3.2)

where $\Sigma=\operatorname{diag}(\gamma_{1},\dots,\gamma_{N_{v}})$. The
formulation in (3.2) extends the one in Section 4.4 of [13], adapting it to a
coupled PDE system with integral constraints. Our approach, which involves
discretizing the integral constraints and incorporating them into the
integration process, circumvents the need for kernel integration during
covariance matrix assembly. Should the solution of (3.1) lack sufficient
smoothness, we recommend employing the vanishing viscosity method [8], or
regularizing the MFG using smooth mollifiers [12]. This approach yields a
system with solutions of higher regularity, to which our numerical methods can
be applied.

The minimization in (3.2) involves optimization within function spaces. For
computational feasibility through numerical algorithms, we reformulate (3.2)
into a finite-dimensional minimization problem. To do that, we make the
following assumptions about the forms of $\boldsymbol{\mathcal{P}}$ and
$\mathcal{B}$.

###### Assumption 1.

Assume that there exists collections of bounded linear operators
$\boldsymbol{L}^{u,\Omega}$, $\boldsymbol{L}^{u,\partial\Omega}$,
$\boldsymbol{L}^{m,\Omega}$, $\boldsymbol{L}^{m,\partial\Omega}$,
$\boldsymbol{L}^{\boldsymbol{{V}},\Omega}$, and
$\boldsymbol{L}^{\boldsymbol{{V}},\partial\Omega}$, and continuous nonlinear
maps $\boldsymbol{{P}}$ and $B$ such that

$\displaystyle\begin{cases}\boldsymbol{\mathcal{P}}(u,m,\overline{H},\boldsymbol{V})(x)=\boldsymbol{{P}}(\boldsymbol{L}^{u,\Omega}(u)(x),\boldsymbol{L}^{m,\Omega}(m)(x),\overline{H},\boldsymbol{L}^{\boldsymbol{V},\Omega}(\boldsymbol{V})(x)),\forall
x\in\operatorname{int}\Omega,\\\
\mathcal{B}(u,m,\boldsymbol{V})(x)=B(\boldsymbol{L}^{u,\partial\Omega}(u)(x),\boldsymbol{L}^{m,\partial\Omega}(m)(x),\boldsymbol{L}^{\boldsymbol{V},\partial\Omega}(\boldsymbol{V})(x)),\forall
x\in\partial\Omega.\end{cases}$

In the following texts, we denote
$\boldsymbol{L}^{u}\mathrel{\mathop{\mathchar
58\relax}}=(\boldsymbol{L}^{u,\Omega},\boldsymbol{L}^{u,\partial\Omega})$,
$\boldsymbol{L}^{m}\mathrel{\mathop{\mathchar
58\relax}}=(\boldsymbol{L}^{m,\Omega},\boldsymbol{L}^{m,\partial\Omega})$, and
$\boldsymbol{L}^{\boldsymbol{{V}}}\mathrel{\mathop{\mathchar
58\relax}}=(\boldsymbol{L}^{\boldsymbol{{V}},\Omega},\boldsymbol{L}^{\boldsymbol{{V}},\partial\Omega})$.

Under Assumption 1, we define a functional vector
$\boldsymbol{\delta}\mathrel{\mathop{\mathchar
58\relax}}=(\boldsymbol{\delta}^{\Omega},\boldsymbol{\delta}^{\partial\Omega})$,
where
$\boldsymbol{\delta}^{\Omega}=(\delta_{x_{1}},\dots,\delta_{x_{M_{\Omega}}})$
and
$\boldsymbol{\delta}^{\partial\Omega}=(\delta_{x_{M_{\Omega}+1}},\dots,\delta_{x_{M}})$.
Let
$\widetilde{\boldsymbol{\phi}}^{u}=\boldsymbol{\delta}\circ\boldsymbol{L}^{u}$,
$\widetilde{\boldsymbol{\phi}}^{m}=\boldsymbol{\delta}\circ\boldsymbol{L}^{m}$,
and
$\widetilde{\boldsymbol{\phi}}^{\boldsymbol{V}}=\boldsymbol{\delta}\circ\boldsymbol{L}^{\boldsymbol{V}}$.
We define the nonlinear map $G$ such that for any $u\in\mathcal{U}$,
$m\in\mathcal{M}$, $\overline{H}\in\mathbb{R}$, and
$\boldsymbol{V}\in\boldsymbol{\mathcal{V}}$, we have

$\displaystyle(G([\widetilde{\boldsymbol{\phi}}^{u},u],[\widetilde{\boldsymbol{\phi}}^{m},m],\overline{H},[\widetilde{\boldsymbol{\phi}}^{\boldsymbol{V}},\boldsymbol{V}]))_{i}=\begin{cases}\boldsymbol{{P}}([\delta_{x_{i}}\circ\boldsymbol{L}^{\Omega,u},u],[\delta_{x_{i}}\circ\boldsymbol{L}^{\Omega,m},m],\overline{H},[\delta_{x_{i}}\circ\boldsymbol{L}^{\Omega,\boldsymbol{V}},\boldsymbol{V}])\\\
\quad\quad\quad\quad\quad\quad\quad\text{if}\ i\in\\{1,\dots,M_{\Omega}\\},\\\
B([\delta_{x_{i}}\circ\boldsymbol{L}^{\partial\Omega,u},u],[\delta_{x_{i}}\circ\boldsymbol{L}^{\partial\Omega,m},m],[\delta_{x_{i}}\circ\boldsymbol{L}^{\partial\Omega,\boldsymbol{V}},\boldsymbol{V}])\\\
\quad\quad\quad\quad\quad\quad\quad\text{if}\
i\in\\{M_{\Omega}+1,\dots,M\\}.\end{cases}$

Then, we rewrite (3.2) as

$\displaystyle\begin{cases}\min\limits_{(u,m,\overline{H},\boldsymbol{V})\in\mathcal{U}\times\mathcal{M}\times\mathbb{R}\times\boldsymbol{\mathcal{V}}}\|u\|_{\mathcal{U}}^{2}+\|m\|_{\mathcal{M}}^{2}+|\overline{H}|^{2}+\|\boldsymbol{V}\|_{\boldsymbol{\mathcal{V}}}^{2}+\frac{1}{\gamma^{2}}|[\boldsymbol{\phi}^{o},m]-\boldsymbol{m}^{o}|^{2}+|\Sigma^{-1}([\Psi,\boldsymbol{V}]-\boldsymbol{V}^{o})|^{2}\\\
\quad\quad\quad\quad\quad\quad\quad+\beta\bigg{|}\frac{1}{M_{\Omega}}\sum\limits_{i=1}^{M_{\Omega}}u(x_{i})\bigg{|}^{2}+\beta\bigg{|}\frac{1}{M_{\Omega}}\sum\limits_{i=1}^{M_{\Omega}}m(x_{i})-1\bigg{|}^{2}\\\
\text{s.t.}\quad
G([\widetilde{\boldsymbol{\phi}}^{u},u],[\widetilde{\boldsymbol{\phi}}^{m},m],\overline{H},[\widetilde{\boldsymbol{\phi}}^{\boldsymbol{V}},\boldsymbol{V}])=0.\end{cases}$
(3.3)

To solve (3.3), we introduce latent variables,
$\boldsymbol{z},\boldsymbol{\rho}$, and $\boldsymbol{v}$, and rewrite (3.3) as

$\displaystyle\begin{cases}\min\limits_{\boldsymbol{z},\boldsymbol{\rho},\boldsymbol{v},\overline{H}}\begin{cases}\min\limits_{(u,m,\boldsymbol{V})\in\mathcal{U}\times\mathcal{M}\times\boldsymbol{\mathcal{V}}}\|u\|_{\mathcal{U}}^{2}+\|m\|_{\mathcal{M}}^{2}+\|\boldsymbol{V}\|_{\boldsymbol{\mathcal{V}}}^{2}\\\
\text{s.t.}\quad[\boldsymbol{\delta}^{\Omega},u]=\boldsymbol{z}^{(1),\Omega},[\boldsymbol{\delta}^{\partial\Omega},u]=\boldsymbol{z}^{(1),\partial\Omega},[\widetilde{\boldsymbol{\phi}}^{u},u]=\boldsymbol{z}^{(2)},\\\
\quad\quad[\boldsymbol{\delta}^{\Omega},m]=\boldsymbol{\rho}^{(1),\Omega},[\boldsymbol{\delta}^{\partial\Omega},m]=\boldsymbol{\rho}^{(1),\partial\Omega},[\boldsymbol{\phi}^{o},m]=\boldsymbol{\rho}^{(2)},[\widetilde{\boldsymbol{\phi}}^{m},m]=\boldsymbol{\rho}^{(3)},\\\
\quad\quad[\widetilde{\boldsymbol{\phi}}^{\boldsymbol{V}},\boldsymbol{V}]=\boldsymbol{v}^{(1)},[\Psi,\boldsymbol{V}]=\boldsymbol{v}^{(2)},\end{cases}\\\
\quad\quad+|\overline{H}|^{2}+\frac{1}{\gamma^{2}}|\boldsymbol{\rho}^{(2)}-\boldsymbol{m}^{o}|^{2}+|\Sigma^{-1}(\boldsymbol{v}^{(2)}-\boldsymbol{V}^{o})|^{2}+\beta\bigg{|}\frac{1}{M_{\Omega}}\sum\limits_{i=1}^{M_{\Omega}}z^{(1),\Omega}_{i}\bigg{|}^{2}+\beta\bigg{|}\frac{1}{M_{\Omega}}\sum\limits_{i=1}^{M_{\Omega}}\rho^{(1),\Omega}_{i}-1\bigg{|}^{2}\\\
\text{s.t.}\quad
G(\boldsymbol{z}^{(2)},\boldsymbol{\rho}^{(3)},\overline{H},\boldsymbol{v}^{(1)})=0,\end{cases}$
(3.4)

where
$\boldsymbol{z}=(\boldsymbol{z}^{(1),\Omega},\boldsymbol{z}^{(1),\partial\Omega},\boldsymbol{z}^{(2)})$,
$\boldsymbol{\rho}=(\boldsymbol{\rho}^{(1),\Omega},\boldsymbol{\rho}^{(1),\partial\Omega},\boldsymbol{\rho}^{(2)},\boldsymbol{\rho}^{(3)})$,
and $\boldsymbol{v}=(\boldsymbol{v}^{(1)},\boldsymbol{v}^{(2)})$. Denote
$\boldsymbol{\phi}^{u}\mathrel{\mathop{\mathchar
58\relax}}=(\boldsymbol{\delta},\widetilde{\boldsymbol{\phi}}^{u})$,
$\boldsymbol{\phi}^{m}\mathrel{\mathop{\mathchar
58\relax}}=(\boldsymbol{\delta},\boldsymbol{\phi}^{o},\widetilde{\boldsymbol{\phi}}^{m})$,
and $\boldsymbol{\phi}^{\boldsymbol{V}}\mathrel{\mathop{\mathchar
58\relax}}=(\widetilde{\boldsymbol{\phi}}^{\boldsymbol{V}},\Psi)$.

For a kernel $K$ and a collection $\boldsymbol{\phi}$ of linear operators in
the dual space of the RKHS associated with $K$, let $K(x,\boldsymbol{\phi})$
be the vector with entries
$\int_{\Omega}K(x,x^{\prime})\phi_{i}(x^{\prime})\operatorname{d\\!}x^{\prime}$
and let $K(\boldsymbol{\phi},\boldsymbol{\phi})$, called the Gram matrix, be
with entries
$\int_{\Omega}\int_{\Omega}K(x,x)\phi_{i}(x^{\prime})\phi_{j}(x^{\prime})\operatorname{d\\!}x\operatorname{d\\!}x^{\prime}$.
Let $K_{u}$, $K_{m}$, and $K_{\boldsymbol{V}}$ be kernels associated with
RKHSs $\mathcal{U}$, $\mathcal{M}$, and $\mathcal{V}$. Let
$(u^{\dagger},m^{\dagger},\boldsymbol{V}^{\dagger})$ be the solution to the
first level minimization problem in (3.4) given
$(\boldsymbol{z},\boldsymbol{\rho},\boldsymbol{v},\overline{H})$. By the
representer formulas in (2.4) and (2.5), we get

$\displaystyle\begin{cases}u^{\dagger}(x)=K_{u}(x,\boldsymbol{\phi}^{u})K_{u}(\boldsymbol{\phi}^{u},\boldsymbol{\phi}^{u})^{-1}\boldsymbol{z},\\\
m^{\dagger}(x)=K_{m}(x,{\boldsymbol{\phi}}^{m})K({\boldsymbol{\phi}}^{m},{\boldsymbol{\phi}}^{m})^{-1}\boldsymbol{\rho},\\\
\boldsymbol{V}^{\dagger}(x)={K}_{\boldsymbol{V}}(x,{\boldsymbol{\phi}}^{\boldsymbol{V}}){K}_{\boldsymbol{V}}({\boldsymbol{\phi}}^{\boldsymbol{V}},{\boldsymbol{\phi}}^{\boldsymbol{V}})^{-1}\boldsymbol{v}.\end{cases}$
(3.5)

Thus, we have

$\displaystyle\begin{cases}\|u^{\dagger}\|_{\mathcal{U}}^{2}=\boldsymbol{z}^{T}K_{u}(\boldsymbol{\phi}^{u},\boldsymbol{\phi}^{u})^{-1}\boldsymbol{z},\\\
\|m^{\dagger}\|_{\mathcal{M}}^{2}=\boldsymbol{\rho}^{T}K_{m}({\boldsymbol{\phi}}^{m},{\boldsymbol{\phi}}^{m})^{-1}\boldsymbol{\rho},\\\
\|\boldsymbol{V}^{\dagger}\|_{\boldsymbol{\mathcal{V}}}^{2}=\boldsymbol{v}^{T}{K}_{\boldsymbol{V}}({\boldsymbol{\phi}}^{\boldsymbol{V}},{\boldsymbol{\phi}}^{\boldsymbol{V}})^{-1}\boldsymbol{v},\end{cases}$

Hence, we can formulate (3.4) as a finite-dimensional optimization problem

$\displaystyle\begin{cases}\min\limits_{\boldsymbol{z},\boldsymbol{\rho},\boldsymbol{v},\overline{H}}\boldsymbol{z}^{T}K_{u}(\boldsymbol{\phi}^{u},\boldsymbol{\phi}^{u})^{-1}\boldsymbol{z}+\boldsymbol{\rho}^{T}K_{m}({\boldsymbol{\phi}}^{m},{\boldsymbol{\phi}}^{m})^{-1}\boldsymbol{\rho}+\boldsymbol{v}^{T}{K}_{\boldsymbol{V}}({\boldsymbol{\phi}}^{\boldsymbol{V}},{\boldsymbol{\phi}}^{\boldsymbol{V}})^{-1}\boldsymbol{v}+|\overline{H}|^{2}\\\
\quad\quad\quad\quad+\frac{1}{\gamma^{2}}|\boldsymbol{\rho}^{(2)}-\boldsymbol{m}^{o}|^{2}+|\Sigma^{-1}(\boldsymbol{v}^{(2)}-\boldsymbol{V}^{o})|^{2}+\beta\bigg{|}\frac{1}{M_{\Omega}}\sum\limits_{i=1}^{M_{\Omega}}z^{(1),\Omega}_{i}\bigg{|}^{2}+\beta\bigg{|}\frac{1}{M_{\Omega}}\sum\limits_{i=1}^{M_{\Omega}}\rho^{(1),\Omega}_{i}-1\bigg{|}^{2}\\\
\text{s.t.}\quad
G(\boldsymbol{z}^{(2)},\boldsymbol{\rho}^{(3)},\overline{H},\boldsymbol{v}^{(1)})=0.\end{cases}$
(3.6)

To deal with the nonlinear constraints in (3.6), we introduce a prescribed
penalization parameter $\alpha>0$ and consider the following relaxation

$\displaystyle\begin{cases}\min\limits_{\boldsymbol{z},\boldsymbol{\rho},\boldsymbol{v},\overline{H}}\boldsymbol{z}^{T}K_{u}(\boldsymbol{\phi}^{u},\boldsymbol{\phi}^{u})^{-1}\boldsymbol{z}+\boldsymbol{\rho}^{T}K_{m}({\boldsymbol{\phi}}^{m},{\boldsymbol{\phi}}^{m})^{-1}\boldsymbol{\rho}+\boldsymbol{v}^{T}{K}_{\boldsymbol{V}}({\boldsymbol{\phi}}^{\boldsymbol{V}},{\boldsymbol{\phi}}^{\boldsymbol{V}})^{-1}\boldsymbol{v}+|\overline{H}|^{2}\\\
\quad\quad\quad\quad+\frac{1}{\gamma^{2}}|\boldsymbol{\rho}^{(2)}-\boldsymbol{m}^{o}|^{2}+|\Sigma^{-1}(\boldsymbol{v}^{(2)}-\boldsymbol{V}^{o})|^{2}+\beta\bigg{|}\frac{1}{M_{\Omega}}\sum\limits_{i=1}^{M_{\Omega}}z^{(1),\Omega}_{i}\bigg{|}^{2}+\beta\bigg{|}\frac{1}{M_{\Omega}}\sum\limits_{i=1}^{M_{\Omega}}\rho^{(1),\Omega}_{i}-1\bigg{|}^{2}\\\
\quad\quad\quad\quad+\alpha|G(\boldsymbol{z}^{(2)},\boldsymbol{\rho}^{(3)},\overline{H},\boldsymbol{v}^{(1)})|^{2}.\end{cases}$
(3.7)

The problem (3.7) is the foundation of the GP method for solving (3.2). We use
the Gauss–Newton method to solve (3.7) (see Section 3 of [13]).

###### Remark 3.1.

Generally, for a kernel $K$ and a vector $\boldsymbol{\phi}$ of linear
operators in the dual space of RKHS associated with $K$, the Gram matrix
$K(\boldsymbol{\phi},\boldsymbol{\phi})$ is ill-conditioned. To compute
$K(\boldsymbol{\phi},\boldsymbol{\phi})^{-1}$, we perform the Cholesky
decomposition on $K(\boldsymbol{\phi},\boldsymbol{\phi})+\eta R$, where
$\eta>0$ is a chosen regularization constant, and $R$ is a block diagonal
nugget constructed using the approach introduced in [13]. In the numerical
experiments, we precompute the Cholesky decomposition of
$K(\boldsymbol{\phi},\boldsymbol{\phi})+\eta R$ and store Cholesky factors for
further uses.

###  3.3. A Concrete Example

To illustrate the general framework described in Subsection 3.2, we apply our
method to solve a one-dimensional stationary MFG, which is analogous to the
examples in Section 5. Let $\mathbb{T}$ be the one-dimensional torus and be
characterized by $[0,1)$. Let $V^{*}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{T}\mapsto\mathbb{R}$ be a smooth solution. We consider the
following stationary MFG

$\displaystyle\begin{cases}\frac{(u_{x}^{*})^{2}}{2}+V^{*}(x)=(m^{*})^{2}+\overline{H}^{*},&\text{on}\
\mathbb{T},\\\ -(m^{*}(u_{x}^{*}))_{x}=0,&\text{on}\ \mathbb{T},\\\
\int_{\mathbb{T}}m^{*}\operatorname{d\\!}x=1,\int_{\mathbb{T}}u^{*}\operatorname{d\\!}x=0,&\end{cases}$
(3.8)

where $u^{*}$ is the value function, $m^{*}$ presents the probability density
of the population, and $\overline{H}^{*}$ is a real number. We note that given
$V^{*}$, (3.8) admits a smooth solution $(u^{*},m^{*},\overline{H}^{*})$ such
that $u^{*}=0$, $m^{*}=\sqrt{V^{*}-\overline{H}^{*}}$, and $\overline{H}^{*}$
is a constant satisfying
$\int_{\mathbb{T}}\sqrt{V^{*}-\overline{H}^{*}}\operatorname{d\\!}x=1$. In
this subsection, our goal is to recover $u^{*},m^{*},V^{*}$, and
$\overline{H}^{*}$ using the framework presented in the previous section under
noisy observations of $m^{*}$ and $V^{*}$, . More precisely, let
$\\{x_{i}^{o}\\}_{i=1}^{N_{m}}$ be a collection of observation points for
$m^{*}$ and let $\\{x_{i}^{o,1}\\}_{i=1}^{N_{v}^{1}}$ be another set of
observation points for $V^{*}$. Let $\boldsymbol{m}^{o}$ be a noisy
observation of $(m^{*}(x_{i}^{o}))_{i=1}^{N_{m}}$ such that
${m}^{o}_{i}-m^{*}(x_{i}^{o})\sim\mathcal{N}(0,\gamma^{2})$. Meanwhile, let
$\boldsymbol{V}^{o}$ be a noisy observation of $V^{*}$ at
$\\{x_{i}^{o,1}\\}_{i=1}^{N_{v}^{1}}$ such that
$V^{*}(\boldsymbol{x}^{o,1})-\boldsymbol{V}^{o}\sim\mathcal{N}(0,\Sigma)$,
where $\Sigma$ is the covariance matrix for noises. Let
$\\{x_{i}\\}_{i=1}^{M}$ be a collection of collocation points on $\mathbb{T}$.
Let $\mathcal{U}$, $\mathcal{M}$, and $\mathcal{V}$ be the RKHSs associated
with positive definite kernels $K_{u}$, $K_{m}$, and $K_{{V}}$, respectively.
Then, we recover $(u^{*},m^{*},V^{*},\overline{H}^{*})$ by a minimizer of the
following optimal recovery problem

$\displaystyle\begin{cases}\min\limits_{(u,m,V,\overline{H})\in\mathcal{U}\times\mathcal{M}\times\mathcal{V}\times\mathbb{R}}\|u\|_{\mathcal{U}}^{2}+\|m\|_{\mathcal{M}}^{2}+\|V\|_{\mathcal{V}}^{2}+|\overline{H}|^{2}+\frac{1}{\gamma^{2}}|[\boldsymbol{\delta}^{o},m]-\boldsymbol{m}^{o}|^{2}+|\Sigma^{-1}(\boldsymbol{V}(\boldsymbol{x}^{o,1})-\boldsymbol{V}^{o})|^{2}\\\
\quad\quad\quad\quad\quad\quad+\beta\bigg{|}\frac{1}{M}\sum\limits_{i=1}^{M}m(x_{i})-1\bigg{|}^{2}+\beta\bigg{|}\frac{1}{M}\sum\limits_{i=1}^{M}u(x_{i})\bigg{|}^{2}\\\
\text{s.t.}\
\frac{u_{x}^{2}(x_{i})}{2}+V(x_{i})=m^{2}(x_{i})+\overline{H},\quad\quad\quad\forall
i=1,\dots,M,\\\ \quad\ \
m_{x}(x_{i})u_{x}(x_{i})+m(x_{i})u_{xx}(x_{i})=0,\quad\forall
i=1,\dots,M,\end{cases}$ (3.9)

where $\beta>0$ is a penalization parameter. Let $K_{u}$, $K_{m}$, and
$K_{{V}}$ be kernels associated with RKHSs $\mathcal{U}$, $\mathcal{M}$, and
${\mathcal{V}}$ such that $u^{*}\in\mathcal{U}$, $m^{*}\in\mathcal{M}$, and
${V}^{*}\in{\mathcal{V}}$. Then, the constraints of (3.9) are well-defined.
Let $(u^{\dagger},m^{\dagger},\overline{H}^{\dagger})$ be a minimizer to
(3.9), whose existence is given by Theorem 3.2 later. Then, (3.9) can be
viewed as MAP estimators of GPs conditioned on the MFG at the sample points
$\\{x_{i}\\}_{i=1}^{M}$ while matching the observations. The conditioned GPs
are not Gaussian since the constraints in (3.9) are nonlinear.

Next, we rewrite (3.9) as a two-level optimization problem

$\displaystyle\begin{cases}\min\limits_{\boldsymbol{z},\boldsymbol{\rho},\boldsymbol{v},\overline{H}}\begin{cases}\min\limits_{(u,m,V)\in\mathcal{U}\times\mathcal{M}\times\mathcal{V}}\|u\|_{\mathcal{U}}^{2}+\|m\|_{\mathcal{M}}^{2}+\|{V}\|_{\mathcal{V}}^{2}\\\
\text{s.t.}\
u(x_{i})=z_{i}^{(1)},u_{x}(x_{i})=z_{i}^{(2)},u_{xx}(x_{i})=z_{i}^{(3)},\quad\quad\forall
i=1,\dots,M,\\\ \quad\ \
m(x_{i})=\rho_{i}^{(1)},m_{x}(x_{i})=\rho_{i}^{(2)},\quad\quad\quad\quad\quad\quad\quad\quad\forall
i=1,\dots,M,\\\ \quad\ \
m(x_{i}^{o})=\rho_{i}^{(3)},\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\
\forall i=1,\dots,N_{m},\\\ \quad\ \
V(x_{i})=v_{i}^{(1)},\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\
\ \forall i=1,\dots,M,\\\ \quad\ \
V(x_{i}^{o,1})=v_{i}^{(2)},\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\forall
i=1,\dots,N_{v}^{1},\end{cases}\\\
\quad\quad+|\overline{H}|^{2}+\frac{1}{\gamma^{2}}|\boldsymbol{\rho}^{(3)}-\boldsymbol{m}^{o}|^{2}+|\Sigma^{-1}(\boldsymbol{v}-\boldsymbol{V}^{o})|^{2}+\beta\bigg{|}\frac{1}{M}\sum\limits_{i=1}^{M}\rho_{i}^{(1)}-1\bigg{|}^{2}+\beta\bigg{|}\frac{1}{M}\sum\limits_{i=1}^{M}z_{i}^{(1)}\bigg{|}^{2}\\\
\text{s.t.}\
\frac{(z^{(2)}_{i})^{2}}{2}+v_{i}^{(1)}=(\rho_{i}^{(1)})^{2}+\overline{H},\quad\forall
i=1,\dots,M,\\\ \quad\ \
\rho_{i}^{(2)}z_{i}^{(2)}+\rho_{i}^{(1)}z_{i}^{(3)}=0,\quad\quad\quad\,\,\forall
i=1,\dots,M,\end{cases}$ (3.10)

where
$\boldsymbol{z}=(z^{(1)}_{1},\dots,z^{(1)}_{M},z^{(2)}_{1},\dots,z^{(2)}_{M},z^{(3)}_{1},\dots,z^{(3)}_{M})$,
$\boldsymbol{\rho}=(\rho^{(1)}_{1},\dots,\rho^{(1)}_{M},\rho^{(2)}_{1},\dots,\rho^{(2)}_{M},\rho_{1}^{(3)},\dots,\rho_{N_{m}}^{(3)})$,
and
$\boldsymbol{v}=(v^{(1)}_{1},\dots,v^{(1)}_{M},v_{1}^{(2)},\dots,v_{N_{v}^{1}}^{(2)})$.
Let $\delta_{x}$ be the Dirac delta function concentrated at $x$. We define
$\boldsymbol{\delta}=(\delta_{x_{i}})_{i=1}^{M}$,
$\boldsymbol{\delta}^{o}=(\delta_{x_{i}^{o}})_{i=1}^{N_{m}}$, and
$\boldsymbol{\delta}^{o,1}=(\delta_{x_{i}^{o,1}})_{i=1}^{N_{v}^{1}}$. Let
$\boldsymbol{\phi}^{u}=(\boldsymbol{\delta},\boldsymbol{\delta}\circ\partial_{x},\boldsymbol{\delta}\circ\partial_{xx})$,
$\boldsymbol{\phi}^{m}=(\boldsymbol{\delta},\boldsymbol{\delta}\circ\partial_{x},\boldsymbol{\delta}^{o})$,
and $\boldsymbol{\phi}^{V}=(\boldsymbol{\delta},\boldsymbol{\delta}^{o,1})$.
Let $(u^{\dagger},m^{\dagger},V^{\dagger})$ be the minimizer of the first
level optimization problem in (3.10). Then, by the representer formulas in
(2.4) (see [49, Sec. 17.8]), we get

$\displaystyle\begin{cases}u^{\dagger}(x)=K_{u}(x,\boldsymbol{\phi}^{u})K_{u}(\boldsymbol{\phi}^{u},\boldsymbol{\phi}^{u})^{-1}\boldsymbol{z},\\\
m^{\dagger}(x)=K_{m}(x,\boldsymbol{\phi}^{m})K_{m}(\boldsymbol{\phi}^{m},\boldsymbol{\phi}^{m})^{-1}\boldsymbol{\rho},\\\
V^{\dagger}(x)=K_{V}(x,\boldsymbol{\phi}^{V})K_{V}(\boldsymbol{\phi}^{V},\boldsymbol{\phi}^{V})^{-1}\boldsymbol{v}.\end{cases}$

Thus, we have

$\displaystyle\begin{cases}\|u^{\dagger}\|_{\mathcal{U}}^{2}=\boldsymbol{z}^{T}K_{u}(\boldsymbol{\phi}^{u},\boldsymbol{\phi}^{u})^{-1}\boldsymbol{z},\\\
\|m^{\dagger}\|_{\mathcal{M}}^{2}=\boldsymbol{\rho}^{T}K_{m}(\boldsymbol{\phi}^{m},\boldsymbol{\phi}^{m})^{-1}\boldsymbol{\rho},\\\
\|V^{\dagger}\|_{\mathcal{V}}^{2}=\boldsymbol{v}^{T}K_{V}(\boldsymbol{\phi}^{V},\boldsymbol{\phi}^{V})^{-1}\boldsymbol{v}.\end{cases}$

Hence, we can formulate (3.10) as a finite-dimensional optimization problem

$\displaystyle\begin{cases}\min\limits_{\boldsymbol{z},\boldsymbol{\rho},\boldsymbol{v},\overline{H}}\boldsymbol{z}^{T}K_{u}(\boldsymbol{\phi}^{u},\boldsymbol{\phi}^{u})^{-1}\boldsymbol{z}+\boldsymbol{\rho}^{T}K_{m}(\boldsymbol{\phi}^{m},\boldsymbol{\phi}^{m})^{-1}\boldsymbol{\rho}+\boldsymbol{v}^{T}K_{V}(\boldsymbol{\phi}^{V},\boldsymbol{\phi}^{V})^{-1}\boldsymbol{v}\\\
\quad\quad+|\overline{H}|^{2}+\frac{1}{\gamma^{2}}|\boldsymbol{\rho}^{(3)}-\boldsymbol{m}^{o}|^{2}+|\Sigma^{-1}(\boldsymbol{v}-\boldsymbol{V}^{o})|^{2}+\beta\bigg{|}\frac{1}{M}\sum\limits_{i=1}^{M}\rho_{i}^{(1)}-1\bigg{|}^{2}+\beta\bigg{|}\frac{1}{M}\sum\limits_{i=1}^{M}z_{i}^{(1)}\bigg{|}^{2}\\\
\text{s.t.}\
\frac{(z^{(2)}_{i})^{2}}{2}+v_{i}^{(1)}=(\rho_{i}^{(1)})^{2}+\overline{H},\forall
i=1,\dots,M,\\\ \quad\ \
\rho_{i}^{(2)}z_{i}^{(2)}+\rho_{i}^{(1)}z_{i}^{(3)}=0,\forall
i=1,\dots,M,\end{cases}$ (3.11)

To deal with the nonlinear constraints in (3.11), we introduce a prescribed
penalization parameter $\alpha>0$ and consider the following relaxation

$\displaystyle\begin{cases}\min\limits_{\boldsymbol{z},\boldsymbol{\rho},\boldsymbol{v},\overline{H}}\boldsymbol{z}^{T}K_{u}(\boldsymbol{\phi}^{u},\boldsymbol{\phi}^{u})^{-1}\boldsymbol{z}+\boldsymbol{\rho}^{T}K_{m}(\boldsymbol{\phi}^{m},\boldsymbol{\phi}^{m})^{-1}\boldsymbol{\rho}+\boldsymbol{v}^{T}K_{V}(\boldsymbol{\phi}^{V},\boldsymbol{\phi}^{V})^{-1}\boldsymbol{v}\\\
\quad\quad+|\overline{H}|^{2}+\frac{1}{\gamma^{2}}|\boldsymbol{\rho}^{(3)}-\boldsymbol{m}^{o}|^{2}+|\Sigma^{-1}(\boldsymbol{v}-\boldsymbol{V}^{o})|^{2}+\beta\bigg{|}\frac{1}{M}\sum\limits_{i=1}^{M}\rho_{i}^{(1)}-1\bigg{|}^{2}+\beta\bigg{|}\frac{1}{M}\sum\limits_{i=1}^{M}z_{i}^{(1)}\bigg{|}^{2}\\\
\quad\quad+\alpha\sum_{i=1}^{M}|\frac{(z^{(2)}_{i})^{2}}{2}+v_{i}^{(1)}-(\rho_{i}^{(1)})^{2}-\overline{H}|^{2}+\alpha\sum_{i=1}^{M}|\rho_{i}^{(2)}z_{i}^{(2)}+\rho_{i}^{(1)}z_{i}^{(3)}|^{2}.\end{cases}$
(3.12)

In numerical implementations, we use the Gauss–Newton method to solve (3.12)
(see also Section 3 of [13]).

###  3.4. Existence of Minimizers

By adapting the proof of Theorem 3.3 in [47], the next theorem shows that
there exists a solution to (3.2). Hence, we establish the solvability of the
optimization problem in (3.2).

###### Theorem 3.2.

Let $\mathcal{U}$, $\mathcal{M}$, and $\boldsymbol{\mathcal{V}}$ be RKHSs
associated with kernels $K_{u}$, $K_{m}$, and $K_{\boldsymbol{V}}$. Assume
that $(u^{*},m^{*},\overline{H}^{*},\boldsymbol{V}^{*})$ satisfies the MFG in
(3.1) such that $u^{*}\in\mathcal{U}$, $m^{*}\in\mathcal{M}$,
$\overline{H}^{*}\in\mathbb{R}$, and
$\boldsymbol{V}^{*}\in\boldsymbol{\mathcal{V}}$. Define the vectors of linear
operators $\boldsymbol{\phi}^{u}$, $\boldsymbol{\phi}^{m}$, and
$\boldsymbol{\phi}^{\boldsymbol{V}}$ as outlined in Subsection 3.2. Under
Assumption 1, the minimization problem in (3.2) admits a minimizer
$(u^{\dagger},m^{\dagger},\overline{H}^{\dagger},\boldsymbol{V}^{\dagger})$
such that

$\displaystyle\begin{cases}u^{\dagger}(x)=K_{u}(x,\boldsymbol{\phi}^{u})K_{u}(\boldsymbol{\phi}^{u},\boldsymbol{\phi}^{u})^{-1}\boldsymbol{z}^{\dagger},\\\
m^{\dagger}(x)=K_{m}(x,{\boldsymbol{\phi}}^{m})K({\boldsymbol{\phi}}^{m},{\boldsymbol{\phi}}^{m})^{-1}\boldsymbol{\rho}^{\dagger},\\\
\boldsymbol{V}^{\dagger}(x)={K}_{\boldsymbol{V}}(x,{\boldsymbol{\phi}}^{\boldsymbol{V}}){K}_{\boldsymbol{V}}({\boldsymbol{\phi}}^{\boldsymbol{V}},{\boldsymbol{\phi}}^{\boldsymbol{V}})^{-1}\boldsymbol{v}^{\dagger}.\end{cases}$
(3.13)

where
$(\boldsymbol{z}^{\dagger},\boldsymbol{\rho}^{\dagger},\overline{H}^{\dagger},\boldsymbol{V}^{\dagger})$
is a minimizer of (3.6).

###### Proof.

By the above arguments, (3.2) is equivalent to (3.6). Hence, the key is to
prove the existence of a minimizer to (3.6). The argument here is similar to
the proof of Theorem 3.3 in [47]. Let
$(u^{*},m^{*},\overline{H}^{*},\boldsymbol{V}^{*})$ be the solution to (3.1).
We define the tuple
$(\boldsymbol{z}_{*},\boldsymbol{\rho}_{*},\boldsymbol{v}_{*})$ such that
$\boldsymbol{z}_{*}=[\boldsymbol{\phi}^{u},u]$,
$\boldsymbol{\rho}_{*}=[\boldsymbol{\phi}^{m},m]$, and
$\boldsymbol{v}_{*}=[\boldsymbol{\phi}^{\boldsymbol{V}},\boldsymbol{V}]$.
Then,
$(\boldsymbol{z}_{*},\boldsymbol{\rho}_{*},\overline{H}^{*},\boldsymbol{v}_{*})$
satisfies the constraints in (3.6). For a little abuse of notations, we denote
by $G(\boldsymbol{z},\boldsymbol{\rho},\ \overline{H},\boldsymbol{v})=0$ the
constraints in (3.6) and define

$\displaystyle\begin{split}\mathcal{C}_{1}=&\bigg{\\{}(\boldsymbol{z},\boldsymbol{\rho},\overline{H},\boldsymbol{v})|\boldsymbol{z}^{T}K_{u}(\boldsymbol{\phi}^{u},\boldsymbol{\phi}^{u})^{-1}\boldsymbol{z}+\boldsymbol{\rho}^{T}K_{m}({\boldsymbol{\phi}}^{m},{\boldsymbol{\phi}}^{m})^{-1}\boldsymbol{\rho}+\boldsymbol{v}^{T}{K}_{\boldsymbol{V}}({\boldsymbol{\phi}}^{\boldsymbol{V}},{\boldsymbol{\phi}}^{\boldsymbol{V}})^{-1}\boldsymbol{v}+|\overline{H}|^{2}\\\
&\quad\quad\quad\quad+\frac{1}{\gamma^{2}}|\boldsymbol{\rho}^{(2)}-\boldsymbol{m}^{o}|^{2}+|\Sigma^{-1}(\boldsymbol{v}^{(2)}-\boldsymbol{V}^{o})|^{2}+\beta\bigg{|}\frac{1}{M_{\Omega}}\sum\limits_{i=1}^{M_{\Omega}}z^{(1),\Omega}_{i}\bigg{|}^{2}+\beta\bigg{|}\frac{1}{M_{\Omega}}\sum\limits_{i=1}^{M_{\Omega}}\rho^{(1),\Omega}_{i}-1\bigg{|}^{2}\\\
&\quad\quad\quad\quad\leqslant\boldsymbol{z}^{T}_{*}K_{u}(\boldsymbol{\phi}^{u},\boldsymbol{\phi}^{u})^{-1}\boldsymbol{z}_{*}+\boldsymbol{\rho}^{T}_{*}K_{m}({\boldsymbol{\phi}}^{m},{\boldsymbol{\phi}}^{m})^{-1}\boldsymbol{\rho}_{*}+\boldsymbol{v}^{T}_{*}{K}_{\boldsymbol{V}}({\boldsymbol{\phi}}^{\boldsymbol{V}},{\boldsymbol{\phi}}^{\boldsymbol{V}})^{-1}\boldsymbol{v}_{*}+|\overline{H}^{*}|^{2}\\\
&\quad\quad\quad\quad+\frac{1}{\gamma^{2}}|\boldsymbol{\rho}^{(2)}_{*}-\boldsymbol{m}^{o}|^{2}+|\Sigma^{-1}(\boldsymbol{v}^{(2)}_{*}-\boldsymbol{V}^{o})|^{2}+\beta\bigg{|}\frac{1}{M_{\Omega}}\sum\limits_{i=1}^{M_{\Omega}}z^{(1),\Omega}_{*,i}\bigg{|}^{2}+\beta\bigg{|}\frac{1}{M_{\Omega}}\sum\limits_{i=1}^{M_{\Omega}}\rho^{(1),\Omega}_{*,i}-1\bigg{|}^{2}\bigg{\\}}.\end{split}$

Then, (3.6) is equivalent to

$\displaystyle\begin{cases}\min\limits_{\boldsymbol{z},\boldsymbol{\rho},\overline{H},\boldsymbol{v}}\boldsymbol{z}^{T}K_{u}(\boldsymbol{\phi}^{u},\boldsymbol{\phi}^{u})^{-1}\boldsymbol{z}+\boldsymbol{\rho}^{T}K_{m}({\boldsymbol{\phi}}^{m},{\boldsymbol{\phi}}^{m})^{-1}\boldsymbol{\rho}+\boldsymbol{v}^{T}{K}_{\boldsymbol{V}}({\boldsymbol{\phi}}^{\boldsymbol{V}},{\boldsymbol{\phi}}^{\boldsymbol{V}})^{-1}\boldsymbol{v}+|\overline{H}|^{2}\\\
\quad\quad\quad\quad+\frac{1}{\gamma^{2}}|\boldsymbol{\rho}^{(2)}-\boldsymbol{m}^{o}|^{2}+|\Sigma^{-1}(\boldsymbol{v}^{(2)}-\boldsymbol{V}^{o})|^{2}+\beta\bigg{|}\frac{1}{M_{\Omega}}\sum\limits_{i=1}^{M_{\Omega}}z^{(1),\Omega}_{i}\bigg{|}^{2}+\beta\bigg{|}\frac{1}{M_{\Omega}}\sum\limits_{i=1}^{M_{\Omega}}\rho^{(1),\Omega}_{i}-1\bigg{|}^{2}\\\
\text{s.t.}\
(\boldsymbol{z},\boldsymbol{\rho},\overline{H},\boldsymbol{v})\in\mathcal{C}\mathrel{\mathop{\mathchar
58\relax}}=\bigg{\\{}(\boldsymbol{z},\boldsymbol{\rho},\overline{H},\boldsymbol{v})|G(\boldsymbol{z},\boldsymbol{\rho},\overline{H},\boldsymbol{v})=0\bigg{\\}}\cap\mathcal{C}_{1}.\end{cases}$
(3.14)

By the continuity of $G$ in
$(\boldsymbol{z},\boldsymbol{\rho},\overline{H},\boldsymbol{v})$ and the fact
that
$(\boldsymbol{z}_{*},\boldsymbol{\rho}_{*},\overline{H}^{*},\boldsymbol{v}_{*})\in\mathcal{C}$,
$\mathcal{C}$ is compact and non-empty. Hence, the objective function (3.14)
achieves a minimum on $\mathcal{C}$. Thus, (3.6) admits a minimizer. We
conclude (3.13) by (3.5). ∎

###  3.5. A Probabilistic Perspective

Here, we give a probabilistic interpretation of the GP method for the MFG
inverse problems. First, we put independent Gaussian priors on the solution
$(u^{*},m^{*},\overline{H}^{*},\boldsymbol{V}^{*})$ to (3.1) such that
$u^{*}\sim\mathcal{N}(0,{K}_{u})$, $m^{*}\sim\mathcal{N}(0,{K}_{m})$,
$\overline{H}^{*}\sim\mathcal{N}(0,1)$, and
$V_{i}^{*}\sim\mathcal{N}(0,{K}_{V_{i}})$, where $K_{u}$, $K_{m}$, $K_{V_{i}}$
are kernels associated with the RKHSs $\mathcal{U}$, $\mathcal{M}$, and
$\mathcal{V}_{i}$, $i=1,\dots,N_{v}$. Let $\boldsymbol{\phi}^{u}$,
$\boldsymbol{\phi}^{m}$, and $\boldsymbol{\phi}^{\boldsymbol{V}}$ be the
vectors of linear operators defined in Subsection 3.2. Let $\boldsymbol{z}$,
$\boldsymbol{\rho}$, and $\boldsymbol{v}$ be observations defined as

$\displaystyle\boldsymbol{z}=[\boldsymbol{\phi}^{u},u^{*}],\boldsymbol{\rho}=[\boldsymbol{\phi}^{m},m^{*}],\text{
and }\boldsymbol{v}=[\boldsymbol{\phi}^{\boldsymbol{V}},\boldsymbol{V}^{*}],$

Thus,
$\boldsymbol{z}\sim\mathcal{N}(0,K_{u}(\boldsymbol{\phi}^{u},\boldsymbol{\phi}^{u}))$,
$\boldsymbol{\rho}\sim\mathcal{N}(0,K_{m}(\boldsymbol{\phi}^{m},\boldsymbol{\phi}^{m}))$
,
$\boldsymbol{v}\sim\mathcal{N}(0,K_{\boldsymbol{V}}(\boldsymbol{\phi}^{\boldsymbol{V}},\boldsymbol{\phi}^{\boldsymbol{V}}))$.
Moreover, the probability densities of $\boldsymbol{z}$, $\boldsymbol{\rho}$,
and $\boldsymbol{V}$ satisfy

$\displaystyle
p(\boldsymbol{z})=\frac{1}{C_{\boldsymbol{z}}\sqrt{\operatorname{det}(K(\boldsymbol{\phi}^{u},\boldsymbol{\phi}^{u}))}}e^{-\frac{1}{2}\boldsymbol{z}^{T}K(\boldsymbol{\phi}^{u},\boldsymbol{\phi}^{u})^{-1}\boldsymbol{z}},$
$\displaystyle
p(\boldsymbol{\rho})=\frac{1}{C_{\boldsymbol{\rho}}\sqrt{\operatorname{det}(K(\boldsymbol{\phi}^{m},\boldsymbol{\phi}^{m}))}}e^{-\frac{1}{2}\boldsymbol{\rho}^{T}K(\boldsymbol{\phi}^{m},\boldsymbol{\phi}^{m})^{-1}\boldsymbol{\rho}},$
$\displaystyle
p(\boldsymbol{v})=\frac{1}{C_{\boldsymbol{v}}\sqrt{\operatorname{det}(K(\boldsymbol{\phi}^{\boldsymbol{V}},\boldsymbol{\phi}^{\boldsymbol{V}}))}}e^{-\frac{1}{2}\boldsymbol{v}^{T}K(\boldsymbol{\phi}^{\boldsymbol{V}},\boldsymbol{\phi}^{\boldsymbol{V}})^{-1}\boldsymbol{v}},$

where $C_{\boldsymbol{z}}$, $C_{\boldsymbol{\rho}}$, and $C_{\boldsymbol{v}}$
are constants guaranteeing the unit mass of the probability density. By
assumptions, $\boldsymbol{m}^{o}$ is a noisy observation of
$\boldsymbol{\rho}^{(2)}$, i.e.
$\boldsymbol{\epsilon}=\boldsymbol{\rho}^{(2)}-\boldsymbol{m}^{0}$,
$\boldsymbol{\epsilon}\sim\mathcal{N}(0,\gamma^{2}I)$, where
$\boldsymbol{\epsilon}$ is independent with other random variables. Thus,
$\gamma^{-1}(\boldsymbol{\rho}^{(2)}-\boldsymbol{m}^{o})\sim\mathcal{N}(0,I)$
and $\gamma^{-1}(\boldsymbol{\rho}^{(2)}-\boldsymbol{m}^{o})$ is independent
with other random variables. Similarly,
$\Sigma^{-1}(\boldsymbol{v}^{(2)}-\boldsymbol{V}^{o})\sim\mathcal{N}(0,I)$ and
is another independent random vector. Then, solving Problem 2 by (3.6) is
equivalent to finding $\boldsymbol{z}$, $\boldsymbol{\rho}$,
$\overline{\boldsymbol{H}}$, and $\boldsymbol{v}$ with the maximum likelihood
while satisfying the PDE system and with penalizations on noisy observations
and the means of $\boldsymbol{z}$ and $\boldsymbol{\rho}$, i.e., (3.6) is
equivalent to

$\displaystyle\begin{cases}\min\limits_{\boldsymbol{z},\boldsymbol{\rho},\boldsymbol{v},\overline{H}}-\log
p(\boldsymbol{z})-\log p(\boldsymbol{\rho})-\log p(\overline{H})-\log
p(\boldsymbol{v})-\log
p(\gamma^{-1}(\boldsymbol{\rho}^{(2)}-\boldsymbol{m}^{o}))\\\
\quad\quad\quad\quad-\log
p(\Sigma^{-1}(\boldsymbol{v}^{(2)}-\boldsymbol{V}^{o}))+\beta\bigg{|}\frac{1}{M_{\Omega}}\sum\limits_{i=1}^{M_{\Omega}}z^{(1),\Omega}_{i}\bigg{|}^{2}+\beta\bigg{|}\frac{1}{M_{\Omega}}\sum\limits_{i=1}^{M_{\Omega}}\rho^{(1),\Omega}_{i}-1\bigg{|}^{2}\\\
\text{s.t.}\quad
G(\boldsymbol{z}^{(2)},\boldsymbol{\rho}^{(3)},\overline{H},\boldsymbol{v}^{(1)})=0.\end{cases}$

##  4\. An Inverse Problem For Time-dependent MFGs

This section considers Problem 1 in settings of time-dependent MFGs and
introduces a GP framework for solving the inverse problem. Differing from
stationary MFGs, time-dependent MFGs typically don’t incorporate mass
constraints, a feature naturally ensured by PDEs. In Subsection 4.1, we
reformulate Problem 1 for time-dependent MFGs. Subsection 4.2 details the GP
framework built upon a forward-backward time-discretization scheme. In
contrast to the approach outlined in [13], our methodology is distinctively
tailored for time-dependent MFGs involving forward-backward PDEs. In this
context, neither a purely forward nor a purely backward Euler scheme is
applicable. Instead, we integrate a mixed forward-backward time discretization
technique, which is proposed in [1], with GPs for the inverse problem.

Another feasible approach for solving the inverse problems of time-dependent
MFGs is to treat the temporal variable as an additional spatial variable,
applying the methods described in the preceding section to address
corresponding challenges. Here, we propose a time-discretized approach, where
the time-dependent PDE system is discretized in time and GPs are employed to
approximate the unknowns in the entire time-independent PDE system. This
discretization offers two main benefits: firstly, it separates the
complexities of time and space, simplifying kernel selection based on spatial
variability; secondly, it allows for distinct kernel choices at different time
slices, catering to the varying spatial regularities of functions over time.

###  4.1. An Inverse Problem for Time-dependent MFGs

Let $\Omega\subset\mathbb{R}^{d}$. In this section, we discuss the inverse
problem for time-dependent MFGs of the following form

$\displaystyle\begin{cases}-\partial_{t}u^{*}+\mathcal{P}_{1}(u^{*},m^{*},\boldsymbol{V}^{*})(t,x)=0,&\forall(t,x)\in(0,T)\times\operatorname{int}\Omega,\\\
\partial_{t}m^{*}+\mathcal{P}_{2}(u^{*},m^{*},\boldsymbol{V}^{*})(t,x)=0,&\forall(t,x)\in(0,T)\times\operatorname{int}\Omega,\\\
\mathcal{B}(u^{*},m^{*},\boldsymbol{V}^{*})(t,x)=0,&\forall(t,x)\in(0,T)\times\partial\Omega,\\\
u^{*}(T,x)=\phi^{*}(x),m^{*}(0,x)=\mu^{*}(x),&\forall
x\in\overline{\Omega}.\end{cases}$ (4.1)

where $\boldsymbol{V}^{*}$ comprises a set of functions that depict
environment configurations, $\phi^{*}$ represents the terminal cost of agents,
and $\mu^{*}$ is the initial distribution of agents. Given
$\boldsymbol{V}^{*}$, $\mu^{*}$, and $\phi^{*}$, $(u^{*},m^{*})$ solves the
time-dependent MFG such that $u^{*}$ stands for the value function and $m^{*}$
represents the density of the population of agents. In this context,
$\mathcal{P}_{1}$ and $\mathcal{P}_{2}$ are two nonlinear operators, while
$\mathcal{B}$ stands for a boundary operator. We proceed under the assumption
that equation (4.1) admits a unique classical solution $(u^{*},m^{*})$ when
the underlying environment configuration $\boldsymbol{V}^{*}$, the initial
distribution $\mu^{*}$, and the terminal cost $\phi^{*}$ are given. We detail
Problem 1 in the setting of time-dependent MFGs of the form (4.1) in the
following problem.

###### Problem 3.

Consider a set of functions on $\Omega$, denoted by
$\boldsymbol{V}^{*}=\\{V^{*}_{i}\\}_{i=1}^{N_{v}}$, with $N_{v}\in\mathbb{N}$.
Let ${\phi}^{*}$ be a function defined over $\overline{\Omega}$, and let
$\mu^{*}$ denote a distribution on $\Omega$. When provided with
$\boldsymbol{V}^{*}$, $\phi^{*}$, and ${\mu}^{*}$, (4.1) yields a unique
classical solution $(u^{*},m^{*})$. Our objective is to estimate the
quantities $u^{*}$, $m^{*}$, $\boldsymbol{V}^{*}$, $\mu^{*}$, and $\phi^{*}$
based on partial noisy observations of $m^{*}$ and of a subset
$\boldsymbol{V}_{o}^{*}$ from $\boldsymbol{V}^{*}$. Specifically, we have

  1. 1.

Partial noisy observations of $m^{*}$. We possess a set of linear operators
$\\{\phi^{o}_{l}\\}_{l=1}^{N_{m}}$, $N_{m}\in\mathbb{N}$ at which we have
imprecise observations of $m^{*}$. These observations are denoted by
$\boldsymbol{m}^{o}$. Formally,
$\boldsymbol{m}^{o}=([\phi_{1}^{o},m^{*}],\dots,[\phi_{N_{m}}^{o},m^{*}])+\boldsymbol{\epsilon}$
where $\boldsymbol{\epsilon}\sim\mathcal{N}(0,\gamma^{2}I)$.

  2. 2.

Partial noisy observations of $\boldsymbol{V}^{*}$. Denote a subset of
$\boldsymbol{V}^{*}$ by
$\boldsymbol{V}_{o}^{*}=\\{V_{j}^{*}\\}_{j\in\mathcal{O}_{v}}$, where
$\mathcal{O}_{v}\subset\\{1,\dots,N_{v}\\}$. For this subset, we have
observations at specific sets of linear operators,
$\\{\\{\psi_{e}^{o,j}\\}_{e=1}^{N^{j}_{v}}\\}_{j\in\mathcal{O}_{v}}$. The
observations are given by

$\displaystyle\boldsymbol{V^{o}}=(\boldsymbol{V}^{o}_{j})_{j\in\mathcal{O}_{v}},\text{
where
}\boldsymbol{V}^{o}_{j}=([\psi_{1}^{o,j},V_{j}^{*}],\dots,[\psi_{N_{v}}^{o,j},V_{j}^{*}]+\boldsymbol{\epsilon}_{j},\boldsymbol{\epsilon}_{j}\sim\mathcal{N}(0,\gamma_{j}^{2}I).$

Inverse Problem Statement

Suppose that agents play the time-dependent MFG specified in (4.1), we infer
$u^{*}$, $m^{*}$, $\boldsymbol{V}^{*}$, $\boldsymbol{\phi}^{*}$, and $\mu^{*}$
based on the data $\boldsymbol{m}^{o}$ and $\boldsymbol{V}^{o}$.

###  4.2. A Time-discretization Framework

In this section, we present a GP framework to solve Problem 3. This is
accomplished by utilizing a time-discretized variant of (4.1), which is given
in (4.2). The key challenge in (4.1) lies in its incorporation of coupled
forward-backward PDEs. To effectively handle this, we employ a mixed forward-
backward time-discretization approach. This methodology is inspired by and
builds upon established strategies found in the MFG literature [1, 3]. The key
difference between our discretization given in (4.2) and the schemes in [1, 3]
is that we only discretize (4.1) in time and approximate solutions in the
resulting time-independent PDEs by GPs but they discretize the MFG both in
time and space by using the adjoint structure inherent in MFGs. Therefore, the
existence and uniqueness of the solution to the discretized PDE system in
(4.2) can be proved using the same fixed-point argument as in [1].

Let $\\{t_{k}\\}_{k=0}^{N_{T}}$ be an equiv-distributed grid points on the
time interval $[0,T]$ with $t_{0}=0$ and $t_{N_{T}}=T$. We approximate $u$,
$m$, and $\boldsymbol{V}$ at the time grid point by sets of functions
$\\{u_{k}\\}_{k=0}^{N_{T}}$, $\\{m_{k}\\}_{k=0}^{N_{T}}$, and
$\\{\boldsymbol{V}_{k}\\}_{k=0}^{N_{T}}$ such that $u_{k}\approx u(t_{k})$,
$m_{k}\approx m(t_{k})$, and $\boldsymbol{V}_{k}\approx\boldsymbol{V}(t_{k})$
for $k=0,\dots,N_{T}$. Then, we discretize (4.1) in time and get the following
system

$\displaystyle\begin{cases}-\frac{u_{k+1}(x)-u_{k}(x)}{\Delta
t}+\mathcal{P}_{1}(u_{k},m_{k+1},\boldsymbol{V}_{k})(x)=0,&\forall
x\in\operatorname{int}\Omega,k\in\\{0,\dots,N_{T}-1\\},\\\
\frac{m_{k+1}(x)-m_{k}(x)}{\Delta
t}+\mathcal{P}_{2}(u_{k},m_{k+1},\boldsymbol{V}_{k})(x)=0,&\forall
x\in\operatorname{int}\Omega,k\in\\{0,\dots,N_{T}-1\\},\\\
\mathcal{B}(u_{k},m_{k+1},\boldsymbol{V}_{k})(x)=0,&\forall
x\in\partial\Omega,k\in\\{0,\dots,N_{T}-1\\},\\\
u_{N_{T}}(x)=\phi^{*}(x),m_{0}(x)=\mu^{*}(x),&\forall
x\in\overline{\Omega}.\end{cases}$ (4.2)

To simplify the presentation, we assume that the time grid points are rich
enough such that noisy observations of the unknown variables occur only at
discrete time slots. If observations are made between these intervals, the
issue can be addressed by either employing a more refined time discretization
or interpolating between the values of adjacent time slices. Next, we impose a
GP prior on each unknown and approximate unknowns by GPs conditioned on PDEs
while penalizing the observations. Let $\\{x_{i}\\}_{i=1}^{M}$ be a set of
samples such that $x_{1},\dots,x_{M_{\Omega}}\in\operatorname{int}\Omega$ and
$x_{M_{\Omega}+1},\dots,x_{M}\in\partial\Omega$ for $1\leqslant
M_{\Omega}\leqslant M$. Let
$(\\{u_{k}^{*}\\}_{k=0}^{N_{T}},\\{m_{k}^{*}\\}_{k=0}^{N_{T}},\\{\boldsymbol{V}_{k}^{*}\\}_{k=0}^{N_{T}})$
be the solution of (4.2). Denote by $\\{\mathcal{U}_{k}\\}_{k=0}^{N_{T}}$ and
$\\{\mathcal{M}_{k}\\}_{k=0}^{N_{T}}$ the collections of RKHSs such that
$u_{k}^{*}\in\mathcal{U}_{k}$ and $m_{k}^{*}\in\mathcal{M}_{k}$ for each
$0\leqslant k\leqslant N_{T}$. Meanwhile, let
$\boldsymbol{\mathcal{V}}\mathrel{\mathop{\mathchar
58\relax}}=\\{\\{\mathcal{V}_{i,k}\\}_{k=1}^{N_{T}}\\}_{i=1}^{N_{v}}$ be a
collection of RKHSs where $V_{i,k}^{*}\in\mathcal{V}_{i,k}$ and let
$\boldsymbol{\mathcal{V}}_{k}=\otimes_{{i=1}}^{N_{v}}\mathcal{V}_{i,k}$. Here,
we assume that the observation operators for $m^{*}$ are fixed for all the
time slots and collect them into a vector $\boldsymbol{\phi}^{o}$. Let
$\boldsymbol{m}^{o}_{k}$ be a noisy observation of
$[\boldsymbol{\phi}^{o},m_{k}^{*}]$ such that
$([\boldsymbol{\phi}^{o},m_{k}^{*}]-\boldsymbol{m}^{o}_{k})\sim\mathcal{N}(0,\gamma_{k}^{2}I)$.
Let $\Psi_{k}$ be a collection of linear functionals in
$\boldsymbol{\mathcal{V}}^{*}_{k}$, the dual of
$\boldsymbol{\mathcal{V}}_{k}$, and let $\boldsymbol{V}_{k}^{o}$ be a noisy
observation of $[\Psi_{k},\boldsymbol{V}^{*}_{k}]$ such that
$([\Psi_{k},\boldsymbol{V}^{*}_{k}]-\boldsymbol{V}^{o}_{k})\sim\mathcal{N}(0,\Sigma_{k})$,
where $\Sigma_{k}$ is the covariance matrix for the noise. Then, to solve
Problem 3, we approximate $(u_{k},m_{k},\boldsymbol{V}_{k})_{k=0}^{N_{T}}$ by
a minimizer of the following optimal recovery problem

$\displaystyle\begin{cases}\min\limits_{\overset{(u_{k},m_{k},\boldsymbol{V}_{k})\in\mathcal{U}_{k}\times\mathcal{M}_{k}\times\boldsymbol{\mathcal{V}}_{k}}{k=0,\dots,N_{T}}}\sum_{k=0}^{N_{T}}\|u_{k}\|_{\mathcal{U}_{k}}^{2}+\sum_{k=0}^{N_{T}}\|m_{k}\|_{\mathcal{M}_{k}}^{2}+\sum_{k=0}^{N_{T}}\|\boldsymbol{V}_{k}\|_{\boldsymbol{\mathcal{V}}_{k}}^{2}\\\
\quad\quad\quad\quad\quad\quad\quad\quad+\sum_{k=0}^{N_{T}}\frac{1}{\gamma^{2}_{k}}|[\boldsymbol{\phi}^{o},m_{k}]-\boldsymbol{m}^{o}_{k}|^{2}+\sum_{k=0}^{N_{T}}|\Sigma^{-1}_{k}([\Psi_{k},\boldsymbol{V}_{k}]-\boldsymbol{V}^{o}_{k})|^{2}\\\
\operatorname{s.t.}\quad-\frac{u_{k+1}(x_{i})-u_{k}(x_{i})}{\Delta
t}+\mathcal{P}_{1}(u_{k},m_{k+1},\boldsymbol{V}_{k})(x_{i})=0,\forall
i\in\\{1,\dots,M_{\Omega}\\},k\in\\{0,\dots,N_{T}-1\\},\\\
\quad\quad\quad\frac{m_{k+1}(x_{i})-m_{k}(x_{i})}{\Delta
t}+\mathcal{P}_{2}(u_{k},m_{k+1},\boldsymbol{V}_{k})(x_{i})=0,\forall
i\in\\{1,\dots,M_{\Omega}\\},k\in\\{0,\dots,N_{T}-1\\},\\\
\quad\quad\quad\mathcal{B}(u_{k},m_{k+1},\boldsymbol{V}_{k})(x_{i})=0,\,\,\,\quad\quad\quad\quad\quad\quad\quad\quad\forall
i\in\\{M_{\Omega}+1,\dots,M\\},k\in\\{0,\dots,N_{T}-1\\},\\\ \quad\quad\quad
u_{N_{T}}(x_{i})=\phi^{*}(x_{i}),m_{0}(x_{i})=\mu^{*}(x_{i}),\,\quad\quad\quad\quad\forall
i\in\\{1,\dots,M\\}.\end{cases}$ (4.3)

Similar arguments to those outlined in the proof of Theorem 3.2 can be applied
to establish the existence of minimizers in the context of (4.3).

Next, we reformulate (4.3) into a finite dimensional minimization problem. For
this transformation, we adopt specific assumptions regarding the structures of
$\mathcal{P}_{1}$, $\mathcal{P}_{2}$, and $\mathcal{B}$ as follows.

###### Assumption 2.

Assume that there exists collections of bounded linear operators
$\boldsymbol{L}^{u,\Omega}$, $\boldsymbol{L}^{u,\partial\Omega}$,
$\boldsymbol{L}^{m,\Omega}$, $\boldsymbol{L}^{m,\partial\Omega}$,
$\boldsymbol{L}^{\boldsymbol{{V}},\Omega}$, and
$\boldsymbol{L}^{\boldsymbol{{V}},\partial\Omega}$, and continuous nonlinear
maps $P_{1}$, $P_{2}$, and $B$ such that

$\displaystyle\begin{cases}\mathcal{P}_{1}(u,m,\boldsymbol{V})(x)=P_{1}(\boldsymbol{L}^{u,\Omega}(u)(x),\boldsymbol{L}^{m,\Omega}(m)(x),\boldsymbol{L}^{\boldsymbol{V},\Omega}(\boldsymbol{V})(x)),\forall
x\in\operatorname{int}\Omega,\\\
\mathcal{P}_{2}(u,m,\boldsymbol{V})(x)=P_{2}(\boldsymbol{L}^{u,\Omega}(u)(x),\boldsymbol{L}^{m,\Omega}(m)(x),\boldsymbol{L}^{\boldsymbol{V},\Omega}(\boldsymbol{V})(x)),\forall
x\in\operatorname{int}\Omega,\\\
\mathcal{B}(u,m,\boldsymbol{V})(x)=B(\boldsymbol{L}^{u,\partial\Omega}(u)(x),\boldsymbol{L}^{m,\partial\Omega}(m)(x),\boldsymbol{L}^{\boldsymbol{V},\partial\Omega}(\boldsymbol{V})(x)),\forall
x\in\partial\Omega.\end{cases}$

In the following texts, we denote
$\boldsymbol{L}^{u}\mathrel{\mathop{\mathchar
58\relax}}=(\boldsymbol{L}^{u,\Omega},\boldsymbol{L}^{u,\partial\Omega})$,
$\boldsymbol{L}^{m}\mathrel{\mathop{\mathchar
58\relax}}=(\boldsymbol{L}^{m,\Omega},\boldsymbol{L}^{m,\partial\Omega})$, and
$\boldsymbol{L}^{\boldsymbol{{V}}}\mathrel{\mathop{\mathchar
58\relax}}=(\boldsymbol{L}^{\boldsymbol{{V}},\Omega},\boldsymbol{L}^{\boldsymbol{{V}},\partial\Omega})$.

Under Assumption 2, we rewrite (4.3) as

$\displaystyle\begin{cases}\min\limits_{\overset{(u_{k},m_{k},\boldsymbol{V}_{k})\in\mathcal{U}_{k}\times\mathcal{M}_{k}\times\boldsymbol{\mathcal{V}}_{k}}{k=0,\dots,N_{T}}}\sum_{k=0}^{N_{T}}\|u_{k}\|_{\mathcal{U}_{k}}^{2}+\sum_{k=0}^{N_{T}}\|m_{k}\|_{\mathcal{M}_{k}}^{2}+\sum_{k=0}^{N_{T}}\|\boldsymbol{V}_{k}\|_{\boldsymbol{\mathcal{V}}_{k}}^{2}\\\
\quad\quad\quad\quad\quad\quad\quad\quad+\sum_{k=0}^{N_{T}}\frac{1}{\gamma^{2}_{k}}|[\boldsymbol{\phi}^{o},m_{k}]-\boldsymbol{m}^{o}_{k}|^{2}+\sum_{k=0}^{N_{T}}|\Sigma^{-1}_{k}([\Psi_{k},\boldsymbol{V}_{k}]-\boldsymbol{V}^{o}_{k})|^{2}\\\
\operatorname{s.t.}\quad-\frac{[\delta_{x_{i}},u_{k+1}-u_{k}]}{\Delta
t}+P_{1}([\delta_{x_{i}}\circ\boldsymbol{L}^{\Omega,u},u_{k}],[\delta_{x_{i}}\circ\boldsymbol{L}^{\Omega,m},m_{k+1}],[\delta_{x_{i}}\circ\boldsymbol{L}^{\Omega,\boldsymbol{V}},\boldsymbol{V}_{k}])=0,\quad\\\
\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\forall
i\in\\{1,\dots,M_{\Omega}\\},k\in\\{0,\dots,N_{T}-1\\},\\\
\quad\quad\quad\frac{[\delta_{x_{i}},m_{k+1}-m_{k}]}{\Delta
t}+P_{2}([\delta_{x_{i}}\circ\boldsymbol{L}^{\Omega,u},u_{k}],[\delta_{x_{i}}\circ\boldsymbol{L}^{\Omega,m},m_{k+1}],[\delta_{x_{i}}\circ\boldsymbol{L}^{\Omega,\boldsymbol{V}},\boldsymbol{V}_{k}])=0,\\\
\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\forall
i\in\\{1,\dots,M_{\Omega}\\},k\in\\{0,\dots,N_{T}-1\\},\\\ \quad\quad\quad
B([\delta_{x_{i}}\circ\boldsymbol{L}^{\partial\Omega,u},u],[\delta_{x_{i}}\circ\boldsymbol{L}^{\partial\Omega,m},m],[\delta_{x_{i}}\circ\boldsymbol{L}^{\partial\Omega,\boldsymbol{V}},\boldsymbol{V}])=0,\\\
\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\forall
i\in\\{M_{\Omega}+1,\dots,M\\},k\in\\{0,\dots,N_{T}-1\\},\\\
\quad\quad\quad[\delta_{x_{i}},u_{N_{T}}]=\phi^{*}(x_{i}),[\delta_{x_{i}},m_{0}]=\mu^{*}(x_{i}),\,\,\quad\quad\forall
i\in\\{1,\dots,M\\},\end{cases}$ (4.4)

Denote $\boldsymbol{\delta}\mathrel{\mathop{\mathchar
58\relax}}=(\boldsymbol{\delta}^{\Omega},\boldsymbol{\delta}^{\partial\Omega})$,
where
$\boldsymbol{\delta}^{\Omega}=(\delta_{x_{1}},\dots,\delta_{x_{M_{\Omega}}})$
and
$\boldsymbol{\delta}^{\partial\Omega}=(\delta_{x_{M_{\Omega}}+1},\dots,\delta_{x_{M}})$.
Let
$\widetilde{\boldsymbol{\phi}}^{u}=\boldsymbol{\delta}\circ\boldsymbol{L}^{u}$,
$\widetilde{\boldsymbol{\phi}}^{m}=\boldsymbol{\delta}\circ\boldsymbol{L}^{m}$,
and
$\widetilde{\boldsymbol{\phi}}^{\boldsymbol{V}}=\boldsymbol{\delta}\circ\boldsymbol{L}^{\boldsymbol{V}}$.
To solve (4.4), we introduce latent variables and rewrite (4.4) as

$\displaystyle\begin{cases}\min\limits_{\boldsymbol{z},\boldsymbol{\rho},\boldsymbol{v}}\begin{cases}\min\limits_{\overset{(u_{k},m_{k},\boldsymbol{V}_{k})\in\mathcal{U}_{k}\times\mathcal{M}_{k}\times\boldsymbol{\mathcal{V}}_{k}}{k=0,\dots,N_{T}}}\sum_{k=0}^{N_{T}}\|u_{k}\|_{\mathcal{U}_{k}}^{2}+\sum_{k=0}^{N_{T}}\|m_{k}\|_{\mathcal{M}_{k}}^{2}+\sum_{k=0}^{N_{T}}\|\boldsymbol{V}_{k}\|_{\boldsymbol{\mathcal{V}}_{k}}^{2}\\\
\text{s.t.}\quad[\boldsymbol{\delta}^{\Omega},u_{k}]=\boldsymbol{z}^{(1),\Omega}_{k},[\boldsymbol{\delta}^{\partial\Omega},u_{k}]=\boldsymbol{z}^{(1),\partial\Omega}_{k},[\widetilde{\boldsymbol{\phi}}^{u},u_{k}]=\boldsymbol{z}^{(2)}_{k},\\\
\quad\quad[\boldsymbol{\delta}^{\Omega},m_{k}]=\boldsymbol{\rho}^{(1),\Omega}_{k},[\boldsymbol{\delta}^{\partial\Omega},m_{k}]=\boldsymbol{\rho}^{(1),\partial\Omega}_{k},[\boldsymbol{\phi}^{o},m_{k}]=\boldsymbol{\rho}^{(2)}_{k},[\widetilde{\boldsymbol{\phi}}^{m},m_{k}]=\boldsymbol{\rho}^{(3)}_{k},\\\
\quad\quad[\widetilde{\boldsymbol{\phi}}^{\boldsymbol{V}},\boldsymbol{V}_{k}]=\boldsymbol{v}^{(1)}_{k},[\Psi,\boldsymbol{V}_{k}]=\boldsymbol{v}^{(2)}_{k},\end{cases}\\\
\quad\quad\quad\quad+\sum_{k=0}^{N_{T}}\frac{1}{\gamma^{2}_{k}}|\boldsymbol{\rho}^{(2)}_{k}-\boldsymbol{m}^{o}_{k}|^{2}+\sum_{k=0}^{N_{T}}|\Sigma^{-1}_{k}(\boldsymbol{v}^{(2)}_{k}-\boldsymbol{V}^{o}_{k})|^{2}\\\
\operatorname{s.t.}\quad-\frac{z_{k+1,i}^{(1),\Omega}-z_{k,i}^{(1),\Omega}}{\Delta
t}+P_{1}(\boldsymbol{z}_{k}^{(2),\Omega},\boldsymbol{\rho}_{k+1}^{(3),\Omega},\boldsymbol{v}_{k}^{(1),\Omega})=0,\quad\forall
i\in\\{1,\dots,M_{\Omega}\\},k\in\\{0,\dots,N_{T}-1\\},\\\
\quad\quad\quad\frac{\rho_{k+1,i}^{(1),\Omega}-\rho_{k,i}^{(1),\Omega}}{\Delta
t}+P_{2}(\boldsymbol{z}_{k}^{(2),\Omega},\boldsymbol{\rho}_{k+1}^{(3),\Omega},\boldsymbol{v}_{k}^{(1),\Omega})=0,\quad\,\,\forall
i\in\\{1,\dots,M_{\Omega}\\},k\in\\{0,\dots,N_{T}-1\\},\\\ \quad\quad\quad
B(\boldsymbol{z}_{k}^{(2),\partial\Omega},\boldsymbol{\rho}_{k+1}^{(3),\partial\Omega},\boldsymbol{v}_{k}^{(1),\partial\Omega})=0,\quad\quad\quad\quad\quad\quad\,\,\,\forall
i\in\\{M_{\Omega}+1,\dots,M\\},k\in\\{0,\dots,N_{T}-1\\},\\\ \quad\quad\quad
z_{N_{T},i}^{(1)}=\phi^{*}(x_{i}),\rho_{0,i}^{(1)}=\mu^{*}(x_{i}),\quad\quad\quad\quad\quad\quad\quad\,\,\forall
i\in\\{1,\dots,M\\},\end{cases}$ (4.5)

where
$\boldsymbol{z}=(\boldsymbol{z}_{k})_{k=0}^{N_{T}},\boldsymbol{z}_{k}=(\boldsymbol{z}^{(1),\Omega}_{k},\boldsymbol{z}^{(1),\partial\Omega}_{k},\boldsymbol{z}^{(2)}_{k})$,
$\boldsymbol{\rho}=(\boldsymbol{\rho}_{k})_{k=0}^{N_{T}},\boldsymbol{\rho}_{k}=(\boldsymbol{\rho}^{(1),\Omega}_{k},\boldsymbol{\rho}^{(1),\partial\Omega}_{k},\boldsymbol{\rho}^{(2)}_{k},\boldsymbol{\rho}^{(3)}_{k})$,
and
$\boldsymbol{v}=(\boldsymbol{v}_{k})_{k=0}^{N_{T}},\boldsymbol{v}_{k}=(\boldsymbol{v}^{(1)}_{k},\boldsymbol{v}^{(2)}_{k})$.
Denote $\boldsymbol{\phi}^{u}\mathrel{\mathop{\mathchar
58\relax}}=(\boldsymbol{\delta},\widetilde{\boldsymbol{\phi}}^{u})$,
$\boldsymbol{\phi}^{m}\mathrel{\mathop{\mathchar
58\relax}}=(\boldsymbol{\delta},\boldsymbol{\phi}^{o},\widetilde{\boldsymbol{\phi}}^{m})$,
and $\boldsymbol{\phi}^{\boldsymbol{V}}\mathrel{\mathop{\mathchar
58\relax}}=(\widetilde{\boldsymbol{\phi}}^{\boldsymbol{V}},\Psi)$.

Let $\\{K_{u,k}\\}_{k=0}^{N_{T}}$, $\\{K_{m,k}\\}_{k=0}^{N_{T}}$,
$\\{K_{\boldsymbol{V},k}\\}_{k=0}^{N_{T}}$ be collections of kernels
corresponding to the sets of RKHSs $\\{\mathcal{U}_{k}\\}_{k=0}^{N_{T}}$,
$\\{\mathcal{M}_{k}\\}_{k=0}^{N_{T}}$, and
$\\{\boldsymbol{\mathcal{V}}_{k}\\}_{k=0}^{N_{T}}$ respectively. Let
$(u_{k}^{\dagger},m_{k}^{\dagger},\boldsymbol{V}_{k}^{\dagger})_{k=0}^{N_{T}}$
be a minimizer to the first level minimization problem in (4.5). Then, by the
representer formulas in (2.4) and (2.5) (also see [49, Sec. 17.8]), we obtain

$\displaystyle\begin{cases}u_{k}^{\dagger}(x)=K_{u,k}(x,\boldsymbol{\phi}^{u})K_{u,k}(\boldsymbol{\phi}^{u},\boldsymbol{\phi}^{u})^{-1}\boldsymbol{z}_{k},\\\
m_{k}^{\dagger}(x)=K_{m,k}(x,{\boldsymbol{\phi}}^{m})K_{m,k}({\boldsymbol{\phi}}^{m},{\boldsymbol{\phi}}^{m})^{-1}\boldsymbol{\rho}_{k},\\\
\boldsymbol{V}_{k}^{\dagger}(x)={K}_{\boldsymbol{V},k}(x,{\boldsymbol{\phi}}^{\boldsymbol{V}}){K}_{\boldsymbol{V},k}({\boldsymbol{\phi}}^{\boldsymbol{V}},{\boldsymbol{\phi}}^{\boldsymbol{V}})^{-1}\boldsymbol{v}_{k}.\end{cases}$

Thus, we have

$\displaystyle\begin{cases}\|u_{k}^{\dagger}\|_{\mathcal{U}_{k}}^{2}=\boldsymbol{z}_{k}^{T}K_{u,k}(\boldsymbol{\phi}^{u},\boldsymbol{\phi}^{u})^{-1}\boldsymbol{z}_{k},\\\
\|m_{k}^{\dagger}\|_{\mathcal{M}_{k}}^{2}=\boldsymbol{\rho}_{k}^{T}K_{m,k}({\boldsymbol{\phi}}^{m},{\boldsymbol{\phi}}^{m})^{-1}\boldsymbol{\rho}_{k},\\\
\|\boldsymbol{V}_{k}^{\dagger}\|_{\boldsymbol{\mathcal{V}}_{k}}^{2}=\boldsymbol{v}_{k}^{T}{K}_{\boldsymbol{V},k}({\boldsymbol{\phi}}^{\boldsymbol{V}},{\boldsymbol{\phi}}^{\boldsymbol{V}})^{-1}\boldsymbol{v}_{k}.\end{cases}$

Hence, we can formulate (4.5) as a finite-dimensional optimization problem

$\displaystyle\begin{cases}\min\limits_{\boldsymbol{z},\boldsymbol{\rho},\boldsymbol{v}}\sum_{k=0}\boldsymbol{z}_{k}^{T}K_{u,k}(\boldsymbol{\phi}^{u},\boldsymbol{\phi}^{u})^{-1}\boldsymbol{z}_{k}+\sum_{k=0}\boldsymbol{\rho}_{k}^{T}K_{m,k}(\boldsymbol{\phi}^{m},\boldsymbol{\phi}^{m})^{-1}\boldsymbol{\rho}_{k}+\sum_{k=0}\boldsymbol{v}_{k}^{T}{K}_{\boldsymbol{V},k}(\boldsymbol{\phi}^{\boldsymbol{V}},\boldsymbol{\phi}^{\boldsymbol{V}})^{-1}\boldsymbol{v}_{k}\\\
\quad\quad\quad\quad+\sum_{k=0}^{N_{T}}\frac{1}{\gamma^{2}_{k}}|\boldsymbol{\rho}^{(2)}_{k}-\boldsymbol{m}^{o}_{k}|^{2}+\sum_{k=0}^{N_{T}}|\Sigma^{-1}_{k}(\boldsymbol{v}^{(2)}_{k}-\boldsymbol{V}^{o}_{k})|^{2}\\\
\operatorname{s.t.}\quad-\frac{z_{k+1,i}^{(1),\Omega}-z_{k,i}^{(1),\Omega}}{\Delta
t}+P_{1}(\boldsymbol{z}_{k}^{(2),\Omega},\boldsymbol{\rho}_{k+1}^{(3),\Omega},\boldsymbol{v}_{k}^{(1),\Omega})=0,\quad\forall
i\in\\{1,\dots,M_{\Omega}\\},k\in\\{0,\dots,N_{T}-1\\},\\\
\quad\quad\quad\frac{\rho_{k+1,i}^{(1),\Omega}-\rho_{k,i}^{(1),\Omega}}{\Delta
t}+P_{2}(\boldsymbol{z}_{k}^{(2),\Omega},\boldsymbol{\rho}_{k+1}^{(3),\Omega},\boldsymbol{v}_{k}^{(1),\Omega})=0,\quad\,\,\forall
i\in\\{1,\dots,M_{\Omega}\\},k\in\\{0,\dots,N_{T}-1\\},\\\ \quad\quad\quad
B(\boldsymbol{z}_{k}^{(2),\partial\Omega},\boldsymbol{\rho}_{k+1}^{(3),\partial\Omega},\boldsymbol{v}_{k}^{(1),\partial\Omega})=0,\quad\quad\quad\quad\quad\quad\,\,\,\forall
i\in\\{M_{\Omega}+1,\dots,M\\},k\in\\{0,\dots,N_{T}-1\\},\\\ \quad\quad\quad
z_{N_{T},i}^{(1)}=\phi^{*}(x_{i}),\rho_{0,i}^{(1)}=\mu^{*}(x_{i}),\quad\quad\quad\quad\quad\quad\quad\,\,\forall
i\in\\{1,\dots,M\\}.\end{cases}$ (4.6)

In the numerical experiments presented in Section 5, we use the Gauss–Newton
method to solve (4.6).

##  5\. Numerical Experiments

In this section, we illustrate the effectiveness of our approach using various
examples. All experiments involve reconstructing complete population profiles,
agent strategies, and entire environmental configurations from partial, noisy
data about populations and environments. Specifically, Subsection 5.1
highlights the determination of the viscosity constant, a critical parameter
measuring uncertainty in agents’ dynamics. Subsections 5.2 and 5.3 are
dedicated to identifying the coupling functions in MFGs, which encode agents’
interactions. Subsection 5.4 addresses a misspecification problem,
demonstrating our method’s robustness by recovering configurations of a MFG
with local coupling using a model designed for non-local coupling. Lastly, in
Subsection 5.5, we apply our methodology to a time-dependent MFG scenario,
showcasing its adaptability to dynamic settings.

In all of our experiments, we employ Python in conjunction with the JAX
package, which facilitates automatic differentiation. The current setup relies
solely on CPU processing. However, notable enhancements in performance can be
realized by incorporating accelerated hardware like Graphics Processing Units
(GPUs), which are designed to significantly expedite computational tasks.

###  5.1. A Stationary MFG

In this example, we study inverse problems related to the following stationary
MFG

$\displaystyle\begin{cases}-\nu\Delta u+H(x,\nabla
u)=m^{2}-V(x)+\overline{H},&\forall x\in\mathbb{T}^{2},\\\ -\nu\Delta
m-\operatorname{div}(D_{p}H(x,\nabla u)m)=0,&\forall x\in\mathbb{T}^{2},\\\
\int_{\mathbb{T}^{2}}u\,dx=0,\int_{\mathbb{T}^{2}}m\,dx=1,&\end{cases}$ (5.1)

where $H(x,p)=\frac{1}{2}|p|^{2}$. We identify $\mathbb{T}^{2}$ by
$[-0.5,0.5)\times[-0.5,0.5)$, and the true potential $V(x)=\sin(2\pi
y)+\sin(2\pi x)+\cos(4\pi x)$.

Our first goal is to recover $(u,m,V,\overline{H})$ given noisy observations
of the population density $m$ under a fixed viscosity $\nu$. To verify the
accuracy of the recovered quantities, we first solve (5.1) using the GP method
in [47] given $\nu$ and $V$. These numerical solutions, alongside the explicit
formulation of $V$, are viewed as reference solutions. This designation
underscores their role as benchmarks against which we evaluate the accuracy of
the corresponding recovered solutions. For simplicity, we denote the reference
solutions by $(u^{*},m^{*},V^{*},\overline{H}^{*})$. To solve the inverse
problem, we randomly select $M=400$ collocation points from $\mathbb{T}^{2}$
and choose the first $I=40$ points as observation points. We treat the
addition of the values of $m^{*}$ at the first $I$ observation points with
independent standard Gaussian noises following $\mathcal{N}(0,\gamma^{2}I)$,
where $\gamma=10^{-3}$, as the observed data for $m$. We employed periodic
kernels for recovering $m$, $u$, and $V$. For each kernel, we set the
lengthscale parameter to $1.41$, and incorporated adaptive diagonal
regularization terms (“nuggets”) with parameters $\eta=10^{-8}$ to the
covariance matrices. Figure 1 presents the results of our experiment in
recovering $m$, $u$, and $V$ with a viscosity $\nu=0.1$. We notice that the
accuracy in recovering the variables $u$ and $m$ is superior to that of $V$
and $\overline{H}$. The reason behind this lies in the nature of the variables
and their interdependencies. Given known values of density $m$ and viscosity
$\nu$, the solution for $u$ can be precisely determined through the FP
equation and integration constraints (see discussions in the introduction
section). In contrast, the recovery of the pair $(V,\overline{H})$ faces a
uniqueness challenge. This challenge arises because if a pair
$(V,\overline{H})$ satisfies the equation in (5.1), then so does the pair
$(V+C,\overline{H}+C)$ for any constant $C$. This implies that without direct
observations of $V$, it’s impossible to uniquely determine the values of both
$V$ and $\overline{H}$. Furthermore, we observe that the pointwise error in
estimating $V-\overline{H}$ is less than that in estimating $V$ alone. The
$L_{2}$ error for the recovered $V$ is $0.311662$, in contrast, the $L_{2}$
error for the recovered $V-\overline{H}$ is significantly lower, recorded at
$0.0035159$. Hence, it is promising to accurately recover $V$ up to a constant
when $\nu$ and only noisy data on $m$ are given.

Figure 2 illustrates the log-log plot of error versus the number of observed
data when $\nu=0.1$. We see that the errors decrease as we increase the number
of observations. The results are obtained by averaging over 50 independent
realizations.

In scenarios where viscosity $\nu$ is unknown, we demonstrate the feasibility
of simultaneously recovering both $\nu$ and $V$ from observations of the
population $m$ and $V$. Table 1 shows the recovered $\nu$ values under various
true viscosities. We choose $M=400$ sample points, out of which $I=20$ data
points of $m$ and $V$ are noisily observed with additive noises following
$\mathcal{N}(0,\gamma^{2}I)$, where $\gamma=10^{-3}$. We see that our GP
framework can accurately recover the uncertainty in the MFG system.

a Samples & observations

b error contour of $V-\bar{H}$

c error contour of $m$

d error contour of $u$

e error contour of $V$

f baseline for $m$

g baseline for $u$

h baseline for $V$

i recovered $m$

j recovered $u$

k recovered $V$

Fig. 1. Numerical results for the MFG in (5.1). $\nu=0.1$, $I=40$ observed
data on $m$, $M=400$ collocation points, ground truth
$\overline{H}=-1.066853$, recovered $\overline{H}^{\dagger}=-0.752241$. The
$L_{2}$ norm of error of $V$ is $0.311662$, and the $L_{2}$ norm of error of
$V-\overline{H}$ is $0.035159$.

a log-log plot of $L_{2}$ error of $m$

b log-log plot of $L_{2}$ error of $u$

c log-log plot of $L_{2}$ error of $V$

Fig. 2. Numerical results for the MFG in (5.1). Log-log plots of errors v.s. number of observed data $I$, when $\nu=0.1$.  Table 1. Recovered viscosity constants for the MFG in (5.1) with $20$ points of observations (both $m$ and $V$ are observed), and $M=400$ collocation points. $\nu$ stands for the true viscosity and $\nu^{\dagger}$ represents the recovered viscosity.  $\nu$ | 1 | 0.5 | 0.1 | 0.05 | 0.01
---|---|---|---|---|---
$\nu^{\dagger}$ | 1.000804 | 0.497811 | 0.099312 | 0.050905 | 0.015023

###  5.2. A Stationary MFG with an Unknown Local Coupling

We investigate a stationary MFG with an unknown local coupling, formulated as:

$\displaystyle\begin{cases}-\Delta u+H(x,\nabla
u)=\Gamma(m)-V(x)+\overline{H},&\forall x\in\mathbb{T}^{2},\\\ -\Delta
m-\operatorname{div}(D_{p}H(x,\nabla u)m)=0,&\forall x\in\mathbb{T}^{2},\\\
\int_{\mathbb{T}^{2}}u\,dx=0,\int_{\mathbb{T}^{2}}m\,dx=1,&\end{cases}$ (5.2)

where $H(x,p)=\frac{1}{2}|p|^{2}$, $\Gamma(m)=m^{\alpha}$, and the exponent
$\alpha$ is a positive constant. We identify $\mathbb{T}^{2}$ with
$[-0.5,0.5)\times[-0.5,0.5)$, and choose the true potential $V(x)=\cos(2\pi
x)+\sin(2\pi y)+\sin(4\pi x)$. Our objective is to recover both $V$ and
$\alpha$ based on noisy observations of the population density $m$ and $V$.

As before, we use the GP method in [47] to compute the numerical solutions for
(5.2) for given $\alpha$ and $V$. These numerical solutions, in conjunction
with the explicit formulation of $V$ and $\alpha$ serve as reference solutions
for accuracy comparisons. We randomly select $M=400$ collocation points from
$\mathbb{T}^{2}$, observing the values of $m$ at the first $I=40$ points and
the values of $V$ at the first $I_{V}=20$ points as data. Independent Gaussian
noises following $\mathcal{N}(0,\gamma^{2}I)$ with a standard deviation of
$\gamma=10^{-3}$ are added to these observations. Periodic kernels are
deployed to recover $m$, $u$, and $V$, with lengthscale parameters being 0.6,
0.6, and 1, respectively. Additionally, adaptive diagonal regularization terms
(’nuggets’) with parameters $\eta=10^{-8}$ are incorporated into the
covariance matrices for regularization purposes.

Figure 3 presents our experimental results in recovering $m$, $u$, and $V$
with $\alpha=2$. The results show that we recover $u,m,V,\overline{H}$, and
the coupling $\Gamma$ within reasonable accuracy with limited observational
data. Table 2 displays the recovered values of $\alpha$ and $\overline{H}$
under various ground truths, with $I_{V}=20$ observed data points for $V$,
$I=40$ observed data points for $m$, and $M=400$ samples.

a Samples & observations

b error contour of $m$

c error contour of $u$

d error contour of $V$

e baseline for $\Gamma$

f baseline for $m$

g baseline for $u$

h baseline for $V$

i recovered $\Gamma$

j recovered $m$

k recovered $u$

l recovered $V$

Fig. 3. Numerical results for the MFG in (5.2) with $\Gamma(m)=m^{\alpha}$. $I=40$ observed data on $m$ and $I_{V}=20$ observed data on $V$, $M=400$ collocation points, ground truth $\alpha=2$, $\overline{H}=-0.975970$, recovered $\alpha^{\dagger}=1.990808$, $\overline{H}^{\dagger}=-0.966401$. Table 2. Numerical results for the MFG in (5.2) with $\Gamma(m)=m^{\alpha}$. Recovered power $\alpha$ and $\overline{H}$ with $I_{V}=20$ data points of $V$, $I=40$ observations on $m$, and $M=400$ sample points. $\alpha$ and $\overline{H}$ stand for the references and $\alpha^{\dagger}$ and $\overline{H}^{\dagger}$ represent the recovered variables.  $\alpha$ | 0.5 | 1 | 2 | 3 | 5
---|---|---|---|---|---
$\alpha^{\dagger}$ | 0.504769 | 1.182533 | 1.990808 | 2.89700 | 4.776276
$\overline{H}$ | -0.970467 | -0.971760 | -0.975970 | -0.985608 | -0.998791
$\overline{H}^{\dagger}$ | -0.970979 | -0.991179 | -0.966402 | -0.985758 | -0.981125

###  5.3. A Stationary MFG with an Unknown Non-local Coupling

In this example, we study the stationary MFG (5.2) with a non-local coupling
$\Gamma$ defined as

$\Gamma(m)=\int_{\mathbb{T}^{2}}e^{-\frac{(x-y)^{2}}{2\sigma^{2}}}m(y)\,\operatorname{d\\!}y,$

where the lengthscale $\sigma$ is unknown. We identify $\mathbb{T}^{2}$ with
$[0,1)\times[0,1)$, and the true potential $V(x)=\sin(2\pi y)+\sin(2\pi
x)+\cos(4\pi x)$. Our objective is to recover $m,u,V,\overline{H}$ and
$\sigma$ based on noisy observations of the population density $m$ and
coefficient $V$.

We use the GP method in [47] to compute the reference solutions given $\sigma$
and $V$. Then, we take $M=950$ sample points in domain, of which $900$ points
are Gauss-Legendre quadrature points and $50$ points are randomly selected
observation points for $m$ and $V$, i.e., $I=I_{V}=50$. Independent Gaussian
noises following $\mathcal{N}(0,\gamma^{2}I)$ with a standard deviation of
$\gamma=10^{-3}$ are added to these observations. The Gauss-Legendre
quadrature rule is utilized to compute the convolution in the non-local
coupling term. Periodic kernels are chosen for $m$, $u$, and $V$ with
lengthscale parameters $0.6$, $0.6$, and $1$, respectively. Additionally,
adaptive diagonal regularization terms (“nuggets”) with parameters
$\eta=10^{-8}$ are incorporated into the covariance matrices.

Figure 4 presents our experimental results in recovering $m$, $u$, and $V$
with $\sigma=0.5$. Table 3 displays the recovered values of $\sigma$ and
$\overline{H}$ under various ground truths, with $I_{V}=I=50$ observed $m,V$
data points and a total of $M=950$ sample points. These results demonstrate
the accuracy of our GP framework.

a Samples & observations

b error contour of $m$

c error contour of $u$

d error contour of $V$

e baseline for $\Gamma$

f baseline for $m$

g baseline for $u$

h baseline for $V$

i recovered $\Gamma$

j recovered $m$

k recovered $u$

l recovered $V$

Fig. 4. Numerical results for the MFG in (5.2) with $\Gamma(m)=\int_{\mathbb{T}^{2}}e^{\frac{(x-y)^{2}}{2\sigma^{2}}}m(y)\operatorname{d\\!}y$. $50$ observation points for $m$ and $V$, $M=950$ sample points, ground truth $\sigma=0.5$, $\overline{H}=-0.545407$, recovered $\sigma^{\dagger}=0.504769$, $\overline{H}^{\dagger}=-0.626463$. Table 3. Numerical results for the MFG in (5.2) with $\Gamma(m)=\int_{\mathbb{T}^{2}}e^{\frac{(x-y)^{2}}{2\sigma^{2}}}m(y)\operatorname{d\\!}y$. Recovered values of $\sigma$ and $\overline{H}$ with $50$ observation points for $m$ and $V$ and $M=950$ sample points.  $\sigma$ | 0.4 | 0.5 | 1 | 2
---|---|---|---|---
$\sigma^{\dagger}$ | 0.541382 | 0.504769 | 0.978857 | 1.921419
$\overline{H}$ | -0.431692 | -0.545408 | -0.819267 | -0.926975
$\overline{H^{\dagger}}$ | -0.584162 | -0.626463 | -0.814232 | -0.921420

###  5.4. A Misspecified Recovery

In this example, we study a misspecification problem, where we recover the
strategies and environment data generated by a MFG with a local coupling by
specifying another MFG model with a non-local coupling. More precisely, we
observe noisy data generated by the MFG system in (5.2) with $V(x)=\cos(2\pi
x)+\cos(2\pi y)+\cos(4\pi x)$ and $\Gamma(m)=m^{2}$, but our goal is to
recover $m$, $u$, $V$, and $\overline{H}$ by assuming that $\Gamma(m)$ takes
the form
$\Gamma(m)=\int_{\mathbb{T}^{2}}e^{-\frac{(x-y)^{2}}{2\sigma^{2}}}m(y)\,dy$,
where the lengthscale $\sigma$ is to be determined.

We use the GP method [47] to solve (5.2) with $\Gamma(m)=m^{2}$ to get the
reference solutions. Then, to solve the inverse problems, we take $625$ Gauss-
Legendre quadrature points as samples and randomly select $I=30$ observation
points for $m$. Independent Gaussian noises following
$\mathcal{N}(0,\gamma^{2}I)$ with a standard deviation of $\gamma=10^{-3}$ are
added to these observations. The Gauss-Legendre quadrature rule are employed
to calculate the convolution in the non-local coupling. In the numerical
experiments, periodic kernels are chosen for $m$, $u$, and $V$ with
lengthscale parameters set to $0.6$, $0.6$, and $1$, respectively.
Additionally, adaptive diagonal regularization terms (“nuggets”) with
parameters $\eta=10^{-8}$ are added into the covariance matrices.

Figure 5 presents our experimental results in recovering $m$, $u$, and $V$. It
is important to highlight that in this scenario, there are no observations
available for $V$. The accuracy of the results imply the robustness of our GP
framework.

a Samples & observations

b error contour of $m$

c error contour of $u$

d error contour of $V$

e baseline for $m$

f baseline for $u$

g baseline for $V$

h recovered $m$

i recovered $u$

j recovered $V$

Fig. 5. Numerical results for the misspecification problem. We choose $M=625$
samples and $I=30$ observations on $m$. The ground truth
$\overline{H}=-0.920221$, and the recovered
$\overline{H}^{\dagger}=-0.524876$, recovered $\sigma=0.619232$.

###  5.5. A Time-dependent MFG

In this subsection, we investigate a time-dependent MFG with local coupling,
represented by the following equations:

$\displaystyle\begin{cases}-\partial_{t}u-\nu\Delta u+H(x,\nabla
u)=m^{2}-V(x),&\forall(t,x)\in(0,T)\times\mathbb{T},\\\
\partial_{t}m-\nu\Delta m-\operatorname{div}(D_{p}H(x,\nabla
u)m)=0,&\forall(t,x)\in(0,T)\times\mathbb{T},\\\ m(0,x)=m_{0}(x),\quad
u(T,x)=u_{T}(x),&\forall x\in\mathbb{T},\end{cases}$

where $T=1$, $H(x,p)=\frac{1}{2}|p|^{2}$. Let $\mathbb{T}$ be identified by
$[-0.5,0.5)$, $m_{0}(\cdot)=1$, and $u_{T}(\cdot)=0$, and the true potential
$V(x)=2\sin(2\pi x)+\cos(2\pi x)$. Our goal is to recover $u,m,V$ in the whole
time and space domain under noisy observations of the population density $m$
and of the potential $V$ for a fixed viscosity $\nu$.

Following the methodology outlined in Section 4.2, we divide the entire time
period into $N_{T}=20$ intervals. At each interval, we randomly sample $M=50$
collocation points from $\mathbb{T}$ and observe the value of $m$ at the first
$I_{m}=5$ points as data, along with observations of $V$ at $I_{V}=3$ points.
Independent Gaussian noises following $\mathcal{N}(0,\gamma^{2}I)$ with a
standard deviation of $\gamma=10^{-3}$ are added to these observations.
Periodic kernels are used for $m$, $u$, and $V$, with lengthscale parameters
$1.5,1.5$, and $1$, separately. Additionally, we set the nugget parameters
$\eta=10^{-8}$ to regularize the covariance matrices.

Figure 6 presents our experimental results in recovering $m$, $u$, and $V$
with a viscosity of $\nu=0.1$. The recovery process demonstrated reasonable
accuracy, effectively recover the quantities of interest with limited
observations.

a error contour of $m$

b error contour of $u$

c error of $V$

d baseline for $m$

e baseline for $u$

f baseline for $V$, the $3$ scatter points are observations on $V$

g recovered $m$

h recovered $u$

i recovered $V$, the $3$ scatter points are observations on $V$

Fig. 6. Numerical results for the time-dependent MFG. We set $\nu=0.1$. At
each time slot, we select $M=50$ samples on $\mathbb{T}$. Within these
samples, we choose $I_{m}=5$ observation points for $m$ and $I_{V}=3$
observation points for $V$.

##  6\. Conclusions and Future Works

In this study, we introduce a GP framework designed to infer strategies and
environmental configurations for agents engaged in MFGs. This framework
operates under the challenge of working with partial and noisy data regarding
the populations of agents and their environmental settings. To tackle this, we
employ GPs to approximate the quantities of interest and address the inverse
problems through Maximum A Posteriori (MAP) estimations conditioned on PDEs
and constraints that are satisfied at specific collocation points and within
the observed data.

Our method’s effectiveness is showcased through a variety of examples,
highlighting its capability to accurately recover diverse quantities in MFG
contexts. An intriguing avenue for future work involves integrating approaches
that address scalability issues, such as Random Fourier Features, sparse GPs,
and mini-batch methods. This integration aims to enhance the scalability of
our framework, making it more adept at handling large datasets.

Furthermore, another promising direction is the application of our methodology
in fields like economics, biology, or finance. Here, the goal would be to
develop MFG models in scenarios where no established models matching the data
currently exist, using our approach to interpret and manage the data
effectively.

###### Acknowledgement A0.

CM is supported by CityU Start-up Grant 7200684, Hong Kong RGC Grant ECS
21302521, Hong Kong RGC Grant GRF 11311422 and Hong Kong RGC Grant GRF
11303223. XY acknowledges support from the Air Force Office of Scientific
Research under MURI award number FA9550-20-1-0358 (Machine Learning and
Physics-Based Modeling and Simulation). CZ is supported by Singapore MOE
(Ministry of Education) AcRF Grant A-8000453-00-00, IoTex Foundation Industry
Grant A-8001180-00-00 and NSFC Grant No. 11871364.

## References

  * [1] Y. Achdou and I. Capuzzo-Dolcetta.  Mean field games: numerical methods.  SIAM Journal on Numerical Analysis, 48(3):1136–1162, 2010.
  * [2] Y. Achdou and M. Laurière.  Mean field games and applications: Numerical aspects.  Mean Field Games, pages 249–307, 2020.
  * [3] Y. Achdou and V. Perez.  Iterative strategies for solving linearized discrete mean field games systems.  Networks and Heterogeneous Media, 7(2):197, 2012.
  * [4] M. A. Alvarez, L. Rosasco, and N. D. Lawrence.  Kernels for vector-valued functions: A review.  Foundations and Trends® in Machine Learning, 4(3):195–266, 2012.
  * [5] L. M. Briceño-Arias, D. Kalise, Z. Kobeissi, M. Laurière, A. M. González, and F. Silva.  On the implementation of a primal-dual algorithm for second order time-dependent mean field games with local couplings.  ESAIM: Proceedings and Surveys, 65:330–348, 2019.
  * [6] L. M. Briceno-Arias, D. Kalise, and F. Silva.  Proximal methods for stationary mean field games with local couplings.  SIAM Journal on Control and Optimization, 56(2):801–836, 2018.
  * [7] P. Cannarsa, W. Cheng, C. Mendico, and K. Wang.  Long-time behavior of first-order mean field games on euclidean space.  Dynamic Games and Applications, 10(2):361–390, 2020.
  * [8] P. Cardaliaguet.  Notes on mean field games.  Technical report, Technical report, 2010.
  * [9] P. Cardaliaguet, J.-M. Lasry, P.-L. Lions, and A. Porretta.  Long time average of mean field games.  Networks and heterogeneous media, 7(2):279–301, 2012.
  * [10] 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.  SIAM Journal on Numerical Analysis, 59(3):1455–1485, 2021.
  * [11] 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.  accepted in Annals of Applied Probability, 2021.
  * [12] A. Cesaroni, M. Cirant, S. Dipierro, M. Novaga, and E. Valdinoci.  On stationary fractional mean field games.  Journal de Mathématiques Pures et Appliquées, 122:1–22, 2019.
  * [13] Y. Chen, B. Hosseini, H. Owhadi, and A. M. Stuart.  Solving and learning nonlinear PDEs with Gaussian processes.  Journal of Computational Physics, 2021.
  * [14] Y. Chen, H. Owhadi, and F. Schäfer.  Sparse Cholesky factorization for solving nonlinear PDEs via Gaussian processes.  arXiv preprint arXiv:2304.01294, 2023.
  * [15] Y. T. Chow, S. W. Fung, S. Liu, L. Nurbekyan, and S. Osher.  A numerical algorithm for inverse problem from partial boundary measurement arising from mean field game problem.  Inverse Problems, 39(1):014001, 2022.
  * [16] M. Cirant and A. Porretta.  Long time behavior and turnpike solutions in mildly non-monotone mean field games.  ESAIM. COCV, 27, 2021.
  * [17] L. Ding, W. Li, S. Osher, and W. Yin.  A mean field game inverse problem.  Journal of Scientific Computing, 92(1):7, 2022.
  * [18] M Ding, H. Liu, and G. Zheng.  Determining a stationary mean field game system from full/partial boundary measurement.  arXiv preprint arXiv:2308.06688, 2023.
  * [19] D. Evangelista, R. Ferreira, D. Gomes, L. Nurbekyan, and V. Voskanyan.  First-order, stationary mean-field games with congestion.  Nonlinear Analysis, 173:37–74, 2018.
  * [20] H. Gao, A. Lin, R. A Banez, W. Li, Z. Han, S. J. Osher, and H. V. Poor.  Belief and opinion evolution in social networks: A high-dimensional mean field game approach.  In ICC 2021-IEEE International Conference on Communications, pages 1–6. IEEE, 2021.
  * [21] D. Gomes, L. Nurbekyan, and E. Pimentel.  Economic models and mean-field games theory.  Publicaoes Matematicas, IMPA, Rio, Brazil, 2015.
  * [22] D. Gomes and J. Saúde.  A mean-field game approach to price formation.  Dynamic Games and Applications, 11(1):29–53, 2021.
  * [23] D. Gomes and X. Yang.  The Hessian Riemannian flow and Newton’s method for effective Hamiltonians and Mather measures.  ESAIM: Mathematical Modelling and Numerical Analysis, 54(6):1883–1915, 2020.
  * [24] O. Guéant, J.-M. Lasry, and P.-L. Lions.  Mean field games and applications.  In Paris-Princeton lectures on mathematical finance 2010, pages 205–266. Springer, 2011.
  * [25] B. Hamzi, H. Owhadi, and Y. Kevrekidis.  Learning dynamical systems from data: A simple cross-validation perspective, part iv: case with partial observations.  Physica D: Nonlinear Phenomena, 454:133853, 2023.
  * [26] M. Huang, P. E. Caines, and R. P. Malhamé.  An invariance principle in large population stochastic dynamic games.  Journal of Systems Science and Complexity, 20(2):162–172, 2007\.
  * [27] M. Huang, P. E. Caines, and R. P. Malhamé.  Large-population cost-coupled lqg problems with nonuniform agents: individual-mass behavior and decentralized $\epsilon$-nash equilibria.  IEEE transactions on automatic control, 52(9):1560–1571, 2007.
  * [28] 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 2007 46th IEEE Conference on Decision and Control, pages 121–126. IEEE, 2007.
  * [29] 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 and Systems, 6(3):221–252, 2006.
  * [30] O. Imanuvilov, H. Liu, and M. Yamamoto.  Lipschitz stability for determination of states and inverse source problem for the mean field game equations.  arXiv preprint arXiv:2304.06673, 2023.
  * [31] M. V. Klibanov, J. Li, and H. Liu.  Hölder stability and uniqueness for the mean field games system via carleman estimates.  arXiv preprint arXiv:2304.00646, 2023.
  * [32] M. V. Klibanov, J. Li, and Z. Yang.  Convexification for a coefficient inverse problem of mean field games.  arXiv preprint arXiv:2310.08878, 2023.
  * [33] J.-M. Lasry and P.-L. Lions.  Jeux à champ moyen. I–le cas stationnaire.  Comptes Rendus Mathématique, 343(9):619–625, 2006.
  * [34] J.-M. Lasry and P.-L. Lions.  Jeux à champ moyen. II–horizon fini et contrôle optimal.  Comptes Rendus Mathématique, 343(10):679–684, 2006.
  * [35] J.-M. Lasry and P.-L. Lions.  Mean field games.  Japanese journal of mathematics, 2(1):229–260, 2007.
  * [36] M. Laurière.  Numerical methods for mean field games and mean field type control.  arXiv preprint arXiv:2106.06231, 2021.
  * [37] W. Lee, S. Liu, W. Li, and S. Osher.  Mean field control problems for vaccine distribution.  arXiv preprint arXiv:2104.11887, 2021.
  * [38] W. Lee, S. Liu, H. Tembine, W. Li, and S. J. Osher.  Controlling propagation of epidemics via mean-field control.  SIAM Journal on Applied Math, 2020.
  * [39] A. T. Lin, S. W. Fung, W. Li, L. Nurbekyan, and S. J. Osher.  Alternating the population and control neural networks to solve high-dimensional stochastic mean-field games.  Proceedings of the National Academy of Sciences, 118(31), 2021.
  * [40] P.-L. Lions.  Cours au Collège de France, 2007-2013.
  * [41] H. Liu, C. Mou, and S. Zhang.  Inverse problems for mean field games.  Inverse Problems, 2023.
  * [42] H. Liu and S. Zhang.  Simultaneously recovering running cost and hamiltonian in mean field games system.  arXiv preprint arXiv:2303.13096, 2023.
  * [43] S. Liu, M. Jacobs, W. Li, L. Nurbekyan, and S. J. Osher.  Computational methods for first-order nonlocal mean field games with applications.  SIAM Journal on Numerical Analysis, 59(5):2639–2668, 2021.
  * [44] S. Liu and L. Nurbekyan.  Splitting methods for a class of non-potential mean field games.  Journal of Dynamics and Games, 2020, 2020.
  * [45] R. Meng and X. Yang.  Sparse Gaussian processes for solving nonlinear PDEs.  Journal of Computational Physics, 490:112340, 2023.
  * [46] C. A. Micchelli and M. Pontil.  On learning vector-valued functions.  Neural computation, 17(1):177–204, 2005.
  * [47] C. Mou, X. Yang, and C. Zhou.  Numerical methods for mean field games based on Gaussian processes and Fourier features.  Journal of Computational Physics, 2022.
  * [48] L. Nurbekyan and J. Saúde.  Fourier approximation methods for first-order nonlocal mean-field games.  Portugaliae Mathematica, 75(3):367–396, 2019.
  * [49] H. Owhadi and C. Scovel.  Operator-Adapted Wavelets, Fast Solvers, and Numerical Homogenization: From a Game Theoretic Approach to Numerical Approximation and Algorithm Design, volume 35.  Cambridge University Press, 2019.
  * [50] M. Raissi, P. Perdikaris, and G. E. Karniadakis.  Machine learning of linear differential equations using Gaussian processes.  Journal of Computational Physics, 348:683–693, 2017.
  * [51] M. Raissi, P. Perdikaris, and G. E. Karniadakis.  Numerical Gaussian processes for time-dependent and nonlinear partial differential equations.  SIAM Journal on Scientific Computing, 40(1):A172–A198, 2018.
  * [52] L. Ruthotto, S. J. Osher, W. Li, L. Nurbekyan, and S. W. Fung.  A machine learning framework for solving high-dimensional mean field game and mean field control problems.  Proceedings of the National Academy of Sciences, 117(17):9183–9193, 2020.
  * [53] C. K. Williams and C. E. Rasmussen.  Gaussian processes for machine learning.  MIT press Cambridge, MA., 2006.
  * [54] L. Yang, S. Liu, T. Meng, and S. J. Osher.  In-context operator learning with data prompts for differential equation problems.  Proceedings of the National Academy of Sciences, 120(39):e2310142120, 2023.
  * [55] L. Yang, X. Sun, B. Hamzi, H. Owhadi, and N. Xie.  Learning dynamical systems from data: A simple cross-validation perspective, part v: Sparse kernel flows for 132 chaotic dynamical systems.  arXiv preprint arXiv:2301.10321, 2023.
  * [56] X. Yang and H. Owhadi.  A mini-batch method for solving nonlinear PDEs with Gaussian processes.  arXiv preprint arXiv:2306.00307, 2023.