Solve High-Dimensional Reﬂected Partial Differential Equations by Neural Network Method

: Reﬂected partial differential equations (PDEs) have important applications in ﬁnancial mathematics, stochastic control, physics, and engineering. This paper aims to present a numerical method for solving high-dimensional reﬂected PDEs. In fact, overcoming the “dimensional curse” and approximating the reﬂection term are challenges. Some numerical algorithms based on neural networks developed recently fail in solving high-dimensional reﬂected PDEs. To solve these problems, ﬁrstly, the reﬂected PDEs are transformed into reﬂected backward stochastic differential equations (BSDEs) using the reﬂected Feyman–Kac formula. Secondly, the reﬂection term of the reﬂected BSDEs is approximated using the penalization method. Next, the BSDEs are discretized using a strategy that combines Euler and Crank–Nicolson schemes. Finally, a deep neural network model is employed to simulate the solution of the BSDEs. The effectiveness of the proposed method is tested by two numerical experiments, and the model shows high stability and accuracy in solving reﬂected PDEs of up to 100 dimensions.


Introduction
As a significant mathematical tool for characterizing singularities, the reflection theory and reflected PDEs have a broad range of applications in various fields, including physics, engineering, and finance. These applications encompass describing materials with memory and heredity, determining system equilibrium points, and pricing financial assets. Solving reflected PDEs can be considered as a Skorohod problem [1]. The reflection term in these equations restricts the solution above an obstacle and represents the minimum external force that prevents the solution from breaching the obstacle [2,3]. Since it is challenging to find an explicit solution to reflected PDEs, the development of efficient numerical algorithms is a crucial area of research.
There are many numerical methods for solving PDEs without reflection. For lowdimensional PDEs, practical algorithms include the neural network-based method [4], finite element method [5,6], finite difference method [7][8][9], Monte Carlo approximation method [10,11], spectral Galerkin approximation method [12,13], and sparse grid approximation method [14,15]. These methods have been proven effective and stable in a large body of literature. However, when it comes to high-dimensional problems, these traditional methods are often not applicable due to the "dimensional curse" [16], which causes an exponential increase in computational complexity and rapid reduction in stability and efficiency as dimension increases. Although the Monte Carlo method can overcome the "dimensional curse", it can only approximate the solution of an isolated point in space. Recently, with the development of deep learning technology, more and more deep learning-based numerical methods have been developed for solving PDEs, successfully overcoming the "dimensional

Approximating Schemes for Reflected PDEs
In this section, we propose a numerical algorithm named the Deep C-N method for solving the obstacle problems for high-dimensional non-linear parabolic PDEs. Generally, an obstacle problem for PDEs can be totally solved via three steps based on the above method. In the first step, a connection will be built between as the solution for obstacle problems for PDEs and that of the corresponding reflected BSDEs by the non-linear reflected Feynman-Kac formula. As a result, the original problems for PDEs are converted to problems for RBSDEs. In the second step, the problems for RBSDEs will be transformed further. In this step, the issue about RBSDEs will be considered as an optimal stopping-time problem via the penalty approach. Consequently, our goal is settling a stochastic control problem. In the final step, the optimal stopping-time problem will be regarded as a deep learning issue based on the neural network model. Specifically, the neural network model will act on the policy function Z, which represents the gradient of the solution for the optimal stopping-time problem.

Nonlinear Parabolic Reflected PDEs
The PDEs with reflection can be presented by nonlinear parabolic PDEs with minimum constraints. The so-called reflection term forces the solution of the equation to be non-negative, and the reflected part is the minimum power that prevents the solution from leaving the non-negative interval. That is, the reflected PDEs are essentially obstacle problems, and the solutions have only two states: located over the obstacle or on the obstacle. In the beginning, we have some assumptions for the basic parameters and functions which are involved in this problem. Let T ∈ (0, ∞), d ∈ N, f(t, x, y, z) is a given non-linear function defined by f : is a given continuous function defined by g : R d → R . In this paper, we describe the continuous obstacle by function h(t, x), and the reflected problems for PDEs can be presented by the following scheme: The terminal value of the solution for Equation (1) satisfies u(T, x) = g(x), x ∈ R d . For Equation (1), there are three further hypotheses as follows: Hypothesis 2 (H2). h : [0, T] × R d → R is a continuous function related to parameters t and x. Meanwhile, it satisfies:

