Actor critic learning algorithms for mean-field control with moment
  neural networks
We develop a new policy gradient and actor-critic algorithm for solving
mean-field control problems within a continuous time reinforcement learning
setting. Our approach leverages a gradient-based representation of the value
function, employing parametrized randomized policies. The learning for both the
actor (policy) and critic (value function) is facilitated by a class of moment
neural network functions on the Wasserstein space of probability measures, and
the key feature is to sample directly trajectories of distributions. A central
challenge addressed in this study pertains to the computational treatment of an
operator specific to the mean-field framework. To illustrate the effectiveness
of our methods, we provide a comprehensive set of numerical results. These
encompass diverse examples, including multi-dimensional settings and nonlinear
quadratic mean-field control problems with controlled volatility.
Huyên Pham, Xavier Warin
# Actor critic learning algorithms for mean-field control

###### Abstract

We develop a new policy gradient and actor-critic algorithm for solving mean-field control problems within a continuous time reinforcement learning setting. Our approach leverages a gradient-based representation of the value function, employing parametrized randomized policies. The learning for both the actor (policy) and critic (value function) is facilitated by a class of moment neural network functions on the Wasserstein space of probability measures, and the key feature is to sample directly trajectories of distributions. A central challenge addressed in this study pertains to the computational treatment of an operator specific to the mean-field framework. To illustrate the effectiveness of our methods, we provide a comprehensive set of numerical results. These encompass diverse examples, including multi-dimensional settings and nonlinear quadratic mean-field control problems with controlled volatility.

**Keywords:** Mean-field control, reinforcement learning, policy gradient, moment neural network, actor-critic algorithms.

## 1 Introduction

This paper is concerned with the numerical resolution of mean-field (a.k.a. McKean Vlasov) control in continuous time in a partially model-free reinforcement learning setting. The dynamics of the controlled mean field stochastic differential equation on \(\mathbb{R}^{d}\) is in the form

\[\mathrm{d}X_{t}=\;b(t,X_{t},\mathbb{P}_{X_{t}},\alpha_{t})\mathrm{d}t+ \sigma(t,X_{t},\mathbb{P}_{X_{t}},\alpha_{t})\mathrm{d}W_{t},\quad 0\leq t \leq T,\;X_{0}\sim\mu_{0},\]

where \(W\) is a standard \(d\)-dimensional Brownian motion on a filtered probability space \((\Omega,\mathcal{F},\mathbb{F}=(\mathcal{F}_{t})_{t},\mathbb{P})\), \(\mu_{0}\in\mathcal{P}_{2}(\mathbb{R}^{d})\), the Wasserstein space of square integrable probability measures, \(\mathbb{P}_{X_{t}}\) denotes the marginal distribution of \(X_{t}\) at time \(t\), and the control process \(\alpha\) is valued in \(A\subset\mathbb{R}^{p}\). The coefficients \(b,\sigma\) are in the separable form:

\[\mathbf{(SC)}\qquad\qquad b(t,x,\mu,a)\;=\;\beta(t,x,\mu)+C(t,a),\qquad \sigma\sigma^{\tau}(t,x,\mu,a)\;=\;\Sigma(t,x,\mu)+\vartheta(t,a),\]

for \((t,x,\mu,a)\in[0,T]\times\mathbb{R}^{d}\times\mathcal{P}_{2}(\mathbb{R}^{d}) \times\mathbb{R}^{p}\), where \(\beta\), \(\Sigma\) depending on the state variable and its probability distribution are unknown functions, while the coefficients \(C\), \(\vartheta\) on the control are known functions on \([0,T]\times\mathbb{R}^{p}\).

The expected total cost associated to a control \(\alpha\) is given by

\[\mathbb{E}\Big{[}\int_{0}^{T}f(X_{t},\mathbb{P}_{X_{t}},\alpha_{t})\mathrm{d }t+g(X_{T},\mathbb{P}_{X_{T}})\Big{]}, \tag{1.1}\]
and we denote by \((t,x,\mu)\mapsto V(t,x,\mu)\) the associated value function defined on \([0,T]\times\mathbb{R}^{d}\times\mathcal{P}_{2}(\mathbb{R}^{d})\). The analytical forms of \(f\), \(g\) are unknown, but it is assumed that given an input \((x,\mu,a)\in\mathbb{R}^{d}\times\mathcal{P}_{2}(\mathbb{R}^{d})\times A\), we obtain (by observation or from a blackbox) the realized output costs \(f(x,\mu,a)\) and \(g(x,\mu)\). Such setting is partially model-free, as the model coefficients \(\beta\), \(\Sigma\), \(f\), and \(g\) are unknown, and we only assume that the action functions \(C\), \(\vartheta\) on the control are known.

The theory and applications of mean-field control (MFC) problems have generated a vast literature in the last decade, and we refer to the monographs [2], [3] for a comprehensive treatment of this topic. From a numerical aspect, the main challenging issue is the infinite dimensional feature of MFC coming from the distribution law state variable. In a model-based setting, i.e., when all the coefficients \(b\), \(\sigma\), \(f\) and \(g\) are known, several recent works have proposed deep learning schemes for MFC, based on neural network approximations of the feedback control and/or the value function solution to the Hamilton-Jacobi-Bellman equation from the dynamic programming or backward stochastic differential equations (BSDEs) from the maximum principle, see [5], [8], [9], [11], [17], [16], [14].

The approximation of solutions to MFC in a (partially) model-free setting is the purpose of reinforcement learning (RL) where one learns in an unknown environment the optimal control (and the value function) by repeatedly trying policies, observing state, receiving and evaluating rewards, and improving policies. RL is a very active branch of machine learning, see the seminal reference monograph [18], and has recently attracted attention in the context of mean-field control in discrete-time, and mostly by \(Q\)-learning methods, see [6], [1], [10].

In this paper, we consider a partially model-free continuous time setting as described above, and adopt a policy gradient approach as in [7]. This relies on a gradient representation of the cost functional associated to randomized policies, which makes appear an additional operator term \(\mathcal{H}\) compared to the classical diffusion setting of [12], and specific to the mean-field setting. The computational treatment of this operator \(\mathcal{H}\) on functions defined on the Wasserstein space is the crucial issue, and has been handled in [7] only for one-dimensional linear quadratic (LQ) models, hence with a very particular dependence of the value function and optimal control on the law of the state. Here, we address the general dependence of the coefficients on the distribution state, and deal with the operator \(\mathcal{H}\) by means of the class of moment neural networks. This class of neural networks consists of functions that depend on the measure via its first \(L\) moments, and satisfies some universal approximation theorem for functions defined on the Wasserstein space, see [15], [20]. We then design an actor-critic algorithm for learning alternately the optimal policy (actor) and value function (critic) with moment neural networks, which provides an effective resolution of MFC control problems in a (partially) model-free setting beyond the LQ setting for multivariate dynamics with control on the drift and the volatility. Our actor-critic algorithm has the structure of general actor-critic algorithms but during gradient iterations, instead of following a single state trajectory by sampling, we follow the evolution of an entire distribution that is initially randomly chosen and described by its empirical measure obtained with a large fixed number of particles. Then, the batch version of the algorithm consists in sampling and following \(m\) distributions together to estimate the gradients.

The outline of the paper is organized as follows. We recall in Section 2 the gradient representation of the functional cost with randomized policies and formulate notably the expression of the operator \(\mathcal{H}\). In Section 3, we consider the class of moment neural networks, and show how it acts on the operator \(\mathcal{H}\). We present in Section 4 the actor-critic algorithm, and Section 5 is devoted to numerical results for illustrating the accuracy and efficiency of our algorithm. We present various examples with control on the drift and on the volatility, and non LQ examples in a multi-dimensional setting.

