Machine learning approximation algorithms for high-dimensional fully nonlinear partial differential equations and second-order backward stochastic differential equations

Motivation

• High-dimensional nonlinear PDEs arise in financially relevant models involving effects such as default risk, transaction costs, trading constraints, and Knightian uncertainty, but are exceedingly difficult to approximate.Their dimension often reflects portfolios containing many financial assets, creating strong practical demand for effective numerical algorithms.Source: §1 · p. 3

Problem

• Existing approximation methods generally require low dimension or restrictive settings such as small nonlinearities, small terminal conditions, or semilinear structure, leaving no practical algorithm known for high-dimensional fully nonlinear parabolic PDEs.Only a limited number of practical high-dimensional PDE algorithms had been developed in the cited literature.Source: §1 · p. 4

Contribution

• The proposed algorithm combines the PDE–2BSDE connection, a merged formulation, temporal discretization, neural-network approximation, and a stochastic optimization procedure to solve fully nonlinear PDEs.The merged formulation converts the PDE problem into a learning problem whose loss compares the prescribed 2BSDE terminal condition with a neural-network forward discretization.Source: §1 · p. 4

Scope

• The method can also sample Shige Peng’s nonlinear expectation in high space-dimensions.The paper reports applications to nonlinear expectations of G-Brownian motions, including the 100-dimensional case.Source: §1 · pp. 4–5

Novelty

• Unlike related BSDE-based deep-learning methods [33] [52], the method applies to fully nonlinear PDEs and nonlinear expectations rather than only semilinear PDEs.The distinction arises from using the fully nonlinear PDE–2BSDE connection of Cheridito et al. instead of the classical PDE–BSDE connection.Source: §1 · pp. 4–5

2 Main ideas of the deep 2BSDE method

• The method combines PDE–2BSDE connections, a merged PDE–2BSDE formulation, temporal forward discretization, deep neural networks, and stochastic-gradient optimization.Its derivation is mainly based on ideas in E, Han, & Jentzen [33] and Cheridito et al..Source: §2 · p. 5

2.1 Fully nonlinear second-order PDEs

• The method targets approximating u(0, ξ) for a scalar terminal-value fully nonlinear second-order PDE, with extensions possible to PDE systems.The terminal-value formulation is motivated by its connection to 2BSDEs and financial applications such as the Black-Scholes-Barenblatt equation.Source: §2.1 · p. 6

2.2 Connection between fully nonlinear second-order PDEs and 2BSDEs

• Under suitable smoothness and regularity hypotheses, the PDE solution satisfies a 2BSDE representation with Y_t = u(t, ξ + W_t) and Z_t = (∇_x u)(t, ξ + W_t).The connection holds almost surely for all t ∈ [0,T].Source: §2.2 · p. 6

2.3 Merged formulation of the PDE and the 2BSDE

• The merged formulation combines the PDE and 2BSDE equations into relations connecting the solution, stochastic processes, and Hessian terms over time intervals.It is derived by substituting the PDE-related identities into the 2BSDE system.Source: §2.3 · p. 7

2.4 Forward-discretization of the merged PDE-2BSDE system

• The merged system is approximated forward on a partition 0 = t_0 < t_1 < … < t_N = T with sufficiently small mesh size.For sufficiently large N, the discretized relations approximate the continuous system at each time step.Source: §2.4 · p. 7

2.5 Deep learning approximations

• Deep neural networks approximate the Hessian-related and derivative terms at each time step, while θ_1 and (θ_2, …, θ_{d+1}) represent u(0, ξ) and (∇_x u)(0, ξ).The described networks use four layers, two hidden layers, and rectifier activation functions.Source: §2.5 · pp. 8–9

2.6 Stochastic gradient descent-type optimization

• Stochastic gradient descent minimizes the expected loss, producing random parameter iterates whose first component approximates u(0, ξ) for sufficiently large discretization and network sizes.The approximation relies on suitable assumptions and sufficiently large iteration counts.Source: §2.6 · pp. 9–10

2.7 Framework for the algorithm in a specific case

• The specific algorithmic framework uses Euler–Maruyama discretization and plain stochastic gradient descent with constant learning rate γ, while more general variants allow batch normalization or Adam.Under suitable assumptions, sufficiently large N, ν, and m and sufficiently small γ yield parameters approximating a local loss minimum and u(0, ξ).Source: §2.7 · pp. 10–11

3.1 Fully nonlinear second-order PDEs

• The method targets approximation of the fully nonlinear PDE solution u(0, x) through its connection with a second-order BSDE.The stated objective is to compute the function x ↦ u(0, x) in high-dimensional space.Source: §3.1 · p. 12, §3.2 · p. 12

3.2 Connection between fully nonlinear second-order PDEs and 2BSDEs