Hypothesis 3 (H3).
Remarks. C((·)) > 0 are constants from (H1) to (H3). ) It has been proven that the solution u(t,x) for Equation (1) is existing and unique. Furthermore, it is a function with up to polynomial growth.

From Reflected PDEs to Related Reflected BSDEs
In this subsection, we build a connection between the solution of reflected PDEs and reflected BSDEs via the nonlinear reflected Feynman-Kac formula.
Let (Ω, F , P) be a given probability space, B : [t, T] × Ω → R d be a d-dimensional standard Brownian motion in this space, F be a normal filtration set generated by B, and continuous functions b(t, x) and σ(t, x) satisfy the hypothesis (H1) in Section 2.1. Now, consider a stochastic process with d-dimensional, {X s } : [t, T] × Ω → R d satisfies: From [1], under the assumptions (H1)-(H3), the reflected BSDE (3) has a unique solution represented by a triple (Y t , Z t , K t ) [1].
Moreover, the solution of Equation (3) satisfies the following three properties:

Proposition 2.
(Section 8 in [2]) Under the assumptions (H1)-(H3), a special relationship between reflected PDEs and reflected BSDEs is established by a classic formula, namely, the Feyman-Kac formula: Consequently, from the above discussions, we transform the reflected problems for PDEs to related reflected problems for BSDEs. That is, u(t, x) can be simulated by Y t . We approximate the reflection term K T − K t via the penalization method [2]: where 0 ≤ t ≤ T, m → +∞ . Therefore, the reflected BSDE (3) can be approximated by the following scheme: Obviously, Equation (5) is a general BSDE without any reflection. We use a binary group (Y m t , Z m t ) to illustrate the solution of Equation (5). Then, Equation (1) can be solved approximately by Equation (5). That is, the solution of reflected PDEs can be presented effectively by the solution of related BSDEs. In this work, we simulate the gradient of the solution u(t, x) via a deep neural network rather than solving ∇u(t, x) directly. The final numerical approximations can be calculated by a series of explicit iteration formats. Details of the whole process will be described in Section 3.

Discretizing via Two Approaches
From Proposition 2, the goal of approximating u(t, x) can be transformed into approximations for Y t . Since the problem we are interested in is the initial value of u(t, x), our goal now is simulating Y 0 by the deep learning algorithm based on the deep neural network. To achieve this goal, we discretize Equations (2) and (5) in time dimension by the Euler and C-N schemes, separately. Firstly, we show some general instructions about these two schemes: the time span is a finite interval from 0 to T, 0 = t 0 < t 1 < · · · < t N = T, with Since this discrete format is implicit, i.e., the value at the latter moment cannot be described by a deterministic formula for the value at the previous moment, then next, we perform a second iteration to obtain the explicit format at any time interval [t n , t n+1 ], n = 0, 1, · · · , N − 1. Let s = 0, 1, 2, · · · , ε = 10 −6 , for any time interval [t n , t n+1 ], we have: It should be noted that the condition for stopping the iteration is

Numerical Experiments
There are three subsections in this part. Firstly, we introduce a general numerical algorithm for solving high-dimensional nonlinear reflected PDEs, then we illustrate the validity and reliability of this method proposed by us with two examples in Sections 3.2 and 3.3. We operate the equations from continuous time to separation time using the two discrete formats mentioned above, respectively, and the comparison of the significant results will be described by means of tables. The numerical experiments we present are performed in Python code using TensorFlow on a Lenovo Pro-13 with a Radeon microprocessor.