**Notations.** The scalar product between two vectors \(x\) and \(y\) is denoted by \(x\cdot y\), and \(|\cdot|\) is the Euclidian norm. Let \(Q=(Q_{i_{1}\ldots i_{q}})\in\mathbb{R}^{d_{1}\times\ldots\times d_{q}}\) be a tensor of order \(q\), and \(P=(P_{i_{1}\ldots i_{p}})\in\mathbb{R}^{d_{1}\times\ldots\times d_{p}}\) be a tensor of order \(p\leq q\). We denote by \(Q\circ P\) the circ product defined as the tensor in \(\mathbb{R}^{d_{p+1}\times\ldots\times d_{q}}\) with components:

\[[Q\circ P]_{i_{p+1}\ldots i_{q}}=\ \sum_{i_{1},\ldots,i_{p}}Q_{i_{1}\ldots i_{ p}i_{p+1}\ldots i_{q}}P_{i_{1}\ldots i_{p}}.\]

When \(q=p=1\), \(\circ\) is the scalar product in \(\mathbb{R}^{d_{1}}\). When \(q=2\), \(p=1\), \(Q\circ P=Q^{\intercal}P\in\mathbb{R}^{d_{2}}\) where \({}^{\intercal}\) is the transpose matrice operator. When \(q=p=2\), \(\circ\) is the inner product \(Q\circ P=\operatorname{tr}(Q^{\intercal}P)\) where \(\operatorname{tr}\) is the trace operator. When \(q=3\), \(Q\circ P\) is a vector in \(\mathbb{R}^{d_{3}}\) for \(p=2\), and a matrix in \(\mathbb{R}^{d_{2}\times d_{3}}\) for \(p=1\).
Preliminaries

We adopt a policy gradient approach by searching optimal control among parametrized randomized policies, i.e., family of probability transition kernels \(\pi_{\theta}\) from \([0,T]\times\mathbb{R}^{d}\times\mathcal{P}_{2}(\mathbb{R}^{d})\) into \(\mathbb{R}^{p}\), with densities \(p_{\theta}(t,x,\mu,.)\) w.r.t. some measure on \(\mathbb{R}^{p}\), and thus by minimizing over the parameters \(\theta\in\mathbb{R}^{D}\) the functional

\[\mathrm{J}(\theta)=\ \mathbb{E}_{\alpha\sim\pi_{\theta}}\Big{[}\int_{0}^{T}f (X_{t},\mathbb{P}_{X_{t}},\alpha_{t})\mathrm{d}t+g(X_{T},\mathbb{P}_{X_{T}}) \Big{]}. \tag{2.1}\]

Here \(\alpha\sim\pi_{\theta}\) means that at each time \(t\), the action \(\alpha_{t}\) is sampled (independently from \(W\)) from the probability distribution \(\pi_{\theta}(.|t,X_{t},\mathbb{P}_{X_{t}})\).

**Remark 2.1**.: _We may include a entropy (e.g. Shannon) regularizer term in the functional cost (2) as proposed in [19] for encouraging exploration of randomized policies. This can slightly help the convergence of the policy gradient algorithms by permitting the use of higher learning rates, but it turns out that it does not really improve the accuracy of the results. Here, we only consider exploration through the randomization of policies._

We have the gradient representation of \(\mathrm{J}\) as derived in [7]:

\[\mathrm{G}(\theta):= \ \nabla_{\theta}\mathrm{J}(\theta)\] \[= \ \mathbb{E}_{\alpha\sim\pi_{\theta}}\Big{[}\int_{0}^{T}\nabla_{ \theta}\log p_{\theta}(t,X_{t},\mathbb{P}_{X_{t}},\alpha_{t})\big{[}\mathrm{d }J_{\theta}(t,X_{t},\mathbb{P}_{X_{t}})+f(X_{t},\mathbb{P}_{X_{t}},\alpha_{t}) \mathrm{d}t\big{]}\] \[\qquad\qquad\qquad\qquad+\ \int_{0}^{T}\mathcal{H}_{\theta}[J_{ \theta}](t,X_{t},\mathbb{P}_{X_{t}})\mathrm{d}t\Big{]}, \tag{2.2}\]

where \(J_{\theta}:[0,T]\times\mathbb{R}^{d}\times\mathcal{P}_{2}(\mathbb{R}^{d})\to \mathbb{R}\) is the dynamic value function associated to (2), hence satisfying the property that

\[(\mathbf{MJ})\quad\{J_{\theta}(t,X_{t},\mathbb{P}_{X_{t}})+\int_{0}^{t}f(X_{ s},\mathbb{P}_{X_{s}},\alpha_{s})\mathrm{d}s,\ 0\leq t\leq T\}\ \ \text{is a martingale},\]

and \(\mathcal{H}_{\theta}\) is the operator specific to the mean-field framework, defined by

\[\mathcal{H}_{\theta}[\varphi](t,x,\mu)=\ \nabla_{\theta}\mathbb{E}_{ \xi\sim\mu}\Big{[}b_{\theta}(t,\xi,\mu)\cdot\partial_{\mu}\varphi(t,x,\mu)( \xi)+\frac{1}{2}\Sigma_{\theta}(t,\xi,\mu)\circ\partial_{\xi}\partial_{\mu} \varphi(t,x,\mu)(\xi)\Big{]}\ \in\mathbb{R}^{D}\]

with \(b_{\theta}(t,x,\mu)=\int_{A}b(t,x,\mu,a)\pi_{\theta}(\mathrm{d}a|t,x,\mu)\), \(\Sigma_{\theta}(t,x,\mu)=\int_{A}\sigma\sigma^{\tau}(t,x,\mu,a)\pi_{\theta}( \mathrm{d}a|t,x,\mu)\). Here \(\partial_{\mu}\varphi(t,x,\mu)(.)\) is the Lions-derivative with respect to \(\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})\), and it is a function from \(\mathbb{R}^{d}\) into \(\mathbb{R}^{d}\), and \(\mathbb{E}_{\xi\sim\mu}[.]\) means that the expectation is taken with respect to the random variable \(\xi\) distributed according to \(\mu\).

Notice that under the structure condition **(SC)**, we have

\[b_{\theta}(t,x,\mu)\ =\ \beta(t,x,\mu)+C_{\theta}(t,x,\mu),\quad\Sigma_{ \theta}(t,x,\mu)\ =\ \Sigma(t,x,\mu)+\vartheta_{\theta}(t,x,\mu)\]

where \(C_{\theta}(t,x,\mu):=\int_{A}C(t,a)\pi_{\theta}(\mathrm{d}a|t,x,\mu)\), \(\vartheta_{\theta}(t,x,\mu):=\int_{A}\vartheta(t,a)\pi_{\theta}(\mathrm{d}a|t,x,\mu)\) are known functions, and thus

\[\mathcal{H}_{\theta}[\varphi](t,x,\mu)=\ \nabla_{\theta}\mathbb{E}_{\xi\sim\mu} \Big{[}C_{\theta}(t,\xi,\mu)\cdot\partial_{\mu}\varphi(t,x,\mu)(\xi)+\frac{1} {2}\vartheta_{\theta}(t,\xi,\mu)\circ\partial_{\xi}\partial_{\mu}\varphi(t,x, \mu)(\xi)\Big{]}.\]

## 3 Parametrization of actor/critic functions with moment neural networks

A moment neural network function on \([0,T]\times\mathbb{R}^{d}\times\mathcal{P}_{2}(\mathbb{R}^{d})\) of order \(L\in\mathbb{N}^{*}\) is a parametric function in the form

