DGM: A deep learning algorithm for solving partial differential equations
Justin Sirignano Konstantinos Spiliopoulos


Efficient Tree Decomposition of High-Rank Tensors
    Adam S. Jermyn
    December 30, 2023
=================================================





High-dimensional PDEs have been a longstanding computational challenge. We propose to solve high-dimensional PDEs by approximating the solution with a deep neural network which is trained to satisfy the differential operator, initial condition, and boundary conditions. Our algorithm 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 algorithm is tested on a class of high-dimensional free boundary PDEs, which we are able to accurately solve in up to 200 dimensions. The algorithm is also tested on a high-dimensional Hamilton-Jacobi-Bellman PDE and Burgers' equation. The deep learning algorithm approximates the general solution to the Burgers' equation for a continuum of different boundary conditions and physical conditions (which can be viewed as a high-dimensional space). We call the algorithm a “Deep Galerkin Method (DGM)" since it is similar in spirit to Galerkin methods, with the solution approximated by a neural network instead of a linear combination of basis functions. In addition, we prove a theorem regarding the approximation power of neural networks for a class of quasilinear parabolic PDEs.





§ 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 𝒪^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 <ref>. Promising numerical results are presented later in Section <ref> for a class of high-dimensional free boundary PDEs. We also accurately solve a high-dimensional Hamilton-Jacobi-Bellman PDE in Section <ref> and Burger's equation in Section <ref>. DGM converts the computational cost of finite difference to a more convenient form: instead of a huge mesh of 𝒪^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 <ref> 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 <cit.>. 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 <ref>), and relevant literature (Section <ref>). The deep learning algorithm for solving PDEs is presented in Section <ref>. An efficient scheme for evaluating the diffusion operator is developed in Section <ref>.  Numerical analysis of the algorithm is presented in Sections <ref>, <ref>, and <ref>. 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 <ref>. Conclusions are in Section <ref>.  For readability purposes, proofs from Section <ref> have been collected in  Appendix <ref>.


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

    ∂_t u(t,x)  + ℒ u(t,x) = 0, ....... (t,x) ∈ [0,T] ×Ω
           u(0,x) = u_0(x), ........ x ∈Ω
           u(t,x) = g(t,x), ........ x ∈ [0,T] ×∂Ω,

where ∂Ω is the boundary of the domain Ω. The solution u(t,x) is of course unknown, but an approximate solution f(t,x) can be found by minimizing the L^2 error

    J( f )   =∂_t f + ℒf ^2_2, [0,T] ×Ω +     f -g^2_2, [0,T]  ×∂Ω +  f (0, ·) - u_0 ^2_2, Ω.

The error function J(f) measures how well the approximate solution f satisfies the differential operator, boundary condition, and initial condition. Note that no knowledge of the actual solution u is assumed; J(f) can be directly calculated from the PDE (<ref>) for any approximation f. The goal is to construct functions f for which  J(f)  is as close to  0 as possible. Define ℭ^n as the class of neural networks with a single hidden layer and n hidden units.[A neural network with a single hidden layer and n hidden units is a function of the form ℭ^n={h(t,x):ℝ^1+d↦ℝ: h(t,x)=∑_i=1^nβ_iψ(α_1,it+∑_j=1^dα_j,ix_j+c_j)} where Ψ: ℝ→ℝ is a nonlinear “activation" function such as a sigmoid or tanh function.] Let f^n be a neural network with n hidden units which minimizes J(f). We prove that, under certain conditions,

    there exists f^n∈ℭ^n such that   J(f^n)→ 0,  as n→∞,  and
           f^n→ u ....as.... n →∞,





strongly in, L^ρ([0,T] ×Ω), with ρ<2, for a class of quasilinear parabolic PDEs; see subsection <ref> and Theorem <ref> therein for the precise statement. That is, the neural network will converge in L^ρ, ρ<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 <ref>. 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) → 0 does not necessarily imply that f^n→ u, given that we only have L^2 control on the approximation error. First, we prove that J(f^n) → 0 as n →∞. We then establish that each neural network { f^n}_n=1^∞ 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→ u as n →∞ in L^ρ([0,T] ×Ω), for ρ<2,    using  the smoothness of the neural network approximations and compactness arguments.

Theorem <ref> 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 <ref>.







 §.§ Relevant Literature
 

Solving PDEs with a neural network as an approximation is a natural idea, and has been considered in various forms previously. <cit.>, <cit.>, <cit.>,  <cit.>, and <cit.> 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, <cit.> 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. <cit.> 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. <cit.> 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 <cit.> further develops the algorithm. The algorithm developed in <cit.> 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 <cit.>, seems to be able to recover accurately the entire solution (at least for the equations that we studied).  <cit.> 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, <cit.> has recently developed a novel partial differential equation approach to optimize deep neural networks.


<cit.> 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 <cit.>, <cit.>, and others.  Sparse grid methods have also been used to solve high-dimensional PDEs; see <cit.>, <cit.>, <cit.>, <cit.>, and <cit.>.

In regards to general results on the approximation power of neural networks we refer the interested reader to classical works <cit.> and we also mention the recent work by <cit.>, 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:

    .   ∂ u/∂ t(t,x) +  ℒ u(t,x)  = 0 , .... (t,x) ∈ [0,T] ×Ω, 
       .    u(t=0,x )  = u_0(x), 
       .    u(t,x) = g(t, x),  .... x ∈∂Ω,

where x ∈Ω⊂ℝ^d. The DGM algorithm approximates u(t,x) with a deep neural network f(t,x; θ) where θ∈ℝ^K are the neural network's parameters. Note that the differential operators ∂ f/∂ t(t, x; θ) and ℒ f(t,x; θ) can be calculated analytically.  Construct the objective function:

    J(f) =   ∂ f/∂ t(t,x; θ) +  ℒ f(t,x; θ)_[0,T] ×Ω, ν_1^2 +     f(t,x; θ ) - g(t,x) _ [0,T] ×∂Ω, ν_2^2 +  f(0, x; θ) - u_0(x) _Ω, ν_3^2.

Here, f(y) ^2_𝒴, ν = ∫_𝒴|f(y)|^2  ν(y) dy where ν(y) is a positive probability density on y ∈𝒴.  J(f) measures how well the function f(t,x; θ) satisfies the PDE differential operator, boundary conditions, and initial condition.  If J(f) = 0, then f(t,x; θ) is a solution to the PDE (<ref>).


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

Estimating θ by directly minimizing J(f) is infeasible when the dimension d is large since the integral over Ω 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 Ω and ∂Ω. This avoids ever forming a mesh.

The DGM algorithm is:
 

  * Generate random points (t_n, x_n) from [0,T] ×Ω and (τ_n, z_n) from [0,T] ×∂Ω according to respective probability densities ν_1 and ν_2.  Also, draw the random point w_n from Ω with probability density ν_3.

  * Calculate the squared error G(θ_n ,s_n) at the randomly sampled points s_n = { (t_n, x_n), ( τ_n, z_n), w_n } where:

    G(θ_n , s_n ) =   (∂ f/∂ t(t_n ,x_n; θ_n) +  ℒ f(t_n,x_n; θ_n) )^2 +   ( f(τ_n ,z_n ; θ_n ) - g(τ_n,z_n) )^2 + (  f(0, w_n; θ_n) - u_0(w_n) )^2.


  * Take a descent step at the random point s_n:

    θ_n+1 = θ_n - α_n ∇_θ G(θ_n, s_n)


  * Repeat until convergence criterion is satisfied.

The “learning rate" α_n decreases with n. The steps ∇_θ G(θ_n, s_n) are unbiased estimates of ∇_θ J( f(·;θ_n) ):

    𝔼[∇_θ G(θ_n, s_n) |θ_n  ]  = ∇_θ J(f(·; θ_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(·;θ_n+1) ) < J( f(·;θ_n) ) ), and θ_n+1 is therefore a better parameter estimate than θ_n.

Under (relatively mild) technical conditions (see <cit.>), the algorithm θ_n will converge to a critical point of the objective function J(f(·;θ)) as n →∞:

    lim_n →∞∇_θ J( f(·;θ_n) )   = 0.


It's important to note that θ_n may only converge to a local minimum when f(t, x; θ) 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 ℒ f(t,x; θ) contains second derivatives ∂^2 f/∂ x_i x_j (t,x; θ) 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 (<ref>)). The computational cost for calculating second derivatives (in both total arithmetic operations and memory) is 𝒪(d^2 × 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 𝒪(d × N). The cost associated with the second derivatives is further increased since we actually need the third-order derivatives ∇_θ∂^2 f/∂ x^2(t,x; θ) 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 ℒ f(t,x,; θ) is of the form 1/2∑_i,j=1^d ρ_i,jσ_i(x) σ_j(x)  ∂^2 f/∂ x_i x_j (t,x; θ), assume [ρ_i,j]_i,j=1^d is a positive definite matrix, and define σ(x) = (σ_1(x), …, σ_d(x) ). For example, such PDEs arise when considering expectations of functions of stochastic differential equations, where the σ(x) represents the diffusion coefficient. See equation (<ref>) 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,

    ∑_i,j=1^d ρ_i,jσ_i(x) σ_j(x)  ∂^2 f/∂ x_i x_j (t,x; θ) =  lim_Δ→ 0𝔼[∑_i=1^d  ∂ f/∂ x_i(t, x + σ(x)W_Δ; θ ) - ∂ f/∂ x_i(t,x; θ)  /Δσ_i(x)  W_Δ^i  ] ,