Deep C-N Algorithm for Solving High-Dimensional Nonlinear Reflected PDEs
In this subsection, we propose a general deep learning framework named the Deep C-N algorithm to calculate high-dimensional nonlinear PDEs with reflection. The issue we are interested in is the initial value of the solutions for Equation (1), i.e., u(0, x). From what we have discussed above, with the intermediary roles of reflected BSDEs, a deterministic relationship is established between the solutions for reflected PDEs and BSDEs. That is: In the whole process, the key step is approximating the gradient term of the solution for reflected PDEs, i.e., (∇uσ)(t n , x), using multilayer feedback deep neural networks. First of all, we define the definitions of some parameters in this neural networks-based algorithm. Let ρ ∈ N denote the dimension of reflected PDEs, θ ∈ R ρ is a set of parameters which need to be learned by neural networks, X θ t n , Y θ t n , Z θ t n , and n ∈ {0, 1, · · · , N} is the solution for BSDEs at time t n . For convenience, we use which is a part of the initial value of the solution for BSDEs. Then, we have: Further, the following relationships hold that: The loss function is defined by the main square error between the terminal output value calculated via the neural network and the real value at that moment. That is, the loss value is given by The parameters in neural networks are updating until the loss function is stable. In addition, the existence and uniqueness of the global minimum for Equation (15) is proven by [6]. We choose the stochastic gradient descent (SGD) approach as the optimization algorithm for calculating the loss value. In the back-propagation simulating process, the Adam optimizer is employed to update the parameters layer by layer.
Next, based on the deep neural network technique, the general algorithm framework for solving high-dimensional nonlinear PDEs with reflection is presented (see Figure 1). By dividing the time interval [0, T] by N, a total of N − 1 sub-neural networks are computed separately. Assuming that each sub-neural network has H hidden layers, all parameters to be learned will be optimized in each hidden layer. In particular, the results of the operations located in the hidden layer are batch normalized before they are passed through the activation function to the next layer.

3.
X , u t , X → h → h → (subnetworks) ⋯ → ∇u t , X is the key step in the whole calculating procedure. Our goal in this step is approximating the spatial gradients, and meanwhile, the weights θ are optimized in the (N − 1) sub-networks. 4. u t , X , ∇u t , X , K − K , B − B → u t , X is a forward iteration procedure that yields the neural network s final output as the unique approximation of u(T, X ), totally characterized by approximating scheme (14).

Allen-Cahn Equation
In this subsection, we consider the solution of PDEs without any boundary. For comparing the results with the Deep BSDE [22] approach, we choose the Allen-Cahn equation as our experimental case. The significant results are presented by figures and tables. The parameters and format of Allen-Cahn equations are as follows: ∀t ∈ [0, T], x, ω ∈ R , y ∈ R, z ∈ R × , suppose that d = 100, N = 20, T = 0.3, μ(t, x) = 0, σ(t, x) = √2ω, X = x + √2ω, f(t, x, y, z) = y − y is the nonlinear part of the equation, g(x) = 2 + ‖x‖ is the value function of the solution at the terminal moment. u(t, x) satisfies: Regarding the settings of the neural network, for comparison with the numerical results in [22], the hyperparameters and the selection of the optimization method are kept the same strategy as that in [22][23][24], except for the change in the discrete format and the Figure 1. Illustration of the neural network framework for solving obstacle problems for PDEs. There are N − 1 sub-networks in total and H hidden layers in each sub-network. Therefore, there exist (H + 1)(N − 1) layers in the whole network with parameters that need to be optimized. We divide the time internal [0, T] for intervals and each column for t = t 1 , t 2 , · · · , t N−1 .
The entire deep neural network-based algorithm consists of four types of operations: 1.
X t n , B t n+1 , B t n → X t n+1 (n = 1, 2, · · · N − 1 and the same settings in 2 to 4) is a forward iterative procedure, which is determined by approximating scheme (6); this procedure does not contain any parameters that need to be optimized.

2.
(K t n , u(t n , X t n ), X t n ) → K t n+1 is a forward iterative procedure too, which is characterized by approximating scheme (7). As in the previous step, no parameters need to be optimized in this operation. 3. X t n , u t n−1 , X t n−1 → h 1 n → h 2 n → (subnetworks) · · · → ∇u(t n , X t n ) is the key step in the whole calculating procedure. Our goal in this step is approximating the spatial gradients, and meanwhile, the weights θ n are optimized in the (N − 1) sub-networks. 4. u(t n , X t n ), ∇u(t n , X t n ), K t n+1 − K t n , B t n+1 − B t n → u t n+1 , X t n+1 is a forward iteration procedure that yields the neural network's final output as the unique approximation of u(T, X T ), totally characterized by approximating scheme (14).

