Neural networks-based backward scheme for fully nonlinear PDEs

TL;DR

This paper introduces a neural network-based backward scheme for efficiently solving high-dimensional fully nonlinear PDEs, demonstrating high accuracy and applicability to complex problems like control and Monge-Ampère equations.

Contribution

It extends neural network methods to fully nonlinear PDEs, enabling simultaneous estimation of solutions, gradients, and Hessians through automatic differentiation.

Findings

Accurate solutions for high-dimensional nonlinear PDEs.

Effective handling of nonlinearity in Hessian terms.

Successful application to control and Monge-Ampère equations.

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 \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.

Tables

Table 1. Table 1: Estimate of u ​ ( 0 , x 0 = 1 5 ) 𝑢 0 subscript 𝑥 0 subscript 1 5 u(0,x_{0}=1_{5}) on the Monge Ampere problem ( 3.1 ) with N = 120 𝑁 120 N=120 . Average and standard deviation observed over 10 independent runs are reported.

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 2. Table 2: Estimate of the solution, its derivative and the optimal control at the initial time t = 0 𝑡 0 t=0 in the Merton problem ( 3.10 ). Average and standard deviation observed over 10 independent runs are reported.

	Averaged value	Standard deviation	Theoretical value	Relative error (%)
$u\left(0,x_0=1\right)$	-0.50561	0.00029	-0.50662	0.20
$D_{x}u\left(0,x_0=1\right)$	0.25081	0.00088	0.25331	0.99
$\alpha \left(0,x_0=1\right)$	0.83552	0.02371	0.80438	3.87

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

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

Table 4. Table 4: Estimate of u ( 0 , x 0 = 1 , θ ) u(0,\mathrm{x}_{0}=1,\theta) on the No Leverage problem ( 3.26 ). Average and standard deviation observed over 10 independent runs are reported.

Dimension $d$	Averaged value	Standard deviation	Relative error (%)	Theoretical solution
2	-0.49980	0.00073	0.35	-0.501566
5	-0.43768	0.00137	0.92	-0.441765
8	-0.38720	0.00363	1.96	-0.394938
10	-0.27920	0.05734	1.49	-0.275092

Table 5. Table 5: Estimate of the solution, its derivative and the optimal control at the initial time t = 0 𝑡 0 t=0 in the Merton problem ( 3.10 ) 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\left(0,x_0=1\right)$	-0.50572	0.00034	-0.50662	0.18
$D_{x}u\left(0,x_0=1\right)$	0.25091	0.00067	0.25331	0.95
$\alpha \left(0,x_0=1\right)$	0.85254	0.01956	0.80438	5.99

Table 6. Table 6: Estimate of the solution, its derivative and the optimal control at the initial time t = 0 𝑡 0 t=0 in the Merton problem ( 3.10 ) with maturity T = 0.1 𝑇 0.1 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\left(0,x_0=1\right)$	5	-0.60667	0.01588	-0.59571	1.84
$u\left(0,x_0=1\right)$	10	-0.59841	0.02892	-0.59571	0.45
$u\left(0,x_0=1\right)$	20	-0.59316	0.04251	-0.59571	0.43
$D_{x}u\left(0,x_0=1\right)$	5	0.09668	0.25630	0.29786	67.54
$D_{x}u\left(0,x_0=1\right)$	10	0.03810	0.44570	0.29786	93.36
$D_{x}u\left(0,x_0=1\right)$	20	0.07557	0.55030	0.29786	74.63
$\alpha \left(0,x_0=1\right)$	5	-0.15243	0.61096	0.80438	118.95
$\alpha \left(0,x_0=1\right)$	10	0.59971	1.97906	0.80438	25.44
$\alpha \left(0,x_0=1\right)$	20	0.28385	0.43775	0.80438	64.71

