DGM: A deep learning algorithm for solving partial differential equations

Deep learning and high-dimensional PDEs

High-dimensional partial differential equations (PDEs) are used in physics, engineering, and finance. Their numerical solution has been a longstanding challenge. Finite difference methods become infeasible in higher dimensions due to the explosion in the number of grid points and the demand for reduced time step size. If there are $d$ space dimensions and $1$ time dimension, the mesh is of size $\mathcal{O}^{d+1}$. This quickly becomes computationally intractable when the dimension $d$ becomes even moderately large. We propose to solve high-dimensional PDEs using a meshfree deep learning algorithm. The method is similar in spirit to the Galerkin method, but with several key changes using ideas from machine learning. The Galerkin method is a widely-used computational method which seeks a reduced-form solution to a PDE as a linear combination of basis functions. The deep learning algorithm, or “Deep Galerkin Method” (DGM), uses a deep neural network instead of a linear combination of basis functions. The deep neural network is trained to satisfy the differential operator, initial condition, and boundary conditions using stochastic gradient descent at randomly sampled spatial points. By randomly sampling spatial points, we avoid the need to form a mesh (which is infeasible in higher dimensions) and instead convert the PDE problem into a machine learning problem.

DGM is a natural merger of Galerkin methods and machine learning. The algorithm in principle is straightforward; see Section 2. Promising numerical results are presented later in Section 4 for a class of high-dimensional free boundary PDEs. We also accurately solve a high-dimensional Hamilton-Jacobi-Bellman PDE in Section 5 and Burger’s equation in Section 6. DGM converts the computational cost of finite difference to a more convenient form: instead of a huge mesh of $\mathcal{O}^{d+1}$ (which is infeasible to handle), many batches of random spatial points are generated. Although the total number of spatial points could be vast, the algorithm can process the spatial points sequentially without harming the convergence rate.

Deep learning has revolutionized fields such as image, text, and speech recognition. These fields require statistical approaches which can model nonlinear functions of high-dimensional inputs. Deep learning, which uses multi-layer neural networks (i.e., “deep neural networks”), has proven very effective in practice for such tasks. A multi-layer neural network is essentially a “stack” of nonlinear operations where each operation is prescribed by certain parameters that must be estimated from data. Performance in practice can strongly depend upon the specific form of the neural network architecture and the training algorithms which are used. The design of neural network architectures and training methods has been the focus of intense research over the past decade. Given the success of deep learning, there is also growing interest in applying it to a range of other areas in science and engineering (see Section 1.2 for some examples).

Evaluating the accuracy of the deep learning algorithm is not straightforward. PDEs with semi-analytic solutions may not be sufficiently challenging. (After all, the semi-analytic solution exists since the PDE can be transformed into a lower-dimensional equation.) It cannot be benchmarked against traditional finite difference (which fails in high dimensions). We test the deep learning algorithm on a class of high-dimensional free boundary PDEs which have the special property that error bounds can be calculated for any approximate solution. This provides a unique opportunity to evaluate the accuracy of the deep learning algorithm on a class of high-dimensional PDEs with no semi-analytic solutions.

This class of high-dimensional free boundary PDEs also has important applications in finance, where it used to price American options. An American option is a financial derivative on a portfolio of stocks. The number of space dimensions in the PDE equals the number of stocks in the portfolio. Financial institutions are interested in pricing options on portfolios ranging from dozens to even hundreds of stocks . Therefore, there is a significant need for numerical methods to accurately solve high-dimensional free boundary PDEs.

We also test the deep learning algorithm on a high-dimensional Hamilton-Jacobi-Bellman PDE with accurate results. We consider a high-dimensional Hamilton-Jacobi-Bellman PDE motivated by the problem of optimally controlling a stochastic heat equation.

Finally, it is often of interest to find the solution of a PDE over a range of problem setups (e.g., different physical conditions and boundary conditions). For example, this may be useful for the design of engineering systems or uncertainty quantification. The problem setup space may be high-dimensional and therefore may require solving many PDEs for many different problem setups, which can be computationally expensive. We use our deep learning algorithm to approximate the general solution to the Burgers’ equation for different boundary conditions, initial conditions, and physical conditions.

In the remainder of the Introduction, we provide an overview of our results regarding the approximation power of neural networks for quasilinear parabolic PDEs (Section 1.1), and relevant literature (Section 1.2). The deep learning algorithm for solving PDEs is presented in Section 2. An efficient scheme for evaluating the diffusion operator is developed in Section 3. Numerical analysis of the algorithm is presented in Sections 4, 5, and 6. We implement and test the algorithm on a class of high-dimensional free boundary PDEs in up to 200 dimensions. The theorem and proof for the approximation of PDE solutions with neural networks is presented in Section 7. Conclusions are in Section 8. For readability purposes, proofs from Section 7 have been collected in Appendix A.

We also prove a theorem regarding the approximation power of neural networks for a class of quasilinear parabolic PDEs. Consider the potentially nonlinear PDE

where $\partial\Omega$ is the boundary of the domain $\Omega$. The solution $u(t,x)$ is of course unknown, but an approximate solution $f(t,x)$ can be found by minimizing the L2 error

strongly in, $L^{\rho}([0,T]\times\Omega)$, with $\rho<2$, for a class of quasilinear parabolic PDEs; see subsection 7.2 and Theorem 7.3 therein for the precise statement. That is, the neural network will converge in $L^{\rho},\rho<2$ to the solution of the PDE as the number of hidden units tends to infinity. The precise statement of the theorem and its proof are presented in Section 7. The proof requires the joint analysis of the approximation power of neural networks as well as the continuity properties of partial differential equations. Note that $J(f^{n})\rightarrow 0$ does not necessarily imply that $f^{n}\rightarrow u$, given that we only have $L^{2}$ control on the approximation error. First, we prove that $J(f^{n})\rightarrow 0$ as $n\rightarrow\infty$. We then establish that each neural network $\{f^{n}\}_{n=1}^{\infty}$ satisfies a PDE with a source term $h^{n}(t,x)$. We are then able to prove, under certain conditions, the convergence of $f^{n}\rightarrow u$ as $n\rightarrow\infty$ in $L^{\rho}([0,T]\times\Omega)$, for $\rho<2$, using the smoothness of the neural network approximations and compactness arguments.