\[\phi_{\eta}(t,x,\mu)=\ \Psi_{\eta}(t,x,\bar{\mathbf{\mu}}_{L}),\]
where \(\bar{\boldsymbol{\mu}}_{L}=(\bar{\mu}^{\ell})_{\ell\in\mathcal{L}}\), with \(\bar{\mu}^{\ell}=\mathbb{E}_{\xi\sim\mu}[\prod_{i=1}^{d}\xi_{i}^{\ell_{i}}]\) for \(\ell=(\ell_{i})_{i\in\mathbb{I}_{1},d\mathbb{I}}\in\mathcal{L}=\{\ell=(\ell_{1},\ldots,\ell_{d})\in\mathbb{N}^{d}:\sum_{i=1}^{d}\ell_{i}\leq L\}\) of cardinality \(L_{d}\), and \((t,x,y)\in[0,T]\times\mathbb{R}^{d}\times\mathbb{R}^{L_{d}}\mapsto\Psi_{\eta} (t,x,y)\) is a classical finite-dimensional feedforward neural network with parameters \(\eta\). Moment neural networks have been considered in [20] as a special case of cylindrical mean-field neural networks, and satisfy a universal approximation theorem for continuous functions on \([0,T]\times\mathbb{R}^{d}\times\mathcal{P}_{2}(\mathbb{R}^{d})\), see [15]. By abuse of notation and language, we identify \(\phi_{\eta}\) and \(\Psi_{\eta}\), and call them indifferently moment neural networks.

We shall parametrize the randomized policy (actor) by a Gaussian probability transition kernel in the form

\[\pi_{\theta}(.|t,x,\mu) = \mathcal{N}(m_{\theta}(t,x,\bar{\boldsymbol{\mu}}_{L}),\lambda \mathbb{I}_{p}),\]

where \(m_{\theta}\) is a moment neural network, hence with log density:

\[\log p_{\theta}(t,x,\mu,a) = -\frac{1}{2}\log(2\pi\lambda)-\frac{|a-m_{\theta}(t,x,\bar{ \boldsymbol{\mu}}_{L})|^{2}}{2\lambda},\]

and \(\lambda>0\) is a parameter for exploration. Notice that in this case, the known functions \(C_{\theta}\), \(\vartheta_{\theta}\) depend on \(\mu\) only though its \(L\) moments \(\bar{\boldsymbol{\mu}}_{L}\), and by misuse of notation we also write: \(C_{\theta}(t,x,\bar{\boldsymbol{\mu}}_{L})\), \(\vartheta_{\theta}(t,x,\bar{\boldsymbol{\mu}}_{L})\).

The value function (critic) is parametrized by a moment neural network \(\mathcal{J}_{\eta}(t,x,\mu)=\mathcal{J}_{\eta}(t,x,\bar{\boldsymbol{\mu}}_{L})\), and we notice that

\[\partial_{\mu}\mathcal{J}_{\eta}(t,x,\mu)(\xi)= D_{1}(\xi)\circ\nabla_{y}\mathcal{J}_{\eta}(t,x,\bar{\boldsymbol{\mu}}_{L})\] \[\partial_{\xi}\partial_{\mu}\mathcal{J}_{\eta}(t,x,\mu)(\xi)= D_{2}(\xi)\circ\nabla_{y}\mathcal{J}_{\eta}(t,x,\bar{\boldsymbol{\mu}}_{L}),\]

where \(D_{1}(\xi)\) is the matrix in \(\mathbb{R}^{L_{d}\times d}\), and \(D_{2}(\xi)\) is the tensor in \(\mathbb{R}^{L_{d}\times d\times d}\) with components

\[[D_{1}(\xi)]_{\ell i}= \ \ell_{i}\xi_{i}^{\ell_{i}-1}\prod_{k\neq i}\xi_{k}^{\ell_{k}}, \quad\text{ for }\xi=(\xi_{i})_{i\in[\mathbb{I},d]},\ \ell=(\ell_{i})_{i\in[\mathbb{I},d]},\] \[[D_{2}(\xi)]_{\ell ij}= \ \begin{cases}\ell_{i}(\ell_{i}-1)\xi_{i}^{\ell_{i}-2}\prod_{k\neq i }\xi_{k}^{\ell_{k}},\quad i=j\\ \ell_{i}\ell_{j}\xi_{i}^{\ell_{i}-1}\xi_{j}^{\ell_{j}-1}\prod_{k\neq i,j}\xi_ {k}^{\ell_{k}},\quad i\neq j.\end{cases}\]

The expression of the operator \(\mathcal{H}_{\theta}\) applied to the moment neural network critic function is then given by

\[\mathcal{H}_{\theta}[\mathcal{J}_{\eta}](t,x,\mu)=\ \nabla_{\theta}\Big{[} \mathbb{E}_{\xi\sim\mu}\big{[}D_{1}(\xi)C_{\theta}(t,\xi,\bar{\boldsymbol{ \mu}}_{L})+\frac{1}{2}D_{2}^{\intercal}(\xi)\circ\vartheta_{\theta}(t,\xi,\bar{ \boldsymbol{\mu}}_{L})\big{]}\cdot\nabla_{y}\mathcal{J}_{\eta}(t,x,\bar{ \boldsymbol{\mu}}_{L})\Big{]}. \tag{3.1}\]

Here \(D_{2}^{\intercal}(\xi)\) is the tensor in \(\mathbb{R}^{d\times d\times L_{d}}\) with components \([D_{2}^{\intercal}(\xi)]_{ij\ell}=[D_{2}(\xi)]_{\ell ij}\).

In the algorithm, we shall use the expectation of \(\mathcal{H}_{\theta}\), which is given from (3) by

\[\overline{\mathcal{H}}_{\theta}[\mathcal{J}_{\eta}](t,\mu):= \ \mathbb{E}_{\xi\sim\mu}\big{[}\mathcal{H}_{\theta}[\mathcal{J}_{ \eta}](t,\xi,\mu)\big{]}\] \[= \ \nabla_{\theta}\Big{[}\mathbb{E}_{\xi\sim\mu}\big{[}D_{1}(\xi)C_ {\theta}(t,\xi,\bar{\boldsymbol{\mu}}_{L})+\frac{1}{2}D_{2}^{\intercal}(\xi) \circ\vartheta_{\theta}(t,\xi,\bar{\boldsymbol{\mu}}_{L})\big{]}\cdot\nabla_{y} \mathbb{E}_{\xi\sim\mu}\big{[}\mathcal{J}_{\eta}(t,\xi,\bar{\boldsymbol{\mu}}_{L })\big{]}\Big{]}.\]

**Remark 3.1**.: **1.** _For complexity argument, it is crucial to rely on the above expression of the operator \(\mathcal{H}_{\theta}\) where the differentiation is taken on the expectation \(\mathbb{E}_{\xi\sim\mu}[.]\), and not the reversal: expectation of the differentiation. Indeed, in the latter case, after empirical approximation of the expectation with \(M\) samples \(\xi^{j}\sim\mu\), \(j=1,\ldots,M\), one should compute by automatic differentiation_

\[\nabla_{\theta}\big{[}D_{1}(\xi^{j})C_{\theta}(t,\xi^{j},\bar{\boldsymbol{\mu}} _{L})+\frac{1}{2}D_{2}^{\intercal}(\xi^{j})\circ\vartheta_{\theta}(t,\xi^{j}, \bar{\boldsymbol{\mu}}_{L})\big{]},\quad\nabla_{y}\mathcal{J}_{\eta}(t,\xi^{j}, \bar{\boldsymbol{\mu}}_{L}),\quad j=1,\ldots,M,\]

_which is very costly as \(M\) is of order \(10^{4}\). In the former case, \(\overline{\mathcal{H}}_{\theta}[\mathcal{J}_{\eta}](t,\mu)\) is approximated by automatic differentiation via_

\[\widehat{\mathcal{H}}_{\theta}^{M}[\mathcal{J}_{\eta}](t,\mu) := \ \nabla_{\theta}\Big{[}\tfrac{1}{M}\sum_{j=1}^{M}D_{1}(\xi^{j})C_{ \theta}(t,\xi^{j},\bar{\boldsymbol{\mu}}_{L})+\tfrac{1}{2}D_{2}^{\intercal}(\xi^ {j})\circ\vartheta_{\theta}(t,\xi^{j},\bar{\boldsymbol{\mu}}_{L})\cdot\nabla_{y} \tfrac{1}{M}\sum_{j=1}^{M}\mathcal{J}_{\eta}(t,\xi^{j},\bar{\boldsymbol{\mu}}_{L })\Big{]}.\]
hence saving an order \(M\) for the complexity cost. In theory, it is also possible to choose other networks for taking into account the dependency on the distribution \(\mu\): the cylindrical network proposed in [15] could be used but some automatic differentiation are then requested to calculate \(\partial_{\mu}\mathcal{J}_{\eta}(t,x,\mu)(\xi)\) and \(\partial_{\xi}\partial_{\mu}\mathcal{J}_{\eta}(t,x,\mu)(\xi)\) for each sample of \(\xi\) leading to an explosion in the computation time._

