Convergence of the Deep BSDE Method for Coupled FBSDEs

Introduction

Forward-backward stochastic differential equations (FBSDEs) and partial differential equations (PDEs) of parabolic type have found numerous applications in stochastic control, finance, physics, etc., as a ubiquitous modeling tool. In most situations encountered in practice the equations cannot be solved analytically but require certain numerical algorithms to provide approximate solutions. On the one hand, the dominant choices of numerical algorithms for PDEs are mesh-based methods, such as finite differences, finite elements, etc. On the other hand, FBSDEs can be tackled directly through probabilistic means, with appropriate methods for the approximation of conditional expectation. Since these two kinds of equations are intimately connected through the nonlinear Feynman–Kac formula , the algorithms designed for one kind of equation can often be used to solve another one.

However, the aforementioned numerical algorithms become more and more difficult, if not impossible, when the dimension increases. They are doomed to run into the so-called “curse of dimensionality” when the dimension is high, namely, the computational complexity grows exponentially as the dimension grows. The classical mesh-based algorithms for PDEs require a mesh of size $O(N^{d})$. The simulation of FBSDEs faces a similar difficulty in the general nonlinear cases, due to the need to compute conditional expectation in high dimension. The conventional methods, including the least squares regression , Malliavin approach , and kernel regression , are all of exponential complexity. There are a limited number of cases where practical high-dimensional algorithms are available. For example, in the linear case, Feynman–Kac formula and Monte Carlo simulation together provide an efficient approach to solving PDEs and associated BSDEs numerically. In addition, methods based on the branching diffusion process and multilevel Picard iteration overcome the curse of dimensionality in their considered settings. We refer for the detailed discussion on the complexity of the algorithms mentioned above. Overall there is no numerical algorithm in literature so far proved to overcome the curse of dimensionality for general quasilinear parabolic PDEs and the corresponding FBSDEs.

A recently developed algorithm, called the deep BSDE method , has shown astonishing power in solving general high-dimensional FBSDEs and parabolic PDEs . In contrast to conventional methods, the deep BSDE method employs neural networks to approximate unknown gradients and reformulates the original equation-solving problem into a stochastic optimization problem. Thanks to the universal approximation capability and parsimonious parameterization of neural networks, in practice the objective function can be effectively optimized in high-dimensional cases, and the function values of interests are obtained quite accurately.

The deep BSDE method was initially proposed for decoupled FBSDEs. In this paper, we extend the method to deal with coupled FBSDEs and a broader class of quasilinear parabolic PDEs. Furthermore, we present an error analysis of the proposed scheme, including decoupled FBSDEs as a special case. Our theoretical result consists of two theorems. Theorem 1 provides a posteriori error estimation of the deep BSDE method. As long as the objective function is optimized to be close to zero under fine time discretization, the approximate solution is close to the true solution. In other words, in practice, the accuracy of the numerical solution is effectively indicated by the value of the objective function. Theorem 2 shows that such a situation is attainable, by relating the infimum of the objective function to the expression ability of neural networks. As an implication of the universal approximation property (in the $L^{2}$ sense), there exist neural networks with suitable parameters such that the obtained numerical solution is approximately accurate. To the best of our knowledge, this is the first theoretical result of the deep BSDE method for solving FBSDEs and parabolic PDEs. Although our numerical algorithm is based on neural networks, the theoretical result provided here is equally applicable to the algorithms based on other forms of function approximations.

The article is organized as follows. In section 2, we precisely state our numerical scheme for coupled FBSDEs and quasilinear parabolic PDEs and give the main theoretical results of the proposed numerical scheme. In section 3, the basic assumptions and some useful results from the literature are given for later use. The proofs of the two main theorems are provided in section 4 and section 5, respectively. Some numerical experiments with the proposed scheme are presented in section 6.

A Numerical Scheme for Coupled FBSDEs and Main Results

Without loss of clarity, here we use the notation $X_{0}^{\pi}$ as $X_{t_{0}}^{\pi}$, $X_{T}^{\pi}$ as $X_{t_{N}}^{\pi}$, etc.

Following the spirit of the deep BSDE method, we employ a stochastic optimizer to solve the following stochastic optimization problem

where $Y_{0}$ is $\mathcal{F}_{0}$-measurable and square-integrable, and $Z_{t}$ is a $\mathcal{F}_{t}$-adapted square-integrable process. The solution of the FBSDEs (2.1)(2.2) is a minimizer of the above problem since the loss function attains zero when it is evaluated at the solution. In addition, the wellposedness of the FBSDEs (under some regularity conditions) ensures the existence and uniqueness of the minimizer. Therefore, we expect (2.4), as a discretized counterpart of (2.5), defines a benign optimization problem and the associated near-optimal solution provides us a good approximate solution of the original FBSDEs. The reason we do not represent $Z_{t_{i}}$ as a function of $X_{t_{i}}$ only is that the process $\{X_{t_{i}}^{\pi}\}_{0\leq i\leq N}$ is not Markovian, while the process $\{X_{t_{i}}^{\pi},Y_{t_{i}}^{\pi}\}_{0\leq i\leq N}$ is Markovian, which facilitates our analysis considerably. If $b$ and $\sigma$ are both independent of $Y$, then the FBSDEs (2.1)(2.2) are decoupled, we can take $\phi_{i}^{\pi}$ as a function of $X_{t_{i}}^{\pi}$ only, as the numerical scheme introduced in .