Table 7. Table 7: Estimate of the solution, its derivative and the optimal control at the initial time t = 0 𝑡 0 t=0 in the Merton problem ( 3.10 ) with maturity T = 0.1 𝑇 0.1 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\left(0,x_0=1\right)$	5	-0.59564	0.01136	-0.59571	0.01
$u\left(0,x_0=1\right)$	10	-0.59550	0.00037	-0.59571	0.04
$u\left(0,x_0=1\right)$	20	-0.59544	0.00054	-0.59571	0.04
$D_{x}u\left(0,x_0=1\right)$	5	0.29848	0.00044	0.29786	0.21
$D_{x}u\left(0,x_0=1\right)$	10	0.29842	0.00084	0.29786	0.19
$D_{x}u\left(0,x_0=1\right)$	20	0.29785	0.00054	0.29786	0.001
$\alpha \left(0,x_0=1\right)$	5	0.82322	0.01014	0.80438	2.34
$\alpha \left(0,x_0=1\right)$	10	0.85284	0.07565	0.80438	6.02
$\alpha \left(0,x_0=1\right)$	20	0.84201	0.09892	0.80438	4.68

Equations

$$\{\begin{array}{rl}{\partial }_{t}u+f(.,.,u,D_{x}u,D_x^{2}u)& =0,\text{mbox}on[0,T)\times R^d,\\ u(T,.)& =g,\text{mbox}on{R}^d,\end{array}$$

$$x\in R^d$$

$$A_{\ell }(x)$$

$$X_t$$

$$Y_t=u(t,X_t),Z_t=D_{x}u(t,X_t),{\Gamma }_t=D_x^{2}u(t,X_t),0\le t\le T.$$

$$Y_t$$

$$-{\int }_t^T{\sigma }^{⊺}(s,X_s)Z_s.d{W}_s,0\le t\le T.$$

$$X_{t_{i+1}}$$

$$Y_{t_i}$$

$$Y_{t_{i+1}}$$

$$F(t,x,y,z,\gamma ,h,\Delta )$$

$$\tilde{f}(t,x,y,z,\gamma )$$

$${\hat{L}}_i(\theta )$$

$$u(t_{i+1},X_{t_{i+1}})$$

$${\hat{L}}_i(\theta )$$

$$\displaystyle\;\;\;-\;\Delta t_{i}\big{[}\tilde{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}}))-\tilde{f}(t_{i},X_{t_{i}},{\cal U}_{i}(X_{t_{i}};\theta),{\cal Z}_{i}(X_{t_{i}};\theta),D\hat{\cal Z}_{i+1}({\cal T}(X_{t_{i+1}})))\big{]}\Big{|}^{2}.$$

$${\hat{L}}_i(\theta )$$