Theorem 7.3 establishes the approximation power of neural networks for solving PDEs (at least within a class of quasilinear parabolic PDEs); however, directly minimizing $J(f)$ is not computationally tractable since it involves high-dimensional integrals. The DGM algorithm minimizes $J(f)$ using a meshfree approach; see Section 2.

2 Relevant Literature

Solving PDEs with a neural network as an approximation is a natural idea, and has been considered in various forms previously. , , , , and propose to use neural networks to solve PDEs and ODEs. These papers estimate neural network solutions on an a priori fixed mesh. This paper proposes using deep neural networks and is meshfree, which is key to solving high-dimensional PDEs.

In particular, this paper explores several new innovations. First, we focus on high-dimensional PDEs and apply deep learning advances of the past decade to this problem (deep neural networks instead of shallow neural networks, improved optimization methods for neural networks, etc.). Algorithms for high-dimensional free boundary PDEs are developed, efficiently implemented, and tested. In particular, we develop an iterative method to address the free boundary. Secondly, to avoid ever forming a mesh, we sample a sequence of random spatial points. This produces a meshfree method, which is essential for high-dimensional PDEs. Thirdly, the algorithm incorporates a new computational scheme for the efficient computation of neural network gradients arising from the second derivatives of high-dimensional PDEs.

Recently, develop physics informed deep learning models. They estimate deep neural network models which merge data observations with PDE models. This allows for the estimation of physical models from limited data by leveraging a priori knowledge that the physical dynamics should obey a class of PDEs. Their approach solves PDEs in one and two spatial dimensions using deep neural networks. uses a deep neural network to model the Reynolds stresses in a Reynolds-averaged Navier-Stokes (RANS) model. RANS is a reduced-order model for turbulence in fluid dynamics. have also recently developed a scheme for solving a class of quasilinear PDEs which can be represented as forward-backward stochastic differential equations (FBSDEs) and further develops the algorithm. The algorithm developed in focuses on computing the value of the PDE solution at a single point. The algorithm that we present here is different; in particular, it does not rely on the availability of FBSDE representations and yields the entire solution of the PDE across all time and space. In addition, the deep neural network architecture that we use, which is different from the ones used in , seems to be able to recover accurately the entire solution (at least for the equations that we studied). use a convolutional neural network to solve a large sparse linear system which is required in the numerical solution of the Navier-Stokes PDE. In addition, has recently developed a novel partial differential equation approach to optimize deep neural networks.

developed an algorithm for the solution of a discrete-time version of a class of free boundary PDEs. Their algorithm, commonly called the “Longstaff-Schwartz method”, uses dynamic programming and approximates the solution using a separate function approximator at each discrete time (typically a linear combination of basis functions). Our algorithm directly solves the PDE, and uses a single function approximator for all space and all time. The Longstaff-Schwartz algorithm has been further analyzed by , , and others. Sparse grid methods have also been used to solve high-dimensional PDEs; see , , , , and .

In regards to general results on the approximation power of neural networks we refer the interested reader to classical works and we also mention the recent work by , where the authors study the necessary and sufficient complexity of ReLU neural networks that is required for approximating classifier functions in the mean square sense.

Algorithm

Consider a parabolic PDE with $d$ spatial dimensions:

Here, $\left\lVert f(y)\right\rVert^{2}_{\mathcal{Y},\nu}=\int_{\mathcal{Y}}\left|f(y)\right|^{2}\nu(y)dy$ where $\nu(y)$ is a positive probability density on $y\in\mathcal{Y}$. $J(f)$ measures how well the function $f(t,x;\theta)$ satisfies the PDE differential operator, boundary conditions, and initial condition. If $J(f)=0$, then $f(t,x;\theta)$ is a solution to the PDE (2.1).

The goal is to find a set of parameters $\theta$ such that the function $f(t,x;\theta)$ minimizes the error $J(f)$. If the error $J(f)$ is small, then $f(t,x;\theta)$ will closely satisfy the PDE differential operator, boundary conditions, and initial condition. Therefore, a $\theta$ which minimizes $J(f(\cdot;\theta))$ produces a reduced-form model $f(t,x;\theta)$ which approximates the PDE solution $u(t,x)$.

Estimating $\theta$ by directly minimizing $J(f)$ is infeasible when the dimension $d$ is large since the integral over $\Omega$ is computationally intractable. However, borrowing a machine learning approach, one can instead minimize $J(f)$ using stochastic gradient descent on a sequence of time and space points drawn at random from $\Omega$ and $\partial\Omega$. This avoids ever forming a mesh.

Generate random points $(t_{n},x_{n})$ from $[0,T]\times\Omega$ and $(\tau_{n},z_{n})$ from $[0,T]\times\partial\Omega$ according to respective probability densities $\nu_{1}$ and $\nu_{2}$. Also, draw the random point $w_{n}$ from $\Omega$ with probability density $\nu_{3}$.

Calculate the squared error $G(\theta_{n},s_{n})$ at the randomly sampled points $s_{n}=\{(t_{n},x_{n}),(\tau_{n},z_{n}),w_{n}\}$ where:

Take a descent step at the random point $s_{n}$:

Repeat until convergence criterion is satisfied.

The “learning rate” $\alpha_{n}$ decreases with $n$. The steps $\nabla_{\theta}G(\theta_{n},s_{n})$ are unbiased estimates of $\nabla_{\theta}J(f(\cdot;\theta_{n}))$:

