[Paper Review] Multilevel Picard approximation algorithm for semilinear partial integro-differential equations and its complexity analysis

Ariel Neufeld, Sizhou Wu|arXiv (Cornell University)|May 19, 2022

Fractional Differential Equations Solutions4 citations

TL;DR

This paper introduces a multilevel Picard (MLP) approximation algorithm for solving high-dimensional semilinear parabolic partial integro-differential equations (PIDEs). By deriving a Feynman-Kac representation for the viscosity solution and leveraging Monte Carlo sampling with time discretization, the method achieves convergence without the curse of dimensionality, with computational complexity bounded polynomially in dimension $d$ and inverse accuracy $
$. The approach is validated numerically in up to 10,000 dimensions.

ABSTRACT

In this paper we introduce a multilevel Picard approximation algorithm for semilinear parabolic partial integro-differential equations (PIDEs). We prove that the numerical approximation scheme converges to the unique viscosity solution of the PIDE under consideration. To that end, we derive a Feynman-Kac representation for the unique viscosity solution of the semilinear PIDE, extending the classical Feynman-Kac representation for linear PIDEs. Furthermore, we show that the algorithm does not suffer from the curse of dimensionality, i.e. the computational complexity of the algorithm is bounded polynomially in the dimension $d$ and the reciprocal of the prescribed accuracy $\varepsilon$. We also provide a numerical example in up to 10'000 dimensions to demonstrate its applicability.

Motivation & Objective

• To develop a numerical method for solving high-dimensional semilinear parabolic PIDEs that overcomes the curse of dimensionality.
• To establish a Feynman-Kac representation for the unique viscosity solution of the PIDE under study.
• To prove convergence of the multilevel Picard approximation to the viscosity solution.
• To analyze the computational complexity of the algorithm and show it scales polynomially in dimension $d$ and accuracy $
  ^{-1}$.
• To demonstrate the method’s practical applicability through a numerical example in up to 10,000 dimensions.

Proposed method

• The method is based on a multilevel Picard iteration scheme derived from a Feynman-Kac representation of the PIDE’s viscosity solution.
• The solution is represented as an expectation involving a jump-diffusion SDE with drift, diffusion, and Lévy jump components.
• Time discretization is applied to the SDE, and expectations are approximated using Monte Carlo sampling.
• The algorithm uses a nested structure where each level refines the approximation using a hierarchy of discretization and sampling parameters.
• The compensator integral for jumps is approximated via Monte Carlo sampling of the Lévy measure, following techniques from prior work on linear PIDEs.
• The convergence and complexity analysis rely on stochastic fixed-point arguments, Gronwall-type inequalities, and moment estimates for the discretized processes.

Experimental results

Research questions

• RQ1 Can a multilevel Picard scheme be constructed for semilinear PIDEs that avoids the curse of dimensionality?
• RQ2 Does the proposed algorithm converge to the unique viscosity solution of the PIDE under suitable regularity and growth conditions on the coefficients?
• RQ3 Can a Feynman-Kac representation be established for semilinear PIDEs with nonlocal (jump) terms, extending the classical linear case?
• RQ4 What is the computational complexity of the algorithm in terms of dimension $d$ and accuracy $
  $?
• RQ5 How well does the method perform in high-dimensional settings, such as in 10,000 dimensions?

Key findings

• The multilevel Picard algorithm converges to the unique viscosity solution of the semilinear PIDE under mild regularity and growth conditions on the coefficients.
• A novel Feynman-Kac representation is derived for the viscosity solution, extending classical results from linear to semilinear PIDEs with jumps.
• The computational complexity of the algorithm is bounded polynomially in both the dimension $d$ and the inverse of the prescribed accuracy $
  ^{-1}$, thus overcoming the curse of dimensionality.
• The method is numerically validated in a test case with up to 10,000 spatial dimensions, confirming its practical feasibility in high dimensions.
• The convergence analysis relies on moment estimates and stochastic fixed-point arguments, with rigorous bounds on the error of the discretized and sampled approximations.
• The error bounds depend on the dimension $d$, the time horizon $T$, and the jump intensity, with explicit control via parameters like $M$ (number of jump samples) and $N$ (time steps).