where W_t ∈ℝ^d is a Brownian motion and Δ∈ℝ_+ is the step-size. The convergence rate for (<ref>) is 𝒪(√(Δ)).[Let f be a three-times differentiable function in x with bounded third-order derivatives in x. Then, it directly follows from a Taylor expansion that |∑_i,j=1^d ρ_i,jσ_i(x) σ_j(x)  ∂^2 f/∂ x_i x_j (t,x; θ) -  𝔼[∑_i=1^d  ∂ f/∂ x_i(t, x + σ(x)W_Δ; θ ) - ∂ f/∂ x_i(t,x; θ)  /Δσ_i(x)  W_Δ^i  ]|≤ C(x) √(Δ). The constant C(x) depends upon ρ, f_xxx(t,x; θ) and σ(x).] Define:

    G_1(θ_n , s_n )     =   (∂ f/∂ t(t_n ,x_n; θ_n) +  ℒ f(t_n,x_n; θ_n) )^2, 
    
    G_2(θ_n , s_n )     =   ( f(τ_n ,z_n ; θ_n ) - g(τ_n ,z_n) )^2 , 
    
    G_3(θ_n , s_n )      =   (  f(0, w_n; θ_n) - u_0(w_n) )^2, 
    
    G(θ_n , s_n )      =    G_1(θ_n , s_n )  + G_2(θ_n , s_n )  + G_3(θ_n , s_n ) .

The DGM algorithm use the gradient ∇_θ G_1 (θ_n , s_n ), which requires the calculation of the second derivative terms in ℒ f(t_n,x_n; θ_n).  Define the first derivative operators as

    ℒ_1 f(t_n,x_n; θ_n)    =ℒ f(t_n,x_n; θ_n) - 1/2∑_i,j=1^d ρ_i,jσ_i(x_n) σ_j(x_n)  ∂^2 f/∂ x_i x_j (t_n,x_n; θ).

  Using (<ref>), ∇_θ G_1 is approximated as G̃_1 with a fixed constant Δ > 0:

    G̃_1(θ_n, s_n)     =   2(∂ f/∂ t(t_n ,x_n; θ_n) +  ℒ_1 f(t_n,x_n; θ_n)  +  1/2∑_i=1^d  ∂ f/∂ x_i(t, x_n + σ(x_n) W_Δ; θ ) - ∂ f/∂ x_i(t,x_n; θ)  /Δσ_i(x_n)   W_Δ^i )
       ×   ∇_θ(∂ f/∂ t (t_n ,x_n; θ_n) +  ℒ_1 f(t_n,x_n; θ_n)  + 1/2∑_i=1^d  ∂ f/∂ x_i(t, x_n + σ(x_n) W̃_Δ; θ ) - ∂ f/∂ x_i(t,x_n; θ)  /Δσ_i(x_n)  W̃_Δ^i  ),

where W_Δ is a d-dimensional normal random variable with 𝔼[W_Δ ] = 0 and Cov[ ( W_Δ )_i, ( W_Δ )_j ] = ρ_i,jΔ.  W̃_Δ has the same distribution as W_Δ.  W_Δ and W̃_Δ are independent.  G̃_1(θ_n, s_n) is a Monte Carlo approximation of ∇_θ G_1(θ_n, s_n).  G̃_1(θ_n, s_n) has 𝒪(√(Δ)) bias as an approximation for   ∇_θ G_1(θ_n, s_n).  This approximation error can be further improved via the following scheme using “antithetic variates":


    G̃_1(θ_n, s_n)      =   G̃_1,a(θ_n, s_n)  +  G̃_1,b(θ_n, s_n)   
    G̃_1,a (θ_n, s_n)     =   (∂ f/∂ t(t_n ,x_n; θ_n) +  ℒ_1 f(t_n,x_n; θ_n)  +  1/2∑_i=1^d ∂ f/∂ x_i(t, x_n + σ(x_n)  W_Δ; θ ) - ∂ f/∂ x_i(t,x_n; θ)  /Δσ_i(x_n)   W_Δ^i )
       ×   ∇_θ(∂ f/∂ t (t_n ,x_n; θ_n) +  ℒ_1 f(t_n,x_n; θ_n)  + 1/2∑_i=1^d ∂ f/∂ x_i(t, x_n + σ(x_n) W̃_Δ; θ ) - ∂ f/∂ x_i(t,x_n; θ)  /Δσ_i(x_n) W̃_Δ^i  ), 
    G̃_1,b (θ_n, s_n)     =   (∂ f/∂ t(t_n ,x_n; θ_n) +  ℒ_1 f(t_n,x_n; θ_n)  - 1/2∑_i=1^d  ∂ f/∂ x_i(t, x_n - σ(x_n)  W_Δ; θ ) - ∂ f/∂ x_i(t,x_n; θ)  /Δσ_i(x_n)  W_Δ^i )
       ×   ∇_θ(∂ f/∂ t (t_n ,x_n; θ_n) +   ℒ_1 f(t_n,x_n; θ_n)  - 1/2∑_i=1^d ∂ f/∂ x_i(t, x_n - σ(x_n) W̃_Δ; θ ) - ∂ f/∂ x_i(t,x_n; θ)  /Δσ_i(x_n)  W̃_Δ^i  ).

The approximation (<ref>) has 𝒪( Δ ) bias as an approximation for ∇_θ G_1(θ_n, s_n). (<ref>)  uses antithetic variates in the sense that G̃_1,a(θ_n, s_n) uses the random variables (W_Δ, W̃_Δ) while G̃_1,b (θ_n, s_n) uses (-W_Δ, -W̃_Δ). See <cit.> for a background on antithetic variates in simulation algorithms. A Taylor expansion can be used to show the approximation error is 𝒪( Δ ). It is important to highlight that there is no computational cost associated with the magnitude of Δ; an arbitrarily small Δ can be chosen with no additional computational cost (although there may be numerical underflow or overflow problems). The modified algorithm using the Monte Carlo approximation for the second derivatives is:



  * Generate random points (t_n, x_n) from [0,T] ×Ω and (τ_n, z_n) from [0,T] ×∂Ω according to respective densities ν_1 and ν_2.  Also, draw the random point w_n from Ω with density ν_3.

  * Calculate the step G̃(θ_n ,s_n) =  G̃_1(θ_n , s_n )  + ∇_θ G_2(θ_n , s_n )  + ∇_θ G_3(θ_n , s_n ) at the randomly sampled points s_n = { (t_n, x_n), ( τ_n, z_n), w_n }. G̃(θ_n ,s_n) is an approximation for ∇_θ G(θ_n ,s_n).

  * Take a step at the random point s_n:

    θ_n+1 = θ_n - α_n G̃(θ_n, s_n)


  * Repeat until convergence criterion is satisfied.


In conclusion, the modified algorithm here is computationally less expensive than the original algorithm in Section <ref> 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 <ref> 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 <ref>. 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
 
We test our algorithm on a class of high-dimensional free boundary PDEs. These free boundary PDEs are used in finance to price American options and are often referred to as “American option PDEs". An American option is a financial derivative on a portfolio of stocks. The option owner may at any time t ∈ [0,T] choose to exercise the American option and receive a payoff which is determined by the underlying prices of the stocks in the portfolio. T is called the maturity date of the option and the payoff function is g(x): ℝ^d →ℝ. Let X_t ∈ℝ^d be the prices of d stocks.
If at time t the stock prices X_t = x, the price of the option is u(t,x).  The price function u(t,x) satisfies a free boundary PDE on [0,T] ×ℝ^d. For American options, one is primarily interested in the solution u(0, X_0) since this is the fair price to buy or sell the option.

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 <ref> 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 <ref> with an iterative method;  see Section <ref>. Section <ref> describes the architecture and implementation details for the neural network. Section <ref> reports numerical accuracy for a case where a semi-analytic solution exists.  Section <ref> reports numerical accuracy for a case where no semi-analytic solution exists.






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

    d X_t^i    =   μ(X_t^i) dt +σ(X_t^i) d W_t^i, 
    
    u(t,x)    =   sup_τ≥ t𝔼[ e^- r( τ∧ T) g(X_τ∧ T) | X_t = x],

where W_t ∈ℝ^d is a standard Brownian motion and Cov[ dW_t^i, dW_t^j] = ρ_i,j dt.  The price of the American option is u(0, X_0).