Therefore, the stochastic gradient descent algorithm will on average take steps in a descent direction for the objective function $J$. A descent direction means that the objective function decreases after an iteration (i.e., $J(f(\cdot;\theta_{n+1}))<J(f(\cdot;\theta_{n}))$ ), and $\theta_{n+1}$ is therefore a better parameter estimate than $\theta_{n}$.

Under (relatively mild) technical conditions (see ), the algorithm $\theta_{n}$ will converge to a critical point of the objective function $J(f(\cdot;\theta))$ as $n\rightarrow\infty$:

It’s important to note that $\theta_{n}$ may only converge to a local minimum when $f(t,x;\theta)$ is non-convex. This is generally true for non-convex optimization and is not specific to this paper’s algorithm. In particular, deep neural networks are non-convex. Therefore, it is well known that stochastic gradient descent may only converge to a local minimum (and not a global minimum) for a neural network. Nevertheless, stochastic gradient descent has proven very effective in practice and is the fundamental building block of nearly all approaches for training deep learning models.

A Monte Carlo Method for Fast Computation of Second Derivatives

This section describes a modified algorithm which may be more computationally efficient in some cases. The term $\mathcal{L}f(t,x;\theta)$ contains second derivatives $\frac{\partial^{2}f}{\partial x_{i}x_{j}}(t,x;\theta)$ which may be expensive to compute in higher dimensions. For instance, $20,000$ second derivatives must be calculated in $d=200$ dimensions.

The complicated architectures of neural networks can make it computationally costly to calculate the second derivatives (for example, see the neural network architecture (4.2)). The computational cost for calculating second derivatives (in both total arithmetic operations and memory) is $\mathcal{O}(d^{2}\times N)$ where $d$ is the spatial dimension of $x$ and $N$ is the batch size. In comparison, the computational cost for calculating first derivatives is $\mathcal{O}(d\times N)$. The cost associated with the second derivatives is further increased since we actually need the third-order derivatives $\nabla_{\theta}\frac{\partial^{2}f}{\partial x^{2}}(t,x;\theta)$ for the stochastic gradient descent algorithm. Instead of directly calculating these second derivatives, we approximate the second derivatives using a Monte Carlo method.

Suppose the sum of the second derivatives in $\mathcal{L}f(t,x,;\theta)$ is of the form $\frac{1}{2}\sum_{i,j=1}^{d}\rho_{i,j}\sigma_{i}(x)\sigma_{j}(x)\frac{\partial^{2}f}{\partial x_{i}x_{j}}(t,x;\theta)$, assume $[\rho_{i,j}]_{i,j=1}^{d}$ is a positive definite matrix, and define \sigma(x)=\bigg{(}\sigma_{1}(x),\ldots,\sigma_{d}(x)\bigg{)}. For example, such PDEs arise when considering expectations of functions of stochastic differential equations, where the $\sigma(x)$ represents the diffusion coefficient. See equation (4.1) and the corresponding discussion. A generalization of the algorithm in this section to second derivatives with nonlinear coefficient dependence on $u(t,x)$ is also possible. Then,

The DGM algorithm use the gradient $\nabla_{\theta}G_{1}(\theta_{n},s_{n})$, which requires the calculation of the second derivative terms in $\mathcal{L}f(t_{n},x_{n};\theta_{n})$. Define the first derivative operators as

Generate random points $(t_{n},x_{n})$ from $[0,T]\times\Omega$ and $(\tau_{n},z_{n})$ from $[0,T]\times\partial\Omega$ according to respective densities $\nu_{1}$ and $\nu_{2}$. Also, draw the random point $w_{n}$ from $\Omega$ with density $\nu_{3}$.

Repeat until convergence criterion is satisfied.

In conclusion, the modified algorithm here is computationally less expensive than the original algorithm in Section 2 but introduces some bias and variance. The variance essentially increases the i.i.d. noise in the stochastic gradient descent step; this noise averages out over a large number of samples though. The original algorithm in Section 2 is unbiased and has lower variance, but is computationally more expensive. We numerically implement the algorithm for a class of free boundary PDEs in Section 4. Future research may investigate other methods to further improve the computational evaluation of the second derivative terms (for instance, multi-level Monte Carlo).

Numerical Analysis for a High-dimensional Free Boundary PDE

Besides the high dimensions and the free boundary, the American option PDE is challenging to numerically solve since the payoff function $g(x)$ (which both appears in the initial condition and determines the free boundary) is not continuously differentiable.

Section 4.1 states the free boundary PDE and the deep learning algorithm to solve it. To address the free boundary, we supplement the algorithm presented in Section 2 with an iterative method; see Section 4.1. Section 4.2 describes the architecture and implementation details for the neural network. Section 4.3 reports numerical accuracy for a case where a semi-analytic solution exists. Section 4.4 reports numerical accuracy for a case where no semi-analytic solution exists.

We now specify the free boundary PDE for $u(t,x)$. The stock price dynamics and option price are:

The model (4.1) for the stock price dynamics is widely used in practice and captures several desirable characteristics. First, the drift $\mu(x)$ measures the “average” growth in the stock prices. The Brownian motion $W_{t}$ represents the randomness in the stock price, and the magnitude of the randomness is given by the coefficient function $\sigma(X_{t}^{i})$. The movement of stock prices are correlated (e.g., if Microsoft’s price increases, it is likely that Apple’s price will also increase). The magnitude of the correlation between two stocks $i$ and $j$ is specified by the parameter $\rho_{i,j}$. An example is the well-known Black-Scholes model $\mu(x)=\mu x$ and $\sigma(x)=\sigma x$. In the Black-Scholes model, the average rate of return for each stock is $\mu$.

The price function $u(t,x)$ in (4.1) is the solution to a free boundary PDE and will satisfy:

