Solving high-dimensional partial differential equations using deep learning

Introduction

Partial differential equations (PDEs) are among the most ubiquitous tools used in modeling problems in nature. Some of the most important ones are naturally formulated as PDEs in high dimensions. Well-known examples include:

The Schrödinger equation in quantum many-body problem. In this case the dimensionality of the PDE is roughly three times the number of electrons or quantum particles in the system.

The nonlinear Black-Scholes equation for pricing financial derivatives, in which the dimensionality of the PDE is the number of underlying financial assets under consideration.

The Hamilton-Jacobi-Bellman equation in dynamic programming. In a game theory setting with multiple agents, the dimensionality goes up linearly with the number of agents. Similarly, in a resource allocation problem, the dimensionality goes up linearly with the number of devices and resources.

As elegant as these PDE models are, their practical use has proven to be very limited due to the curse of dimensionality : the computational cost for solving them goes up exponentially with the dimensionality.

Another area where the curse of dimensionality has been an essential obstacle is machine learning and data analysis, where the complexity of nonlinear regression models, for example, goes up exponentially with the dimensionality. In both cases the essential problem we face is how to represent or approximate a nonlinear function in high dimensions. The traditional approach, by building functions using polynomials, piecewise polynomials, wavelets, or other basis functions, is bound to run into the curse of dimensionality problem.

In recent years a new class of techniques, the deep neural network model, has shown remarkable success in artificial intelligence (see, e.g., ). Neural network is an old idea but recent experience has shown that deep networks with many layers seem to do a surprisingly good job in modeling complicated data sets. In terms of representing functions, the neural network model is compositional: it uses compositions of simple functions to approximate complicated ones. In contrast, the approach of classical approximation theory is usually additive. Mathematically, there are universal approximation theorems stating that a single hidden layer neural network can approximate a wide class of functions on compact subsets (see, e.g., survey and the references therein), even though we still lack a theoretical framework for explaining the seemingly unreasonable effectiveness of multilayer neural networks, which are widely employed nowadays. Despite this, the practical success of deep neural networks in artificial intelligence has been very astonishing and encourages applications to other problems where the curse of dimensionality has been a tormenting issue.

In this paper, we extend the power of deep neural networks to another dimension by developing a strategy for solving a large class of high-dimensional nonlinear PDEs using deep learning. The class of PDEs that we deal with are (nonlinear) parabolic PDEs. Special cases include the Black-Scholes equation and the Hamilton-Jacobi-Bellman equation. To do so, we make use of the reformulation of these PDEs as backward stochastic differential equations (BSDEs) (see, e.g., ) and approximate the gradient of the solution using deep neural networks. The methodology bears some resemblance to deep reinforcement learning with the BSDE playing the role of model-based reinforcement learning (or control theory models) and the gradient of the solution playing the role of policy function. Numerical examples manifest that the proposed algorithm is quite satisfactory in both accuracy and computational cost.

Due to the “curse of dimensionality”, there are only a very limited number of cases where practical high-dimensional algorithms have been developed in literature. For linear parabolic PDEs, one can use the Feynman-Kac formula and Monte Carlo methods to develop efficient algorithms to evaluate solutions at any given space-time locations. For a class of inviscid Hamilton-Jacobi equations, Darbon & Osher have recently developed an effective algorithm in the high-dimensional case (see ), based on the Hopf formula for the Hamilton-Jacobi equations. A general algorithm for nonlinear parabolic PDEs based on the multilevel decomposition of Picard iteration is developed in and has been shown to be quite efficient on a number of examples in finance and physics. The branching diffusion method has been proposed in , which exploits the fact that solutions of semilinear PDEs with polynomial nonlinearity can be represented as an expectation of a functional of branching diffusion processes. This method does not suffer from the curse of dimensionality, but still has limited applicability due to the blow up of approximated solutions in finite time.

The starting point of the present paper is deep learning. It should be stressed that even though deep learning has been a very successful tool for a number of applications, adapting it to the current setting with practical success is still a highly non-trivial task. Here by using the reformulation of BSDEs, we are able to cast the problem of solving PDEs as a learning problem and we design a deep learning framework that fits naturally to that setting. This has proven to be quite successful in practice.

Methodology

We consider a general class of PDEs known as semilinear parabolic PDEs. These PDEs can be represented as follows:

Let $\{W_{t}\}_{t\in[0,T]}$ be a $d$-dimensional Brownian motion and $\{X_{t}\}_{t\in[0,T]}$ be a $d$-dimensional stochastic process which satisfies

Then the solution of (1) satisfies the following BSDE (cf., e.g., ):

We refer to the Section Materials and Methods for further explanation of (3).

