## Neural networks-based backward scheme for fully nonlinear PDEs ∗

Huyˆ en Pham †

Xavier Warin ‡

Maximilien Germain §

This version:

December 7, 2020

## Abstract

We propose a numerical method for solving high dimensional fully nonlinear partial differential equations (PDEs). Our algorithm estimates simultaneously by backward time induction the solution and its gradient by multi-layer neural networks, while the Hessian is approximated by automatic differentiation of the gradient at previous step. This methodology extends to the fully nonlinear case the approach recently proposed in [HPW20] for semi-linear PDEs. Numerical tests illustrate the performance and accuracy of our method on several examples in high dimension with non-linearity on the Hessian term including a linear quadratic control problem with control on the diffusion coefficient, Monge-Amp` ere equation and Hamilton-JacobiBellman equation in portfolio optimization.

Key words: Neural networks, fully nonlinear PDEs in high dimension, backward scheme.

MSC Classification: 60H35, 65C20, 65M12.

## 1 Introduction

This paper is devoted to the resolution in high dimension of fully nonlinear parabolic partial differential equations (PDEs) of the form

$$\dots \begin{cases} \partial _ { t } u + f (.,., u, D _ { x } u, D _ { x } ^ { 2 } u ) = 0, & \text{ on } [ 0, T ) \times \mathbb { R } ^ { d }, \\ u ( T,. ) = g, & \text{ on } \mathbb { R } ^ { d }, \end{cases}$$

with a non-linearity in the solution, its gradient D u x and its hessian D u 2 x via the function f ( t, x, y, z, γ ) defined on [0 , T ] × R d × R × R d × S d (where S d is the set of symmetric d × d matrices), and a terminal condition g .

The numerical resolution of this class of PDEs is far more difficult than the one of classical semi-linear PDEs where the nonlinear function f does not depend on γ . In fact, rather few methods are available to solve fully nonlinear equations even in moderate dimension.

- · First based on the work of [Che+07], an effective scheme developed in [FTW11] using some regression techniques has been shown to be convergent under some ellipticity conditions later removed by [Tan13]. Due to the use of basis functions, this scheme does not permit to solve PDE in dimension greater than 5.
- · A scheme based on nesting Monte Carlo has been recently proposed in [War18]. It seems to be effective in very high dimension for maturities T not too long and linearities not too important.
- · A numerical algorithm to solve fully nonlinear equations has been proposed by [BEJ19] based on the second order backward stochastic differential equations (2BSDE) representation of [Che+07] and global deep neural networks minimizing a terminal objective function, but no test on real fully nonlinear case is given. This extends the idea introduced in the pioneering papers [EHJ17; HJE18], which were the first serious works for using machine learning methods to solve high dimensional PDEs.
- · The Deep Galerkin method proposed in [SS18] based on some machine learning techniques and using some automatic differentiation of the solution seems to be effective on some cases. It has been tested in [AA+18] for example on the Merton problem.

∗ This work is supported by FiME, Laboratoire de Finance des March´s de l'Energie, and the 'Finance and Sustainable Devele opment' EDF - CACIB Chair.

‡ EDF R&amp;D &amp; FiME xavier.warin at edf.fr

† LPSM, Universit´ e de Paris, CREST-ENSAE &amp; FiME pham at lpsm.paris

§ EDF R&amp;D, LPSM, Universit´ e de Paris mgermain at lpsm.paris

In this article, we introduce a numerical method based on machine learning techniques and backward in time iterations, which extends the proposed schemes in [VSS18] for linear problems, and in the recent work [HPW20] for semi-linear PDEs. The approach in these works consists in estimating simultaneously the solution and its gradient by multi-layer neural networks by minimizing a sequence of loss functions defined in backward induction. A basic idea to extend this method to the fully nonlinear case would rely on the representation proposed in [Che+07]: at each time step t n of an Euler scheme, the Hessian D u 2 x at t n is approximated by a neural network minimizing some local L 2 criterion associated to a BSDE involving D u x at date t n +1 and D u 2 x . Then, the pair ( u, D x u ) at date t n is approximated/learned with a second minimization similarly as in the method described by [HPW20]. The first minimization can be implemented with different variations but numerical results show that the global scheme does not scale well with the dimension. Instability on the D u 2 x calculation rapidly propagates during the backward resolution. Besides, the methodology appears to be costly when using two optimizations at each time step. An alternative approach that we develop here, is to combine the ideas of [HPW20] and the splitting method in [Bec+19] in order to derive a new deep learning scheme that requires only one local optimization during the backward resolution for learning the pair ( u, D x u ) and approximating D u 2 x by automatic differentiation of the gradient computed at the previous step.

The outline of the paper is organized as follows. In Section 2, we briefly recall the mathematical description of the classical feedforward approximation, and then derive the proposed neural networks-based backward scheme. We test our method in Section 3 on various examples. First we illustrate our results with a PDE involving a nonlinearity of type uD u 2 x . Then, we consider a stochastic linear quadratic problem with controlled volatility where an analytic solution is available, and we test the performance and accuracy of our algorithm up to dimension 20. Next, we apply our algorithm to a Monge-Amp` ere equation, and finally, we provide numerical tests for the solution to fully nonlinear Hamilton-Jacobi-Bellman equation, with non-linearities of the form | D u x | 2 /D u 2 x , arising in portfolio selection problem with stochastic volatilities.

## 2 The proposed deep backward scheme

↦

Our aim is to numerically approximate the function u : [0 , T ] × R d → R , assumed to be the unique smooth solution to the fully nonlinear PDE (1.1) under suitable conditions. This will be achieved by means of neural networks approximations for u and its gradient D u x , relying on a backward scheme and training simulated data of some forward diffusion process. Approximations of PDE in high dimension by neural networks have now become quite popular, and are supported theoretically by recent results in [Hut+18] and [DLM20] showing their efficiency to overcome the curse of dimensionality.

## 2.1 Feedforward neural network to approximate functions

We denote by d 0 the dimension of the input variables, and d 1 the dimension of the output variable. A (deep) neural network is characterized by a number of layers L +1 ∈ N \ { 1 2 , } with m glyph[lscript] , glyph[lscript] = 0 , . . . , L , the number of neurons (units or nodes) on each layer: the first layer is the input layer with m 0 = d , the last layer is the output layer with m L = d 1 , and the L -1 layers between are called hidden layers, where we choose for simplicity the same dimension m glyph[lscript] = m glyph[lscript] , = 1 , . . . , L -1.

A feedforward neural network is a function from R d 0 to R d 1 defined as the composition

↦

$$x \in \mathbb { R } ^ { d } \longmapsto \, A _ { L } \circ \varrho \circ A _ { L - 1 } \circ \dots \circ \varrho \circ A _ { 1 } ( x ) \, \in \, \mathbb { R } \,.$$

Here A glyph[lscript] , glyph[lscript] = 1 , . . . , L are affine transformations: A 1 maps from R d 0 to R m , A , . . . , A 2 L -1 map from R m to R m , and A L maps from R m to R d 1 , represented by

$$A _ { \ell } ( x ) = \, \mathcal { W } _ { \ell } x + \beta _ { \ell },$$

for a matrix W glyph[lscript] called weight, and a vector β glyph[lscript] called bias term, glyph[rho1] : R → R is a nonlinear function, called activation function, and applied component-wise on the outputs of A glyph[lscript] , i.e., glyph[rho1] ( x , . . . , x 1 m ) = ( glyph[rho1] ( x 1 ) , . . . , glyph[rho1] ( x m )). Standard examples of activation functions are the sigmoid, the ReLu, the Elu, tanh.

All these matrices W glyph[lscript] and vectors β glyph[lscript] , glyph[lscript] = 1 , . . . , L , are the parameters of the neural network, and can be identified with an element θ ∈ R N m , where N m = ∑ L -1 glyph[lscript] =0 m glyph[lscript] (1+ m glyph[lscript] +1 ) = d 0 (1+ m )+ m (1+ m L )( -2)+ m (1+ d 1 ) is the number of parameters. We denote by N d ,d 0 1 ,L,m the set of all functions generated by (2.1) for θ ∈ R N m .

## 2.2 Forward-backward representation

Let us introduce a forward diffusion process

$$X _ { t } = X _ { 0 } + \int _ { 0 } ^ { t } \mu ( s, X _ { s } ) d s + \int _ { 0 } ^ { t } \sigma ( s, X _ { s } ) d W _ { s }, \ \ 0 \leq t \leq T,$$

where µ is a function defined on [0 , T ] × R d with values in R d , σ is a function defined on [0 , T ] × R d with values in M d the set of d × d matrices, and W a d -dimensional Brownian motion on some probability space (Ω , F , P ) equipped with a filtration F = ( F t ) 0 ≤ ≤ t T satisfying the usual conditions. The process X will be used for the simulation of training data in our deep learning algorithm, and we shall discuss later the choice of the drift and diffusion coefficients µ and σ , see Remark 2.3.

Let us next denote by ( Y, Z, Γ) the triple of F -adapted processes valued in R × R d × S d , defined by

$$Y _ { t } \, = \, u ( t, X _ { t } ), \ \ Z _ { t } \, = \, D _ { x } u ( t, X _ { t } ), \ \ \Gamma _ { t } \, = \, D _ { x } ^ { 2 } u ( t, X _ { t } ), \ \ 0 \leq t \leq T.$$

By Itˆ's formula applied to o u t, X ( t ), and since u is solution to (1.1), we see that ( Y, Z, Γ) satisfies the backward equation:

