Paper deep dive
A Posteriori Error Analysis for Decoupled Neural Approximations of Fully Coupled FBSDEs with Control Mismatch
Xichuan Zhang
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 98%
Last extracted: 7/5/2026, 2:04:08 AM
Summary
The paper develops an a posteriori error analysis framework for decoupled neural approximations of fully coupled forward-backward stochastic differential equations (FBSDEs). It addresses the 'control mismatch' that occurs when an auxiliary control process is used in the forward coefficients instead of the actual backward component. The authors establish a continuous-time stability estimate and transfer it to a discrete-time setting, deriving computable error bounds based on terminal defect, pathwise residual, and control mismatch. The framework is validated through numerical experiments on linear-quadratic and Burgers-type FBSDEs.
Entities (7)
Relation Signals (3)
Xichuan Zhang → affiliatedwith → Intelligent Game and Decision Lab
confidence 100% · Xichuan Zhang Intelligent Game and Decision Lab, Beijing 100091, China
Control Mismatch → affects → Deep BSDE Method
confidence 90% · This decoupling is useful in practical deep learning implementations, but it creates a control mismatch that must be included in the error analysis.
Deep BSDE Method → approximates → Fully Coupled FBSDE
confidence 90% · By reformulating the solution of an FBSDE as a stochastic optimization problem and parameterizing the unknown control process with neural networks, the deep BSDE method provides an effective approach...
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:This paper develops an a posteriori error analysis framework for decoupled neural approximations of fully coupled forward--backward stochastic differential equations (FBSDEs). It provides an a posteriori error-analysis for the idealized discrete adapted trajectory. The main feature of the proposed formulation is the use of an auxiliary control process in the forward coefficients, which may differ from the backward component approximated by the neural network. This decoupling is useful in practical deep learning implementations, but it creates a control mismatch that must be included in the error analysis. We first establish a continuous-time stability estimate for fully coupled FBSDEs under perturbations of the drift, diffusion, generator, terminal condition, and auxiliary control input. We then transfer this estimate to the discrete-time setting and derive computable a posteriori error bounds depending only on the terminal defect, the pathwise residual, and the control mismatch. When the auxiliary control is identified with the backward approximation, the mismatch term vanishes and the bound reduces to the standard two-term form. Numerical experiments on a linear--quadratic FBSDE with an explicit reference solution and a multidimensional Burgers-type FBSDE without a reference solution illustrate the diagnostic role of the proposed indicators and the contribution of the mismatch penalty to the consistency and reproducibility of the numerical approximations.
Tags
Links
- Source: https://arxiv.org/abs/2606.29474v1
- Canonical: https://arxiv.org/abs/2606.29474v1
Trouble viewing inline? Open PDF directly →
Full Text
96,441 characters extracted from source content.
Expand or collapse full text
A Posteriori Error Analysis for Decoupled Neural Approximations of Fully Coupled FBSDEs with Control Mismatch Xichuan Zhang Intelligent Game and Decision Lab, Beijing 100091, China Abstract This paper develops an a posteriori error analysis framework for decoupled neural approximations of fully coupled forward–backward stochastic differential equations (FBSDEs). It provides an a posteriori error-analysis for the idealized discrete adapted trajectory. The main feature of the proposed formulation is the use of an auxiliary control process in the forward coefficients, which may differ from the backward component approximated by the neural network. This decoupling is useful in practical deep learning implementations, but it creates a control mismatch that must be included in the error analysis. We first establish a continuous-time stability estimate for fully coupled FBSDEs under perturbations of the drift, diffusion, generator, terminal condition, and auxiliary control input. We then transfer this estimate to the discrete-time setting and derive computable a posteriori error bounds depending only on the terminal defect, the pathwise residual, and the control mismatch. When the auxiliary control is identified with the backward approximation, the mismatch term vanishes and the bound reduces to the standard two-term form. Numerical experiments on a linear–quadratic FBSDE with an explicit reference solution and a multidimensional Burgers-type FBSDE without a reference solution illustrate the diagnostic role of the proposed indicators and the contribution of the mismatch penalty to the consistency and reproducibility of the numerical approximations Keywords fully coupled FBSDE, deep BSDE method, a posteriori error estimate, control mismatch, deep neural networks MSC codes 65C30, 60H35, 68T07, 65M15 1 Introduction Backward stochastic differential equations (BSDEs) were systematically introduced by Pardoux and Peng [P90], who proved the existence and uniqueness of nonlinear BSDEs under Lipschitz conditions. For fully coupled forward–backward stochastic differential equations (FBSDEs), in which the forward drift and diffusion may depend on the backward component Y and the control process Z, fundamental well-posedness results were subsequently established by Peng and Wu [HP95, PW99a] and Ma and Yong [MY99]. These equations provide a probabilistic representation for a wide range of problems in stochastic control, mathematical finance, recursive utility, and quasilinear parabolic partial differential equations (PDEs) through nonlinear Feynman–Kac type formulas [EQ97, EPQ97, PEN90, DE92, PEN91]. Despite their theoretical importance, the numerical solution of FBSDEs remains challenging, especially in high dimensions. Classical grid-based methods, such as finite difference and finite element methods [MPY94, BT04, BET04, BZ08, FZZ16, BS12, GL08], suffer from the curse of dimensionality [BEL57]. Regression-based Monte Carlo methods and multistep schemes have greatly improved the numerical treatment of BSDEs and FBSDEs, but their implementation and accuracy may still deteriorate rapidly as the dimension increases [GLW05, GT16]. A major development in this direction is the deep BSDE method proposed by Han, Jentzen and E [HJE18]. By reformulating the solution of an FBSDE as a stochastic optimization problem and parameterizing the unknown control process with neural networks, the deep BSDE method provides an effective approach for high-dimensional nonlinear PDEs and stochastic systems. Since then, a growing literature has developed deep learning-based methods for BSDEs, FBSDEs, and related PDEs [RPK19, EHJ17, HPB+21, BHL+19, JPP+20, GMW22]. On the theoretical side, Han and Long [HL20] established a posteriori error estimates and convergence results for deep BSDE methods under suitable approximation assumptions. Related developments include extensions to FBSDEs with jumps, optimal stopping problems, multi-FBSDE formulations, and mean-field or McKean–Vlasov type systems [GOP25, GCZ+23, AAO25, RSZ24]. The closest work to the present paper is the a posteriori error analysis of fully coupled McKean–Vlasov FBSDEs by Reisinger, Stockinger and Zhang [RSZ24]. Their framework treats a more general class of law-dependent fully coupled systems and establishes residual-type a posteriori estimates for time-discretized McKean–Vlasov FBSDEs. Since a classical fully coupled FBSDE can be regarded as a special case without law dependence, their results are closely related to the setting considered here; however, they do not address the control mismatch that arises from decoupled neural parametrizations. Previous work [JPP+20] proposed three deep learning algorithms for high-dimensional fully coupled FBSDEs. The second algorithm therein introduced an auxiliary control process U and its training loss already included a penalty on the mismatch between U and Y. The loss therefore consisted only of the terminal defect and the control-mismatch penalty. While that two-term formulation proved effective in practice, no rigorous a posteriori error analysis was available, and the connection between the training objective and a provable error bound remained unclear. In the present paper, we build upon the decoupling idea and develop a comprehensive a posteriori error estimation framework that explicitly incorporates the control mismatch. Crucially, we augment the algorithm with an explicit pathwise residual term ℒRL_R that penalizes deviations from the backward dynamic equation at each time step. Together with the terminal defect ℒTL_T and the control mismatch ℒUL_U, this yields a three‑component loss that fully mirrors the structure of the a posteriori error bound we derive. The main contributions are as follows. First, we establish a continuous-time stability estimate for fully coupled FBSDEs under perturbations in the drift, diffusion, generator, terminal condition, and auxiliary control input. The estimate explicitly contains the mismatch term measuring the distance between the backward component and the auxiliary control. Second, we transfer this stability estimate to the discrete-time setting and obtain a computable a posteriori error bound for arbitrary adapted neural approximations. The bound depends only on the terminal defect, the pathwise residual, and the control mismatch, all of which can be evaluated from the numerical output. Therefore, the estimate is genuinely a posteriori in the sense that it does not require knowledge of the true solution. Third, we show through numerical experiments that the mismatch penalty is not merely a theoretical artifact. The first example is a high-dimensional linear–quadratic FBSDE for which an exact solution can be obtained through the associated Riccati equation. This allows us to compare the computable indicators with the true approximation error. The second example is a multidimensional Burgers-type FBSDE for which no closed-form reference solution is used. In this case, the a posteriori indicators serve as diagnostic tools for comparing different training strategies and for assessing the effect of removing individual loss components. The remainder of the paper is organized as follows. Section 2 introduces the mathematical setting of fully coupled FBSDEs and the neural network approximation framework. Section 3 presents the main continuous stability estimate and the discrete-time a posteriori error bounds. Section 4 contains the proofs of the main theorems. Section 5 reports the numerical experiments and discusses the behavior of the proposed indicators. Section 6 concludes the paper and outlines possible directions for future work. 2 Preliminaries This section introduces the basic concepts needed for this article, including the forward-backward stochastic differential equation (FBSDEs), deep neural networks and universal approximation. This preparatory discussion aims to establish the necessary theoretical foundation and computational framework for our subsequent algorithmic development. 2.1 Mathematical Framework of Fully-Coupled FBSDEs Let ℝnR^n denote the n-dimensional Euclidean space equipped with the standard inner product ⟨x,y⟩ x,y and the Euclidean norm ‖x‖2=⟨x,x⟩\|x\|^2= x,x for all x,y∈ℝnx,y ^n. Let ℝm×nR^m× n be the Hilbert space of all (m×n)(m× n)-matrices endowed with the inner product ⟨A,B⟩=tr(AB⊤),∀A,B∈ℝm×n, A,B =tr(AB ), ∀ A,B ^m× n, and the induced norm ‖A‖2=⟨A,A⟩\|A\|^2= A,A . Let T>0T>0 be fixed, and let (Ω,ℱ,,ℙ)( ,F,F,P) be a complete probability space equipped with a d-dimensional standard Brownian motion Bt0≤t≤T\B_t\_0≤ t≤ T. We denote by =ℱt0≤t≤TF=\F_t\_0≤ t≤ T the natural filtration generated by Bt\B_t\, augmented by all ℙP-null sets in ℱF. Let L2(ℱt;ℝn)L^2(F_t;R^n) be the space of ℱtF_t-measurable ℝnR^n-valued square-integrable random variables. We define the space M2(0,T;ℝn)M^2(0,T;R^n) as the set of all F-progressively measurable processes v:[0,T]×Ω→ℝnv:[0,T]× ^n satisfying [∫0T‖v(t)‖2t]<∞.E [ _0^T\|v(t)\|^2\,dt ]<∞. This space becomes a Hilbert space when endowed with the inner product ⟨u(⋅),v(⋅)⟩M2≔[∫0T⟨u(t),v(t)⟩ℝnt]. u(·),v(·) _M^2 [ _0^T u(t),v(t) _R^n\,dt ]. A fully coupled forward–backward stochastic differential equation defined on (Ω,ℱ,,ℙ)( ,F,F,P) has the following general form dXt=b(t,Xt,Yt,Zt)dt+σ(t,Xt,Yt,Zt)dBt,−dYt=f(t,Xt,Yt,Zt)dt−ZtdBt,X0=a,YT=g(XT), casesdX_t=b(t,X_t,Y_t,Z_t)\,dt+σ(t,X_t,Y_t,Z_t)\,dB_t,\\ -dY_t=f(t,X_t,Y_t,Z_t)\,dt-Z_t\,dB_t,\\ X_0=a, Y_T=g(X_T), cases (2.1) where X(⋅)∈M2(0,T;ℝn)X(·)∈ M^2(0,T;R^n), Y(⋅)∈M2(0,T;ℝm)Y(·)∈ M^2(0,T;R^m), and Z(⋅)∈M2(0,T;ℝm×d)Z(·)∈ M^2(0,T;R^m× d) are F-adapted stochastic processes, a and g(XT)g(X_T) are the initial and terminal conditions, respectively. b and σ denote the forward SDE drift coefficient and diffusion coefficient, and f denotes the generator of the backward SDE. When b(⋅)b(·) and σ(⋅)σ(·) do not depend on the processes Y(⋅),Z(⋅)Y(·),Z(·), the above fully coupled FBSDE can be simplified as follows dXt=b(t,Xt)dt+σ(t,Xt)dBt,−dYt=f(t,Xt,Yt,Zt)dt−ZtdBt,X0=a,YT=g(XT), casesdX_t=b(t,X_t)\,dt+σ(t,X_t)\,dB_t,\\ -dY_t=f(t,X_t,Y_t,Z_t)\,dt-Z_t\,dB_t,\\ X_0=a, Y_T=g(X_T), cases (2.2) which is a basic uncoupled BSDE. Definition 2.1. A triple of processes (X(⋅),Y(⋅),Z(⋅))∈M2(0,T;ℝn×ℝm×ℝm×d)(X(·),Y(·),Z(·))∈ M^2(0,T;R^n×R^m×R^m× d) is called an adapted solution of (2.1) if it satisfies the (2.1) ℙP-almost surely on [0,T][0,T]. Let u=(xyz)∈ℝn×ℝm×ℝm×d,A(t,u)=(−G⊤fGbGσ)(t,u),u= pmatrixx\\ y\\ z pmatrix ^n×R^m×R^m× d, A(t,u)= pmatrix-G f\\ Gb\\ Gσ pmatrix(t,u), where G is given a full rank matrix satisfied Gσ=(Gσ1,⋯,Gσd)Gσ=(G _1,·s,G _d) and ⟨u1,u2⟩=⟨x1,x2⟩+⟨y1,y2⟩+⟨z1,z2⟩. u^1,u^2 = x^1,x^2 + y^1,y^2 + z^1,z^2 . Assumption 1. 1. A(t,u)A(t,u) is uniformly Lipschitz with respect to u, A(⋅,u)A(·,u) is in M2(0,T)M^2(0,T) for each u. g(x)g(x) on x∈ℝnx ^n is uniformly Lipschitz with respect to x, and is in L2(ℱT;ℝn)L^2(F_T;R^n). The Lipschitz constant denotes as K>0K>0. 2. ϕ=b,σ,fφ=b,σ,f is uniformly Lipschitz with respect to u 3. ∀u,u¯∈ℝn×ℝm×ℝm×d∀\,u, u ^n×R^m×R^m× d and x,x¯x, x, there exist constants β1,β2,μ≥0 _1, _2,μ≥ 0 with β1+β2>0 _1+ _2>0 and μ+β2>0μ+ _2>0, such that: ⟨A(t,u)−A(t,u¯),u−u¯⟩ A(t,u)-A(t, u),u- u ≤−β1|G(x−x¯)|2−β2(|G⊤(y−y¯)|2+|G⊤(z−z¯)|2), ≤- _1|G(x- x)|^2- _2 (|G (y- y)|^2+|G (z- z)|^2 ), ⟨g(x)−g(x¯),G(x−x¯)⟩ g(x)-g( x),G(x- x) ≥μ|G(x−x¯)|2, ≥μ|G(x- x)|^2, where β1>0 _1>0 and μ>0μ>0 if m>nm>n, while β2>0 _2>0 if n>mn>m. Proposition 2.1. Under the 1, the above FBSDE (2.1) admits a unique process solution (X(⋅),Y(⋅),Z(⋅))(X(·),Y(·),Z(·)) in M2(0,T;ℝn×ℝm×ℝm×d)M^2(0,T;R^n×R^m×R^m× d). Proof Sketch. The existence follows from a Picard iteration scheme, while uniqueness is guaranteed by the Lipschitz continuity and monotonicity conditions in 1. The proof of the theorem is given in [HP95, PW99b]. ∎ A foundational connection between stochastic processes and deterministic partial differential equations is established via the nonlinear Feynman-Kac formula. This formula provides a probabilistic representation for solutions to a broad class of parabolic PDEs through the lens of (forward-)backward stochastic differential equations. For the more general fully coupled FBSDE system as defined in (2.1), where the drift and diffusion of the forward process depend explicitly on (Yt,Zt)(Y_t,Z_t), the associated PDE becomes quasilinear. The solution u(t,x)u(t,x), if sufficiently smooth, satisfies a quasilinear parabolic PDE whose precise form can be found in [MY99] or [PW99a]. 2.2 Feedforward Neural Networks and Universal Approximation Let din,dout∈ℕd_in,d_out , a feedforward neural network (FNN) is a hierarchical function mapping an input x∈ℝdinx ^d_in to an output fθ(x)∈ℝdoutf_θ(x) ^d_out through L layers of parameterized transformations. Mathematically, it is a function defined by h0(x) h_0(x) =x, =x, (2.3) hℓ(x) h_ (x) =ρ(Aℓhℓ−1+bℓ),for ℓ=1,…,L−1, =ρ(A_ h_ -1+b_ ), =1,…,L-1, fθ(x) f_θ(x) =hL(x)=ALhL−1+bL, =h_L(x)=A_Lh_L-1+b_L, where Aℓ∈ℝdℓ×dℓ−1A_ ^d_ × d_ -1 and bℓ∈ℝdℓb_ ^d_ are learnable parameters, ρ:ℝ→ℝρ:R is a component-wise activation function. The architecture is denoted by =(d0,d1,…,dL)S=(d_0,d_1,…,d_L), where d0=dind_0=d_in and dL=doutd_L=d_out are the input and output dimensions, respectively. The depth of the network is denoted by D()=LD(S)=L, the width is denoted by ‖∞=max1≤i≤Ldi\|S\|_∞= _1≤ i≤ Ld_i, and the total number of the parameters of network is denoted by ||=∑i=0L−1(di×di+1+di+1)|S|= _i=0^L-1(d_i× d_i+1+d_i+1). θ=Aℓ,bℓ0≤ℓ≤Lθ=\A_ ,b_ \_0≤ ≤ L are called the weight parameters, and Θ:=θ :=\θ\ is denoted as the space of parameters. Without loss of clarity, we also use (x;θ)U(x;θ) or (x;θ)V(x;θ) to represent the neural network in this following paper. In this paper, we employ a bounded, smooth activation function, such as the logistic function or the hyperbolic tangent, ρ(x)=11+e−x,orρ(x)=ex−e−xex+e−x,ρ(x)= 11+e^-x, ρ(x)= e^x-e^-xe^x+e^-x, as the activation function ρ for all hidden layers. The output layer uses the linear activation. A fundamental theoretical justification for employing neural networks to solve functional equations is the universal approximation property. Theorem 2.1 (Universal Approximation Theorem). Let K⊂ℝdinK ^d_in be a compact set and let f:K→ℝdoutf:K ^d_out be a continuous function. For any ε>0 >0, there exists a feedforward neural network (⋅;θ)U(·;θ) with architecture =(din,d1,dout)S=(d_in,d_1,d_out), i.e. a single hidden layer containing d1d_1 neurons, and a non-constant, bounded, continuous activation function ρ, such that supx∈K|f(x)−(x;θ)|<ε. _x∈ K |f(x)-U(x;θ) |< . The number of hidden neurons d1d_1 can be chosen finite but sufficiently large. Remark. In our deep learning algorithm we shall use deeper networks to parameterize the backward processes YtY_t, ZtZ_t and the auxiliary control UtU_t. Combined with the stability estimates of Section 3, this provides motivation that a sufficiently well-trained network can yield an accurate solution in the sense of the a posteriori error bounds. The quantitative relationship between the network architecture and the resulting approximation error is beyond the scope of the present paper and will be investigated separately. 3 The Main Results of Error Analysis In this section we establish the core error estimates of the paper. The presentation is organized from the general to the specific. We first prove in §3.1 a continuous-time stability result for fully coupled FBSDEs under the most general perturbations: the forward drift and diffusion may depend on an auxiliary control U that differs from the backward component Y. The classical stability bound for identical dynamics follows immediately as a corollary. In §3.2 we transfer this continuous theory to the discrete-time setting, deriving computable a posteriori error bounds for arbitrary adapted discrete approximations. Again, the general statement allows a mismatch between the backward variable and the control surrogate used in the forward coefficients; the standard case without mismatch is recovered as a special case. 3.1 Continuous Stability under General Perturbations Let (,,)(X,Y,Z) and U be adapted processes satisfying the perturbed dynamics dt=[b(t,t,t,t)+αt]dt+[σ(t,t,t,t)+βt]dBt,−dt=[f(t,t,t,t)+γt]dt−tdBt,0=a,T=g(T)−η, casesdX_t= [b(t,X_t,U_t,Z_t)+ _t ]\,dt+ [σ(t,X_t,U_t,Z_t)+ _t ]\,dB_t,\\ -dY_t= [f(t,X_t,Y_t,Z_t)+ _t ]\,dt-Z_t\,dB_t,\\ X_0=a, _T=g(X_T)-η, cases (3.1) where α,β,γα,β,γ are progressively measurable perturbations, η is a square-integrable terminal defect, and U is an auxiliary process that enters both the forward drift b and the forward diffusion σ. When ≡U , system (3.1) reduces to the standard perturbed FBSDE with identical dynamics; otherwise U models a control surrogate that may deviate from Y in both coefficients. Theorem 3.1 (General continuous stability). Suppose Assumption 1 holds. Let (X,Y,Z)(X,Y,Z) be the solution of (2.1) and (,,)(X,Y,Z) be the solution of (3.1). Then there exists a constant C>0C>0, depending only on the constants in Assumption 1 and T, such that sup0≤t≤T _0≤ t≤ T [|Xt−t|2+|Yt−t|2]+[∫0T|Zt−t|2t] [|X_t-X_t|^2+|Y_t-Y_t|^2 ]+E [ _0^T|Z_t-Z_t|^2\,dt ] (3.2) ≤C([|η|2]+[∫0T(|αt|2+|βt|2+|γt|2)t]+[∫0T|t−t|2t]). ≤ C (E [|η|^2 ]+E [ _0^T(| _t|^2+| _t|^2+| _t|^2)\,dt ]+E [ _0^T|Y_t-U_t|^2\,dt ] ). Corollary 3.1 (Stability under identical dynamics). If ≡U in (3.1), then the control mismatch integral vanishes and estimate (3.2) reduces to sup0≤t≤T _0≤ t≤ T [|Xt−t|2+|Yt−t|2]+[∫0T|Zt−t|2t] [|X_t-X_t|^2+|Y_t-Y_t|^2 ]+E [ _0^T|Z_t-Z_t|^2\,dt ] ≤C([|η|2]+[∫0T(|αt|2+|βt|2+|γt|2)t]). ≤ C (E [|η|^2 ]+E [ _0^T(| _t|^2+| _t|^2+| _t|^2)\,dt ] ). Proof of Corollary 3.1. Setting t=tU_t=Y_t for all t in (3.1) makes the control mismatch term [∫0T|t−t|2t]E [ _0^T|Y_t-U_t|^2\,dt ] identically zero. The stated bound then follows directly from Theorem 3.1. ∎ Remark. Theorem 3.1 shows that the solution map of the fully coupled FBSDE is stable in the mean-square sense with respect to simultaneous perturbations in drift, diffusion, generator, terminal condition, and—when present—forward control mismatch. The control mismatch penalty [∫0T|t−t|2t]E [ _0^T|Y_t-U_t|^2\,dt ] accounts for the discrepancy between the true backward component and the auxiliary control used in both b and σ. 3.2 A Posteriori Error Estimates for Discrete Approximations We now translate the continuous stability theory into computable error bounds for discrete-time schemes. Let π=0=t0<⋯<tN=Tπ=\0=t_0<…<t_N=T\ be a partition of [0,T][0,T] with Δti:=ti+1−ti t_i:=t_i+1-t_i and |π|:=maxiΔti|π|:= _i t_i. Given an ℱtiF_t_i-adapted discrete trajectory (X^i,Y^i,Z^i,U^i)i=0N\( X_i, Y_i, Z_i, U_i)\_i=0^N with X^0=a X_0=a, we denote by (X¯,Y¯,Z¯,U¯)( X, Y, Z, U) its piecewise constant interpolation. Define the local residuals RiX R_i^X :=X^i+1−X^i−b(ti,X^i,U^i,Z^i)Δti−σ(ti,X^i,U^i,Z^i)ΔBi, := X_i+1- X_i-b(t_i, X_i, U_i, Z_i)\, t_i-σ(t_i, X_i, U_i, Z_i)\, B_i, (3.3) RiY R_i^Y :=Y^i+1−Y^i+f(ti,X^i,Y^i,Z^i)Δti−Z^iΔBi, := Y_i+1- Y_i+f(t_i, X_i, Y_i, Z_i)\, t_i- Z_i\, B_i, (3.4) the terminal mismatch η:=g(X^N)−Y^Nη:=g( X_N)- Y_N, and the computable error indicators ℜπ:=∑i=0N−1[|RiX|2Δti+|RiY|2Δti],π:=∑i=0N−1[|Y^i−U^i|2]Δti. R_π:= _i=0^N-1E [ |R_i^X|^2 t_i+ |R_i^Y|^2 t_i ], P_π:= _i=0^N-1E [| Y_i- U_i|^2 ]\, t_i. (3.5) To convert the piecewise constant approximation into a continuous-time perturbed system, we require the following mild time-regularity of the coefficients. Assumption 2 (1/21/2-Hölder continuity in time). There exists LH>0L_H>0 such that for all s,t∈[0,T]s,t∈[0,T] and all (x,y,z)(x,y,z), |b(t,x,y,z)−b(s,x,y,z)|+|σ(t,x,y,z)−σ(s,x,y,z)|+|f(t,x,y,z)−f(s,x,y,z)|≤LH|t−s|1/2.|b(t,x,y,z)-b(s,x,y,z)|+|σ(t,x,y,z)-σ(s,x,y,z)|+|f(t,x,y,z)-f(s,x,y,z)|≤ L_H|t-s|^1/2. Assumption 3. The discrete approximation (X^i,Y^i,Z^i,U^i)i=0N\( X_i, Y_i, Z_i, U_i)\_i=0^N satisfies the uniform moment bound sup0≤i≤N[|X^i|2+|Y^i|2+|Z^i|2+|U^i|2]≤M, _0≤ i≤ NE [| X_i|^2+| Y_i|^2+| Z_i|^2+| U_i|^2 ]≤ M, where M>0M>0 is a constant that may depend on the network architecture and the training outcome, but is independent of the time step Δt t. Theorem 3.2 (General a posteriori error estimate). Suppose Assumptions 1– 3 hold, and let (X,Y,Z)(X,Y,Z) be the solution of (2.1). Then there exists C>0C>0, depending only on the constants in Assumption 1, LHL_H, M and T, such that sup0≤t≤T[|Xt−X¯t|2+|Yt−Y¯t|2]+[∫0T|Zt−Z¯t|2t]≤C([|η|2]+ℜπ+π+|π|). _0≤ t≤ TE [|X_t- X_t|^2+|Y_t- Y_t|^2 ]+E [ _0^T|Z_t- Z_t|^2\,dt ]≤ C (E [|η|^2 ]+ R_π+ P_π+|π| ). (3.6) In particular, on the grid points, max0≤i≤N[|Xti−X^i|2+|Yti−Y^i|2]+[∫0T|Zt−Z¯t|2t]≤C([|η|2]+ℜπ+π+|π|). _0≤ i≤ NE [|X_t_i- X_i|^2+|Y_t_i- Y_i|^2 ]+E [ _0^T|Z_t- Z_t|^2\,dt ]≤ C (E [|η|^2 ]+ R_π+ P_π+|π| ). (3.7) Corollary 3.2 (A posteriori estimate without control mismatch). If U^i≡Y^i U_i≡ Y_i for all i, then π=0 P_π=0 and the bound (3.6) simplifies to sup0≤t≤T[|Xt−X¯t|2+|Yt−Y¯t|2]+[∫0T|Zt−Z¯t|2t]≤C([|η|2]+ℜπ+|π|). _0≤ t≤ TE [|X_t- X_t|^2+|Y_t- Y_t|^2 ]+E [ _0^T|Z_t- Z_t|^2\,dt ]≤ C (E [|η|^2 ]+ R_π+|π| ). Proof of Corollary 3.2. If U^i=Y^i U_i= Y_i for every i, then by definition π=∑i=0N−1[|Y^i−U^i|2]Δti=0 P_π= _i=0^N-1E [| Y_i- U_i|^2 ] t_i=0. Substituting this into (3.6) yields the claimed estimate. ∎ Corollary 3.3 (Euler-type schemes). Assume the discrete trajectory is generated by the forward Euler updates X^i+1 X_i+1 =X^i+b(ti,X^i,U^i,Z^i)Δti+σ(ti,X^i,U^i,Z^i)ΔBi, = X_i+b(t_i, X_i, U_i, Z_i)\, t_i+σ(t_i, X_i, U_i, Z_i)\, B_i, (3.8) Y^i+1 Y_i+1 =Y^i−f(ti,X^i,Y^i,Z^i)Δti+Z^iΔBi+εiY, = Y_i-f(t_i, X_i, Y_i, Z_i)\, t_i+ Z_i\, B_i+ _i^Y, (3.9) where εiY _i^Y is an ℱti+1F_t_i+1-measurable local defect. Then ℜπ=∑i=0N−1[|εiY|2]/Δti R_π= _i=0^N-1E [| _i^Y|^2 ]/ t_i, and (3.6) becomes sup0≤t≤T[|Xt−X¯t|2+|Yt−Y¯t|2] _0≤ t≤ TE [|X_t- X_t|^2+|Y_t- Y_t|^2 ] +[∫0T|Zt−Z¯t|2t] +E [ _0^T|Z_t- Z_t|^2\,dt ] ≤C([|η|2]+∑i=0N−1[|εiY|2]Δti+π+|π|). ≤ C (E [|η|^2 ]+ _i=0^N-1 E [| _i^Y|^2 ] t_i+ P_π+|π| ). If in addition εiY=0 _i^Y=0 and U^i=Y^i U_i= Y_i, we recover the classical O(|π|)O(|π|) consistency bound. Proof of Corollary 3.3. Under the forward Euler update (3.8), the forward residual satisfies RiX=0R_i^X=0 for all i. From (3.9), we obtain RiY=εiYR_i^Y= _i^Y. Consequently, ℜπ=∑i=0N−1[|RiX|2Δti+|RiY|2Δti]=∑i=0N−1[|εiY|2]Δti. R_π= _i=0^N-1E [ |R_i^X|^2 t_i+ |R_i^Y|^2 t_i ]= _i=0^N-1 E [| _i^Y|^2 ] t_i. Inserting this expression for ℜπ R_π into Theorem 3.2 yields the first bound. If additionally εiY=0 _i^Y=0 and U^i=Y^i U_i= Y_i, then ℜπ=0 R_π=0 and π=0 P_π=0, giving the O(|π|)O(|π|) estimate. ∎ Remark. The indicators η, ℜπ R_π, and π P_π are computable directly from the discrete output without knowledge of the true solution. They serve as practical a posteriori criteria for model selection, adaptive refinement, and stopping in deep learning algorithms for fully coupled FBSDEs. The mismatch penalty π P_π is particularly relevant when separate network parametrizations are used for the backward component Y Y and the control surrogate U U that enters both b and σ. If one further sets ℜπ=0 R_π=0 (i.e., discards the pathwise residual penalty) while keeping U U decoupled from Y Y, the scheme reduces to the second case of [JPP+20]. In that case, the general estimate (3.6) still formally applies, but the residual indicator ℜπ R_π may become uninformative because the residual is not explicitly controlled during training. This explains the lack of rigorous error guarantees for that earlier algorithm. Remark. The constant C in Theorems 3.1 and 3.2 originates solely from the stability analysis of the FBSDE and the time discretization. It therefore controls the error under the assumption that the discrete approximation (X^i,Y^i,Z^i,U^i)\( X_i, Y_i, Z_i, U_i)\ is already available. When such an approximation is generated by a deep learning algorithm, additional sources of error—primarily the approximation error due to finite network capacity and the statistical error from finite-sample Monte Carlo estimates—are not accounted for in C. Consequently, the a posteriori bound should be interpreted as a computable error indicator that faithfully reflects the quality of the discrete trajectory, rather than a certified numerical upper bound in the presence of those unmodeled errors. Its practical value lies in guiding the training process and comparing different approximations without reference to the true solution, which is exactly how we employ it in Section 5. 4 Proofs of the Main Theorems In this section we provide the complete and detailed proofs of the two principal theorems stated in Section 3: the general continuous stability estimate (Theorem 3.1) and the general a posteriori error bound for discrete approximations (Theorem 3.2). 4.1 Proof of Theorem 3.1 Proof. Let ΔXt:=Xt−t,ΔYt:=Yt−t,ΔZt:=Zt−t. X_t:=X_t-X_t, Y_t:=Y_t-Y_t, Z_t:=Z_t-Z_t. Also define ΔUt:=Yt−t. U_t:=Y_t-U_t. Then, by (2.1) and (3.1), dΔXt=(b(t,Xt,Yt,Zt)−b(t,t,t,t)−αt)dt+(σ(t,Xt,Yt,Zt)−σ(t,t,t,t)−βt)dBt,−dΔYt=(f(t,Xt,Yt,Zt)−f(t,t,t,t)−γt)dt−ΔZtdBt,ΔX0=0,ΔYT=g(XT)−g(T)+η. casesd X_t= (b(t,X_t,Y_t,Z_t)-b(t,X_t,U_t,Z_t)- _t )\,dt\\ + (σ(t,X_t,Y_t,Z_t)-σ(t,X_t,U_t,Z_t)- _t )\,dB_t,\\ -d Y_t= (f(t,X_t,Y_t,Z_t)-f(t,X_t,Y_t,Z_t)- _t )\,dt- Z_t\,dB_t,\\ X_0=0, Y_T=g(X_T)-g(X_T)+η. cases (4.1) For convenience, write Δbt b_t :=b(t,Xt,Yt,Zt)−b(t,t,t,t), :=b(t,X_t,Y_t,Z_t)-b(t,X_t,Y_t,Z_t), δbt δ b_t :=b(t,t,t,t)−b(t,t,t,t), :=b(t,X_t,Y_t,Z_t)-b(t,X_t,U_t,Z_t), Δσt _t :=σ(t,Xt,Yt,Zt)−σ(t,t,t,t), :=σ(t,X_t,Y_t,Z_t)-σ(t,X_t,Y_t,Z_t), δσt δ _t :=σ(t,t,t,t)−σ(t,t,t,t), :=σ(t,X_t,Y_t,Z_t)-σ(t,X_t,U_t,Z_t), Δft f_t :=f(t,Xt,Yt,Zt)−f(t,t,t,t). :=f(t,X_t,Y_t,Z_t)-f(t,X_t,Y_t,Z_t). Then (4.1) becomes dΔXt=(Δbt+δbt−αt)dt+(Δσt+δσt−βt)dBt,dΔYt=−(Δft−γt)dt+ΔZtdBt,ΔX0=0,ΔYT=g(XT)−g(T)+η. casesd X_t= ( b_t+δ b_t- _t )\,dt+ ( _t+δ _t- _t )\,dB_t,\\ d Y_t=- ( f_t- _t )\,dt+ Z_t\,dB_t,\\ X_0=0, Y_T=g(X_T)-g(X_T)+η. cases (4.2) Forward estimate for ΔX X. Applying Itô’s formula to |ΔXt|2| X_t|^2 and using (4.2), d|ΔXt|2=2⟨ΔXt,Δbt+δbt−αt⟩dt+2⟨ΔXt,Δσt+δσt−βt⟩dBt+|Δσt+δσt−βt|2dt.d| X_t|^2=2 X_t, b_t+δ b_t- _t dt+2 X_t, _t+δ _t- _t dB_t+| _t+δ _t- _t|^2dt. Taking expectations and integrating from 0 to t, |ΔXt|2=∫0t(2⟨ΔXs,Δbs+δbs−αs⟩+|Δσs+δσs−βs|2)s.E| X_t|^2=E _0^t (2 X_s, b_s+δ b_s- _s +| _s+δ _s- _s|^2 )ds. Using the Lipschitz properties of b and σ, we obtain the pointwise bounds 2⟨ΔX,Δb+δb−α⟩ 2 X, b+δ b-α ≤C1|ΔX|2+|ΔY|2+|ΔZ|2+|α|2+|−|2, ≤ C_1| X|^2+| Y|^2+| Z|^2+|α|^2+|Y-U|^2, |Δσ+δσ−β|2 | σ+δσ-β|^2 ≤C2(|ΔX|2+|ΔY|2+|ΔZ|2+|β|2+|−|2). ≤ C_2 (| X|^2+| Y|^2+| Z|^2+|β|^2+|Y-U|^2 ). Thus |ΔXt|2≤C∫0t|ΔXs|2s+C∫0T(|ΔYs|2+|ΔZs|2)s+C∫0T(|αs|2+|βs|2+|s−s|2)s.E| X_t|^2≤ C _0^tE| X_s|^2ds+C _0^TE (| Y_s|^2+| Z_s|^2 )ds+C _0^TE (| _s|^2+| _s|^2+|Y_s-U_s|^2 )ds. Gronwall’s inequality then yields sup0≤t≤T|ΔXt|2≤C(∫0T(|ΔYt|2+|ΔZt|2)t+∫0T(|αt|2+|βt|2+|t−t|2)t). _0≤ t≤ TE| X_t|^2≤ C (E _0^T(| Y_t|^2+| Z_t|^2)dt+E _0^T(| _t|^2+| _t|^2+|Y_t-U_t|^2)dt ). (4.3) Backward estimate for (ΔY,ΔZ)( Y, Z). Apply Itô’s formula to |ΔYt|2| Y_t|^2: d|ΔYt|2=−2⟨ΔYt,Δft−γt⟩dt+2⟨ΔYt,ΔZt⟩dBt+|ΔZt|2dt.d| Y_t|^2=-2 Y_t, f_t- _t dt+2 Y_t, Z_t dB_t+| Z_t|^2dt. Integrating from t to T and taking expectation, |ΔYt|2+∫tT|ΔZs|2s=|ΔYT|2+2∫tT⟨ΔYs,Δfs−γs⟩s.E| Y_t|^2+E _t^T| Z_s|^2ds=E| Y_T|^2+2E _t^T Y_s, f_s- _s ds. The Lipschitz continuity of f gives |Δfs|≤K(|ΔXs|+|ΔYs|+|ΔZs|)| f_s|≤ K(| X_s|+| Y_s|+| Z_s|). Then 2⟨ΔYs,Δfs−γs⟩≤K2|ΔXs|2+(1+2K+2K2)|ΔYs|2+12|ΔZs|2+|γs|2.2 Y_s, f_s- _s ≤ K^2| X_s|^2+(1+2K+2K^2)| Y_s|^2+ 12| Z_s|^2+| _s|^2. Hence |ΔYt|2+∫tT|ΔZs|2s≤|ΔYT|2+C∫tT(|ΔXs|2+|ΔYs|2+|γs|2)s+12∫tT|ΔZs|2s.E| Y_t|^2+E _t^T| Z_s|^2ds | Y_T|^2+CE _t^T (| X_s|^2+| Y_s|^2+| _s|^2 )ds+ 12E _t^T| Z_s|^2ds. Moving the 12∫|ΔZ|2 12 | Z|^2 term to the left, |ΔYt|2+12∫tT|ΔZs|2s | Y_t|^2+ 12E _t^T| Z_s|^2ds ≤|ΔYT|2+C∫tT(|ΔXs|2+|ΔYs|2+|γs|2)s | Y_T|^2+CE _t^T (| X_s|^2+| Y_s|^2+| _s|^2 )ds ≤|ΔYT|2+C∫0T(|ΔXt|2+|γt|2)t+C∫tT|ΔYs|2s. | Y_T|^2+CE _0^T (| X_t|^2+| _t|^2 )dt+CE _t^T| Y_s|^2ds. Since g is Lipschitz, |ΔYT|2≤C(|ΔXT|2+|η|2)E| Y_T|^2≤ C(E| X_T|^2+E|η|^2). Applying Gronwall’s inequality to the function ϕ(t)=|ΔYt|2φ(t)=E| Y_t|^2 (noting that the integral of |ΔY|2| Y|^2 can be absorbed) yields sup0≤t≤T|ΔYt|2≤C(∫0T(|ΔXt|2+|γt|2)t+|ΔXT|2+|η|2). _0≤ t≤ TE| Y_t|^2≤ C (E _0^T(| X_t|^2+| _t|^2)dt+E| X_T|^2+E|η|^2 ). (4.4) Together with the bound for ∫|ΔZ|2 | Z|^2 obtained from the same inequality, we get sup0≤t≤T|ΔYt|2+∫0T|ΔZt|2t≤C(∫0T(|ΔXt|2+|γt|2)t+|ΔXT|2+|η|2). _0≤ t≤ TE| Y_t|^2+E _0^T| Z_t|^2dt≤ C (E _0^T(| X_t|^2+| _t|^2)dt+E| X_T|^2+E|η|^2 ). (4.5) For u=(x,y,z)u=(x,y,z), recall A(t,u):=(−G⊤f(t,x,y,z)Gb(t,x,y,z)Gσ(t,x,y,z)).A(t,u):= pmatrix-G f(t,x,y,z)\\ Gb(t,x,y,z)\\ Gσ(t,x,y,z) pmatrix. Applying Itô’s formula to ⟨GΔXt,ΔYt⟩ G X_t, Y_t gives d⟨GΔXt,ΔYt⟩ d G X_t, Y_t =⟨G(Δbt+δbt−αt),ΔYt⟩dt+⟨G(Δσt+δσt−βt),ΔYt⟩dBt = G( b_t+δ b_t- _t), Y_t dt+ G( _t+δ _t- _t), Y_t dB_t +⟨GΔXt,−(Δft−γt)⟩dt+⟨GΔXt,ΔZt⟩dBt + G X_t,-( f_t- _t) dt+ G X_t, Z_t dB_t +⟨G(Δσt+δσt−βt),ΔZt⟩dt. + G( _t+δ _t- _t), Z_t dt. Integrating from 0 to T and taking expectation (the stochastic integrals are martingales and vanish), we obtain [⟨GΔXT,ΔYT⟩] [ G X_T, Y_T ] =[∫0T⟨A(t,ut)−A(t,u¯t),ut−u¯t⟩t] =E [ _0^T A(t,u_t)-A(t, u_t),u_t- u_t dt ] (4.6) +[∫0T⟨Gδbt,ΔYt⟩t]+[∫0T⟨Gδσt,ΔZt⟩t] +E [ _0^T Gδ b_t, Y_t dt ]+E [ _0^T Gδ _t, Z_t dt ] +[∫0T(⟨GΔXt,γt⟩−⟨Gαt,ΔYt⟩−⟨Gβt,ΔZt⟩)t], +E [ _0^T ( G X_t, _t - G _t, Y_t - G _t, Z_t )dt ], where ut:=(Xt,Yt,Zt),u¯t:=(t,t,t),u_t:=(X_t,Y_t,Z_t), u_t:=(X_t,Y_t,Z_t), and ⟨A(t,ut)−A(t,u¯t),ut−u¯t⟩=−⟨GΔXt,Δft⟩+⟨GΔbt,ΔYt⟩+⟨GΔσt,ΔZt⟩. A(t,u_t)-A(t, u_t),u_t- u_t =- G X_t, f_t + G b_t, Y_t + G _t, Z_t . By the Lipschitz continuity of b and σ (Assumption 1), |δbt|=|b(t,t,t,t)−b(t,t,t,t)|≤K|t−t|,|δ b_t|=|b(t,X_t,Y_t,Z_t)-b(t,X_t,U_t,Z_t)|≤ K|Y_t-U_t|, (4.7) |δσt|=|σ(t,t,t,t)−σ(t,t,t,t)|≤K|t−t|.|δ _t|=|σ(t,X_t,Y_t,Z_t)-σ(t,X_t,U_t,Z_t)|≤ K|Y_t-U_t|. (4.8) Moreover, by the monotonicity condition in Assumption 1, ⟨A(t,ut)−A(t,u¯t),ut−u¯t⟩≤−β1|GΔXt|2−β2(|G⊤ΔYt|2+|G⊤ΔZt|2). A(t,u_t)-A(t, u_t),u_t- u_t ≤- _1|G X_t|^2- _2 (|G Y_t|^2+|G Z_t|^2 ). (4.9) We will repeatedly use Young’s inequality 2ab≤εa2+ε−1b22ab≤ a^2+ ^-1b^2 with different small parameters ε1 _1 (for terms involving |GΔX|2|G X|^2) and ε2 _2 (for |ΔY|2,|ΔZ|2| Y|^2,| Z|^2). The precise values of ε1,ε2 _1, _2 will be chosen later sufficiently small to absorb certain terms. Constants CεC_ depend on ε and may change from line to line. From (4.7) and (4.8) we obtain |⟨Gδbt,ΔYt⟩| | Gδ b_t, Y_t | ≤ε2|ΔYt|2+Cε2|t−t|2, ≤ _2| Y_t|^2+C_ _2|Y_t-U_t|^2, |⟨Gδσt,ΔZt⟩| | Gδ _t, Z_t | ≤ε2|ΔZt|2+Cε2|t−t|2. ≤ _2| Z_t|^2+C_ _2|Y_t-U_t|^2. Hence |⟨Gδbt,ΔYt⟩|+|⟨Gδσt,ΔZt⟩|≤ε2(|ΔYt|2+|ΔZt|2)+Cε2|t−t|2.| Gδ b_t, Y_t |+| Gδ _t, Z_t |≤ _2(| Y_t|^2+| Z_t|^2)+C_ _2|Y_t-U_t|^2. (4.10) Similarly, |⟨Gαt,ΔYt⟩| | G _t, Y_t | ≤ε2|ΔYt|2+Cε2|αt|2, ≤ _2| Y_t|^2+C_ _2| _t|^2, |⟨Gβt,ΔZt⟩| | G _t, Z_t | ≤ε2|ΔZt|2+Cε2|βt|2, ≤ _2| Z_t|^2+C_ _2| _t|^2, |⟨GΔXt,γt⟩| | G X_t, _t | ≤ε1|GΔXt|2+Cε1|γt|2. ≤ _1|G X_t|^2+C_ _1| _t|^2. Thus |⟨Gαt,ΔYt⟩−⟨GΔXt,γt⟩+⟨Gβt,ΔZt⟩| | G _t, Y_t - G X_t, _t + G _t, Z_t | (4.11) ≤ε1|GΔXt|2+ε2(|ΔYt|2+|ΔZt|2)+Cε1|γt|2+Cε2(|αt|2+|βt|2). ≤ _1|G X_t|^2+ _2(| Y_t|^2+| Z_t|^2)+C_ _1| _t|^2+C_ _2(| _t|^2+| _t|^2). Inserting the monotonicity bound (4.9) and the estimates (4.10), (4.11) into (4.6) yields ⟨GΔXT,ΔYT⟩+β1∫0T|GΔXt|2t+β2∫0T(|G⊤ΔYt|2+|G⊤ΔZt|2)t G X_T, Y_T + _1E _0^T|G X_t|^2dt+ _2E _0^T (|G Y_t|^2+|G Z_t|^2 )dt (4.12) ≤∫0T(ε1|GΔXt|2+2ε2(|ΔYt|2+|ΔZt|2))t _0^T ( _1|G X_t|^2+2 _2(| Y_t|^2+| Z_t|^2) )dt +∫0T(Cε2(|αt|2+|βt|2+|t−t|2)+Cε1|γt|2)t. +E _0^T (C_ _2(| _t|^2+| _t|^2+|Y_t-U_t|^2)+C_ _1| _t|^2 )dt. By Assumption 1, the terminal condition satisfies ⟨GΔXT,g(XT)−g(T)⟩≥μ|GΔXT|2. G X_T,g(X_T)-g(X_T) ≥μ|G X_T|^2. Using ΔYT=g(XT)−g(T)+η Y_T=g(X_T)-g(X_T)+η, we have ⟨GΔXT,ΔYT⟩=⟨GΔXT,g(XT)−g(T)⟩+⟨GΔXT,η⟩≥μ|GΔXT|2−|⟨GΔXT,η⟩|. G X_T, Y_T = G X_T,g(X_T)-g(X_T) + G X_T,η ≥μ|G X_T|^2-| G X_T,η |. Substituting this into (4.12) and rearranging, we obtain μ|GΔXT|2+β1∫0T|GΔXt|2t+β2∫0T(|G⊤ΔYt|2+|G⊤ΔZt|2)t |G X_T|^2+ _1E _0^T|G X_t|^2dt+ _2E _0^T (|G Y_t|^2+|G Z_t|^2 )dt (4.13) ≤|⟨GΔXT,η⟩|+∫0T(ε1|GΔXt|2+2ε2(|ΔYt|2+|ΔZt|2))t | G X_T,η |+E _0^T ( _1|G X_t|^2+2 _2(| Y_t|^2+| Z_t|^2) )dt +∫0T(Cε2(|αt|2+|βt|2+|t−t|2)+Cε1|γt|2)t. +E _0^T (C_ _2(| _t|^2+| _t|^2+|Y_t-U_t|^2)+C_ _1| _t|^2 )dt. Case 1: m>nm>n (the X-coercive regime). Here G has full column rank, so there exists cG>0c_G>0 such that |Gx|≥cG|x||Gx|≥ c_G|x| for all x∈ℝnx ^n. In particular, |GΔX|2≥cG2|ΔX|2|G X|^2≥ c_G^2| X|^2 and |GΔXT|2≥cG2|ΔXT|2|G X_T|^2≥ c_G^2| X_T|^2. Young’s inequality gives |⟨GΔXT,η⟩|≤μ2|GΔXT|2+12μ|η|2.| G X_T,η |≤ μ2|G X_T|^2+ 12μ|η|^2. Taking expectations, substituting it into (4.12) and rearranging terms, we obtain μ2|GΔXT|2+β1∫0T|GΔXt|2t+β2∫0T(|G⊤ΔYt|2+|G⊤ΔZt|2)t μ2E|G X_T|^2+ _1E _0^T|G X_t|^2dt+ _2E _0^T (|G Y_t|^2+|G Z_t|^2 )dt (4.14) ≤12μ|η|2+∫0T(ε1|GΔXt|2+2ε2(|ΔYt|2+|ΔZt|2))t ≤ 12μE|η|^2+E _0^T ( _1|G X_t|^2+2 _2(| Y_t|^2+| Z_t|^2) )dt +∫0T(Cε2(|αt|2+|βt|2+|t−t|2)+Cε1|γt|2)t. +E _0^T (C_ _2(| _t|^2+| _t|^2+|Y_t-U_t|^2)+C_ _1| _t|^2 )dt. Choose ε1=β1/2 _1= _1/2 in (4.14). Then the terms ε1∫|GΔX|2 _1 |G X|^2 on the right can be absorbed into the left-hand side β1∫|GΔX|2 _1 |G X|^2. Dropping the non‑negative β2 _2 term, we obtain μ2|GΔXT|2+β12∫0T|GΔXt|2t μ2E|G X_T|^2+ _12E _0^T|G X_t|^2dt ≤Cε1|η|2+2ε2∫0T(|ΔYt|2+|ΔZt|2)t ≤ C_ _1E|η|^2+2 _2E _0^T(| Y_t|^2+| Z_t|^2)dt +Cε2∫0T(|αt|2+|βt|2+|t−t|2)t+Cε1∫0T|γt|2t. +C_ _2E _0^T(| _t|^2+| _t|^2+|Y_t-U_t|^2)dt+C_ _1E _0^T| _t|^2dt. The left-hand side dominates μcG22|ΔXT|2+β1cG22∫0T|ΔXt|2t μ c_G^22E| X_T|^2+ _1c_G^22E _0^T| X_t|^2dt from the full column rank property. Consequently, ∫0T|ΔXt|2t+|ΔXT|2 _0^T| X_t|^2dt+E| X_T|^2 ≤C(|η|2+ε2∫0T(|ΔYt|2+|ΔZt|2)dt ≤ C (E|η|^2+ _2E _0^T(| Y_t|^2+| Z_t|^2)dt +∫0T(|αt|2+|βt|2+|γt|2+|t−t|2)dt). +E _0^T(| _t|^2+| _t|^2+| _t|^2+|Y_t-U_t|^2)dt ). Note that (4.5) holds, we have sup0≤t≤T[|ΔYt|2]+[∫0T|ΔZt|2t] _0≤ t≤ TE [| Y_t|^2 ]+E [ _0^T| Z_t|^2dt ] ≤ ≤ C([∫0T|ΔXt|2t]+[|ΔXT|2]+[∫0T|γt|2t]+[|η|2]) C (E [ _0^T| X_t|^2dt ]+E [| X_T|^2 ]+E [ _0^T| _t|^2dt ]+E [|η|^2 ] ) ≤ ≤ C([∫0Tε2(|ΔYt|2+|ΔZt|2)dt]+[|η|2]+[∫0T(|αt|2+|βt|2+|γt|2)dt] C (E [ _0^T _2 (| Y_t|^2+| Z_t|^2 )dt ]+E [|η|^2 ]+E [ _0^T (| _t|^2+| _t|^2+| _t|^2 )dt ] +[∫0T|t−t|2dt]). +E [ _0^T|Y_t-U_t|^2dt ] ). Choosing a small ε2 _2 (e.g., ε2=12C(1+T) _2= 12C(1+T)), we get the estimate of Y,ZY,Z sup0≤t≤T[|ΔYt|2]+[∫0T|ΔZt|2t] _0≤ t≤ TE [| Y_t|^2 ]+E [ _0^T| Z_t|^2dt ] ≤ ≤ C([|η|2]+[∫0T(|αt|2+|βt|2+|γt|2)t]+[∫0T|t−t|2t]). C (E [|η|^2 ]+E [ _0^T (| _t|^2+| _t|^2+| _t|^2 )dt ]+E [ _0^T|Y_t-U_t|^2dt ] ). Substituting it into (4.3) can derive the estimation of X sup0≤t≤T[|ΔXt|2]≤C([|η|2]+[∫0T(|αt|2+|βt|2+|γt|2)t]+[∫0T|t−t|2t]). _0≤ t≤ TE [| X_t|^2 ]≤ C (E [|η|^2 ]+E [ _0^T (| _t|^2+| _t|^2+| _t|^2 )dt ]+E [ _0^T|Y_t-U_t|^2dt ] ). Thus the desired estimate holds in Case 1. Case 2: m<nm<n (the (Y,Z)(Y,Z)-coercive regime). In this case G has full row rank, so there exists cG′>0c_G >0 such that |G⊤y|≥cG′|y||G y|≥ c_G |y| for all y∈ℝmy ^m, and the same holds columnwise for matrices. Hence |G⊤ΔY|2≥(cG′)2|ΔY|2|G Y|^2≥(c_G )^2| Y|^2 and |G⊤ΔZ|2≥(cG′)2|ΔZ|2|G Z|^2≥(c_G )^2| Z|^2. Similarly, Young’s inequality gives |⟨GΔXT,η⟩|≤ε1|GΔXT|2+14ε1|η|2.| G X_T,η |≤ _1|G X_T|^2+ 14 _1|η|^2. Thus, we will obtain from (4.12) μ|GΔXT|2+β1∫0T|GΔXt|2t+β2∫0T(|G⊤ΔYt|2+|G⊤ΔZt|2)t |G X_T|^2+ _1E _0^T|G X_t|^2dt+ _2E _0^T (|G Y_t|^2+|G Z_t|^2 )dt (4.15) ≤Cε1|η|2+ε1|GΔXT|2+∫0T(ε1|GΔXt|2+2ε2(|ΔYt|2+|ΔZt|2))t ≤ C_ _1E|η|^2+ _1E|G X_T|^2+E _0^T ( _1|G X_t|^2+2 _2(| Y_t|^2+| Z_t|^2) )dt +∫0T(Cε2(|αt|2+|βt|2+|t−t|2)+Cε1|γt|2)t. +E _0^T (C_ _2(| _t|^2+| _t|^2+|Y_t-U_t|^2)+C_ _1| _t|^2 )dt. Choose ε2=β2(cG′)2/4 _2= _2(c_G )^2/4 in (4.15). Dropping the non‑negative μ|GΔXT|2 |G X_T|^2 and β1∫|GΔX|2 _1 |G X|^2 and rearranging terms, we obtain β2(cG′)22∫0T(|ΔYt|2+|ΔZt|2)t _2(c_G )^22E _0^T(| Y_t|^2+| Z_t|^2)dt ≤Cε1|η|2+ε1|GΔXT|2+ε1∫0T|GΔXt|2t ≤ C_ _1E|η|^2+ _1E|G X_T|^2+ _1E _0^T|G X_t|^2dt +Cε2∫0T(|αt|2+|βt|2+|t−t|2)t+Cε1∫0T|γt|2t. +C_ _2E _0^T(| _t|^2+| _t|^2+|Y_t-U_t|^2)dt+C_ _1E _0^T| _t|^2dt. Now use the forward estimate (4.3) (note |GΔX|2≤‖G‖2|ΔX|2|G X|^2≤\|G\|^2| X|^2), we have |GΔXT|2+∫0T|GΔXt|2t≤C(∫0T(|ΔYt|2+|ΔZt|2)t+∫0T(|αt|2+|βt|2+|t−t|2)t).E|G X_T|^2+E _0^T|G X_t|^2dt≤ C (E _0^T(| Y_t|^2+| Z_t|^2)dt+E _0^T(| _t|^2+| _t|^2+|Y_t-U_t|^2)dt ). Substituting this into the previous inequality gives ∫0T(|ΔYt|2+|ΔZt|2)t≤C(|η|2+ε1∫0T(|ΔY|2+|ΔZ|2)+∫0T(|α|2+|β|2+|γ|2+|−|2)).E _0^T(| Y_t|^2+| Z_t|^2)dt≤ C (E|η|^2+ _1E _0^T(| Y|^2+| Z|^2)+E _0^T(|α|^2+|β|^2+|γ|^2+|Y-U|^2) ). Choosing ε1 _1 sufficiently small (e.g., ε1=12C _1= 12C) allows us to absorb the ε1 _1 term, yielding ∫0T(|ΔYt|2+|ΔZt|2)t≤C(|η|2+∫0T(|αt|2+|βt|2+|γt|2+|t−t|2)t).E _0^T(| Y_t|^2+| Z_t|^2)dt≤ C (E|η|^2+E _0^T(| _t|^2+| _t|^2+| _t|^2+|Y_t-U_t|^2)dt ). (4.16) Plugging this bound into (4.3) gives the same estimate for |ΔXt|2| X_t|^2. sup0≤t≤T|ΔXt|2≤C(|η|2+∫0T(|αt|2+|βt|2+|γt|2+|t−t|2)t). _0≤ t≤ TE| X_t|^2≤ C (E|η|^2+E _0^T(| _t|^2+| _t|^2+| _t|^2+|Y_t-U_t|^2)dt ). Finally, (4.5) together with (4.16) and |ΔXT|2≤sup0≤t≤T|ΔXt|2E| X_T|^2≤ _0≤ t≤ TE| X_t|^2 yields sup0≤t≤T|ΔYt|2≤C(|η|2+∫0T(|αt|2+|βt|2+|γt|2+|t−t|2)t). _0≤ t≤ TE| Y_t|^2≤ C (E|η|^2+E _0^T(| _t|^2+| _t|^2+| _t|^2+|Y_t-U_t|^2)dt ). Thus the desired estimate also holds in Case 2. Case 3: m=nm=n (square case). When m=nm=n, the matrix G is invertible. If β2>0 _2>0, the argument of Case 2 applies; if β2=0 _2=0, then Assumption 1 forces β1>0 _1>0 and μ>0μ>0, and Case 1 applies. Hence the estimate follows. Combining all cases completes the proof of Theorem 3.1. ∎ 4.2 Proof of Theorem 3.2 Proof. We only prove (3.6), since (3.7) follows immediately from (3.6) by evaluating the estimate at the grid points. For each i=0,…,N−1i=0,…,N-1, set bi:=b(ti,X^i,U^i,Z^i),σi:=σ(ti,X^i,U^i,Z^i),fi:=f(ti,X^i,Y^i,Z^i).b_i:=b(t_i, X_i, U_i, Z_i), _i:=σ(t_i, X_i, U_i, Z_i), f_i:=f(t_i, X_i, Y_i, Z_i). By the martingale representation theorem, there exist progressively measurable processes ρX,iρ^X,i and ρY,iρ^Y,i on [ti,ti+1][t_i,t_i+1] such that RiX=ti[RiX]+∫titi+1ρsX,iBs,R_i^X=E_t_i[R_i^X]+ _t_i^t_i+1 _s^X,i\,dB_s, and RiY=ti[RiY]+∫titi+1ρsY,iBs.R_i^Y=E_t_i[R_i^Y]+ _t_i^t_i+1 _s^Y,i\,dB_s. Moreover, ∫titi+1|ρsX,i|2s≤|RiX|2,∫titi+1|ρsY,i|2s≤|RiY|2.E _t_i^t_i+1| _s^X,i|^2\,ds |R_i^X|^2, _t_i^t_i+1| _s^Y,i|^2\,ds |R_i^Y|^2. Define, for t∈[ti,ti+1]t∈[t_i,t_i+1], t:=X^i+∫tit(bi+ti[RiX]Δti)s+∫tit(σi+ρsX,i)Bs,X_t:= X_i+ _t_i^t (b_i+ E_t_i[R_i^X] t_i )ds+ _t_i^t ( _i+ _s^X,i )dB_s, and t:=Y^i+∫tit(−fi+ti[RiY]Δti)s+∫tit(Z^i+ρsY,i)Bs.Y_t:= Y_i+ _t_i^t (-f_i+ E_t_i[R_i^Y] t_i )ds+ _t_i^t ( Z_i+ _s^Y,i )dB_s. Set t:=Z^i+ρtY,i,t:=U^i,t∈[ti,ti+1].Z_t:= Z_i+ _t^Y,i, _t:= U_i, t∈[t_i,t_i+1]. Then, by the definitions of RiXR_i^X and RiYR_i^Y, we have ti=X^i,ti+1=X^i+1,X_t_i= X_i, _t_i+1= X_i+1, and ti=Y^i,ti+1=Y^i+1.Y_t_i= Y_i, _t_i+1= Y_i+1. In particular, 0=a,T=Y^N=g(X^N)−η=g(T)−η.X_0=a, _T= Y_N=g( X_N)-η=g(X_T)-η. Equivalently, YT−T=g(XT)−g(T)+η.Y_T-Y_T=g(X_T)-g(X_T)+η. We now write (,,,)(X,Y,Z,U) as a perturbed continuous-time system. On t∈[ti,ti+1]t∈[t_i,t_i+1], define αt:=bi+ti[RiX]Δti−b(t,t,t,t), _t:=b_i+ E_t_i[R_i^X] t_i-b(t,X_t,U_t,Z_t), βt:=σi+ρtX,i−σ(t,t,t,t), _t:= _i+ _t^X,i-σ(t,X_t,U_t,Z_t), and γt:=fi−f(t,t,t,t)−ti[RiY]Δti. _t:=f_i-f(t,X_t,Y_t,Z_t)- E_t_i[R_i^Y] t_i. Then (,,)(X,Y,Z) satisfies dt=[b(t,t,t,t)+αt]dt+[σ(t,t,t,t)+βt]dBt,dX_t= [b(t,X_t,U_t,Z_t)+ _t ]dt+ [σ(t,X_t,U_t,Z_t)+ _t ]dB_t, and −dt=[f(t,t,t,t)+γt]dt−tdBt,-dY_t= [f(t,X_t,Y_t,Z_t)+ _t ]dt-Z_t\,dB_t, with terminal mismatch η. Hence Theorem 3.1 can be applied to compare (X,Y,Z)(X,Y,Z) and (,,)(X,Y,Z). It remains to estimate the perturbation terms. By the Lipschitz continuity in (x,y,z)(x,y,z), the 1/21/2-Hölder continuity in time, and the definitions above, for t∈[ti,ti+1]t∈[t_i,t_i+1], |αt|2 | _t|^2 ≤C(|t−ti|+|t−X^i|2+|t−Z^i|2+|ti[RiX]|2(Δti)2), ≤ C (|t-t_i|+|X_t- X_i|^2+|Z_t- Z_i|^2+ |E_t_i[R_i^X]|^2( t_i)^2 ), |βt|2 | _t|^2 ≤C(|t−ti|+|t−X^i|2+|t−Z^i|2+|ρtX,i|2), ≤ C (|t-t_i|+|X_t- X_i|^2+|Z_t- Z_i|^2+| _t^X,i|^2 ), |γt|2 | _t|^2 ≤C(|t−ti|+|t−X^i|2+|t−Y^i|2+|t−Z^i|2+|ti[RiY]|2(Δti)2). ≤ C (|t-t_i|+|X_t- X_i|^2+|Y_t- Y_i|^2+|Z_t- Z_i|^2+ |E_t_i[R_i^Y]|^2( t_i)^2 ). Furthermore, t−Z^i=ρtY,i.Z_t- Z_i= _t^Y,i. By Assumption 3, together with the linear growth of b,σ,fb,σ,f (which follows from their Lipschitz continuity), this implies supi[|bi|2+|σi|2+|fi|2]≤C. _iE[|b_i|^2+| _i|^2+|f_i|^2]≤ C. Using these bounds and the construction of tX_t, tY_t, a standard computation yields ∑i=0N−1∫titi+1|t−X^i|2t≤C(|π|+ℜπ), _i=0^N-1E _t_i^t_i+1|X_t- X_i|^2\,dt≤ C (|π|+ R_π ), and ∑i=0N−1∫titi+1|t−Y^i|2t≤C(|π|+ℜπ). _i=0^N-1E _t_i^t_i+1|Y_t- Y_i|^2\,dt≤ C (|π|+ R_π ). Also, ∑i=0N−1∫titi+1|ρtX,i|2t≤∑i=0N−1|RiX|2≤|π|ℜπ≤Tℜπ, _i=0^N-1E _t_i^t_i+1| _t^X,i|^2\,dt≤ _i=0^N-1E|R_i^X|^2≤|π| R_π≤ T R_π, and similarly, ∑i=0N−1∫titi+1|ρtY,i|2t≤Tℜπ. _i=0^N-1E _t_i^t_i+1| _t^Y,i|^2\,dt≤ T R_π. Moreover, ∑i=0N−1∫titi+1|ti[RiX]|2(Δti)2t≤∑i=0N−1|RiX|2Δti, _i=0^N-1E _t_i^t_i+1 |E_t_i[R_i^X]|^2( t_i)^2\,dt≤ _i=0^N-1E |R_i^X|^2 t_i, and ∑i=0N−1∫titi+1|ti[RiY]|2(Δti)2t≤∑i=0N−1|RiY|2Δti. _i=0^N-1E _t_i^t_i+1 |E_t_i[R_i^Y]|^2( t_i)^2\,dt≤ _i=0^N-1E |R_i^Y|^2 t_i. Combining the preceding estimates gives ∫0T(|αt|2+|βt|2+|γt|2)t≤C(ℜπ+|π|).E _0^T (| _t|^2+| _t|^2+| _t|^2 )dt≤ C ( R_π+|π| ). Moreover, since t=U^iU_t= U_i and Y¯t=Y^i Y_t= Y_i on [ti,ti+1)[t_i,t_i+1), ∫0T|t−t|2t _0^T|Y_t-U_t|^2dt ≤C∫0T|t−Y¯t|2t+C∫0T|Y¯t−U¯t|2t ≤ CE _0^T|Y_t- Y_t|^2dt+CE _0^T| Y_t- U_t|^2dt ≤C(|π|+ℜπ)+Cπ. ≤ C (|π|+ R_π )+C P_π. Applying Theorem 3.1, we therefore obtain sup0≤t≤T[|Xt−t|2+|Yt−t|2]+∫0T|Zt−t|2t≤C(|η|2+ℜπ+π+|π|). _0≤ t≤ TE [|X_t-X_t|^2+|Y_t-Y_t|^2 ]+E _0^T|Z_t-Z_t|^2dt≤ C (E|η|^2+ R_π+ P_π+|π| ). It remains to pass from the continuous interpolation (,,)(X,Y,Z) to the piecewise constant interpolation (X¯,Y¯,Z¯)( X, Y, Z). From the estimates above, sup0≤t≤T|t−X¯t|2+sup0≤t≤T|t−Y¯t|2≤C(|π|+ℜπ), _0≤ t≤ TE|X_t- X_t|^2+ _0≤ t≤ TE|Y_t- Y_t|^2≤ C (|π|+ R_π ), and ∫0T|t−Z¯t|2t=∑i=0N−1∫titi+1|ρtY,i|2t≤Cℜπ.E _0^T|Z_t- Z_t|^2dt= _i=0^N-1E _t_i^t_i+1| _t^Y,i|^2dt≤ C R_π. Consequently, by the triangle inequality, sup0≤t≤T[|Xt−X¯t|2+|Yt−Y¯t|2]+∫0T|Zt−Z¯t|2t _0≤ t≤ TE [|X_t- X_t|^2+|Y_t- Y_t|^2 ]+E _0^T|Z_t- Z_t|^2dt ≤C(|η|2+ℜπ+π+|π|). ≤ C (E|η|^2+ R_π+ P_π+|π| ). This proves (3.6). The grid-point estimate (3.7) follows since X¯ti=X^i X_t_i= X_i and Y¯ti=Y^i Y_t_i= Y_i. The proof is complete. ∎ 5 Numerical Experiments In this section we present numerical experiments designed to illustrate the a posteriori error estimates developed in Section 3. The purpose of the experiments is not to prove convergence rates for a particular neural network architecture, but rather to examine whether the computable quantities appearing in Theorem 3.2, namely [|η|2],ℜπ,π,E[|η|^2], R_π, P_π, provide meaningful diagnostic information for deep learning approximations of fully coupled FBSDEs. The experiments are organized as follows. In the first example, an explicit solution is available through a Riccati equation. This allows us to compare the neural network approximation with a reference solution and to check whether the a posteriori indicators are consistent with the observed approximation error. In the second example, we consider a Burgers-type fully coupled FBSDE for which no closed-form solution is used. This example is intended to test whether the three components of the a posteriori loss can serve as computable diagnostics and whether removing some of them leads to less stable or less consistent numerical approximations. 5.1 Network Architecture and Training Configuration We now describe how the abstract discrete approximation appearing in Theorem 3.2 is realized by neural network parameterizations. Let π=0=t0<t1<⋯<tN=T,Δt=TN,π=\0=t_0<t_1<·s<t_N=T\, t= TN, be a uniform partition of the time interval. The numerical method constructs a discrete adapted trajectory (X^i,Y^i,Z^i,U^i)i=0N−1,(X^N,Y^N),\( X_i, Y_i, Z_i, U_i)\_i=0^N-1, ( X_N, Y_N), where the auxiliary process U^i U_i is used in the forward coefficients and is allowed to differ from Y^i Y_i. We introduce three families of feedforward neural networks: i(⋅;θiY) _i(·; _i^Y) :ℝn→ℝm,i=0,1,…,N, :R^n ^m, i=0,1,…,N, (5.1) i(⋅;θiZ) _i(·; _i^Z) :ℝn→ℝm×d,i=0,1,…,N−1, :R^n ^m× d, i=0,1,…,N-1, (5.2) i(⋅;θiU) _i(·; _i^U) :ℝn→ℝm,i=0,1,…,N−1. :R^n ^m, i=0,1,…,N-1. (5.3) Here the networks i\Z_i\ and i\U_i\ serve as standard approximations of the control process and the auxiliary control entering the forward dynamics, respectively. The backward component is, however, constructed recursively via residual networks. Specifically, 0Y_0 acts as an initial value network, while for i≥1i≥ 1, each iY_i directly outputs the one-step backward residual Ri−1YR_i-1^Y defined in (3.4). This parameterization exposes the pathwise residual exactly where it is needed for the a posteriori loss. Concretely, given simulated Brownian increments ΔBi B_i, the discrete trajectory is generated by Y^0 Y_0 =0(X^0;θ0Y), =Y_0( X_0; _0^Y), (5.4) Z^i Z_i =i(X^i;θiZ),U^i=i(X^i;θiU). =Z_i( X_i; _i^Z), U_i=U_i( X_i; _i^U). (5.5) X^i+1 X_i+1 =X^i+b(ti,X^i,U^i,Z^i)Δt+σ(ti,X^i,U^i,Z^i)ΔBi, = X_i+b(t_i, X_i, U_i, Z_i) t+σ(t_i, X_i, U_i, Z_i) B_i, (5.6) Y^i+1 Y_i+1 =Y^i−f(ti,X^i,Y^i,Z^i)Δt+Z^iΔBi+i+1(X^i+1;θi+1Y). = Y_i-f(t_i, X_i, Y_i, Z_i) t+ Z_i\, B_i+Y_i+1( X_i+1; _i+1^Y). (5.7) Comparing (5.7) with (3.4) shows that RiY=i+1(X^i+1;θi+1Y),R_i^Y=Y_i+1( X_i+1; _i+1^Y), i.e., the residual network i+1Y_i+1 directly learns the defect in the backward dynamics. The terminal defect is η=g(X^N)−Y^N.η=g( X_N)- Y_N. The empirical training loss is chosen in accordance with Theorem 3.2 and takes the three-component form ℒ=ℒT+λRℒR+λUℒU,L=L_T+ _RL_R+ _UL_U, (5.8) where ℒT _T =^[|g(X^N)−Y^N|2], = E [|g( X_N)- Y_N|^2 ], (5.9) ℒR _R =∑i=0N−1^[|i+1(X^i+1)|2Δt]=∑i=0N−1^[|RiY|2Δt], = _i=0^N-1 E [ |Y_i+1( X_i+1)|^2 t ]= _i=0^N-1 E [ |R_i^Y|^2 t ], (5.10) ℒU _U =∑i=0N−1^[|Y^i−U^i|2]Δt. = _i=0^N-1 E [| Y_i- U_i|^2 ] t. (5.11) The constants λR _R and λU _U are penalty weights; they are not to be confused with the coupling parameters of the FBSDE itself. The loss (5.8) is a Monte Carlo approximation of the computable terms appearing in the a posteriori bound, up to the deterministic discretization term |π||π|. The complete training procedure, which integrates the forward simulation, the computation of the three-component loss (5.8), and the parameter updates, is summarized in Algorithm 1. At each iteration, a batch of Brownian motion paths is generated, the forward SDE is solved using the current networks i,i,i\Y_i,Z_i,U_i\, and the empirical loss ℒL is evaluated by Monte Carlo averaging. All network parameters are then updated by a stochastic gradient descent variant (Adam). The algorithm corresponds exactly to the discrete scheme described in Section 3.2, and every term in the loss function mirrors a computable component of the a posteriori bound in Theorem 3.2. Algorithm 1 Deep FBSDE with control mismatch Input: Time grid tii=0N\t_i\_i=0^N, number of sample paths M, penalty weights λR,λU _R, _U, network architectures for ii=0N\Y_i\_i=0^N, ii=0N−1\Z_i\_i=0^N-1, ii=0N−1\U_i\_i=0^N-1. Initialize: All network parameters θ. for iter = 1 to max_iter do Generate M independent Brownian increments ΔBi(m)i=0N−1\ B_i^(m)\_i=0^N-1. for m=1,…,Mm=1,…,M do X^0(m)←a X_0^(m)← a Y^0(m)←0(X^0(m);θ0Y) Y_0^(m) _0( X_0^(m); _0^Y) ℒR(m)←0L_R^(m)← 0, ℒU(m)←0L_U^(m)← 0 for i=0,…,N−1i=0,…,N-1 do Z^i(m)←i(X^i(m);θiZ) Z_i^(m) _i( X_i^(m); _i^Z), U^i(m)←i(X^i(m);θiU) U_i^(m) _i( X_i^(m); _i^U) X^i+1(m)←X^i(m)+b(ti,X^i(m),U^i(m),Z^i(m))Δt+σ(ti,X^i(m),U^i(m),Z^i(m))ΔBi(m) X_i+1^(m)← X_i^(m)+b(t_i, X_i^(m), U_i^(m), Z_i^(m)) t+σ(t_i, X_i^(m), U_i^(m), Z_i^(m)) B_i^(m) ei(m)←i+1(X^i+1(m);θi+1Y)e_i^(m) _i+1( X_i+1^(m); _i+1^Y) Y^i+1(m)←Y^i(m)−f(ti,X^i(m),Y^i(m),Z^i(m))Δt+Z^i(m)ΔBi(m)+ei(m) Y_i+1^(m)← Y_i^(m)-f(t_i, X_i^(m), Y_i^(m), Z_i^(m)) t+\; Z_i^(m) B_i^(m)+e_i^(m) ℒR(m)←ℒR(m)+|ei(m)|2/ΔtL_R^(m) _R^(m)+|e_i^(m)|^2/ t ℒU(m)←ℒU(m)+|Y^i(m)−U^i(m)|2ΔtL_U^(m) _U^(m)+| Y_i^(m)- U_i^(m)|^2 t X^i(m)←X^i+1(m) X_i^(m)← X_i+1^(m), Y^i(m)←Y^i+1(m) Y_i^(m)← Y_i+1^(m) end for ℒT(m)←|g(X^N(m))−Y^N(m)|2L_T^(m)←|g( X_N^(m))- Y_N^(m)|^2 end for ℒ←1M∑m=1M(ℒT(m)+λRℒR(m)+λUℒU(m))L← 1M _m=1^M (L_T^(m)+ _RL_R^(m)+ _UL_U^(m) ) Update all parameters θ using Adam optimizer on ∇θℒ _θL end for return optimized parameters θ. 5.2 Example 1: An FBSDE with Explicit Solution We first consider a multi-dimensional fully coupled FBSDE with a linear-quadratic structure. The exact solution is available through the associated Riccati equation, which makes this example suitable for checking the relationship between the computable indicators and the true approximation error. The system is dXti=(−2Xti+Yti)dt+(3Xti+Zti)dBti,−dYti=(−Xti−2Yti+3Zti)dt−ZtidBti,X0=a,YT=−QXT, casesdX_t^i= (-2X_t^i+Y_t^i )\,dt+ (3X_t^i+Z_t^i )\,dB_t^i,\\[3.0pt] -dY_t^i= (-X_t^i-2Y_t^i+3Z_t^i )\,dt-Z_t^i\,dB_t^i,\\[3.0pt] X_0=a, Y_T=-QX_T, cases (5.12) where Xt=(Xt1,…,Xtn)⊤,Yt=(Yt1,…,Ytn)⊤,Zt=diag(Zt1,…,Ztn),X_t=(X_t^1,…,X_t^n) , Y_t=(Y_t^1,…,Y_t^n) , Z_t=diag(Z_t^1,…,Z_t^n), and Bt=(Bt1,…,Btn)⊤B_t=(B_t^1,…,B_t^n) is an n-dimensional Brownian motion. In this example, we set the dimension n=100n=100 for the numerical experiments, the initial condition a=na=1_n (the all-ones vector), the terminal time T=0.1T=0.1, and the terminal coefficient Q=5.0InQ=5.0\,I_n. The time horizon is partitioned into N equal subintervals. Unless otherwise stated, we take N=20N=20 as the default grid; variations of N (e.g., N=10,30,40,50N=10,30,40,50) are explicitly indicated when analyzing the discretization effect in Figure 2 and Table 1. We seek a solution of the form Yt=−KtXt,Zt=diag(−MtXt).Y_t=-K_tX_t, Z_t=diag(-M_tX_t). Substitution into (5.12) gives the Riccati system K˙t−Kt2−4Kt+3Mt+In=0,−3Kt+KtMt+Mt=0,KT=Q. cases K_t-K_t^2-4K_t+3M_t+I_n=0,\\ -3K_t+K_tM_t+M_t=0,\\ K_T=Q. cases (5.13) With the above parameters, a Runge–Kutta solver for (5.13) gives Y0=−2.8730 1n.Y_0=-2.8730\,1_n. This value is used as the benchmark for the numerical approximation. Figure 1 reports the training dynamics for a representative run. The terminal defect, the pathwise residual, and the control mismatch all decrease during training. The error in the initial value Y0Y_0 also decreases over the training iterations. This indicates that, for a fixed time grid and a fixed network configuration, the computable loss is consistent with the improvement of the approximation. We emphasize, however, that the absolute value of the loss should not be interpreted as a direct pointwise estimate of the error, because it is affected by the time discretization, the accumulation of local residuals, the penalty weights, and the optimization difficulty. (a) Training loss (b) Y0Y_0 error Figure 1: Training dynamics of the deep FBSDE solver. Panel (a) shows the total loss ℒL and its components ℒTL_T, ℒRL_R, and ℒUL_U on a logarithmic scale. Panel (b) shows the absolute error |Y0N−Y0||Y_0 N-Y_0|. The results suggest that the decrease of the computable a posteriori loss is accompanied by an improvement in the estimated initial value. The influence of the time discretization is examined in Figure 2 and Table 1. The table reports the terminal indicator, the residual indicator, the mismatch indicator, the total loss, and the absolute error in Y0Y_0, averaged over five independent runs. Figure 2: Absolute error |Y0N−Y0||Y_0^N-Y_0| for different numbers of time steps N, averaged over five independent runs. The reference value is obtained from the Riccati equation. After N=20N=20 the error stabilizes, and further grid refinement does not lead to a significant improvement in this experiment. Table 1: A posteriori indicators and Y0Y_0 error for different N. Mean and standard deviation over five runs are reported. N η^2 η^2 ℜ^π R_π ^π P_π Total loss |Y0N−Y0||Y_0^N-Y_0| 10 8.191 ± 0.228 0.454 ± 0.036 33.75 ± 1.061 42.39 ± 1.191 0.017027 ± 0.0074 20 17.53 ± 0.611 1.354 ± 0.083 81.29 ± 1.400 100.2 ± 1.918 0.005288 ± 0.0042 30 31.28 ± 1.563 2.679 ± 0.040 128.6 ± 1.671 162.6 ± 3.090 0.005947 ± 0.0052 40 47.82 ± 2.407 3.909 ± 0.188 164.8 ± 1.464 216.5 ± 3.805 0.005601 ± 0.0018 50 62.52 ± 4.394 5.005 ± 0.263 197.1 ± 3.290 264.6 ± 7.137 0.004788 ± 0.0038 Several observations can be made. First, the total computable loss increases with N in this experiment. This does not contradict the a posteriori estimate, since the empirical loss is a sum of local residual-type quantities and its magnitude is also influenced by the number of time steps, the number of networks to be trained, and the optimization landscape. Second, the error in Y0Y_0 drops markedly from N=10N=10 to N=20N=20 and then remains roughly stable, fluctuating around 0.0050.005 without a clear monotone trend for larger N. This indicates that the time discretization error is already well controlled with a moderate number of time steps, while further grid refinement does not bring a significant additional benefit under the present network and training setup. A finer grid reduces the local discretization error in principle, but it also increases the number of neural network components and may make the optimization problem more difficult. The results also confirm that the a posteriori indicators should be interpreted as diagnostic tools rather than as a simple monotone ranking across different discretizations. For a fixed discretization, their decrease reflects better satisfaction of the terminal condition, the backward dynamics, and the control consistency. Across different discretizations, they help reveal the interplay between discretization accuracy and optimization difficulty. Figure 3 compares the true and approximated trajectories of YtY_t and ZtZ_t along one representative Brownian path for a 5-dim case with N=50N=50. The value process Y is approximated very accurately along this path, whereas the control process Z exhibits a more visible discrepancy. This is consistent with the fact that Z is a gradient-type quantity and is usually more sensitive to local fluctuations of the Brownian path. The result supports the use of residual and mismatch indicators as useful diagnostics, but one should not infer full pathwise accuracy from a single trajectory alone. (a) Trajectory of YtY_t (b) Trajectory of ZtZ_t Figure 3: Comparison of the true trajectories and the neural network approximations for one representative path in the five-dimensional case with N=50N=50. The Y component is accurately reproduced along this path, while the Z component is more difficult to approximate. 5.3 Example 2: FBSDE Associated with a Burgers-Type PDE We next consider a Burgers-type fully coupled FBSDE with a tunable feedback strength. Let Xt∈ℝnX_t ^n, Yt∈ℝY_t , and Zt∈ℝnZ_t ^n. For a viscosity parameter ν>0ν>0, a damping coefficient λ∈ℝλ , and a coupling strength ρ≥0ρ≥ 0, consider dXt=−ρYtndt+2νdBt,−dYt=λYtdt−Zt⊤dBt,X0=x0,YT=g(XT), casesdX_t=-ρ Y_t1_n\,dt+ 2ν\,dB_t,\\[4.0pt] -dY_t=λ Y_t\,dt-Z_t dB_t,\\[4.0pt] X_0=x_0, Y_T=g(X_T), cases (5.14) where BtB_t is an n-dimensional Brownian motion. The parameter ρ controls the strength of the feedback from the backward component YtY_t into the forward dynamics. When ρ=0ρ=0, the forward equation is decoupled from Y; when ρ>0ρ>0, the forward dynamics depend directly on the backward component. Assume that the FBSDE admits a smooth decoupling field u:[0,T]×ℝn→ℝu:[0,T]×R^n such that Yt=u(t,Xt),Zt=2ν∇xu(t,Xt).Y_t=u(t,X_t), Z_t= 2ν\, _xu(t,X_t). By Itô’s formula, u satisfies ∂tu−ρu 1n⋅∇xu+νΔxu+λu=0,(t,x)∈[0,T)×ℝn,u(T,x)=g(x). cases _tu-ρ u\,1_n· _xu+ν _xu+λ u=0, (t,x)∈[0,T)×R^n,\\[4.0pt] u(T,x)=g(x). cases (5.15) Equivalently, after the time reversal v(τ,x)=u(T−τ,x)v(τ,x)=u(T-τ,x), one obtains the viscous damped Burgers-type equation [CS11, BGN14] ∂τv+ρv 1n⋅∇xv=νΔxv+λv. _τv+ρ v\,1_n· _xv=ν _xv+λ v. In the numerical experiment, we set n=10,g(x)=sin(∑i=110xi),T=1,n=10, g(x)= \! ( _i=1^10x_i ), T=1, and x0=,ν=0.2,λ=0.5,ρ=0.5.x_0=0, ν=0.2, λ=0.5, ρ=0.5. The time interval [0,T][0,T] is uniformly partitioned into N subintervals. Unless otherwise stated, we take N=20N=20 as the default grid. Since no closed-form solution is available for this configuration, this example is not intended to measure the true approximation error directly. Instead, it examines whether the three computable components ℒTL_T, ℒRL_R, and ℒUL_U provide useful information about the consistency and stability of the trained solution. We consider the following configurations: • Full loss: all three terms ℒTL_T, ℒRL_R, and ℒUL_U are included. The penalty weights λR _R and λU _U are varied to examine the sensitivity of the method. • w/o R: the pathwise residual term ℒRL_R is removed, while the control mismatch term ℒUL_U is retained. • w/o RU: both ℒRL_R and ℒUL_U are removed, leaving only the terminal loss. The two ablation configurations are used to assess the numerical role of the residual and mismatch terms. Since no reference solution is available, the comparison should be interpreted in terms of stability, consistency, and sensitivity of the estimates, rather than in terms of exact accuracy. Figures 4 and 5 display the training loss curves for different choices of the penalty weights. Each panel reports the total loss and its components on a logarithmic scale. The loss curves show that the three components can be reduced simultaneously under the full-loss configuration. This indicates that the trained networks can approximately satisfy the terminal condition, the discrete backward dynamics, and the control consistency condition at the same time. The final loss levels vary with the penalty weights, which is expected because λR _R and λU _U change the relative strength of the corresponding soft constraints. Figure 4: Training loss curves for different residual penalty weights λR _R. Each panel shows the total loss together with ℒTL_T, ℒRL_R, and ℒUL_U. Figure 5: Training loss curves for different mismatch penalty weights λU _U. The layout is the same as in Figure 4. Tables 2 and 3 report the estimated values of Y0Y_0 and a scalar summary of Z0Z_0 for different configurations. Since Z0∈ℝnZ_0 ^n in this example, the reported value of Z0Z_0 should be understood as the first value of its components in this example, namelyZ¯0:=Z01 Z_0:=Z_0^1. Table 2: Ablation study: estimated Y0Y_0 and scalar summary Z¯0 Z_0 for the full-loss configuration with different residual penalty weights λR _R, and for two ablation configurations. “w/o R” removes the pathwise residual ℒRL_R; “w/o RU” removes both ℒRL_R and the control mismatch ℒUL_U. Mean and standard deviation over five runs are reported. λR _R / Configuration Y0Y_0 (mean ± std) Z¯0 Z_0 (mean ± std) 0.1 0.256436 ± 0.003725 0.110142 ± 0.010282 0.5 0.242612 ± 0.003119 0.115508 ± 0.009206 1.0 0.237267 ± 0.006471 0.122293 ± 0.009611 2.0 0.229739 ± 0.003257 0.120487 ± 0.012230 w/o R 0.227853 ± 0.011740 0.135288 ± 0.009247 w/o RU 0.233889 ± 0.002582 0.135016 ± 0.003326 Table 3: Ablation study: estimated Y0Y_0 and scalar summary Z¯0 Z_0 for the full-loss configuration with different mismatch penalty weights λU _U, and for the baseline “w/o RU”. Mean and standard deviation over five runs are reported. λU _U / Configuration Y0Y_0 (mean ± std) Z¯0 Z_0 (mean ± std) 0.1 0.233716 ± 0.005595 0.138136 ± 0.006634 0.5 0.234612 ± 0.004430 0.119596 ± 0.012485 1.0 0.233662 ± 0.003395 0.120903 ± 0.004734 2.0 0.221752 ± 0.003079 0.086442 ± 0.018901 w/o RU 0.233889 ± 0.002582 0.135016 ± 0.003326 (a) Y0Y_0 vs. λR _R (b) Z¯0 Z_0 vs. λR _R Figure 6: Dependence of the estimated Y0Y_0 and Z¯0 Z_0 on the residual penalty weight λR _R. (a) Y0Y_0 vs. λU _U (b) Z¯0 Z_0 vs. λU _U Figure 7: Dependence of the estimated Y0Y_0 and Z¯0 Z_0 on the mismatch penalty weight λU _U. The ablation results lead to the following observations. First, the full-loss configuration produces stable estimates over a range of penalty weights. The estimated Y0Y_0 and the scalar summary of Z0Z_0 change with λR _R and λU _U, but the variation is moderate for most parameter choices. This suggests that the method is not overly sensitive to a single choice of penalty weight, although excessively large or small weights may still affect the balance between the different constraints. Second, removing the pathwise residual or the mismatch penalty changes the estimated control component more visibly than the estimated value component. In particular, the ablation configurations tend to produce larger values of the reported Z0Z_0 summary. Since no exact solution is available in this example, this should not be described as a direct loss of accuracy. A more appropriate interpretation is that the residual and mismatch penalties influence the consistency of the learned control process with the discrete FBSDE dynamics. Third, the configuration without the residual term shows a larger variability in the estimate of Y0Y_0 in Table 2. This supports the view that the residual term plays a stabilizing role in the training process. The terminal loss alone can enforce the final condition, but it does not directly control the intermediate backward dynamics. The residual term therefore acts as an additional soft constraint that helps prevent the learned trajectory from fitting the terminal condition while deviating from the FBSDE structure in the interior of the time interval. Overall, the Burgers-type experiment provides numerical evidence that the three-term loss suggested by the a posteriori estimate is useful in practice. In the absence of a reference solution, the indicators cannot prove the exact accuracy of the approximation, but they provide computable diagnostics for terminal consistency, dynamic consistency, and control consistency. This is precisely the practical role expected from the a posteriori framework. 6 Conclusion In this paper, we developed an a posteriori error estimation framework for deep learning approximations of fully coupled forward-backward stochastic differential equations. The analysis allows the forward coefficients to depend on an auxiliary control process that may differ from the backward component, a formulation motivated by the decoupled training scheme in the second algorithm of [JPP+20] and by practical stabilization techniques such as target networks. By introducing an explicit pathwise residual constraint ℜπ R_π into this decoupled framework, we obtained a fully computable a posteriori error bound that covers the terminal defect, the pathwise residual, and the control mismatch. The main theoretical contribution is a stability estimate for fully coupled FBSDEs under simultaneous perturbations of the drift, diffusion, generator, terminal condition, and control input. Based on this continuous-time stability result, we derived a discrete-time a posteriori error estimate for arbitrary adapted approximations. The resulting bound depends on three computable quantities: the terminal defect, the pathwise residual, and the control mismatch. These quantities can be evaluated from the numerical solution itself and do not require prior knowledge of the exact solution. In this sense, the estimate provides a genuine a posteriori diagnostic for deep FBSDE solvers. The numerical experiments support the relevance of the theoretical framework. For the linear-quadratic example with an explicit Riccati solution, the computable indicators decrease during training and are consistent with the improvement of the estimated initial value. The experiments with different time steps also show that a finer grid does not automatically lead to a better neural approximation, because the increased number of time steps may make the optimization problem more difficult. For the Burgers-type example without a closed-form reference solution, the ablation studies indicate that the pathwise residual and the control mismatch terms play important stabilizing roles. Although the absence of an exact solution prevents a direct measurement of the true error, the indicators provide useful information about terminal consistency, dynamic consistency, and control consistency. Several questions remain open. The present analysis does not quantify how the choice of neural network architecture, the number of trainable parameters, the optimization algorithm, or the Monte Carlo sample size affects the constants in the error estimates. A systematic investigation of these issues would help connect the present a posteriori theory with more explicit convergence rate results. It would also be interesting to extend the framework to FBSDEs with mean-field interaction, reflecting boundaries, or more general path-dependent coefficients. These topics are left for future research. References [AAO25] K. Andersson, A. Andersson, and C. W. Oosterlee (2025) The deep multi-FBSDE method: a robust deep learning method for coupled FBSDEs. arXiv:2503.13193. Cited by: §1. [BHL+19] A. Bachouch, C. Huré, N. Langrené, and H. Pham (2019) Deep neural networks algorithms for stochastic control problems on finite horizon: numerical applications. Methodology and Computing in Applied Probability 24 (1), p. 143–178. Cited by: §1. [BEL57] R. E. Bellman (1957) Dynamic programming. Princeton University Press. Cited by: §1. [BS12] C. Bender and J. Steiner (2012) Least-squares monte carlo for backward sdes. In Numerical Methods in Finance, R. A. Carmona, P. Del Moral, P. Hu, and N. Oudjane (Eds.), Berlin, Heidelberg, p. 257–289. External Links: ISBN 978-3-642-25746-9 Cited by: §1. [BZ08] C. Bender and J. Zhang (2008) Time discretization and markovian iteration for coupled fbsdes. Annals of Applied Probability 18 (1), p. 143–177. Cited by: §1. [BET04] B. Bouchard, I. Ekeland, and N. Touzi (2004) On the Malliavin approach to Monte Carlo approximation of conditional expectations. Finance and Stochastics 8 (1), p. 45–71. External Links: Document Cited by: §1. [BT04] B. Bouchard and N. Touzi (2004) Discrete-time approximation and Monte-Carlo simulation of backward stochastic differential equations. Stochastic Processes and their Applications 111 (2), p. 175–206. Cited by: §1. [BGN14] Z. Brzeźniak, B. Goldys, and M. Neklyudov (2014) Multidimensional stochastic burgers equation. SIAM Journal on Mathematical Analysis 46 (1), p. 871–889. External Links: Document, Link, https://doi.org/10.1137/120866117 Cited by: §5.3. [CS11] A. B. Cruzeiro and E. Shamarova (2011) On a forward-backward stochastic system associated to the Burgers equation. In Stochastic Analysis with Financial Applications, Progress in Probability, Vol. 65, p. 43–59. External Links: Document Cited by: §5.3. [DE92] D. Duffie and L. G. Epstein (1992) Stochastic differential utility. Econometrica 60 (2), p. 353–394. Cited by: §1. [EHJ17] W. E, J. Han, and A. Jentzen (2017) Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations. Communications in Mathematics and Statistics 5 (4), p. 349–380. Cited by: §1. [EPQ97] N. El Karoui, S. Peng, and M. C. Quenez (1997) Backward stochastic differential equations in finance. Mathematical Finance 7 (1), p. 1–71. External Links: Document, Link, https://onlinelibrary.wiley.com/doi/pdf/10.1111/1467-9965.00022 Cited by: §1. [EQ97] N. El Karoui and M. C. Quenez (1997) Imperfect markets and backward stochastic differential equations. In Numerical Methods in Finance, L. C. G. Rogers and D. Talay (Eds.), p. 181–214. External Links: Document Cited by: §1. [FZZ16] Y. Fu, W. Zhao, and T. Zhou (2016) Multistep schemes for forward backward stochastic differential equations with jumps. Journal of Scientific Computing 69 (2), p. 1–22. Cited by: §1. [GCZ+23] C. Gao, S. Chen, Z. Zhu, and Z. Wang (2023) Convergence of the backward deep BSDE method with applications to optimal stopping problems. SIAM Journal on Financial Mathematics. Cited by: §1. [GMW22] M. Germain, J. Mikael, and X. Warin (2022) Numerical resolution of McKean-Vlasov FBSDEs using neural networks. Methodology and Computing in Applied Probability 24 (4), p. 2557–2586. External Links: Document Cited by: §1. [GOP25] A. Gnoatto, K. Oberpriller, and A. Picarelli (2025) Convergence of a deep BSDE solver with jumps. arXiv:2501.09727. Cited by: §1. [GL08] E. Gobet and C. Labart (2008) Error expansion for the discretization of backward stochastic differential equations. Stochastic Processes and their Applications 118 (5), p. 803–829. Cited by: §1. [GLW05] E. Gobet, J. Lemor, and X. Warin (2005) A regression-based monte carlo method to solve backward stochastic differential equations. Annals of Applied Probability 15 (3), p. 2172–2202. Cited by: §1. [GT16] E. Gobet and P. Turkedjiev (2016) Linear regression MDP scheme for discrete backward stochastic differential equations under general conditions. Mathematics of Computation 85 (299), p. 1359–1391. External Links: Document Cited by: §1. [HJE18] J. Han, A. Jentzen, and W. E (2018) Solving high-dimensional partial differential equations using deep learning. Proceedings of the National Academy of Sciences 115 (34), p. 8505–8510. Cited by: §1. [HL20] J. Han and J. Long (2020) Convergence of the deep BSDE method for coupled FBSDEs. Probability, Uncertainty and Quantitative Risk 5 (1), p. 1–33. Cited by: §1. [HP95] Y. Hu and S. Peng (1995) Solution of forward-backward stochastic differential equations. Probability Theory and Related Fields 103 (2), p. 273–283. Cited by: §1, §2.1. [HPB+21] C. Huré, H. Pham, A. Bachouch, and N. Langrené (2021) Deep neural networks algorithms for stochastic control problems on finite horizon: convergence analysis. SIAM Journal on Numerical Analysis 59 (1), p. 525–557. External Links: Document, Link, https://doi.org/10.1137/20M1316640 Cited by: §1. [JPP+20] S. Ji, S. Peng, Y. Peng, and X. Zhang (2020) Three algorithms for solving high-dimensional fully coupled fbsdes through deep learning. IEEE Intelligent Systems 35 (3), p. 71–84. Cited by: §1, §1, §6, Remark. [MPY94] J. Ma, P. Protter, and J. Yong (1994) Solving forward-backward stochastic differential equations explicitly—a four step scheme. Probability theory and related fields 98 (3), p. 339–359. Cited by: §1. [MY99] J. Ma and J. Yong (1999) Forward-backward stochastic differential equations and their applications. Lecture Notes in Mathematics, Vol. 1702, Springer. Cited by: §1, §2.1. [P90] É. Pardoux and S. Peng (1990) Adapted solution of a backward stochastic differential equation. Systems and Control Letters 14 (1), p. 55–61. Cited by: §1. [PW99a] S. Peng and Z. Wu (1999) Fully coupled forward-backward stochastic differential equations and applications to optimal control. SIAM Journal on Control and Optimization 37 (3), p. 825–843. Cited by: §1, §2.1. [PW99b] S. Peng and Z. Wu (1999) Fully coupled forward-backward stochastic differential equations and applications to optimal control. SIAM Journal on Control and Optimization 37 (3), p. 825–843. Cited by: §2.1. [PEN90] S. Peng (1990) A general stochastic maximum principle for optimal control problems. Siam Journal on Control and Optimization 28 (4), p. 966–979. Cited by: §1. [PEN91] S. Peng (1991) Probabilistic interpretation for systems of quasilinear parabolic partial differential equations. Stochastics and Stochastic Reports 37 (1-2), p. 61–74. Cited by: §1. [RPK19] M. Raissi, P. Perdikaris, and G. E. Karniadakis (2019) Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics 378, p. 686–707. Cited by: §1. [RSZ24] C. Reisinger, W. Stockinger, and Y. Zhang (2024) A posteriori error estimates for fully coupled McKean–Vlasov forward-backward SDEs. IMA Journal of Numerical Analysis 44 (4), p. 2323–2369. External Links: Document Cited by: §1, §1.