The free boundary set is F=\big{\{}(t,x):u(t,x)=g(x)\big{\}}. $u(t,x)$ satisfies a partial differential equation “above” the free boundary set $F$, and $u(t,x)$ equals the function $g(x)$ “below” the free boundary set $F$.

The deep learning algorithm for solving the PDE (4.1) requires simulating points above and below the free boundary set $F$. We use an iterative method to address the free boundary. The free boundary set $F$ is approximated using the current parameter estimate $\theta_{n}$. This approximate free boundary is used in the probability measure that we simulate points with. The gradient is not taken with respect to the $\theta_{n}$ input of the probability density used to simulate random points. For this purpose, define the objective function:

Generate the random batch of points $B^{3}=\{w_{m}\}_{m=1}^{M}$ from $\Omega$ with probability density $\nu_{3}$.

Take a descent step for the random batch $S_{n}$:

Repeat until convergence criterion is satisfied.

The second derivatives in the above algorithm can be approximated using the method from Section 3.

2 Implementation details for the algorithm

This section provides details for the implementation of the algorithm, including the DGM network architecture, hyperparameters, and computational approach.

The architecture of a neural network can be crucial to its success. Frequently, different applications require different architectures. For example, convolution networks are essential for image recognition while long short-term networks (LSTMs) are useful for modeling sequential data. Clever choices of architectures, which exploit a priori knowledge about an application, can significantly improve performance. In the PDE applications in this paper, we found that a neural network architecture similar in spirit to that of LSTM networks improved performance.

The PDE solution requires a model $f(t,x;\theta)$ which can make “sharp turns” due to the final condition, which is of the form $u(T,x)=\max(p(x),0)$ (the first derivative is discontinuous when $p(x)=0$). The shape of the solution $u(t,x)$ for $t<T$, although “smoothed” by the diffusion term in the PDE, will still have a nonlinear profile which is rapidly changing in certain spatial regions. In particular, we found the following network architecture to be effective:

where $\overset{\rightarrow}{x}=(t,x)$, the number of hidden layers is $L+1$, and $\odot$ denotes element-wise multiplication (i.e., z\odot v=\big{(}z_{0}v_{0},\ldots,z_{N}v_{N}\big{)}). The parameters are

The architecture (4.2) is relatively complicated. Within each layer, there are actually many “sub-layers” of computations. The important feature is the repeated element-wise multiplication of nonlinear functions of the input. This helps to model more complicated functions which are rapidly changing in certain time and space regions. The neural network architecture (4.2) is similar to the architecture for LSTM networks (see ) and highway networks (see ).

We emphasize that the only input to the network is $(t,x)$. We do not use any custom-designed nonlinear transformations of $(t,x)$. If properly chosen, such additional inputs might help performance. For example, the European option PDE solution (which has an analytic formula) could be included as an input.

Our computational approach to training the neural network involved several components. The second derivatives are approximated using the method from Section 3. Training is distributed across $6$ GPU nodes using asynchronous stochastic gradient descent (we provide more details on this below). Parameters are updated using the well-known ADAM algorithm (see ) with a decaying learning rate schedule (more details on the learning rate are provided below). Accuracy can be improved by calculating a running average of the neural network solutions over a sequence of training iterations (essentially a computationally cheap approach for building a model ensemble). We also found that model ensembles (of even small sizes of 5) can slightly increase accuracy.

Training of the neural network is distributed across several GPU nodes in order to accelerate training. We use asynchronous stochastic gradient descent, which is a widely-used method for parallelizing training of machine learning models. On each node, i.i.d. space and time samples are generated. Each node calculates the gradient of the objective function with respect to the parameters on its respective batch of simulated data. These gradients are then used to update the model, which is stored on a central node called a “parameter server”. Figure 1 displays the computational setup. Updates occur asynchronously; that is, node $i$ updates the model immediately upon completion of its work, and does not wait for node $j$ to finish its work. The “work” here is calculating the gradients for a batch of simulated data. Before a node calculates the gradient for a new batch of simulated data, it receives an updated model from the parameter server. For more details on asynchronous stochastic gradient descent, see .

During training, we decrease the learning as the number of iterations increases. We use a learning rate schedule where the learning rate is a piecewise constant function of the number of iterations. This is a typical choice. We found the following learning rate schedule to be effective:

We use approximately $100,000$ iterations. An “iteration” involves batches of size $1,000$ on each of the GPU nodes. Therefore, there are $5,000$ simulated time/space points for each iteration. In total, we used approximately $500$ million simulated time/space points to train the neural network.

We implement the algorithm using TensorFlow and PyTorch, which are software libraries for deep learning. TensorFlow has reverse mode automatic differentiation which allows the calculation of derivatives for a broad range of functions. For example, TensorFlow can be used to calculate the gradient of the neural network (4.2) with respect to $x$ or $\theta$. TensorFlow also allows for the training of models on graphics processing units (GPUs). A GPU, which has thousands of cores, can be use to highly parallelize the training of deep learning models. We furthermore distribute our computations across multiple GPU nodes, as described above. The computations in this paper were performed on the Blue Waters supercomputer which has a large number of GPU nodes.

3 A High-dimensional Free Boundary PDE with a Semi-Analytic Solution

We implement our deep learning algorithm to solve the PDE (4.1). The accuracy of our deep learning algorithm is evaluated in up to $200$ dimensions. The results are reported below in Table 1.

The semi-analytic solution used in Table 1 is provided below. Let $\mu(x)=(r-c)x$, $\sigma(x)=\sigma x$, and $\rho_{i,j}=\rho$ for $i\neq j$ (i.e., the Black-Scholes model). If the payoff function in (4.1) is g(x)=\max\big{(}(\prod_{i=1}^{d}x_{i})^{1/d}-K,0\big{)}, then there is a semi-analytic solution to (4.1):

where $v(t,x)$ satisfies the one-dimensional free boundary PDE