$$Y _ { t } = \ g ( X _ { T } ) - \int _ { t } ^ { T } \left [ \mu ( s, X _ { s } ). Z _ { s } + \frac { 1 } { 2 } \text{tr} ( \sigma \sigma ^ { \dagger } ( s, X _ { s } ) \Gamma _ { s } ) - f ( s, X _ { s }, Y _ { s }, Z _ { s }, \Gamma _ { s } ) \right ] d s \\ - \int _ { t } ^ { T } \sigma ^ { \dagger } ( s, X _ { s } ) Z _ { s }, d W _ { s }, \quad 0 \leq t \leq T.$$

Remark 2.1 This BSDE does not uniquely characterize a triple ( Y, Z, Γ) contrarily to the semilinear case (without a non-linearity with respect to Γ) in which proper assumptions on the equation coefficients provide existence and uniqueness for a solution couple ( Y, Z ). In the present case at least two options can be used to estimate the Γ component:

- · Rely on the 2BSDE representation from [Che+07] which extends the probabilistic representation of [PP90] for semilinear equations to the fully nonlinear case. It is the approach used by [BEJ19] with a global large minimization problem, as in [HJE18].
- · Compute the second order derivative by automatic differentiation. This is the point of view we adopt in this paper together with a local approach solving several small optimization problems. In this way, we provide an extension of [HPW20] to cover a broader range of nonlinear PDEs.

glyph[square]

## 2.3 Algorithm

We now provide a numerical approximation of the forward backward system (2.2)-(2.4), and consequently of the solution u (as well as its gradient D u x ) to the PDE (1.1).

We start from a time grid π = { t , i i = 0 , . . . , N } of [0 , T ], with t 0 = 0 &lt; t 1 &lt; . . . &lt; t N = T , and time steps ∆ t i := t i +1 -t i , i = 0 , . . . , N -1. The time discretization of the forward process X on π is then equal (typically when µ and σ are constants) or approximated by an Euler scheme:

$$X _ { t _ { i + 1 } } = \, X _ { t _ { i } } + \mu ( t _ { i }, X _ { t _ { i } } ) \Delta t _ { i } + \sigma ( t _ { i }, X _ { t _ { i } } ) \Delta W _ { t _ { i } }, \quad i = 0, \dots, N - 1,$$

where we set ∆ W t i := W t i +1 -W t i (by misuse of notation, we keep the same notation X for the continuous time diffusion process and its Euler scheme). The backward SDE (2.4) is approximated by the time discretized scheme