The model (<ref>) for the stock price dynamics is widely used in practice and captures several desirable characteristics. First, the drift μ(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 σ(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 ρ_i,j. An example is the well-known Black-Scholes model μ(x) = μ x and σ(x) = σ x. In the Black-Scholes model, the average rate of return for each stock is μ.

An American option is a financial derivative which the owner can choose to “exercise" at any time t ∈ [0,T]. If the owner exercises the option, they receive the financial payoff g(X_t) where X_t is the prices of the underlying stocks. If the owner does not choose to exercise the option, they receive the payoff g(X_T) at the final time T. The value (or price) of the American option at time t is u(t,X_t). Some typical examples of the payoff function g(x): ℝ^d →ℝ are g(x)=max( (∏_i=1^d x_i )^1/d- K, 0 ) and g(x)=max(1/d∑_i=1^d x_i - K, 0 ). The former is referred to as a “geometric payoff function" while the latter is called an “arithmetic payoff function." K is the “strike price" and is a positive number.


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

    0    =   ∂ u/∂ t(t,x)  + μ(x) ·∂ u/∂ x(t,x)  + 1/2∑_i,j =1^d ρ_i,jσ(x_i) σ(x_j) ∂^2 u/∂ x_i ∂ x_j (t,x) -r u(t,x),  ....∀..{ (t,x) : u(t,x) > g(x) }.  
    
    u(t,x)    ≥     g(x),  ....∀.. (t,x).  
    
    u(t,x)    ∈    C^1(ℝ_+×ℝ^d), ....∀..{ (t,x) : u(t,x) = g(x) }. 
    
    u(T,x)    =     g(x),  ....∀.. x.

The free boundary set is F = { (t,x) : u(t,x) = g(x) }.  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 (<ref>) 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 θ_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 θ_n input of the probability density used to simulate random points.  For this purpose, define the objective function:

    J(f;θ, θ̃)     =   ∂ f/∂ t(t,x; θ)  + μ(x) ·∂ f/∂ x(t,x; θ)  + 1/2∑_i,j =1^d ρ_i,jσ(x_i) σ(x_j) ∂^2 f/∂ x_i ∂ x_j (t,x; θ) -r f(t,x; θ) _[0,T] ×Ω, ν_1(θ̃) ^2 
       +   max ( g(x) - f(t,x; θ), 0 ) _ [0,T] ×Ω, ν_2(θ̃) ^2  
       +     f(T,x; θ ) - g(x) _Ω, ν_3 ^2.

Descent steps are taken in the direction -∇_θ J( f;θ, θ̃).  ν_1(θ̃) and ν_2(θ̃) are the densities of the points in B̃^1 and B̃^2, which are defined below.  The deep learning algorithm is:

 

  * Generate the random batch of points B^1 = { t_m, x_m }_m=1^M from [0,T] ×Ω according to the probability density ν_1^0.  Select the points B̃^1 = { (t,x) ∈ B^1 : f(t,x; θ_n) > g(x) }.

  * Generate the random batch of points B^2 = {τ_m, z_m }_m=1^M from [0,T] ×∂Ω according to the probability density ν_2^0.  Select the points B̃^2 = { (τ,z) ∈ B^2 : f(τ, z ;  θ_n) ≤ g(z) }.

  * Generate the random batch of points B^3 = { w_m }_m=1^M from Ω with probability density ν_3.

  * Approximate J(f;θ_n, θ̃_n) as J(f;θ_n, S_n) at the randomly sampled points S_n = {B̃^1, B̃^2, B^3 }:

    J(f;θ_n , S_n)    =   1/|B̃^1 | ∑_(t_m, x_m) ∈B̃^1 (∂ f/∂ t(t_m,x_m; θ_n)  + μ(x_m) ·∂ f/∂ x(t_m,x_m; θ_n)   
       +   1/2∑_i,j =1^d ρ_i,jσ(x_i) σ(x_j) ∂^2 f/∂ x_i ∂ x_j (t_m,x_m; θ_n) -r f(t_m,x_m; θ_n) )^2  
       +   1/ | B̃^2 | ∑_(τ_m, z_m) ∈B̃^2 max( g(z_m) - f(τ_m,z_m; θ_n), 0 )^2 
       +   1/| B^3 |∑_ w_m ∈ B^3( f(T,w_m; θ ) - g(w_m)  )^2.


  * Take a descent step for the random batch S_n:

    θ_n+1 = θ_n - α_n ∇_θ J(f;θ_n, S_n).


  * Repeat until convergence criterion is satisfied.


The second derivatives in the above algorithm can be approximated using the method from Section <ref>.



 §.§ 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; θ) 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:


    S^1   =   σ( W^1 → x + b^1 ), 
    
    Z^ℓ   =   σ ( U^z,ℓ→ x  + W^z,ℓ S^ℓ + b^z,ℓ ), ....ℓ = 1, …, L, 
    
    G^ℓ   =   σ ( U^g,ℓ→ x  + W^g,ℓ S^1 + b^g,ℓ ), ....ℓ = 1, …, L, 
    
    R^ℓ   =   σ ( U^r,ℓ→ x  + W^r,ℓ S^ℓ + b^r, ℓ ), ....ℓ = 1, …, L, 
    
    H^ℓ   =   σ ( U^h,ℓ→ x  + W^h,ℓ  ( S^ℓ⊙ R^ℓ ) + b^h, ℓ ), ....ℓ = 1, …, L, 
    
    S^ℓ+1   =     (1 - G^ℓ ) ⊙ H^ℓ + Z^ℓ⊙ S^ℓ, ....ℓ = 1, …, L, 
    
    f(t,x; θ)    =    W S^L+1 + b,

where → x  = (t,x), the number of hidden layers is L+1, and ⊙ denotes element-wise multiplication (i.e., z ⊙ v = ( z_0 v_0, …, z_N v_N )). The parameters are

    θ = { W^1, b^1, ( U^z, ℓ, W^z, ℓ, b^z, ℓ)_ℓ=1^L,  ( U^g, ℓ, W^g, ℓ, b^g, ℓ)_ℓ=1^L, ( U^r, ℓ, W^r, ℓ, b^r, ℓ)_ℓ=1^L, ( U^h, ℓ, W^h, ℓ, b^h, ℓ)_ℓ=1^L, W, b }.

The number of units in each layer is M and σ: ℝ^M →ℝ^M is an element-wise nonlinearity:

    σ(z) = (ϕ(z_1), ϕ(z_2), …, ϕ(z_M) ),

where ϕ: ℝ→ℝ is a nonlinear activation function such as the tanh function, sigmoidal function e^y/1 + e^y, or rectified linear unit (ReLU) max(y,0). The parameters in θ have dimensions W^1 ∈ℝ^M × (d+1), b^1 ∈ℝ^M, U^z,ℓ∈ℝ^M × (d +1), W^z, ℓ∈ℝ^M × M, b^z, ℓ∈ℝ^M, U^g,ℓ∈ℝ^M × (d +1), W^g, ℓ∈ℝ^M × M, b^g, ℓ∈ℝ^M, U^r,ℓ∈ℝ^M × (d +1), W^r, ℓ∈ℝ^M × M, b^r, ℓ∈ℝ^M, U^h, ℓ∈ℝ^M × (d +1), W^h, ℓ∈ℝ^M × M, b^h, ℓ∈ℝ^M, W ∈ℝ^1 × M, and b ∈ℝ.

The architecture (<ref>) 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 (<ref>) is similar to the architecture for LSTM networks (see <cit.>) and highway networks (see <cit.>). 

The key hyperparameters in the neural network (<ref>) are the number of layers L, the number of units M in each sub-layer, and the choice of activation unit ϕ(y). We found for the applications in this paper that the hyperparameters L = 3 (i.e., four hidden layers), M = 50, and ϕ(y) = tanh(y) were effective. It is worthwhile to note that the choice of ϕ(y) = tanh(y) means that f(t,x; θ) is smooth and therefore can solve for a “classical solution" of the PDE. The neural network parameters are initialized using the Xavier initialization (see <cit.>). The architecture (<ref>) is bounded in the input x (for a fixed choice of parameters θ) if σ(·) is a tanh or sigmoidal function; it may be helpful to allow the network to be unbounded for approximating unbounded/growing functions. We found that replacing the σ(·) in the H^ℓ sub-layer with the identity function can be an effective way to develop an unbounded network.

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.

A regularization term (such as an ℓ^2 penalty) could also be included in the objective function for the algorithm. Such regularization terms are used for reducing overfitting in machine learning models estimated using datasets which have a limited number of data samples. (For example, a model estimated on a dataset of 60,000 images.) However, it's unclear if this will be helpful in the context of this paper's application, since there is no strict upper bound on the size of the dataset (i.e., one can always simulate more time/space points).

Our computational approach to training the neural network involved several components. The second derivatives are approximated using the method from Section <ref>. 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 <cit.>) 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 <ref> 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 <cit.>.



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:

 
    α_n = {[                  10^-4              n ≤ 5,000;              5 × 10^-4     5,000 < n ≤ 10,000;                  10^-5 <  10,000 < n ≤ 20,000;              5 × 10^-6    20,000 < n ≤ 30,000;                  10^-6    30,000 < n ≤ 40,000;              5 × 10^-7    40,000 < n ≤ 45,000;                  10^-7            45,000 <  n;                        ].


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 (<ref>) with respect to x or θ. 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. 








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

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





The semi-analytic solution used in Table <ref> is provided below. Let μ(x) = (r -c) x, σ(x) = σ x, and ρ_i,j = ρ for i ≠ j (i.e., the Black-Scholes model). If the payoff function in (<ref>) is g(x)=max( (∏_i=1^d x_i )^1/d- K, 0 ), then there is a semi-analytic solution to (<ref>):

    u(t,x)    =    v(t, (∏_i=1^d x_i )^1/d- K ),

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

    0    =   ∂ v/∂ t(t,x)  + μ̂x ∂ v/∂ x(t,x)  + 1/2σ̂^2 x ∂^2 v/∂ x^2(t,x) -r v(t,x),  ....∀..{ (t,x) : v(t,x) > ĝ(x) }.  
    
    v(t,x)    ≥   ĝ(x),  ....∀.. (t,x).  
    v̂(t,x)    ∈    C^1(ℝ_+×ℝ^d), ....∀..{ (t,x) : v(t,x) = ĝ(x) }. 
    
    v(T,x)    =   ĝ(x),  ....∀.. x,

where σ̂^2 =  d σ^2 + d(d-1) ρσ^2 /d^2, μ̂= (r-c) - 1/2σ̂^2 + 1/2σ^2, and ĝ(x) = max(x, 0). The one-dimensional PDE (<ref>) can be solved using finite difference methods. If f(t,x; θ) is the deep learning algorithm's estimate for the PDE solution at (t,x), the relative error at the point
(t,x) is | f(t, x; θ) - u(t, x) |/|u(t,x) | × 100% and the absolute error at the point (t,x) is | f(t, x; θ) - u(t, x) |. The relative error and absolute error at the point (t,x) can be evaluated using the semi-analytic solution (<ref>).





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) ∈ [0, T] ×Ω. As an example, we present in Figure <ref> 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:


  * Sample time points t^ℓ uniformly on [0,T] and sample spatial points x^ℓ = (x_1^ℓ, …, x_20^ℓ) from the joint distribution of X_t^1, …, X_t^20 in equation (<ref>). This produces an “envelope" of sampled points since X_t spreads out as a diffusive process from X_0 = 1.

  * Calculate the error E^ℓ at each sampled point (t^ℓ, x^ℓ) for ℓ =1, …, L.

  * Aggregate the error over a two-dimensional subspace ( t^ℓ,  (∏_i=1^20 x_i^ℓ )^1/20, E^ℓ) for ℓ =1, …, L.

  * Produce a contour plot from the data (  t^ℓ, (∏_i=1^20 x_i^ℓ )^1/20, E^ℓ)_ℓ=1^L. The x-axis is t and the y-axis is the geometric average (∏_i=1^20 x_i )^1/20, which corresponds to the final condition g(x).