where $\hat{\sigma}^{2}=\frac{d\sigma^{2}+d(d-1)\rho\sigma^{2}}{d^{2}}$, $\hat{\mu}=(r-c)-\frac{1}{2}\hat{\sigma}^{2}+\frac{1}{2}\sigma^{2}$, and $\hat{g}(x)=\max(x,0)$. The one-dimensional PDE (4.5) can be solved using finite difference methods. If $f(t,x;\theta)$ is the deep learning algorithm’s estimate for the PDE solution at $(t,x)$, the relative error at the point $(t,x)$ is $\frac{|f(t,x;\theta)-u(t,x)|}{|u(t,x)|}\times 100\%$ and the absolute error at the point $(t,x)$ is $|f(t,x;\theta)-u(t,x)|$. The relative error and absolute error at the point $(t,x)$ can be evaluated using the semi-analytic solution (4.4).

Although the solution at $(t,x)=(0,X_{0})$ is of primary interest for American options, most other PDE applications are interested in the entire solution $u(t,x)$. The deep learning algorithm provides an approximate solution across all time and space $(t,x)\in[0,T]\times\Omega$. As an example, we present in Figure 2 contour plots of the absolute error and percent error across time and space for the American option PDE in $20$ dimensions. The contour plot is produced in the following way:

Figure 2 reports both the absolute error and the percent error. The percent error $\frac{|f(t,x;\theta)-u(t,x)|}{|u(t,x)|}\times 100\%$ is reported for points where $|u(t,x)|>0.05$. The absolute error becomes relatively large in a few areas; however, the solution $u(t,x)$ also grows large in these areas and therefore the percent error remains small.

4 A High-dimensional Free Boundary PDE without a Semi-Analytic Solution

We now consider a case of the American option PDE which does not have a semi-analytic solution. The American option PDE has the special property that it is possible to calculate error bounds on an approximate solution. Therefore, we can evaluate the accuracy of the deep learning algorithm even on cases where no semi-analytic solution is available.

We previously only considered a symmetrical case where $\rho_{i,j}=0.75$ and $\sigma=0.25$ for all stocks. This section solves a more challenging heterogeneous case where $\rho_{i,j}$ and $\sigma_{i}$ vary across all dimensions $i=1,2,\ldots,d$. The coefficients are fitted to actual data for the stocks IBM, Amazon, Tiffany, Amgen, Bank of America, General Mills, Cisco, Coca-Cola, Comcast, Deere, General Electric, Home Depot, Johnson & Johnson, Morgan Stanley, Microsoft, Nordstrom, Pfizer, Qualcomm, Starbucks, and Tyson Foods from 2000-2017. This produces a PDE with widely-varying coefficients for each of the $\frac{d^{2}+d}{2}$ second derivative terms. The correlation coefficients $\rho_{i,j}$ range from $-0.53$ to $0.80$ for $i\neq j$ and $\sigma_{i}$ ranges from $0.09$ to $0.69$.

Let $f(t,x;\theta)$ be the neural network approximation. derived that the PDE solution $u(t,x)$ lies in the interval:

where $\tau=\inf\{t\in[0,T]:f(t,X_{t};\theta)<g(X_{t})\}$ and $M_{s}$ is a martingale constructed from the approximate solution $f(t,x;\theta)$

The bounds (4.6) depend only on the approximation $f(t,x;\theta)$, which is known, and can be evaluated via Monte Carlo simulation. The integral for $M_{s}$ must also be discretized. The best estimate for the price of the American option is the midpoint of the interval $[\underline{u}(0,X_{0}),\overline{u}(0,X_{0})]$, which has an error bound of $\frac{\overline{u}(0,X_{0})-\underline{u}(0,X_{0})}{2\underline{u}(0,X_{0})}\times 100\%$. Numerical results are in Table 2.

We present in Figure 3 contour plots of the absolute error bound and percent error bound across time and space for the American option PDE in $20$ dimensions for strike price $K=1$. The contour plot is produced in the following way:

Figure 3 reports both the absolute error and the percent error. The percent error $\frac{|f(t,x;\theta)-u(t,x)|}{|u(t,x)|}\times 100\%$ is reported for points where $|u(t,x)|>0.05$. It should be emphasized that these are error bounds; therefore, the actual error could be lower. The contour plot 3 requires significant computations. For each point at which calculate an error bound, a new simulation of (4.6) is required. In total, a large number of simulations are required, which we distribute across hundreds of GPUs on the Blue Waters supercomputer.

High-dimensional Hamilton-Jacobi-Bellman PDE

We also test the deep learning algorithm on a high-dimensional Hamilton-Jacobi-Bellman (HJB) equation corresponding to the optimal control of a stochastic heat equation. Specifically, we demonstrate that the deep learning algorithm accurately solves the high-dimensional PDE (5.5). The PDE (5.5) is motivated by the problem of optimally controlling the stochastic partial differential equation (SPDE):

The constant $\gamma>0$ is a discount factor. The constant $\lambda>0$ penalizes large values for the control $u(x)$. The goal is to reach the target $\bar{v}(x)$ while expending the minimum amount of energy. The optimal control $u(x)$ satisfies an infinite-dimensional HJB equation. We refer the reader to Theorems 5.3 and 5.4 of as well as and for an analysis of infinite-dimensional HJB equations for the stochastic heat equation.

An example of a problem represented by the SPDE (5.1) is the heating of a rod to a target temperature profile. One can control the heat applied to each portion of the rod along its length. There are also random fluctuations in the temperature of the rod due to other environmental factors, which is represented by the Brownian sheet $W(t,x)$. The goal is to guide the temperature profile of the rod to the target profile while expending the least amount of energy; see the objective function (5.2).

