Deep Learning-Based Numerical Methods for High-Dimensional Parabolic Partial Differential Equations and Backward Stochastic Differential Equations

Abstract

We study a new algorithm for solving parabolic partial differential equations (PDEs) and backward stochastic differential equations (BSDEs) in high dimension, which is based on an analogy between the BSDE and reinforcement learning with the gradient of the solution playing the role of the policy function, and the loss function given by the error between the prescribed terminal condition and the solution of the BSDE. The policy function is then approximated by a neural network, as is done in deep reinforcement learning. Numerical results using TensorFlow illustrate the efficiency and accuracy of the studied algorithm for several 100-dimensional nonlinear PDEs from physics and finance such as the Allen–Cahn equation, the Hamilton–Jacobi–Bellman equation, and a nonlinear pricing model for financial derivatives.

Appendix A: Special Cases of the Proposed Algorithm

In this section, we illustrate the general algorithm in Subsect. 3.2 in several special cases. More specifically, in Subsects. 5.1 and 5.2, we provide special choices for the functions $\psi _m$, $m \in {\mathbb {N}}$, and $\Psi _m$, $m \in {\mathbb {N}}$, employed in (3.14), and in Subsects. 5.3 and 5.4, we provide special choices for the function $\Upsilon$ in (3.9).

1.1 Stochastic Gradient Descent (SGD)

Example 5.1

Assume the setting in Subsect. 3.2, let $( \gamma _m )_{ m \in {\mathbb {N}}} \subseteq (0,\infty )$, and assume for all $m \in {\mathbb {N}}$, $x \in {\mathbb {R}}^{ \varrho }$, $( \varphi _j )_{ j \in {\mathbb {N}}} \in ( {\mathbb {R}}^{ \rho } )^{ {\mathbb {N}}}$ that

$$\begin{aligned} \varrho = \rho , \qquad \Psi _m( x, ( \varphi _j )_{ j \in {\mathbb {N}}} ) = \varphi _1 , \qquad \text {and} \qquad \psi _m( x ) = \gamma _m x .\nonumber \\ \end{aligned}$$

(5.1)

Then it holds for all $m \in {\mathbb {N}}$ that

$$\begin{aligned} \Theta _{ m } = \Theta _{ m - 1 } - \gamma _{ m } \Phi ^{ m - 1, 1 }_{ \mathbb {S}_m }( \Theta _{ m - 1 } ) . \end{aligned}$$

(5.2)

1.2 Adaptive Moment Estimation (Adam) with Mini-Batches

In this subsection, we illustrate how the so-called Adam optimizer (see [25]) can be employed in conjunction with the deep BSDE method in Subsect. 3.2 (cf. also Subsect. 4.1 above).

Example 5.2

Assume the setting in Subsect. 3.2, assume that $\varrho = 2 \rho$, let ${\text {Pow}}_r :{\mathbb {R}}^{ \rho } \rightarrow {\mathbb {R}}^{ \rho }$, $r \in (0,\infty )$, be the functions which satisfy for all $r \in (0,\infty )$, $x = ( x_1, \dots , x_{ \rho } ) \in {\mathbb {R}}^{ \rho }$ that

$$\begin{aligned} {\text {Pow}}_{ r }( x ) = ( | x_1 |^r, \dots , | x_{ \rho } |^r ) , \end{aligned}$$

(5.3)

let $\varepsilon \in (0,\infty )$, $( \gamma _m )_{ m \in {\mathbb {N}}} \subseteq (0,\infty )$, $( J_m )_{ m \in {\mathbb {N}}_0 } \subseteq {\mathbb {N}}$, $\mathbb {X}, \mathbb {Y} \in (0,1)$, let $\mathbf{m} = ( \mathbf{m}^{ (1) } , \dots , \mathbf{m}^{ ( \rho ) } ) :$ ${\mathbb {N}}_0 \times \Omega \rightarrow {\mathbb {R}}^{ \rho }$ and $\mathbb {M} = ( \mathbb {M}^{ (1) } , \dots \mathbb {M}^{ ( \rho ) } ) :{\mathbb {N}}_0 \times \Omega \rightarrow {\mathbb {R}}^{ \rho }$ be the stochastic processes which satisfy for all $m \in {\mathbb {N}}_0$ that $\Xi _m = ( \mathbf{m}_m^{ (1) }, \dots , \mathbf{m}^{ (\rho ) }_m , \mathbb {M}_m^{ (1) }, \dots , \mathbb {M}_m^{ (\rho ) } )$, and assume for all $m \in {\mathbb {N}}$, $x = ( x_1, \dots , x_{ \rho } ) , y = ( y_1, \dots , y_{ \rho } ) \in {\mathbb {R}}^{ \rho }$, $( \varphi _j )_{ j \in {\mathbb {N}}} \in ( {\mathbb {R}}^{ \rho } )^{ {\mathbb {N}}}$ that

