Numerical resolution of McKean-Vlasov FBSDEs using neural networks
Maximilien Germain Joseph Mikael Xavier Warin


The spin-forbidden vacuum-ultraviolet absorption spectrum of ^14N^15N
    W. Ubachs
    December 30, 2023
=====================================================================




We propose several algorithms to solve McKean-Vlasov Forward Backward Stochastic Differential Equations (FBSDEs).
Our schemes rely on the approximating power of neural networks to estimate the solution or its gradient through minimization problems. As a consequence, we obtain methods able to tackle both mean-field games and mean-field control problems in moderate dimension. We analyze the numerical behavior of our algorithms on several examples including non linear quadratic models.





Key words: Neural networks,  McKean-Vlasov FBSDEs, Deep BSDE, mean-field games, machine learning.  




MSC Classification:  65C30, 68T07, 49N80, 35Q89.



§ INTRODUCTION



  

This paper is dedicated to the numerical resolution in moderate dimension of the following McKean-Vlasov Forward Backward Stochastic Differential Equations (MKV FBSDEs)

    X_t     = ξ + ∫_0^t b(s,X_s,Y_s,Z_s,ℒ(X_s),ℒ(Y_s),ℒ(Z_s))  s + ∫_0^t σ(s,X_s,ℒ(X_s))  W_s 
    
    Y_t     = g(X_T,ℒ(X_T))  + ∫_t^T f(s,X_s,Y_s,Z_s,ℒ(X_s),ℒ(Y_s),ℒ(Z_s))  s - ∫_t^T Z_s  W_s

with b: ×^d×^k×^k× d×𝒫_2(^d)×𝒫_2(^k)×𝒫_2(^k× d) ↦^d, σ: ×^d ×𝒫_2(^d) ↦^d, g: ^d×𝒫_2(^d) ↦^k, and f: ×^d×^k×^k× d×𝒫_2(^d)×𝒫_2(^k)×𝒫_2(^k× d) ↦^k.  

W_t is a d-dimensional ℱ_t-Brownian motion where (Ω,𝒜,ℱ_t,) is a given filtered probability space and T>0.

