Deep Backward and Galerkin Methods for the Finite State Master Equation

In 1950, John Nash introduced the concept of a Nash equilibrium, an idea that describes when players trying to minimize an individual cost have no incentive to change their strategy [45]. Nash proved that, under mild conditions, a Nash equilibrium always exists in mixed strategies. The question remained however, even when a Nash equilibrium was known to exist, whether there is an algorithmic way to solve for it. It turns out that, in general, computing a Nash equilibrium is a very difficult task [17].

The origins of limiting models for many-agent games can be traced back to [3] and to [51], both addressing one-shot games. In the realm of stochastic dynamical games, mean field game (MFG) theory was pioneered by [40, 41], along with an independent development by [36]. MFG paradigms study situations in which a very large number of agents interact in a strategic manner. Specifically, MFGs provide one way to approximate Nash equilibria of large, anonymous, symmetric finite-player games.

In characterizing an equilibrium in MFGs, we introduce a representative player who, instead of reacting to $n$ players, responds to a flow of measures. The MFG equilibrium is consequently defined by a fixed point, ensuring that the distribution of the dynamics of the representative player under optimality aligns with the flow of measures to which they respond. It is worth noting that Nash equilibria in finite-player games are similarly described using fixed points.

Returning to the motivating question for $n$-player games, an MFG equilibrium provides an approximate Nash equilibrium in a corresponding finite-player game: if all the $n$ players use the strategy determined by the MFG, each player can be at most $\varepsilon_{n}$ better off by choosing another strategy, and $\varepsilon_{n}$ converges to zero as $n\to\infty$. This approximation is valid for games in which players are homogeneous and the interactions are symmetric, in the sense that they occur through the players’ empirical distribution. While static games are prevalent in game theory for the sake of simplicity, dynamic games are more realistic for applications. In this work, we consider the MFG studied by [5, 15, 30] which approximates an $n$-player stochastic differential game with a finite number of states.

Classically, an MFG equilibrium can be described by a coupled forward-backward system of equations, called the MFG system. The backward equation is a Bellman-type equation describing the evolution of a representative player’s value function (the optimal cost), from which the optimal strategy can be deduced. The forward equation describes the evolution of the population’s distribution, which coincides with the distribution of one representative player over the state space. We refer to the books by [7] and [12] for more details on this approach using partial differential equations (PDEs) or stochastic differential equations (SDEs). In finite-state MFGs as we will consider in this work, the forward-backward system is a system of ordinary differential equations (ODEs); the interested reader is referred to [12, Chapter 7.2].

The solution to the forward-backward system depends on the initial distribution of the population. The representative player’s value function and control depend only implicitly on the population distribution. This dependence can be expressed explicitly through the notion of a decoupling field. The MFG equilibrium can be described by the so-called master equation, introduced by [43]. The interested reader is referred to the monograph by [11] for a detailed analysis of the master equation and the question of convergence of finite-player Nash equilibria to MFG equilibria. The solution to this equation is a function of the time, the representative player’s state, and the population distribution, which will be denoted respectively by $t$, $x$ and $\eta$. The master equation derives from the MFG system (see Proposition 1 below), whose structure is crucial for our analysis. Besides fully characterizing the MFG equilibrium for any initial distribution, the master equation is also used to study the convergence of the finite player equilibrium to its MFG counterpart [11, 19, 20]. In this paper, we are interested in the finite state, finite horizon master equation studied by [5, 15]:

Above, $[d]:=\{1,\dots,d\}$ is the state space, $T$ is the terminal time of the game, $\bar{H}$ is the Hamiltonian (defined in Section 2), and $g$ is a terminal condition. The function $\gamma^{*}$ denotes the minimizer of the Hamiltonian, which under regularity assumptions is equal to the gradient of the Hamiltonian. The variables are the time $t$, the state $x$, and the distribution $\eta$, which lies in the $(d-1)$-dimensional simplex ${\mathcal{P}}([d])$. We use the notation $\Delta_{x}U(t,\cdot,\eta):=(U(t,y,\eta)-U(t,x,\eta))_{y\in[d]}$ for each $x\in[d]$. Derivatives on the simplex, denoted by $D^{\eta}_{yz}U(t,x,\eta)$, can be thought of as usual directional derivatives and will be explicitly defined in Section 2.