**2.** _In order to calculate the term \(\nabla_{y}\mathcal{J}_{\eta}\), it is necessary to explore different initial distributions, otherwise \(\mathcal{J}_{\eta}\) only depends on \(t\) and \(x\) at convergence and the gradient is impossible to estimate._

## 4 Algorithm

The actor-critic method consists in two optimization stages that are performed alternately:

1. _Policy evaluation_: given an actor policy \(\pi_{\theta}\), evaluate its cost functional with the critic function \(\mathcal{J}_{\eta}\) that minimizes the loss function arising from the martingale property **(MJ)** after time discretization of the interval \([0,T]\) with the time grid \(\{t_{k}=k\Delta t,k=0,\ldots,n\}\): \[L^{PE}(\eta)=\ \mathbb{E}_{\alpha\sim\pi_{\theta}}\Big{[}\sum_{k=0}^{n-1} \Big{|}g_{t_{n}}+\sum_{l=k}^{n-1}f_{t_{l}}\Delta t-\mathcal{J}_{\eta}(t_{k},X _{t_{k}},\mu_{t_{k}})\Big{|}^{2}\Delta t\Big{]},\] where we set \(\mu_{t_{l}}=\mathbb{P}_{X_{t_{l}}}\) for the law of \(X_{t_{l}}\), and \(f_{t_{l}}=f(X_{t_{l}},\mu_{t_{l}},\alpha_{t_{l}})\) as the output cost at time \(t_{l}\) for input state \(X_{t_{l}}\), law \(\mu_{t_{l}}\), action \(\alpha_{t_{l}}\sim\pi_{\theta}(.|t_{l},X_{t_{l}},\mu_{t_{l}})\), and \(g_{t_{n}}=g(X_{t_{n}},\mu_{t_{n}})\) the terminal output cost for input \(X_{t_{n}}\), \(\mu_{t_{n}}\).
2. _Policy gradient_: given a critic cost function \(\mathcal{J}_{\eta}\), update the parameter \(\theta\) of the actor by stochastic gradient descent by using the gradient, which is given from (2) and after time discretization by \[G(\theta)=\ \mathbb{E}_{\alpha\sim\pi_{\theta}}\Big{[}\sum_{k=0}^{n- 1}\nabla_{\theta}\log p_{\theta}(t_{k},X_{t_{k}},\mu_{t_{k}},\alpha_{t_{k}}) \big{[}\mathcal{J}_{\eta}(t_{k+1},X_{t_{k+1}},\mu_{t_{k+1}})-\mathcal{J}_{\eta }(t_{k},X_{t_{k}},\mu_{t_{k}})\] \[\qquad\qquad\qquad+\ f_{t_{k}}\Delta t\big{]}+\mathcal{H}_{\theta}[ \mathcal{J}_{\eta}](t_{k},X_{t_{k}},\mu_{t_{k}})\Big{]}.\]

In the practical implementation, we proceed as follows for each epoch e (gradient iteration descent) with a given exploration parameter \(\lambda(e)\) decreasing to \(0\):

* We start with a batch \(N\) (of order \(10\)) of initial distributions \(\mu_{0}^{i}\), \(i=1,\ldots,N\), e.g. Gaussian distributions by varying the mean and std deviations parameters, and sample \(X_{0}^{i,j}\sim\mu_{0}^{i}\), \(j\)\(=1,\ldots,M\) with \(M\) of order \(10^{4}\). If our ultimate goal is to learn the optimal control and function value for other families of initial distributions, the initial distributions should be sampled accordingly.
* We then run by forward induction in time: for \(k=0,\ldots,n-1\):
* Empirical estimate of \(\mu_{t_{k}}^{i}\) from \((X_{t_{k}}^{i,j})_{j\in[\![1,M]\!]}\), for \(i=1,\ldots,N\).
* Sample \(\alpha_{t_{k}}^{i,j}\sim\pi_{\theta}(.|t_{k},X_{t_{k}}^{i,j},\mu_{t_{k}}^{i})\), \(i\in[\![1,N]\!]\), \(j\in[\![1,M]\!]\) using the exploration parameter \(\lambda(e)\)
* Observe running cost \(f_{t_{k}}^{i,j}=f(X_{t_{k}}^{i,j},\mu_{t_{k}}^{i},\alpha_{t_{k}}^{i,j})\), and next state \(X_{t_{k+1}}^{i,j}\), \(i\in[\![1,N]\!]\), \(j\in[\![1,M]\!]\)
* Observe final cost \(g_{t_{n}}^{i,j}=g(X_{t_{n}}^{i,j},\mu_{t_{n}}^{i})\), \(i\in[\![1,N]\!]\), \(j\in[\![1,M]\!]\)
* Compute the empirical mean approximation of \(L^{PE}(\eta)\) on all initial distributions \(\mu_{0}^{i}\), \(i\in[\![1,N]\!]\): \[\widetilde{L}_{M}^{PE}(\eta)=\ \frac{1}{MN}\sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{k=0}^{n- 1}\big{|}g_{t_{n}}^{i,j}+\sum_{l=k}^{n-1}f_{t_{l}}^{i,j}\Delta t-\mathcal{J}_ {\eta}(t_{k},X_{t_{k}}^{i,j},\mu_{t_{k}}^{i})\big{|}^{2}\Delta t,\] and update the critic parameter by \[\eta\longleftarrow\ \eta-\rho^{C}\nabla_{\eta}\tilde{L}_{M}^{PE}(\eta),\] where \(\rho^{C}\) is a learning rate. Notice that the gradient is calculated by automatic differentiation.
* Compute the empirical mean approximation of \(G(\theta)\) on all initial distributions \(\mu_{0}^{i}\), \(i\in[\![1,N]\!]\): \[\widetilde{G}_{M}(\theta)=\ \nabla_{\theta}\frac{1}{MN}\sum_{i=1}^{N}\sum_{j=1}^{M} \Big{\{}\sum_{k=0}^{n-1}\log p_{\theta}(t_{k},X_{t_{k}}^{i,j},\mu_{t_{k}}^{i}, \alpha_{t_{k}}^{i,j})\big{[}\mathcal{J}_{\eta}(t_{k+1},X_{t_{k+1}}^{i,j},\mu_{ t_{k+1}}^{i})-\mathcal{J}_{\eta}(t_{k},X_{t_{k}}^{i,j},\mu_{t_{k}}^{i})\] \[+\ f_{t_{k}}^{i,j}\Delta t\big{]}\Big{\}}+\widetilde{\mathcal{H}}_{ \theta}^{M}[\mathcal{J}_{\eta}]\] where \[\widetilde{\mathcal{H}}_{\theta}^{M}[\mathcal{J}_{\eta}]= \nabla_{\theta}\Big{(}\frac{1}{N}\sum_{i=1}^{N}\Big{[}\sum_{k=0} ^{n-1}\Big{(}\frac{1}{M}\sum_{j=1}^{M}D_{1}(X_{t_{k}}^{i,j})C_{\theta}(t_{k},X _{t_{k}}^{i,j},\mu_{t_{k}}^{i})+\frac{1}{2}D_{2}^{*}(X_{t_{k}}^{i,j})\circ \vartheta_{\theta}(t_{k},X_{t_{k}}^{i,j},\mu_{t_{k}}^{i})\Big{)}\cdot\] \[\nabla_{y}\Big{(}\frac{1}{M}\sum_{j=1}^{M}\mathcal{J}_{\eta}(t_{k},X _{t_{k}}^{i,j},\mu_{t_{k}}^{i})\Big{)}\Big{]}\Big{)}\] and update the actor parameter by \[\theta\longleftarrow\ \theta-\rho^{A}\widetilde{G}_{M}(\theta),\] where \(\rho^{A}\) is a learning rate. Again for efficiency, it is crucial to compute by automatic differentiation the gradient after computing all the different expectations as in Remark 3.1.