Figure <ref> reports both the absolute error and the percent error. The percent error | f(t,x; θ) - u(t,x) |/ | u(t,x) | × 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.
















 §.§ 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 ρ_i,j = 0.75 and σ = 0.25 for all stocks. This section solves a more challenging heterogeneous case where ρ_i,j and σ_i vary across all dimensions i =1, 2, …, 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 d^2 +d/2 second derivative terms.  The correlation coefficients ρ_i,j range from -0.53 to 0.80 for i ≠ j and σ_i ranges from 0.09 to 0.69.



Let f(t,x; θ) be the neural network approximation.  <cit.> derived that the PDE solution u(t,x) lies in the interval:

    u(t,x)    ∈   [ u(t,x), u(t,x) ], 
    u(t,x)    =   𝔼 [ g( X_τ ) | X_t = x, τ > t  ], 
    u(t,x)    =   𝔼 [  sup_s ∈ [t,T] [ e^-r (s-t)  g( X_s)  - M_s  ]  ].

where τ = inf{ t ∈ [0,T]: f(t, X_t; θ) < g(X_t ) } and M_s is a martingale constructed from the approximate solution f(t,x; θ)

    M_s    =    f(s,X_s; θ) - f(t,X_t; θ)  - ∫_t^s [∂ f/∂ t(s',X_s'; θ)  + μ(X_s') ∂ f/∂ x(s',X_s'; θ)  
       +   1/2∑_i,j =1^d σ(X_s', i) σ(X_s', j) ∂^2 f/∂ x_i ∂ x_j (s',X_s'; θ) -r f(s',X_s'; θ) ] ds'.

The bounds (<ref>) depend only on the approximation f(t,x; θ), 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 [u(0, X_0), u(0,X_0) ], which has an error bound of u(0,X_0) - u(0, X_0)  /  2  u(0, X_0) × 100%. Numerical results are in Table <ref>.






We present in Figure <ref> 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:


  * Sample time points t^ℓ uniformly on [0,T] and sample spatial points x^ℓ = (x_1^ℓ, …, x_20^ℓ) from the joint distribution of X_t^1, …, X_t^20 in equation (<ref>).

  * Calculate the error E^ℓ at each sampled point (t^ℓ, x^ℓ) for ℓ =1, …, L.

  * Aggregate the error over a two-dimensional subspace ( t^ℓ,  1/20∑_i=1^20 x_i^ℓ , E^ℓ) for ℓ =1, …, L.

  * Produce a contour plot from the data (  t^ℓ,  1/20∑_i=1^20 x_i^ℓ, E^ℓ)_ℓ=1^L. The x-axis is t and the y-axis is the geometric average 1/20∑_i=1^20 x_i^ℓ, which corresponds to the final condition g(x).

Figure <ref> reports both the absolute error and the percent error. The percent error | f(t,x; θ) - u(t,x) |/ | u(t,x) | × 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 <ref> requires significant computations. For each point at which calculate an error bound, a new simulation of (<ref>) 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 (<ref>). The PDE (<ref>) is motivated by the problem of optimally controlling the stochastic partial differential equation (SPDE):

    ∂ v/∂ t(t,x)    =   α∂^2 v/∂ x^2(t,x)  + u(x)  + σ∂^2W/∂ t∂ x (t,x), .... x ∈ [0,L],
    
    v(t,x=0)    =   v̅(0), 
    
    v(t, x=L)    =   v̅(L), 
    
    v(t=0,x)    =    v_0(x),

where u(x) is the control and W(t,x) is   a Brownian sheet (i.e., ∂^2W/∂ t∂ x (t,x) is space-time white noise) defined on a stochastic basis (Ω,ℱ,ℱ_t,ℙ).  The
 square integrable, adapted to the filtration ℱ_t, control u is a source/sink term which can be used to guide the temperature v(t,x) towards a target profile v̅(x) on [0,L]. As it is discussed in <cit.> such problems admit unique solutions in the appropriate generalized sense, see Theorem 3.1 in <cit.>. The endpoints at x=0,L are held at the target temperatures. Specifically, the optimal control minimizes

    𝔼[∫_0^∞ e^- γ s∫_0^L [ ( v(s,x) - v̅ (x) )^2 + λ u(x)^2  ] dx ds  ].

 The constant γ > 0 is a discount factor. The constant λ > 0 penalizes large values for the control u(x). The goal is to reach the target 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  <cit.> as well as <cit.> and <cit.> for an analysis of infinite-dimensional HJB equations for the stochastic heat equation.

 An example of a problem represented by the SPDE (<ref>) 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 (<ref>).

 (<ref>) can be discretized in space, which yields a system of stochastic differential equations (SDEs). (For example, see Section 3.2 of <cit.>.) 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.

    d X_t^j = α/Δ^2   ( X_t^j+1 - 2 X_t^j + X_t^j-1 ) dt + U^j_t dt + σ/√(Δ)  d W_t^j,  X^j_0=v_0(jΔ),

where Δ is the mesh size, v(t,j Δ) = X_t^j, u(j Δ) = U^j_t, and W_t^j are independent standard Brownian motions (see <cit.>, <cit.>, and <cit.> regarding numerical schemes for stochastic parabolic PDEs of the form considered in this section). The dimension of the SDE system (<ref>) is d = L/Δ - 1. Note that (<ref>) uses a central difference scheme for the diffusion term in (<ref>).

The objective function (<ref>) becomes:

    V( x) = inf_U_t ∈𝒰𝔼[∫_0^∞ e^- γ s∑_j=1^d [ ( X^j_s - v̅ (j Δ) )^2 + λ (U^j_s)^2  ]Δ ds | X_0  = x ].

The value function V(x) satisfies a nonlinear PDE with d spatial dimensions x_1, x_2, …, x_d.

    0    =   Δ (x - v̅)^⊤ (x - v̅) - 1/4 λΔ∑_j=1^d (∂ V/∂ x_j(x) )^2 
       +   σ^2/2 Δ∑_j=1^d  ∂^2 V/∂ x_j^2 (x) + α/Δ^2∑_j=1^d   ( x_j+1 - 2 x_j + x_j-1 ) ∂ V/∂ x_j(x) - γ V(x).

The vector v̅ = ( v̅ ( Δ), v̅ (2 Δ), …, v̅ (d Δ) ). Note that the values x_d+1 = v̅(L) and x_0 = v̅(0) are constants which correspond to the boundary conditions in (<ref>). The PDE (<ref>) is high dimensional since the number of dimensions d = L/Δ - 1. The optimal control is

    U_t^j  = - 1/2 λΔ∂ V/∂ x_j(X_t).

We solve the PDE (<ref>) using the deep learning algorithm for d = 21 dimensions. The size of the domain is L = 10^-1. The coefficients are α = 10^-4, σ = 10^-1/2, λ = 1, and γ = 1. The target profile is v̅(x) = 0.



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


  * Sample spatial points x^ℓ = (x_1^ℓ, …, x_21^ℓ) from the distribution of (<ref>) for ℓ =1, …, L.

  * Calculate the percent error at each sampled point. The percent error is A^ℓ =  | f(x^ℓ; θ) - V(x^ℓ) |/ | V(x^ℓ) | × 100%.

  * Aggregate the accuracy over a two-dimensional subspace ( x_11^ℓ, 1/21∑_i=1^21 x_i^ℓ, A^ℓ) for ℓ =1, …, L.

  * Produce a contour plot from the data ( x_11^ℓ, 1/21∑_i=1^21 x_i^ℓ, A^ℓ)_ℓ=1^L. The x-axis is x_11 and the y-axis is the average 1/21∑_i=1^21 x_i. This corresponds to v(t,x) at the midpoint x = L/2 and the average 1/L∫_0^L v(t,x) dx, respectively.

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



Lastly, we close this section by mentioning that in the recent paper <cit.> (see also <cit.>) 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 𝒫, 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

    ∂ u/∂ t(t,x; p)    =   ℒ_p u(t,x; p), .... (t,x) ∈ [0,T] ×Ω, 
    
    u(t,x;p)    =    g_p(x),  .... (t,x) ∈ [0,T] ×∂Ω, 
    
    u(t=0,x;p)    =    h_p(x), .... x ∈Ω.


A traditional approach would be to discretize the 𝒫-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 𝒫 is typically high-dimensional.

We propose to use the DGM algorithm to approximate the general solution to the PDE (<ref>) 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,



  * Initialize θ.

  * Repeat until convergence:


  * Generate random samples (t, x, p) from [0,T] ×Ω×𝒫,  (t̃, x̃) from [0,T] ×∂Ω, and
x̂ from Ω.

  * Construct the objective function

    J(θ)    =   (∂ f/∂ t(t,x,p; θ) -  ℒ_p f(t,x,p;  θ) )^2 
       +   ( g_p(x̃) - f( t̃,x̃,p; θ) )^2 
       +   ( h_p(x̂) - f( 0,x̂,p; θ) )^2.


  * Update θ with a stochastic gradient descent step

    θ⟶θ - α∇_θ J( θ ),

where α is the learning rate.



If x is low-dimensional (d ≤ 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.


    .   ∂ u/∂ t =  ν∂^2 u/∂ x^2 - α u ∂ u/∂ x, .... (t,x) ∈ [0,1] × [0,1], 
       .     u(t, x =0) = a, 
       .    u(t, x=1)  = b, 
       .    u(t=0, x) = g(x), .... x ∈ [0,1].


The problem setup space is 𝒫 = (ν, α, a, b) ∈ℝ^4. The initial condition g(x) is chosen to be a linear function which matches the boundary conditions u(t, x =0) = a and u(t, x =1) = b. We train a single neural network to approximate the solution of u(t, x; p) over the entire space (t,x, ν, α, a, b) ∈ [0,1] × [0,1] × [10^-2, 10^-1 ] × [10^-2, 1 ] × [-1, 1] × [-1, 1]. We use a larger network (6 layers, 200 units per layer) than in the previous numerical examples. Figure <ref> compares the deep learning solution with the exact solution for several different problem setups p. The solutions are very close; in several cases, the two solutions are visibly indistinguishable. The deep learning algorithm is able to accurately capture the shock layers and boundary layers.

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

























§ NEURAL NETWORK APPROXIMATION THEOREM FOR PDES
  
Let the L^2 error J(f) measure how well the neural network f satisfies the differential operator, boundary condition, and initial condition. Define ℭ^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





    there exists f^n∈ℭ^n such that   J(f^n)→ 0,  as n→∞,  and
           f^n→ u ....as.... n →∞,





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.

    J(f^n) → 0 .....as..... n →∞.

Let u be the solution to the PDE. The statement (<ref>) does not necessarily imply that f^n→ 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^∞ 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→ u as n →∞ 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 <ref> proves that J(f^n) → 0 as n →∞. Section <ref> contains convergence results of f^n to the solution u of the PDE as n →∞. The main result is Theorem <ref>. For readability purposes the corresponding proofs are in Appendix <ref>.




 §.§ Convergence of the L^2 error J(f)
 
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 <cit.> 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 <cit.> apply).

Consider a bounded set Ω⊂ℝ^d with a smooth boundary ∂Ω and denote Ω_T=(0,T]×Ω and ∂Ω_T=(0,T]×∂Ω. In this subsection we consider the class of quasilinear parabolic PDE's of the form 


    ∂_t u(t,x)   -div(α(t,x,u(t,x), ∇ u(t,x)))+γ(t,x,u(t,x),∇ u(t,x))=0,  for (t,x)∈Ω_T
    
    u(0,x)   =u_0(x),  for x∈Ω
    
    u(t,x)   =g(t,x),  for (t,x)∈∂Ω_T


For notational convenience, let us write the operator of (<ref>) as 𝒢. Namely, let us denote

    𝒢[u](t,x)   =∂_t u(t,x)-div(α(t,x,u(t,x), ∇ u(t,x)))+γ(t,x,u(t,x),∇ u(t,x)).


Notice that we can write

    𝒢[u](t,x)   =∂_t u(t,x)-∑_i,j=1^d∂α_i(t,x,u(t,x), ∇ u(t,x))/∂ u_x_j∂_x_i,x_ju(t,x)+γ̂(t,x,u(t,x),∇ u(t,x)),

where

    γ̂(t,x,u,p)=γ(t,x,u,p)-∑_i=1^d∂α_i(t,x,u,p)/∂ u∂_x_iu -∑_i=1^d∂α_i(t,x,u,p)/∂ x_i.


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

In particular we assume that there is a unique u(t,x) solving (<ref>) such that

    u(t,x)∈𝒞(Ω̅_T)⋂𝒞^1+η/2,2+η(Ω_T)  with η∈(0,1)  and that sup_(t,x)∈Ω_T∑_k=1^2|∇_x^(k)u(t,x)|<∞.


We refer the interested reader to Theorems  5.4, 6.1 and 6.2 of Chapter V in <cit.> for specific general conditions on α,γ guaranteeing the validity of the aforementioned statement.

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

    ℭ^n(ψ)={ζ(t,x):ℝ^1+d↦ℝ: ζ(t,x)=∑_i=1^nβ_iψ(α_1,it+∑_j=1^dα_j,ix_j+c_j)}.

 where θ=(β_1,⋯,β_n,α_1,1,⋯,α_d,n,c_1,c_1,⋯, c_n)∈ℝ^2n+n(1+d) compose the elements of the parameter space.
Then we have the following result.


Let ℭ^n(ψ) be given by (<ref>) where ψ is assumed to be in 𝒞^2(ℝ^d), bounded and non-constant. Set ℭ(ψ)=⋃_n≥ 1ℭ^n(ψ). Assume that Ω_T is compact and consider the measures ν_1, ν_2, ν_3 whose support is contained in Ω_T, Ω and ∂Ω_T respectively. In addition, assume that the PDE (<ref>) has a unique classical solution such that (<ref>) holds. Also, assume that the nonlinear terms ∂α_i(t,x,u, p)/∂ p_j and γ̂(t,x,u,p) are 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. Then, for every ϵ>0, there exists a positive constant K>0 that may depend on 
sup_Ω_T|u|, sup_Ω_T|∇_xu| and sup_Ω_T|∇^(2)_xu|
such that there exists a function f∈ℭ(ψ) that satisfies

    J(f)≤ Kϵ.



The proof of this theorem is in the Appendix.





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

    ∂_t u(t,x)   -div(α(t,x,u(t,x), ∇ u(t,x)))+γ(t,x,u(t,x),∇ u(t,x))=0,  for (t,x)∈Ω_T
    
    u(0,x)   =u_0(x),  for x∈Ω
    
    u(t,x)   =0,  for (t,x)∈∂Ω_T,

as n →∞. 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,  for (t,x)∈∂Ω_T, i.e., g=0, to circumvent certain technical difficulties arising due to inhomogeneous boundary conditions. If g≠0 such that g is the trace of some appropriately smooth function, say ϕ, then one can reduce the inhomogeneous boundary conditions on ∂Ω_T to the homogeneous one by introducing in place of u the new function u-ϕ, see Section 4 of Chapter V in <cit.> or Chapter 8 of <cit.> 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.]