$$\begin{aligned}&\Psi _m( x, y, ( \varphi _j )_{ j \in {\mathbb {N}}} ) = \big ( \mathbb {X} x + ( 1 - \mathbb {X} ) \big ( \tfrac{ 1 }{ J_m } \textstyle \sum _{ j = 1 }^{ J_m } \varphi _j \big ) , \mathbb {Y} y + ( 1 - \mathbb {Y} ) \nonumber \\&\quad {\text {Pow}}_2\big ( \frac{ 1 }{ J_m } \textstyle \sum _{ j = 1 }^{ J_m } \varphi _j \big ) \big ) \end{aligned}$$

(5.4)

and

$$\begin{aligned} \psi _m( x, y )\! =\! \left( \left[ \varepsilon + \tfrac{ \sqrt{ | y_1 | } }{ \sqrt{ 1 - \mathbb {Y}^m } } \right] ^{ - 1 } \frac{ \gamma _m x_1 }{ ( 1 - \mathbb {X}^m ) } , \dots , \left[ \varepsilon \!+\! \tfrac{ \sqrt{ | y_{ \rho } | } }{ \sqrt{ 1 \!-\! \mathbb {Y}^m } } \right] ^{ - 1 } \frac{ \gamma _m x_{ \rho } }{ ( 1 - \mathbb {X}^m ) } \right) . \end{aligned}$$

(5.5)

Then it holds for all $m \in {\mathbb {N}}$ that

$$\begin{aligned} \begin{aligned} \Theta _{ m }&= \Theta _{ m - 1 } - \left( \left[ \varepsilon + \tfrac{ \sqrt{ | \mathbb {M}^{ (1) }_m | } }{ \sqrt{ 1 - \mathbb {Y}^m } } \right] ^{ - 1 } \frac{ \gamma _m \mathbf{m}^{ (1) }_m }{ ( 1 - \mathbb {X}^m ) } , \dots , \left[ \varepsilon + \tfrac{ \sqrt{ | \mathbb {M}^{ ( \rho ) }_m | } }{ \sqrt{ 1 - \mathbb {Y}^m } } \right] ^{ - 1 } \frac{ \gamma _m \mathbf{m}^{ (\rho ) }_m }{ ( 1 - \mathbb {X}^m ) } \right) , \\ \mathbf{m}_m&= \mathbb {X} \, \mathbf{m}_{ m - 1 } + \frac{ ( 1 - \mathbb {X} ) }{ J_m } \left( \sum _{ j = 1 }^{ J_m } \Phi ^{ m - 1 , j }_{ \mathbb {S}_m }( \Theta _{ m - 1 } ) \right) , \\ \mathbb {M}_m&= \mathbb {Y} \, \mathbb {M}_{ m - 1 } + \left( 1 - \mathbb {Y} \right) {\text {Pow}}_{ 2 }\left( \frac{ 1 }{ J_m } \sum _{ j = 1 }^{ J_m } \Phi ^{ m - 1 , j }_{ \mathbb {S}_m }( \Theta _{ m - 1 } ) \right) . \end{aligned} \end{aligned}$$

(5.6)

1.3 Euler–Maruyama Scheme

Example 5.3

Assume the setting in Subsect. 3.2, let $\mu :[0,T] \times {\mathbb {R}}^d \rightarrow {\mathbb {R}}^d$ and $\sigma :[0,T] \times {\mathbb {R}}^d \rightarrow {\mathbb {R}}^d$ be functions, and assume for all $s, t \in [0,T]$, $x, w \in {\mathbb {R}}^d$ that

$$\begin{aligned} \Upsilon ( s, t, x, w ) = x + \mu ( s, x ) \, ( t - s ) + \sigma ( s, x ) \, w . \end{aligned}$$

(5.7)

Then it holds for all $m, j \in {\mathbb {N}}_0$, $n \in \{ 0, 1, \dots , N - 1 \}$ that

$$\begin{aligned} \mathcal {X}^{ m, j }_n = \mathcal {X}^{ m, j }_n + \mu \left( t_n, \mathcal {X}^{ m, j }_n \right) \, \left( t_{ n + 1 } - t_n \right) + \sigma \left( t_n, \mathcal {X}^{ m, j }_n \right) \, \left( W_{ t_{ n + 1 } } - W_{ t_n } \right) . \end{aligned}$$

