Source-linked AI summary
On multilevel Picard numerical approximations for high-dimensional nonlinear parabolic partial differential equations and high-dimensional nonlinear backward stochastic differential equations
Weinan E, Martin Hutzenthaler, Arnulf Jentzen, Thomas Kruse
TL;DR
High-dimensional nonlinear PDEs and BSDEs often lack explicit solutions, while available approximation approaches face limitations in their applicability or computational scaling. This paper studies multilevel Picard approximations through simulations on 100-dimensional nonlinear PDEs and a review of other methods. The simulations indicate favorable accuracy and speed, with complexity growing at most linearly in dimension and quartically in reciprocal accuracy for the tested examples.
Problem
High-dimensional nonlinear PDEs and BSDEs often lack explicit solutions, while methods for nonlinear problems face difficulty achieving polynomial computational growth in both dimension and reciprocal accuracy.
Source: §1 · pp. 2–3Method
The paper tests multilevel Picard approximations based on fixed-point iterations and multilevel Monte Carlo, and reviews other approximation methods for nonlinear PDEs and BSDEs.
Source: §1 · p. 3, §2.1 · p. 4Results
Simulations on several 100-dimensional nonlinear PDEs indicate at most linear growth in dimension and quartic growth in reciprocal accuracy, with satisfactory accuracy and speed.
Source: §1 · p. 3Takeaways & Limitations
The simulations support the applicability of multilevel Picard approximations to the tested high-dimensional PDE examples, including cases with non-smooth or non-Lipschitz nonlinearities or terminal conditions.
Source: §1 · p. 3The review notes that quantization-tree effort grows exponentially with dimension, while cubature error bounds may have constants depending on dimension.
Source: §4.2 · p. 20, §4.5 · p. 21
Abstract
from arXiv · showhide
Parabolic partial differential equations (PDEs) and backward stochastic differential equations (BSDEs) are key ingredients in a number of models in physics and financial engineering. In particular, parabolic PDEs and BSDEs are fundamental tools in the state-of-the-art pricing and hedging of financial derivatives. The PDEs and BSDEs appearing in such applications are often high-dimensional and nonlinear. Since explicit solutions of such PDEs and BSDEs are typically not available, it is a very active topic of research to solve such PDEs and BSDEs approximately. In the recent article [E, W., Hutzenthaler, M., Jentzen, A., and Kruse, T. Linear scaling algorithms for solving high-dimensional nonlinear parabolic differential equations. arXiv:1607.03295 (2017)] we proposed a family of approximation methods based on Picard approximations and multilevel Monte Carlo methods and showed under suitable regularity assumptions on the exact solution for semilinear heat equations that the computational complexity is bounded by $O( d \, \epsilon^{-(4+\delta)})$ for any $\delta\in(0,\infty)$, where $d$ is the dimensionality of the problem and $\epsilon\in(0,\infty)$ is the prescribed accuracy. In this paper, we test the applicability of this algorithm on a variety of $100$-dimensional nonlinear PDEs that arise in physics and finance by means of numerical simulations presenting approximation accuracy against runtime. The simulation results for these 100-dimensional example PDEs are very satisfactory in terms of accuracy and speed. In addition, we also provide a review of other approximation methods for nonlinear PDEs and BSDEs from the literature.
1 Introduction and main results
The introduction motivates approximating high-dimensional nonlinear PDEs and BSDEs, where explicit solutions are typically unavailable and existing numerical approaches have important limits. It presents multilevel Picard approximations and reports simulations on 100-dimensional examples alongside a review of other methods.
1 Introduction and main results
- High-dimensional nonlinear PDEs arising in derivative pricing often lack explicit solutions, creating demand for approximate computation.Models may incorporate features such as borrowing and lending rates, default risk, transaction costs, or uncertain parameters.Source: §1 · p. 2
- Deterministic approximation methods can be efficient in low or moderate dimensions, but their computational effort grows exponentially with dimension and can prevent termination even at low accuracy.Examples include finite element, finite difference, spectral Galerkin, and sparse grid methods.Source: §1 · p. 2
- Multilevel Monte Carlo methods yield polynomial dimension scaling for linear PDEs, but the classical Feynman-Kac representation does not directly extend this approach to nonlinear PDEs.For linear PDEs, the effort grows up to an arbitrarily small order quadratically in reciprocal precision.Source: §1 · p. 2
- BSDE discretization methods require approximating nested conditional expectations, and the reviewed approaches do not achieve polynomial growth in both dimension and reciprocal accuracy.The text treats solving semilinear PDEs at individual space-time points and nonlinear BSDEs at individual times as essentially equivalent via the generalized nonlinear Feynman-Kac formula.Source: §1 · pp. 2–3
- Branching diffusion methods offer a Monte Carlo approach with complexity O(d ε^-2), but require the PDE's terminal or initial condition to be quite small.The cited representations apply to some nonlinear PDEs with polynomial or more general analytic nonlinearities.Source: §1 · p. 3
- The simulations indicate at most linear growth in dimension and quartic growth in reciprocal accuracy for several 100-dimensional nonlinear PDEs, including cases with non-smooth or non-Lipschitz terms.The earlier complexity bound O(d ε^-(4+δ)) applies under suitable regularity assumptions on the exact solution; the simulations complement that analysis.Source: §1 · p. 3
1.1 Notation
- The notation section defines conventions for inner products, measurable-function spaces, L2 norms, Lipschitz functions, matrix operations, and positive parts.These include B(E) for the Borel sigma-algebra and diag(x) for the diagonal matrix with entries of x.Source: §1.1 · p. 3
2 Multilevel Picard approximations for high-dimensional semilinear PDEs
The section develops multilevel Picard approximations by recasting a PDE solution as a fixed point and approximating its iterations with multilevel numerical estimates. It then specializes the scheme and describes its stochastic, sampling, and quadrature components.
2.1 Abstract picture for multilevel Picard approximations
- The abstract construction approximates a fixed point through recursively computed Picard iterates and computable approximations to the fixed-point map.The method uses multilevel stages corresponding to fixed-point iterations, requiring the full history of earlier approximations; it aims to retain the iterates' convergence speed at moderate computational cost.Source: §2.1 · pp. 3–4
2.2 A fixed-point equation for semilinear PDEs
- The PDE solution and its spatial derivative are formulated as a fixed point using Feynman-Kac and Bismut-Elworthy-Li formulas.A zero-expectation term is included to slightly reduce variance in Monte Carlo approximation; expectations and time integrals are then discretized at different levels.Source: §2.2 · p. 5
2.3 The approximation scheme
- The scheme defines multilevel Picard approximations for semilinear PDEs in a stochastic setting with indexed Brownian motions and approximations of the forward processes.The formulation specifies the model and process inputs under its stated assumptions.Source: §2.3 · p. 5
2.4 Remarks on the approximation scheme
- The implementation combines Picard iteration, multilevel techniques, SDE discretization, Monte Carlo sampling, and quadrature for time integrals.The recursion requires all earlier approximation levels; sample counts and quadrature rules are indexed by iteration, level, and parameter choices.Source: §2.4 · pp. 5–6
2.5 Special case: semilinear heat equations
- The section states a specialization of the approximation scheme to semilinear heat equations under the listed assumptions.A proposition gives the resulting specialized form, with its proof omitted.Source: §2.5 · p. 6
2.6 Special case: geometric Brownian motion
- The section also specializes the scheme to a geometric Brownian motion, a case noted as common in financial engineering.A proposition gives the specialized form under its stated assumptions; the proof is omitted.Source: §2.6 · p. 6
3 Numerical simulations of high-dimensional nonlinear PDEs
This section uses numerical simulations to evaluate algorithm (8) on nonlinear PDE examples across dimensions up to 100, examining approximation errors or increments, empirical convergence, and runtime.
3 Numerical simulations of high-dimensional nonlinear PDEs
- The simulations test five PDEs: solutions are unknown in Subsections 3.1–3.4 and explicitly known in 3.5; dimensions 1 and 100 are used in the first four, and runtimes are recorded for dimensions 5–100.For the one-dimensional cases, errors are measured against finite-difference approximations; in 100 dimensions, increments are examined for Subsections 3.1–3.4 and error against the exact solution for 3.5.Source: §3 · p. 7
- One-dimensional errors are compared with high-precision finite-difference approximations, while 100-dimensional increments are plotted against runtime for Subsections 3.1–3.4.For Subsection 3.5, the 100-dimensional error is instead measured against the explicitly known exact solution; the figures plot these quantities against average runtime.Source: §3 · pp. 7, 10
- Results in Subsections 3.1–3.3 average 10 independent simulation runs, whereas Subsections 3.4–3.5 use one run with an averaging parameter set to 10.The first two runs in 3.1–3.3 use different pseudorandom seeds.Source: §3 · p. 7
- Subsections 3.1–3.3 use fewer Monte Carlo samples for the nonlinearity because its variance is lower than the terminal-condition variance, allowing 7 Picard iterations rather than 5 in Subsections 3.4–3.5.The paper states that reducing averaging for the nonlinearity saves computational effort while maintaining Monte Carlo error of the same magnitude as for the terminal condition.Source: §3 · p. 9
- The simulations report runtimes and approximation statistics across parameter settings, while figures examine runtime against dimension from d = 5 to 100.Runtime plots cover the examples in Subsections 3.1–3.5; some average 20 realizations for smoother results.Source: §3 · p. 11
3.4 Allen-Cahn equation
- For the Allen-Cahn equation, empirical convergence is close to 1/4 at d = 1 and 1/3 at d = 100.These rates are suggested by the left and right sides of Figure 5, respectively.Source: §3.4 · p. 16
3.5 An example with an explicit solution
- A three-dimensional version of this explicit-solution example was considered by Chassagneux.Source: §3.5 · p. 17
- The explicit-solution example shows an empirical convergence rate close to 1/4 in Figure 7; the simulations are also reported in Table 9.The example uses T = 0.5 and dimension d = 100.Source: §3.5 · p. 17
4 Discussion of approximation methods from the literature
The discussion excludes deterministic methods because their effort grows exponentially in high dimensions, then reviews stochastic approaches and compares their computational effort under selected parameter choices.
4 Discussion of approximation methods from the literature
- Deterministic methods are excluded because their effort grows exponentially, making them unsuitable for high-dimensional nonlinear second-order parabolic PDEs.The discussion therefore focuses on stochastic approximation methods.Source: §4 · p. 17
- The reviewed stochastic methods often discretize BSDEs into nested conditional expectations, then approximate those expectations with Monte Carlo, quantization, regression, projection, or cubature methods.The review also covers Wiener chaos, branching diffusion, and density-estimation methods.Source: §4 · p. 18
- For the discussed algorithms other than branching diffusion, computational effort fails to grow polynomially in dimension and reciprocal best-known error bound under the chosen parameter settings.The settings optimize the best-known upper error bound; the truly optimal parameters are unknown because lower error estimates are unavailable.Source: §4 · p. 18
4.1 The ’straight-forward’ Monte Carlo method
The straightforward Monte Carlo method approximates nested conditional expectations in backward Euler-type BSDE schemes, but its computational effort grows at least exponentially in the reciprocal of the error bound.
4.1 The ’straight-forward’ Monte Carlo method
- Straightforward Monte Carlo approximates each conditional expectation in a backward Euler-type scheme with Monte Carlo simulations, producing nested averages.With time grid π and M samples per conditional expectation, computing Y^π,M requires (Md)^|π| scalar standard-normal realizations.Source: §4.1 · p. 19
- The required computational effort grows at least exponentially in the reciprocal of the right-hand side of (27), suggesting at most logarithmic convergence.The cited discussion reports no known literature statement claiming polynomial convergence for this method.Source: §4.1 · p. 19
4.2 The quantization tree method
The quantization tree method discretizes time and approximates conditional expectations using Monte Carlo-derived transition probabilities on a space-time grid. The grid must grow with the number of time steps, making the computational effort exponential in dimension.
4.2 The quantization tree method
- The method uses an explicit backward Euler-type time discretization and Monte Carlo-computed transition probabilities on a space-time grid to approximate conditional expectations.The underlying forward diffusion is projected onto the chosen grid.Source: §4.2 · p. 20
- At least Nd|π| scalar standard normal realizations are required for a grid with N space nodes and |π| time points.Source: §4.2 · p. 20
- Keeping the stated upper bound from exploding as |π| grows requires N ≥ |π|^(d+1), which makes computational effort exponential in dimension.The authors report finding no literature statement that some parameter choice yields polynomial effort in both dimension and reciprocal accuracy.Source: §4.2 · p. 20
4.3 The Malliavin calculus based regression method
The Malliavin calculus based regression method uses an implicit backward Euler-type scheme, but its iterated Skorohod integrals lead to exponential computational effort in dimension. A kernel-estimation alternative lacks a reported upper error estimate for choosing its bandwidth.
4.3 The Malliavin calculus based regression method
- The method is based on the implicit backward Euler-type method and involves iterated Skorohod integrals.Source: §4.3 · p. 20
- Numerically computing the iterated Skorohod integrals requires 2d independent standard normal variables, and the resulting effort grows exponentially with dimension.The authors report no known approximation method for these integrals whose effort does not grow exponentially in dimension.Source: §4.3 · p. 20
- A kernel-estimation alternative for the conditional expectations has no reported upper error estimate, leaving bandwidth selection from the number of time-grid points undetermined.Source: §4.3 · p. 20
4.4 The projection on function spaces method
The projection-on-function-spaces method estimates conditional expectations by projecting random variables onto a finite-dimensional basis and then using Monte Carlo. The cited basis choices leave high-dimensional computational guarantees unresolved.
4.4 The projection on function spaces method
- The method estimates conditional expectations by projecting random variables onto a finite-dimensional function space and estimating those projections by Monte Carlo.Projection error and computational effort generally depend on the basis choice.Source: §4.4 · p. 20
- With hypercube-indicator basis functions, achieving projection error ε requires simulating ⌊cε^-(3+2d)|log(ε)|⌋ forward-diffusion paths, giving exponential effort in dimension at fixed accuracy.Source: §4.4 · p. 20
- The cosine-function basis was specified only in one dimension, leaving its computational effort in high dimensions unclear.Source: §4.4 · p. 20
- The authors report no literature claim that a choice of function spaces and other parameters yields polynomial effort in both dimension and reciprocal accuracy.Source: §4.4 · p. 20
4.5 The cubature on Wiener space method
Cubature on Wiener space constructs finitely supported approximations to the forward SDE distribution and combines them with an implicit backward Euler-type scheme. Complexity-reduction methods yield polynomial convergence rates, but joint polynomial complexity in dimension and reciprocal accuracy is not established in the cited literature.
4.5 The cubature on Wiener space method
- The method combines an implicit backward Euler-type scheme with cubature paths that approximate the distribution of the forward SDE solution.The path construction is independent of f and g and can be computed once and tabularized.Source: §4.5 · p. 20
- Without complexity reduction, calculating Y^πn requires at least (N_m,d)^n paths, with effort growing exponentially in the reciprocal of the stated error bound.Source: §4.5 · p. 21
- Tree-based branching and recombination reduce the number of nodes and yield approximations with polynomial convergence rates, but the error-bound constant may depend on dimension.Source: §4.5 · p. 21
- Reported cubature simulations were in dimensions 1 and 5, and the authors found no claim that complexity-reduced cubature is polynomial in both dimension and reciprocal accuracy.Source: §4.5 · p. 21
4.6 The Wiener chaos decomposition method
The Wiener chaos decomposition method combines Picard iterations with chaos-decomposition formulas for nested conditional expectations. Its complexity does not remain polynomial in both dimension and reciprocal accuracy as the chaos order increases.
4.6 The Wiener chaos decomposition method
- Increasing the chaos-decomposition order is required for the approximation error to converge to zero, so the stated complexity does not grow at most polynomially in both dimension and reciprocal accuracy.The complexity is O(K × M × p × (N × d)^(p+1)) for p chaos order, K Picard iterations, M Monte Carlo samples per conditional expectation, and N time steps.Source: §4.6 · p. 21
- The method uses Picard iterations of the associated BSDE and Wiener chaos decomposition formulas to evaluate its nested conditional expectations.It does not discretize time because Wiener integrals with explicitly known antiderivatives can be simulated exactly.Source: §4.6 · p. 21
- The paper reports no known literature result establishing a polynomial convergence rate for this method.The passage cites Remark 4.8 in Briand & Labart.Source: §4.6 · p. 21
4.7 The branching diffusion method
The section reviews branching diffusion approximations for nonlinear PDEs and their dimensional scope, convergence, and applicability conditions. In an Allen–Cahn example, approximation quality worsens as the terminal value increases and condition (30) ceases to hold.
4.7 The branching diffusion method
- The method represents polynomial-nonlinearity PDE solutions using branching diffusion expectations, and works in all dimensions without suffering from the curse of dimensionality by construction.Its convergence rate is 0.5 when the forward diffusion can be simulated exactly.Source: §4.7 · p. 21
- Polynomial approximations can be unavailable or require a priori bounds on the exact solution, and certain polynomial nonlinearities can cause finite-time blow-up.The method is known to work well when the approximating polynomial, time horizon, and terminal condition satisfy a stated condition preventing blow-up until time T.Source: §4.7 · pp. 21–22
- Convergence at rate 0.5 requires finite second moments for the stochastic representation, imposing a more restrictive condition than (30).Source: §4.7 · p. 21
- The section illustrates possible problems beyond finite-time blow-up with simulations for an Allen–Cahn equation, noting that theoretical results in this direction are unavailable.The example assumes T = 1, constant coefficients and terminal value, and nonlinearity f(t, x, y, z) = y − y^3.Source: §4.7 · p. 21
- Branching diffusion approximations become poor as g(0) increases from 0.1 to 0.7 and fail to give good approximations once condition (30) is not satisfied.Source: §4.7 · p. 22
4.8 Approximations based on density representations
The section discusses methods representing the PDE solution as the density of a measure obtained from a McKean–Vlasov-type SDE. It identifies high-dimensional density estimation as a difficulty for these approaches.
4.8 Approximations based on density representations
- The methods proposed in Chang, Liu, & Xiong and Le Cavil, Oudjane, & Russo [60] use a McKean–Vlasov-type SDE representation in which the solution is a measure density, creating a density-estimation difficulty in high dimensions.Source: §4.8 · p. 22
- For the approximation in [60], the cited upper bound becomes small only with a small bandwidth and with both n and N growing exponentially in dimension.Source: §4.8 · p. 22
- The estimates in the proofs of imply that the initial particle count in its branching-particle approximations must grow exponentially in dimension.Source: §4.8 · p. 22