(5.1) can be discretized in space, which yields a system of stochastic differential equations (SDEs). (For example, see Section 3.2 of .) This system of SDEs can be used to derive a finite, high-dimensional PDE for the value function and optimal control. That is, we first approximate the SPDE with a finite-dimensional system of SDEs, and then we solve the high-dimensional PDE corresponding to the finite-dimensional system of SDEs.

where $\Delta$ is the mesh size, $v(t,j\Delta)=X_{t}^{j}$, $u(j\Delta)=U^{j}_{t}$, and $W_{t}^{j}$ are independent standard Brownian motions (see , , and regarding numerical schemes for stochastic parabolic PDEs of the form considered in this section). The dimension of the SDE system (5.3) is $d=\frac{L}{\Delta}-1$. Note that (5.3) uses a central difference scheme for the diffusion term in (5.1).

The value function $V(x)$ satisfies a nonlinear PDE with $d$ spatial dimensions $x_{1},x_{2},\ldots,x_{d}$.

The vector $\bar{v}=(\bar{v}(\Delta),\bar{v}(2\Delta),\ldots,\bar{v}(d\Delta))$. Note that the values $x_{d+1}=\bar{v}(L)$ and $x_{0}=\bar{v}(0)$ are constants which correspond to the boundary conditions in (5.1). The PDE (5.5) is high dimensional since the number of dimensions $d=\frac{L}{\Delta}-1$. The optimal control is

We solve the PDE (5.5) using the deep learning algorithm for $d=21$ dimensions. The size of the domain is $L=10^{-1}$. The coefficients are $\alpha=10^{-4}$, $\sigma=10^{-\frac{1}{2}}$, $\lambda=1$, and $\gamma=1$. The target profile is $\bar{v}(x)=0$.

The deep learning algorithm’s accuracy can be evaluated since a semi-analytic solution is available for (5.5).The PDE (5.5) has a semi-analytic solution which satisfies a Riccati equation. The Riccati equation can be solved using an iterative method. Figure 4 shows a contour plot of the percent error over space. The contour plot is produced in the following way:

The average percent error over the entire space is $0.1\%$.

Lastly, we close this section by mentioning that in the recent paper (see also ) the authors develop a machine learning algorithm that provides the value at a single point in time and space of the solution to a class of HJB equations which admit explicit solution that can be obtained through the Cole-Hopf transformation. Their method relies on characterizing the solution via backward stochastic differential equations (BSDE). In contrast, the current work (a) does not rely on BSDE type representations through nonlinear Feynman-Kac formulas, and (b) allows to recover the whole object (i.e. the solution across all points in time and space).

Burgers’ equation

It is often of interest to find the solution of a PDE over a range of problem setups (e.g., different physical conditions and boundary conditions). For example, this may be useful for the design of engineering systems or uncertainty quantification. The problem setup space may be high-dimensional and therefore may require solving many PDEs for many different problem setups, which can be computationally expensive.

Let the variable $p$ represent the problem setup (i.e., physical conditions, boundary conditions, and initial conditions). The variable $p$ takes values in the space $\mathcal{P}$, and we are interested in the solution of the PDE $u(t,x;p)$. (This is sometimes called a “parameterized class of PDEs”.) In particular, suppose $u(t,x;p)$ satisfies the PDE

A traditional approach would be to discretize the $\mathcal{P}$-space and re-solve the PDE many times for many different points $p$. However, the total number of grid points (and therefore the number of PDEs that must be solved) grows exponentially with the number of dimensions, and $\mathcal{P}$ is typically high-dimensional.

We propose to use the DGM algorithm to approximate the general solution to the PDE (6.1) for different boundary conditions, initial conditions, and physical conditions. The deep neural network is trained using stochastic gradient descent on a sequence of random time, space, and problem setup points $(t,x,p)$. Similar to before,

Update $\theta$ with a stochastic gradient descent step

If $x$ is low-dimensional ($d\leq 3$), which is common in many physical PDEs, the first and second partial derivatives of $f$ can be calculated via chain rule or approximated by finite difference. We implement our algorithm for Burgers’ equation on a finite domain.

Figure 6 presents the accuracy of the deep learning algorithm for different times $t$ and different choices of $\nu$. As $\nu$ becomes smaller, the solution becomes steeper. It also shows the shock layer forming over time. The contour plot (7) reports the absolute error of the deep learning solution for different choices of $b$ and $\nu$.

Neural Network Approximation Theorem for PDEs

Let the L2 error $J(f)$ measure how well the neural network $f$ satisfies the differential operator, boundary condition, and initial condition. Define $\mathfrak{C}^{n}$ as the class of neural networks with $n$ hidden units and let $f^{n}$ be a neural network with $n$ hidden units which minimizes $J(f)$. We prove that

in the appropriate sense, for a class of quasilinear parabolic PDEs with the principle term in divergence form under certain growth and smoothness assumptions on the nonlinear terms. Our theoretical result only covers a class of quasilinear parabolic PDEs as described in this section. However, the numerical results of this paper indicate that the results are more broadly applicable.

The proof requires the joint analysis of the approximation power of neural networks as well as the continuity properties of partial differential equations. First, we show that the neural network can satisfy the differential operator, boundary condition, and initial condition arbitrarily well for sufficiently large $n$.

Let $u$ be the solution to the PDE. The statement (7.1) does not necessarily imply that $f^{n}\rightarrow u$. One challenge to proving convergence is that we only have $L^{2}$ control of the error. We prove convergence for the case of homogeneous boundary data, i.e., $g(t,x)=0$, by first establishing that each neural network $\{f^{n}\}_{n=1}^{\infty}$ satisfies a PDE with a source term $h^{n}(t,x)$. Importantly, the source terms $h^{n}(t,x)$ are only known to be vanishing in $L^{2}$. We are then able to prove that the convergence of $f^{n}\rightarrow u$ as $n\rightarrow\infty$ in the appropriate space holds using compactness arguments.