The objective function is

    J(f)   =𝒢[f]^2_2,Ω_T +   f^2_2, ∂Ω_T +  f(0, ·) - u_0 ^2_2, Ω






Recall that the norms above are L^2(X) norms in the respective space X=Ω_T, ∂Ω_T and Ω respectively. From Theorem <ref>, we have that

    J(f^n) → 0  as n→∞.


Each neural network f^n satisfies the PDE

    𝒢[f^n](t,x)   =h^n(t,x),  for (t,x)∈Ω_T
    
    f^n(0,x)   =u^n_0(x),  for x∈Ω
    
    f^n(t,x)   =g^n(t,x),  for (t,x)∈∂Ω_T

for some  h^n, u^n_0, and g^n  such that

    h^n^2_2,Ω_T+   g^n^2_2, ∂Ω_T +  u^n_0 - u_0 ^2_2, Ω→ 0  as n→∞.


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



  * There is a constant  μ>0 and positive functions κ(t,x),λ(t,x) such that for all (t,x)∈Ω_T we have

    α(t,x,u,p)≤μ(κ(t,x)+p),  and  |γ(t,x,u,p)|≤λ(t,x)p,

with κ∈ L^2(Ω_T), λ∈ L^d+2+η(Ω_T) for some η>0.



  * α(t,x,u,p) and γ(t,x,u,p) are Lipschitz continuous in (t,x,u,p)∈Ω_T×ℝ×ℝ^d uniformly on compacts of the form {(t,x)∈Ω̅_T, |u|≤ C, |p|≤ C}.

  * α(t,x,u,p) is differentiable with respect to (x,u,p) with continuous derivatives. 






  * There is a positive constant ν>0 such that

    α(t,x,u,p)p≥ν |p|^2

and

    <α(t,x,u,p_1)-α(t,x,u,p_2), p_1-p_2> >0,  for every   p_1,p_2∈ℝ^d, p_1≠ p_2.

   


  *  u_0(x)∈𝒞^0,2+ξ(Ω̅) for some ξ>0[In general, the Hölder space 𝒞^0,ξ(Ω̅) is the Banach space of continuous functions in Ω̅ having continuous derivatives up to order [ξ] in Ω̅ with finite corresponding uniform norms and finite uniform ξ-[ξ] Hölder norm. Analogously, we also define the Hölder space 𝒞^0,ξ,ξ/2(Ω̅_T) which in addition has finite [ξ]/2 and (ξ-[ξ])/2 regular and Hölder derivatives norms in time respectively. These spaces are denoted by H^ξ(Ω̅) and H^ξ,ξ/2(Ω̅_T) respectively in <cit.>. ] with itself and its first derivative bounded in Ω̅.

  * Ω is a bounded, open subset of ℝ^d with boundary ∂Ω∈𝒞^2.

  * For every n∈ℕ,  f^n∈𝒞^1,2(Ω̅_T). In addition, (f^n)_n∈ℕ∈ L^2(Ω_T).