Our two main theorems regarding the deep BSDE method are the following, mainly on the justification and property of the objective function (2.4) in the general coupled case, regardless of the specific choice of parametric function spaces. An important assumption for the two theorems is the so-called weak coupling or monotonicity condition, which will be explained in detail in section 3. The precise statement of the theorems can be found in Theorem 1′ (section 4) and Theorem 2′ (section 5), respectively.

Under some assumptions, there exists a constant C, independent of h, d, and m, such that for sufficiently small h,

where $\hat{X}_{t}^{\pi}=X_{t_{i}}^{\pi}$, $\hat{Y}_{t}^{\pi}=Y_{t_{i}}^{\pi}$, $\hat{Z}_{t}^{\pi}=Z_{t_{i}}^{\pi}$ for $t\in[t_{i},t_{i+1})$.

Under some assumptions, there exists a constant C, independent of h, d and m, such that for sufficiently small h,

Briefly speaking, Theorem 1 states that the simulation error (left side of equation (2.6)) can be bounded through the value of the objective function (2.4). To the best of our knowledge, this is the first result for the error estimation of the coupled FBSDEs, concerning both time discretization error and terminal distance. Theorem 2 states that the optimal value of the objective function can be small if the approximation capability of the parametric function spaces ($\mathcal{N}^{\prime}_{0}$ and $\mathcal{N}_{i}$ above) is high. Neural networks are a promising candidate for such a requirement, especially in high-dimensional problems. There are numerous results, dating back to the 90s (see, e.g., ), in regard to the universal approximation and complexity of neural networks. There are also some recent analysis on approximating the solutions of certain parabolic partial differential equations with neural networks. However, the problem is still far from resolved. Theorem 2 implies that if the involved conditional expectations can be approximated by neural networks whose numbers of parameters growing at most polynomially both in the dimension and the reciprocal of the required accuracy, then the solutions of the considered FBSDEs can be represented in practice without the curse of dimensionality. Under what conditions this assumption is true is beyond the scope of this work and remains for further investigation.

The above-mentioned scheme in (2.3)(2.4) is for solving FBSDEs. The so-called nonlinear Feynman–Kac formula, connecting FBSDEs with the quasilinear parabolic PDEs, provides an approach to numerically solve quasilinear parabolic PDEs (2.7) below through the same scheme. We recall a concrete version of the nonlinear Feynman–Kac formula in Theorem 3 below and refer interested readers to e.g., for more details. According to this formula, the term $E|Y_{0}-Y_{0}^{\pi}|^{2}$ can be interpreted as $E|u(0,\xi)-\mu_{0}^{\pi}(\xi)|^{2}$. Therefore, we can choose the random variable $\xi$ with a delta distribution, a uniform distribution in a bounded region, or any other distribution we are interested in. After solving the optimization problem, we obtain $\mu_{0}^{\pi}(\xi)$ as an approximation of $u(0,\xi)$. See for more details.

$m=d$ and $b(t,x,y)$, $\sigma(t,x,y)$, $f(t,x,y,z)$ are smooth functions with bounded first-order derivatives with respect to $x,y,z$.

There exist a positive continuous function $\nu$ and a constant $\mu$, satisfying that

Then the following quasilinear PDE has a unique classical solution $u(t,x)$ that is bounded with bounded $u_{t}$, $\nabla_{x}u$, and $\nabla^{2}_{x}u$,

The associated FBSDEs (2.1)(2.2) have a unique solution $(X_{t},Y_{t},Z_{t})$ with $Y_{t}=u(t,X_{t})$, $Z_{t}=\sigma^{\operatorname{T}}(t,X_{t},u(t,X_{t}))\nabla_{x}u(t,X_{t})$, and $X_{t}$ is the solution of the following SDE

The statement regarding FBSDEs (2.1)(2.2) in Theorem 3 is developed through a PDE-based argument, which requires $m=d$, uniform ellipticity of $\sigma$, and high-order smoothness of $b,\sigma,f$, and $g$. An analogous result through probabilistic argument is given below in Theorem 4 (point 4). In that case, we only need the Lipschitz condition for all of the involved functions, in addition to some weak coupling or monotonicity conditions demonstrated in Assumption 3. Note that the Lipschitz condition alone does not guarantee the existence of a solution to the coupled FBSDEs, even in the situation when $b,f,\sigma$ are linear (see for a concrete counterexample).