ξ is a given random variable in L^2(Ω,ℱ,;^d) and
𝒫_2(^n) stands for the space of square integrable probability measures over ^n endowed with the 2-Wasserstein distance 

    𝒲_2(μ,ν) = inf{√([(X-X')^2]) | X,X'∈ L^2(Ω,ℱ,;^n), ℒ(X) = μ, ℒ(X') = ν}.

At last ℒ(.) is a generic notation for the law of a random variable.


  

This kind of equation is linked to non local PDEs kwown as master equations. We refer to  <cit.>, chapter 4 and 5 of volume 2, for an introduction on the subject.
In <cit.>, it is shown for example that under regularity conditions, when the drift b is independent of  time and the law ℒ(Z_t), the driver f does not depend on time and the law ℒ(Z_t), σ does not depend on time, then the resolution of equation (<ref>) provides a way to estimate the solution of the equation

    ∂_t 𝒰(t,x, μ) +     b(x, 𝒰(t,x, μ),∂_x 𝒰(t,x, μ) σ(x, μ),μ, η).∂_x 𝒰(t,x, μ)  
        + 1/2[ ∂^2_xx𝒰(t,x, μ) σσ^⊤(x, μ)] + f(x, 𝒰(t,x, μ), ∂_x 𝒰(t,x, μ) σ(x, μ),μ, η)  
        + ∫_^d∂_μ𝒰(t,x, μ)(y). b(y,𝒰(t,y, μ),∂_x 𝒰(t,x, μ) σ(x, μ),μ, η)  μ(y)  
       
     + ∫_^d1/2[ ∂_x ∂_μ𝒰(t,x, μ)(y) 
    σσ^⊤(y, μ)] μ(y) =0 on [0,T]×^d×_2(^d)
 with initial condition 𝒰(0,x, μ) = g(x,μ) on ^d×_2(^d), and
where  η is a notation for the image of the probability measure μ by the mapping x∈^d ⟶𝒰(t,x,μ) ∈^k.  
Under suitable assumptions, <cit.> proves that the solution to (<ref>) admits the so-called decoupling field representation Y_s = U(s,X_s,(X_s)) where (t,x,μ)∈[0,T]×^d×_2(^d)↦ U(T-t,x,μ) is a classical solution to (<ref>). See Theorem 2.9 and equation (2.12) in <cit.>. Their results are stated in the more general of dynamics depending in the joint law of (X_t,Y_t) but we state them here with the particular case of dependence in the marginal laws (X_t),(Y_t) to be coherent with (<ref>). This MKV FBSDE representation is used in <cit.> to build a numerical scheme for the approximation of (<ref>) when the drift b is independent of Z.
The equation above is non local due to the integral terms
 and the term ∂_μ𝒰(t,x, μ)(y) stands
 for the Wasserstein derivative of 𝒰 in the direction of the measure at point (t,x, μ) and evaluated at the continuous coordinate y. 



  

Equations (<ref>)   appear as well  as probabilistic formulations of mean-field games or mean-field controls characterizing the value function V of the game.
Mean-field games are introduced by <cit.> and <cit.> to model games with interactions between many similar players. In this theory, each player's 
dynamics and cost take into account the empirical distribution of all agents. At the limit of an infinite number of players, the search for a Nash equilibrium with close loop controls boils down to a control problem concerning a representative player whose law enters in the cost and dynamics.

Two probabilistic approaches based on Forward Backward Stochastic Differential Equations can be used to solve these problems:

    
  * A first approach  called the Pontryagin approach consists as shown in <cit.> in applying  the strong Pontryagin principle to these control problems. Under regularity and convexity conditions,  Y_t appears to be  a stochastic representation of the gradient of the value function V. In this case the coefficients b,f of the related MKV FBSDE (<ref>) do not depend on Z_t, ℒ(Y_t),ℒ(Z_t) and k = d. 
    
  * Another approach called the Weak approach permits to solve the optimization problem by estimating directly Y_t as the value function V of the problem as shown in  <cit.>. In this case the coefficients b,f of the related MKV FBSDE (<ref>) do not depend on Y_t, ℒ(Y_t),ℒ(Z_t) and k = 1. 



  

The numerical resolution of equations (<ref>) is rather difficult since:

    
  * The dynamics are coupled through both the drift and the driver of the BSDE.
    
  * The McKean-Vlasov structure of the problem requires to solve a fixed point in probability spaces. 

In the linear-quadratic setting (with quadratic cost to minimize but linear dynamics) the weak approach applied to mean-field problems gives a problem in low dimension but with quadratic coupling in Z_t appearing in the backward dynamic. In contrast, the Pontryagin approach exhibits a problem in potentially high dimension (the Z component is a d× d matrix in this case) but with a linear coupling in Y_t which is easier to solve numerically.

In the case of mean-field games, only the law of X_t is present in the dynamic of (<ref>).
In mean-field games of controls (also called extended mean-field games, see <cit.> and the references therein), individuals interact through their controls instead of their states as in the model of trade crowding in <cit.>. The law of the control thus appears in the dynamic of (<ref>) and may give rise to some FBSDE depending on  the law of Z_t in the weak approach or the law of Y_t in the Pontryagin approach.


  

Existence and uniqueness of a solution to the fully coupled system (<ref>) are studied by <cit.> when the drift does not depend on the law of Z. Their Theorem 4.29 gives existence of a solution under a non-degeneracy condition. However, uniqueness is a priori only expected to hold in small time, as stated in Theorem 4.24 from <cit.>.
In the latter we will assume that existence and uniqueness hold for the MKV FBSDE we aim to solve.



  

In <cit.> and <cit.>, tree and grid algorithms are proposed and tested in dimension 1. It is worth mentioning that these techniques suffer from the so-called curse of dimensionality and cannot be applied when the dimension describing a player state is high (typically greater than 3 or 4). This is due to the discretization of the state space.


  

However, new approaches using machine learning are developed since 2017 to solve non linear parabolic PDEs through a BSDE representation.
Two kinds of methods have emerged:

    
  * The first to appear are global methods first proposed in <cit.> to solve semi-linear PDEs. They rely on a single high-dimensional optimization problem whose resolution is difficult. It consists in the training of as many neural networks as time steps by solving in a forward way the backward representation of the PDE solution.  The Z_t process is represented by a different neural network Z^θ_i with parameters θ at each date t_i. Instead of solving the BSDE starting from the terminal condition, the method writes it down as a forward equation and an optimization problem aiming to reach the terminal condition g(X_T) by minimizing a mean squared error |Y_T - g(X_T)|^2. The approach is extended to fully nonlinear equations (nonlinear in the solution, its gradient and hessian) in <cit.> and the authors show that the methodology can solve some equations in high dimension.
    <cit.> showed that it is more effective to used a single network for all dates and besides proposed an original fixed point algorithm to solve semilinear PDEs.
    
  * A second kind of local methods first proposed in <cit.> is based on local optimization problems solved at each time step in a backward way.  Contrarily to the global method, the successive optimization problems are here in moderate dimension. Each optimization step at date t_i consists in the training of only two local neural networks Y^θ_i, Z^θ_i with parameters θ. For instance, instead of solving 1 optimization problem with N neural networks in the global method, <cit.> solves N learning problems with 2 neural networks. Moreover the resolution is simplified by the initialization of the neural networks at time t_i to their previously computed values at time t_i+1, which provides a good approximation for the current value. This strategy is inspired by the standard backward resolution of BSDE with conditional expectations from <cit.> and <cit.>. 
    The methodology is extended to the much more challenging case of fully nonlinear PDEs in <cit.> by combining it with some ideas proposed in <cit.>. Extensive tests performed in <cit.> show that the local method gives better results than the global one, such as <cit.> in the case of fully nonlinear dynamics. Especially, these papers show that local methods can be used with a larger time horizon T than the global method.



  

Machine learning techniques to solve coupled FBSDEs are investigated by several authors in <cit.> and <cit.>, and a first method for McKean-Vlasov FBSDEs with delay is studied by <cit.> for a linear quadratic equation in dimension one. Similar and more general ideas are presented and tested in dimension one in <cit.> alongside convergence results, together with an additional method directly solving mean-field control problems by minimizing the cost without writing down optimality conditions.
The resulting algorithms proposed all rely on the global approach first initiated in <cit.>.


  

Our paper aims to extend these methods and to propose new ones for the resolution of McKean-Vlasov FBSDEs in moderate dimension, and go beyond one dimensional examples for which standard methods are already available (see <cit.>). In fact, one major advantage for the use of neural networks for solving control problems is their ability to efficiently represent high-dimensional functions without using space grids. We also study the influence of the maturity T on the algorithms.
We first propose to modify the previously proposed algorithm to
stabilize its convergence. 
Our modification allows us to reduce the variance of the estimators used in the dynamic of X_t and Y_t.
Then we propose a second algorithm relaxing the fixed point iteration algorithm by adding a neural network learning the distribution of the solution thanks to a penalization in the loss function.
At last we propose a resolution scheme based on some local resolution as in <cit.>.

  

To simplify the presentation, we consider first order interaction models, that is to say that the dependency of the drift and cost function with respect to the laws ℒ(X_t),ℒ(Y_t),ℒ(Z_t) only concerns expectations in the form

    u_t = (u_t^X,u_t^Y,u_t^Z) := ([φ_1(X_t)],  [φ_2(Y_t)], [φ_3(Z_t)]),
 for some continuous functions φ_1,φ_2,φ_3 with adequate domains and codomains. In this framework we can rewrite by abuse of notation 

    b(s,X_s,Y_s,Z_s,ℒ(X_s),ℒ(Y_s),ℒ(Z_s))     = b(s,X_s,Y_s,Z_s,u_s)
    σ(s,X_s,ℒ(X_s))     = σ(s,X_s,u_s^X) 
    
    g(X_T,ℒ(X_T))     = g(X_T,u_T^X)
    
    f(s,X_s,Y_s,Z_s,ℒ(X_s),ℒ(Y_s),ℒ(Z_s))     = f(s,X_s,Y_s,Z_s,u_s).

For instance when φ_1,φ_2,φ_3 are power functions with positive integers as exponents we recover probability distribution moments.
See Remark <ref> for more general cases beyond first order interaction.


  

We provide multidimensional tests to show how these machine learning approaches can overcome the curse of dimensionality on some test cases first coming from a mean-field game of controls: we solve the FBSDE derived from the weak approach and the Pontryagin approach. We also consider an example arising from a non linear quadratic mean-field game.
Then we compare all the methods on some general test cases of FBSDE involving linear or quadratic dependence on the processes X_t, Y_t, Z_t and on their distributions.


  

The structure of the paper is the following: in sections <ref> and <ref> we describe the proposed schemes, and in section <ref> we provide a numerical study of our methods in dimension 10 (except for the one-dimensional Example of Section <ref>). We show that our algorithms can solve non linear-quadratic models with small maturities.


§ MACHINE LEARNING GLOBAL SOLVERS

In this section we propose three global   algorithms based on the approach in <cit.>. 



 §.§ Algorithm principle




  

We propose a generalized and refined version of the Algorithm 2 from <cit.>. We recall that a similar technique with additional networks is used in <cit.> for delayed McKean-Vlasov equations but is tested only on a one dimensional linear quadratic example. Our methods also take advantage of different expectation computation methods, introduced in section <ref>. 
We present in section <ref> several tests in dimension 10 where the laws of X, Y, Z are involved. 


  

We consider the Euler-Maruyama discretized FBSDE system (<ref>) on a regular time grid t_k = kT/N for k∈0,N:

    X_t_i+1    = X_t_i + b(t_i,X_t_i,Y_t_i,Z_t_i,u_t_i) Δ t+ σ(t_i,X_t_i,u_t_i^X ) Δ W_i
    
    Y_t_i+1    = Y_t_i - f(t_i,X_t_i,Y_t_i,Z_t_i,u_t_i) Δ t + Z_t_iΔ W_i,

with terminal condition Y_t_N = g(X_t_N,u_t_N^X) and initial condition X_0 = ξ. We recall that u_t is defined in (<ref>). We note Δ t := t_i+1 - t_i = T/N and (Δ W_i)_i=0,⋯,N-1 := (W_t_i+1 - W_t_i)_i=0,⋯,N-1 the Brownian increments.
In the FBSDE theory, one requires the processes (X_t_i,Y_t_i,Z_t_i) to be ℱ_t_i-adapted. Therefore the backward part of the system can also be written in the conditional expectation form

    Y_t_i = [Y_t_i+1 + f(t_i,X_t_i,Y_t_i,Z_t_i,u_t_i) Δ t | ℱ_t_i]
    
           Z_t_i =  [Y_t_i+1Δ W_i/Δ t | ℱ_t_i],
 where we see how the process Z is defined. This process, specific to the stochastic case, allows the Y component to be ℱ_t-adapted, even though we fix its terminal condition. It is a major difference between backward ordinary differential equations and backward stochastic differential equations.
  We see that the whole system is coupled therefore we need to design a method allowing to solve simultaneously both equations of (<ref>).


  

We solve the system by the Merged Deep BSDE method introduced in <cit.>. Z_t_i is approximated by a single feedforward neural network _θ^z(t_i, X_t_i) and  Y_0 by a neural network _θ^y(X_0) with parameters θ = (θ^y,θ^z). 
With this point of view, the discretized Brownian motion W_t acts as training data in the language of machine learning, so that we can generate a training set as large as desired. 
Extensive tests conducted in <cit.> show that the use of a Merged network improves the training in comparison with the Deep BSDE method of <cit.>. Indeed it lowers the number of parameters to learn hence reduces the complexity of the problem. It empirically improves the accuracy of the method but also makes the training faster. That is why we focus on this architecture. It is also used by <cit.> which considers controls in the form (ζ(t_i,X_i))_i=1,…,N for a neural network ζ. The use of  recurrent networks such as Long Short Term Memory networks as in <cit.> is possible  but tests achieved in <cit.> seem to show that is does not bring more accuracy on Markovian problems. Other alternatives may include the GroupSort network <cit.> for a better control of the  Lipschitz constant of the approximation, or some special networks preserving some properties of the solution but they have not been tested.

 The motivation for such an approximation comes from the notion of decoupling field, also used for numerical purposes in <cit.> or <cit.>, which gives the existence of functions u,v (see the paragraph below (<ref>)) such that

    Y_t = u(t,X_t,ℒ(X_t)), Z_t = v(t,X_t,ℒ(X_t)).

Numerically, it is enough to consider Y and Z as a function of the couple (t,X_t). In fact, the law of the solution (and therefore its moments) can be seen as a function of t. That's why we search for a representation 

    Y_t = u(t,X_t), Z_t = v(t,X_t).
 The forward-backward system is transformed into a forward system and an optimization problem aiming to satisfy the terminal condition of the BSDE through the loss function [(Y_T - g(X_T, u^X_t_N)
 )^2]. To simplify notations, X_i := X_t_i and similarly for Y and Z.



  

In practice the loss function is minimized with the Adam gradient descent method <cit.>. 
In any case, the goal of our scheme is to learn both the optimal control and the distribution of X_t, Y_t, Z_t. In the following, B is the batch size, that is the number of particles we will sample to approximate the loss function in expectation form, see for instance (<ref>). N is the number of time steps and M is the number of previous batches expectations to keep in memory. 


  

We use for  and  feedforward neural networks with 3 hidden layers (d + 10 neurons in each) with hyperbolic tangent function as activation functions and an output layer with identity as activation. It is worth noticing that because the merged neural network takes the couple (t,X) as inputs, we cannot use batch normalization since the distribution of X_i is not stationary over time.



 §.§ Estimation of the expectation



  

A key step for the methods is to estimate the mean-field parameter u. It has a significant effect on the algorithms performances. We note θ_m=(θ^y_m,θ^z_m) the neural network parameters at optimization iteration m and u_i = (u_i^X,u_i^Y,u_i^Z) the estimation of 
u_t_i. 
In the algorithms described below, the approximated processes are considered as functions of the parameters θ of the neural network.

Several methods can be used to approximate the moments of the solution involved in the stochastic McKean-Vlasov dynamics:

    
  * Direct: use the empirical mean of the current batch of particles
    
    u_i = 1/B(∑_j=1^Bφ_1(X^j_i(θ_m)), ∑_j=1^Bφ_2(Y^j_i(θ_m)), ∑_j=1^Bφ_3(Z^j_i(θ_m))), i=0, ⋯ , N-1.

    
    Alternatively one could use instead the last batch particles
    to estimate the law 
    
    u_i = 1/B(∑_j=1^Bφ_1(X^j_i(θ_m-1)), ∑_j=1^Bφ_2(Y^j_i(θ_m-1)), ∑_j=1^Bφ_3(Z^j_i(θ_m-1))), i=0, ⋯ , N-1.
 The difference lies in the fact that in one case the optimization of the parameters θ_m at iteration m modifies the current estimation of the law whereas using the previously computed parameters θ_m-1 fixes the law and simplifies the optimization problem.
     In practice, for the numerical tests of Section <ref> we use the formula (<ref>). This approach requires to handle very large batches, typically of the order of B=10,000 sample paths get a reasonable approximation of the laws. This is the approach used by <cit.> and <cit.>. 
    
    We solve the following optimization problem
    
    min_θ = (θ^y,θ^z)    1/B∑_k=1^B |Y_N^k(θ) - g(X_N^k(θ),1/B∑_j=1^B φ_1(X_N^j(θ))) |^2
    
             X_i+1^j(θ) =     X_i^j(θ) + b(t_i,X_i^j(θ),Y_i^j(θ),_θ^z(t_i,X_i^j(θ)) ,u_i) Δ t 
             +  σ(t_i,X_i^j(θ) ,u_i^X) Δ W_i^j
    
             Y_i+1^j(θ) =     Y_i^j(θ) - f(t_i,X_i^j(θ),Y_i^j(θ),_θ^z(t_i,X_i^j (θ)) ,u_i) Δ t 
             +  _θ^z(t_i,X_i^j(θ)) Δ W_i^j 
    
             u_i =    1/B(∑_j=1^Bφ_1(X^j_i(θ)), ∑_j=1^Bφ_2(Y^j_i(θ)), ∑_j=1^Bφ_3(_θ^z(t_i,X_i^j(θ))))
    
             X_0^j =    ξ^j ∼ξ, j=1,⋯,B 
    
             Y_0^j(θ) =    _θ^y(X_0^j), 
             
              i =      1,⋯,N-1.

    
    The Global Direct solver leads to algorithm <ref>. 
    

    
  * Dynamic: a method which dynamically updates the estimation on (M+1) B samples. The expectations from the last M batches are kept in memory in an array 
    (ζ_i,r)_[   i=0, …, N-1,; r= 0, … ,  M-1 ]
 
    initialized with values ([φ_1(ξ)], φ_2(0),φ_3(0))^N× M.

    At iteration m - 1, ν_i^(m-1) is defined as the empirical mean on these previous sample paths. On a new batch, the expectation is computed by averaging the previous estimation ν_i^(m-1) and the current batch empirical mean  by the following algorithm  used for i=0, ⋯ , N-1:
    
    ν_i^(m-1)   = 1/M∑_r=0^M-1ζ_i,r, 
    
             u_i    = M ν_i^(m-1) + 1/B(∑_j=1^Bφ_1(X^j_i(θ_m)),∑_j=1^Bφ_2(Y^j_i(θ_m)), ∑_j=1^Bφ_3(Z^j_i(θ_m)))/M + 1, 
    ζ_i,m% M    = 1/B(∑_j=1^Bφ_1(X^j_i(θ_m)),∑_j=1^Bφ_2(Y^j_i(θ_m)), ∑_j=1^Bφ_3(Z^j_i(θ_m))).
 The notation m%M refers to the remainder of the Euclidian division of m by M. This technique allows to use smaller batches of size 100 or 1000. Thus it is more efficient in terms of convergence speed in comparison with the direct approach. This method can be seen as a dynamic fixed point approach. 
    

  

    The idea behind this update rule comes from online learning in machine learning. 1/M can be interpreted as a learning rate quantifying the updating speed. From the current estimation of the particles law, we introduce a small correction related to the new observed samples. Therefore the estimation is much more stable through iterations compared to the instantaneous update of the law used by the Direct method. After M batches, the older samples are forgotten, since they don't represent anymore the current law. Indeed we expect the convergence for a good choice of M. If this parameter is too small the stabilization would be inefficient and on the contrary a too large M would slow down the learning process by introducing a bias in the law.
    For instance in our numerical experiments of Section <ref> we use M = 100 for a total of 2000 gradient descent iterations.
    
    
    	If the law dependence is more general than a first order interaction 
    	and is given by a continuous function F: μ∈𝒫_2(^d) ↦^k then the Direct method can be straightforwardly applied to the equation by estimating F(ℒ(X_t)) by the so-called empirical projection F(1/B∑_j=1^B δ_X^j_t_i) for identically distributed particles (X^j_t_i)_j=1,…,B on a time grid t_0,⋯,t_N. Concerning the Dynamic approach, it would require to keep in memory the previously computed particles from the last M batches which is costly.
    
    
    
    The fixed point approach is known to be convergent theoretically only for small maturities. In practice, the theoretical bound on the maturity found on the simple example given for example in paragraph 3.1 in <cit.> is far too pessimistic.  We will see that the restriction is not relevant on all our test cases.
    
    
   
    For a given iteration m, given the estimations (ζ_i,r)_[   i=0, …, N-1,; r= 0, … ,  M-1 ] of u_i on the last M iterations, we perform one gradient descent step for the following optimization problem
    
    min_θ_m=(θ^y_m,θ^z_m)    1/B∑_k=1^B |Y_N^k(θ_m) - g(X_N^k(θ_m),1/B∑_j=1^B φ_1(X_N^j(θ_m))) |^2
    
             X_i+1^j(θ_m) =     X_i^j(θ_m) + b(t_i,X_i^j(θ_m),Y_i^j(θ_m),_θ^z_m(t_i,X_i^j(θ_m)),u_i) Δ t 
        +  σ(t_i,X_i^j(θ_m),u_i^X(θ_m)) Δ W_i^j
    
             Y_i+1^j(θ_m) =     Y_i^j(θ_m) - f(t_i,X_i^j(θ_m),Y_i^j(θ_m),_θ^z_m(t_i,X_i^j (θ_m)),u_i) Δ t 
        +  _θ^z_m(t_i,X_i^j(θ_m)) Δ W_i^j 
    
             u_i =    1/B(∑_j=1^Bφ_1(X^j_i(θ_m)), ∑_j=1^Bφ_2(Y^j_i(θ_m)), ∑_j=1^Bφ_3(_θ^z_m(t_i,X_i^j(θ_m))))
    u_i =    ∑_r=0^M-1ζ_i,r + u_i/M + 1
    
             X_0^j =    ξ^j ∼ξ, j=1,⋯,B
    
             Y_0^j(θ_m) =    _θ^y_m(X_0^j)
    
             i =      1,⋯,N-1.
 Then we update (ξ_i)_i by forgetting the oldest estimation and keeping in memory the new one,  (u_i)_i (see (<ref>)). The Global Dynamic solver is given more explicitly in algorithm <ref>.


   
    
  * Expectation: estimate u_t by a neural network Ψ_θ^Ψ with input t and parameters θ^Ψ. 
    
    u_i (θ^Ψ) = Ψ_θ^Ψ(t_i) = ( Ψ_θ^Ψ^X(t_i),Ψ_θ^Ψ^Y(t_i) ,Ψ_θ^Ψ^Z(t_i) ), i=0, ⋯ , N.

    A penalization term
    
    [λ/N∑_i=0^N-1‖Ψ_θ^Ψ(t_i) - 1/B( ∑_j=1^Bφ_1(X^j_i(θ)), ∑_j=1^Bφ_2(Y^j_i(θ)), ∑_j=1^Bφ_3(Z^j_i(θ)))‖_2^2],

    is added to the loss function. 
    We will see that in practice this method is quite involved to use because the performances heavily depend upon the choice of the parameter λ. This approach provides a relaxation of the fixed point method.
    
    We solve the following optimization problem
    
    min_θ = (θ^y,θ^z,θ^Ψ)    1/B∑_k=1^B |Y_N^k(θ) - g(X_N^k(θ),1/B∑_j=1^B φ_1(X_N^j(θ))) |^2
        + λ/N∑_i=0^N-1‖ u_i(θ^Ψ) - Ψ_θ^Ψ(t_i)‖^2
    
             X_i+1^j(θ) =     X_i^j(θ) + b(t_i,X_i^j(θ),Y_i^j(θ),_θ^z(t_i,X_i^j(θ)),Ψ_θ^Ψ(t_i)) Δ t 
        +  σ(t_i,X_i^j(θ),Ψ_θ^Ψ(t_i)^X) Δ W_i^j
    
             Y_i+1^j(θ) =     Y_i^j(θ) - f(t_i,X_i^j(θ),Y_i^j(θ),_θ^z(t_i,X_i^j (θ)),Ψ_θ^Ψ(t_i)) Δ t 
        +  _θ^z(t_i,X_i^j(θ)) Δ W_i^j 
    
             u_i =    1/B(∑_j=1^Bφ_1(X^j_i(θ)), ∑_j=1^Bφ_2(Y^j_i(θ)), ∑_j=1^Bφ_3(_θ^z(t_i,X_i^j(θ))))
    
             X_0^j =    ξ^j ∼ξ, j=1,⋯,B
    
             Y_0^j(θ) =    _θ^y(X_0^j)
    
             i =      1,⋯,N-1.

    The Global Expectation solver is described in algorithm <ref>. 
    The parameter λ is chosen by trial and error.
    

We will compare the performances of these techniques on several examples in section <ref>.



§ A LOCAL SOLVER


We also propose a local method inspired by the Deep Backward Dynamic Programming introduced by <cit.> and <cit.>. It considers local minimization problems between contiguous time steps. In this case there are as many networks as time steps. We replace a global optimization setting by a set of smaller problems.


  

In this method for i∈ 0, N-1, Z_i and Y_i are approximated by a neural network (^i_θ_i^z(·) ,^i_θ_i^y(·)) with parameters θ = (θ_0^y,θ_0^z,⋯,θ_N-1^y, θ_N-1^z). At iteration m, with θ_m = (θ_m,0^y , θ_m,0^z,⋯, θ_m,N-1^y, θ_m,N-1^z) , we simulate X_i(θ_m) with the previously computed parameters θ_m. 

    X_i+1^j(θ_m) =     X_i^j(θ_m) + b(t_i,X_i^j(θ_m),^i_θ_m,i^y(X_i^j(θ_m)),^i_θ_m,i^z(X_i^j(θ_m)),u_i) Δ t 
     +    σ(t_i,X_i^j(θ_m),u_i^X) Δ W_i^j.

This first step allows to find the areas visited by the controlled process. Using R samples, we compute the empirical mean m_i^m = 1/R( ∑_j=1^R X^j_i(θ_m)) and variance V_i^m = 1/R∑_j=1^R(X^j_i(θ_m))^2 - 1/R^2(∑_j=1^R X^j_i(θ_m))^2 of X_i(θ_m). We estimate u_t as in the Dynamic method (<ref>): 

    u_i =    1/R(∑_j=1^Rφ_1(X^j_i(θ_m)), ∑_j=1^Rφ_2(^i_θ_m,i^y(t_i,X_i^j(θ_m))), ∑_j=1^Rφ_3(^i_θ_m,i^z(t_i,X_i^j(θ_m))))
    u_i  =    ∑_r=0^M-1ζ_i,r + u_i/M + 1.

A priori R can be different from the batch size B. It is the batch size used for the estimation of the law with the previously computed parameters. In the numerical tests of Section <ref> we use B=100 or B=300 for the backward optimization and a larger value R= 50000 for the Monte-Carlo forward estimation of the law. Then we solve backward problems to find the θ_m+1 by sampling B independent copies of X_i through a Gaussian distribution 𝒩( m_i^m,V_i^m ) with frozen parameters θ_m : 

    
  * First sample B independent copies X^1_N,⋯,X^B_N of X_N following a Gaussian distribution 𝒩( m_N^m,V_N^m). Y^N_θ_m+1,N^y(X_N^j) is set to the terminal condition g(X_N^j ,u_N^X).

  * For i from N-1 to 0:

Sample B independent copies X^1_i,⋯,X^B_i of X_i following a Gaussian distribution 𝒩(m_i^m, V_i^m). Diffuse according to the dynamics (<ref>) of X, starting from X^1_i,⋯,X^B_i, to obtain X^1_i+1,⋯,X^B_i+1.

Solve the local optimization problem 

    min_θ=(θ^y,θ^z)1/B∑_j=1^B    |Y^i+1_θ_m+1,i+1^y(X_i+1^j) - Y^i_θ^y(X_i^j)+ f(t_i,X_i^j,^i_θ^y(X_i^j),^i_θ^z(X_i^j),u_i) Δ t . 
       . - ^i_θ^z(X_i^j)Δ W_i^j |^2,
starting from the parameter value θ_m+1,i+1. We can then update the θ value by denoting as  θ_m+1,i the argmin value of the minimization problem.

  * Repeat the previous steps for the iteration m+1 until reaching K iterations.


In the version of the Local Dynamic solver given in algorithm <ref>, we use the dynamic update of the expectations introduced previously in the dynamic solver of section <ref>. In this algorithm H stands for the number of gradient steps to perform at each step of the algorithm and  R is the number of samples for the laws estimation.


Because we have to learn the dynamic of the forward process, the use of a backward resolution is not as obvious as in <cit.>. We have to alternate between forward dynamic estimations and backward resolutions. 

More precisely, here we solve fully coupled FBSDEs, whereas the works <cit.> consider decoupled FBSDEs, where the forward process X can be simulated independently of Y,Z. Estimating the law of X and sampling from a normal distribution allows us to decouple and solve locally the FBSDEs in areas visited by X and its law. However, because of this freezing of the forward dynamics, another fixed point problem has to be solved. Other approaches for fully coupled FBSDEs like <cit.> and <cit.> rely on the global machine learning method initiated by <cit.>.

Notice that we estimate the law dependent part of the dynamics in the first loop but for the resolution loop, we use a normal distribution as training measure. 

Another natural training measure could be the empirical law of the process. Here, we employ an hybrid method between the two methods proposed in Remark 2.1 from <cit.>. We exploit the approximated control to estimate the law of X and use a Gaussian law with the same two first moments as X as training measure. In the numerical examples of <cit.>, the authors also use Gaussian distributions as training measures, but with fixed parameters through the iterations. Here we adapt the training measures during the training to have the same expectation and variance as the controlled process thanks to an exploration step.







§ NUMERICAL RESULTS

The algorithms are implemented in Python with the Tensorflow library <cit.>. Each numerical experiment is conducted using a node composed of 2 Intel® Xeon® Gold 5122 Processors, 192 Go of RAM, and 2 GPU nVidia® Tesla® V100 16Go. The multi-GPU parallelization on the global solver is conducted using the Horovod library <cit.>.
The methods we test are:


    
  * Global Direct: algorithm <ref> at page algo: direct. Batch size B = 10000.
    
  * Global Dynamic: algorithm <ref> at page algo: dynamic. Batch size B = 200 and M = 100.
    
  * Global Expectation: algorithm <ref> at page algo: expectation. Batch size B = 2000.
    
  * Local Dynamic: algorithm <ref> at page algo: local. Batch size B = 300 (Weak), B = 100 (Pontryagin) and M = 20, R=50000.

The name depends on the two choices made regarding the algorithm: the use a Global or Local method but also what law estimation technique is applied.
If the algorithm is applied to equations coming from the Pontryagin (abbreviated in Pont.) or the Weak approach, it is specified in its name. In the following tables, the incorrect results are highlighted in red. In Section <ref> the three best results are highlighted in green whereas in Section <ref> only the best value is.



 §.§ Linear price impact model


We use a linear-quadratic mean-field game of controls model studied in <cit.> and <cit.> for comparison. This model is useful for numerical tests since the analytic solution is known. The MFG of controls model for the representative player is given by:

    min_α∈𝔸      [∫_0^T (c_α/2α_t^2 + c_X/2X_t^2 - γ X_t · u_t )  t + c_g/2X_T^2 ]
       subject to      X_t = x_0 + ∫_0^t α_s  s + σ W_t
    .
 
and the fixed point [α_t] =u_t. In this case, the mean-field interaction is exerted through the law of the control process.

The Pontryagin optimality principle gives the system:

    X_t    = - 1/c_αY_t  t + σ  W_t
    
    X_0    = x_0
     Y_t    = - (c_X X_t + γ/c_α[Y_t])  t + Z_t  W_t
    
    Y_T    = c_g X_T.

In this case, the output Z of the neural network is a matrix of size d× d and Y is a vector of size d. 

The weak representation of the value function gives:

    X_t    = - 1/c_ασ^-1Z_t  t + σ  W_t
    
    X_0    = x_0
     Y_t    = - (c_X/2X_t^2 + γ/c_αX_t·σ^-1[Z_t] + 1/2c_ασ^-1 Z_t^2)  t + Z_t  W_t
    
    Y_T    = c_g/2X_T^2.

 In this case, the output Z of the neural network is a vector of size d and Y is a scalar. Therefore we may be able work in higher dimensions.

With LQ models, the dynamics of Y is linear in the Pontryagin approach and quadratic in the Weak approach. Thus the potentially high dimension of one method is counterbalanced by the complex dynamics of the other technique.


For our numerical experiments we take c_X = 2, x_0 = 1, σ = 0.7, γ =2, c_α = 2/3, c_g = 0.3. If not stated otherwise, the simulations are conducted with T = 1, d = 10, Δ t = 0.01. 

















  

All methods except Global Expectation Weak and Local Dynamic Weak converge to the exact solution for small maturities.
These two solvers do not converge to the right solution for any time horizon.

The choice of the parameter λ influences a lot the output of the Global Expectation Pontryagin scheme, as observed in Table <ref>. The best results are obtained when λ is of order 10 but we notice that a large range of values seems to work fine for very small time horizon T. However the Global Expectation Weak scheme never works. 
We will not test the expectation methods on the other test cases since they are less efficient than the other methods.

We see in Figure <ref> that the Local Dynamic Pontryagin  method needs more iterations for the loss to stabilize than the Global methods. We cannot hope for more iterations to help the convergence in the Global Weak methods since the loss in the learning curves of Figure <ref> reaches a plateau. The algorithms solving the system coming from the Pontryagin principle perform better than the others. The dynamic estimation of the expectation allows to gain training speed and to smooth the loss, as seen in Figure <ref> and Table <ref>. As another accuracy test, we can also plot the optimal control for which we have an analytical expression. We see in Figures <ref>, <ref> that Global Pontryagin and Local Dynamic Pontryagin methods perform well but that the Local Dynamic Weak method does not seem to converge, which confirms what is observed in Table <ref>.



 §.§ A one-dimensional mixed model

We consider the following one-dimensional example from <cit.>:

    X_t     = - ρ Y_t  t + σ  W_t, X_0 = x_0
     Y_t     = arctan([X_t])  t + Z_t  W_t, Y_T = arctan(X_T).
 This model comes from the Pontryagin principle applied to the mean-field game problem

    min_α    [∫_0^T (1/2ρα_s^2 - X_s arctan( u_s))  s + g(X_T)]
        X_t  = α_t  t +  W_t, X_0 = x_0,
 with the fixed point u_s = [X_s], and where g is an antiderivative of arctan. 
We take the same model parameters as in <cit.> (T=1 and x_0=1) and obtain in Figure <ref> with all our methods the same results as in their Figure 4. For the numerical resolution we choose 100 time steps. Notice that we use 3 hidden layers with 11 neurons in each when <cit.> uses 100 neurons by layer. We see that our smaller number of neurons is enough for this example resolution.

The duration of the methods are given in Table <ref> which illustrates again the speed gain in using the dynamic method.





 §.§ Beyond the mean-field games case



  

In this Section we design non Linear Quadratic models in order to test the limitations of our methods. We construct general MKV FBSDES with explicit solutions following a log-normal distribution. Let X_t^i be defined by


    X_t^i     = a X^i_t  t + σ X^i_t  W^i_t, 
    
      X^i_0     = ξ,
 with a, σ, ξ∈. We obtain explicitly

    X^i_t     = ξ e^(a-σ^2/2) t +σ W^i_t, 
    
       h^i_t     :=  [ X^i_t] = ξ e^a t,  
    
        k^i_t     := [ (X^i_t)^2]  = ξ e^(2a +σ^2) t .


We choose  u: (t,x) ↦ e^α tlog( ∏_i=1^n x^i )  and the following dynamic for Y_t

    Y_t     = e^α tlog(∏_i X^i_t) = e^α t∑_i [ log(ξ)+ (a-σ^2/2) t + σ W^i_t ],

  such that
  
    c_t     :=[ Y_t ]  = e^α t∑_i [log(ξ)+ (a-σ^2/2) t],
    
      d_t     := [Y_t^2]  = e^2 α t( [∑_i  ( log(ξ)+ a-σ^2/2) t]^2 + ∑_i σ^2 t ).

As we have Y_t = u(t,X_t), we obtain Z_t^i = σ X^i_t ∂_x u(t,X_t) following:

    Z_t^i     = σ e^α t
    
      e_t^i     :=[Z_t^i ]  = σ e^α t,
    
      f_t^i     := [(Z_t^i)^2 ]  = σ^2 e^2α t.

Introducing 

    ϕ(t,x)     := ∂_t u + ∑_i a x_i ∂_x_i u + 
     ∑_i(σ x_i)^2/2 ∂_x_i^2^2 u  
        = e^α t( αlog(∏_i x^i)  + ∑_i (a - σ^2/2)),

u(t,X_t) solves the PDE

    ∂_t u + ∑_i a x^i ∂_x_i u + 
    ∑_iσ^2/2 ∂_x_i^2^2 u - ϕ(t,x)= 0.

This semilinear PDE is related to the BSDE associated with the driver f(t,x) = -ϕ(t,x) for forward dynamics (<ref>).

Using some chosen ^d valued functions ψ and ^k valued functions κ,
 we express all dynamics in a McKean-Vlasov setting:

    X_t^i      = (a X^i_t +  ψ(Y_t, Z_t^i, [X_t^i],[(X_t^i)^2], [Y_t], [Y_t^2],
    [Z_t^i], [(Z^i_t)^2])
        - ψ(e^α tlog(∏_i X^i_t), σ e^α t, h_t^i,k_t^i,c_t,d_t,e_t^i,f_t^i)  t + σ X^i_t  W^i_t  
    
      X^i_0     = ξ
     Y_t     = - f(t,X_t,Y_t,Z_t,[X_t],[X_t^2], [Y_t], [Y_t^2],
    [Z_t], [Z^2_t])  t 
        + Z_t  W_t
    
    Y_T     = e^α Tlog(∏_i X^i_T)

with

    f(t,X_t,Y_t,Z_t,x_1,x_2,y_1,y_2,z_1,z_2 ) 
        = -ϕ(t,x) + κ(Y_t, Z_t,x_1,x_2,y_1,y_2,z_1,z_2 ) 
       -
    κ(e^α tlog(∏_i X^i_t), σ e^α t, h_t^i,k_t^i,c_t,d_t,e_t^i,f_t^i).

and f: ×^d××^1× d×^d×^d×××^1× d×^1× d↦.


We consider two models of this kind for numerical tests.



  §.§.§ A fully coupled linear example 

We consider a linear McKean-Vlasov FBSDE  in Y_t, Z_t and their law dynamics for X_t and Y_t:

    X_t^i      = (a X^i_t +  b(Y_t + Z_t^i + [X_t^i] + [Y_t] +
    [Z_t^i]) 
         - b(e^α tlog(∏_i=1^d X^i_t)+ σ e^α t + h_t^i + c_t + e_t^i)  t  + σ X^i_t  W^i_t  
    
      X^i_0     = ξ
      Y_t     =  (ϕ(t,X_t) + b(Y_t + 1/d∑_i=1^d Z_t^i + 1/d∑_i=1^d[X_t^i] + [Y_t] 
       + 1/d∑_i=1^d [Z_t^i])   - b(e^α tlog(∏_i=1^d X^i_t) + 1/d∑_i=1^d σ e^α t + 1/d∑_i=1^d h_t^i . 
       + . c_t +  1/d∑_i=1^d e_t^i) )  t + Z_t  W_t
    
    Y_T     = e^α Tlog(∏_i=1^d X^i_T).















  

We take a=b=0.1, α = 0.5, σ = 0.4, ξ = 1. The three algorithms demonstrate good performances on this test case. Both processes Y, Z are well represented by the neural network. However the Local Dynamic method is less precise than the global methods when the maturity grows. We see in Table <ref>, Figure <ref>, Figure <ref> that the Local Dynamic method is biased when T=1 when the Global methods achieve a great accuracy. It looks like the results of the Local Dynamic method cannot be improved since the loss flattens in Figure <ref>.



  §.§.§ A fully coupled quadratic example 

We consider a quadratic McKean-Vlasov FBSDE in Y_t, Z_t and their law dynamics for X_t and Y_t:
We consider a quadratic McKean-Vlasov FBSDE in Y_t, Z_t and their law dynamics for X_t and Y_t:

    X_t^i      = (a X^i_t +  b(Y_t + Z_t^i + [X_t^i] + [Y_t] +
    [Z_t^i]) + σ X^i_t  W^i_t  
         - b(e^α tlog(∏_i=1^d X^i_t)+ σ e^α t + h_t^i + c_t + e_t^i) 
         +  c (Y_t^2 + (Z_t^i)^2 + [(X_t^i)^2] + [Y_t^2] +
    [(Z_t^i)^2]) 
         - c(e^2α tlog(∏_i=1^d X^i_t)^2+ σ^2 e^2α t + (h_t^i)^2 + c_t^2 + (e_t^i)^2) )  t   
    
      X^i_0     = ξ
     Y_t     = (ϕ(t,X_t) + b(Y_t + 1/d∑_i=1^d Z_t^i + 1/d∑_i=1^d[X_t^i] + [Y_t]) 
        + b/d∑_i=1^d [Z_t^i])   - b(e^α tlog(∏_i=1^d X^i_t) + 1/d∑_i=1^d σ e^α t)  + c [Y_t^2]
        -b( 1/d∑_i=1^d h_t^i + c_t +  1/d∑_i=1^d e_t^i) + c(Y_t^2 + 1/d∑_i=1^d (Z_t^i)^2 + 1/d∑_i=1^d[(X_t^i)^2] ) 
       + c/d∑_i=1^d [(Z_t^i)^2])  - c(e^2α tlog(∏_i=1^d X^i_t)^2 + 1/d∑_i=1^d σ^2 e^2α t) 
       - c( 1/d∑_i=1^d (h_t^i)^2 + c_t^2+ 1/d∑_i=1^d (e_t^i)^2))  t + Z_t  W_t
    
    Y_T     = e^α Tlog(∏_i=1^d X^i_T).
  




















  

We take a=b=c=0.1, α = 0.5, σ = 0.4, ξ = 1. We observe in Table <ref> convergence of the methods for small maturities and divergence beyond T=1. Note that the dynamic estimation of the expectation prevents the algorithm to explode for T=1, contrarily to the Global Direct method. However, it does not converge to the true solution in this case. Indeed the loss plateaus at the value 2 in Figure <ref> (right), so the terminal condition of the BSDE is not properly respected. The Global Dynamic method also produces better result than the Global Direct one (see Figure <ref> and Table <ref>). We notice from Figure <ref> and Figure <ref> that the estimated Y,Z processes have the good shape but some errors are still present after convergence.

Concerning the Local Dynamic method, we see in Figure <ref> that the estimated expectations are stable around zero for a few iterations but then become negative. It may be due to the lack of a contraction for the fixed point problem. The loss explodes for T=1.5, as seen on the learning curve from Figure <ref>. For T=1, we see that it stays above 10.



§ CONCLUSION

We have shown that neural network methods can solve some moderate dimensional FBSDE of McKean-Vlasov type.
Comparing the different algorithms we find out that:

    
  * The dynamic update of the expectation is efficient in terms of computation speed (about 30% faster than direct method) and seems to smooth the learning curve.
    
  * For the mean-field games of controls example, the Pontryagin approach performs better than the Weak one for large maturities. On the contrary, the Weak approach is the best for small maturities.
    
  * For the fully coupled linear model we observe no convergence problem whereas for the fully coupled quadratic one we can solve only the equation on a small time horizon. However the Local Dynamic method is not very accurate for larger maturities.
    
  * The Local Dynamic method faces more difficulties for quadratic problems than the global methods do. It also requires more iterations, hence more time, to converge.
    
  * The expectation law estimation method does not work well and requires to empirically choose a proper penalization parameter, which is troublesome.
    
  * The methods can be used in dimension 10, thus applied to more realistic problems than usually. For instance, in the price impact model, the number of dimensions corresponds to the number of assets involved in the trading. Thus, developing methods able to deal with problems in high dimensions can help us to handle large portfolios.

We recommend the use of the global method combined with dynamic moment estimation which offers the best accuracy and training speed among all the tested methods. For linear quadratic mean-field games of controls, it appears to be better to use the Weak approach for small maturities and the Pontryagin method for larger time horizons. The use of a Local Dynamic method is possible but requires too many iterations to converge hence it is not competitive in terms of computation time.



