# Convergence of the deep BSDE method for stochastic control problems
formulated through the stochastic maximum principle

Zhipeng Huang  Balint Negyesi  Cornelis W. Oosterlee

###### Abstract

It is well-known that decision-making problems from stochastic control can be
formulated by means of forward-backward stochastic differential equation
(FBSDE). Recently, the authors of [1] proposed an efficient deep learning-
based algorithm which was based on the stochastic maximum principle (SMP). In
this paper, we provide a convergence result for this deep SMP-BSDE algorithm
and compare its performance with other existing methods. In particular, by
adopting a similar strategy as in [2], we derive a posteriori error estimate,
and show that the total approximation error can be bounded by the value of the
loss functional and the discretization error. We present numerical examples
for high-dimensional stochastic control problems, both in case of drift- and
diffusion control, which showcase superior performance compared to existing
algorithms.

###### keywords:

stochastic control, deep SMP-BSDE, stochastic maximum principle, vector-valued
FBSDE

††journal: Mathematics and Computers in Simulation

[label1] organization=Mathematical Institute, Utrecht University

[label2] organization=Delft Institute of Applied Mathematics (DIAM), Delft
University of Technology

##  1 Introduction

Stochastic control theory is a powerful paradigm for modelling and analyzing
decision-making problems that are subject to some random dynamics. Classical
approaches for solving these kinds of problems include methods based on the
dynamic programming principle (DP) [3], the stochastic maximum principle (SMP)
[4] [5] and other techniques, see e.g. [6] [7] [8] [9]. However, these
approaches cannot easily handle high-dimensional problems and suffer from the
“curse of dimensionality”. A recent candidate solution technique for
stochastic control problems is formed by deep learning-based approaches, due
to their remarkable performance in high-dimensional settings. The paper [10]
developed a deep learning algorithm that directly approximates the optimal
control process by a neural network at each step in time, and by training a
terminal loss functional for all time steps simultaneously. Similar approaches
have been explored in the control theory community [11], [12] and [13], before
the rise of computing power and machine learning.

Inspired by the remarkable performance in [10], the research community has
developed several neural network-based algorithms for stochastic control
problems, see e.g. [14] [15] [16]. Many of these algorithms build upon
deriving a forward-backward stochastic differential equation (FBSDE),
associated to the control problem. Pioneered by the well-known deep BSDE
method, initially proposed by [17] and later extended by [2] to the coupled
setting, several solution approaches have been proposed, showing outstanding
empirical performance in high-dimensional frameworks. These methods derive the
FBSDE through a stochastic representation of the solution of the Hamilton-
Jacobi-Bellman (HJB) equation, motivated by the non-linear Feyman-Kac formula.
However, such techniques become infeasible when the diffusion of the state
process is also controlled, as they do not solve for the value function’s
second derivative, which is necessary to compute the optimal (diffusion)
control. As a remedy, enabling diffusion control, the authors in [1] proposed
a deep BSDE algorithm where the associated FBSDE is derived from the
stochastic maximum principle (SMP), which we call the deep SMP-BSDE. Other
SMP-based algorithms include [18] [19] in the context of mean-field control
and mean-field games. We refer to [20] and [21] for a detailed overview of
deep learning algorithms for stochastic control problems.

Furthermore, the authors in [22] found that deep BSDE methods may fail to
converge for FBSDEs stemming from stochastic control problems via DP, due to
local minima. They proposed a robust counterpart by adding a regularization
component to the loss function which resolved this issue in the case of drift
control. Despite the fact that there are many studies concerning stochastic
control, only a few theoretical derivations are available regarding the
convergence of machine learning-based approaches for FBSDEs stemming from
stochastic control. Theoretical work in this direction includes the study of
convergence of the deep BSDE method by [2], which provides a posteriori error
estimate. The authors of [22] prove the convergence of a robust deep BSDE
method by exploiting the special structure of the FBSDE. For works in the
decoupled framework see e.g. [23], [24] and [25].

In this paper, we provide convergence results for the deep SMP-BSDE algorithm
and compare the method with existing algorithms supported by numerical
results. Unlike in [17] [2] and [22], we consider FBSDEs that come from the
SMP, instead of the HJB equation, and therefore the results and standard
estimates for FBSDE stemming from DP can not be directly applied.
Nevertheless, with some extra effort, we are able to adopt a similar strategy
as in [2], and derive a posteriori error estimate for the numerical solutions
of the deep SMP-BSDE algorithm. By this, we are able to tackle diffusion
control problems.

This paper is organized as follows: In Section 2, we briefly review the
theoretical foundations related to the SMP. In Section 3, we formulate the
deep SMP-BSDE algorithm for diffusion control problems. We carry out a
convergence analysis in Section 4, and, in particular, a posteriori error
estimate is derived for the deep SMP-BSDE algorithm. In Section 5, we
demonstrate the performance of the algorithm through numerical examples both
in the case of drift- and diffusion control.

##  2 Background

In this section, we review some basic results from stochastic control theory
and show how to reformulate a stochastic control problem into an FBSDE through
the SMP.

Let $\left(\Omega,\mathcal{F},\left\\{\mathcal{F}_{t}\right\\}_{0\leq t\leq
T},\mathbb{P}\right)$ be a filtered probability space, supporting an
$m$-dimensional Brownian motion $W$ and its natural filtration
$\mathcal{F}=\left\\{\mathcal{F}_{t}\right\\}_{0\leq t\leq T}$, augmented by
all $\mathbb{P}$-null sets. Fixing $0<T<\infty$, we consider the following
finite horizon stochastic control problem