Note that (1) is, in general, a nonlinear PDE, for which there is no analytical solution and the question of well-posedness is challenging (for an overview, see the book by [22]). In fact, even computing numerically an approximate solution is a daunting problem. When $d\leq 3$, traditional methods such as finite difference schemes can be employed. However, the computational cost of these methods becomes prohibitive for dimensions greater than $3$. For example, in grid-based methods, the number of points grows exponentially quickly in $d$. Other methods include Markov chain approximations, see [4]. We refer to [38] for an overview of classical methods for PDEs and [1, 42] for an overview of classical methods for MFGs. Besides, many applications require $d>3$. And additionally, [8] recently demonstrated that MFGs with a continuous state space can be approximated by MFGs with finite, but large, state spaces. This approximation relies on the convergence of the solution of the associated master equations as the number of states tends to infinity. Neural networks provide one avenue to mitigate the curse of dimensionality, and this has triggered interest for machine learning-based methods, as we discuss in the next section.

The value function under the unique Nash equilibrium in the $n$-player game converges to the solution to the master equation up to an error of order ${\mathcal{O}}(n^{-1})$ (and so do the optimal controls) [5, 15]. Then, we prove that the two algorithms we propose, the deep backward master equation (DBME) and deep Galerkin master equation (DGME) methods, approximate the master equation which, as discussed in the prior paragraph, is not easy to solve analytically. So while the primary object of interest for this paper is the master equation, our main contribution completes the picture for how to provably approximate the solution to the $n$-player game.

In this work, we propose and rigorously analyze deep learning methods to tackle the master equation (1) for finite-state MFGs beyond what is possible with traditional numerical methods. Our numerical methods build upon the recent development of deep learning approaches for PDEs. In the past few years, numerous methods have been proposed along these lines, such as the deep BSDE method [21], the deep Galerkin method [53], physics informed neural networks [49] and deep backward dynamic programming [48]. In the context of MFGs and mean field control problems, deep learning methods have been applied using the probabilistic approach [14, 24, 28] or the analytical approach [2, 10, 13, 50]. The aforementioned works consider neural networks whose inputs are only the individual state and possibly time. A few works have considered neural networks taking as inputs a representation of the mean field. The paper by [27] proposed a deep backward method for mean field control problems. As with the first method (called DBME) we present below, it is based on a discrete-time, backward scheme. The rate of convergence in terms of number of particles was studied in [29]. However, this method solves the Bellman equation for mean field control problems while we solve the master equation for MFGs, which fails to satisfy a comparison principle. Furthermore, [27] considered a continuous space model using particle-based approximations while we focus on a finite space model. Work by [18] solves continuous space mean field control problems using population-dependent controls by directly learning the control function instead of using a backward dynamic programming equation. As for MFG with population-dependent controls, [47] proposed a deep reinforcement learning method for finite-state MFGs based on a fictitious play and deep Q-networks. In contrast, our two methods do not rely on reinforcement learning and we provide a detailed analysis of the algorithms’ errors. We refer to the overview of [35] for a review of deep learning methods in optimal control and games.

In this work, we will focus on two types of methods. The first approach finds its roots in the pioneering work of  [21], which exploits the connection between SDEs and PDEs: using Feynman–Kac type formulas, the PDE solution satisfies a backward SDE, that is, they satisfy a terminal condition. Corresponding to the SDE, the PDE has a terminal condition. Such SDEs and PDEs crop up often in mathematical finance, most iconically in the Black–Scholes equation for pricing European options [9]. While the original method of [21] consists in replacing the BSDE by a shooting method for a forward SDE, subsequent works exploit the backward structure differently. In particular, deep backward methods were explored in control problems in the nonlinear and fully-nonlinear cases [37, 48]. The Deep Backward Dynamic Programming algorithm (DBDP) from [37] (more specifically their DBDP2 algorithm) uses a sequence of neural networks on a discretization of $[0,T]$. Training is done by starting from the terminal time and, going backward in time, the neural network for a given time step is trained to minimize a loss function which involves the neural network at the following time step.

The second method we propose is the DGME, introduced as the deep Galerkin method by [53]. The PDE solution is approximated by a neural network, which is a function of the time and space variables. Training is done by minimizing an empirical loss which captures the residual of the PDE, as well as possible boundary conditions. The loss is computed over points sampled inside the domain and, if needed, on its boundary. In this way, the method does not require any discretization of space or time. The solution is learned over the whole domain thanks to the generalization capability of neural networks.

For these methods, most theoretical results rely on universal approximation theorems for suitable classes of functions related to the regularity of PDE solutions. In that sense, other classes of approximators could be used instead of neural networks. However, neural networks seem to be a very suitable choice to solve PDEs. Two of the main limitations of neural networks are the lack of explainability and the large amount of data required for training, but these are not issues when solving PDEs with the aforementioned methods. Indeed, explainability is not a priority when approximating PDE solutions, and data points are obtained by sampling over the PDE domain as many times as desired.