Allen-Cahn Equation
In this subsection, we consider the solution of PDEs without any boundary. For comparing the results with the Deep BSDE [22] approach, we choose the Allen-Cahn equation as our experimental case. The significant results are presented by figures and tables. The parameters and format of Allen-Cahn equations are as follows: ∀t ∈ [0, T], x, ω ∈ R d , y ∈ R, z ∈ R 1×d , suppose that d = 100, N = 20, x, y, z) = y − y 3 is the nonlinear part of is the value function of the solution at the terminal moment. u(t, x) satisfies: Regarding the settings of the neural network, for comparison with the numerical results in [22], the hyperparameters and the selection of the optimization method are kept the same strategy as that in [22][23][24], except for the change in the discrete format and the input information of the neural network. The Relu function is chosen as the activation function, r = 5 × 10 −4 is the learning rate, and the optimization method uses the Adam optimizer and batch normalization, with a total of N − 1 subneural networks for stacking, each containing two hidden layers and each hidden layer contains d + 10 neurons.
The numerical results of Equation (16) are illustrated by two tables. Specifically, the results calculated by the Deep BSDE method are described in Table 1, and the results calculated by the Deep C-N method are described in Table 2. When solving Equation (16) using the Deep C-N algorithm, we tried appropriate values of iteration steps repeatedly in the interval [2000, 15,000], and found that 10,000 is a suitable number. Too-small values will make the loss value unstable, and too-large values will encounter an overfitting problem. The numerical experiment results shown in the table are based on the mean values of five independent tests. The true value of Equation (16) is 0.052802, which is derived from the Branching diffusion algorithm in [22]. The two tables above illustrate that the numerical results of the Deep C-N algorithm are better than the Deep BSDE algorithm. The most significant aspect is that the computational accuracy of the Deep C-N algorithm is better than that of the Deep BSDE algorithm. This is demonstrated by the fact that the loss value of the Deep BSDE algorithm is approaching 2 × 10 −4 , while the value of the Deep C-N algorithm is approaching 3 × 10 −5 , which means that the accuracy is improved by nearly an order of magnitude. In other aspects, the standard deviation of Y 0 of the Deep C-N method is smaller than the Deep BSDE method, which demonstrates that the approximating stability of the Deep C-N method is superior to the Deep BSDE method. In addition, the smaller relative L 1 -approximate error, to a certain extent, reflects the strengths of the model.
We now discuss the stability and accuracy of these two numerical models via eight equations with 50 to 400 dimensions. In Figures 2 and 3, it can be clearly seen that the numerical results of the Deep C-N algorithm are more compact. In particular, the loss function curves of the Deep C-N algorithm are close to overlapping when the dimensionality is greater than 100, while the Deep BSDE algorithm shows a larger difference. This indicates that the Deep C-N model has better stability than the Deep BSDE model when dealing with high-dimensional problems. Regarding the approximation accuracy of the model, Table 3 shows that in tests for different dimensions, the loss values of the Deep C-N model are smaller than the corresponding results of the Deep BSDE model. This indicates that the computational accuracy of the Deep C-N algorithm is higher in the interval of 50 to 400 dimensions. model, Table 3 shows that in tests for different dimensions, the loss values of th N model are smaller than the corresponding results of the Deep BSDE model cates that the computational accuracy of the Deep C-N algorithm is higher in t of 50 to 400 dimensions.