$$Y _ { t _ { i } } \simeq \, Y _ { t _ { i + 1 } } - [ \mu ( t _ { i }, X _ { t _ { i } } ). Z _ { t _ { i } } + \frac { 1 } { 2 } \text{tr} ( \sigma \sigma ^ { r } ( t _ { i }, X _ { t _ { i } } \Gamma _ { t _ { i } } ) - f ( t _ { i }, X _ { t _ { i } }, Y _ { t _ { i } }, Z _ { t _ { i } }, \Gamma _ { t _ { i } } ) ] \Delta t _ { i } - \sigma ^ { r } ( t _ { i }, X _ { t _ { i } } ) Z _ { t _ { i } }. \Delta W _ { t _ { i } },$$

that is written in forward form as

$$Y _ { t _ { i + 1 } } \simeq \ F ( t _ { i }, X _ { t _ { i } }, Y _ { t _ { i } }, Z _ { t _ { i } }, \Gamma _ { t _ { i } }, \Delta t _ { i }, \Delta W _ { t _ { i } } ), \ \ i = 0, \dots, N - 1,$$

with

$$F ( t, x, y, z, \gamma, h, \Delta ) & \coloneqq \ y - \tilde { f } ( t, x, y, z, \gamma ) h \ + \ z ^ { \dagger } \sigma ( t, x ) \Delta, \\ \tilde { f } ( t, x, y, z, \gamma ) & \coloneqq \ f ( t, x, y, z, \gamma ) - \mu ( t, x ). z - \frac { 1 } { 2 } \text{tr} ( \sigma \sigma ^ { \dagger } ( t, x ) \gamma ).$$

↦

The idea of the proposed scheme is the following. Similarly as in [HPW20], we approximate at each time t i , u t , . ( i ) and its gradient D u t , . x ( i ), by neural networks x ∈ R d → U ( i ( x θ , ; ) Z i ( x θ ; )) with parameter θ that are learned optimally by backward induction: suppose that ˆ U i +1 := U i +1 ( ; . θ ∗ i +1 ), ˆ Z i +1 := Z i +1 ( ; . θ ∗ i +1 ) is an

approximation of u t ( i +1 , . ) and D u t x ( i +1 , . ) at time t i +1 , then θ ∗ i is computed from the minimization of the quadratic loss function:

$$\hat { L } _ { i } ( \theta ) = \, \mathbb { E } \left | \hat { U } _ { i + 1 } - F ( t _ { i }, X _ { t _ { i } }, \mathcal { U } _ { i } ( X _ { t _ { i } } ; \theta ), \mathcal { Z } _ { i } ( X _ { t _ { i } } ; \theta ), D \hat { \mathcal { Z } } _ { i + 1 } ( \mathcal { T } ( X _ { t _ { i + 1 } } ) ), \Delta t _ { i }, \Delta W _ { t _ { i } } ) \right | ^ { 2 }$$

where T is a truncation operator such that T ( X ) is bounded for example by a quantile of the diffusion process and D ˆ Z i +1 stands for the automatic differentiation of ˆ Z i +1 . The idea behind the truncation is the following. During one step resolution, the estimation of the gradient is less accurate at the edge of the explored domain where samples are rarely generated. Differentiating the gradient gives a very oscillating Hessian at the edge of the domain. At the following time step resolution, these oscillations propagate to the gradient and the solution even if the domain where the oscillations occur is rarely attained. In order to avoid these oscillations, a truncation is achieved, permits to avoid that the oscillations of the neural network fit in zone where the simulations propagate scarcely to areas of importance. This truncation may be necessary to get convergence on some rather difficult cases. Of course this truncation is only valid if the real Hessian does not varies too much.

The intuition for the relevance of this scheme to the approximation of the PDE (1.1) is the following. From (2.3) and (2.5), the solution u to (1.1) should approximately satisfy

$$u ( t _ { i + 1 }, X _ { t _ { i + 1 } } ) \simeq \ F ( t _ { i }, X _ { t _ { i } }, u ( t _ { i }, X _ { t _ { i } } ), D _ { x } u ( t _ { i }, X _ { t _ { i } } ), D _ { x } ^ { 2 } u ( t _ { i }, X _ { t _ { i } } ), \Delta t _ { i }, \Delta W _ { t _ { i } } ).$$

Suppose that at time t i +1 , ˆ U i +1 is an estimation of u t ( i +1 , . ). Recalling the expression of F in (2.6), the quadratic loss function at time t i is then approximately equal to

$$\hat { L } _ { i } ( \theta ) & \simeq \mathbb { E } \Big | u ( t _ { i }, X _ { t _ { i } } ) - \mathcal { U } _ { i } ( X _ { t _ { i } } ; \theta ) + ( D _ { z } u ( t _ { i }, X _ { t _ { i } } ) - Z _ { i } ( X _ { t _ { i } } ; \theta ) \Big ) ^ { \tau } \sigma ( t _ { i }, X _ { t _ { i } } ) \Delta W _ { t _ { i } } \\ & \quad - \ \Delta t _ { i } \Big [ \tilde { f } ( t _ { i }, X _ { t _ { i } }, u ( t _ { i }, X _ { t _ { i } } ), D _ { z } u ( t _ { i }, X _ { t _ { i } } ), D _ { z } ^ { 2 } u ( t _ { i }, X _ { t _ { i } } ) ) - \tilde { f } ( t _ { i }, X _ { t _ { i } }, \mathcal { U } ( X _ { t _ { i } } ; \theta ), \mathcal { Z } _ { i } ( X _ { t _ { i } } ; \theta ), D \hat { \mathcal { Z } } _ { i + 1 } ( \mathcal { T } ( X _ { t _ { i + 1 } } ) ) ) \Big ] \Big | ^ { 2 }.$$

By assuming that ˜ has small non-linearities in its arguments ( f y, z, γ ), say Lipschitz, possibly with a suitable choice of µ , σ , the loss function is thus approximately equal to

$$\hat { L } _ { i } ( \theta ) \simeq \ ( 1 + O ( \Delta t _ { i } ) ) \mathbb { E } | u ( t _ { i }, X _ { t _ { i } } ) - \mathcal { U } _ { i } ( X _ { t _ { i } } ; \theta ) | ^ { 2 } + O ( \Delta t _ { i } ) \mathbb { E } | D _ { x } u ( t _ { i }, X _ { t _ { i } } ) - \mathcal { Z } _ { i } ( X _ { t _ { i } } ; \theta ) | ^ { 2 } + O ( | \Delta t _ { i } | ^ { 2 } ).$$

Therefore, by minimizing over θ this quadratic loss function, via stochastic gradient descent (SGD) based on simulations of ( X ,X t i t i +1 , ∆ W t i ) (called training data in the machine learning language), one expects the neural networks U i and Z i to learn/approximate better and better the functions u t , . ( i ) and D u t , x ( i ) in view of the universal approximation theorem for neural networks. The rigorous convergence of this algorithm is postponed to a future work.

To sum up, the global algorithm is given in Algo 1 in the case where g is Lipschitz and the derivative can be analytically calculated almost everywhere. If the derivative of g is not available, it can be calculated by automatic differentiation of the neural network approximation of g .

## Algorithm 1 Algorithm for fully non-linear equations.

- 1: Use a single deep neural network ( U N ( ; . θ , ) Z N ( ; . θ )) ∈ N d, 1+ d,L,m and minimize (by SGD)

$$\begin{cases} \hat { L } _ { N } ( \theta ) \, \coloneqq \, \mathbb { E } \left | \mathcal { U } _ { N } ( X _ { t _ { N } } ; \theta ) - g ( X _ { t _ { N } } ) \right | ^ { 2 } + \frac { \Delta t _ { N - 1 } } { d } \mathbb { E } \left | Z _ { N } ( X _ { t _ { N } } ; \theta ) - D g ( X _ { t _ { N } } ) \right | ^ { 2 } \\ \theta _ { N } ^ { * } \in \arg \min _ { \theta \in \mathbb { R } ^ { N _ { m } } } \hat { L } _ { N } ( \theta ). \end{cases}$$

- 2: ̂ U N = U N ( ; . θ ∗ N ), and set ̂ Z N = Z N ( ; . θ ∗ N )
- 3: for i = N -1 , . . . , 0 do
- 4: Use a single deep neural network ( U i ( ; . θ , ) Z i ( ; . θ )) ∈ N d, 1+ d,L,m for the approximation of ( u t , . ( i ) , D x u t , . ( i )), and compute (by SGD) the minimizer of the expected quadratic loss function

$$\begin{cases} \hat { L } _ { i } ( \theta ) \coloneqq \mathbb { E } \Big | \widehat { \mathcal { U } } _ { i + 1 } ( X _ { t _ { i + 1 } } ) \ - \ F ( t _ { i }, X _ { t _ { i } }, \mathcal { U } _ { i } ( X _ { t _ { i } } ; \theta ), \mathcal { Z } _ { i } ( X _ { t _ { i } } ; \theta ), D \hat { \mathcal { Z } } _ { i + 1 } ( \mathcal { T } ( X _ { t _ { i + 1 } } ) ), \Delta t _ { i }, \Delta W _ { t _ { i } } ) \Big | ^ { 2 } \\ \quad \theta _ { i } ^ { * } \in \arg \min _ { \theta \in \mathbb { R } ^ { N _ { m } } } \hat { L } _ { i } ( \theta ). \end{cases}$$

- 5: Update: ̂ U i = U i ( ; . θ ∗ i ), and set ̂ Z i = Z i ( ; . θ ∗ i ).

Remark 2.2 Several alternatives can be implemented for the computation of the second order derivative. A natural candidate would consist in choosing to approximate the solution u at time t i by a neural network U i

and estimate Γ i as the iterated automatic differentiation D 2 x U i . However, it is shown in [HPW20] that choosing only a single neural network for u and using its automatic derivative to estimate the Z component degrades the error in comparison to the choice of two neural networks U Z , . A similar behavior has been observed during our tests for this second order case and the most efficient choice was to compute the derivative of the Z network. This derivative can also be estimated at the current time step t i instead of t i +1 . However this method leads to an additional cost for the neural networks training by complicating the computation of the automatic gradients performed by Tensorflow during the backpropagation. It also leads numerically to worse results on the control estimation, as empirically observed in Table 5 and described in the related paragraph 'Comparison with an implicit version of the scheme'. For this reason, we decided to apply a splitting method and evaluate the Hessian at time t i +1 . For this reason, we decided to apply a splitting method and evaluate the Hessian at time t i +1 . glyph[square]

Remark 2.3 The diffusion process X is used for the training simulations in the stochastic gradient descent method for finding the minimizer of the quadratic loss function in (2.7), where the expectation is replaced by empirical average for numerical implementation. The choice of the drift and diffusion parameters are explained in Section 3.1. glyph[square]

## 3 Numerical results

We first construct an example with different non-linearities in the Hessian term and the solution. We graphically show that the solution is very well calculated in dimension d = 1 and then move to higher dimensions. We then use an example derived from a stochastic optimization problem with an analytic solution and show that we are able to accurately calculate the solution. Next, we consider the numerical resolution of the MongeAmp` ere equation, and finally, give some tests for a fully nonlinear Hamilton-Jacobi-Bellman equation arising from portfolio optimization with stochastic volatilities.

## 3.1 Choice of the algorithm hyperparameters

We describe in this paragraph how we choose the various hyperparameters of the algorithm and explain the learning strategy.

- · Parameters of the training simulations : the choice of the drift coefficient is typically related to the underlying probabilistic problem associated to the PDE (for example a stochastic control problem), and should drive the training process to regions of interest, e.g.., that are visited with large probability by the optimal state process in stochastic control. In practice, we can take a drift function µ . ( ) equal to the drift associated to some a priori control. This choice of control could be an optimal control for a related problem for which we know the solution, or could be the control obtained by the first iteration of the algorithm. The choice of the diffusion coefficient σ is also important: large σ induces a better exploration of the state space, but as we will see in most of examples below, it gives a scheme slowly converging to the solution with respect to the time discretization and it generates a higher variance on the results. Moreover, for the applications in stochastic control, we might explore some region that are visited with very small probabilities by the optimal state process, hence representing few interest. On the other hand, small σ means a weak exploration, and we might lack information and precision on some region of the state space: the solution calculated at each time step is far more sensitive to very local errors induced by the neural network approximation and tends to generate a bias. Therefore a trade off has to be found between rather high variance with slow convergence in time and fast convergence in time with a potential bias. We also refer to [NR20] for a discussion on the role of the diffusion coefficient.

In practice and for the numerical examples in the next section, we test the scheme for different σ and by varying the number of time steps, and if it converges to the same solution, one can consider that we have obtained the correct solution. We also show the impact of the choice of the diffusion coefficient σ .

- · Parameters of truncation : Given the training simulations X , we choose a truncation operator T p indexed by a parameter p close to 1, so that T p ( X t ) corresponds to a truncation of X t at a given quantile φ p . In the numerical tests, we shall vary p between 0 95 and 0 999. . .
- · Parameters of the optimization algorithm over neural networks : In the whole numerical part, we use a classical Feedforward network using layers with m neurons each and a tanh activation function, the output layer uses an identity activation function. At each time step the resolution of equation (2.7) is achieved using a mini-batch with 1000 training trajectories. The training and learning rate adaptation procedure is the following:
- · Every 40 inner gradient descent iterations, the loss is checked on 10000 validation trajectories.
- · This optimization sequence is repeated with 200 outer iterations for the first optimization step at date t N = T and only 100 outer iterations at the dates t i with i &lt; N .

- · An average of the loss calculated on 10 successive outer iterations is performed. If the decrease of the average loss every 10 outer iterations is less than 5% then the learning rate is divided by 2.

The optimization is performed using the Adam gradient descent algorithm, see [KB14]. Notice that the adaptation of the learning rate is not common with the Adam method but in our case it appears to be crucial to have a steady diminution of the loss of the objective function. The procedure is also described in [CWNMW19] and the chosen parameters are similar to this article. At the initial optimization step at time t N = T , the learning rate is taken equal to 1 E -2 and at the following optimization steps, we start with a learning rate equal to 1 E -3.

During time resolution, it is far more effective to initialize the solution of equations (2.7) with the solution ( U Z , ) at the next time step. Indeed the previously computed values at time step t i +1 are good approximations of the processes at time step t i if the PDE solution and its gradient are continuous. All experiments are achieved using Tensorflow [Aba+15]. In the sequel, the PDE solutions on curves are calculated as the average of 10 runs. We provide the standard deviation associated to these results. We also show the influence of the number of neurons on the accuracy of the results.

## 3.2 A non-linearity in uD u 2 x

We consider a generator in the form

$$f ( t, x, y, z, \gamma ) = \ y t r ( \gamma ) + \frac { y } { 2 } + 2 y ^ { 2 } - 2 y ^ { 4 } e ^ { - ( T - t ) },$$

and g x ( ) = tanh ( ∑ d i =1 x i √ d ) , so that an analytical solution is available:

$$u ( t, x ) = \, \tanh \left ( \frac { \sum _ { i = 1 } ^ { d } x _ { i } } { \sqrt { d } } \right ) e ^ { - \frac { T - t } { 2 } }.$$

We fix the horizon T = 1, and choose to evaluate the solution at t = 0 and x = 0 5 . 1 I d √ d (here 1 I d denotes the vector in R d with all components equal to 1), for which u t, x ( ) = 0 761902 while its derivative is equal to . 1 2966. . This initial value x is chosen such that independently of the dimension the solution is varying around this point and not in a region where the tanh function is close to -1 or 1.

The coefficients of the forward process used to solve the equation are (here I d is the identity d × d -matrix)

$$\sigma = \frac { \hat { \sigma } } { \sqrt { d } } \mathbf I _ { d }, \quad \mu = 0,$$

and here the truncation operator is chosen equal to

$$\mathcal { T } _ { p } ( X _ { t } ^ { 0, x } ) = \min \left \{ \max [ x - \sigma \sqrt { t } \phi _ { p }, X _ { t } ^ { 0, x } ], x + \sigma \sqrt { t } \phi _ { p } \right \},$$

where φ p = N -1 ( p ), with N is the CDF of a unit centered Gaussian random variable.

In the numerical results, we take p = 0 999 and . m = 20 neurons. We first begin in dimension d = 1, and show in Figure 1 how u , D u x and D u 2 x are well approximated by the resolution method.

On Figure 2, we check the convergence, for different values of ˆ of both the solution σ u and its derivative at point x and date 0. Standard deviation of the function value is very low and the standard deviation of the derivative still being low.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-35910177-afda-4196-92c2-4fcdb146befd.jpg)

Y at date t = 0 5. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-3a7c25fd-bdd0-41bc-8dfc-1a5f52251f41.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-ea9a8ad8-c34a-43f2-a4f4-9f83e727f852.jpg)

Z at date t = 0 5 .

Γ at date t = 0 5 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-b1f990aa-2af7-46ee-a7d9-8010d0b23c41.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-4a5b6bce-9408-4db7-9c94-a94fdff0eaf4.jpg)

Y at date t = 0 006125. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-84ebd0aa-a78d-4a9d-9de6-87791db21d1f.jpg)

Z

at date

t

= 0 006125

.

Γ at date

t

= 0 006125

.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-1a3369cf-2c53-4ad9-b46e-e70949b3ebc0.jpg)

Figure 1: A single valuation run for test case one 1D using 160 time steps, ˆ = 2 , σ . p = 0 999, 20 neurons, 2 . layers.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-7ae51068-2a4a-4e54-9e47-7f15ebb6e186.jpg)

Convergence of u depending on ˆ σ

Standard deviation of u

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-12ae882f-251a-4294-8968-d0718f6d3dd4.jpg)

Convergence of D u x depending on ˆ σ

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-9427a79f-50d1-49e0-92a5-f4764c9a3f56.jpg)

Standard deviation of D u x

Figure 2: Convergence in 1D of the case one, number of neurons par layer equal to 20, 2 layers, p = 0 999. .

As the dimension increases, we have to increase the value of ˆ of the forward process. σ In dimension 3, the value ˆ = 0 5 gives high standard deviation in the result obtained as shown on Figure 3, while in dimension 10, σ . see Figure 4, we see that the value ˆ = 1 is too low to give good results. σ We also clearly notice that in 10D, a smaller time step should be used but in our test cases we decided to consider a maximum number of time steps equal to 160.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-eec995d8-0c6c-4c2b-ba9a-bf800009c54d.jpg)

Convergence of u depending on ˆ σ

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-029d1001-9fbb-4d77-b8a5-e4379f739ce5.jpg)

Convergence of D u x (first component) depending on ˆ σ

Convergence of u depending on ˆ σ

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-4b7219f3-fee1-4024-a429-ea7a3a9ccade.jpg)

Figure 3: Convergence in 3D of the case one, number of neurons par layer equal to 20, 2 layers, p = 0 999. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-40fff7e1-0184-440e-a3cb-720b5ffef569.jpg)

Convergence of D u x depending on ˆ (first component) σ

Figure 4: Convergence in 10D of the case one, number of neurons par layer equal to 20, 2 layers, p = 0 999. .

On this simple test case, the dimension is not a problem and very good results are obtained in dimension 20 or above with only 20 neurons and 2 layers.

## 3.3 A linear quadratic stochastic test case.

In this example, we consider a controlled process X = X α with dynamics in R d according to

$$d \mathcal { X } _ { t } = \, ( A \mathcal { X } _ { t } + B \alpha _ { t } ) d t + D \alpha _ { t } d W _ { t }, \ \ 0 \leq t \leq T, \ \ \mathcal { X } _ { 0 } \ = \ x,$$

where W is a real Brownian motion, the control process α is valued in R , and the constant coefficients A ∈ M d , B ∈ R d , D ∈ R d . The quadratic cost functional to be minimized is

$$J ( \alpha ) = & \mathbb { E } \left [ \int _ { 0 } ^ { T } \left ( \mathcal { X } _ { t } ^ { \tau } Q \mathcal { X } _ { t } + \alpha _ { t } ^ { 2 } N \right ) d t + \mathcal { X } _ { T } ^ { \tau } P \mathcal { X } _ { T } \right ],$$

where P , Q are non negative d × d symmetric matrices and N ∈ R is strictly positive. The Bellman equation associated to this stochastic control problem is:

$$\frac { \partial u } { \partial t } + \inf _ { a \in \mathbb { R } } \left [ ( A x + B a ). D _ { x } u + \frac { a ^ { 2 } } { 2 } \text{tr} ( D D ^ { \tau } D _ { x } ^ { 2 } u ) + x ^ { \tau } Q x + N a ^ { 2 } \right ] = \, 0, \ \ ( t, x ) \in [ 0, T ) \times \mathbb { R } ^ { d },$$

$$u ( T, x ) = \ x ^ { \tau } P x, \ \ x \in \mathbb { R } ^ { d },$$

which can be rewritten as a fully nonlinear equation in the form (1.1) with

$$f ( t, x, y, z, \gamma ) = \ x ^ { \tau } Q x + A x. z - \frac { 1 } { 2 } \frac { | B ^ { \tau } z | ^ { 2 } } { \text{tr} ( D D ^ { \tau } \gamma ) + 2 N }.$$

An explicit solution to this PDE is given by

$$u ( t, x ) = \, x ^ { \tau } K ( t ) x,$$

where K t ( ) is non negative d × d symmetric matrix function solution to the Riccati equation:

$$\dot { K } + A ^ { \top } K + K A + Q - \frac { K B B ^ { \top } K } { N + D ^ { \top } K D } = 0, \ \ K ( T ) = P.$$

We take T = 1. The coefficients of the forward process used to solve the equation are

$$\sigma = \frac { \hat { \sigma } } { \sqrt { d } } \mathbf I _ { d }, \ \mu ( t, x ) = A x.$$

In our numerical example we take the following parameters for the optimization problem:

$$A = \mathbf I _ { d }, \ B = D = \mathbf I _ { d }, \ Q = P = \frac { 1 } { d } \mathbf I _ { d }, \ \ N = d$$

and we want to estimate the solution at x = 1 I d .

In this example, the truncation operator (indexed by p between 0 and 1 and close to 1) is as follows:

$$\mathcal { T } _ { p } ( X _ { t } ^ { x } ) = \, \min \left \{ \, \max \, [ x e ^ { \hat { A } t } - \sigma \sqrt { \frac { e ^ { 2 \hat { A } t } - \hat { 1 } } { 2 \hat { A } } } \phi _ { p }, X _ { t } ^ { x } ], x e ^ { \hat { A } t } + \sigma \sqrt { \frac { e ^ { 2 \hat { A } t } - \hat { 1 } } { 2 \hat { A } } } \phi _ { p } \right \},$$

where φ p = N -1 ( p ), ˆ is a vector so that A ˆ A i = A ii , i = 1 , ..., d , ˆ is a unit vector, and the square root is taken 1 componentwise.

On Figure 5 we give the solution of the PDE with d = 1 using ˆ = 1 5 obtained for two dates: σ . at t = 0 5 . and at t close to zero. We observe that we have a very good estimation of the function value and a correct one of the Γ value at date t = 0 5. . The precision remains good for Γ close to t = 0 and very good for u and D u x .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-a34e3a96-2d8a-459f-825f-8ad1c4bae8f6.jpg)

Y at date t = 0 5. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-aa9b1f63-03a9-4ceb-bb00-49409807bbe9.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-c1c3b0b2-54dc-44ad-8422-69d4574a429f.jpg)

Γ at date t = 0 5 .

Z at date t = 0 5 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-a64d0e3f-4b82-4007-b037-95b9dfd01835.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-618e9267-528c-4274-89c1-ac0563018ad6.jpg)

Y at date t = 0 006125. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-b169763e-6d8d-4d21-adb4-c18a25e5ec78.jpg)

Z

at date

t

= 0 006125

.

Γ at date t = 0 006125 .

Figure 5: Test case linear quadratic 1D using 160 time steps, ˆ = 1 5, σ . p = 0 999, 100 neurons. .

On Figure 6, we give the results obtained in dimension d = 1 by varying ˆ. σ For a value of ˆ = 2, the σ standard deviation of the result becomes far higher than with ˆ = 0 5 or 1 σ . .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-f0a1c6f3-ece8-4a73-b24b-79096e990200.jpg)

Convergence of u depending on ˆ. σ

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-914e745a-c311-4130-a1d4-d6e78b05f145.jpg)

Standard deviation of u

Figure 6: Convergence in 1D of the linear quadratic case, number of neurons par layer equal to 50, 2 layers, p = 0 999. .

On Figure 7, for d = 3, we take a quite low truncation factor p = 0 95 and observe that the number of . neurons to take has to be rather high. We have also checked that taking a number of hidden layers equal to 3 does not improve the results.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-d621f4cf-f81e-42c3-b876-a927f8d6e34d.jpg)

10 neurons

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-98f632e6-8c6f-4039-a1f4-4e228eaed750.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-0e0bed95-273d-495f-b077-a0cda166a7a9.jpg)

20 neurons

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-63b70642-a3a6-4ef6-ac47-befb2548469d.jpg)

30 neurons

50 neurons

Figure 7: Convergence in 3D of the linear quadratic case, 2 layers, testing the influence of the number of neurons, truncation p = 0 95. .

On Figure 8, for d = 3, we give the same graphs for a higher truncation factor. As we take a higher truncation factor, the results are improved by taking a higher number of neurons (100 in the figure below).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-ee5bf09d-4a46-4cfa-bdc7-51f056f40069.jpg)

10 neurons

20 neurons

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-737876b6-7d2f-46a3-8f90-a460585e75fd.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-6e88682f-c0e0-43ff-8b40-7a55074f6f8a.jpg)

50 neurons

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-3e9bb0bb-02ca-41b4-ac12-9c0aa42e41a6.jpg)

100 neurons

Figure 8: Convergence in 3D of the linear quadratic case, 2 layers, testing the influence of the number of neurons, truncation p = 0 99. .

On Figure 9, we observe in dimension 7 the influence of the number of neurons on the result for a high truncation factor p = 0 999. . We clearly have a bias for a number of neurons equal to 50. This bias disappears when the number of neurons increases to 100.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-4e8f64fa-8e25-413e-943f-e2232d17f79f.jpg)

Convergence with 50 neurons

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-6c188818-ecd2-448b-8ee7-9b22fcdf3876.jpg)

Convergence with 100 neurons

Figure 9: Convergence in 7D of the linear quadratic case, 2 layers, p = 0 999. .

On Figure 10, for d = 7, we check that influence of the truncation factor appears to be slow for higher dimensions.

Figure 10: Function value convergence in 7D of the linear quadratic case with 2 layers, 100 neurons, testing p , using ˆ = 2 σ

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-9be6a082-7fbc-4035-a2b7-8f255b67f112.jpg)

.

Finally, we give results in dimension 10, 15 and 20 for p = 0 999 on Figures 11, 12. . We observe that the number a neurons with 2 hidden layers has to increase with the dimension but also that the increase is rather slow in contrast with the case of one hidden layer as theoretically shown in [Pin99]. For ˆ = 5 we had to take σ 300 neurons to get very accurate results.

Figure 11: Function value convergence in 10D and 15D of the linear quadratic case with 2 layers, p = 0 999. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-0ae8c7f6-3cc3-4a3f-89a5-50d63b05051b.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-c50ee068-9cd2-44ea-85b5-57cbb8489080.jpg)

Figure 12: Function value convergence in 20D of the linear quadratic case with 2 layers, p = 0 999. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-a25fc6b6-a2e5-4f00-b52c-6a52528f954e.jpg)

## 3.4 Monge-Amp` ere equation

Let us consider the parabolic Monge-Amp` ere equation

$$\begin{cases} \partial _ { t } u + \det ( D _ { x } ^ { 2 } u ) \, = \, h ( x ), & ( t, x ) \in [ 0, T ] \times \mathbb { R } ^ { d }, \\ u ( T, x ) \, = \, g ( x ), \end{cases}$$

where det( D u 2 x ) is the determinant of the Hessian matrix D u 2 x . It is in the form (1.1) with

$$f ( t, x, \gamma ) = \, \det ( \gamma ) - h ( x ).$$

We test our algorithm by choosing a C 2 function g , then compute G = det( D g 2 x ), and set h := G -1. Then, by construction, the function

$$u ( t, x ) = \ g ( x ) + T - t,$$

is solution to the Monge-Amp` ere equation (3.1). We choose g x ( ) = cos( ∑ d i =1 x / i √ d ), and we shall train with the forward process X = x 0 + W , where W is a d -dimensional Brownian motion. On this example, we use neural networks with 3 hidden layers, d + 10 neurons per layer, and we do not need to apply any truncation to the forward process X . Actually, we observe that adding a truncation worsens the results. For choosing the truncation level, we first test the method with no truncation before decreasing the quantile parameter p . In the Monge-Amp` ere case the best results are obtained without any truncation. It may be caused by the oscillation of the Hessian.

The following table gives the results in dimension d = 5, 15, and for T = 1.

|   Dimension d |   Averaged value |   Standard deviation |   Relative error (%) |   Theoretical solution |
|---------------|------------------|----------------------|----------------------|------------------------|
|             5 |          0.37901 |              0.00312 |                 0.97 |               0.382727 |
|            15 |          0.25276 |              0.00235 |                 1.17 |               0.255754 |

Table 1: Estimate of u (0 , x 0 = 1 ) on the Monge Ampere problem (3.1) with 5 N = 120. Average and standard deviation observed over 10 independent runs are reported.

## 3.5 Portfolio selection

We consider a portfolio selection problem formulated as follows. There are n risky assets of uncorrelated price process P = ( P , . . . , P 1 n ) with dynamics

$$\mathrm d P _ { t } ^ { i } = \, P _ { t } ^ { i } \sigma ( V _ { t } ^ { i } ) \left [ \lambda _ { i } ( V _ { t } ^ { i } ) \mathrm d t + \mathrm d W _ { t } ^ { i } \right ], \ \ i = 1, \dots, n,$$

where W = ( W ,...,W 1 n ) is a n -dimensional Brownian motion, b = ( b 1 , . . . , b n ) is the rate of return of the assets, λ = ( λ , . . . , λ 1 n ) is the risk premium of the assets, σ is a positive function (e.g. σ v ( ) = e v corresponding to the Scott model), and V = ( V 1 , . . . , V n ) is the volatility factor modeled by an Ornstein-Uhlenbeck (O.U.) process

$$d V _ { t } ^ { i } = \kappa _ { i } [ \theta _ { i } - V _ { t } ^ { i } ] { \mathrm d } t + \nu _ { i } { \mathrm d } B _ { t } ^ { i }, \ \ i = 1, \dots, n,$$

with κ , θ , ν i i i &gt; 0, and B = ( B ,. . . , B 1 n ) a n -dimensional Brownian motion, s.t. d &lt; W ,B i j &gt; = δ ij ρ ij dt , with ρ i := ρ ii ∈ ( -1 1). , An agent can invest at any time an amount α t = ( α , . . . , α 1 t n t ) in the stocks, which generates a wealth process X = X α governed by

$$d \mathcal { X } _ { t } = \, \sum _ { i = 1 } ^ { n } \alpha _ { t } ^ { i } \sigma ( V _ { t } ^ { i } ) \left [ \lambda _ { i } ( V _ { t } ^ { i } ) \mathrm d t + \mathrm d W _ { t } ^ { i } \right ].$$

The objective of the agent is to maximize her expected utility from terminal wealth:

$$\mathbb { E } \left [ U ( \mathcal { X } _ { T } ^ { \alpha } ) \right ] \ \leftarrow \ \maximize \, \text{over} \, \alpha$$

It is well-known that the solution to this problem can be characterized by the dynamic programming method (see e.g. [Pha09]), which leads to the Hamilton-Jacobi-Bellman for the value function on [0 , T ) × R × R n :

$$\begin{cases} \partial _ { i } u + \sum _ { i = 1 } ^ { n } \left [ \kappa _ { i } ( \theta _ { i } - v _ { i } ) \partial _ { v _ { i } } u + \frac { 1 } { 2 } \nu _ { i } ^ { 2 } \partial _ { v _ { i } } ^ { 2 } u \right ] = \, \frac { 1 } { 2 } R ( v ) \frac { ( \partial _ { x } u ) ^ { 2 } } { \partial _ { x x } ^ { 2 } u } + \sum _ { i = 1 } ^ { n } \left [ \rho _ { i } \lambda _ { i } ( v _ { i } ) \nu _ { i } \frac { \partial _ { x } u \partial _ { x v _ { i } } ^ { 2 } u } { \partial _ { x x } ^ { 2 } u } + \frac { 1 } { 2 } \rho _ { i } ^ { 2 } \nu _ { i } ^ { 2 } \frac { ( \partial _ { x v _ { i } } ^ { 2 } u ) ^ { 2 } } { \partial _ { x x } ^ { 2 } u } \right ] \\ u ( T, x, v ) \, = \, U ( x ), \quad x \in \mathbb { R }, \, v \in \mathbb { R } ^ { n }, \end{cases}$$

with a Sharpe ratio R v ( ) := | λ v ( ) | 2 , for v = ( v , . . . , v 1 n ) ∈ (0 , ∞ ) n . The optimal portfolio strategy is then given in feedback form by α ∗ t = ˆ( a t, X ∗ t , V t ), where ˆ = (ˆ a a , . . . , a 1 ˆ n ) is given by

$$\hat { a } _ { i } ( t, x, v ) = \, - \frac { 1 } { \sigma ( v _ { i } ) } \left ( \lambda _ { i } ( v _ { i } ) \frac { \partial _ { x } u } { \partial _ { x x } ^ { 2 } u } + \rho _ { i } \nu _ { i } \frac { \partial _ { x v _ { i } } ^ { 2 } u } { \partial _ { x x } ^ { 2 } u } \right ), \ \ ( t, x, v = ( v _ { 1 }, \dots, v _ { n } ) ) \in [ 0, T ) \times \mathbb { R } \times \mathbb { R } ^ { n },$$

for i = 1 , . . . , n . This Bellman equation is in the form (1.1) with

$$f ( t, x, y, z, \gamma ) = \sum _ { i = 1 } ^ { n } \left [ \kappa _ { i } ( \theta _ { i } - v _ { i } ) z _ { i } + \frac { 1 } { 2 } \nu _ { i } ^ { 2 } \gamma _ { i i } \right ] - \frac { 1 } { 2 } R ( v ) \frac { z _ { 0 } ^ { 2 } } { \gamma _ { 0 0 } } - \sum _ { i = 1 } ^ { n } \left [ \rho _ { i } \lambda _ { i } ( v _ { i } ) \nu _ { i } \frac { z _ { 0 } \gamma _ { 0 i } } { \gamma _ { 0 0 } } + \frac { 1 } { 2 } \rho _ { i } ^ { 2 } \nu _ { i } ^ { 2 } \frac { ( \gamma _ { 0 i } ) ^ { 2 } } { \gamma _ { 0 0 } } \right ],$$

for x = (x , v ) ∈ R n +1 , z = ( z , . . . , z 0 n ) ∈ R n +1 , γ = ( γ ij ) 0 ≤ i,j ≤ n ∈ S n +1 , and displays a high non-linearity in the Hessian argument γ .

The truncation operator indexed by a parameter p is chosen equal to

$$\mathcal { T } _ { p } ( X _ { t } ^ { 0, x } ) = \min \left \{ \max [ x + \mu t - \sigma \sqrt { t } \phi _ { p }, X _ { t } ^ { 0, x } ], x + \mu t + \sigma \sqrt { t } \phi _ { p } \right \},$$

where φ p = N -1 ( p ), N is the CDF of a unit centered Gaussian random variable. We use neural networks with 2 hidden layers and d + 10 neurons per layer. We shall test this example when the utility function U is of exponential form: U x ( ) = -exp( -ηx ), with η &gt; 0, and under different cases for which closed-form solutions are available:

- (1) Merton problem. This corresponds to a degenerate case where the factor V , hence the volatility σ and the risk premium λ are constant, so that the generator of Bellman equation reduces to

$$f ( t, x, y, z, \gamma ) \, = \, - \frac { 1 } { 2 } | \lambda | ^ { 2 } \frac { z ^ { 2 } } { \gamma }, \quad ( t, x, y, z ) \in [ 0, T ] \times \mathbb { R } \times \mathbb { R },$$

with explicit solution given by:

$$u ( t, x ) = \ e ^ { - ( T - t ) \frac { | \lambda | ^ { 2 } } { 2 } } U ( x ), \ \hat { a } _ { i } \ = \ \frac { \lambda _ { i } } { \eta \sigma }.$$

We train with the forward process

$$X _ { k + 1 } = \, X _ { k } + \lambda \Delta t _ { k } + \Delta W _ { k }, \ \ k = 0, \dots, N, \, X _ { 0 } \, = \, x _ { 0 }.$$

- (2) One risky asset: n = 1. A quasi-explicit solution is provided in [Zar01]:

$$u ( t, x, v ) = \ U ( x ) w ( t, v ), \quad \text{with} \ w ( t, v ) \ = \ \left \| \exp \left ( - \, \frac { 1 } { 2 } \int _ { t } ^ { T } R ( \hat { V } _ { s } ^ { t, v } ) d s \right ) \right \| _ { _ { L ^ { 1 - \rho ^ { 2 } } } }$$

where ˆ V t,v s is the solution to the modified O.U. model

$$\mathrm d \hat { V } _ { s } = \ [ \kappa ( \theta - \hat { V } _ { s } ) - \rho \nu \lambda ( \hat { V } _ { s } ) ] \mathrm d s + \nu \mathrm d B _ { s }, \ \ s \geq t, \, \hat { V } _ { t } = v.$$

We test our algorithm with λ v ( ) = λv , λ &gt; 0, for which we have an explicit solution:

$$w ( t, v ) = \, \exp \left ( \, - \phi ( t ) \frac { v ^ { 2 } } { 2 } - \psi ( t ) v - \chi ( t ) \right ), \ \ ( t, v ) \in [ 0, T ] \times \mathbb { R },$$

where ( φ, ψ, χ ) are solutions of the Riccati system of ODEs:

$$\dot { \phi } - 2 \bar { \kappa } \phi - \nu ^ { 2 } ( 1 - \rho ^ { 2 } ) \phi ^ { 2 } + \lambda ^ { 2 } & = \ 0, \quad \phi ( T ) \ = \ 0, \\ \dot { \psi } - ( \bar { \kappa } + \nu ^ { 2 } ( 1 - \rho ^ { 2 } ) \phi ) \psi + \kappa \theta \phi & = \ 0, \quad \psi ( T ) \ = \ 0, \\ \dot { \chi } + \kappa \theta \psi - \frac { \nu ^ { 2 } } { 2 } ( - \phi + ( 1 - \rho ^ { 2 } ) \psi ^ { 2 } ) & = \ 0, \quad \chi ( T ) \ = \ 0,$$

$$2$$

with ¯ = κ κ + ρνλ , and explicitly given by (see e.g. Appendix in [SZ99])

$$\text{with $\bar{ }k=\kappa+\rho\nu\lambda$, and explicitly given by (see e.g. Appendix in [SZ99])} \\ \phi ( t ) & = \lambda ^ { 2 } \frac { \sinh ( \hat { \kappa } ( T - t ) ) } { \hat { \kappa } \cosh ( \hat { \kappa } ( T - t ) ) + \bar { \kappa } \sinh ( \hat { \kappa } ( T - t ) ) } \\ \psi ( t ) & = \lambda ^ { 2 } \frac { \kappa \theta } { \hat { \kappa } } \frac { \cosh ( \hat { \kappa } ( T - t ) ) - 1 } { \hat { \kappa } \cosh ( \hat { \kappa } ( T - t ) ) + \bar { \kappa } \sinh ( \hat { \kappa } ( T - t ) ) } \\ \chi ( t ) & = \frac { 1 } { 2 ( 1 - \rho ^ { 2 } ) } \ln \left [ \cosh ( \hat { \kappa } ( T - t ) ) + \frac { \bar { \kappa } } { \hat { \kappa } } \sinh ( \hat { \kappa } ( T - t ) ) \right ] - \frac { 1 } { 2 ( 1 - \rho ^ { 2 } ) } \bar { \kappa } ( T - t ) \\ & \quad - \lambda ^ { 2 } \frac { ( \kappa \theta ) ^ { 2 } } { \hat { \kappa } ^ { 2 } } \left [ \frac { \sinh ( \hat { \kappa } ( T - t ) ) } { \hat { \kappa } \cosh ( \hat { \kappa } ( T - t ) ) + \bar { \kappa } \sinh ( \hat { \kappa } ( T - t ) ) } - ( T - t ) \right ] \\ & \quad - \lambda ^ { 2 } \frac { ( \kappa \theta ) ^ { 2 } \bar { \kappa } } { \hat { \kappa } ^ { 3 } } \frac { \cosh ( \hat { \kappa } ( T - t ) ) - 1 } { \hat { \kappa } \cosh ( \hat { \kappa } ( T - t ) ) + \bar { \kappa } \sinh ( \hat { \kappa } ( T - t ) ) }, \\ \text{with $\bar{ }k=\sqrt{\kappa^2+2\rho\nu\lambda\kappa+\gamma^2\lambda^2}$. We train with the forward process}$$

with ˆ = κ √ κ 2 +2 ρνλκ + γ λ 2 2 . We train with the forward process

$$\mathcal { X } _ { k + 1 } & = \ \mathcal { X } _ { k } + \lambda \theta \Delta t _ { k } + \Delta W _ { k }, \quad k = 0, \dots, N - 1, \ \mathcal { X } _ { 0 } = x _ { 0 }, \\ V _ { k + 1 } & = \ V _ { k } + \nu \Delta B _ { k }, \quad k = 0, \dots, N - 1, \ V _ { 0 } \ = \ \theta.$$

- (3) No leverage effect, i.e., ρ i = 0 , i = 1 , . . . , n . In this case, there is a quasi-explicit solution given by

$$u ( t, x, v ) = \ U ( x ) w ( t, v ), \ \text{ with } w ( t, v ) \ = \ \mathbb { E } \left [ \exp \left ( - \, \frac { 1 } { 2 } \int _ { t } ^ { T } R ( V _ { s } ^ { t, v } ) d s \right ) \right ], \ ( t, v ) \in [ 0, T ] \times \mathbb { R } ^ { n }, \ ( 3. 4 )$$

where V t,v is the solution to (3.2), starting from v at time t . We test our algorithm with λ i ( v ) = λ v i i , λ i &gt; 0, i = 1 , . . . , n , v = ( v , . . . , v 1 n ), for which we have an explicit solution given by

$$\text{$\text{$r\colon$} } \text{$v\colon$} \text{$i=1$, $d$, $n$, $v=(v_{1},\dots,v_{n})$, for which we have an explicit solution given by} \\ w ( t, v ) = \, \exp \left ( - \sum _ { i = 1 } ^ { n } \left [ \phi _ { i } ( t ) \frac { v _ { 2 } ^ { 2 } } { 2 } + \psi _ { i } ( t ) v _ { i } + \chi _ { i } ( t ) \right ] \right ), \quad ( t, v ) \in [ 0, T ] \times \mathbb { R } ^ { n }, \\ \phi _ { i } ( t ) = \, \lambda _ { i } ^ { 2 } \frac { \sinh ( \tilde { x } _ { i } ( T - t ) ) } { \tilde { x } _ { i } \sinh ( \tilde { x } _ { i } ( T - t ) ) + \tilde { x } _ { i } \cosh ( \tilde { x } _ { i } ( T - t ) ) } \\ \psi _ { i } ( t ) = \, \lambda _ { i } ^ { 2 } \frac { \kappa _ { i } \theta _ { i } } { \tilde { x } _ { i } \sinh ( \tilde { x } _ { i } ( T - t ) ) + \tilde { x } _ { i } \cosh ( \tilde { x } _ { i } ( T - t ) ) } \\ \chi _ { i } ( t ) = \, \frac { 1 } { 2 } \ln \left [ \cosh ( \tilde { x } _ { i } ( T - t ) ) + \frac { \kappa _ { i } } { \tilde { x } _ { i } } \sinh ( \tilde { x } _ { i } ( T - t ) ) \right ] - \frac { 1 } { 2 } \kappa _ { i } ( T - t ) \\ & - \, \lambda _ { i } ^ { 2 } \frac { ( \kappa _ { i } \theta _ { i } ) ^ { 2 } } { \tilde { x } _ { i } ^ { 2 } } \left [ \frac { \sinh ( \tilde { x } _ { i } ( T - t ) ) } { \tilde { x } _ { i } \cosh ( \tilde { x } _ { i } ( T - t ) ) + \kappa _ { i } \sinh ( \tilde { x } _ { i } ( T - t ) ) } - ( T - t ) \right ] \\ & - \, \lambda _ { i } ^ { 2 } \frac { ( \kappa _ { i } \theta _ { i } ) ^ { 2 } \kappa _ { i } } { \tilde { x } _ { i } ^ { 2 } } \frac { \cosh ( \tilde { x } _ { i } ( T - t ) ) - 1 } { \tilde { x } _ { i } \cosh ( \tilde { x } _ { i } ( T - t ) ) + \kappa _ { i } \sinh ( \tilde { x } _ { i } ( T - t ) ) }, \\ \hbar { \hat { x } } _ { i } = \, \sqrt { \kappa _ { i } ^ { 2 } + \nu _ { i } ^ { 2 } \lambda _ { i } ^ { 2 } }. \, \text{We train with the forward process}$$

with ˆ κ i = √ κ 2 i + ν 2 i λ 2 i . We train with the forward process

$$\mathcal { X } _ { k + 1 } & = \ \mathcal { X } _ { k } + \Delta W _ { k }, \quad k = 0, \dots, N - 1, \ \mathcal { X } _ { 0 } \ = \ x _ { 0 }, \\ V _ { k + 1 } ^ { i } & = \ V _ { k } ^ { i } + \nu _ { i } \Delta B _ { k } ^ { i }, \quad k = 0, \dots, N - 1, \ V _ { 0 } ^ { i } \ = \ \theta _ { i },$$

with &lt; W,B i &gt; t = 0.

Merton Problem. We take η = 0 5, . λ = 0 6, . T = 1, N = 120, and σ v ( ) = e v . We plot the neural networks approximation of u, D x u, D 2 x u, α (in blue) together with their analytic values (in orange). For comparison with Figures 6 and 7, we report the error on the gradient and the initial control. In practice, after empirical tests, we choose p = 0 98 for the truncation. .

|                     |   Averaged value |   Standard deviation |   Theoretical value |   Relative error (%) |
|---------------------|------------------|----------------------|---------------------|----------------------|
| u (0 , x 0 = 1)     |         -0.50561 |              0.00029 |            -0.50662 |                 0.2  |
| D x u (0 , x 0 = 1) |          0.25081 |              0.00088 |             0.25331 |                 0.99 |
| α (0 , x 0 = 1)     |          0.83552 |              0.02371 |             0.80438 |                 3.87 |

Table 2: Estimate of the solution, its derivative and the optimal control at the initial time t = 0 in the Merton problem (3.3). Average and standard deviation observed over 10 independent runs are reported.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-de3de506-212a-4cac-9b18-e4e40adbf7a4.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-13db6526-a3c8-4c80-81b7-981f53fe9cd0.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-a89b8482-51d3-4a5d-a2af-d7077edaa0c2.jpg)

Y at date t = 0 5042. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-0ef8ac41-8062-4ca0-8616-edc4b616c901.jpg)

Z at date t = 0 5042 .

Γ at date t = 0 5042 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-9785cfaf-7671-4d4f-b6cf-60d028e92edb.jpg)

Y at date t = 0 0084. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-651f7994-4345-4dcb-b56c-6c01e7f9dc0c.jpg)

Z at date t = 0 0084 .

Γ at date t = 0 0084 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-3c41e2aa-6671-40b9-a779-14a8f8d7a3f6.jpg)

Figure 13: Estimates of the solution and its derivatives on the Merton problem (3.3) using 120 time steps.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-83567c9b-4720-4799-ab39-49b845284e0e.jpg)

α at date t = 0 5042. .

α at date t = 0 0084. .

Figure 14: Estimates of the optimal control α on the Merton problem (3.3).

One asset ( n = 1 ) in Scott volatility model. We take η = 0 5, . λ = 1 5, . θ = 0 4, . ν = 0 4, . κ = 1, ρ = -0 7. . For all tests we choose T = 1, N = 120, and σ v ( ) = e v . In practice, after empirical tests, we choose

p = 0 98 for the truncation. .

Table 3: Estimate of u (0 x , 0 = 1 , θ ) on the One Asset problem with stochastic volatility ( d = 2). Average and standard deviation observed over 10 independent runs are reported. The exact solution is -0 53609477. .

|   Averaged value |   Standard deviation |   Relative error (%) |
|------------------|----------------------|----------------------|
|         -0.53431 |               0.0007 |                 0.34 |

No Leverage in Scott model. In the case with one asset ( n = 1), we take η = 0 5, . λ = 1 5, . θ = 0 4, . ν = 0 2, . κ = 1. For all tests we choose T = 1, N = 120, and σ v ( ) = e v . In practice, after empirical tests, we choose p = 0 95 for the truncation. .

Table 4: Estimate of u (0 , x 0 = 1 , θ ) on the No Leverage problem (3.4). Average and standard deviation observed over 10 independent runs are reported.

|   Dimension d |   Averaged value |   Standard deviation |   Relative error (%) |   Theoretical solution |
|---------------|------------------|----------------------|----------------------|------------------------|
|             2 |         -0.4998  |              0.00073 |                 0.35 |              -0.501566 |
|             5 |         -0.43768 |              0.00137 |                 0.92 |              -0.441765 |
|             8 |         -0.3872  |              0.00363 |                 1.96 |              -0.394938 |
|            10 |         -0.2792  |              0.05734 |                 1.49 |              -0.275092 |

In the case with four assets ( n =4, d = 5), we take η = 0 5, . λ = 1 5 ( . 1 1 . 2 . 0 8 . ) , θ = 0 1 ( . 0 2 . 0 3 . 0 4 . ) , ν = 0 2 ( . 0 15 . 0 25 . 0 31 . ) , κ = 1 ( . 0 8 . 1 1 . 1 3 . ) .

- In the case with seven assets ( n = 7, d = 8) we take η = 0 5, . λ = 1 5 ( . 1 1 . 2 . 0 8 . 0 5 . 1 7 . 0 9 . ) ,
- θ = 0 1 ( . 0 2 . 0 3 . 0 4 . 0 25 . 0 15 . 0 18 . ) , ν = 0 2 ( . 0 15 . 0 25 . 0 31 . 0 4 . 0 35 . 0 22 . ) ,

κ

= 1

(

.

0 8

.

1 1

.

1 3

.

0 95

.

0 99

.

1 02

.

)

.

In the case with nine assets ( n =9, d = 10), we take η = 0 5, . λ = 1 5 ( . 1 1 . 2 . 0 8 . 0 5 . 1 7 . 0 9 . 1 . 0 9 . ) , θ = 0 1 ( . 0 2 . 0 3 . 0 4 . 0 25 . 0 15 . 0 18 . 0 08 . 0 91 . ) , ν = 0 2 ( . 0 15 . 0 25 . 0 31 . 0 4 . 0 35 . 0 22 . 0 4 . 0 15 . ) , κ = 1 ( . 0 8 . 1 1 . 1 3 . 0 95 . 0 99 . 1 02 . 1 06 . 1 6 . ) .

Hamilton-Jacobi-Bellman equation from portfolio optimization is a typical example of full-nonlinearity in the second order derivative, and the above results show that our algorithm performs quite well up to dimension d = 8, but gives a high variance in dimension d = 10.

Comparison with an implicit version of the scheme. As explained in Remark 2.2, an alternative option for the estimation of the Hessian is to approximate it by the automatic differentiation of the current neural network for the Z component. It corresponds to the replacement of D ˆ Z i +1 ( T ( X t i +1 )) by D Z T i ( ( X t i )); θ ) in (2.7). An additional change has to be made to the method for it to work. At the last optimization step (for time step t 0 = 0), we notice empirically that the variable Γ 0 is not able to properly learn the initial Hessian value at all. Therefore for this last step we use variables Y , Z 0 0 and an explicit estimation of the second order derivative given by D ˆ ( Z 1 T ( X t 1 )). We see in Table 5 that the results for the Merton problem are very similar to the ones from Table 2 for the splitting scheme but with a worse estimation of the Hessian and optimal control (the error is multiplied by around 1.5). When we tested this implicit scheme on the Monge Ampere problem we also faced computational problems during the optimization step of Tensorflow. The numerical computation of the gradient of the objective function for the backpropagation step, more precisely for the determinant part, often gives rise to matrix invertibility errors which stops the algorithm execution. For these two reasons, we focused our study on the explicit scheme.

Table 5: Estimate of the solution, its derivative and the optimal control at the initial time t = 0 in the Merton problem (3.3) with implicit estimation of the Hessian. Average and standard deviation (Std) observed over 10 independent runs are reported

|                     |   Average |     Std |   True value |   Relative error (%) |
|---------------------|-----------|---------|--------------|----------------------|
| u (0 , x 0 = 1)     |  -0.50572 | 0.00034 |     -0.50662 |                 0.18 |
| D x u (0 , x 0 = 1) |   0.25091 | 0.00067 |      0.25331 |                 0.95 |
| α (0 , x 0 = 1)     |   0.85254 | 0.01956 |      0.80438 |                 5.99 |

Comparison with the 2BSDE scheme of [BEJ19]. We conclude this paper with a comparison of our algorithm with the global scheme of [BEJ19], called Deep 2BDSE. The tests below concern the Merton

problem (3.3) but similar behavior happens on the other examples with stochastic volatilities. This scheme was implemented in the original paper only for small number of time steps (e.g. N = 30). Thus we tested this algorithm on two discretizations, respectively with N = 20 and N = 120 time steps, as shown in Figure 15, for T = 1 where we plotted the learning curve of the Deep BSDE method. These curves correspond to the values taken by the loss function during the gradient descent iterations. For this algorithm the loss function to minimize in the training of neural networks is defined as the mean L 2 error between the generated Y N value and the true terminal condition g X ( N ). We observe that for this choice of maturity T = 1 the loss function oscillates during the training process and does not vanish. As a consequence the Deep 2BSDE does not converge in this case. Even when decreasing the learning rate, we noticed that we cannot obtain the convergence of the scheme.

However, the Deep 2BSDE method does converge for small maturities T , as illustrated in Table 6 with T = 0 1 and different values for the number of time steps . N . Nevertheless, even if the value function is well approximated, the estimation of the gradient and control did not converge (the corresponding variance is very large), in comparison with our scheme whereas the gradient is very well approximated and the control is quite precise. We also have a much smaller variance in the results. Table 7 shows the results obtained by our method with T = 0 1 in order to compare it with the performance of [BEJ19]. It illustrates the limitations of the global . approach and justifies our introduction of a local method.

Figure 15: Learning curve in logarithmic scale for the scheme [BEJ19] on the Merton problem (3.3) with N = 20 times steps on the left and N = 120 time steps on the right. The maturity is T = 1 . 10000 gradient descent iterations were conducted.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-f86dbb91-c658-458e-a71d-15cf7df36011.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00412v3-figure-ec1d981e-10f2-4918-89d2-83d909943bd9.jpg)

Table 6: Estimate of the solution, its derivative and the optimal control at the initial time t = 0 in the Merton problem (3.3) with maturity T = 0 1 for the [BEJ19] scheme. . Average and standard deviation observed over 10 independent runs are reported.

|                     |   N |   Averaged value |   Standard deviation |   Theoretical value |   Relative error (%) |
|---------------------|-----|------------------|----------------------|---------------------|----------------------|
| u (0 , x 0 = 1)     |   5 |         -0.60667 |              0.01588 |            -0.59571 |                 1.84 |
| u (0 , x 0 = 1)     |  10 |         -0.59841 |              0.02892 |            -0.59571 |                 0.45 |
| u (0 , x 0 = 1)     |  20 |         -0.59316 |              0.04251 |            -0.59571 |                 0.43 |
| D x u (0 , x 0 = 1) |   5 |          0.09668 |              0.2563  |             0.29786 |                67.54 |
| D x u (0 , x 0 = 1) |  10 |          0.0381  |              0.4457  |             0.29786 |                93.36 |
| D x u (0 , x 0 = 1) |  20 |          0.07557 |              0.5503  |             0.29786 |                74.63 |
| α (0 , x 0 = 1)     |   5 |         -0.15243 |              0.61096 |             0.80438 |               118.95 |
| α (0 , x 0 = 1)     |  10 |          0.59971 |              1.97906 |             0.80438 |                25.44 |
| α (0 , x 0 = 1)     |  20 |          0.28385 |              0.43775 |             0.80438 |                64.71 |

Table 7: Estimate of the solution, its derivative and the optimal control at the initial time t = 0 in the Merton problem (3.3) with maturity T = 0 1 for our scheme. Average and standard deviation observed over 10 . independent runs are reported.

|                     |   N |   Averaged value |   Standard deviation |   Theoretical value |   Relative error (%) |
|---------------------|-----|------------------|----------------------|---------------------|----------------------|
| u (0 , x 0 = 1)     |   5 |         -0.59564 |              0.01136 |            -0.59571 |                0.01  |
| u (0 , x 0 = 1)     |  10 |         -0.5955  |              0.00037 |            -0.59571 |                0.04  |
| u (0 , x 0 = 1)     |  20 |         -0.59544 |              0.00054 |            -0.59571 |                0.04  |
| D x u (0 , x 0 = 1) |   5 |          0.29848 |              0.00044 |             0.29786 |                0.21  |
| D x u (0 , x 0 = 1) |  10 |          0.29842 |              0.00084 |             0.29786 |                0.19  |
| D x u (0 , x 0 = 1) |  20 |          0.29785 |              0.00054 |             0.29786 |                0.001 |
| α (0 , x 0 = 1)     |   5 |          0.82322 |              0.01014 |             0.80438 |                2.34  |
| α (0 , x 0 = 1)     |  10 |          0.85284 |              0.07565 |             0.80438 |                6.02  |
| α (0 , x 0 = 1)     |  20 |          0.84201 |              0.09892 |             0.80438 |                4.68  |

## References

| [AA+18]   | A. Al-Aradi et al. 'Solving Nonlinear and High-Dimensional Partial Differential Equations via Deep Learning'. In: arXiv:1811.08782 (2018).                                                                                                                       |
|-----------|------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [Aba+15]  | M. Abadi et al. TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems . Software available from tensorflow.org. 2015. url : https://www.tensorflow.org/ .                                                                                            |
| [Bec+19]  | C. Beck et al. 'Deep splitting method for parabolic PDEs'. In: arXiv:1907.03452 (2019).                                                                                                                                                                          |
| [BEJ19]   | C. Beck, W. E, and A. Jentzen. 'Machine Learning Approximation Algorithms for High- Dimensional Fully Nonlinear Partial Differential Equations and Second-order Backward Stochas- tic Differential Equations'. In: J. Nonlinear Sci. 29.4 (2019), pp. 1563-1619. |
| [Che+07]  | P. Cheridito et al. 'Second-order backward stochastic differential equations and fully nonlinear parabolic PDEs'. In: Communications on Pure and Applied Mathematics 60.7 (2007), pp. 1081- 1110.                                                                |
| [CWNMW19] | Quentin Chan-Wai-Nam, Joseph Mikael, and Xavier Warin. 'Machine learning for semi linear PDEs'. In: Journal of Scientific Computing 79.3 (2019), pp. 1667-1712.                                                                                                  |
| [DLM20]   | J. Darbon, G. Langlois, and T. Meng. 'Overcoming the curse of dimensionality for some Hamilton-Jacobi partial differential equations via neural network architectures'. In: Res Math Sci 7 (2020).                                                               |
| [EHJ17]   | W. E, J. Han, and A. Jentzen. 'Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations'. In: Com- munications in Mathematics and Statistics 5.4 (2017), pp. 349-380.  |
| [FTW11]   | A. Fahim, N. Touzi, and X. Warin. 'A probabilistic numerical method for fully nonlinear parabolic PDEs'. In: The Annals of Applied Probability (2011), pp. 1322-1364.                                                                                            |
| [HJE18]   | J. Han, A. Jentzen, and W. E. 'Solving high-dimensional partial differential equations using deep learning'. In: Proceedings of the National Academy of Sciences of the United States of America 115.34 (2018), pp. 8505-8510.                                   |
| [HPW20]   | Cˆme o Hur´, e Huyˆn e Pham, and Xavier Warin. 'Deep backward schemes for high-dimensional nonlinear PDEs'. In: Mathematics of Computation 89.324 (2020), pp. 1547-1579.                                                                                         |
| [Hut+18]  | M. Hutzenthaler et al. 'Overcoming the curse of dimensionality in the numerical approximation of semilinear parabolic partial differential equations'. In: arXiv:1807.01212 (2018).                                                                              |
| [KB14]    | D. P. Kingma and J. Ba. Adam: A Method for Stochastic Optimization . Published as a conference paper at the 3rd International Conference for Learning Representations, San Diego, 2015. 2014. url : http://arxiv.org/abs/1412.6980 .                             |
| [NR20]    | N. N¨sken u and L. Richter. 'Solving high-dimensional Hamilton-Jacobi-Bellman PDEs using neural networks: perspectives from the theory of controlled diffusions and measures on path space'. In: arXiv:2005.05409 (2020).                                        |
| [Pha09]   | H. Pham. Continuous-time Stochastic Control and Optimization with Financial Applications . Vol. 61. SMAP. Springer, 2009.                                                                                                                                        |
| [Pin99]   | A. Pinkus. 'Approximation theory of the MLP model in neural networks'. In: Acta numerica                                                                                                                                                                         |

| [PP90]   | E. Pardoux and S. Peng. 'Adapted solution of a backward stochastic differential equation'. In: Systems & Control Letters 14.1 (1990), pp. 55-61.                                 |
|----------|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [SS18]   | J. Sirignano and K. Spiliopoulos. 'DGM: A deep learning algorithm for solving partial differ- ential equations'. In: Journal of Computational Physics 375 (2018), pp. 1339-1364. |
| [SZ99]   | R. Sch¨bel o and J. Zhu. 'Stochastic volatility with an Ornstein-Uhlenbeck process and exten- sion'. In: Review of Finance 3.1 (1999), pp. 23-46.                                |
| [Tan13]  | X. Tan. 'A splitting method for fully nonlinear degenerate parabolic PDEs'. In: Electronic Journal of Probability 18 (2013).                                                     |
| [VSS18]  | M. Sabate Vidales, D. Siska, and L. Szpruch. 'Unbiased deep solvers for parametric PDEs'. In: arXiv:1810.05094v2 (2018).                                                         |
| [War18]  | X. Warin. 'Monte Carlo for high-dimensional degenerated Semi Linear and Full Non Linear PDEs'. In: arXiv:1805.05078 (2018).                                                      |
| [Zar01]  | T. Zariphopoulou. 'A solution approach to valuation with unhedgeable risks'. In: Finance and Stochastics 5 (2001), pp. 61-82.                                                    |
Huyen Pham Xavier Warin Maximilien Germain
Neural networks-based backward scheme for fully nonlinear PDEs
We propose a numerical method for solving high dimensional fully nonlinear
partial differential equations (PDEs). Our algorithm estimates simultaneously
by backward time induction the solution and its gradient by multi-layer neural
networks, while the Hessian is approximated by automatic differentiation of the
gradient at previous step. This methodology extends to the fully nonlinear case
the approach recently proposed in \cite{HPW19} for semi-linear PDEs. Numerical
tests illustrate the performance and accuracy of our method on several examples
in high dimension with nonlinearity on the Hessian term including a linear
quadratic control problem with control on the diffusion coefficient,
Monge-Amp{\`e}re equation and Hamilton-Jacobi-Bellman equation in portfolio
optimization.