Theorem 3 also implies that the assumption that the drift function $b$ only depends on $x,y$ is general. If $b$ depends on $z$ as well, one can move the associated term in (2.7) into the nonlinearity $f$ and apply the nonlinear Feynman–Kac formula back to obtain an equivalent system of coupled FBSDEs, in which the new drift function is independent of $z$.

Preliminaries

In this section, we introduce our assumptions and two useful results in . We use the notation $\Delta x=x_{1}-x_{2}$, $\Delta y=y_{1}-y_{2}$, $\Delta z=z_{1}-z_{2}$.

There exist (possibly negative) constants $k_{b}$, $k_{f}$ such that

b, $\sigma$, f, g are uniformly Lipschitz continuous with respect to (x,y,z). In particular, there are non-negative constants K, $b_{y}$, $\sigma_{x}$, $\sigma_{y}$, $f_{x}$, $f_{z}$, and $g_{x}$ such that

$b(t,0,0)$, $f(t,0,0,0)$, and $\sigma(t,0,0)$ are bounded. In particular, there are constants $b_{0}$, $\sigma_{0}$, $f_{0}$, and $g_{0}$ such that

We note here $b_{y}$ et al. are all constants, not partial derivatives. For convenience, we use $\mathscr{L}$ to denote the set of all the constants mentioned above and assume K is the upper bound of $\mathscr{L}$.

$b,\sigma,f$ are uniformly Hölder-$\frac{1}{2}$ continuous with respect to $t$. We assume the same constant K to be the upper bound of the square of the Hölder constants as well.

Small time duration, that is, T is small.

Weak coupling of Y into the forward SDE (2.1), that is, $b_{y}$ and $\sigma_{y}$ are small. In particular, if $b_{y}=\sigma_{y}=0$, then the forward equation does not depend on the backward one and, thus, equations (2.1)(2.2) are decoupled.

Weak coupling of X into the backward SDE (2.2), that is, $f_{x}$ and $g_{x}$ are small. In particular, if $f_{x}=g_{x}=0$, then the backward equation does not depend on the forward one and, thus, equations (2.1)(2.2) are also decoupled. In fact, in this case, Z = 0 and (2.2) reduces to an ODE.

f is strongly decreasing in y, that is, $k_{f}$ is very negative.

b is strongly decreasing in x, that is, $k_{b}$ is very negative.

The assumptions stated above are usually called weak coupling and monotonicity conditions in literature . To make it more precise, we define

Then, a specific quantitative form of the above five conditions can be summarized as:

In other words, if any of the five conditions of the weak coupling and monotonicity conditions holds to certain extent, the two inequalities in (3.1) hold. Below, we refer to (3.1) as Assumption 3 and the five general qualitative conditions described above as the weak coupling and monotonicity conditions.

The above three assumptions are basic assumptions in , which we need in order to use the results from , as stated in Theorems 4 and 5 below. Theorem 4 gives the connections between coupled FBSDEs and quasilinear parabolic PDEs under weaker conditions. Theorem 5 provides the convergence of the implicit scheme for coupled FBSDEs. Our work primarily uses the same set of assumptions except that we assume some further quantitative restrictions related to the weak coupling and monotonicity conditions, which will be transparent through the extra constants we define in proofs. Our aim is to provide explicit conditions on which our results hold and more clearly present the relationship between these constants and the error estimates. As will be seen in the proof, roughly speaking, the weaker the coupling (resp., the stronger the monotonicity, the smaller the time horizon) is, the easier the condition is satisfied, and the smaller the constant $C$ related with error estimates are.

$|u(t,x_{1})-u(t,x_{2})|^{2}\leq L_{1}|x_{1}-x_{2}|^{2}$.

$|u(s,x)-u(t,x)|^{2}\leq C(1+|x|^{2})|s-t|$ with some constant C depending on $\mathscr{L}$ and $T$.

u is a viscosity solution of the PDE (2.7).

The FBSDEs (2.1)(2.2) have a unique solution $(X_{t},Y_{t},Z_{t})$ and $Y_{t}=u(t,X_{t})$. Thus, $(X_{t},Y_{t},Z_{t})$ satisfies decoupled FBSDEs

Furthermore, the solution of the FBSDEs satisfies the path regularity with some constant C depending on $\mathscr{L}$ and T

Several conditions can guarantee $Z_{t}$ admits a càdlàg version, such as $m=d$ and $\sigma\sigma^{\operatorname{T}}\geq\delta I$ with some $\delta>0$, see e.g., .

Under Assumptions 1, 2, and 3, for sufficiently small h, the following discrete-time equation ($0\leq i\leq N-1$)