While the well-posedness of (1) was studied by [5] and [15], no analytical solution is known. Therefore we provide two algorithms, we prove their correctness, and numerically solve the master equation (1) for some examples.

We begin by defining notation in Section 2 before providing a description of the MFG model and recalling the major results concerning the MFG and the master equation in Section 3. In Section 4 we outline the architecture of the neural networks we consider. Then in Sections 5 and 6, we describe and analyze the DBME and DGME, respectively. Finally, in Section 7, we present numerical results for two examples using these algorithms.

We consider a representative player that uses Markovian controls taking values in a set of rates $\mathbb{A}\subseteq\mathbb{R}_{+}:=[0,+\infty)$. We will consider rates bounded in a compact set; that is, $\mathbb{A}:=[\mathfrak{a}_{l},\mathfrak{a}_{u}]$, with $0\leq\mathfrak{a}_{l}\leq\mathfrak{a}_{u}<+\infty$. A Markovian control is a measurable function $\alpha:\mathbb{R}_{+}\times[d]\to\mathbb{A}^{d}_{[d]}$ where $\mathbb{A}^{d}_{[d]}=\cup_{x\in[d]}\mathbb{A}^{d}_{-x}$, and:

Here, for any $y\neq x$, $\alpha_{y}(t,x):=\alpha(t,x)_{y}\in\mathbb{A}$ is the rate of transition at time $t$ to move from state $x$ to state $y$ and as usual with continuous-time Markov chains, $\alpha_{x}(t,x)=-\sum_{y,y\neq x}\alpha_{y}(t,x)$. We denote by ${\mathcal{Q}}[\mathbb{A}]$ the set of $d\times d$ transition-rate matrices with rates in $\mathbb{A}$. We denote by ${\mathcal{A}}$ the set of measurable mappings $\alpha:[0,T]\to{\mathcal{Q}}[\mathbb{A}]$.

Fix $t_{0}\in[0,T]$ to be the initial time. The jump dynamics of the representative player’s state ${\mathcal{X}}$ are given by:

where ${\mathcal{N}}$ is a Poisson random measure with intensity measure $\nu$ given by

with $\mathbb{A}^{d}_{y}:=\{a\in\mathbb{A}^{d}\mid a_{x}=0\text{ for all }x\neq y\}$, where Leb is the Lebesgue measure on $\mathbb{R}$, and $E$ is any Borel measurable set. We can think of $\mathbb{A}^{d}_{y}$ as the $y$-projection of $\mathbb{A}^{d}$, the set of rates for each state transition; thus, $\mathbb{A}^{d}_{y}$ is the set of transition rates from state $y\in[d]$.

In order to set up the MFG equilibrium, we describe the cost structure of the game with respect to a representative player, who selects a Markovian control $\alpha$ in order to play against a smooth flow of measures $\mu\in{\mathcal{C}}^{1}([t_{0},T],{\mathcal{P}}([d]))$. We emphasize that the representative player’s cost includes the flow of measures $\mu$. In this work, we make use of the regularity result from [15], which assumes that the total running cost is separable in the form $f+F$. The running cost is $f:[d]\times\mathbb{A}^{d}_{[d]}\to\mathbb{R}$, a function of the current state and control, while the mean field cost is $F:[d]\times{\mathcal{P}}([d])\to\mathbb{R}$, a function of the state and mean field distribution. Recall that $\alpha_{x}(t,x)$ does not represent a rate, thus $f(x,a)$ is assumed independent of $a_{x}$ as a technicality. Given a flow $\mu\in{\mathcal{C}}^{1}([t_{0},T],{\mathcal{P}}([d]))$ of distributions, the representative agent’s control problem is to minimize:

over all $\alpha\in{\mathcal{A}}$ where $({\mathcal{X}}_{t_{0}},\mu(t_{0}))=(x,\eta)\in[d]\times{\mathcal{P}}([d])$. Given the flow of measure $\mu$ and an initial state ${\mathcal{X}}_{t_{0}}=x\in[d]$, the distribution (or law) of the representative player, denoted $({\mathcal{L}}({\mathcal{X}}_{t}))_{t\in[t_{0},T]}$, is a fully deterministic flow of measures on ${\mathcal{P}}([d])$. Formally, when $({\mathcal{L}}({\mathcal{X}}_{t}))_{t\in[t_{0},T]}=(\mu(t))_{t\in[t_{0},T]}$ and $\alpha$ minimizes the cost:

we say that we have a mean field game equilibrium and we refer to this optimal cost as the value function. To be mathematically precise, the MFG equilibrium is a fixed point, as described next.

