Hybrid PDE solver for data-driven problems and modern branching

TL;DR

This paper introduces the Probabilistic Domain Decomposition (PDD) method, a scalable parallel approach for solving large-scale PDEs in data-driven applications, leveraging stochastic representations and Monte Carlo techniques.

Contribution

It presents the PDD method, combining probabilistic and deterministic PDE solvers, and discusses recent advances in stochastic representations for nonlinear PDEs using branching diffusions.

Findings

PDD exhibits excellent scalability for large-scale PDEs.

Recent stochastic representations enable solving nonlinear PDEs with branching diffusions.

The paper provides a comprehensive overview of advanced parallel algorithms for PDEs.

Abstract

The numerical solution of large-scale PDEs, such as those occurring in data-driven applications, unavoidably require powerful parallel computers and tailored parallel algorithms to make the best possible use of them. In fact, considerations about the parallelization and scalability of realistic problems are often critical enough to warrant acknowledgement in the modelling phase. The purpose of this paper is to spread awareness of the Probabilistic Domain Decomposition (PDD) method, a fresh approach to the parallelization of PDEs with excellent scalability properties. The idea exploits the stochastic representation of the PDE and its approximation via Monte Carlo in combination with deterministic high-performance PDE solvers. We describe the ingredients of PDD and its applicability in the scope of data science. In particular, we highlight recent advances in stochastic representations for nonlinear PDEs using branching diffusions, which have significantly broadened the scope of PDD. We envision this work as a dictionary giving large-scale PDE practitioners references on the very latest algorithms and techniques of a non-standard, yet highly parallelizable, methodology at the interface of deterministic and probabilistic numerical methods. We close this work with an invitation to the fully nonlinear case and open research questions.

Figures

Tables

Table 1. Table 2.1: Acceleration of the most current version of PDD (IterPDD), with improved stochastic numerics, over the previous PDD algorithm. a 0 subscript 𝑎 0 a_{0} is the error tolerance on interfacing nodes. The BVP being solved is ( 2.10 ).

$a_0$	theoretical speedup	observed speedup (approx.)
.04	13.93	15
.02	28.34	29
.01	57.42	60
.005	116.00	125
.0025	233.84	250

Equations

$$\Delta u(\mathbf{x})=0\text{ if x∈Ω⊂Rd},u(\mathbf{x})=g(\mathbf{x})\text{ if x∈∂Ω}.$$

$$\displaystyle u(\mathbf{x})=\mathbb{E}_{\mathbf{x},0}\big{[}\,g({\bf X}_{\tau})\,\big{]}$$

$$\text{for }t\ge 0\mathrm{d}{\mathbf{X}}_t$$