where $\overline{X}_{t}^{\pi}=\overline{X}_{t_{i}}^{\pi}$, $\overline{Y}_{t}^{\pi}=\overline{Y}_{t_{i}}^{\pi}$, $\overline{Z}_{t}^{\pi}=\overline{Z}_{t_{i}}^{\pi}$ for $t\in[t_{i},t_{i+1})$, and C is a constant depending on $\mathscr{L}$ and T.

In , the above result (existence and convergence) is proved for the explicit scheme, which is formulated as replacing $f(t_{i},\overline{X}_{t_{i}}^{\pi},\overline{Y}_{t_{i}}^{\pi},\overline{Z}_{t_{i}}^{\pi})$ with $f(t_{i},\overline{X}_{t_{i}}^{\pi},\overline{Y}_{t_{i+1}}^{\pi},\overline{Z}_{t_{i}}^{\pi})$ in the last equation of (3.3). The same techniques can be used to prove the implicit scheme, as we state in Theorem 5.

Finally, to make sure the system in (2.3) is well-defined, we restrict our parametric function spaces $\mathcal{N}^{\prime}_{0}$ and $\mathcal{N}_{i}$ as in Assumption 4 below. Note that neural networks with common activation functions, including ReLU and sigmoid function, satisfy this assumption. Under Assumption 1 and 4, one can easily prove by induction that $\{X_{t_{i}}^{\pi}\}_{0\leq i\leq N}$, $\{Y_{t_{i}}^{\pi}\}_{0\leq i\leq N}$ and $\{Z_{t_{i}}^{\pi}\}_{0\leq i\leq N-1}$ defined in (2.3) are all measurable and square-integrable random variables.

A Posteriori Estimation of the Simulation Error

We prove Theorem 1 in this section. Comparing the statements of Theorem 1 and Theorem 5, we wish to bound the differences between $(X_{t_{i}}^{\pi},Y_{t_{i}}^{\pi},Z_{t_{i}}^{\pi})$ and $(\overline{X}_{t_{i}}^{\pi},\overline{Y}_{t_{i}}^{\pi},\overline{Z}_{t_{i}}^{\pi})$ with the objective function $E|g(X_{T}^{\pi})-Y_{T}^{\pi}|^{2}$. Recalling the definition of the system of equations (2.3), we have

Taking the expectation $E[\cdot|\mathcal{F}_{t_{i}}]$ on both sides of (4.2), we obtain

Right multiplying $(\Delta W_{i})^{\operatorname{T}}$ on both sides of (4.2) and taking the expectation $E[\cdot|\mathcal{F}_{t_{i}}]$ again, we obtain

The above observation motivates us to consider the following system of equations

Note that (4.3) is defined just like the FBSDEs (2.1)(2.2), where the $X$ component is defined forwardly and the $Y,Z$ components are defined backwardly. However, since we do not specify the terminal condition of $Y_{T}^{\pi}$, the system of equations (4.3) has infinitely many solutions. The following lemma gives an estimate of the difference between two such solutions.

Let $\delta X_{i}=X_{t_{i}}^{\pi,1}-X_{t_{i}}^{\pi,2},\delta Y_{i}=Y_{t_{i}}^{\pi,1}-Y_{t_{i}}^{\pi,2}$, then we have, for $0\leq n\leq N$,

To prove Lemma 1, we need the following lemma to handle the $Z$ component.

By the martingale representation theorem, there exists an $\mathscr{F}_{t}$-adapted square-integrable process $\{\delta Z_{t}\}_{t_{i}\leq t\leq t_{i+1}}$ such that

From equation (4.9), by Assumptions 1, 2 and the root-mean square and geometric mean inequality (RMS-GM inequality), for any $\lambda_{1}>0$, we have

Recall $A_{1}=2k_{b}+\lambda_{1}+\sigma_{x}+Kh$, $A_{2}=(\lambda_{1}^{-1}+h)b_{y}+\sigma_{y}$, $E|\delta X_{0}|^{2}=0$. By induction we can obtain that, for $0\leq n\leq N$,

Similarly, from equation (4.10), for any $\lambda_{2}>0$, we have

To deal with the integral term in (4.11), we apply Lemma 2 to (4.6)(4.8) and get

where $(\cdot)_{k}$ denotes the $k$-th component of the vector. Plugging it into (4.11) gives us

Then for any $\lambda_{2}\geq f_{z}$ and sufficiently small $h$ satisfying $(2k_{f}+\lambda_{2})h<1$, we have

Recall $A_{3}=-h^{-1}\ln[1-(2k_{f}+\lambda_{2})h]$, $A_{4}=f_{x}\lambda_{2}^{-1}[1-(2k_{f}+\lambda_{2})h]^{-1}$. By induction we obtain that, for $0\leq n\leq N$,

Now we are ready to prove Theorem 1, whose precise statement is given below. Note that its conditions are satisfied if any of the five cases in the weak coupling and monotonicity conditions holds.