Crucially, notice that the equilibrium is defined for a fixed initial distribution $\mu_{0}$. For different initial distributions, we generally obtain different equilibrium controls and mean field flows. The master equation characterizes MFG equilibria for any initial distribution.

As explained below in Section 3.3, the equilibrium can be characterized using optimality conditions which involve the Hamiltonian, defined as:

where,

At times, using $\bar{H}$ rather than $H+F$ will simplify the presentation of the PDE.

In this section, we list the assumptions made on our MFG model. Such assumptions are standard and appear in previous work by [5, 15, 16].

Recall that the action set is $\mathbb{A}:=[\mathfrak{a}_{l},\mathfrak{a}_{u}]$ where $0\leq\mathfrak{a}_{l}\leq\mathfrak{a}_{u}<\infty$.¹¹1In the papers by [5, 16] it is assumed that $\mathfrak{a}_{l}>0$; however in his thesis, Cecchin managed to relax the assumption to allow for $\mathfrak{a}_{l}=0$ [15].

We note that since $\gamma^{*}$ is given by (7), one sufficient condition to guarantee this is when $f$ is strictly convex in $a$.

As for the mean field cost $F$, and the terminal cost $g$, we assume

Intuitively, the Lasry–Lions monotonicity condition means that a representative player would prefer to avoid congestion to decrease their cost.

In the following assumption, we specify the regularity of $H$, which, since it is defined through $f$ implicitly defines most of the necessary regularity for $f$. Moreover, since $F$ is defined on a compact set, its regularity assumptions imply it is bounded.

Note that when $H$ is differentiable, [31, Proposition 1] proved:

In order to keep this level of generality for the regularity of $H$ (namely, that these additional assumptions hold on $[-W,W]$ but not necessarily globally), we will on occasion prove that a particular $p$-argument of $H$ is uniformly bounded by $W$. For the purpose of this paper, the assumption on $W$ in (9) is not too restrictive. Computations with the Hamiltonian only appear with the argument $\Delta_{x}U$ or $\Delta_{x}{\mathcal{U}}$, where ${\mathcal{U}}$ is a neural network specified in the following sections that approximates $U$. As we explain in detail in Remark 5, both of these arguments are bounded by $W$ from (9) and hence the Hamiltonian will always be regular enough when it is required.

Throughout the paper $C$ denotes a generic, positive constant that depends only on the problem’s parameters (that is, the parameters introduced above). To alleviate excessive notation, $C$’s value may change from one line to another.

Next up, we recall some established results from the study of finite-state, finite-horizon MFGs and master equations.

The functions $u^{t_{0},\eta}$ and $\mu^{t_{0},\eta}$ are defined in a moment but we first describe their meanings. The measure $\mu^{t_{0},\eta}$ is the evolving MFG equilibrium as in Definition 1. The value function $u^{t_{0},\eta}$ is the value function (optimal cost) of the MFG starting at time $t_{0}$ with initial distribution $\eta\in{\mathcal{P}}([d])$ along the MFG equilibrium. So, the value $u^{t_{,}\eta}(t,x)$ is the remaining optimal cost to a player in state $x$ at time $t\geq t_{0}$ until the game ends, at time $T$.

The MFG system is:

Note that the system is composed of a forward equation and a backward equation. The above ODE system characterizes the MFG equilibrium for a given initial condition $(t_{0},\eta)$. However, our goal is to compute the MFG equilibrium for every initial condition. This would require solving a continuum of ODE systems, which is not feasible. Instead, we will focus on the master equation. Informally, it encompasses the continuum of ODE systems and its solution captures the dependence of the value function on the mean field. More precisely, the master equation’s solution $U$ is defined through its MFG system. The following proposition is a combination of the results from [15, Proposition 1, Proposition 5, Theorem 6]²²2In fact, we reuse the argument from Proposition 1, but applied to the MFG system. and [11, Section 1.2.4].

First, we introduce some notation relevant to the algorithm. Then, we justify some assumptions made on the neural networks under consideration, and introduce an auxiliary Kolmogorov forward ODE, before introducing the DBME algorithm itself.

The time interval $[0,T]$ will be partitioned into $\pi:=\{t_{0}=0,t_{1},\dots,t_{N-1},t_{N}=T\}$, with increments denoted $\Delta t_{i}:=t_{i+1}-t_{i}$. Define

Throughout we assume that there exist $\check{c},\hat{c}>0$ such that:

that is, $|\pi|^{-1}$ is roughly proportional to $N$. For functions and random variables at time $t_{i}\in\pi$, we will often use only a subscript $i$ as opposed to $t_{i}$ for brevity.