The output \((\theta^{*},\eta^{*})\) are the optimal parameters obtained at convergence of the algorithm.

**Remark 4.1**.: _Compared to classical actor critic algorithm where one samples a trajectory for a given distribution, here the batch version of the algorithm consists in sampling and following \(N\) distributions together to estimate the gradients._

**Remark 4.2**.: _In order to check that the algorithm has effectively converged to the solution, we can use the calculated control \(m_{\theta^{*}}(t,x,\mu)\) and apply it from different initial distributions \(\mu_{0}\) sampled as \((X_{0}^{j})_{j\in[\![1,M]\!]}\) in a time discretized version of \((\ref{eq:1})\). Taking discrete expectation, we can compare the result obtained to \(\frac{1}{M}\sum_{j=1}^{M}\mathcal{J}_{\eta^{*}}(0,X_{0}^{j},\mu_{0})\). When results are very close, we can suppose that the algorithm has effectively converged to the right solution._

**Remark 4.3**.: _In the case where we know a priori that the running cost and terminal cost functions depend on the probability distribution \(\mu\) only via its moments \(\bar{\mu}_{L}\), then we only need to estimate the moments of \(\mu_{t_{k}}^{i}\) from \((X_{t_{k}}^{i,j})_{j\in[\![1,M]\!]}\), since all the other coefficients in the algorithm depend upon the measure via its moments._