Suppose Assumptions 1, 2, 3, and 4 hold true and there exist $\lambda_{1}>0,\lambda_{2}\geq f_{z}$ such that $\overline{A_{0}}<1$, where

Then there exists a constant $C>0$, depending on $E|\xi|^{2}$, $\mathscr{L}$, $T$, $\lambda_{1}$, and $\lambda_{2}$, such that for sufficiently small $h$,

where $\hat{X}_{t}^{\pi}=X_{t_{i}}^{\pi}$, $\hat{Y}_{t}^{\pi}=Y_{t_{i}}^{\pi}$, $\hat{Z}_{t}^{\pi}=Z_{t_{i}}^{\pi}$ for $t\in[t_{i},t_{i+1})$.

The above theorem also implies the coercivity of the objective function (2.4) used in the deep BSDE method. Formally speaking, the coercivity means that if $\sum_{i=0}^{N-1}E|Z_{t_{i}}^{\pi}|^{2}+E|Y_{0}^{\pi}|^{2}\rightarrow+\infty$, we have $E|g(X_{T}^{\pi})-Y_{T}^{\pi}|^{2}\rightarrow+\infty$, which is a direct result from Theorem 1′.

If any of the weak coupling and monotonicity conditions introduced in Assumption 3 holds to a sufficient extent, there must exist $\lambda_{1},\lambda_{2}$ satisfying the conditions in Theorem 1′. We discuss the 5 cases in what follows.

Suppose all other constants and $\lambda_{1}>0,\lambda_{2}\geq f_{z}$ are fixed, if $T>0$ is sufficiently small, then the second factor of $\overline{A_{0}}$ could be sufficiently close to 0 such that $\overline{A_{0}}<1$.

Suppose all other constants and $\lambda_{1}>0,\lambda_{2}\geq f_{z}$ are fixed, if $b_{y}\geq 0$ and $\sigma_{y}\geq 0$ are sufficiently small, then $\overline{A_{2}}\geq 0$ could be sufficiently small such that $\overline{A_{0}}<1$.

Suppose all other constants and $\lambda_{1}>0,\lambda_{2}\geq f_{z}$ are fixed, if $f_{x}\geq 0$ and $g_{x}\geq 0$ are sufficiently small, then $\overline{A_{4}}$ and thus the last factor in $\overline{A_{0}}$ could be sufficiently close to 0 such that $\overline{A_{0}}<1$.

Suppose all constants except $k_{f}$ and $\lambda_{2}>0$ are fixed. Let $\overline{A_{1}}^{\prime}\coloneqq\overline{A_{1}}+\overline{A_{3}}=2k_{b}+2k_{f}+\sigma_{x}+\lambda_{1}+\lambda_{2}$ and rewrite $\overline{A_{0}}$ as

It is straightforward to check that there exists a negative constant $C_{1}$ such that when $\overline{A_{1}}^{\prime}\leq C_{1}$, $(e^{\overline{A_{1}}^{\prime}T}-1)/\overline{A_{1}}^{\prime}<1/(2\overline{A_{2}}g_{x})$. By the definition of $\overline{A_{1}}^{\prime}$, if $k_{f}$ is sufficiently negative, there exists $\lambda_{2}\geq f_{x}$ such that $\overline{A_{1}}^{\prime}=C_{1}$ and $\lambda_{2}$ is sufficiently large to ensure

Combining these two estimates gives $\overline{A_{0}}<1$.

Noting that $k_{b}$ and $k_{f}$ play the same role in $\overline{A_{1}}^{\prime}$, we use the same argument as above to show that when $k_{b}$ is sufficiently negative, there exists $\lambda_{2}\geq f_{x}$ such that $\overline{A_{0}}<1$.

From the proof of this theorem and throughout the remainder of the paper, we use $C$ to generally denote a constant that only depends on $E|\xi|^{2}$, $\mathscr{L}$, and $T$, whose value may change from line to line when there is no need to distinguish. We also use $C(\cdot)$ to generally denote a constant depending on $E|\xi|^{2}$, $\mathscr{L}$, $T$ and the constants represented by $\cdot$.

We use the same notations as Lemma 1. Let $X_{t_{i}}^{\pi,1}=X_{t_{i}}^{\pi}$, $Y_{t_{i}}^{\pi,1}=Y_{t_{i}}^{\pi},Z_{t_{i}}^{\pi,1}=Z_{t_{i}}^{\pi}$ (defined in system (2.3)) and $X_{t_{i}}^{\pi,2}=\overline{X}_{t_{i}}^{\pi}$, $Y_{t_{i}}^{\pi,2}=\overline{Y}_{t_{i}}^{\pi}$, $Z_{t_{i}}^{\pi,2}=\overline{Z}_{t_{i}}^{\pi}$ (defined in system (3.3)). It can be easily checked that both $(\{X_{t_{i}}^{\pi,j}\}_{0\leq i\leq N}$, $\{Y_{t_{i}}^{\pi,j}\}_{0\leq i\leq N}$, $\{Z_{t_{i}}^{\pi,j}\}_{0\leq i\leq N-1})$, $j=1,2$ satisfy the system of equations (4.3). Our proof strategy is to use Lemma 1 to bound the difference between two solutions through the objective function $E|g(X_{T}^{\pi})-Y_{T}^{\pi}|^{2}$. This allows us to apply Theorem 5 to derive the desired estimates.

