{'fraction_non_alphanumeric': 0.0518992625573953, 'fraction_numerical': 0.017949074718241267, 'mean_word_length': 4.326581027667984, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Solving high-dimensional partial differential equations using deep learning

Developing algorithms for solving high-dimensional partial differential equations (PDEs) has been an exceedingly difficult task for a long time, due to the notoriously difficult problem known as the"curse of dimensionality". This paper introduces a deep learning-based approach that can handle general high-dimensional parabolic PDEs. To this end, the PDEs are reformulated using backward stochastic differential equations and the gradient of the unknown solution is approximated by neural networks, very much in the spirit of deep reinforcement learning with the gradient acting as the policy function. Numerical results on examples including the nonlinear Black-Scholes equation, the Hamilton-Jacobi-Bellman equation, and the Allen-Cahn equation suggest that the proposed algorithm is quite effective in high dimensions, in terms of both accuracy and cost. This opens up new possibilities in economics, finance, operational research, and physics, by considering all participating agents, assets, resources, or particles together at the same time, instead of making ad hoc assumptions on their inter-relationships.


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: 1. 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.
2. 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.
3. 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 [1]: 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., [2][3][4][5][6]). 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 [7] 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., [8,9]) 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 [10]), 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 [11] 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 [12,13], 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: Tr σσ T (t, x)(Hess x u)(t, x) + ∇u(t, x) · µ(t, x) + f t, x, u(t, x), σ T (t, x)∇u(t, x) = 0 (1) with some specified terminal condition u(T, x) = g(x). Here t and x represent the time and d-dimensional space variable respectively, µ is a known vector-valued function, σ is a known d × d matrix-valued function, σ T denotes the transpose associated to σ, ∇u and Hess x u denote the gradient and the Hessian of function u respect to x, Tr denotes the trace of a matrix, and f is a known nonlinear function. To fix ideas, we are interested in the solution at t = 0, x = ξ for some vector ξ ∈ R d .
Then the solution of (1) satisfies the following BSDE (cf., e.g., [8,9]): 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 ) ≈ θ u 0 , ∇u(0, X 0 ) ≈ θ ∇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 ∇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 < . . . < t N = T , we consider the simple Euler scheme for n = 1, . . . , N − 1: and where Given this temporal discretization, the path {X tn } 0≤n≤N can be easily sampled using (4). Our key step next is to approximate the function x → σ T (t, x) ∇u(t, x) at each time step t = t n by a multilayer feedforward neural network for n = 1, . . . , N − 1, where θ n denotes parameters of the neural network approximating 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, . . . , N − 1. Specifically, this network takes the paths {X tn } 0≤n≤N and {W tn } 0≤n≤N as the input data and gives the final output, denoted byû({X tn } 0≤n≤N , {W tn } 0≤n≤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: We can now use a stochastic gradient descent-type (SGD) algorithm to optimize the parameter θ, just as in the standard training of deep neural networks. In our numerical examples, we use the Adam optimizer [14]. 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 Nonlinear Black-Scholes Equation with Default Risk
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 [15]. 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., [16][17][18][19][20]). Each of these effects results in a nonlinear contribution in the pricing model (see, e.g., [17,21,22]). 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 [21].
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 [16,17]. 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 δ ∈ [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 [16]), 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 , γ h > γ l ): (see [17]). The associated nonlinear Black-Scholes equation We choose T = 1, δ = 2/3, R = 0.02,μ = 0.02,σ = 0.2, v h = 50, v l = 70, γ h = 0.2, γ l = 0.02, and the terminal condition g(x) = min{x 1 , . . . , x 100 } for x = (x 1 , . . . , x 100 ) ∈ R 100 . Fig. 1 shows the mean and the standard deviation of θ u 0 as an approximation of u(t=0, x=(100, . . . , 100)), with the final relative error being 0.46%. The not explicitly known "exact" solution of (11) at t = 0, x = (100, . . . , 100) has been approximately computed by means of the multilevel Picard method [11]: u(t=0, x=(100, . . . , 100)) ≈ 57.300. In comparison, if we do not consider the default risk, we getũ(t=0, x=(100, . . . , 100)) ≈ 60.781. In this case, the model becomes linear and can be solved using straightforward Monte Carlo methods. However, neglecting default risks results in a considerable error in the pricing, as illustrated above. The deep BSDE method allows us to rigorously incorporate default risks into pricing models. This in turn makes it possible to evaluate financial derivatives with substantial lower risks for the involved parties and the societies.