$\left\\{\begin{aligned}
&\inf_{u\in\mathcal{U}[0,T]}J\left(0,x_{0};u(\cdot)\right):=\inf_{u\in\mathcal{U}[0,T]}\mathbb{E}\left(\int_{0}^{T}\bar{f}(s,X_{s},u_{s})\,ds+g(X_{T})\right),\\\
&X_{t}=x_{0}+\int_{0}^{t}\bar{b}(s,X_{s},u_{s})\,ds+\int_{0}^{t}\bar{\sigma}(s,X_{s},u_{s})\,dW_{s},\quad
t\in[0,T],\end{aligned}\right.$ (1)

where
$\bar{b}:[0,T]\times\mathbb{R}^{d}\times\mathbb{R}^{\ell}\to\mathbb{R}^{d}$,
$\bar{\sigma}:[0,T]\times\mathbb{R}^{d}\times\mathbb{R}^{\ell}\to\mathbb{R}^{d\times
m}$, $\bar{f}:[0,T]\times\mathbb{R}^{d}\times\mathbb{R}^{\ell}\to\mathbb{R}$
and $g:\mathbb{R}^{d}\to\mathbb{R}$ are deterministic functions, and
$X_{t},u_{t}$ are $\mathbb{R}^{d},\mathbb{R}^{\ell}$-valued stochastic
processes, respectively. The set of admissible controls, $\mathcal{U}[0,T]$,
is defined as

$\mathcal{U}[0,T]\coloneqq\left\\{u:[0,T]\times\Omega\rightarrow U\mid u\in
L_{\mathcal{F}}^{2}\left(0,T;\mathbb{R}^{\ell}\right)\right\\},$

with

$L_{\mathcal{F}}^{2}\left(0,T;\mathbb{R}^{\ell}\right)\coloneqq\left\\{x:[0,T]\times\Omega\rightarrow\mathbb{R}^{\ell}\mid
x\text{ is }\mathcal{F}\text{-adapted and
}E\left[\int_{0}^{T}\left\|x_{t}\right\|^{2}\,dt\right]<\infty\right\\},$

where we shall denote $\|\cdot\|$ for both the usual Euclidean norm and the
Frobenius norm for matrices. We assume that the control domain $U$ is a convex
body in $\mathbb{R}^{\ell}$.

Any process $u_{t}\in\mathcal{U}[0,T]$ is called an admissible control of (1),
and the $(X_{t},u_{t})$ consisting of the corresponding state process is
called an admissible pair. Furthermore, $(X_{t},u_{t})$ is an optimal pair
whenever the infimum of (1) is achieved, and accordingly, we define the value
function $V(\cdot,\cdot)$ of (1) as follows

$\left\\{\begin{aligned}
&V(t,x)\coloneqq\inf_{u\in\mathcal{U}[t,T]}J(t,x;u(\cdot)),\quad\forall(t,x)\in[0,T)\times\mathbb{R}^{d},\\\
&V(T,y)\coloneqq g(y),\quad\forall y\in\mathbb{R}^{d},\end{aligned}\right.$
(2)

and denote its partial derivatives as $V_{x}\coloneqq\partial_{x}V(t,x)$,
$V_{xx}\coloneqq\partial_{xx}V(t,x)$ and $\partial_{t}V(t,x)$ without
specifying their dependencies in $(t,x)$.

######  Remark 1.

We shall further distinguish two important classes of problem (1). We call (1)
a drift control problem if the drift coefficient depends on the control
$u_{t}$ but not the diffusion coefficient, and it is called a diffusion
control problem if the diffusion coefficient also includes $u_{t}$ as an
argument.

Associated with the stochastic control problem (1), we introduce the adjoint
equation, as follows

$\left\\{\begin{aligned}
&P_{t}=P_{0}-\int_{0}^{t}\nabla_{x}\bar{H}(s,X_{s},u_{s},P_{s},Q_{s})\,ds+\int_{0}^{t}Q_{s}\,dW_{s},\quad
t\in[0,T]\\\ &P_{T}=-\nabla_{x}g(X_{T}),\end{aligned}\right.$ (3)

where the $\nabla_{x}\bar{H}$ is the derivative of Hamiltonian $\bar{H}$,
which is defined by

$\displaystyle\bar{H}(t,x,u,p,q):=p^{\top}\bar{b}(t,x,u)+\operatorname{Tr}\left(q^{\top}\bar{\sigma}(t,x,u)\right)-\bar{f}(t,x,u),$
(4) $\displaystyle\quad\forall(t,x,u,p,q)\in[0,T]\times\mathbb{R}^{d}\times
U\times\mathbb{R}^{d}\times\mathbb{R}^{d\times m}.$

Equation (3) is a BSDE whose solution is formed by a pair of processes,
$(P(\cdot),Q(\cdot))\in
L_{\mathcal{F}}^{2}\left(0,T;\mathbb{R}^{d}\right)\times\left(L_{\mathcal{F}}^{2}\left(0,T;\mathbb{R}^{d}\right)\right)^{m}$.
Concerning the wellposedness of BSDE (3), we denote the maps
$\bar{\varphi}\coloneqq\\{\bar{b},\bar{\sigma},\bar{f},g\\}$, and state the
following assumption first.

######  Assumption 1.

The map $\bar{\varphi}$ is $C^{2}$ in $x$, and $\bar{\varphi}(t,0,u)$ is
bounded for any $u\in U$. Moreover, $\bar{\varphi}$, $\bar{\varphi}_{x}$ and
$\bar{\varphi}_{xx}$ are uniformly Lipschitz in $x$ and $u$.

######  Remark 2.

With Assumption 1, the adjoint equation, or BSDE (3), admits a unique solution
for every admissible pair $(X_{t},u_{t})$. However, there is only one of the
admissible 4-tuples $(X_{t},u_{t},P_{t},Q_{t})$ that also minimizes the
objective function in (1). That is to say, to have an FBSDE that admits a
unique solution without including the objective function, we need extra effort
in the reformulation. For this purpose, we recall the SMP.

######  Theorem 1 (Stochastic maximum principle).

Let Assumption 1 hold, and
$\left(X_{t}^{*},u_{t}^{*},P_{t}^{*},Q_{t}^{*}\right)$ be an admissible
4-tuple. Suppose that $g(\cdot)$ is convex,
$\bar{H}\left(t,\cdot,\cdot,P_{t}^{*},Q_{t}^{*}\right)$ defined by (4) is
concave for all $t\in[0,T]$ $\mathbb{P}$ almost surely, and the maximum
condition

$\bar{H}\left(X_{t}^{*},u_{t}^{*},P_{t}^{*},Q_{t}^{*}\right)=\max_{u\in
U}\bar{H}\left(X_{t}^{*},u,P_{t}^{*},Q_{t}^{*}\right),\quad\text{a.e.}\
t\in[0,T],\quad\mathbb{P}\text{-a.s.}$ (5)

holds. Then, $\left(X_{t}^{*},u_{t}^{*}\right)$ is an optimal pair of problem
(1).

######  Remark 3.

The proof of this theorem is found in [26, pp. 149-150], or, in a more general
setting, [27, pp. 138-140]. Theorem 1 provides a sufficient condition for the
optimal control $u^{*}_{t}$, when certain concavity and convexity conditions
hold, which are crucial in general, see for instance, [27, Example 3.1, pp.
138-140]. On the other hand, Theorem 1 itself does not constitute a necessary
condition unless there is no diffusion control in problem (1). In the rest of
this paper, we shall use superscript $*$ to indicate that a process is
associated with the optimal control $u_{t}^{*}$.

Under sufficient smoothness and concavity assumptions of $\bar{H}$ in Theorem
1, the optimization (5) is uniquely solved by the following first-order
conditions,

$\displaystyle\nabla_{u}\bar{H}(t,X^{*}_{t},u^{*}_{t},P^{*}_{t},Q^{*}_{t})=$
$\displaystyle(\nabla_{u}\bar{b}(t,X^{*}_{t},u^{*}))^{\top}P^{*}_{t}$ (6)
$\displaystyle+\nabla_{u}\operatorname{Tr}(\bar{\sigma}^{\top}(t,X^{*}_{t},u^{*}_{t})Q^{*}_{t})-\nabla_{u}\bar{f}(t,X^{*}_{t},u^{*}_{t})=0.$
(7)

Similar to algorithm 3 in [1], we assume that the first-order conditions
yields an explicit formula for the mapping
$\mathcal{M}:(t,X_{t}^{*},P_{t}^{*},Q_{t}^{*})\mapsto u_{t}^{*}$,

$u_{t}^{*}=\mathcal{M}(t,X^{*}_{t},P^{*}_{t},Q^{*}_{t}),\quad\forall
t\in[0,T].$ (8)

We will call this function the _feedback map_ , and remark that for a rather
wide range of interesting problems such an expression is available in closed-
form. Let us denote the composition
$\varphi\coloneqq\bar{\varphi}(t,x,\mathcal{M}(t,x,p,q))$ for
$\bar{\varphi}\setminus\\{g\\}$, where we define
$\varphi=\\{b,\sigma,f\\}\cup\\{g\\}$. Similarly, we define the function
$\bar{F}(t,x,u,p,q)\coloneqq\nabla_{x}\bar{H}(t,x,u,p,q)$ and write

$F(t,x,p,q)\coloneqq\bar{F}(t,x,\mathcal{M}(t,x,u,p,q),p,q).$ (9)

With this in hand, we reformulate (3) and the controlled SDE of (1) as a
fully-coupled FBSDE, for $t\in[0,T]$,

$\left\\{\begin{aligned}
X_{t}&=x_{0}+\int_{0}^{t}b(s,X_{s},P_{s},Q_{s})\,ds+\int_{0}^{t}\sigma(s,X_{s},P_{s},Q_{s})\,dW_{s},\\\
P_{t}&=P_{0}-\int_{0}^{t}\nabla_{x}H(s,X_{s},P_{s},Q_{s})\,ds+\int_{0}^{t}Q_{s}\,dW_{s},\\\
P_{T}&=-\nabla_{x}g(X_{T}).\end{aligned}\right.$ (10)

We call the BSDE part in (10), subject to its boundary condition, the SMP-
BSDE.

######  Remark 4.

An important feature of (10) is that its solution is equivalent to the
solution of (1). Let $(X^{*}_{t},u^{*}_{t})$ be an optimal pair obtained by
(1), then there exists a unique pair $(P^{*}_{t},Q^{*}_{t})$ which follows
from Remark 2. Therefore, we have the optimal
$(X^{*}_{t},P^{*}_{t},Q^{*}_{t})$ that also solves (10). On the other hand,
solving (10) gives us $(\tilde{X}_{t},\tilde{P}_{t},\tilde{Q}_{t})$, and we
obtain $\tilde{u}_{t}$ by the feedback map (8). Consequently, we must have
$\tilde{u}_{t}=u^{*}_{t}$ and
$(\tilde{X}_{t},\tilde{P}_{t},\tilde{Q}_{t})=(X^{*}_{t},P^{*}_{t},Q^{*}_{t})$,
due to the uniqueness of feedback map (8) and Theorem 1.

##  3 Existing algorithms and the deep SMP-BSDE algorithm

In this section, we first give a comparison between the existing deep BSDE
algorithm and the deep SMP-BSDE algorithm, and show why the former can not be
used to solve diffusion control problems.

Suppose the DP holds, then the value function $V(t,x)$ of problem (1) solves
the following HJB equation

$\left\\{\begin{aligned} &-\partial_{t}V+\sup_{u\in
U}\mathcal{G}\left(t,x,u,-V_{x},-V_{xx}\right)=0,\\\
&V(T,x)=g(x),\end{aligned}\right.$ (11)

where the generalized Hamiltonian $\mathcal{G}$ is defined by

$\mathcal{G}\left(t,x,u,-V_{x},-V_{xx}\right)\coloneqq-\frac{1}{2}\operatorname{Tr}\left(\bar{\sigma}(t,x,u)^{\top}V_{xx}\bar{\sigma}(t,x,u)\right)-V_{x}^{\top}\bar{b}(t,x,u)-\bar{f}(t,x,u)$
(12)

By the stochastic verification theorem, see [27, pp. 268-269], a given
admissible pair $(X_{t}^{*},u_{t}^{*})$ is optimal if and only if

$\displaystyle\mathcal{G}\left(t,X_{t}^{*},u_{t}^{*},-V_{x},-V_{xx}\right)=\max_{u\in
U}\mathcal{G}\left(t,X_{t}^{*},u,-V_{x},-V_{xx}\right),$ (13)

provided that $V$ is a classical solution to (11) with sufficient smoothness.

Suppose (11) admits a unique classical solution given by the value function
$V(t,x)$. Then, it is straightforward to obtain a new feedback map
$\tilde{\mathcal{M}}$ by maximizing $\mathcal{G}$ over $u$, which allows us to
write the optimal control $u^{*}_{t}$ in terms of $X^{*}_{t}$, $V_{x}$ and
$V_{xx}$

$u^{*}_{t}=\tilde{\mathcal{M}}(t,X^{*}_{t},-V_{x},-V_{xx}),\qquad t\in[0,T].$

Notice that whenever the diffusion includes the control variable,
$\tilde{\mathcal{M}}$ will always depend on $V_{xx}$. From the stochastic
representation of $V(t,x)$, we obtain for $t\in[0,T]$

$V(t,x)=g(X^{*}_{T})+\int_{t}^{T}\bar{f}(s,X^{*}_{s},u^{*}_{s})\,ds-\int_{t}^{T}V_{x}^{\top}\bar{\sigma}(t,X_{t}^{*},u_{t}^{*})\,dW_{s}.$
(14)

Equation (14) defines a BSDE whose unique solution coincides with
$(Y_{t},Z_{t})=\left(V(t,X_{t}^{*}),V_{x}^{\top}\bar{\sigma}(t,X^{*}_{t},u^{*}_{t})\right)$,
whenever the feedback map $\tilde{\mathcal{M}}$ is only a function of $V_{x}$
but not of $V_{xx}$. However, this is only the case when the diffusion
coefficient does not depend on $u$, i.e., in the case of drift control. In
particular, equation (14) coincides with the deep BSDE method of [17] [2],
[22], applied to stochastic control problems. In numerical settings where one
does not have direct access to $V_{xx}$, this makes the dynamic programming
approaches described above infeasible, whenever the diffusion coefficient
depends on $u$.

Moreover, the robust counterpart in [22] heavily builds upon the fact that the
corresponding DP approach results in a BSDE whose driver is independent of
$Y_{t}=V(t,X_{t}^{*})$, and uses the fact that $Y_{0}=V(t_{0},x_{0})$ is a
deterministic quantity coinciding with the value function of the control
problem. Hence, both its mean and variance should be minimized, and based on
these terms a new loss functional can be formulated. However, this robustness
technique does not apply to the SMP-BSDE setting, as in general $F$ in (9)
depends on $P_{t}$, and since $P_{0}=-V_{x}(t_{0},x_{0})$ its mean is not
necessarily minimized at $t_{0}$. From this point of view, the deep SMP-BSDE
algorithm seems more promising for solving general stochastic control
problems.

For the reasons above, in order to be able to treat diffusion control
problems, in what follows we focus on the SMP-based FBSDE formulation. Let us
state the discrete scheme of the deep SMP-BSDE algorithm originating from
Algorithm 3 of [1], but instead of using one single neural network for both
$P_{0}$ and the process $Q_{t}$, as in the aforementioned paper, we
approximate $P_{0}$ and $Q_{t}$ at each step in time by a separate neural
network, respectively. We consider the following classical Euler scheme
corresponding to (10)

$\displaystyle\inf_{\mu_{0}^{\pi}\in\theta_{0}^{P},\phi_{i}^{\pi}\in\theta_{i}^{Q}}$
$\displaystyle
E\left\|-\nabla_{x}g\left(X_{T}^{\pi}\right)-P_{T}^{\pi}\right\|^{2},$ (15a)
$\displaystyle\text{s.t.}\left\\{\begin{aligned} &X_{0}^{\pi}=x_{0},\\\
&P_{0}^{\pi}=\mu_{0}^{\pi}(x_{0}),\\\ &X_{t_{i+1}}^{\pi}=\begin{aligned}
X_{t_{i}}^{\pi}&+b\left(t_{i},X_{t_{i}}^{\pi},P_{t_{i}}^{\pi},Q_{t_{i}}^{\pi}\right)\Delta
t\\\
&+\sigma\left(t_{i},X_{t_{i}}^{\pi},P_{t_{i}}^{\pi},Q_{t_{i}}^{\pi}\right)\Delta
W_{i},\end{aligned}\\\
&Q_{t_{i}}^{\pi}=\phi_{i}^{\pi}\left(X_{t_{i}}^{\pi}\right),\\\
&P_{t_{i+1}}^{\pi}=P_{t_{i}}^{\pi}-F\left(t_{i},X_{t_{i}}^{\pi},P_{t_{i}}^{\pi},Q_{t_{i}}^{\pi}\right)\Delta
t+\left(Q_{t_{i}}^{\pi}\right)^{\top}\Delta W_{i},\end{aligned}\right.$ (15b)

with a time partition, $\pi:0=t_{0}<t_{1}<\cdots<t_{N}=T$, $h=T/N$, $t_{i}=ih$
and $\Delta W_{i}:=W_{t_{i+1}}-W_{t_{i}}$ for $i=0,1,\ldots,N-1$. We recall
the definitions of $b$, $\sigma$ and $F$ in Section 2. Moreover, we let
$\theta_{0}^{P}$ and $\theta_{i}^{Q}$ be the corresponding parametric spaces
for the neural networks
$\theta_{0}^{P}\ni\mu_{0}^{\pi}:\mathbb{R}^{d}\to\mathbb{R}^{d}$ and
$\theta_{i}^{Q}\ni\phi_{i}^{\pi}:\mathbb{R}^{d}\to\mathbb{R}^{d\times m}$,
respectively. The objective function (15a) serves as the loss function in the
machine learning algorithm, and through the training of the neural networks we
wish to find appropriate functions $\mu_{0}^{\pi}(x_{0})$ and
$\phi_{i}^{\pi}\left(X_{t_{i}}^{\pi}\right)$ that can approximate $P_{0}$ and
$Q_{t_{i}}$ sufficiently well. The complete pseudo-code for the deep SMP-BSDE
algorithm is given in Algorithm 1, which will be used in the later Section 5
for numerical experiments.

Algorithm 1 deep SMP-BSDE algorithm

1:Input: Initial parameters
$\left(\theta^{P}_{0},\theta^{Q}_{0},\ldots,\theta^{Q}_{N-1}\right)$, learning
rate $\eta$; batch size $M$; number of iteration $K$.

2:Data: Simulated Brownian increments $\left\\{\Delta
W_{t_{i},k}\right\\}_{0\leq i\leq N-1,1\leq k\leq K}$

3:Output: The triple $(X_{t_{i}},P_{t_{i}},Q_{t_{i}})$

4:for $k=1$ to $K$ do

5: $X_{t_{0},k}^{\pi}=x_{0}$,
$P_{t_{0},k}^{\pi}=\mu_{0}^{\pi}(x_{0};\theta^{P}_{0})$

6: for $i=0$ to $N-1$ do

7:
$Q_{t_{i},k}^{\pi}=\phi_{i}^{\pi}\left(X_{t_{i}}^{\pi};\theta^{Q}_{i}\right)$

8:
$u_{t_{i},k}^{\pi}=\mathcal{M}\left(t_{i},X_{t_{i},k}^{\pi},P_{t_{i},k}^{\pi},Q_{t_{i},k}^{\pi}\right)$

9:
$X_{t_{i+1},k}^{\pi}=X_{t_{i},k}^{\pi}+\bar{b}\left(t_{i},X_{t_{i},k}^{\pi},u_{t_{i},k}^{\pi}\right)\Delta
t_{i}+\bar{\sigma}\left(t_{i},X_{t_{i},k}^{\pi},u_{t_{i},k}^{\pi}\right)\Delta
W_{t_{i},k}$

10:
$P_{t_{i+1},k}^{\pi}=P_{t_{i},k}^{\pi}-\bar{F}\left(t_{i},X_{t_{i},k}^{\pi},u_{t_{i},k}^{\pi}\right)\Delta
t_{i}+Q_{t_{i},k}^{\pi}\Delta W_{t_{i},k}$

11: end for

12:
$\text{Loss}=\frac{1}{M}\sum_{j=1}^{M}\left\|-\nabla_{x}g\left(X_{t_{N},k}^{\pi}\right)-P_{t_{N},k}^{\pi}\right\|^{2}$

13:
$\left(\theta^{P}_{0},\theta^{Q}_{0},\ldots,\theta^{Q}_{N-1}\right)\longleftarrow\left(\theta^{P}_{0},\theta^{Q}_{0},\ldots,\theta^{Q}_{N-1}\right)-\eta\nabla\text{Loss}$

14:end for

##  4 Convergence Analysis

This section is dedicated to the convergence analysis for the deep SMP-BSDE
method, reviewed in Section 2, and the discrete scheme (15b). In particular,
we will show that the errors of the numerical solution to the FBSDE are
bounded by the simulation error of the objective function, and such errors in
theory could be made arbitrarily small, due to the universal approximation
theorem. Our analysis follows a similar strategy as [2]. For the sake of this
section, in order to be able to use techniques established therein, we need
the following restriction.

######  Assumption 2.

The coefficients $\bar{b},\bar{\sigma}$ in (1) imply a feedback map such that
the compositions $\bar{b}(t,x,\mathcal{M}(t,x,p,q))$ and
$\bar{\sigma}(t,x,\mathcal{M}(t,x,p,q))$ do not depend on $q$.

Even though this condition may seem abstract, note that it in particular
covers the important subclass of drift control problems, i.e. whenever
$\bar{\sigma}$ does not depend on $u$. In that case, under the convexity
conditions of Theorem 1, the feedback map in (8) is only a function of
$(t,X_{t},P_{t})$ and not of $Q_{t}$ – see the first-order conditions in (6).
Consequently, the solution pair of the backward equation in the FBSDE (10)
only couples into the forward component via $P_{t}$, similarly as in [2].
Nevertheless, in what follows we follow a more general presentation with
Assumption 2, such that special cases of diffusion control problems can also
be treated.

Our main contribution is establishing a convergence result for the numerical
solution of the stochastic control problem (1) by Algorithm 1, where the
corresponding FBSDE is vector-valued, derived from the SMP. In this regard,
our results can be viewed as an extension to [2], where the backward equation
in (10) was only scalar-valued and relied on the HJB equation. Additionally,
we require the following technical, standing assumptions to hold.

######  Assumption 3.

Suppose

  1. 1.

The maps $\bar{b},\bar{\sigma},\bar{F}$ are uniformly
Hölder-$\frac{\sqrt{2}}{2}$ continuous in $t$.

  2. 2.

The feedback map $\mathcal{M}(t,x,p,q)$ is uniformly
Hölder-$\frac{\sqrt{2}}{2}$ continuous in $t$, and Lipschitz continuous in
$x$, $p$ and $q$, respectively.

######  Assumption 4.

There exist constants $k^{b}$ and $k^{F}$, that are possibly negative, such
that

$\displaystyle\left(b(t,x_{1},p)-b(t,x_{2},p)\right)^{\top}\Delta x$
$\displaystyle\leq k^{b}\|\Delta x\|^{2},$ (16)
$\displaystyle\left(F(t,x,p_{1},q)-F(t,x,p_{2},q)\right)^{\top}\Delta p$
$\displaystyle\leq k^{F}\|\Delta p\|^{2}.$

######  Remark 5.

With Assumptions 1, 2 and 3, this implies that $b,\sigma,F$ and $g$ are
uniformly Lipschitz continuous in all spatial variables, and therefore we may
write

$\displaystyle\left\|b(t,x_{1},p_{1})-b(t,x_{2},p_{2})\right\|^{2}$
$\displaystyle\leq L^{b}_{x}\|\Delta x\|^{2}+L^{b}_{p}\|\Delta p\|^{2},$ (17)
$\displaystyle\left\|\sigma(t,x_{1},p_{1})-\sigma(t,x_{2},p_{2})\right\|^{2}$
$\displaystyle\leq L^{\sigma}_{x}\|\Delta x\|^{2}+L^{\sigma}_{p}\|\Delta
p\|^{2},$
$\displaystyle\left\|F\left(t,x_{1},p_{1},q_{1}\right)-F\left(t,x_{2},p_{2},q_{2}\right)\right\|^{2}$
$\displaystyle\leq L^{F}_{x}\|\Delta x\|^{2}+L^{F}_{p}\|\Delta
p\|^{2}+L^{F}_{q}\|\Delta q\|^{2},$
$\displaystyle\left\|\nabla_{x}g(x_{1})-\nabla_{x}g(x_{2})\right\|^{2}$
$\displaystyle\leq L^{g_{x}}_{x}\|\Delta x\|^{2}.$

Similarly, the Hölder-continuity of $b,F$ and $\sigma$ follow from the same
set of assumptions, which imply that $b(t,0,0),F(t,0,0,0)$ and $\sigma(t,0,0)$
are bounded in $t$, and the boundedness of $\nabla_{x}g(x)$ directly follows
from Assumption 1. For convenience, we use $\mathscr{L}$ to denote the set of
all constants mentioned above and denote its upper bound by $L$.

Next, we introduce the following system of quasi-linear parabolic PDEs,
associated with (10),

$\left\\{\begin{aligned}
&\partial_{t}\nu^{i}+\frac{1}{2}\operatorname{Tr}\left(\partial_{xx}\nu^{i}\sigma\sigma^{\top}(t,x,\nu)\right)+\partial_{x}\nu^{i}b\left(t,x,\nu,\partial_{x}\nu\sigma(t,x,\nu)\right)\\\
&\qquad+F^{i}\left(t,x,\nu,\partial_{x}\nu\sigma(t,x,\nu)\right)=0,\quad\forall
i=1,\cdots,d;\\\ &\nu(T,x)=-\nabla_{x}g(x),\end{aligned}\right.$ (18)

where $\nu^{i}$ denotes the $i$-th component of the vector $\nu$. This can be
obtained by Ito’s lemma and the nonlinear Feyman-Kac theorem, for sufficiently
smooth $\nu$.

######  Assumption 5.

The PDE (18) has a classical solution $\nu$ with bounded derivatives
$\partial_{x}\nu$ and $\partial_{xx}^{2}\nu$, and $\sigma$ is bounded.

######  Theorem 2.

Suppose Assumptions 1-5 hold, then the FBSDE (10) has a unique solution
$(X_{t},P_{t},Q_{t})$, and it holds that for $t\in[0,T]$,

$P_{t}=\nu\left(t,X_{t}\right),\quad
Q_{t}=\partial_{x}\nu\left(t,X_{t}\right)\sigma\left(t,X_{t},\nu\left(t,X_{t}\right)\right).$
(19)

######  Remark 6.

This is a direct consequence of the nonlinear Feyman-Kac formula and the proof
can be found in [28, pp. 185-186]. One should note that (18) is not the HJB
equation (11). Instead, the solution of (18) is given by
$\nu(t,X_{t}^{*})=-V_{x}(t,X_{t}^{*})$, when it reduces to the drift control
case, and, in general,

$P_{t}^{*}=-V_{x}(t,X_{t}^{*}),\quad
Q_{t}^{*}=-V_{xx}(t,X_{t}^{*})\bar{\sigma}(t,X_{t}^{*},u_{t}^{*})$ (20)

follows directly from [26, pp. 151-152] or [28, pp. 250-253], see also Remark
4.

With Theorem 2, we can decouple (10) with $\nu$, and therefore standard
arguments and results for decoupled FBSDEs can be used. In particular,
whenever the compositions $b(t,x,\nu(t,x))$ and $\sigma(t,x,\nu(t,x))$ are
Lipschitz in $x$ and Hölder-$1/2$ continuous in $t$, we can derive the
following convergence result for the implicit scheme of the FBSDE.

######  Theorem 3 (Convergence of the implicit scheme).

Suppose Assumptions 1-5 hold, then, for sufficiently small $h$, the following
discrete-time equation $(0\leq i\leq N-1)$

$\left\\{\begin{aligned} &\bar{X}_{0}^{\pi}=x_{0}\\\
&\bar{X}_{t_{i+1}}^{\pi}=\bar{X}_{t_{i}}^{\pi}+b\left(t_{i},\bar{X}_{t_{i}}^{\pi},\bar{P}_{t_{i}}^{\pi}\right)h+\sigma\left(t_{i},\bar{X}_{t_{i}}^{\pi},\bar{P}_{t_{i}}^{\pi}\right)\Delta
W_{i},\\\ &\bar{P}_{T}^{\pi}=-\nabla_{x}g\left(\bar{X}_{T}^{\pi}\right),\\\
&\bar{Q}_{t_{i}}^{\pi}=\frac{1}{h}E\left(\bar{P}_{t_{i+1}}^{\pi}\Delta
W_{i}^{\top}\mid\mathcal{F}_{t_{i}}\right),\\\
&\bar{P}_{t_{i}}^{\pi}=E\left(\bar{P}_{t_{i+1}}^{\pi}+F\left(t_{i},\bar{X}_{t_{i}}^{\pi},\bar{P}_{t_{i}}^{\pi},\bar{Q}_{t_{i}}^{\pi}\right)h\mid\mathcal{F}_{t_{i}}\right),\end{aligned}\right.$
(21)

has a solution,
$\left(\bar{X}_{t_{i}}^{\pi},\bar{P}_{t_{i}}^{\pi},\bar{Q}_{t_{i}}^{\pi}\right)$,
such that $\bar{X}_{t_{i}}^{\pi}\in
L^{2}\left(\Omega,\mathcal{F}_{t_{i}},\mathbb{P}\right)$ and

$\displaystyle\begin{aligned}
\sup_{t\in[0,T]}\Big{(}E\left\|X_{t}-\bar{X}_{t}^{\pi}\right\|^{2}&+E\left\|P_{t}-\bar{P}_{t}^{\pi}\right\|^{2}\Big{)}\\\
&+\int_{0}^{T}E\left\|Q_{t}-\bar{Q}_{t}^{\pi}\right\|^{2}\,dt\end{aligned}\leq
C\left(1+E\|x_{0}\|^{2}\right)h,$ (22)

where
$\bar{X}_{t}^{\pi}=\bar{X}_{t_{i}}^{\pi},\bar{P}_{t}^{\pi}=\bar{P}_{t_{i}}^{\pi},\bar{Q}_{t}^{\pi}=\bar{Q}_{t_{i}}^{\pi}$,
for $t\in\left[t_{i},t_{i+1}\right)$, $\left(X_{t},P_{t},Q_{t}\right)$ is the
solution to (10), and $C$ is a constant depending on $\mathscr{L}$ and $T$.

######  Remark 7.

This result follows from Theorem 5.3.1 and Theorem 5.3.3 of [28], once we
decouple the FBSDE (10) through Theorem 2. In particular, the author has
provided proofs of these results for the explicit scheme in [28], but the same
convergence result holds for the implicit scheme, as stated in their Remark
5.3.2 (ii).

######  Remark 8.

We point out that Theorem 3 here requires conditions that are different from
[2] or [29], as we do not require the weak and monotonicity conditions for the
convergence of the implicit scheme and we rely on a PDE argument instead.

Now, recall the classical Euler scheme (15b). For the discrete equation of
$P_{t_{i+1}}^{\pi}$ in (15b), we take the conditional expectation,
$E\left(\cdot\mid\mathcal{F}_{t_{i}}\right)$, at both sides and obtain the
conditional expectation representation for $P_{t_{i}}^{\pi}$. For the same
discrete equation of $P_{t_{i+1}}^{\pi}$, we multiply by $\left(\Delta
W_{i}\right)^{\top}$ and take $E\left(\cdot\mid\mathcal{F}_{t_{i}}\right)$
afterwards. Therefore, we obtain a formulation, as follows

$\left\\{\begin{aligned} &X_{0}^{\pi}=x_{0},\\\
&X_{t_{i+1}}^{\pi}=X_{t_{i}}^{\pi}+b\left(t_{i},X_{t_{i}}^{\pi},P_{t_{i}}^{\pi}\right)h+\sigma\left(t_{i},X_{t_{i}}^{\pi},P_{t_{i}}^{\pi}\right)\Delta
W_{i},\\\ &Q_{t_{i}}^{\pi}=\frac{1}{h}E\left(P_{t_{i+1}}^{\pi}\Delta
W_{i}^{\top}\mid\mathcal{F}_{t_{i}}\right),\\\
&P_{t_{i}}^{\pi}=E\left(P_{t_{i+1}}^{\pi}+F\left(t_{i},X_{t_{i}}^{\pi},Y_{t_{i}}^{\pi},Z_{t_{i}}^{\pi}\right)h\mid\mathcal{F}_{t_{i}}\right).\end{aligned}\right.$
(23)

######  Remark 9.

A key feature of this formulation is that it does not include the boundary
condition for $P_{T}^{\pi}$, and therefore there are infinitely many
solutions. In particular, it is easy to see that both the classic Euler scheme
(15b) and the implicit scheme (21) are solutions to this formulation. Such a
feature is particularly suitable for the analysis of the algorithm, as we set
our loss function (15a) to measure the distance between $P_{T}^{\pi}$ and
$-\nabla_{x}g(X_{T}^{\pi})$, and one may expect that the closer the two
solutions of (23) are, the smaller loss we will have, and vice versa.

######  Lemma 1.

For $j=1,2$, suppose
$\big{\\{}(X_{t_{i}}^{\pi,j},P_{t_{i}}^{\pi,j},Q_{t_{i}}^{\pi,j})\big{\\}}_{0\leq
i\leq N-1}$ are two solution triples of (23), with
$X_{t_{i}}^{\pi,j},P_{t_{i}}^{\pi,j}\in
L^{2}\left(\Omega,\mathcal{F}_{t_{i}},\mathbb{P}\right),0\leq i\leq N$. For
any $\lambda_{1}>0,\lambda_{2}\geq L^{F}_{q}$, and sufficiently small $h$,
denote

$\displaystyle A_{1}(h)\coloneqq
2k^{b}+\lambda_{1}+L^{\sigma}_{x}+L^{b}_{x}h,$ (24) $\displaystyle
A_{2}(h)\coloneqq\left(\lambda_{1}^{-1}+h\right)L^{b}_{p}+L^{\sigma}_{p},$
$\displaystyle
A_{3}(h)\coloneqq-\frac{\ln\left(1-\left(2k^{F}+\lambda_{2}\right)h\right)}{h},$
$\displaystyle
A_{4}(h)\coloneqq\frac{L^{F}_{x}}{\left(1-\left(2k^{F}+\lambda_{2}\right)h\right)\lambda_{2}}.$

Then, we have, for $0\leq n\leq N$,

$\displaystyle E\left\|X_{t_{n}}^{\pi,1}-X_{t_{n}}^{\pi,2}\right\|^{2}$
$\displaystyle\leq
A_{2}\sum_{i=0}^{n-1}e^{A_{1}(n-i-1)h}E\left\|P_{t_{i}}^{\pi,1}-P_{t_{i}}^{\pi,2}\right\|^{2}h,$
(25) $\displaystyle E\left\|P_{t_{n}}^{\pi,1}-P_{t_{n}}^{\pi,2}\right\|^{2}$
$\displaystyle\leq\begin{aligned}
&e^{A_{3}(N-n)h}E\left\|P_{t_{N}}^{\pi,1}-P_{t_{N}}^{\pi,2}\right\|^{2}\\\
&+A_{4}\sum_{i=n}^{N-1}e^{A_{3}(i-n)h}E\left\|X_{t_{i}}^{\pi,1}-X_{t_{i}}^{\pi,2}\right\|^{2}h.\end{aligned}$
(26)

###### Proof.

We remark that even though $A_{1},A_{2},A_{3},A_{4}$ are all functions of $h$,
in order to ease the presentation, we do not make this dependence explicit.
Let us define

$\displaystyle\begin{aligned} \delta X_{i}&\coloneqq
X_{t_{i}}^{\pi,1}-X_{t_{i}}^{\pi,2},\\\ \delta P_{i}&\coloneqq
P_{t_{i}}^{\pi,1}-P_{t_{i}}^{\pi,2},\\\ \delta b_{i}&\coloneqq
b\left(t_{i},X_{t_{i}}^{\pi,1},P_{t_{i}}^{\pi,1}\right)-b\left(t_{i},X_{t_{i}}^{\pi,2},P_{t_{i}}^{\pi,2}\right),\\\
\delta\sigma_{i}&\coloneqq\sigma\left(t_{i},X_{t_{i}}^{\pi,1},P_{t_{i}}^{\pi,1}\right)-\sigma\left(t_{i},X_{t_{i}}^{\pi,2},P_{t_{i}}^{\pi,2}\right),\\\
\delta F_{i}&\coloneqq
F\left(t_{i},X_{t_{i}}^{\pi,1},P_{t_{i}}^{\pi,1},Q_{t_{i}}^{\pi,1}\right)-F\left(t_{i},X_{t_{i}}^{\pi,2},P_{t_{i}}^{\pi,2},Q_{t_{i}}^{\pi,2}\right).\end{aligned}$
(27)

Then, we have

$\delta X_{i+1}=\delta X_{i}+\delta b_{i}h+\delta\sigma_{i}\Delta W_{i},$ (28)
$\delta P_{i}=E\left(\delta P_{i+1}+\delta
F_{i}h\mid\mathcal{F}_{t_{i}}\right),$ (29)

and motivated by (23), we define

$\delta Q_{i}\coloneqq\frac{1}{h}E\left(\delta P_{i+1}\Delta
W_{i}^{\top}\mid\mathcal{F}_{t_{i}}\right).$ (30)

The martingale representation theorem implies the existence of an
$\mathcal{F}_{t}$-adapted square-integrable process $\left\\{\delta
Q_{t}\right\\}_{t_{i}\leq t\leq t_{i+1}}$, such that

$\delta P_{i+1}=E\left(\delta
P_{i+1}\mid\mathcal{F}_{t_{i}}\right)+\int_{t_{i}}^{t_{i+1}}\delta
Q_{t}\,dW_{t},$ (31)

which, together with (29), gives us

$\delta P_{i+1}=\delta P_{i}-\delta F_{i}h+\int_{t_{i}}^{t_{i+1}}\delta
Q_{t}\,dW_{t}.$ (32)

From (28) and (32), noting that $\delta X_{i},\delta P_{i},\delta
b_{i},\delta\sigma_{i}$, and $\delta F_{i}$ are all $\mathcal{F}_{t_{i}}$
measurable, and $E\left[\Delta W_{i}\mid\mathcal{F}_{t_{i}}\right]=0$,
$E\left(\int_{t_{i}}^{t_{i+1}}Q_{t}\,dW_{t}\mid\mathcal{F}_{t_{i}}\right)=0$,
we have

$\displaystyle E\left\|\delta X_{i+1}\right\|^{2}$
$\displaystyle=E\left\|\delta X_{i}+\delta
b_{i}h\right\|^{2}+hE\left\|\delta\sigma_{i}\right\|^{2},$ (33) $\displaystyle
E\left\|\delta P_{i+1}\right\|^{2}$ $\displaystyle=E\left\|\delta P_{i}-\delta
F_{i}h\right\|^{2}+\int_{t_{i}}^{t_{i+1}}E\left\|\delta
Q_{t}\right\|^{2}\,dt.$ (34)

We first establish the proclaimed upper bound for the forward diffusion part
in (25). Using Assumption 4 and Remark 5, we can apply the root-mean-square
and geometric mean inequality (RMS-GM) and derive, for any $\lambda_{1}>0$,

$E\left\|\delta X_{i+1}\right\|^{2}\begin{aligned} &=\begin{aligned}
&E\left\|\delta X_{i}\right\|^{2}+E\left\|\delta
b_{i}\right\|^{2}h^{2}+hE\left\|\delta\sigma_{i}\right\|^{2}\\\
&+2hE\left(\left(b\left(t_{i},X_{t_{i}}^{\pi,1},P_{t_{i}}^{\pi,1}\right)-b\left(t_{i},X_{t_{i}}^{\pi,2},P_{t_{i}}^{\pi,1}\right)\right)^{\top}\delta
X_{i}\right)\\\
&+2hE\left(\left(b\left(t_{i},X_{t_{i}}^{\pi,2},P_{t_{i}}^{\pi,1}\right)-b\left(t_{i},X_{t_{i}}^{\pi,2},P_{t_{i}}^{\pi,2}\right)\right)^{\top}\delta
X_{i}\right)\\\ \end{aligned}\\\ &\leq\begin{aligned} &E\left\|\delta
X_{i}\right\|^{2}+\left(L^{b}_{x}E\left\|\delta
X_{i}\right\|^{2}+L^{b}_{p}E\left\|\delta P_{i}\right\|^{2}\right)h^{2}\\\
&+\left(L^{\sigma}_{x}E\left\|\delta
X_{i}\right\|^{2}+L^{\sigma}_{p}E\left\|\delta
P_{i}\right\|^{2}\right)h+2k^{b}hE\left\|\delta X_{i}\right\|^{2}\\\
&+\left(\lambda_{1}E\left\|\delta
X_{i}\right\|^{2}+\lambda_{1}^{-1}L^{b}_{p}E\left\|\delta
P_{i}\right\|^{2}\right)h\\\ \end{aligned}\\\ &=\begin{aligned}
&\left(1+\left(2k^{b}+\lambda_{1}+L^{\sigma}_{x}+L^{b}_{x}h\right)h\right)E\left\|\delta
X_{i}\right\|^{2}\\\
&+\left(\left(\lambda_{1}^{-1}+h\right)L^{b}_{p}+L^{\sigma}_{p}\right)E\left\|\delta
P_{i}\right\|^{2}h,\end{aligned}\end{aligned}$ (35)

where we recall the definitions of $A_{1},A_{2}$. Notice that $E\left\|\delta
X_{0}\right\|^{2}=0$, and, thus, by induction, we have that, for $1\leq n\leq
N$,

$\displaystyle E\left\|\delta X_{n}\right\|^{2}$
$\displaystyle\leq\prod_{i=0}^{n-1}(1+A_{1}h)E\left\|\delta
X_{0}\right\|^{2}+\sum_{i=0}^{n-1}(1+A_{1}h)^{n-1-i}A_{2}E\left\|\delta
P_{i}\right\|^{2}h$ (36)
$\displaystyle=\sum_{i=0}^{n-1}(1+A_{1}h)^{n-1-i}A_{2}E\left\|\delta
P_{i}\right\|^{2}h$ $\displaystyle\leq
A_{2}\sum_{i=0}^{n-1}e^{A_{1}(n-i-1)h}E\left\|\delta P_{i}\right\|^{2}h,$

where we used the inequality $(1+x)\leq e^{x}$, $\forall x\in\mathbb{R}$.

In order to show (26), we use a similar approach and apply the RMS-GM
inequality for any $\lambda_{2}>0$, which yields

$\displaystyle E\left\|\delta P_{i+1}\right\|^{2}\begin{aligned} &=E\|\delta
P_{i}\|^{2}-2hE\left(\delta F_{i}^{\top}\delta P_{i}\right)+h^{2}E\|\delta
F_{i}\|^{2}+\int_{t_{i}}^{t_{i+1}}E\left\|\delta Q_{t}\right\|^{2}\,dt\\\
&\geq\begin{aligned} &E\left\|\delta
P_{i}\right\|^{2}+\int_{t_{i}}^{t_{i+1}}E\left\|\delta
Q_{t}\right\|^{2}\,dt\\\
&-2hE\left(\left(F\left(t_{i},X_{i}^{1,\pi},P_{i}^{1,\pi},Q_{i}^{1,\pi}\right)-F\left(t_{i},X_{i}^{1,\pi},P_{i}^{2,\pi},Q_{i}^{1,\pi}\right)\right)^{\top}\delta
P_{i}\right)\\\
&-2hE\left(\left(F\left(t_{i},X_{i}^{1,\pi},P_{i}^{2,\pi},Q_{i}^{1,\pi}\right)-F\left(t_{i},X_{i}^{2,\pi},P_{i}^{2,\pi},Q_{i}^{2,\pi}\right)\right)^{\top}\delta
P_{i}\right)\end{aligned}\\\ &\geq\begin{aligned} &E\left\|\delta
P_{i}\right\|^{2}+\int_{t_{i}}^{t_{i+1}}E\left\|\delta
Q_{t}\right\|^{2}\,dt-2k^{F}hE\left\|\delta P_{i}\right\|^{2}\\\
&-\left(\lambda_{2}E\left\|\delta
P_{i}\right\|^{2}+\lambda_{2}^{-1}\left(L^{F}_{x}E\left\|\delta
X_{i}\right\|^{2}+L^{F}_{q}E\left\|\delta
Q_{i}\right\|^{2}\right)\right)h.\end{aligned}\end{aligned}$ (37)

To deal with the integral term in the last inequality, we derive the following
relation via Ito’s isometry, (32) and (30),

$\delta Q_{i}=\frac{1}{h}E\left(\int_{t_{i}}^{t_{i+1}}\delta
Q_{t}\,dt\mid\mathcal{F}_{t_{i}}\right).$ (38)

Thereafter, the Jensen and the Cauchy-Schwartz inequalities imply the
following lower bound for the integral term

$\displaystyle E\left\|\delta Q_{i}\right\|^{2}h$
$\displaystyle=\sum_{j=1}^{d}\sum_{k=1}^{m}E\left(\left(\delta
Q_{i}\right)_{j,k}\right)^{2}h$ (39)
$\displaystyle=\sum_{j=1}^{d}\sum_{k=1}^{m}\frac{1}{h}E\left(E\left(\int_{t_{i}}^{t_{i+1}}\left(\delta
Q_{t}\right)_{j,k}\,dt\mid\mathcal{F}_{t_{i}}\right)\right)^{2}$
$\displaystyle\leq\sum_{j=1}^{d}\sum_{k=1}^{m}\frac{1}{h}E\left(\int_{t_{i}}^{t_{i+1}}\left(\delta
Q_{t}\right)_{j,k}\,dt\right)^{2}$
$\displaystyle\leq\sum_{j=1}^{d}\sum_{k=1}^{m}\int_{t_{i}}^{t_{i+1}}E\left(\left(\delta
Q_{t}\right)_{j,k}\right)^{2}\,dt$
$\displaystyle=\int_{t_{i}}^{t_{i+1}}E\left\|\delta Q_{t}\right\|^{2}\,dt,$

where $(\cdot)_{j,k}$ denotes the $(j,k)$-entry of the matrix. Combining (37)
with (39), subsequently gives

$E\left\|\delta P_{i+1}\right\|^{2}\geq\begin{aligned}
&\left(1-\left(2k^{F}+\lambda_{2}\right)h\right)E\left\|\delta
P_{i}\right\|^{2}\\\ &+\left(1-L^{F}_{q}\lambda_{2}^{-1}\right)E\left\|\delta
Q_{i}\right\|^{2}h-L^{F}_{x}\lambda_{2}^{-1}E\left\|\delta
X_{i}\right\|^{2}h.\end{aligned}$ (40)

For any $\lambda_{2}\geq L^{F}_{q}$ and sufficiently small $h$ satisfying
$\left(2k^{F}+\lambda_{2}\right)h<1$, we get the following inequality

$E\left\|\delta
P_{i}\right\|^{2}\leq\left(1-\left(2k^{F}+\lambda_{2}\right)h\right)^{-1}\left(E\left\|\delta
P_{i+1}\right\|^{2}+L^{F}_{x}\lambda_{2}^{-1}E\left\|\delta
X_{i}\right\|^{2}h\right).$ (41)

Finally, recalling the definitions of $A_{3},A_{4}$, by induction, we obtain
that, for $0\leq n\leq N-1$,

$E\left\|\delta P_{n}\right\|^{2}\leq e^{A_{3}(N-n)h}E\left\|\delta
P_{N}\right\|^{2}+A_{4}\sum_{i=n}^{N-1}e^{A_{3}(i-n)h}E\left\|\delta
X_{i}\right\|^{2}h.$ (42)

∎

With the help of the apriori result above, we can now state the main result of
this section, which establishes a posteriori error bound for the convergence
of the deep BSDE algorithm for the SMP formulation of the stochastic control
problem under the aforementioned assumptions.

######  Theorem 4.

Suppose Assumptions 1-5 hold true and there exist $\lambda_{1}>$
$0,\lambda_{2}\geq L^{F}_{q}$ and the constants in $\mathscr{L}$ are such that
$\overline{A_{0}}<1$, where

$\displaystyle\overline{A_{1}}$ $\displaystyle:=\lim_{h\to
0}A_{1}(h)=2k^{b}+\lambda_{1}+L^{\sigma}_{x},$ (43)
$\displaystyle\overline{A_{2}}$ $\displaystyle:=\lim_{h\to
0}A_{2}(h)=L^{b}_{p}\lambda_{1}^{-1}+L^{\sigma}_{p},$ (44)
$\displaystyle\overline{A_{3}}$ $\displaystyle:=\lim_{h\to
0}A_{3}(h)=2k^{F}+\lambda_{2},$ (45) $\displaystyle\overline{A_{4}}$
$\displaystyle:=\lim_{h\to 0}A_{4}(h)=L^{F}_{x}\lambda_{2}^{-1},$ (46)
$\displaystyle\overline{A_{0}}$
$\displaystyle:=\overline{A_{2}}\frac{1-e^{-\left(\overline{A_{1}}+\overline{A_{3}}\right)T}}{\overline{A_{1}}+\overline{A_{3}}}\left\\{L^{g_{x}}_{x}e^{\left(\overline{A_{1}}+\overline{A_{3}}\right)T}+\overline{A_{4}}\frac{e^{\left(\overline{A_{1}}+\overline{A_{3}}\right)T}-1}{\overline{A_{1}}+\overline{A_{3}}}\right\\}.$
(47)

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

$\begin{aligned}
\sup_{t\in[0,T]}\Big{(}E\left\|X_{t}-\hat{X}_{t}^{\pi}\right\|^{2}+E\left\|P_{t}-\hat{P}_{t}^{\pi}\right\|^{2}\Big{)}\\\
+\int_{0}^{T}E\left\|Q_{t}-\hat{Q}_{t}^{\pi}\right\|^{2}\,dt\end{aligned}\leq
C\big{(}h+E\left\|-\nabla_{x}g\left(X_{T}^{\pi}\right)-P_{T}^{\pi}\right\|^{2}\big{)},$
(48)

where
$\hat{X}_{t}^{\pi}=X_{t_{i}}^{\pi},\hat{P}_{t}^{\pi}=P_{t_{i}}^{\pi},\hat{Q}_{t}^{\pi}=Q_{t_{i}}^{\pi}$
for $t\in\left[t_{i},t_{i+1}\right)$, and $\left(X_{t},P_{t},Q_{t}\right)$ is
the solution to (10).

###### Proof.

Let
$X_{t_{i}}^{\pi,1}=X_{t_{i}}^{\pi},P_{t_{i}}^{\pi,1}=P_{t_{i}}^{\pi},Q_{t_{i}}^{\pi,1}=$
$Q_{t_{i}}^{\pi}$, given by the Euler scheme (15b), and
$X_{t_{i}}^{\pi,2}=\bar{X}_{t_{i}}^{\pi},P_{t_{i}}^{\pi,2}=\bar{P}_{t_{i}}^{\pi},Q_{t_{i}}^{\pi,2}=\bar{Q}_{t_{i}}^{\pi}$,
given by the implicit scheme (21). Since both of these schemes solve (23), we
can apply Lemma 1. In what follows, we adopt the notation of Lemma 1.

First, using the triangular and Cauchy-Schwartz inequalities, and the RMS-GM
inequality, for any $\lambda_{3}>0$,

$E\left\|\delta
P_{N}\right\|^{2}=E\left\|-\nabla_{x}g\left(\bar{X}_{T}^{\pi}\right)-P_{T}^{\pi}\right\|^{2}\leq\begin{aligned}
&\left(1+\lambda_{3}^{-1}\right)E\left\|-\nabla_{x}g\left(X_{T}^{\pi}\right)-P_{T}^{\pi}\right\|^{2}\\\
&+L^{g_{x}}_{x}\left(1+\lambda_{3}\right)E\left\|\delta
X_{N}\right\|^{2},\end{aligned}$ (49)

where we also used that $\nabla_{x}g$ is Lipschitz continuous.

Let

$\mathcal{X}:=\max_{0\leq n\leq N}e^{-A_{1}nh}E\left\|\delta
X_{n}\right\|^{2},\quad\mathcal{P}:=\max_{0\leq n\leq
N}e^{A_{3}nh}E\left\|\delta P_{n}\right\|^{2}.$ (50)

From the upper bound in (25), it follows that

$e^{-A_{1}nh}E\left\|\delta X_{n}\right\|^{2}\leq
A_{2}\sum_{i=0}^{n-1}e^{-A_{1}(i+1)h}E\left\|\delta P_{i}\right\|^{2}h\leq
A_{2}\mathcal{P}\sum_{i=0}^{n-1}e^{-A_{1}(i+1)h-A_{3}ih}h.$ (51)

Similarly, (26) and (50) imply

$\displaystyle e^{A_{3}nh}E\left\|\delta P_{n}\right\|^{2}$ $\displaystyle\leq
e^{A_{3}T}E\left\|\delta
P_{N}\right\|^{2}+A_{4}\sum_{i=n}^{N-1}e^{A_{3}ih}E\left\|\delta
X_{i}\right\|^{2}h$ (52) $\displaystyle\leq\begin{aligned}
&e^{A_{3}T}\left(1+\lambda_{3}^{-1}\right)E\left\|-\nabla_{x}g\left(X_{T}^{\pi}\right)-P_{T}^{\pi}\right\|^{2}\\\
&+e^{A_{3}T}L^{g_{x}}_{x}\left(1+\lambda_{3}\right)E\left\|\delta
X_{N}\right\|^{2}+A_{4}\sum_{i=n}^{N-1}e^{A_{3}ih}E\left\|\delta
X_{i}\right\|^{2}h\end{aligned}$ $\displaystyle\leq\begin{aligned}
&e^{A_{3}T}\left(1+\lambda_{3}^{-1}\right)E\left\|-\nabla_{x}g\left(X_{T}^{\pi}\right)-P_{T}^{\pi}\right\|^{2}\\\
&+\left(L^{g_{x}}_{x}\left(1+\lambda_{3}\right)e^{\left(A_{1}+A_{3}\right)T}+A_{4}\sum_{i=n}^{N-1}e^{\left(A_{1}+A_{3}\right)ih}h\right)\mathcal{X},\end{aligned}$

where we used (49) to obtain the second inequality. Combining (52), (51) with
(50) thus yields

$\displaystyle\mathcal{X}$ $\displaystyle\leq
A_{2}he^{-A_{1}h}\frac{e^{-\left(A_{1}+A_{3}\right)T}-1}{e^{-\left(A_{1}+A_{3}\right)h}-1}\mathcal{P},$
(53) $\displaystyle\mathcal{P}$ $\displaystyle\leq\begin{aligned}
&e^{A_{3}T}\left(1+\lambda_{3}^{-1}\right)E\left\|-\nabla_{x}g\left(X_{T}^{\pi}\right)-P_{T}^{\pi}\right\|^{2}\\\
&+\left(L^{g_{x}}_{x}\left(1+\lambda_{3}\right)e^{\left(A_{1}+A_{3}\right)T}+A_{4}h\frac{e^{\left(A_{1}+A_{3}\right)T}-1}{e^{\left(A_{1}+A_{3}\right)h}-1}\right)\mathcal{X}.\end{aligned}$
(54)

Now, we define

$A(h)\coloneqq
A_{2}he^{-A_{1}h}\frac{e^{-\left(A_{1}+A_{3}\right)T}-1}{e^{-\left(A_{1}+A_{3}\right)h}-1}\left(L^{g_{x}}_{x}\left(1+\lambda_{3}\right)e^{\left(A_{1}+A_{3}\right)T}+A_{4}h\frac{e^{\left(A_{1}+A_{3}\right)T}-1}{e^{\left(A_{1}+A_{3}\right)h}-1}\right),$
(55)

by which (54) can be rewritten, as follows

$\mathcal{P}\leq
e^{A_{3}T}\left(1+\lambda_{3}^{-1}\right)E\left\|-\nabla_{x}g\left(X_{T}^{\pi}\right)-P_{T}^{\pi}\right\|^{2}+A(h)\mathcal{P}.$
(56)

Hence, whenever $A(h)<1$, the estimates in (53) and (54) take the following
form

$\displaystyle\mathcal{X}$
$\displaystyle\leq\frac{e^{A_{3}T}\left(1+\lambda_{3}^{-1}\right)A_{2}he^{-A_{1}h}\frac{e^{-\left(A_{1}+A_{3}\right)T}-1}{e^{-\left(A_{1}+A_{3}\right)h}-1}E\left\|-\nabla_{x}g\left(X_{T}^{\pi}\right)-P_{T}^{\pi}\right\|^{2}}{1-A(h)},$
(57) $\displaystyle\mathcal{P}$
$\displaystyle\leq\frac{e^{A_{3}T}\left(1+\lambda_{3}^{-1}\right)E\left\|-\nabla_{x}g\left(X_{T}^{\pi}\right)-P_{T}^{\pi}\right\|^{2}}{1-A(h)}.$
(58)

Similarly to (50), we define

$\overline{\mathcal{X}}\coloneqq\max_{0\leq n\leq
N}e^{-\overline{A_{1}}nh}E\left\|\delta
X_{n}\right\|^{2},\quad\overline{\mathcal{P}}\coloneqq\max_{0\leq n\leq
N}e^{\overline{A_{3}}nh}E\left\|\delta P_{n}\right\|^{2}.$ (59)

Recall $\lim_{h\to 0}A_{i}(h)=\overline{A_{i}}$, from (55), it then follows
that

$\lim_{h\rightarrow
0}A(h)=\overline{A_{2}}\frac{1-e^{-\left(\overline{A_{1}}+\overline{A_{3}}\right)T}}{\overline{A_{1}}+\overline{A_{3}}}\left(L^{g_{x}}_{x}\left(1+\lambda_{3}\right)e^{\left(\overline{A_{1}}+\overline{A_{3}}\right)T}+\overline{A_{4}}\frac{e^{\left(\overline{A_{1}}+\overline{A_{3}}\right)T}-1}{\overline{A_{1}}+\overline{A_{3}}}\right).$
(60)

Recall the definition of $\overline{A_{0}}$ in (47). Whenever
$\overline{A_{0}}<1$, comparing (60) with $\overline{A_{0}}$ and using the
definitions in (59), we notice that for any $\epsilon>0$, there exists a
$\lambda_{3}>0$ and sufficiently small $h$, such that

$\displaystyle\overline{\mathcal{X}}$
$\displaystyle\leq(1+\epsilon)\frac{\overline{A_{2}}e^{\overline{A_{3}}T}}{1-\overline{A_{0}}}\left(1+\lambda_{3}^{-1}\right)\frac{1-e^{-\left(\overline{A_{1}}+\overline{A_{3}}\right)T}}{\overline{A_{1}}+\overline{A_{3}}}E\left\|-\nabla_{x}g\left(X_{T}^{\pi}\right)-P_{T}^{\pi}\right\|^{2},$
(61) $\displaystyle\overline{\mathcal{P}}$
$\displaystyle\leq(1+\epsilon)\frac{e^{\overline{A_{3}}T}}{1-\overline{A_{0}}}\left(1+\lambda_{3}^{-1}\right)E\left\|-\nabla_{x}g\left(X_{T}^{\pi}\right)-P_{T}^{\pi}\right\|^{2}.$
(62)

Consequently, by fixing $\epsilon=1$ and choosing a suitable $\lambda_{3}>0$,
we obtain the following error estimates

$\displaystyle\max_{0\leq n\leq N}E\left\|\delta X_{n}\right\|^{2}$
$\displaystyle\leq e^{\overline{A_{1}}T\vee 0}\bar{\mathcal{X}}\leq
C\left(\lambda_{1},\lambda_{2}\right)E\left\|-\nabla_{x}g\left(X_{T}^{\pi}\right)-P_{T}^{\pi}\right\|^{2},$
(63) $\displaystyle\max_{0\leq n\leq N}E\left\|\delta P_{n}\right\|^{2}$
$\displaystyle\leq e^{\left(-\overline{A_{3}}T\right)\vee
0}\bar{\mathcal{P}}\leq
C\left(\lambda_{1},\lambda_{2}\right)E\left\|-\nabla_{x}g\left(X_{T}^{\pi}\right)-P_{T}^{\pi}\right\|^{2}.$
(64)

Finally, in order to estimate $E\left\|\delta Q_{n}\right\|^{2}$, we consider
estimate (40) from the proof of Lemma 1, in which $\lambda_{2}$ can take any
value such that $\lambda_{2}\geq L^{F}_{q}$. In particular, when
$L^{F}_{q}\neq 0$, we choose $\lambda_{2}=2L^{F}_{q}$ and obtain

$\frac{1}{2}E\left\|\delta
Q_{i}\right\|^{2}h\leq\frac{L^{F}_{x}}{2L^{F}_{q}}E\left\|\delta
X_{i}\right\|^{2}h+E\left\|\delta
P_{i+1}\right\|^{2}-\left(1-\left(2k^{F}+2L^{F}_{q}\right)h\right)E\left\|\delta
P_{i}\right\|^{2}.$ (65)

Summing from 0 to $N-1$ therefore gives

$\displaystyle\sum_{i=0}^{N-1}E\left\|\delta Q_{i}\right\|^{2}h$
$\displaystyle\leq\begin{aligned} &\frac{L^{F}_{x}T}{L^{F}_{q}}\max_{0\leq
n\leq N}E\left\|\delta X_{n}\right\|^{2}\\\
&+\left(4\left(k^{F}+L^{F}_{q}\right)T\vee 0+2\right)\max_{0\leq n\leq
N}E\left\|\delta P_{n}\right\|^{2}.\end{aligned}$ (66)

Using the estimates established by (63) and (64), we collect

$\displaystyle\sum_{i=0}^{N-1}E\left\|\delta Q_{i}\right\|^{2}h$
$\displaystyle\leq
C(\lambda_{1},\lambda_{2})E\left\|-\nabla_{x}g\left(X_{T}^{\pi}\right)-P_{T}^{\pi}\right\|^{2}.$
(67)

The case $L^{F}_{q}=0$ can be dealt with similarly by choosing $\lambda_{2}=1$
and the same type of estimate can be derived.

Finally, combining estimates (63), (64) and (67) with Theorem 3, we prove our
statement. ∎

######  Corollary 1.

Let $\left(X_{t},u_{t},P_{t},Q_{t}\right)$ be the optimal 4-tuple that solves
problem (1). Under the setting of Theorem 4, there exists a constant, $C>0$,
depending on $E\|x_{0}\|^{2},\mathscr{L},T,\lambda_{1}$, and $\lambda_{2}$,
such that for sufficiently small $h$,

$\int_{0}^{T}E\|u_{t}-\hat{u}_{t}^{\pi}\|^{2}\,dt\leq
C\big{(}h+E\left\|-\nabla_{x}g\left(X_{T}^{\pi}\right)-P_{T}^{\pi}\right\|^{2}\big{)},$
(68)

where
$\hat{u}_{t_{i}}^{\pi}\coloneqq\mathcal{M}(t,\hat{X}_{t_{i}}^{\pi},\hat{P}_{t_{i}}^{\pi},\hat{Q}_{t_{i}}^{\pi})$,
and $\hat{u}_{t}^{\pi}=\hat{u}_{t_{i}}^{\pi}$ for
$t\in\left[t_{i},t_{i+1}\right)$.

###### Proof.

Recall the optimal feedback map obtained by using the SMP,

$u_{t}=\mathcal{M}(t,X_{t},P_{t},Q_{t}),\quad\forall t\in[0,T],$ (69)

where $\left(X_{t},u_{t},P_{t},Q_{t}\right)$ is the optimal 4-tuple that
solves stochastic control problem (1) and the FBSDE (10), see Remark 4. Then,
by Theorem 4 and the Lipschitz continuity of $\mathcal{M}$ in Assumption 3, we
immediately have

$\displaystyle\int_{0}^{T}E\left\|u_{t}-\hat{u}_{t}^{\pi}\right\|^{2}\,dt$
$\displaystyle=\int_{0}^{T}E\left\|\mathcal{M}(t,X_{t},P_{t},Q_{t})-\mathcal{M}(t,\hat{X}_{t_{i}}^{\pi},\hat{P}_{t_{i}}^{\pi},\hat{Q}_{t_{i}}^{\pi})\right\|^{2}\,dt$
(70) $\displaystyle\leq\begin{aligned}
L\int_{0}^{T}\Big{(}E\left\|X_{t}-\hat{X}_{t}^{\pi}\right\|^{2}&+E\left\|P_{t}-\hat{P}_{t}^{\pi}\right\|^{2}\\\
&+E\left\|Q_{t}-\hat{Q}_{t}^{\pi}\right\|^{2}\Big{)}\,dt\end{aligned}$
$\displaystyle\leq
C\big{(}h+E\left\|-\nabla_{x}g\left(X_{T}^{\pi}\right)-P_{T}^{\pi}\right\|^{2}\big{)}.$

∎

##  5 Numerical Results

In this section, we demonstrate the accuracy and robustness of the proposed
scheme. We consider stochastic control problems, both in the case of drift and
diffusion control. Algorithm 1 has been implemented in TensorFlow 2.15, and
the experiments were run on an Dell Alienware Aurora R10 machine, equipped
with an AMD Ryzen 9 3950X CPU (16 cores, $64\mathrm{Mb}$ cache,
$4.7\mathrm{GHz}$) and an Nvidia GeForce RTX 3090 GPU (24 Gb). The library
used in this paper will be publicly accessible under the github repository of
the second author. In all numerical experiments presented below, we use
standard fully-connected feedforward neural networks to parametrize
$\mu_{0}^{\pi}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ and each
$\phi_{i}^{\pi}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times m}$,
$i=0,\dots,N-1$. Similarly to [22], each neural network has two hidden layers
of width $20$ and a ReLU activation function. The hidden layers have no
biases, as we apply batch normalization after the non-linearity in each hidden
layer. All parameters are initialized to constant zeros. We employ the Adam
optimizer with default settings to minimize the empirical loss function and
use a batch size of $1024$ independent paths of the Brownian motion for each
SGD step during training. An exponential learning rate schedule is applied
starting from $10^{-2}$ with a decay rate of $10^{-2}$. A fixed number of
$2^{14}=16384$ iteration steps are taken. The Monte Carlo errors are computed
over an independent validation sample of $1024$ paths. In order to assess the
inherent randomness of the deep BSDE methods and the SGD iterations, we run
each experiment 5 times and report on the mean and standard deviation of the
resulting independent approximations of Algorithm 1. Computations were carried
out with single floating-point precision. We denote the total approximation
errors by $\delta\hat{X}^{\pi}_{n}$ at time $t_{n}$ and similarly for other
processes. In the comparison with [22], we use $\lambda=1$ as described in
their paper. We denote by $Y_{n}^{\pi}$ the discrete time counterpart of the
objective functional in (1), which we compute similarly to [22, eq. (3.1)] via
a backward summation, using the relations in Remark 6.

###  5.1 Analytic solution for general LQ problems with $m>1$

Our numerical examples are given for so-called linear quadratic (LQ) type
problems, where a reference solution can be derived semi-analytically through
numerically solving a system of ODEs, induced by the HJB equation. In order to
obtain this reference solution, we first extend a result that is well-known in
the literature for the case $m=1$, see e.g. [27]. A stochastic linear
quadratic problem is a special case of (1), where

$\displaystyle\bar{b}(t,x,u)=Ax+Bu+\beta,\
\bar{\sigma}_{j}(t,x,u)=C_{j}x+D_{j}u+\Sigma_{j},$ (71)
$\displaystyle\bar{f}(t,x,u)=\frac{1}{2}(x^{\top}R_{x}x+u^{\top}R_{xu}x+u^{\top}R_{u}u),\
g(x)=\frac{1}{2}x^{\top}Gx,$ (72)

with $A\in\mathbb{R}^{d\times d}$, $B\in\mathbb{R}^{d\times\ell}$,
$\beta\in\mathbb{R}^{d}$, $R_{x}\in\mathbb{R}^{d\times d}$,
$R_{xu}\in\mathbb{R}^{l\times d}$, $R_{u}\in\mathbb{R}^{\ell\times\ell}$ and
$G\in\mathbb{R}^{d\times d}$, $R_{x}$, $R_{u}$ and $G$ are symmetric, and, for
$j=1,\dots,m$, $C_{j}\in\mathbb{R}^{d\times d}$,
$D_{j}\in\mathbb{R}^{d\times\ell}$, $\Sigma_{j}\in\mathbb{R}^{d}$. For
$\bar{\sigma},\Sigma,q\in\mathbb{R}^{d\times m}$, we denote the $j$-th column
by $\bar{\sigma}_{j},\Sigma_{j},q_{j}$. Utilizing the SMP, we obtain the
feedback map for the optimal control

$u=\mathcal{M}(t,x,p,q)=-R_{u}^{-1}\big{(}R_{xu}x-B^{\top}p-\sum_{j=1}^{m}D_{j}^{\top}q_{j}\big{)},$
(73)

and therein we obtain the corresponding FBSDE (10), by substituting (73) back
into the forward diffusion and adjoint equation.

Next, we derive the reference solution for this general formulation, following
the same approach as in [27], but extending to the case $m>1$. By the DP, the
value function $V(t,x)$ for a LQ problem should satisfy the HJB equation (11)
with a boundary condition $V(T,x)=\frac{1}{2}x^{\top}Gx$, and, in the LQ case,
the generalized Hamiltonian $\mathcal{G}$ is given by

$\displaystyle\mathcal{G}\left(t,x,u,-V_{x},-V_{xx}\right)=$
$\displaystyle-\frac{1}{2}u^{\top}R_{u}u-u^{\top}R_{xu}x-\frac{1}{2}x^{\top}R_{x}x-(Ax+Bu+\beta)^{\top}V_{x}$
(74)
$\displaystyle-\sum_{j=1}^{m}\frac{1}{2}(C_{j}x+D_{j}u+\Sigma_{j})^{\top}V_{xx}(C_{j}x+D_{j}u+\Sigma_{j}).$

We conjecture that the value function takes the form
$V(t,x)=\frac{1}{2}x^{\top}\Gamma(t)x+x^{\top}\gamma(t)+\kappa(t)$, for some
unknowns symmetric $\Gamma(t)\in\mathbb{R}^{d\times d}$,
$\gamma(t)\in\mathbb{R}^{d}$ and $\kappa(t)\in\mathbb{R}$, and for simplicity
we omit the argument $t$ in the following. Notice that $\mathcal{G}$ is
quadratic in $u$, then, by a completion of square technique,

$\displaystyle\mathcal{G}(t,x,u,-V_{x},-V_{xx})=\begin{aligned}
&-\frac{1}{2}(u+\Psi x+\psi)^{\top}\widehat{R}(u+\Psi
x+\psi)+\frac{1}{2}x^{\top}\widehat{S}^{\top}\widehat{R}^{-1}\widehat{S}x\\\
&+\frac{1}{2}x^{\top}\left(-R_{x}-\sum_{j=1}^{m}C_{j}^{\top}\Gamma
C_{j}-\Gamma^{\top}A-A^{\top}\Gamma\right)x\\\
&-x^{\top}\left(A^{\top}\gamma+\sum_{j=1}^{m}C_{j}^{\top}\Gamma\Sigma_{j}-\Psi^{\top}\widehat{R}\psi+\Gamma\beta\right)\\\
&+\frac{1}{2}\psi^{\top}\widehat{R}\psi-\beta^{\top}\gamma-\sum_{j=1}^{m}\frac{1}{2}\Sigma_{j}^{\top}\Gamma\Sigma_{j},\end{aligned}$
(75)

where we define the following

$\displaystyle\widehat{R}\coloneqq R_{u}+\sum_{j=1}^{m}D_{j}^{\top}\Gamma
D_{j},\quad\widehat{S}\coloneqq
B^{\top}\Gamma+R_{xu}+\sum_{j=1}^{m}D_{j}^{\top}\Gamma C_{j},$ (76)
$\displaystyle\Psi\coloneqq\widehat{R}^{-1}\widehat{S},\quad\psi\coloneqq\widehat{R}^{-1}\left(B^{\top}\gamma+\sum_{j=1}^{m}D_{j}^{\top}\Gamma\Sigma_{j}\right),$

provided that $\widehat{R}$ is positive definite. Due to the first quadratic
term, we immediately see that the optimal control should take the following
feedback form,

$u^{*}=-\Big{(}R_{u}+\sum_{j=1}^{m}D_{j}^{\top}V_{xx}D_{j}\Big{)}^{-1}\begin{aligned}
\Big{(}&B^{\top}V_{x}+\big{(}R_{xu}+\sum_{j=1}^{m}D_{j}^{\top}V_{xx}C_{j}\big{)}x\\\
&+\sum_{j=1}^{m}D_{j}^{\top}V_{xx}\Sigma_{j}\Big{)}.\end{aligned}$ (77)

Substituting $\mathcal{G}(t,x,u^{*},-V_{x},-V_{xx})$ back to (11), we find
that $V(t,x)=\frac{1}{2}x^{\top}\Gamma(t)x+x^{\top}\gamma(t)+\kappa(t)$ solves
the HJB equation whenever the following relations hold, at each $t\in[0,T]$

$\left\\{\begin{aligned} &0=\dot{\Gamma}+\Gamma
A+A^{\top}\Gamma+\sum_{j=1}^{m}C_{j}^{\top}\Gamma
C_{j}+R_{x}-\widehat{S}^{\top}\widehat{R}^{-1}\widehat{S},\\\
&0=\dot{\gamma}+A^{\top}\gamma+\sum_{j=1}^{m}C_{j}^{\top}\Gamma\Sigma_{j}-\Psi^{\top}\widehat{R}\psi+\Gamma\beta,\\\
&0=\dot{\kappa}-\frac{1}{2}\psi^{\top}\widehat{R}\psi+\beta^{\top}\gamma+\sum_{j=1}^{m}\frac{1}{2}\Sigma_{j}^{\top}\Gamma\Sigma_{j},\\\
&\Gamma(T)=G,\quad\gamma(T)=0,\quad\kappa(T)=0.\end{aligned}\right.$ (78)

This system of ODEs can be integrated numerically with practically arbitrary
accuracy. With the thereby obtained numerical solution of (78), one can
subsequently compute the value function and all of its derivatives using the
conjecture, and obtain the optimal control $u^{*}_{t}$ through the feedback
map (77). The reference solution to the corresponding BSDE (10), derived via
the SMP, can similarly be computed using the relations (20) in Remark 6. In
what follows, we compare our approximations to a reference solution obtained
through the numerical integration of the Ricatti ODEs in (77), using an
equdistant time grid with $N_{\text{ode}}=10^{5}$, and simulating the coupled
FBSDE system in (10) on the same time grid with a discrete Euler-Maruyama
approximation.

###  5.2 Example 1 – drift control

The following stochastic control problem can be found in [22, sec. 5.1.2]. The
control space is $\ell=2$ dimensional, whereas for the state space $d=m=6$.
The coefficients in (1) are of LQ type and the corresponding matrices in (71)
are defined, as follows,

$\displaystyle A$ $\displaystyle=-\text{diag}([1,2,3,1,2,3]),$ (79)
$\displaystyle B$ $\displaystyle=\begin{pmatrix}&1,&1,&0.5,&1,&0,&0\\\
&-1,&1,&1,&-1,&-1,&1\end{pmatrix}^{\top},$ (80) $\displaystyle\beta$
$\displaystyle=-A([-0.2,-0.1,0,0,0.1,0.2])^{\top},\quad C_{j}=0,\quad
D_{j}=0,$ (81) $\displaystyle\Sigma$
$\displaystyle=\text{diag}([0.05,0.25,0.05,0.25,0.05,0.25]),\quad
R_{xu}=0,\quad R_{u}=2I_{l},$ (82) $\displaystyle R_{x}$
$\displaystyle=2\text{diag}([25,1,25,1,25,1]),\quad
G=2\text{diag}([1,25,1,25,1,25]),$ (83)

for all $j=1,\dots,m$. We consider a deterministic initial condition,
$x_{0}=(0.1,\dots,0.1)$, and a terminal time of $T=1/2$. We remark that this
problem falls under the class of drift control problems, and, recalling the
discussion following Assumption 2, it in particular satisfies the conditions
of Theorem 4 and Corollary 1 over any compact subset of the state space.

(a) Algorithm 1.

(b) Andersson et al. (2023) in [22].

Figure 1: Example 1, $N=50$. Errors of approximations of the optimally
controlled state space and control strategy. On top: results obtained through
Algorithm 1 and the SMP. On the bottom: reference method from [22] using
dynamic programming. Lines correspond to the mean of $5$ independent runs of
the algorithm, shaded areas to the standard deviation. Graphs computed over an
independent Monte Carlo sample of size $M=2^{10}$.

In Figure 1, the relative approximation errors of the optimal control strategy
and the optimally controlled state space are compared to the method proposed
in [22], explained in Section 3. We remark that, as found in the
aforementioned paper, the deep BSDE method on the dynamic programming equation
fails to converge to the true solution. As it can be seen our deep BSDE
formulation, via the SMP in Algorithm 1, leads to more accurate and robust
approximations of both the controlled diffusion and the optimal control. The
relative approximation errors of both processes are below
$\mathcal{O}(10^{-2})$ for over $95\%$ of the sampling paths. Even though,
standard to the FBSDE literature, our theory only establishes error bounds in
the $L^{2}$ sense, see (48) and (68), the tight quantile estimates in Figure 1
suggest that convergence is achieved in a stronger sense as well. Increasing
the number of discretization points $N$ further reduces the errors.

The convergence of proposed algorithm is collected in Figure 2(a), where the
absolute errors of all corresponding processes are plotted for different
values of $N$. We recall that the insights of Theorem 4 suggest that these
errors admit to an upper bound depending on $h=T/N$ and the value of the
objective functional
$E||-\nabla_{x}g(\hat{X}_{N}^{\pi})-\hat{P}_{N}^{\pi}||^{2}$. Indeed, Figure
2(a) confirms our theoretical findings. The approximation errors of all
processes decay with a rate of at least $\mathcal{O}(h)$, whenever the
dominating term of the total error is the discretization part. In particular,
as can be seen from the right end of the curves, once $N$ is large, and the
optimization problem in Algorithm 1 also gets more complex, our fixed
optimization strategy only manages to yield a terminal loss of the same order
as the discretization, leading to stalled convergence. In order to preserve
convergence for very fine time grids, one needs to make sure the terminal loss
is sufficiently small, as implied by Theorem 4.

The relative $L^{2}$ approximation errors over time are depicted in Figure
3(a), for $N=100$. The method yields highly accurate approximations for all
time steps, and of a similar magnitude for all processes. It is worth to
mention the accuracy at $t=0$, where no discretization error from the
reference solution is present. All error measures are collected in Table 2. In
line with the figures, we see that convergence is achieved in the natural
norms of all processes. The resulting approximations are robust regardless of
the underlying randomness of the Monte Carlo machinery, standard deviations of
the errors are orders of magnitude smaller than their corresponding means. As
shown in the last column, Algorithm 1 also offers a competitive runtime for
high-dimensional problems, as an SGD iteration step takes approximately $0.16$
seconds on average when $N=100$. In our experiments, this beats the runtime of
[22] due to the lack of the extra backward summation step considered therein.

(a) Example 1.

(b) Example 2.

Figure 2: Convergence and empirical convergence rates of Algorithm 1 over $N$.
Lines correspond to the mean of 5 independent runs of the algorithm, error
bars to the standard deviation. Graphs compute over an independent Monte Carlo
sample of size $M=2^{10}$.

(a) Example 1.

(b) Example 2.

Figure 3: Relative $L^{2}$ approximation errors over time, $N=100$. Lines
correspond to the mean of 5 independent runs of the algorithm, shaded areas to
the standard deviation. Graphs computed over an independent Monte Carlo sample
of size $M=2^{10}$.

Table 1: Example 1. Error measures for different values of $N$. Figures
correspond to the mean (standard deviation) of $5$ independent runs of the
algorithm. Expectations computed over an independent Monte Carlo sample of
size $M=2^{10}$.

| max. $E||\delta\hat{X}_{n}^{\pi}||^{2}$ | max. $E||\delta\hat{P}_{n}^{\pi}||^{2}$ | avg. $E||\delta\hat{Q}_{n}^{\pi}||^{2}$ | avg. $E||\delta\hat{u}_{n}^{\pi}||^{2}$ | $E||-\nabla_{x}g(\hat{X}_{N}^{\pi})-\hat{P}_{N}^{\pi}||^{2}$ | $|\delta\hat{Y}_{0}^{\pi}|$ | runtime (s)
---|---|---|---|---|---|---|---
1 | $$8.660010-02$($510-06$)$ | $$1.727\,5410+02$($710-03$)$ | $$7.510310+01$($310-03$)$ | $$1.636010-01$($210-05$)$ | $$310-05$($110-05$)$ | $$6.410+01$($310+00$)$ | $$7.810+01$($110+00$)$
5 | $$9.11710-03$($510-06$)$ | $$1.538410+01$($110-03$)$ | $$2.224610+01$($210-03$)$ | $$2.78710-02$($510-05$)$ | $$1.310-04$($810-05$)$ | $$6.810-01$($510-02$)$ | $$2.29810+02$($510-01$)$
10 | $$3.53510-03$($310-06$)$ | $$5.177210+00$($710-04$)$ | $$9.90110+00$($210-03$)$ | $$1.18410-02$($110-05$)$ | $$110-04$($110-04$)$ | $$1.2410-01$($710-03$)$ | $$4.1610+02$($410+00$)$
20 | $$1.63310-03$($410-06$)$ | $$1.499110+00$($410-04$)$ | $$3.68410+00$($110-03$)$ | $$5.0410-03$($210-05$)$ | $$210-04$($210-04$)$ | $$2.410-02$($310-03$)$ | $$7.0610+02$($510+00$)$
50 | $$1.20010-03$($510-06$)$ | $$2.47610-01$($210-04$)$ | $$7.77110-01$($210-04$)$ | $$1.45010-03$($710-06$)$ | $$510-04$($310-04$)$ | $$3.410-03$($510-04$)$ | $$1.30510+03$($710+00$)$
100 | $$1.23910-03$($510-06$)$ | $$6.4010-02$($910-04$)$ | $$2.12010-01$($610-04$)$ | $$4.7910-04$($410-06$)$ | $$4.310-03$($910-04$)$ | $$810-04$($110-04$)$ | $$2.6810+03$($210+01$)$

Table 2: Example 2. Error measures for different values of $N$. Figures
correspond to the mean (standard deviation) of $5$ independent runs of the
algorithm. Expectations computed over an independent Monte Carlo sample of
size $M=2^{10}$.

| max. $E||\delta\hat{X}_{n}^{\pi}||^{2}$ | max. $E||\delta\hat{P}_{n}^{\pi}||^{2}$ | avg. $E||\delta\hat{Q}_{n}^{\pi}||^{2}$ | avg. $E||\delta\hat{u}_{n}^{\pi}||^{2}$ | $E||-\nabla_{x}g(\hat{X}_{N}^{\pi})-\hat{P}_{N}^{\pi}||^{2}$ | $|\delta\hat{Y}_{0}^{\pi}|$ | runtime (s)
---|---|---|---|---|---|---|---
1 | $$9.217310-02$($910-06$)$ | $$1.822\,2910+02$($310-03$)$ | $$7.984910+01$($210-03$)$ | $$1.550510-01$($210-05$)$ | $$310-05$($210-05$)$ | $$6.610+01$($310+00$)$ | $$7.8310+01$($610-01$)$
5 | $$1.012610-02$($210-06$)$ | $$1.70510+01$($110-02$)$ | $$2.39010+01$($110-02$)$ | $$2.84110-02$($410-05$)$ | $$4.06310-01$($310-04$)$ | $$7.910-01$($510-02$)$ | $$2.2910+02$($210+00$)$
10 | $$3.76710-03$($710-06$)$ | $$5.73810+00$($710-03$)$ | $$1.05410+01$($110-02$)$ | $$1.27710-02$($410-05$)$ | $$4.9110-01$($210-03$)$ | $$1.510-01$($310-02$)$ | $$4.20010+02$($910-01$)$
20 | $$1.91410-03$($910-06$)$ | $$1.81810+00$($310-03$)$ | $$3.94610+00$($810-03$)$ | $$6.1910-03$($210-05$)$ | $$6.1310-01$($210-03$)$ | $$3.310-02$($210-03$)$ | $$7.0610+02$($610+00$)$
50 | $$1.60710-03$($410-06$)$ | $$5.4710-01$($710-03$)$ | $$1.05210+00$($610-03$)$ | $$2.79610-03$($110-05$)$ | $$8.0210-01$($510-03$)$ | $$410-03$($110-03$)$ | $$1.29910+03$($910+00$)$
100 | $$1.7110-03$($210-05$)$ | $$4.5810-01$($710-03$)$ | $$6.3810-01$($510-03$)$ | $$2.2010-03$($210-05$)$ | $$9.510-01$($110-02$)$ | $$2.410-03$($410-04$)$ | $$2.6910+03$($110+01$)$

###  5.3 Example 2 – drift and diffusion control

In the next section, we provide an extension to the previous problem by
introducing diffusion control to the LQ problem in Section 5.2. The
coefficients read,

$\displaystyle C_{j}$
$\displaystyle=\frac{1}{60}\text{diag}([1,2,3,1,2,3]),\quad D_{j}$
$\displaystyle=\frac{1}{60}\begin{pmatrix}&1,&0,&1,&0,&1,&0\\\
&0,&-1,&0,&-1,&0,&-1\end{pmatrix}^{T},$ (84)

for all $j=1,\dots,m$, and the rest as in (79). Note that due to the presence
of diffusion control, this equation can not be treated via standard DP, such
as [22], therefore no reference method is provided.

In Figure 4, the relative approximation errors of the controlled state process
and the optimal control strategy are depicted. Even with the introduction of
diffusion control, we obtain accurate approximations up to
$\mathcal{O}(10^{-2})$ relative accuracy in both processes. The tight quantile
bounds suggest that our approximations remain accurate even in stronger senses
than $L^{2}$.

As can easily be seen from (73), $Q$ couples back into the forward diffusion
and hence the conditions of Assumption 2 are not satisfied. Therefore, this
example does not fall under the setting of Theorem 4. Nevertheless, the
convergence drawn in Figure 2(b) gives hope that a similar theoretical bound
may be established, even in the case of diffusion control. As we see, the
empirical convergence rates of $P,Q,u$ are all $\mathcal{O}(h)$. However, once
the loss functional, corresponding to the second term in (48), starts to
dominate the total approximation error, the convergence stalls. This
observation is comparable to Example 1.

The relative $L^{2}$ approximation errors are depicted in Figure 3(b) for
$N=100$. As we see, even in the case of diffusion control, we preserve a
similar relative approximation error of $\mathcal{O}(10^{-2})$, as in Example
1. The slight increase in the relative error of $u$ is due to the feedback map
(73), also accumulating the errors in $Q$. Finally, all error measures are
collected in Table 2. We find that even though convergence in $P,Q$ and $u$
are achieved, the value of the loss functional does not decrease sufficiently
as the discretization is refined. This motivates the fact that improvements to
our fixed optimization strategy need to be made, in order to be able to tackle
more complicated diffusion control problems. Nevertheless, the lack of
difference in total runtimes in Table 2, compared to the drift control case in
2, highlights the great potential of such deep BSDE formulations in the
framework of high-dimensional diffusion control problems.

Figure 4: Example 2, $N=100$. Errors of the approximations of the optimally
controlled state space and control strategy through Algorithm 1 and the SMP.
Lines correspond to the mean of $5$ independent runs of the algorithm, shaded
areas to the standard deviation. Graphs computed over an independent Monte
Carlo sample of size $M=2^{10}$.

## Acknowledgement

The first author would like to thank the China Scholarship Council (CSC) for
the financial support. The second author acknowledges financial support from
the Peter Paul Peterich Foundation via the TU Delft University Fund.

## References

  * [1] Shaolin Ji, Shige Peng, Ying Peng, and Xichuan Zhang.  Solving stochastic optimal control problem via stochastic maximum principle with deep learning method.  Journal of Scientific Computing, 93(1):30, 2022.
  * [2] Jiequn Han and Jihao Long.  Convergence of the deep BSDE method for coupled FBSDEs.  Probability, Uncertainty and Quantitative Risk, 5:1–33, 2020.
  * [3] Richard Bellman.  Dynamic programming and stochastic control processes.  Information and control, 1(3):228–239, 1958.
  * [4] Jean-Michel Bismut.  An introductory approach to duality in optimal stochastic control.  SIAM review, 20(1):62–78, 1978.
  * [5] Lev Semenovich Pontryagin.  Mathematical theory of optimal processes.  Routledge, 2018.
  * [6] Harold J Kushner.  Numerical methods for stochastic control problems in continuous time.  SIAM Journal on Control and Optimization, 28(5):999–1048, 1990\.
  * [7] Nicolai V. Krylov.  The Rate of Convergence of Finite-Difference Approximations for Bellman Equations with Lipschitz Coefficients.  Applied Mathematics and Optimization, 52(3):365–399, October 2005\.
  * [8] Hongjie Dong and Nicolai V Krylov.  The rate of convergence of finite-difference approximations for parabolic bellman equations with lipschitz coefficients in cylindrical domains.  Applied Mathematics and Optimization, 56:37–66, 2007.
  * [9] Espen Robstad Jakobsen.  On the rate of convergence of approximation schemes for Bellman equations associated with optimal stopping time problems.  Mathematical Models and Methods in Applied Sciences, 13(05):613–644, 2003.
  * [10] Jiequn Han et al.  Deep learning approximation for stochastic control problems.  arXiv preprint arXiv:1611.07422, 2016.
  * [11] Kenneth J Hunt, D Sbarbaro, R Żbikowski, and Peter J Gawthrop.  Neural networks for control systems—a survey.  Automatica, 28(6):1083–1112, 1992.
  * [12] Charles-Albert Lehalle and Robert Azencott.  Piecewise affine neural networks and nonlinear control.  In ICANN 98: Proceedings of the 8th International Conference on Artificial Neural Networks, Skövde, Sweden, 2–4 September 1998 8, pages 633–638. Springer, 1998.
  * [13] D. Psaltis, A. Sideris, and A.A. Yamamura.  A multilayered neural network controller.  IEEE Control Systems Magazine, 8(2):17–21, 1988.
  * [14] Achref Bachouch, Côme Huré, Nicolas Langrené, and Huyen Pham.  Deep neural networks algorithms for stochastic control problems on finite horizon: numerical applications.  Methodology and Computing in Applied Probability, 24(1):143–178, 2022.
  * [15] Côme Huré, Huyên Pham, Achref Bachouch, and Nicolas Langrené.  Deep neural networks algorithms for stochastic control problems on finite horizon: convergence analysis.  SIAM Journal on Numerical Analysis, 59(1):525–557, 2021.
  * [16] Marcus Pereira, Ziyi Wang, Tianrong Chen, Emily Reed, and Evangelos Theodorou.  Feynman-Kac neural network architectures for stochastic control using second-order FBSDE theory.  In Learning for Dynamics and Control, pages 728–738. PMLR, 2020\.
  * [17] Jiequn Han, Arnulf Jentzen, and Weinan E.  Solving high-dimensional partial differential equations using deep learning.  Proceedings of the National Academy of Sciences, 115(34):8505–8510, 2018.
  * [18] Jean-Pierre Fouque and Zhaoyu Zhang.  Deep learning methods for mean field control problems with delay.  Frontiers in Applied Mathematics and Statistics, 6:11, 2020.
  * [19] René Carmona and Mathieu Laurière.  Convergence analysis of machine learning algorithms for the numerical solution of mean field control and games: II—the finite horizon case.  The Annals of Applied Probability, 32(6):4065–4105, 2022.
  * [20] Ruimeng Hu and Mathieu Laurière.  Recent developments in machine learning methods for stochastic control and games.  arXiv preprint arXiv:2303.10257, 2023.
  * [21] Maximilien Germain, Huyên Pham, Xavier Warin, et al.  Neural networks-based algorithms for stochastic control and PDEs in finance.  arXiv preprint arXiv:2101.08068, 2021.
  * [22] Kristoffer Andersson, Adam Andersson, and Cornelis W Oosterlee.  Convergence of a robust deep FBSDE method for stochastic control.  SIAM Journal on Scientific Computing, 45(1):A226–A255, 2023.
  * [23] Côme Huré, Huyên Pham, and Xavier Warin.  Deep backward schemes for high-dimensional nonlinear PDEs.  Mathematics of Computation, 89(324):1547–1579, 2020.
  * [24] Balint Negyesi, Kristoffer Andersson, and Cornelis W Oosterlee.  The One Step Malliavin scheme: new discretization of BSDEs implemented with deep learning regressions.  to appear in: IMA Journal of Numerical Analysis, 2024.
  * [25] Chengfan Gao, Siping Gao, Ruimeng Hu, and Zimu Zhu.  Convergence of the backward deep BSDE method with applications to optimal stopping problems.  SIAM J. Financial Mathematics, 14(4):1290–1303, 2023.
  * [26] Huyên Pham.  Continuous-time stochastic control and optimization with financial applications, volume 61.  Springer Science & Business Media, 2009.
  * [27] Jiongmin Yong and Xun Yu Zhou.  Stochastic controls: Hamiltonian systems and HJB equations, volume 43.  Springer Science & Business Media, 1999.
  * [28] Jianfeng Zhang and Jianfeng Zhang.  Backward stochastic differential equations.  Springer, 2017.
  * [29] Christian Bender and Jianfeng Zhang.  Time discretization and Markovian iteration for coupled FBSDEs.  The Annals of Applied Probability, 18(1):143–177, February 2008\.