To begin with, note that for any $\lambda_{3}>0$, the RMS-GM inequality yields

Therefore by definition of $P$ and $S$, we have

When $\overline{A_{0}}<1$, comparing $\lim_{h\rightarrow 0}A(h)$ and $\overline{A_{0}}$, we know that, for any $\epsilon>0$, there exists $\lambda_{3}>0$ and sufficiently small $h$ such that

By fixing $\epsilon=1$ and choosing suitable $\lambda_{3}$, we obtain our error estimates of $E|\delta X_{n}|^{2}$ and $E|\delta Y_{n}|^{2}$ as

To estimate $E|\delta Z_{n}|^{2}$, we consider estimate (4.12), in which $\lambda_{2}$ can take any value no smaller than $f_{z}$. If $f_{z}\neq 0$, we choose $\lambda_{2}=2f_{z}$ and obtain

The case $f_{z}=0$ can be dealt with similarly by choosing $\lambda_{2}=1$ and the same type of estimate can be derived. Finally, combining estimates (4.18)(4.19)(4.20) with Theorem 5, we prove the statement in Theorem 1′. ∎

An Upper Bound for the Minimized Objective Function

We prove Theorem 2 in this section. We first state three useful lemmas. Theorem 2′, as a detailed statement of Theorem 2, and Theorem 6, as an variation of Theorem 2′ under stronger conditions, are then provided, followed by their proofs. The proofs of three lemmas are given at the end of the section.

The main process we analyze is (2.3). Lemma 3 gives an estimate of the final distance $E|g(X_{T}^{\pi})-Y_{T}^{\pi}|^{2}$ provided by (2.3) in terms of the deviation between the approximated variables $Y_{0}^{\pi},Z_{t_{i}}^{\pi}$ and the true solutions.

for $0\leq i\leq N-1$, $j=1,2$. Let $\delta X_{i}=X_{t_{i}}^{\pi,1}-X_{t_{i}}^{\pi,2}$, $\delta Z_{i}=Z_{t_{i}}^{\pi,1}-Z_{t_{i}}^{\pi,2}$, then for any $\lambda_{7}>f_{z}$, and sufficiently small $h$, we have

where $A_{5}\coloneqq-h^{-1}\ln[1-(2k_{f}+\lambda_{7})h]$.

Lemma 5 shows that, similar to the nonlinear Feynman–Kac formula, the discrete stochastic process defined in (2.3) can also be linked to some deterministic functions.

Now we are ready to prove Theorem 2, with a precise statement given below. Like Theorem 1′, the conditions below are satisfied if any of the five cases of the weak coupling and monotonicity conditions holds to certain extent.

Suppose Assumptions 1, 2, 3, and 4 hold true. Given any $\lambda_{1},\lambda_{3}>0$, $\lambda_{2}\geq f_{z}$, and $\lambda_{7}>f_{z}$, let $\overline{A_{i}},(i=1,2,3,4)$ be defined in (4.4) and

If there exist $\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{7}$ satisfying $\overline{A_{0}}^{\prime}<1$ and $\overline{B_{0}}<1$, then there exists a constant C depending on $E|\xi|^{2}$, $\mathscr{L}$, $T$, $\lambda_{1}$, $\lambda_{2}$, $\lambda_{3}$, and $\lambda_{7}$, such that for sufficiently small $h$,

If we take the infimum within the domains of $Y_{0}^{\pi}$ and $Z_{t_{i}}^{\pi}$ on both sides, we recover the original statement in Theorem 2.

If any of the weak coupling and monotonicity conditions introduced in Assumption 3 holds to a sufficient extent, there must exist $\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{7}$ satisfying the conditions in Theorem 2′. The arguments are very similar to those provided in Remark Remark. Hence, we omit the details here for the sake of brevity.

Using Lemma 3 with $\lambda_{4}>0$, we obtain

Plugging estimates (5.6)(5.7) into (5.5) gives us