To derive a numerical algorithm to compute $u(0,X_{0})$, we treat $u(0,X_{0})\approx\theta_{u_{0}},\nabla u(0,X_{0})\approx\theta_{\nabla u_{0}}$ as parameters in the model and view (3) as a way of computing the values of $u$ at the terminal time $T$, knowing $u(0,X_{0})$ and $\nabla u(t,X_{t})$. We apply a temporal discretization to (2)–(3). Given a partition of the time interval $[0,T]$: $0=t_{0}<t_{1}<\ldots<t_{N}=T$, we consider the simple Euler scheme for $n=1,\dots,N-1$:

Given this temporal discretization, the path $\{X_{t_{n}}\}_{0\leq n\leq N}$ can be easily sampled using (4). Our key step next is to approximate the function $x\mapsto\sigma^{\operatorname{T}}(t,x)\,\nabla u(t,x)$ at each time step $t=t_{n}$ by a multilayer feedforward neural network

for $n=1,\dots,N-1$, where $\theta_{n}$ denotes parameters of the neural network approximating $x\mapsto\sigma^{\operatorname{T}}(t,x)\,\nabla u(t,x)$ at $t=t_{n}$.

Thereafter, we stack all the sub-networks in (7) together to form a deep neural network as a whole, based on the summation of (5) over $n=1,\dots,N-1$. Specifically, this network takes the paths $\{X_{t_{n}}\}_{0\leq n\leq N}$ and $\{W_{t_{n}}\}_{0\leq n\leq N}$ as the input data and gives the final output, denoted by $\hat{u}(\{{X_{t_{n}}}\}_{0\leq n\leq N},\{W_{t_{n}}\}_{0\leq n\leq N})$, as an approximation of $u(t_{N},X_{t_{N}})$. We refer to the Section Materials and Methods for more details on the architecture of the neural network. The difference in the matching of given terminal condition can be used to define the expected loss function

The total set of parameters are: $\theta=\{\theta_{u_{0}},\theta_{\nabla u_{0}},\theta_{1},\dots,\theta_{N-1}\}$.

We can now use a stochastic gradient descent-type (SGD) algorithm to optimize the parameter $\theta$, just as in the standard training of deep neural networks. In our numerical examples, we use the Adam optimizer . See the Section Materials and Methods for more details on the training of the deep neural networks. Since the BSDE is used as an essential tool, we call the methodology introduced above deep BSDE method.

Examples

A key issue in the trading of financial derivatives is to determine an appropriate fair price. Black & Scholes illustrated that the price $u$ of a financial derivative satisfies a parabolic PDE, nowadays known as the Black-Scholes equation . The Black-Scholes model can be augmented to take into account several important factors in real markets, including defaultable securities, higher interest rates for borrowing than for lending, transactions costs, uncertainties in the model parameters, etc. (see, e.g., ). Each of these effects results in a nonlinear contribution in the pricing model (see, e.g., ). In particular, the credit crisis and the ongoing European sovereign debt crisis have highlighted the most basic risk that has been neglected in the original Black-Scholes model, the default risk .

Ideally the pricing models should take into account the whole basket of underlyings that the financial derivatives depend on, resulting in high-dimensional nonlinear PDEs. However, existing pricing algorithms are unable to tackle these problems generally due to the curse of dimensionality. To demonstrate the effectiveness of the deep BSDE method, we study a special case of the recursive valuation model with default risk . We consider the fair price of a European claim based on 100 underlying assets conditional on no default having occurred yet. When default of the claim’s issuer occurs, the claim’s holder only receives a fraction $\delta\in[0,1)$ of the current value. The possible default is modeled by the first jump time of a Poisson process with intensity $Q$, a decreasing function of the current value, i.e., the default becomes more likely when the claim’s value is low. The value process can then be modeled by (1) with the generator

(see ), where $R$ is the interest rate of the risk-free asset. We assume that the underlying asset price moves as a geometric Brownian motion and choose the intensity function $Q$ as a piecewise-linear function of the current value with three regions ($v^{h}<v^{l},\,\gamma^{h}>\gamma^{l}$):

Hamilton-Jacobi-Bellman (HJB) Equation

The term “curse of dimensionality” was first used explicitly by Richard Bellman in the context of dynamic programming , which has now become the cornerstone in many areas such as economics, behavioral science, computer science, and even biology, where intelligent decision making is the main issue. In the context of game theory where there are multiple players, each player has to solve a high-dimensional HJB type equation in order to find his/her optimal strategy. In a dynamic resource allocation problem involving multiple entities with uncertainty, the dynamic programming principle also leads to a high-dimensional HJB equation for the value function. Until recently these high-dimensional PDEs have basically remained intractable. We now demonstrate below that the deep BSDE method is an effective tool for dealing with these high-dimensional problems.

We consider a classical linear-quadratic-Gaussian (LQG) control problem in 100 dimension:

(see e.g., Yong & Zhou [24, Chapter 4]). The value of the solution $u(t,x)$ of (13) at $t=0$ represents the optimal cost when the state starts from $x$. Applying Itô’s formula, one can show that the exact solution of (13) with the terminal condition $u(T,x)=g(x)$ admits the explicit formula