Hamilton-Jacobi-Bellman (HJB) Equation
The term "curse of dimensionality" was first used explicitly by Richard Bellman in the context of dynamic programming [1], 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 [23] 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: with t ∈ [0, T ], X 0 = x, and with the cost functional J({m t } 0≤t≤T ) = E T 0 m t 2 dt + g(X T ) .
Here {X t } t∈[0,T ] is the state process, {m t } t∈[0,T ] is the control process, λ is a positive constant representing the "strength" of the control, and {W t } t∈[0,T ] is a standard Brownian motion. Our goal is to minimize the cost functional through the control process. The HJB equation for this problem is given by ∂u ∂t (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.
We solve the Hamilton-Jacobi-Bellman equation (13) in the 100-dimensional case with g(x) = ln (1 + x 2 )/2 for x ∈ R 100 . Fig. 2 (a) shows the mean and the standard deviation of the relative error for u(t=0, x=(0, . . . , 0)) in the case λ = 1. The deep BSDE method achieves a relative error of 0.17% in a runtime of 330 seconds on a Macbook Pro. We also use the BSDE method to approximatively calculate the optimal cost u(t=0, x=(0, . . . , 0)) against different values of λ; see Fig. 2 (b). The curve in Fig. 2 (b) clearly confirms the intuition that the optimal cost decreases as the control strength increases.

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., [25]). Here we consider a typical Allen-Cahn equation with the "double-well potential" in 100-dimensional space with the initial condition u(0, x) = g(x), where g(x) = 1/ 2 + 0.4 x 2 for x ∈ R 100 . By applying a transformation of the time variable t → T − t (T > 0), we can turn (15) into the form of (1) such that the deep BSDE method can be used. Fig. 3 (a) shows the mean and the standard deviation of the relative error of u(t=0.3, x=(0, . . . , 0)). The not explicitly known "exact" solution of (15) at t = 0.3, x = (0, . . . , 0) has been approximatively computed by means of the branching diffusion method (see, e.g., [12,13]

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.
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
In this subsection we briefly illustrate the architecture of the deep BSDE method. To simplify the presentation we restrict ourselves in these illustrations to the case where the diffusion coefficient Fig. 4 illustrates the network architecture for the deep BSDE method. Note that ∇u(t n , X tn ) denotes the variable we approximate directly by sub-networks and u(t n , X tn ) denotes the variable we compute iteratively in the network. There are three types of connections in this network: 1. X tn → h 1 n → h 2 n → · · · → h H n → ∇u(t n , X tn ) is the multilayer feedforward neural network approximating the spatial gradients at time t = t n . The weights θ n of this sub-network are the parameters we aim to optimize.
2. (u(t n , X tn ), ∇u(t n , X tn ), W t n+1 − W tn ) → 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.
3. (X tn , W t n+1 − W tn ) → 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. It should be pointed out that the proposed deep BSDE method can also be employed if we are interested in values of the PDE solution u in a region D ⊂ R d at time t = 0 instead of at a single space-point ξ ∈ R d . In this case we choose X 0 = ξ to be a non-degenerate D-valued random variable and we employ two additional neural networks parameterized by {θ u 0 , θ ∇u 0 } for approximating the functions D x → u(0, x) ∈ R and D x → ∇u(0, x) ∈ R d . Upper and lower bounds for approximation errors of stochastic approximation algorithms for PDEs and BSDEs, respectively, can be found in e.g., [27][28][29] and the references therein.

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 [30] 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 [31] to implement our algorithm with the Adam optimizer [14] 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 [14] 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 [32] (d = 2): ∂u ∂t (t, x) + 1 2 ∆u(t, x) + min 1, u(t, x) − u * (t, x) 2 = 0, (19) in which u * (t, x) is the explicit oscillating solution Parameters are chosen in the same way as in [32]: κ = 1.6, λ = 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 0 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. mean of relative error 2.29% 0.90% 0.60% 0.56% 0.53% std. of relative error 0.0026 0.0016 0.0017 0.0017 0.0014 † The PDE is solved until convergence with 30 equidistant time steps (N =30) and 40000 iteration steps. Learning rate is 0.01 for the first half of iterations and 0.001 for the second half. † † We only count the layers that have free parameters to be optimized.