Assume that Condition <ref> and (<ref>) hold. Then, problem (<ref>) has a  unique bounded solution in 
𝒞^0,δ,δ/2(Ω̅_T)∩ L^2(0,T; W^1,2_0(Ω))∩ W^(1,2),2_0(Ω_T^') for some δ>0 and any interior subdomain Ω_T^' of Ω_T[Here W^(1,2),2_0(Ω_T^') denotes the Banach space which is the closure of C_0^∞(Ω_T^') with elements from L^2(Ω_T^') having generalized derivatives of the form D^r_tD^s_x with r,s such that 2r+s≤ 2 with the usual Sobolev norm.]. 

 In addition,  f^n converges  to u, the unique solution to (<ref>), strongly in L^ρ(Ω_T) for every ρ<2. If, in addition, the sequence {f^n(t,x)}_n∈ℕ is uniformly bounded in n and equicontinuous then the convergence to u is uniform in Ω_T.


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 α(t,x,u,p) and γ(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 <cit.> for details.  As a matter of fact the results of Chapter V.6 in <cit.> show that with assuming a little bit more on the growth of the derivatives of the nonlinear functions α(t,x,u,p),γ(t,x,u,p) will lead to ∇_xu∈𝒞^0,δ',δ'/2(Ω_T) for some δ'>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.



The uniform, in n,  L^2 bound for the sequence {f^n}_n∈ℕ is easily satisfied for a bounded neural network approximation sequence f^n(t,x). However, we believe that it is true for a wider class of models, after all one expects that to be true if f^n indeed converges in L^ρ for ρ<2. The condition on equicontinuity for {f^n(t,x)}  allows to both simplify the proof and make a stronger claim as well. However, it is only a sufficient condition and not necessary. The paper, <cit.>, see Theorems 19 and 20 therein, discusses structural restrictions (a-priori boundedness and summability) that can be imposed on the unknown weights of feedforward neural networks, belonging in the class ℭ(ψ)=⋃_n≥ 1ℭ^n(ψ) as defined by (<ref>), which then guarantee both equicontinuity and universal approximation properties of the neural network for continuous and bounded functions. As it is also discussed in <cit.>, equicontinuity is also related to fault-tolerance properties of neural networks, a subject worthy of further study in the context of PDEs. However, we do not discuss this further here as this would be a topic for a different paper.



Let us present the case of linear parabolic PDEs  in Example <ref> below.
[Linear case]
Let us assume that the operator 𝒢 is linear in u and ∇ u. In particular, let us set

    α_i(t,x,u,p)=∑_j=1^n(σσ^T)_i,j(t,x)p_j, i=1,⋯ d

and

    γ(t,x,u,p)=-<b(t,x),p> +∑_i,j=1^d∂/∂ x_i(σσ^T)_i,j(t,x)p_j-c(t,x) u.


Assume that there are positive constants ν,μ>0 such that for every ξ∈ℝ^d the matrix [(σσ^T)_i,j(t,x)]_i,j=1^d satisfies

    ν |ξ|^2≤∑_i,j=1^d(σσ^T)_i,j(t,x)ξ_iξ_j≤μ |ξ|^2

and that the coefficients b and c are such that

    ∑_i=1^db_i^2_q,r,Ω_T+c_q,r,Ω_T≤μ,  for some μ>0

where we recall for example c_q,r,Ω_T=(∫_0^T(∫_Ω|c(t,x)|^qdx)^r/q)^1/r and  r,q satisfy the relations

    1/r+d/2q=1
    
    q∈(d/2,∞], r∈[1,∞),  for  d≥ 2,
    
    q∈[1,∞], r∈[1,2],  for  d=1.


In particular, the previous bounds always hold in the case of coefficients b and c that are bounded in Ω_T. Under these conditions, standard results for linear PDE's, see for instance Theorem 4.5 of Chapter III of <cit.> for a related result, show that approximation results analogous to that of Theorem <ref> 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 <ref>, <ref> and <ref> 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 <ref> 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.







§ PROOFS OF CONVERGENCE RESULTS

In this section we have gathered the proofs of the theoretical results of Section  <ref>.







By Theorem 3 of <cit.> we know that there is a function f∈ℭ(ψ) that is uniformly 2-dense on compacts of 𝒞^2(ℝ^1+d).  This means that for u∈𝒞^1,2([0,T]×ℝ^d) and ϵ>0, there is f∈ℭ(ψ) such that

    sup_(t,x)∈Ω_T|∂_tu(t,x)-∂_t f(t,x;θ)|+max_|a|≤ 2sup_(t,x)∈Ω̅_T |∂^(a)_xu(t,x)-∂^(a)_xf(t,x;θ)|<ϵ


We have assumed that (u,p)↦γ̂(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

    |γ̂(t,x,u,p)-γ̂(t,x,v,s)|   ≤(|u|^q_1/2+|p|^q_2/2+|v|^q_3/2+|s|^q_4/2)(|u-v|+|p-s|).

for some constants 0≤ q_1,q_2,q_3,q_4<∞. Therefore we obtain, using Hölder inequality with exponents r_1,r_2,

    ∫_Ω_T|γ̂(t,x,f,∇_xf)-γ̂(t,x,u,∇_xu)|^2dν_1(t,x)≤
       ≤∫_Ω_T(|f(t,x;θ)|^q_1+|∇_xf(t,x;θ)|^q_2+|u(t,x)|^q_3+|∇_xu(t,x)|^q_4) 
           ×(|f(t,x;θ)-u(t,x)|^2+|∇_xf(t,x;θ)-∇_xu(t,x)|^2)dν_1(t,x)
       ≤(∫_Ω_T(|f(t,x;θ)|^q_1+|∇_xf(t,x;θ)|^q_2+|u(t,x)|^q_3+|∇_xu(t,x)|^q_4)^r_1dν_1(t,x))^1/r_1
           ×(∫_Ω_T(|f(t,x;θ)-u(t,x)|^2+|∇_xf(t,x;θ)-∇_xu(t,x)|^2)^r_2dν_1(t,x))^1/r_2
       ≤ K(∫_Ω_T(|f(t,x;θ)-u(t,x)|^q_1+|∇_xf(t,x;θ)-∇_xu(t,x)|^q_2+|u(t,x)|^q_1∨ q_3+|∇_xu(t,x)|^q_2∨ q_4)^r_1dν_1(t,x))^1/r_1
           ×(∫_Ω_T(|f(t,x;θ)-u(t,x)|^2+|∇_xf(t,x;θ)-∇_xu(t,x)|^2)^r_2dν_1(t,x))^1/r_2
       ≤ K(ϵ^q_1+ϵ^q_2+sup_Ω_T|u|^q_1∨ q_3+sup_Ω_T|∇_xu|^q_2∨ q_4)ϵ^2

where the unimportant constant K<∞ may change from line to line and for two numbers q_1∨ q_3=max{q_1,q_3}. In the last step we used (<ref>).


In addition, we have also assumed that for every i,j∈{1,⋯ d}, the mapping (u,p)↦∂α_i(t,x,u, p)/∂ 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

    |∂α_i(t,x,u, p)/∂ p_j-∂α_i(t,x,v, s)/∂ s_j|   ≤(|u|^q_1/2+|p|^q_2/2+|v|^q_3/2+|s|^q_4/2)(|u-v|+|p-s|).

for some constants 0≤ q_1,q_2,q_3,q_4<∞. Denote for convenience

    ξ(t,x,u,∇ u, ∇^2u)=∑_i,j=1^d∂α_i(t,x,u(t,x), ∇ u(t,x))/∂ u_x_j∂_x_i,x_ju(t,x).

Then, similarly to (<ref>) we have after an application of Hölder inequality, for some constant K<∞ that may change from line to line,

    ∫_Ω_T|ξ(t,x,f,∇_xf, ∇^2_xf)-ξ(t,x,u,∇_xu,∇^2_xu)|^2dν_1(t,x)≤
       ≤∫_Ω_T|∑_i,j=1^d(∂α_i(t,x,f(t,x;θ), ∇ f(t,x;θ))/∂ f_x_j-∂α_i(t,x,u(t,x), ∇ u(t,x))/∂ u_x_j)∂_x_i,x_ju(t,x)|^2dν_1(t,x)
         +∫_Ω_T|∑_i,j=1^d∂α_i(t,x,f(t,x;θ), ∇ f(t,x;θ))/∂ f_x_j(∂_x_i,x_jf(t,x;θ)-∂_x_i,x_ju(t,x))|^2dν_1(t,x)
       ≤ K ∑_i,j=1^d(∫_Ω_T|∂_x_i,x_ju(t,x)|^2pdν_1(t,x))^1/p×
         ×(
    ∫_Ω_T|∂α_i(t,x,f(t,x;θ), ∇ f(t,x;θ))/∂ f_x_j-∂α_i(t,x,u(t,x), ∇ u(t,x))/∂ u_x_j|^2qdν_1(t,x))^1/q+
         + K ∑_i,j=1^d(∫_Ω_T|∂α_i(t,x,f, ∇ f)/∂ f_x_j|^2p dν_1(t,x))^1/p(∫_Ω_T|∂_x_i,x_jf(t,x;θ)-∂_x_i,x_ju(t,x)|^2qdν_1(t,x))^1/q


    ≤ K ∑_i,j=1^d(∫_Ω_T|∂_x_i,x_ju(t,x)|^2pdν_1(t,x))^1/p×
         ×(∫_Ω_T(|f(t,x;θ)-u(t,x)|^q_1+|∇_xf(t,x;θ)-∇_xu(t,x)|^q_2+|u(t,x)|^q_1∨ q_3+|∇_xu(t,x)|^q_2∨ q_4)^qr_1dν_1(t,x))^1/(qr_1)
         ×(∫_Ω_T(|f(t,x;θ)-u(t,x)|^2+|∇_xf(t,x;θ)-∇_xu(t,x)|^2)^qr_2dν_1(t,x))^1/(qr_2)
         + K ∑_i,j=1^d(∫_Ω_T|∂α_i(t,x,f, ∇ f)/∂ f_x_j|^2p dν_1(t,x))^1/p(∫_Ω_T|∂_x_i,x_jf(t,x;θ)-∂_x_i,x_ju(t,x)|^2qdν_1(t,x))^1/q
       ≤ Kϵ^2,

where in the last step we followed the computation in (<ref>) and used  (<ref>).


Using (<ref>) and (<ref>)-(<ref>) we subsequently obtain for the objective function (note that 𝒢[u](t,x)=0 for u that solves the PDE)
 
    J(f)   =𝒢[f](t,x)_Ω_T, ν_1^2 +     f(t,x; θ ) - g(t,x) _∂Ω_T, ν_2^2 +  f(0, x; θ) - u_0(x) _Ω, ν_3^2
       =𝒢[f](t,x)-𝒢[u](t,x) _Ω_T, ν_1^2 +     f(t,x; θ ) - g(t,x) _∂Ω_T, ν_2^2 +  f(0, x; θ) - u_0(x) _Ω, ν_3^2
       ≤∫_Ω_T|∂_tu(t,x)-∂_t f(t,x;θ)|^2dν_1(t,x)+
    ∫_Ω_T|ξ(t,x,f,∇ f, ∇^2f)-ξ(t,x,u,∇ u, ∇^2u)|^2dν_1(t,x)
       +∫_Ω_T|γ(t,x,f,∇_xf)-γ(t,x,u,∇_xu)|^2dν_1(t,x)+
     ∫_∂Ω_T|f(t,x;θ)-u(t,x)|^2dν_2(t,x)
     +
           +∫_Ω|f(0,x;θ)-u(0,x)|^2dν_3(t,x)
       ≤ Kϵ^2

for an appropriate constant K<∞. The last step completes the proof of the Theorem after rescaling ϵ.





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

Let us first consider problem (<ref>) with g^n(t,x)=0 and let us denote the solution to this problem by f̂^n(t,x). 
Due to Condition <ref>, Lemma 4.1 of <cit.>  applies and gives that {f̂^n}_n∈ℕ is uniformly bounded with respect to n in at least L^∞(0,T;L^2(Ω))∩ L^2(0,T; W^1,2_0(Ω)) (in regards to such uniform energy bound results we also refer the reader to Theorem 2.1 and Remark 2.14 of <cit.> for the case γ=0 and to <cit.>  for related results in more general cases). As a  matter of fact f̂^n is more regular than stated, see Section 6, Chapter V of <cit.>, but we will not make use of this fact in the convergence proof of f̂^n to u. These uniform energy bounds imply that we can extract a subsequence, denoted also by {f̂^n}_n∈ℕ, which converges to some u in the weak-* sense in L^∞(0,T;L^2(Ω)) and weakly in L^2(0,T; W^1,2_0(Ω)) and to some v weakly in L^2(Ω) for every fixed t∈(0,T]. 

Next let us set q=1+d/d+4∈ (1,2) and note that for conjugates, r_1,r_2>1 such that 1/r_1+1/r_2=1

    ∫_Ω_T|γ(t,x,f̂^n,∇_xf̂^n)|^qdtdx    ≤∫_Ω_T|λ(t,x)|^q|∇_xf̂^n(t,x)|^qdtdx
       ≤(∫_Ω_T|λ(t,x)|^r_1qdtdx)^1/r_1(∫_Ω_T|∇_xf̂^n(t,x)|^r_2q dtdx)^1/r_2.



Let us choose r_2=2/q>1. Then we calculate r_1=r_2/r_2-1=2/2-q. Hence, we have that r_1q=d+2. Recalling the assumption λ∈ L^d+2(Ω_T) and the uniform bound on the ∇_xf̂^n_2 we subsequently obtain that for q=1+d/d+4, there is a constant C<∞ such that

    ∫_Ω_T|γ(t,x,f̂^n,∇_xf̂^n)|^qdtdx    ≤ C.


The latter estimate together with the  growth assumptions on α(·)  from Condition <ref>, imply that {∂_tf̂^n}_n∈ℕ is bounded uniformly with respect to n in L^1+d/(d+4)(Ω_T) and in L^2(0,T; W^-1,2(Ω)). Consider the conjugates 1/δ_1+1/δ_2=1  with δ_2> max{2,d}. Due to the embedding

    W^-1,2(Ω)⊂ W^-1,δ_1(Ω),   L^q(Ω)⊂ W^-1,δ_1(Ω),  and  L^2(Ω)⊂ W^-1,δ_1(Ω),

we have that {∂_tf̂^n}_n∈ℕ is bounded uniformly with respect to n in L^1(0,T; W^-1,δ_1(Ω)). Define now the spaces X=W^1,2_0(Ω), B=L^2(Ω) and Y=W^-1,δ_1(Ω), and notice that

    X⊂ B⊂ Y

with the first embedding being compact.  Then, Corollary 4 of <cit.> yields relative compactness of {f̂^n}_n∈ℕ in L^2(Ω_T), which means that {f̂^n}_n∈ℕ converges strongly to u in that space. Thus, up to subsequences, {f̂^n}_n∈ℕ converges almost everywhere to u in Ω_T.



The nonlinearity of the α and γ functions with respect to the gradient prohibits us from passing to the limit directly in the respective weak formulation. However, the uniform boundedness of {f̂_n}_n∈ℕ in L^σ(0,T; W^1,σ_0(Ω)) with σ>1 (in fact here σ=2) and its weak convergence to u in that space,  allows us to conclude, as in Theorem 3.3 of <cit.>,  that

    ∇f̂^n→∇ u  almost everywhere in Ω_T.


Hence, we obtain that {f̂^n}_n∈ℕ converges  to u strongly also in L^ρ(0,T; W^1,ρ_0(Ω)) for every ρ<2. 

In preparation to passing to the limit as n→∞ in the weak formulation, we need to study the behavior of the nonlinear terms. Recalling the assumptions on α(t,x,u,p) we have for ρ<2 and for a measurable set A⊂Ω_T (the constant K<∞ may change from line to line)






    ∫_A|α(t,x,f̂^n,∇f̂^n)|^ρdtdx≤ K [∫_A|κ(t,x)|^ρdtdx+∫_A|∇f̂^n(t,x)|^ρdtdx]
           ≤ K [∫_A|κ(t,x)|^ρdtdx+(∫_Ω_T|∇f̂^n(t,x)|^2dtdx)^ρ/2|A|^1-ρ/2]
           ≤ K [∫_A|κ(t,x)|^ρdtdx+|A|^1-ρ/2].

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

    α(t,x,f̂^n,∇f̂^n)→α(t,x,u,∇ u)  strongly in  L^ρ(Ω_T)

as n→∞, for every 1<ρ<2. For the same reason, an analogous estimate to (<ref>), gives

    ∫_A|γ(t,x,f̂^n,∇_xf̂^n)|^qdtdx ≤  K (∫_A|λ(t,x)|^d+2dtdx)^(2-q)/2≤ K |A|^η/d+2+η

 implying, via Vitali's theorem,  that

    γ(t,x,f̂^n,∇f̂^n)→γ(t,x,u,∇ u)  strongly in  L^q(Ω_T)

as n→∞, for q=1+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(Ω). The weak formulation of the PDE (<ref>) with g^n=0 reads as follows. For every t_1∈(0,T]

    ∫_Ω_t_1[-f̂^n∂_tϕ +<α(t,x,f̂^n,∇f̂^n), ∇ϕ>+(γ(t,x,f̂^n,∇f̂^n)-h^n)ϕ](t,x)dxdt
          +∫_Ωf̂^n(t_1,x)ϕ(t_1,x)dx-∫_Ωu^n_0(x)ϕ(0,x)dx=0

for every ϕ∈ C^∞_0(Ω_T). Using the above convergence results, we then obtain that the limit point u satisfies for every t_1∈(0,T] the equation

    ∫_Ω_t_1[-u∂_tϕ +<α(t,x,u,∇ u), ∇ϕ>+γ(t,x,u,∇ u)ϕ](t,x)dxdt
          +∫_Ωu(t_1,x)ϕ(t_1,x)dx-∫_Ωu_0(x)ϕ(0,x)dx=0,

which is the weak formulation of the equation (<ref>).

It remains to discuss the convergence of f^n-f̂^n to zero, where we recall that f^n is the neural network approximation satisfying (<ref>) and f̂^n satisfies (<ref>) with g^n=0. The functions f^n∈𝒞^1,2(Ω̅_T) and Ω̅_T is compact. We have also assumed that {f^n}_n is uniformly bounded in L^2(Ω_T). This implies that, up to a subsequence, f^n will converge at least weakly in L^2(Ω_T). Moreover,  the boundary values g^n(t,x) (which is nothing else by f^n(t,x) evaluated at the smooth boundary ∂Ω_T) in (<ref>)  converge to zero strongly in L^2. It is then a standard result that  g^n, i.e., f^n evaluated at the boundary, converges to zero, at least, almost uniformly along a subsequence, see for example Lemma 2.1 in Chapter II of <cit.>. As it then follows, for example, by the proof of Theorems 6.3-6.4-6.5  of Chapter V.6 in <cit.>, using smoothness and uniqueness,  f^n will differ from the solution to the PDE (<ref>) with g^n=0,  f̂^n(t,x), by a negligible amount as n→∞ in the almost everywhere sense. The assumed uniform L^2 bound for {f^n}_n∈ℕ together with the previously derived uniform L^2(Ω_T) bound for {f̂^n}_n∈ℕ yield uniform L^2(Ω_T) boundedness for {f^n-f̂^n}_n∈ℕ. Then, by Vitali's theorem again, we get that {f^n-f̂^n}_n∈ℕ goes to zero strongly in L^ρ(Ω_T) for every ρ<2. 


The previously derived strong convergence of {f^n-f̂^n}_n∈ℕ  to zero  in L^ρ(Ω_T) for every ρ<2, together with the strong L^2(Ω_T) convergence of {f̂^n}_n∈ℕ  to u, conclude the proof of the convergence in L^ρ(Ω_T) for every ρ<2 using triangle inequality.


If {f^n(t,x)} is equicontinuous then, Lemma 3.2 of <cit.> gives uniform convergence of g^n to zero. Hence, by the previous analysis, it will certainly be true that  {f^n}_n∈ℕ  converges to u in L^ρ(Ω_T) for every ρ<2. The L^ρ convergence to zero together with boundedness and equicontinuity of the sequence {f^n(t,x)} results then in uniform convergence due to the well known Arzelà-Ascoli theorem. 





99






Glynn S. Asmussen and P. Glynn, Stochastic Simulation: Algorithms and Analysis, Springer, 2007.


Jentzen  C. Beck, W. E., A. Jentzen.  Machine learning approximation algorithms for high-dimensional fully nonlinear partial differential equations and second-order backward stochastic differential equations, arXiv:1709.05963, 2017.




BertsekasThitsiklis2000 D. Bertsekas and J. Tsitsiklis, Gradient convergence in gradient methods via errors, SIAM Journal of Optimization, Vol.10, No. 3, pgs. 627-642, 2000.


BDGO1997 L. Boccardo, A. Dall`Aglio, T. Gallouët and L. Orsina, Nonlinear parabolic equations with measure data, Journal of Functional Analysis, Vol. 147, pp. 237-258, (1997).

BoccardoPorzioPrimo L. Boccardo, M.M. Porzio and A. Primo, Summability and existence results for nonlinear parabolic equations, Nonlinear Analysis: Theory, Methods and Applications, Vol. 71, Issue 304, pp. 1-15, 2009.

Bungartz2 H. Bungartz, A. Heinecke, D. Pfluger, and S. Schraufstetter, Option pricing with a direct adaptive sparse grid approach, Journal of Computational and Applied Mathematics, Vol. 236, No. 15, pgs. 3741-3750, 2012.

Griebel H. Bungartz and M Griebel, Sparse Grids, Acta numerica, Vol. 13, pgs. 174-269, 2004.

ChandraSingh2004
P. Chandra and Y. Singh, Feedforward sigmoidal networks - equicontinuity and fault-tolerance properties,  IEEE Transactions on Neural Networks,   Vol. 15, Issue 6,  pp. 1350-1366, 2004.

Chaudhari  P. Chaudhari, A. Oberman, S. Osher, S. Soatto, and G. Carlier.  Deep relaxation: partial differential equations for optimizing deep neural networks, 2017.

Cerrai S. Cerrai, Stationary Hamilton-Jacobi Equations in Hilbert Spaces and Applications to a Stochastic Optimal Control Problem, SIAM Journal on Control and Optimization, Vol. 40, No. 3, pgs. 824-852, 2001.

Cybenko G. Cybenko, Approximation by superposition of a sigmoidal function, Mathematics of Control, Signals and Systems, Vol. 2, pp. 303-314, 1989.


Gaines  A. Davie and J. Gaines, Convergence of Numerical Schemes for the Solution of Parabolic Stochastic Partial Differential Equations, Mathematics of Computation, Vol. 70, No. 233, pgs.  121-134, 2000.



Dean J. Dean, G. Corrado, R. Monga, K. Chen, M. Devin, M. Mao, A. Senior, P. Tucker, K. Yang, Q. Le, and A. Ng. Large scale distributed deep networks. In Advances in Neural Information Processing Systems, pp. 1223-1231, 2012.



Debussche A. Debussche, M. Fuhrman, and G. Tessitore, Optimal control of a stochastic heat equation with boundary-noise and boundary-control, ESAIM: Control, Optimisation and Calculus of Variations, Vol. 13, No. 1, pgs. 178-205, 2007.




Weinan  W. E., J. Han, and A. Jentzen.  Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations, Communications in Mathematics and Statistics, Springer, 2017.


Masaaki  M. Fujii, A. Takahashi, and M. Takahashi.  Asymptotic Expansion as Prior Knowledge in Deep Learning Method for high dimensional BSDEs, arXiv:1710.07030, 2017.

Garcia1970 A. M. Garcia, E. Rodemich, H. Rumsey Jr. and M. Rosenblatt,  A real variable lemma and the continuity of paths of some Gaussian processes, Indiana University Mathematics Journal,
Vol. 20, No. 6, pp. 565-578, December, 1970.

Xavier  X. Glorot and Y. Bengio.  Understanding the difficulty of training deep feedforward neural networks. In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, pp. 249-256, 2010.

GainesLectureNotes J. Gaines. Numerical experiments with SPDEs. London Mathematical Society Lecture Note Series, 55-71, 1995.


GilbargTrudinger
D. Gilbarg and N.S. Trudinger. Elliptic partial differential equations of second order, second edition. Springer-Verlang, Berlin Heidelberg, 1983.

Gyongy I. Gyöngy,  Lattice Approximations for Stochastic Quasi-Linear Parabolic Partial Differential Equations Driven by Space-Time White Noise I, Potential Analysis, Vol. 9, Issue 1, pgs. 1–25, 1998.



Bungartz1 A. Heinecke, S. Schraufstetter, and H. Bungartz, A highly parallel Black-Scholes solver based on adaptive sparse grids, International Journal of Computer Mathematics, Vol. 89, No. 9, pgs. 1212-1238, 2012.



Haugh M. Haugh and L. Kogan, Pricing American Options: A Duality Approach, Operations Research, Vol. 52, No. 2, pgs. 258-270, 2004.

SchmidhuberLSTM S. Hochreiter and J. Schmidhuber, Long short-term memory, Neural Computation, Vol. 9, No. 8, pgs. 1735-1780, 1997.



HornikEtAl K. Hornik, M. Stinchcombe and H. White, Universal approximation of an unknown mapping and its derivatives using multilayer feedforward networks, Neural Networks, Vol. 3, Issue 5, pp. 551-560, 1990.

Hornik K. Hornik, Approximation capabilities of multilayer feedforward networks, Neural Networks, Vol. 4, pp. 251-257, 1991.

ADAM D. Kingma and J. Ba, ADAM: A method for stochastic optimization. arXiv: 1412.6980, 2014.

Ladyzhenskaya  O. A. Ladyzenskaja, V. A. Solonnikov and N. N. Ural'ceva, Linear and Quasi-linear Equations of Parabolic Type (Translations of Mathematical Monographs Reprint), American Mathematical Society, Vol. 23, 1988.

Lagaris  I. Lagaris, A. Likas, and D. Fotiadis, Artificial neural networks for solving ordinary and partial differential equations, IEEE Transactions on Neural Networks, Vol. 9, No. 5, pgs. 987-1000, 1998.

Likas I. Lagaris, A. Likas, and D. Papageorgiou, Neural-network methods for boundary value problems with irregular boundaries, IEEE Transactions on Neural Networks, Vol. 11, No. 5, pgs. 1041-1049, 2000.

Lee  H. Lee, Neural Algorithm for Solving Differential Equations, Journal of Computational Physics, Vol. 91, pgs. 110-131, 1990.

Ling J. Ling, A. Kurzawski, and J. Templeton. Reynolds averaged turbulence modelling using deep neural networks with embedded invariance. Journal of Fluid Mechanics, Vol. 807, pgs. 155-166, 2016.


LongstaffSchwartz F. Longstaff and E. Schwartz, Valuing American Options by Simulation: A Simple Least-Squares Approach, Review of Financial Studies, Vol. 14, pgs.  113-147, 2001.


Magliocca M. Magliocca,  Existence results for a Cauchy-Dirichlet parabolic problem with a repulsive gradient term, Nonlinear Analysis, Vol. 166, pp. 102-143, 2018


Malek  A. Malek and R. Beidokhti, Numerical solution for high order differential equations using a hybrid neural network-optimization method, Applied Mathematics and Computation, Vol. 183, No. 1, pgs.  260-271, 2006.


Masiero F. Masiero, HJB equations in infinite dimensions, Journal of Evolution Equations, Vol. 16, No. 4, pgs. 789-824, 2016.



NardoEtAll R. Di Nardo, F. Feo, O. Guibé, Existence result for nonlinear parabolic equations with lower order terms, Anal. Appl.(Singap.),  Vol. 09, No. 02,  pp. 161-186, 2011.



Petersen2017
P. Petersen and  F. Voigtlaender, Optimal approximation of piecewise smooth functions using deep ReLU neural networks, (2017), arXiv: 1709.05289v4.








Pinkus
A. Pinkus. Approximation theory of the MLP model in neural networks. Acta Numerica, Vol. 8, pp. 143–195, 1999.






Porzio M. M. Porzio, Existence of solutions for some “noncoercive" parabolic equations, Discrete and Continuous Dynamical Systems, Vol. 5, Issue 3,  pp. 553-568, 1999.


Karniadakis1 M. Raissi, P. Perdikaris, and G. Karniadakis. Physics Informed Deep Learning (Part I): Data-driven Solutions of Nonlinear Partial Differential Equations, arXiv:1711.10561, 2017.

Karniadakis2 M. Raissi, P. Perdikaris, and G. Karniadakis. Physics Informed Deep Learning (Part II): Data-driven Discovery of Nonlinear Partial Differential Equations, arXiv:1711.10566, 2017.

Reisinger C. Reisinger and G. Wittum, Efficient hierarchical approximation of high-dimensional option pricing problems, SIAM Journal on Scientific Computing, Vol. 29, No. 1, pgs. 440-458, 2007.

Reisinger2 C. Reisinger, Analysis of linear difference schemes in the sparse grid combination technique, IMA Journal of Numerical Analysis, Vol. 33, No. 2, pgs. 544-581, 2012.



ChrisRogers L.C.G. Rogers, Monte-Carlo Valuation of American Options, Mathematical Finance, Vol. 12, No. 3, pgs. 271-286, 2002.

Rudd  K. Rudd, Solving Partial Differential Equations using Artificial Neural Networks, PhD Thesis, Duke University, 2013.

SchmidhuberHighway R. Srivastava, K. Greff, and J. Schmidhuber, Training very deep networks, In Advances in Neural Information Processing Systems, pp. 2377-2385, 2015.


Simon J. Simon, Compact sets in the space L^p(0,T;B), Annali di Matematica Pura ed Applicata, Vol. 146, pp. 65-96, 1987.

Tompson  J. Tompson, K. Schlachter, P. Sprechmann, and K. Perlin.  Accelerating Eulerian Fluid Simulation with Convolutional Networks, Proceedings of Machine Learning Research, Vol. 70, pgs. 3424-3433, 2017.