For the time $t_{i}\in\pi$, we approximate the master equation’s solution $U(t_{i},x,\eta)$ by ${\mathcal{U}}_{i}\in{\mathcal{N}}{\mathcal{N}}$ with input dimension $\delta_{1}=d+2$. Its parameters are denoted by $\theta^{i}$ and we use the notation ${\mathcal{U}}_{i}(x,\eta;\theta^{i})$. We will sometimes use the shorthand notation $\hat{\mathcal{U}}_{i}(x,\eta):={\mathcal{U}}_{i}(x,\eta;\hat{\theta}^{i})$, where the vector of parameters $\hat{\theta}^{i}$ optimizes a loss function defined in the sequel. We will not use a neural network at the last time in the partition and simply set $\hat{\mathcal{U}}_{N}(x,\eta):=g(x,\eta)$.

Let $L=2$ for instance and $\bar{\delta}=(\delta_{2}(d+3))$. For the vector of parameters $\theta\in\mathbb{R}^{\bar{\delta}}$, denote the parameters from the scaling inside the activation function by $\check{\theta}$. We can write a neural network ${\mathcal{U}}_{i}(x,\eta;\theta)\in{\mathcal{N}}{\mathcal{N}}$ as:

where $a^{k}:=(a^{y,k})_{y\in[d]}\in\mathbb{R}^{d}$, $a^{0,k}\in\mathbb{R}$, $b\in\mathbb{R}^{\delta_{2}}$, $c^{k}\in\mathbb{R}$ for all $k=1,\dots,\delta_{2}$, and these parameters make up the entries in $\theta\in\mathbb{R}^{\bar{\delta}}$ and $\check{\theta}:=(a^{k})_{k=1}^{\delta_{2}}$. In the case $L=2$, we note that $\check{\theta}$ is a $(\delta_{2}\times 2)$-matrix with real entries and denote $\|\cdot\|_{2}$ as the matrix operator $2$-norm. Per [32, Section 4], we know that when $\|\check{\theta}\|_{2}<K$, some positive real number, the Lipschitz constant of the network with parameters $\theta\in\mathbb{R}^{\bar{\delta}}$ is no more than $K$ (possibly after modifying $K$ to account for the Lipschitz constant of the activation function $\phi$). So, for the rest of this section, which focuses on the DBME algorithm, we define the following set of parameters:

where recall that $C_{L,U}$ is the Lipschitz constant of $U$ depending only on the problem parameters that can be computed from [15, Section 3]. Thus any network ${\mathcal{U}}(\cdot,\cdot;\theta)$ where $\theta\in\Theta$ is Lipschitz with constant no more than $2C_{L,U}$.

One final element to introduce before presenting the algorithm is the solution of a forward Kolmogorov’s equation. This equation mirrors the flow of measures $\mu$ in the MFG equilibrium. However, instead of utilizing the master equation’s solution $U$ in the dynamics, we employ a neural network. This solution is a key element in the algorithm’s loss function. In details: for any $\kappa\in{\mathcal{P}}([d])$, any $i\in\{0,\dots,N-1\}$, and any $\theta\in\Theta$, define $M^{\theta}_{i}(\kappa):=\rho_{i}^{\theta,\kappa}(t_{i+1})$, where $\rho_{i}^{\theta}(t):=(\rho_{i}^{\theta}(t,x))_{x\in[d]}$ uniquely solves the (forward) ODE:

The previous display is Kolmogorov’s forward equation over the time interval $[t_{i},t_{i+1}]$ with initial distribution $\kappa\in{\mathcal{P}}([d])$, and which corresponds to a non-linear Markov chain with rate matrix $\gamma^{*}_{x}(y,\Delta_{y}{\mathcal{U}}_{i}(\cdot,\rho^{\theta,\kappa}_{i}(t);\theta^{i}))$. Note that the control’s input is the difference $\Delta_{y}{\mathcal{U}}_{i}$ of the neural network ${\mathcal{U}}_{i}$, rather than that of $U$; in light of (14), this is similar to the MFG system (LABEL:eqn:mfg). Moreover, $\rho^{\theta,\kappa}_{i}$ appears within the rate matrix, which makes the Markov chain non-linear. The equation (LABEL:eqn:M_theta_i_def) is well-posed because the neural network and $\gamma^{*}$ are Lipschitz. The notation $M^{\theta}_{i}(\kappa)$ can be thought of as the terminal distribution of this non-linear Markov chain at time $t_{i+1}$; that is, the effect of propagating the initial distribution $\kappa$ by a control using neural network values, parameterized by $\theta\in\Theta$.

The next remark motivates the algorithm’s structure and the form of the loss function (23) in the algorithm.