It remains to estimate the term $\sum_{i=0}^{N-1}E|\overline{Z}_{t_{i}}^{\pi}-E[\overline{Z}_{t_{i}}^{\pi}|X_{t_{i}}^{\pi},Y_{t_{i}}^{\pi}]|^{2}h$, to which we intend to apply Lemma 4. Let $X_{t_{i}}^{\pi,1}=X_{t_{i}}^{\pi}$ and $X_{t_{i}}^{\pi,2}=\overline{X}_{t_{i}}^{\pi}$. The associated $Z_{t_{i}}^{\pi,1}$ and $Z_{t_{i}}^{\pi,2}$ are then defined according to equation (5.1). Note that $Z_{t_{i}}^{\pi,2}=\overline{Z}_{t_{i}}^{\pi}$ but $Z_{t_{i}}^{\pi,1}$ is not necessarily equal to $Z_{t_{i}}^{\pi}$, due to the possible violation of the terminal condition. From Lemma 5, we know $Z_{t_{i}}^{\pi,1}$ can be represented as $V_{i}^{\pi}(X_{t_{i}}^{\pi},Y_{t_{i}}^{\pi})$ with $V_{i}^{\pi}$ being a deterministic function. By the property of conditional expectation, we have

for any $V_{i}$. Therefore we have the estimate

Recall that $\delta X_{i}=X_{t_{i}}^{\pi}-\overline{X}_{t_{i}}^{\pi},\delta Z_{i}=Z_{t_{i}}^{\pi,1}-\overline{Z}_{t_{i}}^{\pi}$. Similar to the derivation of estimate (4.16) (using a given $\lambda_{3}>0$ without final specification) in the proof of Theorem 1′, when $\overline{A_{0}}^{\prime}<1$, we have

in which $\overline{P}=\max_{0\leq n\leq N}e^{-\overline{A_{1}}nh}E|\delta X_{i}|^{2}$. Plugging (5.10) into (5.9), and then into (5.8), we get

for sufficiently small $h$. Here $B(h)$ is defined as

The forms of inequalities (5.4) and (5.11) are already very close. When $\displaystyle{\lim_{h\rightarrow 0}B(h)}=\overline{B_{0}}<1$, there exists $\lambda_{4}>0$ such that for sufficiently small $h$, we have $1-(1+\lambda_{4})^{3}B(h)>\frac{1}{2}(1-\overline{B_{0}})$. Rearranging the term $E|g(X_{T}^{\pi})-Y_{T}^{\pi}|^{2}$ in inequality (5.11) yields our final estimate. ∎

Suppose Assumptions 1, 2, 3, 4 and the assumptions in Theorem 3 hold true. Let $u$ be the solution of corresponding quasilinear PDEs (2.7) and $L$ be the squared Lipschitz constant of $\sigma^{\operatorname{T}}(t,x,u(t,x))\nabla_{x}u(t,x)$ with respect to x. With the same notations of Theorem 2′, when $\overline{A_{0}}^{\prime}<1$ and

there exists a constant $C>0$ depending on $E|\xi|^{2}$, $T$, $\mathscr{L}$, $L$, $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{3}$, such that for sufficiently small $h$,

where $f_{i}(x)=\sigma^{\operatorname{T}}(t_{i},x,u(t_{i},x))\nabla_{x}u(t_{i},x)$.

By Theorem 3, we have $Z_{t_{i}}=f_{i}(X_{t_{i}})$, in which $X_{t}$ is the solution of

Using Lemma 3 again with $\lambda_{4}>0$ gives us

Similar to the arguments in inequalities (5.6)(5.7), we have

where the last equality uses the convergence result (3.4). Plugging it into (5.13), we have

We employ the estimate (5.10) again to rewrite inequality (5.14) as

The Lipschitz constant used in Theorem 6 may be further estimated a priori. Denote the Lipschitz constant of function $f$ with respect to $x$ as $L_{x}(f)$, and the bound of function $f$ as $M(f)$. Loosely speaking, we have

Here $L_{x}(u)=M(\nabla_{x}u(t,x))$ can be estimated from the first point of Theorem 4 and $L(\nabla_{x}u(t,x))=M(\nabla_{xx}u)$ can be estimated through the Schauder estimate (see, e.g., [32, Chapter 4, Lemma 2.1]). Note that the resulting estimate may depend on the dimension $d$.

We construct continuous processes $X_{t}^{\pi},Y_{t}^{\pi}$ as follows. For $t\in[t_{i},t_{i+1})$, let

From system (2.3), we see this definition also works at $t_{i+1}$. We are interested in again the estimates of the following terms

For any $\lambda_{5},\lambda_{6}>0$, using Assumptions 1, 2 and the RMS-GM inequality, we have

in which we choose $\epsilon_{1}=\lambda_{6}(K\lambda_{5}^{-1}+\sigma_{x})^{-1}$ and $\epsilon_{2}=\lambda_{6}(b_{y}\lambda_{5}^{-1}+\sigma_{y})^{-1}$. The path regularity in Theorem 4 tells us

Plugging inequalities (5.17)(5.18)(5.19) into (5.16) with simplification, we obtain