This can be used to test the accuracy of the proposed algorithm.

Allen-Cahn Equation

The Allen-Cahn equation is a reaction-diffusion equation that arises in physics, serving as a prototype for the modeling of phase separation and order-disorder transition (see, e.g., ). Here we consider a typical Allen-Cahn equation with the “double-well potential” in 100-dimensional space

Conclusions

The algorithm proposed in this paper opens up a host of new possibilities in several different areas. For example in economics one can consider many different interacting agents at the same time, instead of using the “representative agent” model. Similarly in finance, one can consider all the participating instruments at the same time, instead of relying on ad hoc assumptions about their relationships. In operational research, one can handle the cases with hundreds and thousands of participating entities directly, without the need to make ad hoc approximations.

It should be noted that although the methodology presented here is fairly general, we are so far not able to deal with the quantum many-body problem due to the difficulty in dealing with the Pauli exclusion principle.

Materials and Methods

(cf., e.g., ). Therefore, we can compute the quantity $u(0,X_{0})$ associated to (1) through $Y_{0}$ by solving the BSDE (16)–(17). More specifically, we plug the identities in (18) into (17) and rewrite the equation forwardly to obtain the formula in (3).

Then we discretize the equation temporally and use neural networks to approximate the spacial gradients and finally the unknown function, as introduced in the Section Methodology of the paper.

Neural Network Architecture

$X_{t_{n}}\to h_{n}^{1}\rightarrow h_{n}^{2}\rightarrow\cdots\rightarrow h_{n}^{H}\rightarrow\nabla u(t_{n},X_{t_{n}})$ is the multilayer feedforward neural network approximating the spatial gradients at time $t=t_{n}$. The weights $\theta_{n}$ of this sub-network are the parameters we aim to optimize.

$(u(t_{n},X_{t_{n}}),\nabla u(t_{n},X_{t_{n}}),W_{t_{n+1}}-W_{t_{n}})\rightarrow u(t_{n+1},X_{t_{n+1}})$ is the forward iteration giving the final output of the network as an approximation of $u(t_{N},X_{t_{N}})$, completely characterized by (5)–(6). There are no parameters to be optimized in this type of connection.

$(X_{t_{n}},W_{t_{n+1}}-W_{t_{n}})\rightarrow X_{t_{n+1}}$ is the shortcut connecting blocks at different time, which is characterized by (4) and (6). There are also no parameters to be optimized in this type of connection.

If we use $H$ hidden layers in each sub-network, as illustrated in Fig. 4, then the whole network has $(H+1)(N-1)$ layers in total that involve free parameters to be optimized simultaneously.

Implementation

We describe in detail the implementation for the numerical examples presented in the paper. Each sub-network is fully connected and consists of $4$ layers (except the example in the next subsection), with $1$ input layer ($d$-dimensional), $2$ hidden layers (both $d+10$-dimensional), and $1$ output layer ($d$-dimensional). We choose the rectifier function (ReLU) as our activation function. We also adopted the technique of batch normalization in the sub-networks, right after each linear transformation and before activation. This technique accelerates the training by allowing a larger step size and easier parameter initialization. All the parameters are initialized through a normal or a uniform distribution without any pre-training.

We use TensorFlow to implement our algorithm with the Adam optimizer to optimize parameters. Adam is a variant of the SGD algorithm, based on adaptive estimates of lower-order moments. We set the default values for corresponding hyper-parameters as recommended in and choose the batch size as 64. In each of the presented numerical examples the means and the standard deviations of the relative $L^{1}$-approximation errors are computed approximatively by means of 5 independent runs of the algorithm with different random seeds. All the numerical examples reported are run on a Macbook Pro with a 2.9GHz Intel Core i5 processor and 16 GB memory.

Effect of Number of Hidden Layers

The accuracy of the deep BSDE method certainly depends on the number of hidden layers in the sub-network approximation (7). To test this effect, we solve a reaction-diffusion type PDE with different number of hidden layers in the sub-network. The PDE is a high-dimensional version ($d=100$) of the example analyzed numerically in Gobet & Turkedjiev ($d=2$):

in which $u^{*}(t,x)$ is the explicit oscillating solution

Parameters are chosen in the same way as in : $\kappa=1.6,\,\lambda=0.1,\,T=1$. A residual structure with skip connection is used in each sub-network with each hidden layer having $d$ neurons. We increase the number of hidden layers in each sub-network from to $4$ and report the relative error in Table 1. It is evident that the approximation accuracy increases as the number of hidden layers in the sub-network increases.

Acknowledgement

The work of Han and E is supported in part by Major Program of NNSFC under grant 91130005, DOE grant DE-SC0009248 and ONR grant N00014-13-1-0338.

References