• Lemma 3.1 identifies the 2BSDE variables with the PDE solution, gradient, Hessian, and associated drift along a general Itô process.Specifically, Y_t = u(t, X_t), Z_t = (∇_x u)(t, X_t), and Γ_t = (Hess_x u)(t, X_t), with the resulting processes satisfying the 2BSDE system.Source: §3.2 · pp. 12–14

3.3 Merged formulation of the PDE and the 2BSDE

• The PDE and 2BSDE are merged by representing Y, Z, Γ, and A through u, its gradient, Hessian, and generator applied along X.This produces the merged formulation in equations (55)–(56), extending the corresponding construction beyond the Brownian-motion specialization.Source: §3.3 · p. 14

3.4 Forward-discretization of the merged PDE-2BSDE system

• The merged PDE-2BSDE system is forward-discretized on a sufficiently fine time grid 0 = t0 < t1 < ... < tN = T.The discretization initializes X_t0, Y_t0, and Z_t0 using ξ, u(0, ξ), and (∇_x u)(0, ξ).Source: §3.4 · p. 15

3.5 Deep learning approximations

• Neural approximations represent the initial value and gradient, each time-step’s Y and Z variables, and the Hessian and drift-related terms.The construction does not approximate u(t_n, x) or (∇_x u)(t_n, x) separately at intermediate times, instead propagating corresponding stochastic processes.Source: §3.5 · pp. 16–17

3.6 Stochastic gradient descent-type optimization

• Stochastic gradient descent minimizes the constructed loss and yields random parameter iterates treated, under suitable assumptions, as approximations of a local minimum.For sufficiently large discretization, network, and iteration parameters, the resulting random function approximates the target PDE-related function.Source: §3.6 · p. 17

3.7 Framework for the algorithm in the general case

• Framework 3.2 generalizes the method to broader Itô dynamics, alternative optimizers such as Adam, and machine-learning techniques such as batch normalization.It provides a unified framework covering the preceding subsections and is illustrated by examples in Section 4.Source: §3.7 · pp. 17, 19

4 Examples

• The examples evaluate deep 2BSDE approximations for Allen-Cahn, Black-Scholes-Barenblatt, Hamilton-Jacobi-Bellman equations, and nonlinear G-Brownian expectations in dimensions up to 100.Implementations use Framework 2.1 or Framework 3.2 in Python with TensorFlow.Source: §4 · p. 20

4.1 Allen-Cahn equation with plain gradient descent and no batch normalization

• Table 1 reports deep2BSDE numerical simulations in Framework 2.1 for the 20-dimensional Allen-Cahn equation, based on 10 independent realizations and relative L1-approximation errors.Figure 1 depicts relative L1-approximation errors and approximations of the mean loss function associated with Θ(1).Source: §4.1 · pp. 21–22, Table 1 · p. 22

4.2 Setting for the deep 2BSDE method with batch normalization and the Adam optimizer

• The batch-normalized setting uses Framework 3.2 with the Adam optimizer and evaluates the 20-dimensional Allen-Cahn problem over 5,000 iterations.The Adam parameter ε is 10^-8 by default, while Subsection 4.5 uses ε = 1.Source: §4.2 · pp. 22–23

4.3 A 100-dimensional Black-Scholes-Barenblatt equation

• The Framework 3.2 method approximates a 100-dimensional Black-Scholes-Barenblatt equation, with errors, empirical loss, and runtime evaluated over 400 iterations and ten realizations.Relative-error calculations use the analytic reference value 77.1049 from Lemma 4.2.Source: §4.3 · pp. 23–25, Table 2 · p. 24

4.4 A 100-dimensional Hamilton-Jacobi-Bellman equation

• The method also approximates a 100-dimensional Hamilton-Jacobi-Bellman equation with a quadratic gradient nonlinearity using ten independent runs through 2,000 iterations.The reported relative errors use the reference value 4.5901 calculated with Lemma 4.2 in [33] and classical Monte Carlo.Source: §4.4 · pp. 26–27, Table 3 · p. 28, Figure 3 · p. 28

4.5 A 50-dimensional Allen-Cahn equation

• A Framework 3.2 implementation approximates a 50-dimensional Allen-Cahn equation with cubic nonlinearity, evaluating errors, losses, and runtime over 2,000 iterations.This experiment uses Adam with ε = 1 and the reference value 0.09909 from the branching diffusion method.Source: §4.5 · pp. 27, 29, Table 4 · p. 29, Figure 4 · p. 30

4.6 G-Brownian motions in 1 and 100 space-dimensions

