The MLP algorithm

MLP(; M=4, L=4, K=10, mc_sample = NoSampling())

Multi level Picard algorithm.


  • L: number of Picard iterations (Level),
  • M: number of Monte Carlo integrations (at each level l, M^(L-l)integrations),
  • K: number of Monte Carlo integrations for the non local term
  • mc_sample::MCSampling : sampling method for Monte Carlo integrations of the non local term.

Can be UniformSampling(a,b), NormalSampling(σ_sampling), or NoSampling (by default).


Returns a PIDESolution object.


  • multithreading : if true, distributes the job over all the threads available.
  • verbose: print information over the iterations.

The MLP, for Multi-Level Picard iterations, reformulates the PDE problem as a fixed point equation through the Feynman Kac formula.

The MLP algorithm overcomes the curse of dimensionality, with a computational complexity that grows polynomially in the number of dimension (see M. Hutzenthaler et al. 2020).

`MLP` can only approximate the solution on a single point

MLP only works for PIDEProblem with x0_sample = NoSampling(). If you want to solve over an entire domain, you definitely want to check the DeepSplitting algorithm.

The general idea 💡

Consider the PDE

\[\partial_t u(t,x) = \mu(t, x) \nabla_x u(t,x) + \frac{1}{2} \sigma^2(t, x) \Delta_x u(t,x) + f(x, u(t,x)) \tag{1}\]

with initial conditions $u(0, x) = g(x)$, where $u \colon \R^d \to \R$.

Picard Iterations

The MLP algorithm observes that the Feynman Kac formula can be viewed as a fixed point equation, i.e. $u = \phi(u)$. Introducing a sequence $(u_k)$ defined as $u_0 = g$ and

\[u_{l+1} = \phi(u_l),\]

the Banach fixed-point theorem ensures that the sequence converges to the true solution $u$. Such a technique is known as Picard iterations.

The time integral term is evaluated by a Monte-Carlo integration

\[u_L = \frac{1}{M}\sum_i^M \left[ f(X^{x,(i)}_{t - s_{(l, i)}}, u_{L-1}(T-s_i, X^{x,( i)}_{t - s_{(l, i)}})) + u(0, X^{x,(i)}_{t - s_{(l, i)}}) \right].\]

But the MLP uses an extra trick to lower the computational cost of the iteration.

Telescope sum

The MLP algorithm uses a telescope sum

\[\begin{aligned} u_L = \phi(u_{L-1}) &= [\phi(u_{L-1}) - \phi(u_{L-2})] + [\phi(u_{L-2}) - \phi(u_{L-3})] + \dots \\ &= \sum_{l=1}^{L-1} [\phi(u_{l-1}) - \phi(u_{l-2})] \end{aligned}\]

As $l$ grows, the term $[\phi(u_{l-1}) - \phi(u_{l-2})]$ becomes smaller - and demands more calculations. The MLP algorithm uses this fact by evaluating the integral term at level $l$ with $M^{L-l}$ samples.

  • L corresponds to the level of the approximation, i.e. $u \approx u_L$
  • M characterises the number of samples for the monte carlo approximation of the time integral

Overall, MLP can be summarised by the following formula

\[\begin{aligned} u_L &= \sum_{l=1}^{L-1} \frac{1}{M^{L-l}}\sum_i^{M^{L-l}} \left[ f(X^{x,(l, i)}_{t - s_{(l, i)}}, u(T-s_{(l, i)}, X^{x,(l, i)}_{t - s_{(l, i)}})) + \mathbf{1}_\N(l) f(X^{x,(l, i)}_{t - s_{(l, i)}}, u(T-s_{(l, i)}, X^{x,(l, i)}_{t - s_{(l, i)}}))\right] \\ &\qquad + \frac{1}{M^{L}}\sum_i^{M^{L}} u(0, X^{x,(l, i)}_t)\\ \end{aligned}\]

Note that the superscripts $(l, i)$ indicate the independence of the random variables $X$ across levels.

Nonlocal PDEs

MLP can solve for non-local reaction diffusion equations of the type

\[\partial_t u = \mu(t, x) \nabla_x u(t, x) + \frac{1}{2} \sigma^2(t, x) \Delta u(t, x) + \int_{\Omega}f(x, y, u(t,x), u(t,y))dy\]

The non-localness is handled by a Monte Carlo integration.

\[\begin{aligned} u_L &= \sum_{l=1}^{L-1} \frac{1}{M^{L-l}}\sum_{i=1}^{M^{L-l}} \frac{1}{K}\sum_{j=1}^{K} \bigg[ f(X^{x,(l, i)}_{t - s_{(l, i)}}, Z^{(l,j)}, u(T-s_{(l, i)}, X^{x,(l, i)}_{t - s_{(l, i)}}), u(T-s_{l,i}, Z^{(l,j)})) + \\ &\qquad \mathbf{1}_\N(l) f(X^{x,(l, i)}_{t - s_{(l, i)}}, u(T-s_{(l, i)}, X^{x,(l, i)}_{t - s_{(l, i)}}))\bigg] + \frac{1}{M^{L}}\sum_i^{M^{L}} u(0, X^{x,(l, i)}_t)\\ \end{aligned}\]


In practice, if you have a non-local model you need to provide the sampling method and the number $K$ of MC integration through the keywords mc_sample and K.

  • K characterises the number of samples for the Monte Carlo approximation of the last term.
  • mc_sample characterises the distribution of the Z variables


  • Boussange, V., Becker, S., Jentzen, A., Kuckuck, B., Pellissier, L., Deep learning approximations for non-local nonlinear PDEs with Neumann boundary conditions. arXiv (2022)
  • Becker, S., Braunwarth, R., Hutzenthaler, M., Jentzen, A., von Wurstemberger, P., Numerical simulations for full history recursive multilevel Picard approximations for systems of high-dimensional partial differential equations. arXiv (2020)