**Remark 4.4**.: _When \(C_{\theta}(t,x,\mu):=\int_{A}C(t,a)\pi_{\theta}(\mathrm{d}a[t,x,\mu)\), \(\vartheta_{\theta}(t,x,\mu):=\int_{A}\vartheta(t,a)\pi_{\theta}(\mathrm{d}a[t, x,\mu)\) are not analytically explicit, it is always possible to estimate them numerically for example using a numerical quadrature or a quasi Monte carlo/Monte-Carlo method but with some non negligible extra costs._

## 5 Numerical results

Throughout this section, we use moment neural networks with 3 hidden layers and 20 neurons on each layer, and choose the activation function \(\tanh\). The exploration parameter \(\lambda\) is a function of the number of gradient descent iterations (epoch number \(e\leq\hat{N}\)):

\[\lambda(e)=(\bar{\lambda}-\underline{\lambda})\Big{(}1-S\big{(}\frac{20e-10 \hat{N}}{\hat{N}}\big{)}\Big{)}+\underline{\lambda},\]

where \(\underline{\lambda}=0.0001\) and \(\bar{\lambda}=0.1\) and \(S\) is the sigmoid function: \(S(x)=\frac{1}{1+\exp(-x)}\). In other words, it is chosen so that the exploration period with \(\lambda\) close to \(0.1\) is long enough, then \(\lambda\) slowly decreases to \(0.0001\) and stays close to that value long enough. This fonction is plotted on Figure 1.
At a gradient descent iteration \(e\), we sample the control as:

\[\pi_{\theta}(.|t,x,\mu)\sim\mathcal{N}(m_{\theta}(t,x,\mu),\lambda(e)\mathbb{I}_{ p}).\]

During the gradient descent algorithm, we use the ADAM optimizer [13]. We point out that it is crucial to use two timescales approach (see [1]), for the learning rates \(\rho^{C}\), \(\rho^{A}\) of the critic and actor updates: \(\rho^{C}\) should be at least one order of magnitude higher than \(\rho^{A}\) to get good convergence, hence the approximate critic function should evolve faster. We take a batch size \(N=10\) while the number of samples to estimate distributions is taken equal to \(M=10000\) or \(M=20000\) depending on the examples.

In the tables and figures below, we give the average analytic solution "Anal" at \(t=0\), i.e. \(\mathbb{E}_{X_{0}\sim\mu_{0}}[V(0,X_{0},\mu_{0})]\), and the average calculated value function "Calc": \(\mathbb{E}_{X_{0}\sim\mu_{0}}[\mathcal{J}_{\eta^{*}}(0,X_{0},\mu_{0})]\) obtained by the algorithm at \(t=0\), by varying the initial distributions \(\mu_{0}\). The MSE is the mean square error between the analytic and the critic value computed at \(t=0\), i.e. \(\mathbb{E}_{X_{0}\sim\mu_{0}}\big{|}\mathcal{J}_{\eta^{*}}(0,X_{0},\mu_{0})- V(0,X_{0},\mu_{0})\big{|}^{2}\), and the relative error is

\[\text{RelError}=\ \frac{\mathbb{E}_{X_{0}\sim\mu_{0}}\big{[}\mathcal{J}_{ \eta^{*}}(0,X_{0},\mu_{0})-V(0,X_{0},\mu_{0})\big{]}}{\mathbb{E}_{X_{0}\sim\mu _{0}}[V(0,X_{0},\mu_{0})]}\]

We shall also plot in the one-dimensional case \(d=1\), and \(A=\mathbb{R}\), the trajectories of the optimal control \(t\mapsto\alpha_{t}^{*}\) vs the ones obtained from moment neural networks, i.e., \(t\mapsto m_{\theta^{*}}(t,X_{t},(\mathbb{E}[X_{t}^{\ell}])_{\ell\in[1,L]})\).

All training times are calculated on a on GPU NVidia V100 32Go graphic card.

We consider four examples with control on the drift, including multidimensional setting and nonlinear quadratic mean-field control, and one example with controlled volatility, for which we have analytic solutions to be compared with the approximations calculated from our actor critic algorithm.

### Examples with controlled drift

In the four examples of this paragraph, we take \(\vartheta(t,a)\equiv 0\), \(C(t,a)=a\) and so \(\vartheta_{\theta}\equiv 0\), \(C_{\theta}=m_{\theta}\).

#### 5.1.1 Systemic risk model in one dimension

We consider the model in [4]:

\[\begin{cases}b(x,\mu,a)&=\ \kappa(\bar{\mu}-x)+a,\quad\sigma\ \text{positive constant}\\ f(x,\mu,a)&=\ \frac{1}{2}a^{2}-qa(\bar{\mu}-x)+\frac{p}{2}(\bar{\mu}-x)^{2}, \qquad g(x,\mu)\ =\ \frac{c}{2}(x-\bar{\mu})^{2},\end{cases}\]

for \((x,\mu,a)\in\mathbb{R}\times\mathcal{P}_{2}(\mathbb{R})\times\mathbb{R}\), with some positive constants \(\kappa\), \(q\), \(p\), \(c>0\), \(q^{2}\leq p\). Here we denote by \(\bar{\mu}:=\mathbb{E}_{\xi\sim\mu}[\xi]\).

In this linear quadratic (LQ) model, the value function is explicitly given by

\[V(t,x,\mu)=\ K(t)(x-\bar{\mu})^{2}+\sigma^{2}R(t), \tag{5.1}\]

Figure 1: Exploration function \(\lambda\) for a number of epochs \(\hat{N}=6000\).
where

\[K(t)=-\frac{1}{2}\Big{[}\kappa+q-\sqrt{\Delta}\frac{\sqrt{\Delta} \sinh(\sqrt{\Delta}(T-t))+(\kappa+q+c)\cosh(\sqrt{\Delta}(T-t))}{\sqrt{\Delta} \cosh(\sqrt{\Delta}(T-t))+(\kappa+q+c)\sinh(\sqrt{\Delta}(T-t))}\Big{]},\]

with \(\sqrt{\Delta}=\sqrt{(\kappa+q)^{2}+p-q^{2}}\), and

\[R(t)=\ \frac{\sigma^{2}}{2}\ln\Big{[}\cosh(\sqrt{\Delta}(T-t))+ \frac{\kappa+q+c}{\sqrt{\Delta}}\sinh(\sqrt{\Delta}(T-t))\Big{]}-\frac{\sigma ^{2}}{2}(\kappa+q)(T-t),\]

while the optimal control is given by

\[\alpha_{t}^{*}=\ (2K(t)+q)(\mathbb{E}[X_{t}]-X_{t}),\quad 0\leq t\leq T.\]

The parameters of the model are fixed to the following values: \(\kappa=0.6\), \(\sigma=1\), \(p=c=2\), \(q=0.8\), \(T=1\). We take a number of time steps \(n=100\), \(M=10000\). At each gradient iteration, the initial distribution is sampled with

\[X_{0}\sim\ \mu_{0}\ =\ \upsilon_{0}\ \mathcal{N}(0,1),\]

where \(\upsilon_{0}^{2}\) is sampled at each iteration according to the uniform distribution on \([0,1]\) for each element of the batch. In the table, we give the results obtained in simulation with \(\upsilon_{0}^{2}\in\{0,\frac{1}{10},\ldots,1\}\) after training with \(L=2\) moments, using \(\hat{N}=6000\) gradient iterations.

On Figure 2, we plot 3 trajectories of the optimal control and the ones calculated with moment neural networks, and we observe that the control is very well estimated.

In this LQ example, we know that the suitable number of moments to take is \(L=2\). In the real test case, we do not know which \(L\) to take. In Table 2, we give the results with \(L=4\), which shows that the results are also very accurate and do not depend on \(L\) being small.

\begin{table}
\begin{tabular}{|c|c|c|c|c|c|} \hline \(\upsilon_{0}^{2}\) & 0. & 0.1 & 0.2 & 0.3 & 0.4 \\ \hline Anal & 0.3870 & 0.4095 & 0.4321 & 0.4546 & 0.4772 \\ \hline Calc & 0.3958 & 0.4198 & 0.4421 & 0.4642 & 0.4875 \\ \hline MSE & 0.0001 & 0.0001 & 0.0002 & 0.0002 & 0.0004 \\ \hline \hline \(\upsilon_{0}^{2}\) & 0.5 & 0.6 & 0.7 & 0.8 & 0.9 \\ \hline Anal & 0.4997 & 0.5223 & 0.5448 & 0.5674 & 0.5900 \\ \hline Calc & 0.5112 & 0.5341 & 0.5593 & 0.5858 & 0.6082 \\ \hline MSE & 0.0005 & 0.0005 & 0.0006 & 0.0005 & 0.0007 \\ \hline \end{tabular}
\end{table}
Table 1: Results for the systemic model using \(L=2\), \(\rho^{A}=0.0005\), \(\rho^{C}=0.01\). Training time is 106863s.

Figure 2: Trajectories of control with \(\upsilon_{0}^{2}=0.9\)
#### 5.1.2 An optimal trading example

We consider an optimal trading model taken from [7]:

\[\begin{cases}b(x,\mu,a)&=\ a,\quad\sigma\ \text{positive constant}\\ f(x,\mu,a)&=\ a^{2}+2Pa,\qquad g(x,\mu)\ =\ \gamma(x-\bar{\mu})^{2},\end{cases}\]

for \((x,\mu,a)\in\mathbb{R}\times\mathcal{P}_{2}(\mathbb{R})\times\mathbb{R}\), with \(P>0\) the constant transaction price per trading, and \(\gamma>0\) the risk aversion parameter.

In this LQ framework, the value function has the form as in (5.1.1) with

\[K(t)\ =\ \frac{\gamma}{1+\gamma(T-t)},\qquad R(t)\ =\ \sigma^{2}\log(1+\gamma(T-t))-P^{2}(T-t),\]

while the optimal control is given by

\[\alpha_{t}^{\star}=\ -K(t)(X_{t}-\mathbb{E}[X_{t}])-P,\quad 0\leq t\leq T.\]

We take the following parameters : \(P=3,\gamma=3,\sigma=1,T=0.5\), and \(n=100\), \(M=10000\). At each gradient iteration, the initial distribution is sampled with

\[X_{0}\sim\ \mu_{0}=\bar{\mu}_{0}+\upsilon_{0}\mathcal{N}(0,1),\]

where \((\bar{\mu}_{0},\upsilon_{0}^{2})\) are sampled from \((0.4\mathcal{U}([0,1]),0.5\mathcal{U}([0,1]))\) for each element of the batch.

The relative error is plotted on Figure 3 by varying \((\bar{\mu}_{0},\upsilon_{0})\), while the trajectories of the control (optimal vs moment neural netwok) are plotted on Figure 4. Again, in this LQ example, the suitable number of moments to take is \(L=2\), and when we increase \(L\), convergence is more difficult to achieve and results become more instable and may depend on the run with the same hyper-parameters. Nevertheless we manage to obtain very good results with \(L=4\) as shown on Figure 3.

\begin{table}
\begin{tabular}{|c|c|c|c|c|c|} \hline \(\upsilon_{0}^{2}\) & 0. & 0.1 & 0.2 & 0.3 & 0.4 \\ \hline Anal & 0.3869 & 0.4095 & 0.4320 & 0.4546 & 0.4771 \\ \hline Calc & 0.3917 & 0.4159 & 0.4376 & 0.4588 & 0.4801 \\ \hline MSE & 0.0000 & 0.0000 & 0.0001 & 0.0001 & 0.0001 \\ \hline \hline \(\upsilon_{0}^{2}\) & 0.5 & 0.6 & 0.7 & 0.8 & 0.9 \\ \hline Anal & 0.4997 & 0.5222 & 0.5448 & 0.5674 & 0.5899 \\ \hline Calc & 0.5023 & 0.5249 & 0.5471 & 0.5688 & 0.5891 \\ \hline MSE & 0.0000 & 0.0001 & 0.0001 & 0.0002 & 0.0003 \\ \hline \end{tabular}
\end{table}
Table 2: Results for the systemic case using \(L=4\), \(\rho^{A}=0.0005\), \(\rho^{C}=0.01\). Training time is 115183s.
#### 5.1.3 A non linear quadratic mean-field control

We construct an ad-hoc mean-field control model with

\[\begin{cases}b(t,x,\mu,a)&=\ \beta(t,x,\mu)+a,\quad\sigma\ \mbox{positive constant}\\ f(x,\mu,a)&=F(t,x,\mu)+\frac{1}{2}|a|^{2},\quad g(x,\mu)\ =\mathbb{E}_{\xi\sim\mu}[w(x-\xi)] \end{cases}\]

for some smooth \(C^{2}\) even function \(w\) on \(\mathbb{R}\), e.g. \(w(x)=\cos(x)\), and \(F\) is a function to be chosen later.

In this case, the optimal feedback control valued in \(A=\mathbb{R}\) is given by

\[\mathfrak{a}^{\star}(t,x,\mu) = \hat{\mathrm{a}}(t,x,\mathcal{U}(t,x,\mu))\ =\ -\mathcal{U}(t,x,\mu)\ =\ -\partial_{\mu}v(t,\mu)(x)\] \[\quad\mbox{with }v(t,\mu)\ =\ \mathbb{E}_{\xi\sim\mu}[V(t,\xi,\mu)],\]

and \(V\) is solution to the Master Bellman equation (see section 6.5.2 in [2]):

\[\partial_{t}V(t,x,\mu)+\big{(}\beta(t,x,\mu)-\mathcal{U}(t,x,\mu )\big{)}\partial_{x}V(t,x,\mu)+\frac{\sigma^{2}}{2}\partial_{xx}^{2}V(t,x,\mu)\] \[+\ \mathbb{E}_{\xi\sim\mu}\Big{[}\big{(}\beta(t,\xi,\mu)- \mathcal{U}(t,\xi,\mu)\big{)}\partial_{\mu}V(t,x,\mu)(\xi)+\frac{\sigma^{2}}{2 }\partial_{x^{\prime}}\partial_{\mu}V(t,x,\mu)(\xi)\Big{]}\] \[\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+F(t,x,\mu)+\frac{1}{2}| \mathcal{U}(t,x,\mu)|^{2}=\ 0, \tag{5.2}\]

with the terminal condition \(V(T,x,\mu)=g(x,\mu)\).

We look for a solution to the Master equation in the form: \(V(t,x,\mu)=e^{T-t}\mathbb{E}_{\xi\sim\mu}[w(x-\xi)]\). For such function \(V\), we have \(\partial_{t}V(t,x,\mu)=-V\),

\[\partial_{x}V(t,x,\mu) = e^{T-t}\mathbb{E}_{\xi\sim\mu}[w^{\prime}(x-\xi)],\quad\partial_ {xx}^{2}V(t,x,\mu)\ =\ e^{T-t}\mathbb{E}_{\xi\sim\mu}[w^{\prime\prime}(x-\xi)]\] \[\partial_{\mu}V(t,x,\mu)(\xi) = -e^{T-t}w^{\prime}(x-\xi),\quad\partial_{x^{\prime}}\partial_{ \mu}V(t,x,\mu)(\xi)\ =\ e^{T-t}w^{\prime\prime}(x-\xi),\]

and

\[\mathcal{U}(t,x,\mu) = e^{T-t}\mathbb{E}_{\xi\sim\mu}[w^{\prime}(x-\xi)-w^{\prime}(\xi- x)]\ =\ 2e^{T-t}\mathbb{E}_{\xi\sim\mu}[w^{\prime}(x-\xi)]\ =\ 2 \partial_{x}V(t,x,\mu).\]

since \(w\) is even. By plugging these derivatives expressions of \(V\) into the l.h.s. of (5.2), we then see that by choosing \(F\) equal to

\[F(t,x,\mu) = e^{T-t}\mathbb{E}_{\xi\sim\mu}\Big{[}(w-\sigma^{2}w^{\prime\prime })(x-\xi)+(\beta(t,\xi,\mu)-\beta(t,x,\mu))w^{\prime}(x-\xi)\Big{]}\] \[-\ 2e^{2(T-t)}\mathbb{E}_{(\xi,\xi^{\prime})\sim\mu\otimes\mu} \big{[}w^{\prime}(x-\xi)w^{\prime}(\xi-\xi^{\prime})\big{]},\]

the function \(V\) satisfies the Master Bellman equation, hence is the value function to the mean-field control problem.
Actually, with the choice of \(w(x)=\cos(x)\), and using trigonometric relations, we have

\[F(t,x,\mu) = \cos(x)\Big{[}e^{T-t}\left((1+\sigma^{2})\mathbb{E}_{\xi\sim\mu}[ \cos(\xi)]+\mathbb{E}_{\xi\sim\mu}[\sin(\xi)\beta(t,\xi,\mu)]-\beta(t,x,\mu) \mathbb{E}_{\xi\sim\mu}[\sin(\xi)]\right)\] \[-\ 2e^{2(T-t)}\left(\mathbb{E}_{\xi\sim\mu}[\sin(\xi)\xi)[ \mathbb{E}_{\xi\sim\mu}[\sin(\xi)]-\mathbb{E}_{\xi\sim\mu}[\sin^{2}(\xi)] \mathbb{E}_{\xi\sim\mu}[\cos(\xi)]\right)\Big{]}\] \[+\ \sin(x)\Big{[}e^{T-t}\left((1+\sigma^{2})\mathbb{E}_{\xi\sim\mu}[ \sin(\xi)]-\mathbb{E}_{\xi\sim\mu}[\beta(t,\xi,\mu)\cos(\xi)]+\beta(t,x,\mu) \mathbb{E}_{\xi\sim\mu}[\cos(\xi)]\right)\] \[-\ 2e^{2(T-t)}\left(\mathbb{E}_{\xi\sim\mu}[\sin(\xi)\cos(\xi)] \mathbb{E}_{\xi\sim\mu}[\cos(\xi)]-\mathbb{E}_{\xi\sim\mu}[\cos^{2}(\xi)] \mathbb{E}_{\xi\sim\mu}[\sin(\xi)]\right)\Big{]}.\]

For the test case, we take

\[\beta(t,x,\mu)=\ \kappa(\bar{\mu}-x), \tag{5.3}\]

with the parameters \(\kappa=\sigma=1\), \(T=0.4\), \(n=40\), \(M=20000\). At each gradient iteration, the initial distribution is sampled with

\[X_{0}\sim\ \mu_{0}=\bar{\mu}_{0}+\upsilon_{0}\mathcal{N}(0,1)\]

where \((\bar{\mu}_{0},\upsilon_{0}^{2})\) is sampled from \((0.2\mathcal{U}([0,1]),0.5\mathcal{U}([0,1]))\) for each element of the batch. On Figure 5, we give the analytic solution depending on \((\bar{\mu}_{0},\upsilon_{0}^{2})\) and the relative error obtained by the algorithm. We observe that the results are quite accurate.

We plot in Figure 6 the trajectories of the optimal control vs the moment neural network, and observe that they are very close.
Figure 7 shows that using \(L=2\) or \(3\) is optimal in terms of relative error, while the convergence is more difficult to achieve for high values of \(L\).

#### 5.1.4 A multi-dimensional LQ example

We consider a multi-dimensional extension of the LQ systemic risk model in Section 5.1.1, by supposing that on each dimension, the dynamics satisfies the same equation with independent Brownian motions, and that the cost functions are the sum over each component of the cost function in the univariate case. In this case, the value function is given by \(V(t,x,\mu)=\sum_{i=1}^{d}V_{1}(t,x_{i},\mu_{i})\), for \(t\in[0,T]\), \(x=(x_{i})_{i\in[1,d]}\)\(\in\mathbb{R}^{d}\), \(\mu_{i}\) is the \(i\)-th marginal law of \(\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})\), and \(V_{1}\) is the value function in the univariate model given by (5.1.1).

We keep the same parameters as in Section 5.1.1, with a number of time steps \(n=50\), \(M=10000\), \(L=2\), and test in dimension \(d=2\) and \(3\). At each gradient iteration, the initial distribution is sampled from

\[X_{0}\sim\ \mu_{0}=\mathcal{N}(0,\upsilon_{0}),\]

where \(\upsilon_{0}\) is the diagonal \(d\times d\)-matrix with diagonal elements \(\upsilon_{0,i}\) sampled from \(u_{i}\mathcal{U}([0,1])\), with constants \(u_{i}\in[0,1]\), \(i=1,\ldots,d\), for each element of the batch.

We plot in Figure 8 the relative error in dimension \(d=2\) by varying \((\upsilon_{0,1},\upsilon_{0,2})\), and for \(L=2\) and \(L=4\). In Figure 9, we plot the relative error in dimension \(d=3\) by varying \((\upsilon_{0,1},\upsilon_{0,2})\), and for \(\upsilon_{0,3}\)\(=0\).

Figure 7: Relative error in a non LQ model for different values of \(L\). \(\hat{N}=9000\) gradient iterations, \(\rho^{A}=0.0005\), \(\rho^{C}=0.02\).
### A non LQ example with controlled volatility

We consider a one-dimensional model with

\[\begin{cases}b(t,x,\mu,a)&=\ \beta(t,x,\mu),\quad\sigma(t,x,\mu,a)\ =\ a,\\ f(x,\mu,a)&=\ F(t,x,\mu)+\frac{1}{2}P|a|^{2}-a,\quad g(x,\mu)\ =\ \mathbb{E}_{\xi \sim\mu}[w(x-\xi)],\end{cases}\]

where \(P\) is a positive constant, \(w\) is a smooth \(C^{2}\) even function on \(\mathbb{R}\), e.g. \(w(x)=\cos(x)\), and \(F\) is a function to be chosen later. Notice that \(C_{\theta}\equiv 0\), and

\[\vartheta_{\theta}(t,x,\mu) = m_{\theta}(t,x,\mu)^{2}+\lambda,\]

where \(\lambda=\lambda(e)\) (depending on the epoch \(e\)) is the exploration parameter of \(\pi_{\theta}(.|t,x,\mu)=\mathcal{N}(m_{\theta}(t,x,\mu),\lambda)\).

In this model, the optimal feedback control valued in \(A=\mathbb{R}\) is given by

\[\mathfrak{a}^{*}(t,x,\mu)=\ \hat{\mathrm{a}}(t,x,\partial_{x} \mathcal{U}(t,x,\mu))\ =\ \frac{1}{P+\partial_{x}\mathcal{U}(t,x,\mu)},\] \[\text{with}\ \ \mathcal{U}(t,x,\mu)\ =\ \partial_{\mu}v(t,\mu)(x),\quad v(t,\mu)\ =\ \mathbb{E}_{\xi\sim\mu}[V(t,\xi,\mu)],\]

Figure 8: Relative error in dimension \(d\) = 2, with \(\rho^{C}=0.01\).

Figure 9: Relative error for \(L=2\) in dimension \(d=3\), with \(\hat{N}=9000\) gradient iterations, \(\rho^{A}=0.0005\), \(\rho^{C}=0.01\). Training time is 144000s.
and \(V\) is solution to the Master Bellman equation:

\[\partial_{t}V(t,x,\mu)+\beta(t,x,\mu)\partial_{x}V(t,x,\mu)+\frac{1} {2}\frac{1}{(P+\partial_{x}\mathcal{U}(t,x,\mu))^{2}}\partial_{xx}^{2}V(t,x,\mu)\] \[+\,\mathbb{E}_{\xi\sim\mu}\Big{[}\beta(t,\xi,\mu)\partial_{\mu}V (t,x,\mu)(\xi)+\frac{1}{2}\frac{1}{(P+\partial_{x}\mathcal{U}(t,\xi,\mu))^{2}} \partial_{x^{\prime}}\partial_{\mu}V(t,x,\mu)(\xi)\Big{]}\] \[+F(t,x,\mu)+\frac{1}{2}\frac{P}{(P+\partial_{x}\mathcal{U}(t,x, \mu))^{2}}-\frac{1}{P+\partial_{x}\mathcal{U}(t,x,\mu)}=\;0, \tag{5.4}\]

with the terminal condition \(V(T,x,\mu)=g(x,\mu)\).

We look for a solution to the Master equation in the form: \(V(t,x,\mu)=e^{T-t}\mathbb{E}_{\xi\sim\mu}[w(x-\xi)]\), and by similar calculations as in Section 5.1.3, we would have \(\mathcal{U}(t,x,\mu)=2\partial_{x}V(t,x,\mu)=2e^{T-t}\mathbb{E}_{\xi\sim\mu} [w^{\prime}(x-\xi)]\). Therefore, with \(w(x)=\cos(x)\), and by choosing \(F\) equal to

\[F(t,x,\mu)= -\frac{P}{2(P-2e^{T-t}\mathbb{E}_{\xi\sim\mu}[\cos(x-\xi)])^{2}}+ \frac{1}{P-2e^{T-t}\mathbb{E}_{\xi\sim\mu}[\cos(x-\xi)]}+\] \[\mathbb{E}_{\xi\sim\mu}[\cos(x-\xi)]e^{T-t}(1+\frac{1}{2}\frac{1} {(P-2e^{T-t}\mathbb{E}_{\xi\sim\mu}[\cos(x-\xi)])^{2}})+\] \[e^{T-t}\mathbb{E}_{\xi\sim\mu}[(\beta(t,x,\mu)-\beta(t,\xi,\mu) )\sin(x-\xi)]+\] \[e^{T-t}\mathbb{E}_{\xi\sim\mu}[\cos(x-\xi)\frac{1}{2}\frac{1}{(P -2e^{T-t}\mathbb{E}_{\xi^{\prime}\sim\mu}[\cos(\xi-\xi^{\prime})])^{2}}],\]

the function \(V\) satisfies the Master Bellman equation, hence is the value function to the mean-field control problem.

To be easily computable, the function \(F\) can be rewritten using trigonometric relations as

\[F(t,x,\mu)= -\frac{P}{2}\frac{1}{(P-2e^{T-t}[\cos(x)\overline{\cos}_{\mu}+ \sin(x)\overline{\sin}_{\mu}])^{2}}+\frac{1}{P-2e^{T-t}[\cos(x)\overline{\cos} _{\mu}+\sin(x)\overline{\sin}_{\mu}]}\] \[+\;e^{T-t}[\cos(x)\overline{\cos}_{\mu}+\sin(x)\overline{\sin}_{ \mu}](1+\frac{1}{2}\frac{1}{(P-2e^{T-t}[\cos(x)\overline{\cos}_{\mu}+\sin(x) \overline{\sin}_{\mu}])^{2}})\] \[+\;e^{T-t}[\beta(t,x,\mu)\sin(x)\overline{\cos}_{\mu}-\beta(t,x, \mu)\cos(x)\overline{\sin}_{\mu}\] \[\quad-\;\mathbb{E}_{\xi\sim\mu}[\beta(t,\xi,\mu)\cos(\xi)]\sin(x )+\mathbb{E}_{\xi\sim\mu}[\beta(t,\xi,\mu)\sin(\xi)]\cos(x)]\] \[+\;\frac{e^{T-t}}{2}\cos(x)\mathbb{E}_{\xi\sim\mu}[\frac{\cos( \xi)}{2[P-2e^{T-t}(\cos(\xi)\overline{\cos}_{\mu}+\sin(\xi)\overline{\sin}_{ \mu})]^{2}}]\] \[+\;\frac{e^{T-t}}{2}\sin(x)\mathbb{E}_{\xi\sim\mu}[\frac{\sin( \xi)}{2[P-2e^{T-t}(\cos(\xi)\overline{\cos}_{\mu}+\sin(\xi)\overline{\sin}_{ \mu})]^{2}}]\]

with the notations: \(\overline{\cos}_{\mu}:=\mathbb{E}_{\xi\sim\mu}[\cos(\xi)]\), \(\overline{\sin}_{\mu}:=\mathbb{E}_{\xi\sim\mu}[\sin(\xi)]\).

We take \(P=2.2e^{T}\) so that the control is bounded, the same trend \(\beta\) as in (5.1.3), and parameters as in section 5.1.3: \(\kappa=\sigma=1\), \(T=0.4\), \(n=40\), \(M=20000\). At each gradient iteration, the initial distribution is sampled with

\[x_{0}\sim\;\mu_{0}=\bar{\mu}_{0}+\upsilon_{0}\mathcal{N}(0,1)\]

where \((\bar{\mu}_{0},\upsilon_{0}^{2})\) is sampled from \((0.2U([0,1]),0.5U([0,1]))\). Controlling the volatility is more difficult than controlling the trend, and it is crucial for the method that \(\rho^{A}\) is very small. We take \(\hat{N}=9000\) gradient iterations. Training time with \(L=3\) is 67228s, while it takes 69073s for \(L=4\).

On Figure 10, we give the relative error obtained with \(L=3\) and \(L=4\). Notice that with \(L=3\), \(\rho^{A}\) is small and that we have to take \(\rho^{A}\) even smaller with \(L=4\). The relative error is small, but the control is not as well approximated as in the controlled drift example, as shown in Figure 11.
## 6 Conclusion

In this study, we have presented a robust effective resolution to the challenging problem of mean-field control within a partially model-free continuous-time framework. Leveraging policy gradient techniques and actor-critic algorithms, our approach has demonstrated the valuable role of moment neural networks in the sampling of distributions. We have illustrated the significance of maintaining a low number of moments (typically two or three) while underscoring the critical role played by fine-tuning learning rates for actor and critic updates.

Subsequent developments could encompass the extension to non-separable forms within the state and control components of drift and diffusion coefficients. Furthermore, a compelling direction for further investigation could involve mean-field dynamics governed by jump diffusion processes, where the intensities of the jumps remain unknown.