The precise statement of the theorem and the presentation of the proof is in the next two sections. Section 7.1 proves that $J(f^{n})\rightarrow 0$ as $n\rightarrow\infty$. Section 7.2 contains convergence results of $f^{n}$ to the solution $u$ of the PDE as $n\rightarrow\infty$. The main result is Theorem 7.3. For readability purposes the corresponding proofs are in Appendix A.

In this section, we present a theorem guaranteeing the existence of multilayer feed forward networks $f$ able to universally approximate solutions of quasilinear parabolic PDEs in the sense that there is $f$ that makes the objective function $J(f)$ arbitrarily small. To do so, we use the results of on universal approximation of functions and their derivatives and make appropriate assumptions on the coefficients of the PDEs to guarantee that a classical solution exists (since then the results of apply).

For notational convenience, let us write the operator of (7.2) as $\mathcal{G}$. Namely, let us denote

For the purposes of this section, we consider equations of the type (7.2) that have classical solutions.

In particular we assume that there is a unique $u(t,x)$ solving (7.2) such that

We refer the interested reader to Theorems 5.4, 6.1 and 6.2 of Chapter V in for specific general conditions on $\alpha,\gamma$ guaranteeing the validity of the aforementioned statement.

Universal approximation results for single functions and their derivatives have been obtained under various assumptions in . In this paper, we use Theorem 3 of . Let us recall the setup appropriately modified for our case of interest. Let $\psi$ be an activation function, e.g., of sigmoid type, of the hidden units and define the set

The proof of this theorem is in the Appendix.

2 Convergence of the neural network to the PDE solution

We now prove, under stronger conditions, the convergence of the neural networks $f^{n}$ to the solution $u$ of the PDE

as $n\rightarrow\infty$. Notice that we have restricted the discussion to homogeneous boundary data. We do this for both presentation and mathematical reasons. We set $u(t,x)=0,\textrm{ for }(t,x)\in\partial\Omega_{T}$, i.e., $g=0$, to circumvent certain technical difficulties arising due to inhomogeneous boundary conditions. If $g\neq 0$ such that $g$ is the trace of some appropriately smooth function, say $\phi$, then one can reduce the inhomogeneous boundary conditions on $\partial\Omega_{T}$ to the homogeneous one by introducing in place of $u$ the new function $u-\phi$, see Section 4 of Chapter V in or Chapter 8 of for details on such considerations. We do not explore this here, because our goal is not to prove the most general result possible, but to provide a concrete setup in which we can prove the validity of the approximation results of interest.

Recall that the norms above are $L^{2}(X)$ norms in the respective space $X=\Omega_{T},\partial\Omega_{T}$ and $\Omega$ respectively. From Theorem 7.1, we have that

Each neural network $f^{n}$ satisfies the PDE

for some $h^{n},u^{n}_{0},$ and $g^{n}$ such that

For the purposes of this section, we make the following set of assumptions.

There is a constant $\mu>0$ and positive functions $\kappa(t,x),\lambda(t,x)$ such that for all $(t,x)\in\Omega_{T}$ we have

with $\kappa\in L^{2}(\Omega_{T})$, $\lambda\in L^{d+2+\eta}(\Omega_{T})$ for some $\eta>0$.

$\alpha(t,x,u,p)$ is differentiable with respect to $(x,u,p)$ with continuous derivatives.

There is a positive constant $\nu>0$ such that

$u_{0}(x)\in\mathcal{C}^{0,2+\xi}(\bar{\Omega})$ for some $\xi>0$In general, the Hölder space $\mathcal{C}^{0,\xi}(\bar{\Omega})$ is the Banach space of continuous functions in $\bar{\Omega}$ having continuous derivatives up to order $[\xi]$ in $\bar{\Omega}$ with finite corresponding uniform norms and finite uniform $\xi-[\xi]$ Hölder norm. Analogously, we also define the Hölder space $\mathcal{C}^{0,\xi,\xi/2}(\bar{\Omega}_{T})$ which in addition has finite $[\xi]/2$ and $(\xi-[\xi])/2$ regular and Hölder derivatives norms in time respectively. These spaces are denoted by $H^{\xi}(\bar{\Omega})$ and $H^{\xi,\xi/2}(\bar{\Omega}_{T})$ respectively in . with itself and its first derivative bounded in $\bar{\Omega}$.

The proof of this theorem is in the Appendix. We conclude this section with some remarks and an example.

Despite the restriction made to the zero boundary data case, we do expect that our results are also valid for reasonably smooth inhomogeneous boundary data. In addition, if we make further assumptions on the nonlinearities $\alpha(t,x,u,p)$ and $\gamma(t,x,u,p)$ and on the initial data $u_{0}(x)$, then one can establish existence and uniqueness of classical solutions, see for example Section 6 of Chapter V in for details. As a matter of fact the results of Chapter V.6 in show that with assuming a little bit more on the growth of the derivatives of the nonlinear functions $\alpha(t,x,u,p),\gamma(t,x,u,p)$ will lead to $\nabla_{x}u\in\mathcal{C}^{0,\delta^{\prime},\delta^{\prime}/2}(\Omega_{T})$ for some $\delta^{\prime}>0$. Furthermore, we remark here that stronger claims can be made if more properties are known in regards to the given approximating family $\{f^{n}\}$ such as, for example, a-priori bounds on appropriate Sobolev norms, but we do not explore this further here.

Let us present the case of linear parabolic PDEs in Example 7.6 below.

Let us assume that the operator $\mathcal{G}$ is linear in $u$ and $\nabla u$. In particular, let us set

and that the coefficients $b$ and $c$ are such that

where we recall for example $\left\lVert c\right\rVert_{q,r,\Omega_{T}}=\left(\int_{0}^{T}\left(\int_{\Omega}|c(t,x)|^{q}dx\right)^{r/q}\right)^{1/r}$ and $r,q$ satisfy the relations