(5.8)

In the setting of Example 5.3, we consider under suitable further hypotheses for every sufficiently large $m \in {\mathbb {N}}_0$ the random variable $\mathcal {U}^{ \Theta _m }$ as an approximation of $u(0,\xi )$ where $u :[0,T] \times {\mathbb {R}}^d \rightarrow {\mathbb {R}}^k$ is a suitable solution of the PDE

$$\begin{aligned}&\frac{ \partial u}{ \partial t } ( t, x ) + \frac{ 1 }{ 2 } \sum \limits _{ j = 1 }^d \left( \frac{ \partial ^2 u}{ \partial x^2 } \right) ( t, x )\left[ \sigma ( t, x ) \, e^{ (d) }_j , \sigma ( t, x ) \, e^{ (d) }_j \right] + \left( \frac{ \partial u}{ \partial x } \right) ( t, x ) \, \mu ( t, x ) \nonumber \\&\quad + f\left( t, x, u(t,x), \left( \frac{ \partial u}{ \partial x } \right) ( t, x ) \, \sigma ( t, x ) \right) = 0 \end{aligned}$$

(5.9)

with $u(T,x) = g(x)$, $e^{ (d) }_1 = (1,0,\dots ,0)$, $\dots$, $e^{ (d) }_d = (0,\dots ,0,1) \in {\mathbb {R}}^d$ for $t \in [0,T]$, $x = ( x_1, \dots , x_d ) \in {\mathbb {R}}^d$ (cf. (PDE) in Sect. 2 above).

1.4 Geometric Brownian Motion

Example 5.4

Assume the setting in Subsect. 3.2, let $\bar{\mu }, \bar{\sigma } \in {\mathbb {R}}$, and assume for all $s, t \in [0,T]$, $x = ( x_1, \dots , x_d )$, $w = ( w_1, \dots , w_d ) \in {\mathbb {R}}^d$ that

$$\begin{aligned} \Upsilon ( s, t, x, w ) \!=\! \exp \left( \left( \bar{\mu }\! -\! \frac{ \bar{\sigma }^2 }{ 2 } \right) ( t\! -\! s ) \right) \exp \left( \bar{\sigma } {\text {diag}}_{ {\mathbb {R}}^{ d \times d } }( w_1, \dots , w_d ) \right) x .\quad \end{aligned}$$

(5.10)

Then it holds for all $m, j \in {\mathbb {N}}_0$, $n \in \{ 0, 1, \dots , N \}$ that

$$\begin{aligned} \mathcal {X}^{ \theta , m, j }_n = \exp \left( \left( \bar{\mu } - \frac{ \bar{\sigma }^2 }{ 2 } \right) t_n {\text {Id}}_{ {\mathbb {R}}^d } + \bar{\sigma } {\text {diag}}_{ {\mathbb {R}}^{ d \times d } }\left( W_{ t_n }^{ m, j } \right) \right) \xi . \end{aligned}$$

(5.11)

In the setting of Example 5.4 we view under suitable further hypotheses (cf. Subsect. 4.4 above) for every sufficiently large $m \in {\mathbb {N}}_0$ the random variable $\mathcal {U}^{ \Theta _m }$ as an approximation of $u(0,\xi )$ where $u :[0,T] \times {\mathbb {R}}^d \rightarrow {\mathbb {R}}^k$ is a suitable solution of the PDE

$$\begin{aligned}&\tfrac{ \partial u}{ \partial t } ( t, x ) + \tfrac{\bar{\sigma }^2}{ 2 } \textstyle \sum \limits _{i=1}^d | x_i |^2 \, \big ( \tfrac{ \partial ^2 u}{ \partial x^2_i } \big )(t,x) + \bar{\mu } \sum \limits _{i=1}^d x_i \, \big (\tfrac{\partial u}{\partial x_i}\big )(t,x) \nonumber \\&\quad + f\big ( t, x, u(t,x), \bar{\sigma } \, ( \tfrac{ \partial u}{ \partial x } )( t, x ) {\text {diag}}_{ {\mathbb {R}}^{ d \times d } }(x_1, \dots , x_d) \big ) = 0 \end{aligned}$$

(5.12)

with $u(T,x) = g(x)$ for $t \in [0,T]$, $x = ( x_1, \dots , x_d ) \in {\mathbb {R}}^d$.