• The method approximates nonlinear expectations for 1- and 100-dimensional G-Brownian motions, using an explicit solution in 100 dimensions and finite-difference comparisons in one dimension.The 100-dimensional reference value is 162.5, obtained from Lemma 4.3.Source: §4.6 · p. 30, Table 5 · p. 31
• The experiments measure means, standard deviations, relative L1-approximation errors, empirical losses, and runtimes across training iterations using 10 independent runs.The 100-dimensional study varies m from 0 to 2000, while the one-dimensional study varies m from 0 to 500.Source: §4.6 · pp. 31, 33
• For the 100-dimensional G-Brownian motion, Table 5 and Figure 5 track relative L1-approximation errors and mean empirical losses for m up to 2000.The error calculations use the reference solution value 162.5, obtained from Lemma 4.3 with c = 1, σmax = 1, T = 1, σmin = 1/2, and d = 100.Source: §4.6 · pp. 31–32, Figure 5 · p. 32
• Lemma 4.3 states that, for the functions defined on [0,T] × R^d, u is infinitely differentiable and satisfies the terminal condition u(T,x)=g(x).Source: §4.6 · p. 32
• For the 1-dimensional G-Brownian motion, Table 6 and Figure 6 report relative L1-approximation errors and mean empirical loss for m ∈ {0, 1, . . . , 500}, based on 10 independent runs.The setup uses d = 1, g(x) = 1/(1 + exp(−x^2)), γ_m = 1/100, and ξ(ω) = −2; u(0, −2) is approximated as 0.90471 by finite differences.Source: Table 6 · p. 33, §4.6 · pp. 33–34, Figure 6 · p. 34

A Source codes

• The cited passage presents Python code described as a simplified version of Python code 3 in Subsection A.3.Source: §A.1 · p. 34

A.1 A Python code for the deep 2BSDE method used in Subsection 4.1

• The code is a plain deep2BSDE solver with a hard-coded Allen–Cahn equation.Source: §A.1 · p. 34
• The implementation defines the PDE right-hand side through f0(t,X,Y,Z,Gamma), including the term -Y+tf.pow(Y, 3).Source: §A.1 · p. 35
• Helper functions construct time-specific neural networks for the stochastic variables, including Z and Gamma.Source: §A.1 · p. 35
• A forward discretization advances the stochastic system across N−1 time steps using simulated Brownian increments and updates Y with the PDE driver.Source: §A.1 · pp. 36–37
• Training computes gradients of the loss function, applies stochastic-gradient updates in a training loop, and records loss, Y0, learning rate, and running time.Source: §A.1 · pp. 37–38
• The code uses plain stochastic gradient descent without batch normalization and writes the recorded results to a CSV file.Source: §A.1 · p. 38

• The Matlab implementation is a slightly modified version of the code from E, Han, & Jentzen [33, Subsection 6.2].Source: §A.2 · p. 38
• The code evaluates the branching PDE solver through MCBranchingEvaluation and reports an estimate of u(0,x0) with its standard deviation.Source: §A.2 · p. 39
• Its branching-process routine advances states using drift, diffusion, and Gaussian increments, while tracking the selected nonlinearity and recursively generating descendant processes.Source: §A.2 · p. 40
• The appendix also labels a second Matlab listing as the branching diffusion method used in Subsection 4.5.Source: §A.2 · p. 40

• The code is based on E, Han, & Jentzen [33, Subsection 6.1] and implements a deep 2BSDE solver for a hard-coded Black-Scholes-Barenblatt equation.Source: §A.3 · pp. 40–41
• The implementation defines the Black-Scholes-Barenblatt dynamics and driver through volatility, diffusion, drift, and nonlinear f_tf functions.Source: §A.3 · p. 41
• It represents the solution components with time-indexed neural networks built from multilayer layers, ReLU activation, and batch normalization.Source: §A.3 · p. 42
• A forward loop updates X, Y, and Z across time steps using simulated Brownian increments, the driver f_tf, and Gamma-dependent terms.Source: §A.3 · pp. 43–44
• Training uses exponentially decaying learning rates and records losses, Y0 values, learning rates, running times, and CSV output.Source: §A.3 · pp. 44–45

A.4 A Matlab code for the classical Monte Carlo method used in Subsection 4.4

• The code slightly modifies the Matlab implementation from E, Han, & Jentzen [33, Subsection 6.3].Source: §A.4 · p. 46
• The supplied Matlab implementation is the classical Monte Carlo method used in Subsection 4.4.Source: §A.4 · p. 46

Implementation

• The Matlab listing is explicitly based on code by E et al. [35].The section states that it is inspired by MATLAB code 7 in Section 3 of that work.Source: §A.5 · p. 46
• The code defines the nonlinear function f using a gamma-dependent expression involving sigma_max.The displayed implementation begins the definition as f = @(t,x,y,z,gamma) 0.5 * ((gamma >0).* sigma_max ...Source: §A.5 · p. 46
• The implementation constructs sparse tridiagonal matrices with gallery('tridiag') and specified diagonal entries for the finite-difference operators.The displayed lines include diagonal patterns built from ones, -2*ones, and zeros, followed by boundary-row entries.Source: §A.5 · pp. 46–47
• The appendix supplies the Matlab implementation of the finite-differences method used in Subsection 4.6.The listing is identified as Matlab code 5 and includes the GBrownianMotion function.Source: §A.5 · pp. 46–47