$$\left\{\begin{array}[]{ll}\frac{\partial u}{\partial t}=\mathcal{L}u+c({\bf x},t)u+f({\bf x},t)&\textrm{ if }t>0,{\bf x}\in\Omega,\\ u=p({\bf x})&\textrm{ if }t=0,{\bf x}\in\Omega,\\ u=g({\bf x},t)&\textrm{ if }t>0,{\bf x}\in\partial\Omega_{A},\\ \frac{\partial u}{\partial N}={\varphi}({\bf x},t)u+\psi({\bf x},t)&\textrm{ if }t>0,{\bf x}\in\partial\Omega_{R}.\end{array}\right.$$

$$\mathcal{L}:=\frac{1}{2}\sum _{i,j=1}^d{A}_{ij}(\mathbf{x},t)\frac{{\partial }^2}{\partial x_i\partial x_j}+\sum _{i=1}^d{b}_i(\mathbf{x},t)\frac{\partial }{\partial x_i},$$

$$\varphi$$

$$+\underset{0}{\overset{\min (T,\tau )}{\int }}f({\mathbf{X}}_t,T-t)\Phi (t)dt+\underset{0}{\overset{\min (T,\tau )}{\int }}\psi ({\mathbf{X}}_t,T-t)\Phi (t)d{\xi }_t,$$

$$\Phi(t)=\exp{\Big{(}\,\int_{0}^{t}c({\bf X}_{s},s)ds+\int_{0}^{t}{\varphi}({\bf X}_{s},s)d\xi_{s}\,\Big{)}}.$$

$${\mathbb{E}}\big{[}\phi\big{|}{\bf X}_{0}={\bf x}_{0}\big{]}={\mathbb{E}}\Big{[}\,q({\bf X}_{\tau})Y_{\tau}+Z_{\tau}\,\Big{]},\ \textrm{s.\leavevmode\nobreak\ th. }\ q({\bf X}_{\tau})=\left\{\begin{array}[]{ll}g({\bf X}_{\tau},T-\tau),&\textrm{ if }\tau<T,\\ p({\bf X}_{T}),&\textrm{ if }\tau\geq T,\end{array}\right.$$

$$\left\{\begin{array}[]{ll}d{\bf X}_{t}={\bf b}({\bf X}_{t},T-t)dt+\sigma({\bf X}_{t},T-t)d{\bf W}_{t}-{\bf N}({\bf X}_{t})d\xi_{t}&{\bf X}_{0}={\bf x}_{0},\\ dY_{t}=c({\bf X}_{t},T-t)Y_{t}dt+{\varphi}({\bf X}_{t},T-t)Y_{t}d\xi_{t}&Y_{0}=1,\\ dZ_{t}=f({\bf X}_{t},T-t)Y_{t}dt+\psi({\bf X}_{t},T-t)Y_{t}d\xi_{t}&Z_{0}=0,\\ d\xi_{t}={\mathbbm{1}}_{\{{\bf X}_{t}\in\partial\Omega_{R}\}}dt&\xi_{0}=0.\end{array}\right.$$

$$u(t,{\bf x}_{0})={\mathbb{E}}\big{[}\phi\big{|}{\bf X}_{0}={\bf x}_{0}\big{]}={\mathbb{E}}\Big{[}\,q({\bf X}_{\tau})Y_{\tau}+Z_{\tau}\,\Big{]}.$$

$$\displaystyle\nabla^{2}u+\frac{\cos{(x+y)}}{1.1+\sin{(x+y)}}\big{(}\frac{\partial u}{\partial x}+\frac{\partial u}{\partial y}\big{)}-\frac{x^{2}+y^{2}}{1.1+\sin{(x+y)}}u+f(x,y)=0,$$

$$S(A,B)=\frac{\text{Time taken by method B}}{\text{Time taken by method A}}.$$

$$\{\begin{array}{l}\frac{\partial u}{\partial t}-\frac{{\partial }^{2}u}{\partial x^2}-u(u-1)=0,x\in (x_1,x_2),\text{leavevmode}{x}_1,x_2\in \mathbb{R},\text{leavevmode}t>0,\\ u(x,0)=\psi (x),u(x_1,t)=g(x_1,t),u(x_2,t)=g(x_2,t).\end{array}$$

$$\frac{\partial u}{\partial t}-\mathcal{L}u-c\left(\sum _{i=0}^{\infty }{\alpha }_i{u}^i-u\right)=0,$$

$$\displaystyle u(\mathbf{x},t)=\mathbb{E}_{\mathbf{x},0}\left[\prod_{i=1}^{N_{t}}\left(\psi\big{(}X^{(i)}_{t}\big{)}\mathbbm{1}_{\{t<\tau_{\partial\Omega}\}}+g\big{(}X^{(i)}_{\tau_{\partial\Omega}},t-\tau_{\partial\Omega}\big{)}\mathbbm{1}_{\{t\geq\tau_{\partial\Omega}\}}\right)\right]\,,$$

$$\frac{\partial u}{\partial t}-\frac{{\partial }^{2}u}{\partial x^2}-\frac{1}{2}\left(1+u^2\right)+u=0,u(x,0)=\psi (x),x\in \mathbb{R},\text{leavevmode}t>0.$$

$$\frac{\partial u}{\partial t}-D\frac{{\partial }^{2}u}{\partial x^2}+ru(1-u)=0,(x,t)\in \mathbb{R}\times [0,\infty )$$

$$\displaystyle u(x,0)=\psi(x)=1-\Big{(}1+\exp\big{(}x/\sqrt{6}\big{)}\Big{)}^{-2}\,,$$

$$\displaystyle u(x,t)=1-\Big{(}1+\exp\big{(}\dfrac{x}{\sqrt{6}}-\frac{5t}{{6}}\big{)}\Big{)}^{-2}\,.$$

$$\displaystyle u(x,t)=\mathbb{E}_{\mathbf{x},0}\left[\prod_{i=1}^{N_{t}}u\big{(}W^{(i)}_{t}+x,0\big{)}\right]=\mathbb{E}_{\mathbf{x},0}\left[\prod_{i=1}^{N_{t}}\left(1-\left(1+\exp\left(\big{(}W^{(i)}_{t}+x\big{)}/\sqrt{6}\right)\right)^{-2}\right)\right]\,,$$

$$\frac{\partial u}{\partial t}+\mathcal{L}u+c\left(\sum _{i=0}^L{\alpha }_i(\mathbf{x},t)u^i-u\right)=0,u(\mathbf{x},T)=\psi (\mathbf{x}).$$

$$\displaystyle l_{0}(s):=\sum_{k\geq 0}\|\alpha_{k}\|_{\infty}s^{k}\quad\text{and}\quad l(s):=c\Big{(}\,\frac{l_{0}\big{(}s\|\psi\|_{\infty}\big{)}}{\|\psi\|_{\infty}}-s\,\Big{)},$$

$$u(\mathbf{x},t)={\mathbb{E}}_{\mathbf{x},t}\left[\prod _{k\in {\mathcal{K}}_{T-t}}\psi ({\mathbf{X}}_T^{(i)})\prod _{n=1}^{M_{T-t}}\left(\frac{{\alpha }_{I_n}({\mathbf{X}}_{t+T_n}^{(K_n)},t+T_n)}{q_{I_n}}\right)\right].$$

$$\displaystyle\max\Big{(}\frac{\partial u}{\partial t}+\mathcal{L}u,\varphi(\mathbf{x})-u\Big{)}=0,\qquad u(\mathbf{x},T)=\varphi(\mathbf{x}),\qquad(\mathbf{x},t)\in\mathbb{R}^{d}\times[0,T].$$

$$\frac{\partial u}{\partial t}+\mathcal{L}u=\text{mathbbm}{1}_{\{\phi (\mathbf{x})\ge u\}}\mathcal{L}\phi ,u(\mathbf{x},T)=\phi (\mathbf{x}),(\mathbf{x},t)\in {\mathbb{R}}^d\times [0,T].$$

$$\frac{\partial u}{\partial t}+\mathcal{L}u+r_{0}u+r_1{u}^{+}=0,u(x,T)=\psi (x),x\in \mathbb{R}$$

$$\frac{\partial v}{\partial t}+\mathcal{L}v+r_{0}v+r_1{v}^{+}=0$$

$$\frac{\partial v}{\partial t}+\mathcal{L}v+c\left(F(v)-v\right)=0,v(x,T)=\psi /\Vert \psi {\Vert }_{\infty }$$

$$\frac{\partial v}{\partial t}+\mathcal{L}v+c\left(0.0586+0.5v+0.8199v^2-0.4095v^4-v\right)=0,\text{leavevmode}v(x,T)=\psi /\Vert \psi {\Vert }_{\infty }.$$

$$f(t,\mathbf{x},y,z)=\sum _{l=(l_0,\dots ,l_m)\in L}{c}_l(t,\mathbf{x})y^{l_0}\prod _{i=1}^m(h_i(t,\mathbf{x})\cdot z{)}^{l_i}.$$

Full text

We envision this work as a dictionary giving large-scale PDE practitioners references on the very latest algorithms and techniques of a non-standard, yet highly parallelizable, methodology at the interface of deterministic and probabilistic numerical methods. We close this work with an invitation to the fully nonlinear case and open research questions.

Keywords: Probabilistic Domain Decomposition, high-performance parallel computing, scalability, nonlinear PDEs, Marked branching diffusions, hybrid PDE solvers, Monte-Carlo methods.

**2010 AMS subject classifications:

**Primary: 65C05, 65C30,  Secondary: 65N55, 60H35, 91-XX, 35CXX

1 Introduction

Partial Differential Equations (PDEs) are ubiquitous in modelling, appearing in image analysis and processing, inverse problems, shape analysis and optimization, filtering, data assimilation and optimal control. They are used in Math-Biology to model population dynamics with competition or growth of tumours; or to model complex dynamics of movement of persons in crowds or to model (ir)rational decisions of players in games and they feature in many complex problems in Mathematical Finance. Underpinning all these applications is the necessity of solving numerically such equations either in bounded or unbounded domains.

Deterministic Domain Decomposition. The standard example of a Boundary Value Problem (BVP) is Laplace’s equation with Dirichlet Boundary Conditions (BCs):

[TABLE]

The large data sets involved in realistic applications nearly always imply that the discretization of a BVP such as (1.1) leads to algebraic systems of equations that can only be solved on a parallel computer with a large number (say $p>>1$) of processors. Not only does parallelization require multiple processors but also parallel algorithms. The classical Schwarz’s alternating method was the first and remains the paradigm of such algorithms which we refer to as “Deterministic Domain Decomposition” (DDD) [59]. While state-of-the-art DDD algorithms outperform Schwarz’s alternating method in every respect, the latter nonetheless serves to illustrate the crucial difficulty they all face. The idea of Schwarz’s algorithm is to divide $\Omega$ into a set of $p$ overlapping subdomains, and have processor $j=1,\cdots,p$ solve the restriction of the PDE to the subdomain, $\Omega_{j}$ see Fig. 1.1.

Since the solution is not known in the first place, the BCs on the fictitious interfaces of $\Omega_{j}$ are also unknown, therefore, an initial guess has to be made in order to give processor $j$ a well-posed (yet incorrect) problem. The BCs along the fictitious interfaces of $\Omega_{j}$ are then updated from the solution of the surrounding subdomains in an iterative way until convergence.

Inter-processor communication and the scalability limit of DDD. Since the inter-processor communication involved in DDD’s updating procedure is intrinsically sequential, it sets a limit to the scalability of the algorithm by virtue of the well known Amdahl’s law.

A simple illustrative example is as follows. A fully scalable algorithm would take half the time (say $T/2$) to run if the number of processors was doubled. If there is a fraction $\nu<1$ of the algorithm which is sequential, then the completion time will not drop below $T\nu$, regardless of how many processors are added. For instance, in Schwarz’s method, if $5\%$ (i.e. $\nu=.05$) of the execution time of one given processor is lost by waiting for the artificial BCs to be ready, then the execution time could be shortened by at most a factor of 20 (with infinitely many processors; a factor of about 19 would already take over 1000 processors). In other words, Schwarz’s alternating algorithm, or any DDD algorithm for that matter, cannot exploit the full capabilities of a parallel computer due to the idle time wasted in (essential) communication.

We emphasize this point further by borrowing an example from David Keyes111See http://www.mcs.anl.gov/research/projects/petsc-fun3d/Talks/bellTalk.ppt. The Gordon Bell prizes are annually awarded to numerical schemes which achieve a breakthrough in performance when solving a realistic problem. In 1999, one such problem was the simulation of the compressible Navier-Stokes equations around the wing of an airplane. With $128$ processors, the winning code took $43$ minutes to solve the task. On the other hand, with $3072$ processors it took it $2.5$ minutes instead of $1.79$ as would have been the case with a fully parallelizable algorithm. The remaining $28\%$ of computer time were lost to interprocessor communication. At this point, adding more processors would have led to a faster loss of scalability.

The Probabilistic Domain Decomposition (PDD) method. A conceptual breakthrough was achieved by Acebrón et al. with the PDD method (or rather, the PDD framework) [2], based on a previous, unpublished idea by Renato Spigler. PDD is the only domain decomposition method potentially free of communication, and thus potentially fully scalable. It does so by splitting the simulation in two separate stages, the first stage recasts the BVP into a stochastic formulation (via the so-called Feynman-Kac formula) which allows to compute the solution of the PDE at certain specific points in time/space. Thus, we can compute the “true” solution values of the DDD’s fictitious interfaces for the $\Omega_{j}$’s. Therefore, the fictitious boundaries that were previously unknown are now known!

Consequently, the subdomains are now completely independent of each other, the second stage then involves solving for the solution over the subdomains in a full parallel way. PDD calculations will be affected by two independent sources of numerical error: the subdomain solver and the statistical error of Monte Carlo simulations.

PDD is currently well understood for the linear case, although recent advances in stochastic representations have opened the door for using PDD with nonlinear PDEs. A class of nonlinear PDEs strongly amenable to PDD are those whose solution can be represented by the so-called branching processes, as introduced by [62], [58], [52] and recently extended by [44] and [43]. This methodology avoids backward regressions, the so-called “Monte Carlo of Monte Carlo simulation” problem, see Section 4 below or [42, Section 3.1]. To illustrate the potential of branching in PDD, numerical examples are worked out in Section 3. For general mildly nonlinear PDEs, the straightforward (and often efficient) approach of linearization and solution of each of the linear iterates, would be easily implementable with PDD, without resorting to nonlinear representations.

The general case of systems of nonlinear 2nd order parabolic/elliptic PDEs or fully nonlinear PDEs are, to the best of our knowledge, yet been addressed in the framework of PDD. These types of PDEs admit a probabilistic representation in terms of Forward Backward SDEs (FBSDEs [29]) and numerical methods for FBSDEs is currently the subject of extensive research.

The remainder of the paper is organized as follows.

• •

Section 2 is an overview of PDD. We illustrate the connection between stochastic and BVPs in the linear case, and provide pointers for the numerical methods required. We stress the notion of balancing the various aspects of PDD in order to speed up the overall algorithm. The section closes with a survey of reported results on PDD.

• •

In Section 3, nonlinear equations that can be represented probabilistically in terms of branching processes are discussed, and an account of recent developments is given. We give examples with convincing numerical tests using branching.

• •

Section 4 is an invitation to new, more challenging problems than those tackled so far. In particular, the stochastic representation of systems of nonlinear parabolic equations with Forward Backward SDEs (FBSDEs) is succinctly discussed. We highlight open research problems.

2 An Overview of Probabilistic Domain Decomposition

An alternative to deterministic methods which is specifically designed to circumvent the scalability issue of DDD is the PDD method [2], [3]. First, the domain $\Omega$ under consideration is divided into non-overlapping subdomains. (Rather than overlapping ones, such as with Schwarz’s alternating algorithm.) PDD then consists of two stages, firstly, the solution is calculated only on a set of interfacing nodes along the artificial interfaces, by solving the stochastic representation of the BVP with the Monte Carlo method, using the so-called Feynman-Kac representations. The stochastic representation is the crux of PDD and can be highly non-trivial. Nonetheless, it can also be extremely simple, for instance in (1.1), it is well known (see [46]) that the solution $u$ can be represented as

[TABLE]

where ${\bf X}$ is the solution to the simplest Stochastic Differential Equation (SDE)

[TABLE]

where $({\bf W}_{t})_{t\geq 0}$ is a $d$-dimensional Brownian motion; ${\bf X}_{\tau}$ is the point on $\partial\Omega$ where the trajectory of ${\bf X}_{t}$ first hits the boundary and ${\bf X}$ starts at $t=0$ from $\mathbf{x}$; and lastly, $\mathbb{E}[\cdot]$ is the usual expectation operator (in (2.1), $\mathbb{E}_{\mathbf{x},0}[\cdot]$ emphasizes that the diffusion $({\bf X}_{t})_{t\geq 0}$ starts at time $t=0$ from position $\mathbf{x}$.) See Fig.2.1 for an illustration.

The expected value above is approximated via a Monte Carlo method (the mean over many independent realizations of the SDE (2.2), which can be carried out in a fully parallelizable way), introducing some statistical error. A numerical scheme is required to solve the SDE, we refer to it as “stochastic numerics” [47].

By using the stochastic representation (2.1), one can compute the equation’s solution at each point at the artificial boundaries $\Omega_{i}\cap\Omega_{j}$. It is then possible to reconstruct (approximately) the solution on the interfaces (interpolating the values at interfacing nodes with Chebyshev polynomials, say), so that the PDE restricted to each of the subdomains is now well posed. Now, the $p$ Laplace’s equations on each of the $p$ subdomains are separate problems, and can be independently solved “deterministically”, the second stage of PDD.

Note that both stages in the PDD are embarrassingly parallel by construction. We stress that this has been made possible only by the unique capability of stochastic representations of PDEs for solving a BVP at an isolated point (time-space) of the domain. On the other hand, DDD methods are matrix-based and deliver the solution everywhere; but by doing this are bound to move around information.

Remark 2.1** (Further features: Fault Tolerance and boundary design)**

Another feature of PDD is that it is, by construction, fault-tolerant: the outage of an individual processor does not lead to substantial — let alone fatal — damage to the overall simulation. Given that modern parallel computers boast millions of processors, this is a definite plus (which stands in stark contrast to DDD algorithms).

Lastly, it is up to the user to choose the artificial boundaries. These can be chosen in such a way as to lessen effects of kinks, sharp corners in the numerics. For instance, by better interplaying with the solutions being provided by the stochastic representations and Monte Carlo simulation.

2.1 The ingredients of PDD

PDD requires four ingredients: an interpolator, a subdomain solver, a stochastic representation of the PDE at hand, and an efficient numerical method for exploiting the latter with the Monte Carlo method. Let us briefly comment on them.

$\triangleright$ The interpolator scheme joins the (approximate) solutions at the interface nodes (computed with Monte Carlo) in order to yield an (approximate) Dirichlet BC on the whole interface. (We remark that, in time-dependent problems, the interface runs along the time component, as well as the spatial ones.) A flexible and accurate interpolation scheme well suited to PDD is the RBF interpolation — see [13] and references therein. We note that less stringent approximation schemes than interpolation, such as “denoising” splines or least squares, are also possible.

$\triangleright$ The subdomain solver is any state-of-the-art numerical scheme to solve a PDE (or a BVP) on a finite domain, without parallelization. It takes the approximation to the solution on the boundary delivered by the interpolator as Dirichlet BCs and produces a solution everywhere within the subdomain. All the subdomain solutions are later “glued together” to yield a complete solution of the large-scale original problem. Common examples of solvers are finite differences, finite elements, and meshless methods.

$\triangleright$ The stochastic (or probabilistic) representation of the PDE, like that in (2.1), is the pointwise characterization of the solution of a PDE (stationary or time-dependent) as the expected value of the functional of a system of stochastic differential equations (SDEs). The paradigm of stochastic representations is the well-known Feynman-Kac formula. More details will follow in Section 2.2.

$\triangleright$ Efficient stochastic numerical methods to exploit the stochastic representation. With the Monte Carlo approach, the expected value which yields the solution at an interfacing node is approximated as the average over many independent samples and this has to be repeated for every interfacing node. Contrary to deterministic methods, in the stochastic setting the errors are twofold, brought about by both discretization of time and by the substitution of the expected value by a sample mean. (One can say that the estimator of the expected value is biased due to the discretization of the time variable.) In order to attain a target accuracy $a$, both sources of errors must be balanced and this leads to lengthy simulations: the computational effort is ${\cal O}(a^{-4})$ if naive numerics are used (such as the Euler-Maruyama scheme plus boundary test for linear BVPs, see [45]). This cost can be substantially reduced (to around ${\cal O}(a^{-2})$ in some cases) with sophisticated, novel schemes. The new integrators for bounded stopped diffusions put forward by Gobet [36] and Gobet and Menozzi [40] have doubled the convergence rate of the weak error associated to the interaction with the boundary, see also [12]. For reflected diffusions, Lèpingle’s method [36] or Milstein’s bounded methods [54] are recommended. Two further promising directions are the incorporation of Giles’ Multilevel method [45], [35] and the use of pathwise control variates for variance reduction [13]. In the later sections, pointers will be given to the numerics of nonlinear probabilistic representations. In any case, the numerics of SDEs (and generalizations thereof) are underdeveloped compared to those of deterministic equations. A positive aspect is that many new ideas being currently developed for SDEs in financial mathematics can be transplanted to PDD.

Ideally, a numerical method tailored to the probabilistic representation of the given PDE exists, resulting in huge gains in efficiency. An example is the Walk on Spheres algorithm for the stopped Brownian motion, which is the probabilistic representation of Laplace’s equation with Dirichlet BCs. This and similar, specific numerical schemes should be used whenever possible (see [14] for references).

2.2 Probabilistic representation of linear BVPs

We proceed now to give the probabilistic representation of linear BVPs (a generalization of the well known Feynman-Kac formula). Let us start by considering the general linear parabolic BVP of second order with mixed BCs:

[TABLE]

The notation $\partial\Omega_{A}$ stands for a portion of the boundary where Dirichlet BCs are imposed, while on $\partial\Omega_{R}=\partial\Omega\backslash\partial\Omega_{A}$, BCs involving the normal derivative (Neumann or Robin) hold. In SDE literature $\partial\Omega_{A}$ and $\partial\Omega_{R}$ are known as absorbing and reflecting boundaries respectively. In (2.3), $\mathcal{L}$ is the 2nd order differential operator defined by

[TABLE]

where the matrix $A$ is positive definite and hence can be decomposed as $A:=[A_{ij}]=\sigma\sigma^{T}$ (for instance by Cholesky factorization). For future reference, we point out the importance of this operator $\mathcal{L}$ in obtaining the (main) driving SDE, one can compare the coefficients of $\mathcal{L}$ with $\bf X$ in (2.8) below. It is assumed that $\varphi\leq 0$ and that the remaining coefficients and boundary in (2.3) are smooth enough (in particular, the outward normal vector ${\bf N}$ to $\partial\Omega_{R}$ is well defined) so that an unique solution exists [22], [34].

In order to introduce the stochastic representation of (2.3), let ${\bf X}_{t},\,0\leq t\leq T$ be a diffusion (see (2.8) below) starting at ${\bf X}_{0}={\bf x}_{0}\in\Omega\subset{\mathbb{R}}^{d},\,(d\geq 1)$, and consider the linear functional, in fact, a random variable,

[TABLE]

where

[TABLE]

In (2.2) and (2.6), ${\mathbbm{1}}_{\{H\}}$ is the indicator function ($1$ if $H$ is true and [math] otherwise); $\tau=\inf_{t}\{{\bf X}_{t}\in\partial\Omega_{A}\}$ is the “first exit (or first passage) time” from $\Omega$; which occurs at the “first exit point” ${\bf X}_{\tau}\in\partial\Omega_{A}$; $\xi_{t}$ is the “local time” (the amount of time the trajectory spends infinitesimally close to $\partial\Omega_{R}$); and all functions are assumed continuous and bounded. We want to calculate the expectation of $\phi$ in (2.2), which can also be expressed as

[TABLE]

where the processes $({\bf X}_{t},Y_{t},Z_{t},\xi_{t})$ are governed by a set of stochastic differential equations (SDEs) driven by a $d-$dimensional Wiener process ${\bf W}_{t}$:

[TABLE]

Then, formulas (2.2) and (2.7) are the solution of the parabolic linear BVP, i.e.

[TABLE]

The representation in terms of the SDE system (2.8), favoured by Milstein, is not the most usual one, but it can be programmed in a straightforward manner. If $c({\bf x})\leq 0$, elliptic BVPs can be formally derived from (2.3), the Feynman-Kac formulas are then known as Dynkin’s formulas. By taking $T\to\infty$, removing $p$, and dropping the time dependence from the surviving coefficients, the SDEs in (2.8) are now autonomous as long as $\tau$ is finite. (This fails to happen, for instance, with a purely reflecting boundary $(\partial\Omega=\partial\Omega_{R})$, when the solution $u({\bf x})$ is defined up to an arbitrary constant and requires one more compatibility condition [34].) In the case of purely Dirichlet BCs, $\xi_{\cdot}$ and the last equation in (2.8) drop out. The representation of linear BVPs can be simplified in many cases, such as (2.1) and (2.2) for (1.1).

2.3 An illustration of the effect of improved numerics

Remarkable research into the numerical methods to solve representations of linear BVPs given above during the last fifteen years has meant that they can be solved today at a fraction of the cost required when PDD was introduced. To highlight the effect of improved stochastic numerics on the PDD methodology, let us take an example from [13]. There the authors study the BVP

[TABLE]

with $f(x,y)$ such that u(x,y)=2\cos\big{(}2(y-2)x\big{)}+\sin\big{(}3(x-2)y\big{)}+3.1 is the exact solution, as well as the Dirichlet boundary condition. The BVP domain can be seen in Figure 2.1.

This problem was solved with the most current version of PDD (called IterPDD in [13]). IterPDD is a numerical suite which includes variance reduction techniques with iterative control variates based on nested, increasingly accurate global solutions.

The speedup of a method A over a method B to solve a BVP is defined as

[TABLE]

Table 2.1 shows the speedup of IterPDD over the previous versions of PDD. Note that IterPDD retains all the advantages of any previous version of the PDD algorithm and hence it is not necessary to compare it against DDD solvers, since $S(IterPDD,DDD)=S(IterPDD,PDD)\times S(PDD,DDD)$. The theoretical speedup in Table 2.1 is the optimal speedup predicted by a sensitivity algorithm presented in [13]. Importantly, the latter is designed in such a way that it relies on fast warm up Monte Carlo sampling, and runs in parallel so as not to defeat the purpose of PDD. The agreement with the experimentally observed speedup is consistently conservative.

The acceleration (speedup) grows larger as the target nodal accuracy $a_{0}$ (the largest admissible error of the numerical solution on the interface nodes, obtained via the probabilistic representation) decreases, approximately in an inversely proportional way. In summary, nearly every previously reported result using PDD, which were already faster than deterministic solvers, could be additionally accelerated by one to two orders of magnitude thanks to improved numerics. Further acceleration improvements are expected to follow in combination with Multilevel formulations [35].

Some applications of PDD

The initial undertaking of the PDD programme was due to J. A. Acebrón and his research group [2]-[6]. Some realistic large-scale simulations carried out in supercomputers at the Barcelona Supercomputer Center and Rome’s CASPUR proved the superiority of PDD over ScaLAPACK in terms of both total time and observed scalability when solving a nonlinear equations including KPP [8], [9], [5]. Independently, Gobet and Mairé proposed a PDD-like domain decomposition scheme for the Poisson equation in [39]. Bihlo and Haynes have applied PDD to the generation of meshes for finite elements, a PDE-based task very well suited to stochastic representations [16], [15], [17]. Moreover, in [6], the Vlasov-Poisson equations were tackled, which are of unquestionable interest in plasma physics. The PDD treatment of nonlinear PDEs like this one is further discussed in Section 3.

3 Branching Diffusions and the KPP equation

The PDEs covered in Section 2, are all linear, unfortunately though PDEs arising in applications are typically not linear and standard SDE arguments are insufficient to address the general settings. That being said, one can derive stochastic representations for more general PDEs using so-called Forward Backward SDEs (FBSDEs), but FBSDEs are computationally expensive and difficult to handle (see Section 4). However, [57], [42], [44] and [43] have further developed the “branching diffusions” methodology as an efficient method to solve wider classes of nonlinear PDEs than the originally proposed by McKean.

Branching diffusions allow one to tackle PDEs which are not linear without the FBSDE machinery. Although the classical results were somewhat restrictive on the form of the PDE (see (3.1) below), the recent developments allow one to address in an efficient way more general classes of PDEs. For example, when nonlinearities appear in the solution, or even gradients of the solution (see (3.3) & Section 3.3 below). Throughout this section we supplement the theory with examples of PDEs that are amenable to branching and whose representations open the door for PDD as a viable algorithm for a far larger class of problems than previously available.

In Section 2 we discussed stochastic representations for non-Dirichlet boundary conditions, unfortunately such boundary conditions have not been considered in this more general setting, but are future work for the authors. For the interested reader we mention recent analysis of the Monte Carlo Branching methodology to elliptic PDEs [10].

3.1 Reaction-Diffusion: KPP Equation

The branching diffusion idea is based on the works [58], [62], [52], to solve the KPP equation (in one spatial dimension with Dirichlet boundary condition),

[TABLE]

This was later generalized to allow for nonlinear PDEs such as,

[TABLE]

where $c$ is a positive constant. Solutions to such nonlinear PDEs can be represented through so-called branching diffusion processes. There has also been work for non integer powers of $u$, typically $u^{\alpha}$ for $\alpha\in[0,2]$, such process are referred to as super diffusions, although we will not discuss these here to focus on more recent work, see [27], [44] and references therein for further details. Again the SDE process is governed by $\mathcal{L}$, see (2.4).

Assumption 3.1** (Classical Branching Diffusions [52])**

In the classical setting the stochastic representation of the above PDE relies on the following two conditions for the coefficients $\alpha$; $\alpha_{i}\geq 0$ and $\sum_{i=0}^{\infty}\alpha_{i}=1$.

To explain the stochastic representation to (3.1), for point $(\mathbf{x},t)$ we need to consider a set of particles. We start a particle at $\mathbf{x}$ at time [math], this particle has a life time $\tau$ which is exponentially distributed with intensity $c$ (sometimes known as the “branching rate”), we then simulate this particle until $\min(t,\tau,\tau_{\partial\Omega})$, where $\tau_{\partial\Omega}$ is the particle hitting the spatial boundary, if $\tau$ is this smallest time, then the particle “branches” into $i$ particles with probability $\alpha_{i}$. These particles all start at the same point in space, however, they are all equipped with the own independent exponential random variable (“life time”) and driving Brownian motion. We then continue to simulate every particle until it has hit the space boundary or is alive at time $t$. We provide a detailed algorithm for branching diffusions (without space boundaries) in Section 3.1.1.

Consider the same initial and boundary value as above. The solution $u$ can then be written as the expected value of the product of the surviving particles or particles that hit the space boundary (see [8]),

[TABLE]

where $X^{(i)}_{s}$ is the location of the $i$th particle surviving at time $s$ and $N_{t}$ is the number of surviving particles at time $t$ (a particle that hits the boundary is still “alive” at time $t$).

Although this set up allows for more general PDEs, it is not difficult to see that the computational complexity is greater for these types of problems (though still smaller than that for FBSDEs). Moreover, because of the nature of the branching, we are mainly confined to cases which are not “overly” nonlinear with time horizons that are not too long. Otherwise, we may end up with “explosions” whereby we are required to keep track of an extremely large number of particles, which will also tend to increase the number of branchings. We postpone discussion of this to Section 3.2.

To help give an idea of branching diffusions we illustrate a sample path for a given PDE. Let us consider the PDE (one spatial dimension),

[TABLE]

From the work above there are two possible outcomes at a branching time (compare with (3.1)), either the particle dies, or splits into two descendants, the probability of each of these events is $\frac{1}{2}$. A possible trajectory (to compute $u(x_{0},T)$) is shown in Figure 3.1.

3.1.1 A PDD numerical example with the KPP Equation

Since we already have enough machinery to consider some interesting PDEs let us show an example which comes from the area of reaction diffusion equations, this class of equations are used in many areas such as physics and biology [32], [33], [49]. For illustration purposes we consider the most well known of these, the Kolmogorov-Petrovsky-Piskunov (KPP) equation [48], linked to applications in plasma physics and ecology. The general form of the KPP is the following nonlinear parabolic PDE,

[TABLE]

where $D$ and $r$ are constants. Following [8], taking $r=D=1$ and the initial value as,

[TABLE]

we can write the true solution of this problem as,

[TABLE]

The advantage of such an example is that we can also compare the error of each method. Following (3.2), the stochastic representation of the solution at $(x,t)$ is simply,

[TABLE]

where $N_{t}$ is the number of particles at time $t$, and $W^{(i)}_{t}$ is the Brownian motion for particle $i$ at time $t$. Since the operator $\mathcal{L}$ is $\partial^{2}/\partial x^{2}$, the process $X$ is Brownian motion started at point $x$ (compare with (2.3) and (2.8)).

Algorithm for the problem

We solve this problem222MATLAB was used for the implementation. The simulations ran on a Dell PowerEdge R430 with four intel xeon E5-2680 processors. All polynomial interpolation were carried out using MATLAB’s “polyfit” and the PDE is solved using MATLAB’s “pdepe” function. To keep this consistent with the error we expect from Monte Carlo, we have set the relative and absolute error in the solver as $10^{-3}$. over the domain $x\in[-2000,2000]$ for $t\in[0,1]$. For the PDD algorithm this corresponds to choosing points in space and solving them over a sequence of times to construct the artificial space time boundary, we create $p+1$ boundaries to obtain $p$ subdomains. Hence with five processors we split the domain in four subdomains $[-2000,-1000]$, $[-1000,0]$, $[0,1000]$ and $[1000,2000]$. This may not be the most optimal approach, but our goal is to show how even a simple PDD approach can dramatically improve the computational time required. We further calculate the solution at $11$ equally spaced time points (including $t=0$), which we denote by the set of $T_{i}$, satisfying $0=T_{0}<T_{1}<\dots<T=1$. Denote by $\Gamma$ the set of nodes at which we approximate the boundaries (true and artificial), hence for $D$ domain points, with $\Theta$ time points, $\Gamma$ is an $D$-by-$\Theta$ matrix. We denote by $\Gamma_{k}$ for $k\in\{1,\dots,D\}$ the column vectors of $\Gamma$, i.e. the approximation of the boundary at the $k$th domain point through time.

Remark 3.2** (Pruning - Controlling the growth of the branching)**

A useful technique to control the branching is the “pruning”, which is a method whereby we truncate the number of branches and comes from the observation that the pruned branches contribute very little to the expectation and are extremely computationally expensive. We do not discuss this further here, but direct the interested reader to [7].

It is clear from the PDE that all branching types will be two (the nonlinear part is a squared - compare with (3.1)). Moreover the driving process for the branching diffusions is a Brownian motion, this allows for us to take large step sizes since the simulation of Brownian motion is unbiased. Algorithm 3.1 provides a summary of the method

Error Analysis

Remark 3.3** (Straightforward implementation for Monte Carlo)**

For this example we used only a standard Monte Carlo implementation. Therefore, the results presented should be seen as the lower bound for what is achievable. It possible to improve the speed and accuracy of the algorithm by using methods touched on in Section 2.

Figure 3.2 compares the absolute error of Monte Carlo versus using only the PDE solver over the domain $[-5,5]$. We observe the maximum error of the the pure PDE solver is approximately $1/3$ that of PDD. However, the largest error for PDD is concentrated in one region to the left hand side of the boundary and outside this small region the two errors are comparable.

Scalability and running times of the algorithm

We conclude by considering how the time taken to solve the problem changes as more cores are used. We assume that we have access to the $p+1$ processors for $p$ subdomains. Since PDD consists of two main steps (Monte Carlo then PDE solver), we want to consider how these change as we increase the number cores. These results are presented in Figure 3.3. Note, to check the scaling we use the number of domains rather than the number of processors and PDD is not used when the number of domains is one.

Figure 3.3, clearly shows the scaling capabilities of PDD. As one increases the number of processors the Monte Carlo time stays close to constant (as expected) and the PDE solver time drops dramatically as the number of processors increases. In fact, as is presented on the right graph, this decrease is constant with the number of subdomains with no apparent sign of levelling off; this follows from the no interprocessor communication.

3.2 Marked Branching Diffusions

The idea of marked branching diffusions was developed by [42] as a solution to pricing credit valuation adjustments (CVAs) (see Section 3.2.2), but since then has been generalized by [44] and [43]; similar ideas appeared in [57]. Following [44], we consider a terminal value PDE of the form (recall the time-inversion argument in Section 2.2)

[TABLE]

the notation follows that in Section 3 with $L$ a positive integer.

We now discuss sufficient conditions for our stochastic representation to be the unique viscosity solution of this PDE. The technique in [44] is to use the fact that a large class of PDEs can be represented by FBSDEs, see [55]. The idea is to derive sufficient conditions for the FBSDE to be the unique bounded viscosity solution to a PDE of the form (3.3) and use this to obtain the results on branching diffusions; here, FBSDEs are used only as a theoretical method to enable the branching argument. The key assumption is the following [44, Assumption 2.2] denote the two functions for $s\in[0,\infty)$

[TABLE]

where $\|\cdot\|_{\infty}$ is the $L^{\infty}$-norm and the constants/functions follow from (3.3), then,

Assumption 3.4** (Assumptions for Marked branching)**

**

The power series $l_{0}$ has a radius of convergence $0<R\leq\infty$. Moreover, the function $l$ satisfies one of the following conditions:

1. (i)

$l(1)\leq 0$, 2. (ii)

$l(1)>0$* and for some $\hat{s}>1$, $l(s)>0$, $\forall s\in[1,\hat{s})$ and $l(\hat{s})=0$,* 3. (iii)

$l(s)>0$* $\forall s\in[1,\infty)$ and \int_{1}^{\bar{s}}\big{(}l(s)\big{)}^{-1}ds=T, for some constant $\bar{s}\in(1,{R}/{\|\psi\|_{\infty}})$.* 2. 2.

The terminal function satisfies $\|\psi\|_{\infty}<R$.

Under the above assumption, the stochastic equations associated to (3.3) have a unique bounded solution and yield the unique viscosity solution for PDE (3.3) (see Proposition 2.3 and Theorem 2.13 of [44]). Thus the above assumption provides sufficient conditions for the branching diffusions to be used.

We are now most of the way to stating the marked branching diffusion representation. Of course, $\alpha$ needs not be a probability distribution any more, so let us introduce a probability distribution $q$, where $q$ is chosen so $q_{i}>0$ if $\alpha_{i}\neq 0$. It will beneficial for us to introduce notation to keep track of the particles. Let $T_{n}$ the $n$th branching time of the system where at time $T_{n}$ one of the particles branches into $k$ particles, with probability $q_{k}$. We denote by $I_{n}$ the number of descendants created at time $T_{n}$, hence $I_{n}\in\{0,\dots,L\}$. Denote by $M_{T-t}:=\sup\{n:t+T_{n}\leq T\}$ the number of branches that occurred between time $t$ and $T$. Further denote by $\mathcal{K}_{t}$ the index of all particles alive at time $t$, therefore $\mathcal{K}_{T-t}$ corresponds to the index of all particles alive at time $T-t$. Finally, denote by $K_{n}$ the index of the particle that has branched at time $T_{n}$ (hence $K_{n}\in\mathcal{K}_{T_{n}}$). The representation for the solution of the PDE (3.3) at $(\mathbf{x},t)$ is given as

[TABLE]

Remark 3.5** (Variance of the estimator)**

It is also of note that the specific choice of probability distribution $q$ (subject to the conditions above) does not change the representation, however, from the view point of variance reduction an optimal choice exists, see [42].

Although the set up here has been carried out using exponential random variables as the lifetime of particles, more recent work suggest that other distributions yield smaller variances, see [43], [26], [61] and [20].

3.2.1 More general PDEs through approximation

The stochastic representations shown in this section can be used to deal with many different types of PDEs. Namely, these representations allow us to consider PDEs whose source term $f$ can be well approximated by polynomials. An example of this is considered in Section 3.2.2.

Although we mention CVAs, many other financial contracts rely on maximums between two entries (so called obstacle PDEs), these appear in American options for example. Computing the solution to certain obstacle PDEs in stochastics is deeply related to optimal stopping problems. We write such PDEs in their variational formulation: for a exercise payoff $\varphi$ the price function $u$ solves

[TABLE]

In [11], the authors showed how this PDE can be converted into a semilinear PDE

[TABLE]

A stochastic representation for these equations is well-understood and the methodology we have presented so far can be applied to them, see [42, Section 2.3]. Moreover, [18] use the polynomial representation of trigonometric functions to calculate the nonlinear Poisson-Boltzmann equation.

3.2.2 Numerical example with Option Pricing: Credit Valuation Adjustment (CVA)

The next example we consider is a different problem than the KPP above, but highlights the marked branching diffusion and its application scope very well. This example concerns option (derivative) pricing with Credit Valuation Adjustments (CVAs). Since the financial crisis in 2008, risk, especially default risk (the risk that a counter party fails to meet future obligations) has been at the forefront of many policy decisions and regulation changes for financial firms. CVAs play a crucial role here by adjusting the price of a derivative (a financial contract) with the knowledge that the counter party may default. [41] give a derivation of the PDEs arising in this problem. The specific form of the PDE depends on the way in which one chooses to model the mark-to-market value (the fair value) of the derivative at the time of default, but we obtain PDEs of the form,

[TABLE]

where the $r_{i}$ are functions of $x$ and $t$ and we denote by $y^{+}:=\max(0,y)$.

Although at first glance such a PDE seems beyond the scope of the stochastic representations in Section 3.2, [42] provided a simple yet powerful trick to deal with such PDEs (see also [57]). For simplicity we will consider only the one dimensional case. The first trick, is to rescale the problem, by assuming the payoff is bounded. Then we may consider the following function $v:=u/||\psi||_{\infty}$, such a function satisfies the PDE

[TABLE]

with terminal value bounded above by $1$. Let us consider instead a PDE of the form

[TABLE]

where $F$ is a polynomial, hence this equation is of the type considered previously. The great insight in [42] was to use a polynomial to approximate the semilinear component $v^{+}$. Therefore we have transformed the PDE from one outside the capabilities of marked branching diffusions, into one which can be. [42] uses a polynomial of order $4$ to do this, which provides a good approximation as can be seen in Fig. 3.4.

Thus the PDE we want to calculate is,

[TABLE]

Since we have constant coefficients and a bounded terminal condition, this PDE easily satisfies Assumption 3.4, hence we can use the stochastic representation. Remark the negative coefficient associated to the 4th-power term; the classical Assumption 3.1 does not hold here and hence the motivation behind Marked Branching Diffusions.

Of course, there does not seem much scope for PDD here since the domain is rather small. However, in higher dimensional problems (basket options), the PDD method can very easily provide the computational advantages highlighted in earlier sections.

3.3 Further Generalisations: Aged marked branching diffusions

In [43], the authors generalize the representation to a wider class of semilinear PDEs through so called age-marked branching diffusions, still building on the paradigm of semilinear PDEs with polynomial source terms $f$. Consider for some integer $m\geq 0$ a set $L\subset\mathbb{N}^{m+1}$ and a sequence of functions $(c_{l})_{l\in L}$ and $(h_{i})_{i=1,\cdots,m}$, where $c_{l}:[0,T]\times\mathbb{R}^{d}\rightarrow\mathbb{R}$ and $h_{i}:[0,T]\times\mathbb{R}^{d}\rightarrow\mathbb{R}$. Let the source term function $f:[0,T]\times\mathbb{R}^{d}\times\mathbb{R}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$ be of the form

[TABLE]

With such a function in mind (under some regularity assumptions) [43] provides a stochastic representation for PDEs of the form,

[TABLE]

where $\psi:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is a bounded Lipschitz function. The assumptions required for a stochastic representation are more technical than the previous case considered. We thus omit them here and point the reader to [43] for the details. Under this set up we can handle systems of PDE and also popular PDEs such as Burger’s equation. In fact, using the polynomial approximation trick discussed in Section 3.2.1, we can handle PDEs containing $f(\cdot,u,\nabla u)=u|\nabla u|^{2}$ terms, which appear in many applications for instance, to construct harmonic homotopic maps [60]; in the theory of ferromagnetic materials through the Landau-Lifschitz-Gilbert equation or models of magnetostriction.

Remark 3.6** (Fully nonlinear case)**

Recently, [61] constructed an algorithm from techniques developed in [43] to approximate a fully non-linear PDEs. This work is still experimental (at the present time) and thus we only reference it here for interested readers.

3.3.1 Some other PDE examples in data science solvable by branching diffusions

We focus on branching diffusions, since they are amenable to PDD. In mathematical biology one often wants to consider the affect on a population of a species given the presence of another species and the interaction between them. Such models have been considered in works such as [25], see [38] as well. In the paper by Escher and Matioc [30] they present a model to describe the growth of tumours. The model is a nonlinear two dimensional system with two decoupled Dirichlet problems.

As highlighted in Section 3.2.1, the methodology described in this work can also be used to deal with PDE obstacle problems. Apart from the well-known problem of pricing American type options we point out to the so-called Travel agency problem [38, Section 5.2] where a travel agency needs to decide when to offer travels depending on currency and weather forecast.

If one returns to finance, the approximation method highlighted in Section 3.2.1 combined with the age-marked branching diffusions can be used to tackle the problems in actuarial science, namely the “Reinsurance and Variable Annuity” problem (see [41]) or the pricing of contracts with different borrowing/lending rates [23, Section 3.3].

4 An Invitation to More Challenging Problems

In this Section, we shortly review probabilistic representations of nonlinear problems based on Forward-Backward SDE (FBSDEs) for which branching arises as a particular case. FBSDEs provide stochastic representations for a wide class of PDEs, wider than branching techniques allow, including systems of (fully nonlinear, parabolic or elliptic) PDEs (or BVPs) combining any type of mixed boundary condition.

Consider the system of $K\geq 1$ quasilinear coupled PDEs in ${\mathbb{R}}^{d}\times[0,T]$, $d\geq 1$:

[TABLE]

The associated SDE to $\mathcal{L}$, ${\bf X}_{t}$, is as in equation (2.8) (with the same requirements on the coefficients as there). Under the proper conditions on the coefficients in (4.1) and adequate smoothness of the $u^{(1)}({\bf x},t),\ldots,u^{(K)}({\bf x},t)$ the pointwise solution of the system at $({\bf x}_{0},0)$ is given by $(1\leq k\leq K)$

[TABLE]

where the processes $Y^{(1)}_{t},\ldots,Y^{(K)}_{t},{\bf Z}_{t}^{(1)},\ldots,{\bf Z}_{t}^{(K)}$ are the solution to the the system of FBSDEs (see [37, chapter VII])

[TABLE]

One can show via Itô’s formula that $Y_{t}^{(k)}=u^{(k)}({\bf X}_{t},t)$ and ${\bf Z}^{(k)}_{t}=\nabla u^{(k)}({\bf X}_{t},t)\sigma({\bf X}_{t},t)$, $(1\leq k\leq K)$ — this lifts the apparent inconsistency of having more processes than equations in (4.3). As with linear parabolic equations, the terminal condition can be transformed into an initial one by reversing time, i.e. by letting ${\hat{t}}\shortrightarrow T-t$ and $\partial/\partial{\hat{t}}\shortrightarrow-\partial/\partial t$. System (4.1) accommodates many important applications in biochemistry, stochastic control, finance and physics [37, chapter VII]. Multiple extensions of the above system are possible, including different generators for each of the equations in the system (i.e. different ${\cal L}^{(1)},\ldots,{\cal L}^{(K)}$), as well as BCs of various types. The complexity of the probabilistic representation grows accordingly.

In order to simplify the discussion and focus on difficulties, let $K=1$ in the remainder of this Section (thus $Y_{t}=Y^{(1)}_{t}$, $\psi=\psi^{(1)}$, $f=f^{(1)}$), and $f=f(t,{\bf x},u)$. It can be shown

[TABLE]

The manipulation of the conditional expectation (4.4) requires the introduction of filtered probability spaces, which lies beyond the scope of this paper. The numerical treatment is equally delicate. After a canonical time discretization over a uniform time partition with step $h=T/N$, we have

[TABLE]

Here we emphasize that to compute $u({\bf x},0)=Y_{0}$ one needs, at each time step $t_{i}$, to compute an approximation of $u({\bf x},t_{i})$ for $x$ over the domain of the PDE. Comparing with the result when branching is possible, see equation (3.2), the gain of branching over the general FBSDE machinery is clear. Nonetheless, branching techniques have their limits and the general case must, in principle, be tackled through FBSDEs and its variants. How to carry out a PDD methodology within the framework of FBSDEs that still retains PDD’s characteristic scalability is an open question.

In general, the timestep-per-timestep backward iteration and the computation and/or approximation of the conditional expectation in (4.5) can be done in many ways, for instance, by resorting to Bellman’s dynamic programming principle, whereby the solution is projected on a functional basis by regression (see [37, chapter VIII]). Overviews on numerics of FBSDE can be found in [19], [23] and PDE inspired problems dealt by FBSDEs can be found in [33], [49], [50].

A brief remark on FBSDEs and stochastic representations

We close this section with comments on various probabilistic representations not mentioned elsewhere in this paper.

There is theoretical work relating FBSDEs to the Navier-Stokes equations [24]. While, in general, transport equations do not have a stochastic representation, there are some important exceptions, of which we cite two: the one-dimensional telegraph equation [4]; and the the Vlasov-Poisson system of equations. For the latter, it was proposed in [53] a probabilistic representation in the Fourier space.

In [18] the authors deal with representations (and Monte Carlo simulation) for nonlinear divergence-form elliptic Poisson-Boltzmann PDE over the whole $\mathbb{R}^{3}$. Quasilinear parabolic PDEs (multidimensional and systems) admit a stochastic representation in terms of FBSDEs [56], [55] and many extensions exist ranging from representations for obstacle PDE problems [28] to representations of fully nonlinear elliptic and parabolic PDEs [21]. See as well [33], [49] for stochastic representations for reaction-diffusion PDEs, PDEs with quadratic gradients terms and of Burger’s type nonlinearities. Complex valued PDEs and the FBSDE connection are handled in [63].

Systems of fully nonlinear systems of second-order PDEs (i.e. including possible nonlinearities on the highest derivatives i.e. $f$ can depend on “$\Delta u$” terms) have recently found a representation in terms so-called 2BSDEs [21]; numerical schemes have already been proposed [31]. Lastly, much of the theory (on stochastic representations) can be extended to frameworks where the coefficients $b$ and $\sigma$ depend on $u$ and $\nabla u$ (see [51]).

Acknowledgements.

F. Bernal acknowledges funding from Centre de Mathématiques Appliquées (CMAP), École Polytechnique. G. dos Reis gratefully thanks the partial support by the Fundação para a Ciência e a Tecnologia (Portuguese Foundation for Science and Technology) through the project UID/MAT/00297/2013 (Centro de Matemática e Aplicações). G. Smith was supported by The Maxwell Institute Graduate School in Analysis and its Applications, a Centre for Doctoral Training funded by the UK Engineering and Physical Sciences Research Council (grant [EP/L016508/01]), the Scottish Funding Council, Heriot-Watt University and the University of Edinburgh.