where $A_{6}\coloneqq K\lambda_{5}^{-1}+\sigma_{x}+\lambda_{5}+\lambda_{6}$, $A_{7}\coloneqq b_{y}\lambda_{5}^{-1}+\sigma_{y}+2\lambda_{6}$, and $h$ is sufficiently small.

Similarly, with the same type of estimates in (5.16)(5.20), for any $\lambda_{5},\lambda_{6}>0$, we have

Arguing in the same way of (5.21), by Grönwall inequality, for sufficiently small $h$, we have

with $A_{8}\coloneqq K\lambda_{5}^{-1}+\lambda_{5}+\lambda_{6}$, $A_{9}\coloneqq f_{x}\lambda_{5}^{-1}+2\lambda_{6}$. Choosing $\epsilon_{3}=(1+f_{z}\lambda_{5}^{-1}+\lambda_{6})^{-1}\lambda_{6}$ and using

with $A_{10}\coloneqq 1+f_{z}\lambda_{5}^{-1}+2\lambda_{6}$.

Combining inequalities (5.21)(5.22) together yields

Letting $A_{11}\coloneqq\max\{A_{6},A_{8}\}+\max\{A_{7},A_{9}\}$, we have

We start from $M_{0}=E|Y_{0}-Y_{0}^{\pi}|^{2}$ and apply inequality (5.23) repeatedly to obtain

in which for the last term we use the fact $\sum_{i=0}^{N-1}E_{z}^{i}\leq Ch$ from inequality (3.2). Note that

Given any $\lambda_{4}>0$, we can choose $\lambda_{6}$ small enough such that

This condition and inequality (5.24) together give us

Finally, by decomposing the objective function, we have

We use the same notations as in the proof of Lemma 1. As derived in (4.12), for any $\lambda_{7}>f_{z}\geq 0$, we have

Multiplying both sides of (5.27) by $e^{A_{5}ih}(e^{-A_{5}T}\vee 1)/(1-f_{z}\lambda_{7}^{-1})$ gives us

Summing (5.28) up from $i=0$ to $N-1$, we obtain

Note that $E|\delta Y_{N}|^{2}\leq g_{x}E|\delta X_{N}|^{2}$ by Assumption 1. Plugging it into (5.29), we arrive at the desired result. ∎

Consider the following map defined on $H_{k}$:

By Assumption 3, $\Phi_{k}(Y)$ is square-integrable. Furthermore, following the same argument for $Z_{t_{k}}^{\pi,^{\prime}}$, $\Phi_{k}(Y)$ can also be represented as a deterministic function of $X_{t_{k}}^{\pi},Y_{t_{k}}^{\pi}$. Hence, $\Phi_{k}(Y)\in H_{k}$. Note that Assumption 1 implies $E|\Phi_{k}(Y_{1})-\Phi_{k}(Y_{2})|^{2}\leq Kh^{2}E|Y_{1}-Y_{2}|^{2}$. Therefore $\Phi_{k}$ is a contraction map on $H_{k}$ when $h<1/\sqrt{K}$. By the Banach fixed-point theorem, there exists a unique fixed-point $Y^{*}=\phi_{k}^{*}(X_{t_{k}}^{\pi},Y_{t_{k}}^{\pi})\in H_{k}$ satisfying $Y^{*}=\Phi_{k}(Y^{*})$. We choose $U_{k}^{\pi}=\phi_{k}^{*}$ to validate the statement for $Y_{t_{k}}^{\pi,^{\prime}}$.

When $b$ and $\sigma$ are independent of $y$, all of the arguments above can be made similarly with $U_{i}^{\pi},V_{i}^{\pi}$ also being independent of $Y$. ∎

Numerical Examples

2 Example 1

The first problem is adapted from , in which the original spatial dimension of the problem is 1. We consider the following coupled FBSDEs

where $X_{j,t},Z_{j,t},W_{j,t}$ denote the $j$-th components of $X_{t},Y_{t},W_{t}$, and the coefficient functions are given as

It can be verified by Itô’s formula that the $Y$ part of the solution of (6.1) is given by

Let $\xi=(1,1,\dots,1)$ ($100$-dimensional), $T=5,N=160$. The initial guess of $Y_{0}$ is generated from a uniform distribution on the interval $$while the true value of$Y_{0}\approx 0.81873$. We train 25000 steps with an exponential decay learning rate that decays every 100 steps, with the starting learning rate being 1e-2 and ending learning rate being 1e-5. Figure 1 illustrates the mean of the loss function and relative approximation error of$Y_{0}$against the number of iteration steps. All runs converged and the average final relative error of$Y_{0}$is$0.39\%$.

3 Example 2

The second problem is adapted from , in which the spatial dimension is originally tested up to 10. The coupled FBSDEs are given by

where $\sigma>0,r,D$ are constants. One can easily check by Itô’s formula that the $Y$ part of the solution of (6.2) is

References