American Options
In this subsection, we consider the solution of reflected PDEs with continuo ary. We select the case of pricing American options in order to make a compa the test results in [4]. The significant results are presented as well. The param format of the American options are as follows: ∀t ∈ [0, T], x, ω ∈ R , y ∈ R, z ∈ R × , to be consistent with [4], suppose 5, 10, 20, 40, N = 20 , T = 1, K = 1, r = 0.05, is the risk-free interest rate, σ( dX = μ(t, x)X dt + σ(t, x)X dW , where W is a d-dimension Brownian motion mensions, g(x) = K − ∏ X is the payoff of option. At time t, the value of can option u(t, x) satisfies: N model are smaller than the corresponding results of the Deep BSDE model cates that the computational accuracy of the Deep C-N algorithm is higher in of 50 to 400 dimensions.

American Options
In this subsection, we consider the solution of reflected PDEs with continu ary. We select the case of pricing American options in order to make a compa the test results in [4]. The significant results are presented as well. The para format of the American options are as follows: ∀t ∈ [0, T], x, ω ∈ R , y ∈ R, z ∈ R × , to be consistent with [4], suppos 5, 10, 20, 40, N = 20 , T = 1, K = 1, r = 0.05, is the risk-free interest rate, σ dX = μ(t, x)X dt + σ(t, x)X dW , where W is a d-dimension Brownian motion mensions, g(x) = K − ∏ X is the payoff of option. At time t, the value of can option u(t, x) satisfies:

American Options
In this subsection, we consider the solution of reflected PDEs with continuous boundary. We select the case of pricing American options in order to make a comparison with the test results in [4]. The significant results are presented as well. The parameters and format of the American options are as follows: ∀t ∈ [0, T], x, ω ∈ R d , y ∈ R, z ∈ R 1×d , to be consistent with [4], suppose that d = 5, 10,20,40 , N = 20, T = 1, K = 1, r = 0.05, is the risk-free interest rate, is the payoff of option. At time t, the value of an American option u(t, x) satisfies: where T t,T is defined as the set of stopping times; it is also a solution to the PDEs following with a boundary where x) Now we display the numerical results for pricing American options by three models: The results in Table 4 show that both the Deep C-N method and the RDBDP method have great approximation accuracy up to 40 dimensions, while the Deep BSDE method has exploding loss values or eventually leads to the wrong results. Therefore, it is concluded that both the Deep C-N method and the RDBDP method have excellent performance when dealing with reflected PDE problems up to 40 dimensions, while the Deep BSDE method cannot solve the high-dimensional reflected PDE problems. It has been proved in [4] that the RDBDP algorithm has an obvious limitation when solving PDEs with reflection over 40 dimensions. Next, we demonstrate through Figure 4 that the Deep C-N method still has great stability when dealing with reflected PDE problems with more than 40 dimensions.
have great approximation accuracy up to 40 dimensions, while the Deep BSDE method has exploding loss values or eventually leads to the wrong results. Therefore, it is concluded that both the Deep C-N method and the RDBDP method have excellent performance when dealing with reflected PDE problems up to 40 dimensions, while the Deep BSDE method cannot solve the high-dimensional reflected PDE problems. It has been proved in [4] that the RDBDP algorithm has an obvious limitation when solving PDEs with reflection over 40 dimensions. Next, we demonstrate through Figure 4 that the Deep C-N method still has great stability when dealing with reflected PDE problems with more than 40 dimensions.

Discussion and Conclusions
In this paper, we proposed a deep learning-based numerical algorithm (named the Deep C-N method) for solving high-dimensional nonlinear reflected PDEs. Through numerical experiments, we found that the Deep C-N method has great approximating accuracy and computational stability in solving reflected PDEs of up to 100 dimensions. To obtain this result, we created a hybrid discrete format to solve some specific reflection problems. Compared with a single format, we found that this hybrid discretization strategy effectively overcomes the problem of loss function explosion or computational paralysis during the neural network computation. In addition, we used the penalization method to approximate the reflection terms instead of simulating them directly with numerical algorithms. The results show that the relative errors of the Deep C-N method are all below 4% when solving up to 40-dimensional reflected PDEs. As the dimension of the equation increases, the value of the neural network loss function shows an increasing trend, and the value of the neural network loss function stabilizes between 0.04 and 0.06 when the dimension rises to 100 dimensions. The Deep C-N method encountered challenges when trying to solve the reflection PDE above 100 dimensions. The first challenge is the computation time rising significantly with the increase in dimensions. The second is the computational accuracy decreasing as the dimensions increase. We support that these two problems are the most common difficulties encountered when dealing with highdimensional issues, and we can say that the key to developing numerical algorithms for high-dimensional problems is to overcome these two problems.
The method proposed in this paper contributes in two ways. First, in solving highdimensional reflection PDEs, one of the successful deep learning algorithms is the RDBDP algorithm proposed by Huré et al. [4], which can solve up to 40-dimensional reflected PDEs. The method proposed in this paper successfully extends the dimensions of the equation up to 100 dimensions. Second, in solving high-dimensional PDEs without reflection terms, the method proposed in this paper possesses higher approximating accuracy compared with the Deep BSDE methods proposed by Han et al. [22], Beck et al. [23], and Han et al. [24]. Therefore, the Deep C-N method proposed in this paper can be used to solve both reflected PDEs of up to 100 dimensions and high-dimensional nonlinear PDEs (tested up to 400 dimensions).
Since the discrete strategy we use has more steps, this inevitably leads to an increase in training difficulty and computational cost. One of the most direct evidence is that when approximating the 100-dimensional Allen-Cahn equation, the results can be converged in 426 s with 4000 training sessions using the Deep BSDE method, while with the Deep C-N method proposed in this paper, 10,000 training sessions are required to make the results converge in 3508 s (the exact running time varies with the machine).
Funding: This research received no external funding.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Conflicts of Interest:
The authors declare no conflict of interest.