In particular, the previous bounds always hold in the case of coefficients $b$ and $c$ that are bounded in $\Omega_{T}$. Under these conditions, standard results for linear PDE’s, see for instance Theorem 4.5 of Chapter III of for a related result, show that approximation results analogous to that of Theorem 7.3 hold.

Conclusion

We believe that deep learning could become a valuable approach for solving high-dimensional PDEs, which are important in physics, engineering, and finance. The PDE solution can be approximated with a deep neural network which is trained to satisfy the differential operator, initial condition, and boundary conditions. We prove that the neural network converges to the solution of the partial differential equation as the number of hidden units increases.

Our deep learning algorithm for solving PDEs is meshfree, which is key since meshes become infeasible in higher dimensions. Instead of forming a mesh, the neural network is trained on batches of randomly sampled time and space points. The approach is implemented for a class of high-dimensional free boundary PDEs in up to $200$ dimensions with accurate results. We also test it on a high-dimensional Hamilton-Jacobi-Bellman PDE with accurate results.

The DGM algorithm can be easily modified to apply to hyperbolic, elliptic, and partial-integral differential equations. The algorithm remains essentially the same for these other types of PDEs. However, numerical performance for these other types of PDEs remains to be be investigated.

It is also important to put the numerical results in Sections 4, 5 and 6 in a proper context. PDEs with highly non-monotonic or oscillatory solutions may be more challenging to solve and further developments in architecture will be necessary. Further numerical development and testing is therefore required to better judge the usefulness of deep learning for the solution of PDEs in other applications. However, the numerical results of this paper demonstrate that there is sufficient evidence to further explore deep neural network approaches for solving PDEs.

In addition, it would be of interest to establish results analogous to Theorem 7.3 for PDEs beyond the class of quasilinear parabolic PDEs considered in this paper. Stability analysis of deep learning and machine learning algorithms for solving PDEs is also an important question. It would certainly be interesting to study machine learning algorithms that use a more direct variational formulation of the involved PDEs. We leave these questions for future work.

Appendix A Proofs of convergence results

In this section we have gathered the proofs of the theoretical results of Section 7.

We have assumed that $(u,p)\mapsto\hat{\gamma}(t,x,u,p)$ is locally Lipschitz continuous in $(u,p)$ with Lipschitz constant that can have at most polynomial growth in $u$ and $p$ , uniformly with respect to $t,x$. This means that

for some constants $0\leq q_{1},q_{2},q_{3},q_{4}<\infty$. Therefore we obtain, using Hölder inequality with exponents $r_{1},r_{2}$,

where the unimportant constant $K<\infty$ may change from line to line and for two numbers $q_{1}\vee q_{3}=\max\{q_{1},q_{3}\}$. In the last step we used (A.1).

In addition, we have also assumed that for every $i,j\in\{1,\cdots d\}$, the mapping $(u,p)\mapsto\frac{\partial\alpha_{i}(t,x,u,p)}{\partial p_{j}}$ is locally Lipschitz in $(u,p)$ with Lipschitz constant that can have at most polynomial growth on $u$ and $p$, uniformly with respect to $t,x$. This means that

for some constants $0\leq q_{1},q_{2},q_{3},q_{4}<\infty$. Denote for convenience

Then, similarly to (A.2) we have after an application of Hölder inequality, for some constant $K<\infty$ that may change from line to line,

where in the last step we followed the computation in (A.2) and used (A.1).

Using (A.1) and (A.2)-(A.3) we subsequently obtain for the objective function (note that $\mathcal{G}[u](t,x)=0$ for $u$ that solves the PDE)

for an appropriate constant $K<\infty$. The last step completes the proof of the Theorem after rescaling $\epsilon$. ∎

Existence, regularity and uniqueness for (7.5) follows from Theorem 2.1 combined with Theorems 6.3-6.5 of Chapter V.6 in (see also Theorem 6.6 of Chapter V.6 of ). Boundedness follows from Theorem 2.1 in and Chapter V.2 in . The convergence proof follows by the smoothness of the neural networks together with compactness arguments as we explain below.

Next let us set $q=1+\frac{d}{d+4}\in(1,2)$ and note that for conjugates, $r_{1},r_{2}>1$ such that $1/r_{1}+1/r_{2}=1$

Let us choose $r_{2}=2/q>1$. Then we calculate $r_{1}=\frac{r_{2}}{r_{2}-1}=\frac{2}{2-q}$. Hence, we have that $r_{1}q=d+2$. Recalling the assumption $\lambda\in L^{d+2}(\Omega_{T})$ and the uniform bound on the $\left\lVert\nabla_{x}\hat{f}^{n}\right\rVert_{2}$ we subsequently obtain that for $q=1+\frac{d}{d+4}$, there is a constant $C<\infty$ such that

In preparation to passing to the limit as $n\rightarrow\infty$ in the weak formulation, we need to study the behavior of the nonlinear terms. Recalling the assumptions on $\alpha(t,x,u,p)$ we have for $\rho<2$ and for a measurable set $A\subset\Omega_{T}$ (the constant $K<\infty$ may change from line to line)

In the latter display we used Höder inequality with exponent $2/\rho>1$. By Vitali’s theorem we then conclude that

as $n\rightarrow\infty$, for every $1<\rho<2$. For the same reason, an analogous estimate to (A.4), gives

as $n\rightarrow\infty$, for $q=1+\frac{d}{d+4}$.

Notice also that by construction we have that the initial condition $u^{n}_{0}$ converges to $u_{0}$ strongly in $L^{2}(\Omega)$. The weak formulation of the PDE (7.6) with $g^{n}=0$ reads as follows. For every $t_{1}\in(0,T]$

for every $\phi\in C^{\infty}_{0}(\Omega_{T})$. Using the above convergence results, we then obtain that the limit point $u$ satisfies for every $t_{1}\in(0,T]$ the equation

which is the weak formulation of the equation (7.5).

References