$$\left\{\begin{aligned} \hat{L}_{N}(\theta)&:=\;\mathbb{E}\Big{|}{\cal U}_{N}(X_{t_{N}};\theta)-g(X_{t_{N}})\Big{|}^{2}+\frac{\Delta t_{N-1}}{d}\mathbb{E}\Big{|}{\cal Z}_{N}(X_{t_{N}};\theta)-Dg(X_{t_{N}})\Big{|}^{2}\\ \theta_{N}^{*}&\in{\rm arg}\min_{\theta\in\operatorname{\mathbb{R}}^{N_{m}}}\hat{L}_{N}(\theta).\end{aligned}\right.$$

$$\left\{\begin{aligned} \hat{L}_{i}(\theta)&:=\;\mathbb{E}\Big{|}\widehat{\cal U}_{i+1}(X_{t_{i+1}})\;-\;F(t_{i},X_{t_{i}},{\cal U}_{i}(X_{t_{i}};\theta),{\cal Z}_{i}(X_{t_{i}};\theta),D\hat{\cal Z}_{i+1}({\cal T}(X_{t_{i+1}})),\Delta t_{i},\Delta W_{t_{i}})\Big{|}^{2}\\ \theta_{i}^{*}&\in{\rm arg}\min_{\theta\in\operatorname{\mathbb{R}}^{N_{m}}}\hat{L}_{i}(\theta).\end{aligned}\right.$$

$$f(t,x,y,z,\gamma )$$

$$u(t,x)$$

$$\sigma =$$

$$\displaystyle{\cal T}_{p}(X_{t}^{0,x})=\min\big{\{}\max[x-\sigma\sqrt{t}\phi_{p},X_{t}^{0,x}],x+\sigma\sqrt{t}\phi_{p}\big{\}},$$

$$d{\mathcal{X}}_t$$

$$J(\alpha )=$$

$$\displaystyle\frac{\partial u}{\partial t}+\inf_{a\in\operatorname{\mathbb{R}}}\big{[}(Ax+Ba).D_{x}u+\frac{a^{2}}{2}{\rm tr}(DD^{\scriptscriptstyle{\intercal}}D_{x}^{2}u)+x^{\scriptscriptstyle{\intercal}}Qx+Na^{2}\big{]}$$

$$u(T,x)$$

$$f(t,x,y,z,\gamma )$$

$$u(t,x)$$

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

$$\sigma =$$

$$A={\mathbf{I}}_d,B=D=1{\mathrm{I}}_d,$$

Full text

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

[TABLE]

with a non-linearity in the solution, its gradient $D_{x}u$ and its hessian $D_{x}^{2}u$ via the function $f(t,x,y,z,\gamma)$ defined on $[0,T]\times\operatorname{\mathbb{R}}^{d}\times\operatorname{\mathbb{R}}\times\operatorname{\mathbb{R}}^{d}\times\mathbb{S}^{d}$ (where $\mathbb{S}^{d}$ is the set of symmetric $d\times 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 $\gamma$. 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.

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_{x}^{2}u$ at $t_{n}$ is approximated by a neural network minimizing some local $L_{2}$ criterion associated to a BSDE involving $D_{x}u$ at date $t_{n+1}$ and $D_{x}^{2}u$. 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_{x}^{2}u$ 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_{x}^{2}u$ 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 non-linearity of type $uD_{x}^{2}u$. 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ère equation, and finally, we provide numerical tests for the solution to fully nonlinear Hamilton-Jacobi-Bellman equation, with non-linearities of the form $|D_{x}u|^{2}/D_{x}^{2}u$, 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]\times\operatorname{\mathbb{R}}^{d}$ $\mapsto$ $\operatorname{\mathbb{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_{x}u$, 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$ $\in$ $\mathbb{N}\setminus\{1,2\}$ with $m_{\ell}$, $\ell$ $=$ $0,\ldots,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_{\ell}$ $=$ $m$, $\ell$ $=$ $1,\ldots,L-1$.

A feedforward neural network is a function from $\operatorname{\mathbb{R}}^{d_{0}}$ to $\operatorname{\mathbb{R}}^{d_{1}}$ defined as the composition

[TABLE]

Here $A_{\ell}$, $\ell$ $=$ $1,\ldots,L$ are affine transformations: $A_{1}$ maps from $\operatorname{\mathbb{R}}^{d_{0}}$ to $\operatorname{\mathbb{R}}^{m}$, $A_{2},\ldots,A_{L-1}$ map from $\operatorname{\mathbb{R}}^{m}$ to $\operatorname{\mathbb{R}}^{m}$, and $A_{L}$ maps from $\operatorname{\mathbb{R}}^{m}$ to $\operatorname{\mathbb{R}}^{d_{1}}$, represented by

[TABLE]

for a matrix ${\cal W}_{\ell}$ called weight, and a vector $\beta_{\ell}$ called bias term, $\varrho$ $:$ $\operatorname{\mathbb{R}}$ $\rightarrow$ $\operatorname{\mathbb{R}}$ is a nonlinear function, called activation function, and applied component-wise on the outputs of $A_{\ell}$, i.e., $\varrho(x_{1},\ldots,x_{m})$ $=$ $(\varrho(x_{1}),\ldots,\varrho(x_{m}))$. Standard examples of activation functions are the sigmoid, the ReLu, the Elu, $\tanh$.

All these matrices ${\cal W}_{\ell}$ and vectors $\beta_{\ell}$, $\ell$ $=$ $1,\ldots,L$, are the parameters of the neural network, and can be identified with an element $\theta$ $\in$ $\operatorname{\mathbb{R}}^{N_{m}}$, where $N_{m}$ $=$ $\sum_{\ell=0}^{L-1}m_{\ell}(1+m_{\ell+1})$ $=$ $d_{0}(1+m)+m(1+m)(L-2)+m(1+d_{1})$ is the number of parameters. We denote by ${\cal N}_{d_{0},d_{1},L,m}$ the set of all functions generated by (2.1) for $\theta\in\operatorname{\mathbb{R}}^{N_{m}}$.

2.2 Forward-backward representation

Let us introduce a forward diffusion process

[TABLE]

where $\mu$ is a function defined on $[0,T]\times\operatorname{\mathbb{R}}^{d}$ with values in $\operatorname{\mathbb{R}}^{d}$, $\sigma$ is a function defined on $[0,T]\times\operatorname{\mathbb{R}}^{d}$ with values in $\mathbb{M}^{d}$ the set of $d\times d$ matrices, and $W$ a $d$-dimensional Brownian motion on some probability space $(\Omega,{\cal F},\mathbb{P})$ equipped with a filtration $\mathbb{F}$ $=$ $({\cal F}_{t})_{0\leq t\leq 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 $\mu$ and $\sigma$, see Remark 2.3.

Let us next denote by $(Y,Z,\Gamma)$ the triple of $\mathbb{F}$-adapted processes valued in $\operatorname{\mathbb{R}}\times\operatorname{\mathbb{R}}^{d}\times\mathbb{S}^{d}$, defined by

[TABLE]

By Itô’s formula applied to $u(t,X_{t})$, and since $u$ is solution to (1.1), we see that $(Y,Z,\Gamma)$ satisfies the backward equation:

[TABLE]

Remark 2.1

This BSDE does not uniquely characterize a triple $(Y,Z,\Gamma$) contrarily to the semilinear case (without a non-linearity with respect to $\Gamma$) 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 $\Gamma$ 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.

*** $\Box$*

2.3 Algorithm

We now provide a numerical approximation of the forward backward system (2.3)-(2.5), and consequently of the solution $u$ (as well as its gradient $D_{x}u$) to the PDE (1.1).

We start from a time grid $\pi$ $=$ $\{t_{i},i=0,\ldots,N\}$ of $[0,T]$, with $t_{0}$ $=$ [math] $<$ $t_{1}$ $<$ $\ldots$ $<$ $t_{N}$ $=$ $T$, and time steps $\Delta t_{i}$ $:=$ $t_{i+1}-t_{i}$, $i$ $=$ $0,\ldots,N-1$. The time discretization of the forward process $X$ on $\pi$ is then equal (typically when $\mu$ and $\sigma$ are constants) or approximated by an Euler scheme:

[TABLE]

where we set $\Delta 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.5) is approximated by the time discretized scheme

[TABLE]

that is written in forward form as

[TABLE]

with

[TABLE]

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_{x}u(t_{i},.)$, by neural networks $x$ $\in$ $\operatorname{\mathbb{R}}^{d}$ $\mapsto$ $({\cal U}_{i}(x;\theta),{\cal Z}_{i}(x;\theta))$ with parameter $\theta$ that are learned optimally by backward induction: suppose that $\hat{\cal U}_{i+1}$ $:=$ ${\cal U}_{i+1}(.;\theta_{i+1}^{*})$, $\hat{\cal Z}_{i+1}$ $:=$ ${\cal Z}_{i+1}(.;\theta_{i+1}^{*})$ is an approximation of $u(t_{i+1},.)$ and $D_{x}u(t_{i+1},.)$ at time $t_{i+1}$, then $\theta_{i}^{*}$ is computed from the minimization of the quadratic loss function:

[TABLE]

where ${\cal T}$ is a truncation operator such that ${\cal T}(X)$ is bounded for example by a quantile of the diffusion process and $D\hat{\cal Z}_{i+1}$ stands for the automatic differentiation of $\hat{\cal 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.4) and (2.6), the solution $u$ to (1.1) should approximately satisfy

[TABLE]

Suppose that at time $t_{i+1}$, $\hat{\cal U}_{i+1}$ is an estimation of $u(t_{i+1},.)$. Recalling the expression of $F$ in (2.7), the quadratic loss function at time $t_{i}$ is then approximately equal to

[TABLE]

By assuming that $\tilde{f}$ has small non-linearities in its arguments $(y,z,\gamma)$, say Lipschitz, possibly with a suitable choice of $\mu$, $\sigma$, the loss function is thus approximately equal to

[TABLE]

Therefore, by minimizing over $\theta$ this quadratic loss function, via stochastic gradient descent (SGD) based on simulations of $(X_{t_{i}},X_{t_{i+1}},\Delta W_{t_{i}})$ (called training data in the machine learning language), one expects the neural networks ${\cal U}_{i}$ and ${\cal Z}_{i}$ to learn/approximate better and better the functions $u(t_{i},.)$ and $D_{x}u(t_{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$.

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 ${\cal U}_{i}$ and estimate $\Gamma_{i}$ as the iterated automatic differentiation $D_{x}^{2}{\cal 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 ${\cal U},{\cal 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 ${\cal 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}$. $\Box$**

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.10), 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. * $\Box$*

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 Monge-Ampère 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.

$\bullet$ 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 $\mu(.)$ 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 $\sigma$ is also important: large $\sigma$ 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 $\sigma$ 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 $\sigma$ 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 $\sigma$.

$\bullet$ Parameters of truncation: Given the training simulations $X$, we choose a truncation operator ${\cal T}_{p}$ indexed by a parameter $p$ close to $1$, so that ${\cal T}_{p}(X_{t})$ corresponds to a truncation of $X_{t}$ at a given quantile $\phi_{p}$. In the numerical tests, we shall vary $p$ between $0.95$ and $0.999$.

$\bullet$ 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.10) 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<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 $1E-2$ and at the following optimization steps, we start with a learning rate equal to $1E-3$.

During time resolution, it is far more effective to initialize the solution of equations (2.10) with the solution $({\cal U},{\cal 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_{x}^{2}u$

We consider a generator in the form

[TABLE]

and $g(x)$ $=$ \tanh{\big{(}\frac{\sum_{i=1}^{d}x_{i}}{\sqrt{d}}\big{)}}, so that an analytical solution is available:

[TABLE]

We fix the horizon $T$ $=$ $1$, and choose to evaluate the solution at $t=0$ and $x=0.5\frac{1{\rm I}_{d}}{\sqrt{d}}$ (here $1{\rm I}_{d}$ denotes the vector in $\operatorname{\mathbb{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 ${\bf I}_{d}$ is the identity $d\times d$-matrix)

[TABLE]

and here the truncation operator is chosen equal to

[TABLE]

where $\phi_{p}$ $=$ $\mathcal{N}^{-1}(p)$, with $\mathcal{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_{x}u$ and $D_{x}^{2}u$ are well approximated by the resolution method.

On Figure 2, we check the convergence, for different values of $\hat{\sigma}$ of both the solution $u$ and its derivative at point $x$ and date [math]. Standard deviation of the function value is very low and the standard deviation of the derivative still being low.

As the dimension increases, we have to increase the value of $\hat{\sigma}$ of the forward process. In dimension 3, the value $\hat{\sigma}=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 $\hat{\sigma}=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.

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 ${\cal X}$ $=$ ${\cal X}^{\alpha}$ with dynamics in $\operatorname{\mathbb{R}}^{d}$ according to

[TABLE]

where $W$ is a real Brownian motion, the control process $\alpha$ is valued in $\operatorname{\mathbb{R}}$, and the constant coefficients $A\in\mathbb{M}^{d}$, $B\in\operatorname{\mathbb{R}}^{d}$, $D\in\operatorname{\mathbb{R}}^{d}$. The quadratic cost functional to be minimized is

[TABLE]

where $P$, $Q$ are non negative $d\times d$ symmetric matrices and $N\in\operatorname{\mathbb{R}}$ is strictly positive.

The Bellman equation associated to this stochastic control problem is:

[TABLE]

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

[TABLE]

An explicit solution to this PDE is given by

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

and we want to estimate the solution at $x=1{\rm I}_{d}$.

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

[TABLE]

where $\phi_{p}$ $=$ $\mathcal{N}^{-1}(p)$, $\hat{A}$ is a vector so that $\hat{A}_{i}=A_{ii}$, $i=1,...,d$, $\hat{1}$ is a unit vector, and the square root is taken componentwise.

On Figure 5 we give the solution of the PDE with $d=1$ using $\hat{\sigma}=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 $\Gamma$ value at date $t=0.5$. The precision remains good for $\Gamma$ close to $t=0$ and very good for $u$ and $D_{x}u$.

On Figure 6, we give the results obtained in dimension $d=1$ by varying $\hat{\sigma}$. For a value of $\hat{\sigma}=2$, the standard deviation of the result becomes far higher than with $\hat{\sigma}=0.5$ or $1.$

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.

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).

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$.

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

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 $\hat{\sigma}=5$ we had to take 300 neurons to get very accurate results.

3.4 Monge-Ampère equation

Let us consider the parabolic Monge-Ampère equation

[TABLE]

where ${\rm det}(D_{x}^{2}u)$ is the determinant of the Hessian matrix $D_{x}^{2}u$. It is in the form (1.1) with

[TABLE]

We test our algorithm by choosing a $C^{2}$ function $g$, then compute $G$ $=$ ${\rm det}(D_{x}^{2}g)$, and set $h$ $:=$ $G-1$. Then, by construction, the function

[TABLE]

is solution to the Monge-Ampère equation (3.1). We choose $g(x)=\cos(\sum_{i=1}^{d}x_{i}/\sqrt{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ère 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$.

3.5 Portfolio selection

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

[TABLE]

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

[TABLE]

with $\kappa_{i},\theta_{i},\nu_{i}$ $>$ [math], and $B$ $=$ $(B^{1},\ldots,B^{n})$ a $n$-dimensional Brownian motion, s.t. $d<W^{i},B^{j}>$ $=$ $\delta_{ij}\rho_{ij}dt$, with $\rho_{i}$ $:=$ $\rho_{ii}$ $\in$ $(-1,1)$. An agent can invest at any time an amount $\alpha_{t}$ $=$ $(\alpha_{t}^{1},\ldots,\alpha_{t}^{n})$ in the stocks, which generates a wealth process ${\cal X}$ $=$ ${\cal X}^{\alpha}$ governed by

[TABLE]

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

[TABLE]

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)\times\operatorname{\mathbb{R}}\times\operatorname{\mathbb{R}}^{n}$:

[TABLE]

with a Sharpe ratio $R(v)$ $:=$ $|\lambda(v)|^{2}$, for $v$ $=$ $(v_{1},\ldots,v_{n})$ $\in$ $(0,\infty)^{n}$. The optimal portfolio strategy is then given in feedback form by $\alpha_{t}^{*}$ $=$ $\hat{a}(t,{\cal X}_{t}^{*},V_{t})$, where $\hat{a}$ $=$ $(\hat{a}_{1},\ldots,\hat{a}_{n})$ is given by

[TABLE]

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

[TABLE]

for $x$ $=$ $(\mathrm{x},v)$ $\in$ $\operatorname{\mathbb{R}}^{n+1}$, $z$ $=$ $(z_{0},\ldots,z_{n})$ $\in$ $\operatorname{\mathbb{R}}^{n+1}$, $\gamma$ $=$ $(\gamma_{ij})_{0\leq i,j\leq n}$ $\in$ $\mathbb{S}^{n+1}$, and displays a high non-linearity in the Hessian argument $\gamma$.

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

[TABLE]

where $\phi_{p}$ $=$ $\mathcal{N}^{-1}(p)$, $\mathcal{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(-\eta x)$, with $\eta$ $>$ [math], 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 $\sigma$ and the risk premium $\lambda$ are constant, so that the generator of Bellman equation reduces to

[TABLE]

with explicit solution given by:

[TABLE]

We train with the forward process

[TABLE]

• (2)

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

[TABLE]

where $\hat{V}_{s}^{t,v}$ is the solution to the modified O.U. model

[TABLE]

We test our algorithm with $\lambda(v)$ $=$ $\lambda v$, $\lambda$ $>$ [math], for which we have an explicit solution:

[TABLE]

where $(\phi,\psi,\chi)$ are solutions of the Riccati system of ODEs:

[TABLE]

with $\bar{\kappa}$ $=$ $\kappa+\rho\nu\lambda$, and explicitly given by (see e.g. Appendix in [SZ99])

[TABLE]

with $\hat{\kappa}$ $=$ $\sqrt{\kappa^{2}+2\rho\nu\lambda\kappa+\gamma^{2}\lambda^{2}}$. We train with the forward process

[TABLE]

• (3)

No leverage effect, i.e., $\rho_{i}$ $=$ [math], $i$ $=$ $1,\ldots,n$. In this case, there is a quasi-explicit solution given by

[TABLE]

where $V^{t,v}$ is the solution to (3.5), starting from $v$ at time $t$. We test our algorithm with $\lambda_{i}(v)$ $=$ $\lambda_{i}v_{i}$, $\lambda_{i}$ $>$ [math], $i$ $=$ $1,\ldots,n$, $v$ $=$ $(v_{1},\ldots,v_{n})$, for which we have an explicit solution given by

[TABLE]

with $\hat{\kappa}_{i}$ $=$ $\sqrt{\kappa_{i}^{2}+\nu_{i}^{2}\lambda_{i}^{2}}$. We train with the forward process

[TABLE]

with $<W,B^{i}>_{t}$ $=$ [math].

Merton Problem. We take $\eta=0.5$, $\lambda=0.6$, $T$ $=$ $1$, $N=120$, and $\sigma(v)=e^{v}$. We plot the neural networks approximation of $u,D_{x}u,D^{2}_{x}u,\alpha$ (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.

One asset ($n$ $=$ $1$) in Scott volatility model. We take $\eta=0.5$, $\lambda=1.5$, $\theta=0.4$, $\nu=0.4$, $\kappa=1$, $\rho=-0.7$. For all tests we choose $T$ $=$ $1$, $N=120$, and $\sigma(v)=e^{v}$. In practice, after empirical tests, we choose $p=0.98$ for the truncation.

No Leverage in Scott model. In the case with one asset ($n$ $=$ $1$), we take $\eta=0.5$, $\lambda=1.5$, $\theta=0.4$, $\nu=0.2$, $\kappa=1$. For all tests we choose $T$ $=$ $1$, $N=120$, and $\sigma(v)=e^{v}$. In practice, after empirical tests, we choose $p=0.95$ for the truncation.

In the case with four assets ($n$ $=$ $4$, $d=5$), we take $\eta=0.5$, $\lambda=\begin{pmatrix}1.5&1.1&2.&0.8\end{pmatrix}$, $\theta=\begin{pmatrix}0.1&0.2&0.3&0.4\end{pmatrix}$, $\nu=\begin{pmatrix}0.2&0.15&0.25&0.31\end{pmatrix}$, $\kappa=\begin{pmatrix}1.&0.8&1.1&1.3\end{pmatrix}$.

In the case with seven assets ($n$ $=$ $7$, $d=8$) we take $\eta=0.5$, $\lambda=\begin{pmatrix}1.5&1.1&2.&0.8&0.5&1.7&0.9\end{pmatrix}$,

$\theta=\begin{pmatrix}0.1&0.2&0.3&0.4&0.25&0.15&0.18\end{pmatrix}$, $\nu=\begin{pmatrix}0.2&0.15&0.25&0.31&0.4&0.35&0.22\end{pmatrix}$,

$\kappa=\begin{pmatrix}1.&0.8&1.1&1.3&0.95&0.99&1.02\end{pmatrix}$.

In the case with nine assets ($n$ $=$ $9$, $d=10$), we take $\eta=0.5$, $\lambda=\begin{pmatrix}1.5&1.1&2.&0.8&0.5&1.7&0.9&1.&0.9\end{pmatrix}$, $\theta=\begin{pmatrix}0.1&0.2&0.3&0.4&0.25&0.15&0.18&0.08&0.91\end{pmatrix}$, $\nu=\begin{pmatrix}0.2&0.15&0.25&0.31&0.4&0.35&0.22&0.4&0.15\end{pmatrix}$, $\kappa=\begin{pmatrix}1.&0.8&1.1&1.3&0.95&0.99&1.02&1.06&1.6\end{pmatrix}$.

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\hat{\cal Z}_{i+1}({\cal T}(X_{t_{i+1}}))$ by $D{\cal Z}_{i}({\cal T}(X_{t_{i}}));\theta)$ in (2.10). 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 $\Gamma_{0}$ is not able to properly learn the initial Hessian value at all. Therefore for this last step we use variables $Y_{0},Z_{0}$ and an explicit estimation of the second order derivative given by $D\hat{\cal Z}_{1}({\cal 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.

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.10) 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.