Integration Matters for Learning PDEs with Backward SDEs
Authors: Sungje Park
Stephen Tu
Venue: NeurIPS 2025 spotlight ()
Keywords: Backward Stochastic Differential Equations
Partial Differential Equations
Deep Learning
Scientific Machine Learning
Primary area: deep_learning
TL;DR: We show that the performance gap between BSDE and PINNs methods stems from discretization bias in Euler-Maruyama, and propose a Stratonovich-Heun BSDE method that eliminates the bias and restores BSDE performance.

Abstract: Backward stochastic differential equation (BSDE)-based deep learning methods provide an alternative to Physics-Informed Neural Networks (PINNs) for solving high-dimensional partial differential equations (PDEs), offering potential algorithmic advantages in settings such as stochastic optimal control, where the PDEs of interest are tied to an underlying dynamical system. However, standard BSDE-based solvers have empirically been shown to underperform relative to PINNs in the literature. In this paper, we identify the root cause of this performance gap as a discretization bias introduced by the standard Euler-Maruyama (EM) integration scheme applied to one-step self-consistency BSDE losses, which shifts the optimization landscape off target. We find that this bias cannot be satisfactorily addressed through finer step-sizes or multi-step self-consistency losses. To properly handle this issue, we propose a Stratonovich-based BSDE formulation, which we implement with stochastic Heun integration.  We show that our proposed approach completely eliminates the bias issues faced by EM integration.  Furthermore, our empirical results show that our Heun-based BSDE method consistently outperforms EM-based variants and achieves competitive results with PINNs across multiple high-dimensional benchmarks. Our findings highlight the critical role of integration schemes in BSDE-based PDE solvers, an algorithmic detail that has received little attention thus far in the literature.

Integration Matters for Learning PDEs with Backward SDEs
Sungje Park
Stephen Tu University of Southern California stephen. tu@usc.
University of Southern California sungjepa@usc.eem
Abstract
Backward stochastic differential equation (BSDE)-based deep learning methods provide an alternative to Ihysice-Informed Neural Networks (PINNs) for solving high-dimensional partial differential equations (PDEs), offering potential algorith- mic advantages in settings such as stochastic optimal control, where the PDEs interest tied to underlying dynamical system. However, standard BSDE- based solvers have empirically been shown to underperform relative to PINNs the literature. In this paper, we identify the root cause of this performance gap a discretization bias introduced by the standard Fuler-Maruyane (EM) integra- tion scheme applied one-step self-consistency BSDE losses, which shifts the optimization landscape off target. We find that this bias cannot be satisfactorily addressed through finer step-sizes multi-step self-consistency losses. To properly handle this issue, we propose a Stratonovich-based BSDE formulation, which we implement with stochastic Heun integration. We show that our proposed approach completely eliminates the bias issues faced by EM integration. Furthermore, our empirical results show that our Heun-basse BSDE method consistently outperforms E.M-ased variants and achieves competitive results with PINNs across multiple digh-dimensional benchmarks. Our findings highlight the critical role of integration schemes in BSDE-based PDE solvers, an algorithmic detail that has received little attention thus far in the literature.
Introduction
Numerical solutions partial differential equations (PDEs) are foundational to modeling problems across a diverse set of fields in science and engineering. However, due to the curse dimensionality of traditional numerical methods, application of classic solvers to high dimensional PDEs is com- putationally intractable. In recent years, motivated by the success of deep learning methods, both Physics-Informed Neural Networks (PINNs) [1. 2] and backward stochastic differential equation (BSDE) methods [3-5] have emerged as promising alternatives to classic techniques.
Despite the widespread popularity of PINNs methods, in this paper we focus the use of BSDE- based methods for solving high-dimensional PDEs. The key difference between PINNs and BSDE methods is that while PINNs minimize the PDE residual directly randomly sampled collocation points, BSDE methods reformulate PDEs as serward-backweer minimize the discrepancy between predieted and terminal
SDEs (FBSDEs) and simulate
the resulting stochastic processes conditions at the end of the forward SDE trajectory [3], or across an intermediate time-horizon
via self-consistency losses [4, 5]. BSDE methods especially well-suited for high-dimensional problems where there is underlying fynamics-suchs as in stochastic optimal control or quantitative finance-as the crux of these methods involving sampling over stochastic trajectories rather than over bounded spatial domains. Furthermore, BSDE methods offers a significant advantage in problems where the governing equations of the PDE are unknown and can only be accessed through simulation [6], as in model-free optimal control (cf. Appendix C for more details). In contrast, PINNs
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

methods require explicit knowledge of the PDE equations, which may be either be impractical to obtain for various tasks or require separate model learning step within the training pipeline.
Surprisingly, despite the aforementioned benefits of BSDE methods compared with PINNs, a thorough comparison between the two techniques remains largely absent from the literature. One notable exception is recent work by [5] which finds that on several benchmark problems. PDE solutions found by BSDE-based approaches significantly underperform the corresponding PINNS solutions. To address this gap, they propose a hybrid interpolating loss between the PINNs and BSDE losses. While promising, their result has two key disadvantages. First, the underlying cause of the performance gap between BSDE and PINNs methods is not elucidated. Second, their method introduces new hyper-parameter (the horizon- length controlling the level of interpolation) which must be tuned for optimal performance, adding complexity to the already delicate training process [7].
In this work, identify the key source of the performance gap between BSDE and PINNs methods as the standard Fuler-Maruyama (EM) scheme used in BSDE methods for stochastic integration. Although simple to implement, we show that the EM scheme introduces a significant discretization bias in one-step BSDE losses, resulting in a discrepancy between the optimization objective and the true solution. We furthermore show that the EM discretization bias can only be made arbitrarily small by using multi-step BSDE losses, which we show both theoretically and empirically comes at a significant cost in performance. Our analysis reveals the interpolating loss of [5] as a method find (via
tuning) the best suitable horizon length (i.e., number of self-consistency steps).
As alternative, we propose interpreting both the forward and backward SDEs Stratonovich SDEs-as as opposed to Itô and utilizing the stochastic Heun integration scheme for numerical integration. We prove that the of the stochastic Heun method completely eliminates the non- vanishing bias issues which occur in the EM formulation for one-step BSDE losses. This removes all performance tradeoffs in the horizon-length, allowing us to utilize single-step self-consistency losses. The result is a practical BSDE-based algorithm that is competitive with PINNs methods without the need for interpolating losses. Surprisingly, prior to our work the role of stochastic integration has received little attention in the BSDE literature; we hope that our results inspire further algorithmic and implementation level improvements for BSDE solvers.
2 Related Work
In recent years, PINNs [1, 2, 8-10] has emerged as a popular method for solving high dimensional PDEs. PIÑNs methods parameterize the PDE solution as neural network and directly minimize the PDE residual as a loss function, provides a mesh-free free method that can easily incorporate complex boundary conditions and empirical data. However, the PINNS approach remains an incomplete solution and still suffers from notable issues including optimization challenges [1 1-5]], despite concerted effort to remedy these difficulties [7, 11, 12, 16-20). Hence, the application of PINNs as a general purpose solver for complex high-dimensionas PDEs remains an active area of research.
On the other hand, a complementary line of work proposes methods based on BSDEs to solve high-dimensional PDEs [3-5, 21-25]. These approaches reformulate PDEs as forward-back backward SDEs to derive a trajectory- based loss. While the original deep BSDE methods [3, 21, 22] learn separate neural networks to predict both the value and gradient at each discrete time-step. follow up work [4, 5] uses self-consistency. i.e., the residual of stochastic integration along BSDE trajectories. to form loss. In this work, we exclusively focus on self-consistency BSDE losses. as they generalize the original method while allowing for a single network to parameterize the PDE solution for all space and time, similarly to PINNs. Similar self-consistency losses have also been recently utilized to learn solutions to Fokker-Plann PDEs [26-28]. Further discussions on related works and the paper's s relationship to recent BSDE- methods can be found in Appendix B.
The main purpose of our work is understand the performance differences between PINNs and BSDE methods on high- -igeesimensioaal PDEs. The most relevant work to ours is [5], which to the best of our knowledge the only work in the literature that directly compares PINNs and BSDE methods in head-to-head evaluation. As discussed in Section 1, one of our main contributions shows that the gap in performance between PINNs and BSDE methods observed in [5] is due to the choice of stochastic integration. Previous work [29] studying stochastic Runge- Kutta discretizations for BSDEs methods considers only the original BSDE losses instead of self-consistency methods, and hence does not uncover the issues identified in our work.
2

---

3 Background and Problem Setup
We consider learning approximate solutions to the following non-linear boundary value PDEs
1 tr( (x, v'u(z,t) + (f(a,
tx, () == t,
(u(z,)) h1, =0, (3.1)
over domain z E Q Rd E := [0, with boundary conditions (a) u(r = ((x) Vr E Ω and (b) t(x,t) (z, Vr E ä9, E I. Here, : Ω x I + R is candidate PDE solution, f H Rd is a vector-field, h[u] 1,(,,,,u(a 1), Vu(x, ()) for some ho QxIx IR R2 + I captures the non-linear terms, (x, ((,,()g(() 1) Red for g(x,t) Reddd is positive definite matrix-valued function, o: 12 I and di : 0Q I are boundary conditions, and both v and 2² denote spatial gradients and Hessians, respectively. For expositional clarity, we will assume that Q H and drop the second boundary condition (b), noting that all subsequent arguments can be extended in a straightforward manner for bounded domains 52.
Physics-Informed Neural Networks (PINNs). Under the assumption of knowledge of the operator
R[u] and the boundary condition o, the standard PINNs methodology [1. 2, 8- 8-10] for solving (3.1) works by parameterizing the solution u(x, t) in a function class == (2, t) 0 E Q) (c.g.. represents the weights of a neural network), and minimizing the PINNS loss over U1:
LPINNs (0; A) : (=.)))) ((] [z](2, +)) X Ez-p' [1)
(3.2)
where j measure over Ω x I and p' is measure over Ω. The choice of measures u, in addition to the relative weight À are
which must be carefully selected by the user. To simplify exposition further, we will assume that each uo (Ix, ()) lll satisfies u(:, = (e.g., as in [30, 3 31D, and hence the PINNs loss simplifies further to LPINNs(0) (0) = (z..))) (R[up](z) +)
Back ward SDEs and self-consistenge losses. While the PINNs loss has received much attention in the literature, a separate line of work has advocated for an alternative approach to solving PDEs The key idea is that given the forward (Itô) SDE:
based on backward SDEs [3, 4, 21-23].
dx,- fXx, t)dt g(x,,t)dBB
Xo To,
(3.3)
where (Bo) is standard Brownian motion in Rdd, the corresponding backward (11) SDE:
dY, = t,,,, Y,, Zt)t+ + 7 ()d,
((X),
(3.4)
solved by setting Y. = (X,,,) 1) and Z, Vu( ((,X,,1) t), where u is solution to the PDE (3.1); this equivalence is readily shown with Itô's lemma. The relationship between the forward and backward SDE has motivated several different types of BSDE loss functions for solving (3.1). In this work, we focus on BSDE losses based on self-consistency [4, 5J, which uses the residual of stochastic integration along the BSDE trajectories as supervision. Self-consistency losses are more practical than other BSDE variants only one network is required and the weights can be shared across time (unlike e.g., the original BSDE losses [3, 21] which learn N models to predict both Yú and Z, at N discretization points, and require retraining for every new initial condition zn). Specifically, consider the following H-horizon (for N = TT/H E N. and l, = nH) self-consistency BSDE loss.?
LISDE, ( (0) == Erg..H N
Σ n=
(uo (Un(XXm11) - us (X,,,,, (n) - So (ob, 1))
(3.5)
where So t1) == S(" he t),t S ()) t)B, with ho = M[uo](x, t).
and zi ~ Mo is drawn from a distribution pmo over initial conditions for the forward SDE (3.3).
Euler-Maruyama integration. Unlike the PINNS loss (3.2), the BSDE loss (3.5) must be discretized with an appropriate stochastic integrator. The standard choice is to use the Euler-Maruyama (EM) method, selecting N e N. to define a step-size the forward and backward SDEs as follows:
T/N and time-points in := nt, and integrating
X X,, + T (8n, (n) + + (n) Way Wn ~ La), X. = Io.
'We leave consideration of the PINNs loss with non square losses (e.g.. [18) to future work.
We note that the most general form of self-comsistency losses are due to [5] and take on the form FRSDE((:pp) E ma-Ro. (... A7 (u - (o(XX,,,,) So(t,,,))) involving time-pair distribution
p over I², where A. ty ts. We choose to present (3.5) more closely aligns with the discrete losses.
3

---

n-1 1 +ho (T ' "(ggg(Xns tn Yo = us ()),
(3.6)
With this discretization, the k-step EM-BSDE loss (for N/k 1)4 for LBSDE( is:
(8) =
k Σ n=00
(n=(X(n=1))
us (Knk tnk) (#t+1)) -
- 22 (3.7)
2
In the one-step (or single-step) setting 1, (3.7) reduces to the self-consistent BSDE loss from [4], also discussed in [32, 33];
the shorthand LEM.r(0) (0) to denote this setting. We refer the > case generally as multi-step, which is form interpolating loss from [5]. Another notable case is when k = N, which recovers the full-horizon losses used in the original BSDE works [3, 21].
4 Analysis of One-Step Self-Consistency Losses
In this section we conduct analysis of both EM and Heun stochastic integration applied to the one-step self-consistency BSDE loss.
The Hölder space Chh.1_ Let f : M + M', where both M, M,M' M' are subsets of Euclidean space (with possibly different dimension). We say that is (M, M') (C when M, M' are clear) if is both bounded and k-times continuously differentiable on M, and all j-th derivatives of for j é {O.....kt k} are Lipschitz continuous. The Hölder norm |[ (fMMMM) is the smallest bound possible on I and all the Lipschitz constants for D f.j {O,..., k).
4.1 Analysis of Maruya-a for BSDE
We first illustrate the bias when using EM to integrate the single-step consistency loss. To do this, we define the point-wise EM loss at resolution T for a fixed (x, t) E IR x I as:
fEM,
E. (uo (ārer + uz(2,t) - (0, (0,,
(4.1)
Pesr == xf t) (Fg(I,t)ur w ...NOO,1 Ia).
The point-wise EM loss (4.1) is related to the one-step EM-BSDE loss via LEM, (0)
Nr Enzo Ex,, [fEM,- (0, (0,X,, (n,tt)) Our first result shows that the dominant error term (in r) of the loss (4.1) suffers from additive bias term that is introduced as result of the EM integration.
Lemma 4.1. Suppose that f.g are bounded and ug is C We have that
(xH(I,t) V ())( --2 2. REM, (0, (R[ue](z, ())+ 1 + (AAu
(4.2)
where the O()) hides factors depending on d, the bounds on S.g, and |lu Ica..
Lemma can further be used, in conjunction with standard results on the order 1/2 strong con- vergence of EM integration [34], to show the following statement regarding the full loss (0).
Theorem 4.2. Suppose that f.g. ho e CO.1, ug E C and r E 1. We have that:
1 1 [" E [((R(u)(x.)))) 1 ( (H(((())))) v'uo dt 0(7113), (4.3)
LEM.- (0) T
where the 0()) hides constants that depend on d, T, and the Hölder norms of f. 9. hg, and ug-
Theorem 4.2 implies that even if the function class is expressive enough to contain a PDE solution 249 to (3.1) satisfying R[ up. = 0, in general LEM.- (6.) ) inf LEM. (f), and hence optimizing LEM,- (0) can lead to sub-optimal solutions even in the limit of infinite simulated data. Furthermore, this bias cannot be resolved by simply reducing the step-size since the PDE residual term and the bias term are both the same order (cf. (4.3)). We illustrate this with a one-dimensional example in Figure 1a. Proofs for both Lemma 4.1 and Theorem 4.2 are given in Appendix D.2.

---

r=1x10
r=1x10
r=1x10- True
109
Lox, (8)r"
Lien [B]mT
20
T=5x10-1
204
r=1x x=1x10-
x
20
03
(a) Plot (EM,-(0) at (10"a ,000 2,10" ,,10") levels of discretization.
(b) Plot of LHemm
(5 10"1 1001,
?} levels discretization.
Figure 1: A plot of both (e) and THemm (0) at various levels of discretization. The PDE is one dimensional Linear Quadratic Regulator HJB equation, where parameterizes quadratic value function.
4.2 Stratonovich BSDEs and Stochastic Heun Integration
Our next step is to derive new BSDE loss based on Heun integration. The starting point is to interpret the forward SDE as a Stratonovich SDE (in contrast to (3.3) which is an Itô SDE):
(4.4) For any u that satisfies (3.1), we have du(X;, ["), tu(Xi [X; t) dB, with hr[ (r[u[(x, t) := Mu][ t) - tr(H(z, () by the Stratonovich chain rule, which motivates the following H-horizon seff-consisteney Stratonovich BSDE loss:
:x; = f(X;,
g(G;,) t) dB,,
Xi - 20-
N-1
(S.BSDE,H(O) Exg, B. N
Σ n=0
(us
uo( (n) S (mn))))=)))) (4.5)
with (to, : 1.th (X; leh (X dB, where hig t) == (<[ f). As (4.5) utilizes Stratonovich integration, the as it converges to the Itô solution. Hence. will consider the stochastic Heun integrator [34. 35] which has the favorable property of converging to the Stratonovich solution. We proceed first by defining the augmented forward and backward SDE process Z
scheme cannot be used for integration,
d
X; y
g(Xi!l) f(X,,t) dt (vuo(Xi,)) hg(X;,
((xi.t) (xxxx))
odB
X yo
20 uo (00
o))=
(4.6)
Fo t) + Go( oB.
The augmented SDE is discretized as follows using the stochastic Heun scheme:
2 2.00 FF 2 (n) /FG 22.0 2..0 Wnn w. N(O,
(4.7)
2.0
7.
T
Fo 2::0 tt)) Fo 2+1 (:1 (nt1))
(Z
VT (Gol2n"tth) tw.) :) Go GZ +
wn.
2
with 2,000 (xo: (x 0)). This gives rise to the k-step Heun-BSDE loss defined as:
THHnmn, (0) = 3x
E W
Σ n=0
(no( (no(Xin=)))) (m+1)k) uo &i (ink)(n)) bnnk ) (V (n+1)R - na )² (4.8)
with 3x N For the step t case. use the shorthand LHeun, (0). We now show that
the one-step one-step EM-BSDE loss LEM,- (cf. Section 4.1), To this, analogous to (4.1), we define the point-wise Heun loss resolution for a fixed (æ,t)
loss LHeun. (0) avoids the undesirable bias term which appears for the
r (0, E. (up (wo (Ze.r, t 1o. )
5

---

=lt = 2, +Fo( VTGGoZ2, t)w, - (z,
vT (Gol + Go Col t T))
2 21 2 (F6((2) Fo Fo
0
2
noting that 26.r (Seur, Be Similar to L (0). we have the following identity Leunn, (0) 1 N En-o E xx [HHeum, (0, X;, x [. Our next result illustrates that the point-wise Heun loss avoids the issues identified with the point point-wise EM loss in Lemma 4.1.
Lemma 4.3. Suppose that f, 9, and he are all in Che,, and is in C3.1, We have that
=2. Heun, r (0,x, x, t)) (HQ ue (e]( ((())))))())
0(7))))
(4.9)
where the O()) hides factors depending d and the Hölder norms of f. g. hg. and ug.
Furthermore, analogously to Theorem 4.2, we can utilize Lemma 4.3 in conjunction with the order 1/2 strong convergence of stochastic Heun (ef. Appendix F) to the Stratonovich solution to show the following relationship for the full loss I
Theorem 4.4. Suppose that J. 9, and ha are all in C 6 C and < 1. We have that
L (0) = T L." E [X (()))))
(4.10)
where the O() hides factors depending d, T, and the Hölder norms of f. g, he, and ug.
Therefore, unlike the situation with EM integration in (4.3), any additional bias terms only enter through a O(r term which is of higher order than the leading PDE residue term (cf. (4.10). In Figure 1b, we show the plot of [Heun,-(f)) on the same HJB PDE problem in Figure 1a, and show that the bias issue in LEM. (#) is now resolved. The proofs of both Lemma 4.3 and Theorem 4.4 are given in Appendix D.3.
5 Trade-offs for Multi-Step BSDE Losses
In Section 4, we conducted a thorough analysis of the one-step self-consistency losses LEM. (0) and LHeua,+ -0). We now consider the other extreme: the full-horizon (k N) losses LEMN - (0) and Leeeax.7 (0). The intermediate multi-step regime [5]. where 1 < k < N, serves as an extension to the cited method and studied experimentally in Section 6. Due to space constraints, we defer the precise theorem statements arising from our analysis, in addition to the proofs, to Appendix E.
BSDE loss and Euler-Maruyamo discretization. We start with the H-horizon BSDE loss (3.5). Using Jensen's 's inequality, we show (Proposition E.3) the relationship LBSDE, & L 1/3 Thus, at the SDE level, the benefits of using the full-horizon loss LBSDETT over the --horizon loss LBSDEE, (f) are not clear, given that (a) the full-horizon loss dominated by the latter T-horizon loss (up to an order 1/2 term), and (b) LESDE, (0) does indeed vanish for an optimal 6..
The situation becomes more complex when factoring in EM discretization. Using the order 1/2 strong convergence of FMM, we show (Proposition E.4) that (EM LBSDETT()) O On the other hand, from Theorem 4.2 we also show (Proposition E.5) that LEM,- LISDDE, Bias(0) (/2/22) with the bias term (0 == 2r So" E Etr(((((X,))))) (X,,)) Vuug (x,, 2) dd not vanishing
O. Hence, the loss LEM LEMN.~(O) (0) presents an advantage over LEM.- (0) for sufficiently small discretization sizes in terms of bias. However, the inequality LBSDE,T (0) < LBSDE,-
(271/2) still holds, meaning that while LEMM.r(8 (f) does not suffer from the bias issues identified in LEM,- (f), the trade-off is that the loss LBSDE,T(()) it approximates without bias is nearly dominated by another loss LESDE,- (f): this is precisely the loss that (EM,-(() attempts to approximate, but it does so in way that introduces losses for = nor ie N provides a completely satisfactory solution. In Section 6.2, we illustrate these issues empirically. Furthermore, in light of this analysis, we can interpret the interpolating loss of [5] as attempting to resolve this trade-off by finding the best intermediate multi-step & E {I,...,N)
irreducible bias term Bias(f). Thus, neither of the EM-BSDE
Stratonovich BSDE and Heun discretization. In the setting of the Stratonovich BSDE and the Heun- BSDE loss, we first show (Proposition E.8) that Le-BSDE. (0) & LS-BSDE,r (0) O(71/2) holds at

---

the SDE level, analogous to the relationship between LESDE,T (0) and LBSDE.> (0)) Next, we use the order 1/2 strong convergence Heun to show (Proposition E.9) that _e (0) = LS..E,
0(71/2); again analogous to the relationship between LEMy.r (0) and LBSDE,T(O), However, unlike the one-step EM case, using Theorem 4.4 show (Proposition E.1 that LHeum,=
L (0) + ((71/2), from which we conclude L (0) & Heuu, (0) + 0(71/2). Thus- unlike the EM setting-the the relationship between the full-horizon and one-step case at both the level the same, suggesting questionable benefits of LeeuNN (0) over (0). In Section 6.2, we show that this conclusion is indeed reflected in practice.
SDE level and
Experiments
In this section, we compare the proposed Heun-based BSDE method against both standard PINNs, variant of PINNs which uses the forward SDE to sample collocation points, and standard EM-based BSDE solvers on various high-dimensional PDE problems. Specifically, we compare the methods:
(a) PINNs:
The standard PINNs loss LPINNs (0) from (3.2) is minimized. Since we consider un- bounded domains, the collocation measure μ over (Ir, is chosen by fitting a normal distribution over the spatial dimensions of the forward SDE trajectories prior to training.
(b) FS-PINNS: The standard PINNs loss (3.2) is again minimized, where the measure μ over space-time is chosen by directly sampling trajectories from the forward SDE (3.3).
(C) EM-BSDE: The self-consistency loss discretized with the standard Euler-Maruyama (EM) scheme, i.e., LEM., , (0) from (3.7) as described in Section 3.
(d) EM-BSDE (NR): A variant of EM-BSDE where use the BSDE to propagate Y, instead of setting it directly to us (X,. t) [cf. 4, 5]. We refer to this variant as the no-reser (NR) variant. Specifically, the backward SDEs in (3.6) is integrated as (starting from G0 = us (xo, 0)):
y: n-1 8 + -ho Xn, In, yo Vuo &x,,,,,) bn)) /FV us (xn, (n) g( (xn, Ln Wn,
(6.1)
with Wn ~ N(O, Ia). The loss LEM..-(f) remains the same except replacing (3.6) with (6.1).
(e) Heun- BSDD (Ours): The self-consistency loss discretized with stochastic Heun integration, i.e., L (0) from (4.8) as described in Section 4.2.
We evaluate these methods three PDE benchmark problems: (i) a Hamilton-Incobi-Bellme (HJB) equation [4], (ii)
Shock--Sholes--ee -Breenblaa (BSB) equation [4], and (iii) a fully-coupled FBSDE from Bender & Zhang (BZ) [36]: the PDEs are detailed in Appendix G.1. addition, we evaluate the methods on an optimal control pendulum swing-up problem to demonstrate application to non-linear control problem (see Appendix G.6). To evaluate model performance, the analytical solution (available for all PDEs under consideration) is compared with the model output along 5 forward SDE trajectories, using the relative L2 error (RL2) metric:
N (Uret (iee((X,,, (x,44)) Uprod (X,, 44)
RL2 ==
(6.2)
ΣNn de (X,,
Unless otherwise noted, we set T 1 and N 50 (i.e., == 0.02). Model architectures and training details are described in Appendix G. Additionally, the code reproduce our experiments is available at: https: : / /gi tiub. com/ sungje- e park//eunbaeeeee heunbade. For what follows, we report two main sets of results on (i) one-step self-consistency losses (Section 6. and (ii) multi-step self-consistency losses (Section 6.2).
Efficient sub-sampling BSDE implementation. In experimental results, we consider both full FSDE rollout algorithm, in addition to a batched, sub-sampled FSDE rollout variation of FS-PINNs, EM-BSDE, and Heun-BSDE, which we find performs similarly the original algorithm while providing significant speed improvements. We sketch the details of the sub-sampled FSDE variation Algorithm I; further details on the two algorithms can be found in Appendix G.4.
6.1 One-Step Self-Consistency Losses
For our first set of results, we compare the PINNs baselines and the one-step (k: 1) EM-BSDE baseline with our one-step Heun method. We solve each PDE instantiated with a state space
7

---

Algorithm Batched, Sub-sampling BSDE Algorithm (Simplified)
Input: Neural network uo t), parameters 6, terminal function Ø. time step At, trajectory length N,
evaluation batch size B.
1: Sample initial state: (=:[0], [[]]] (zo, 0), with zo -
2: Sample Brownian noise: <[0: . 1] N(0, Ia)
3: Evaluate network at initial state: (u, u.)
4: /* Forward SDE rollout */
5: for === N do
6: /* Use either EM Heun integration */
7: Propagate forward state (with (u, Uz) coupled): =[]]]] Ar
8: Propagate time: t[i + 1]] 1]= t[$] + At
9: Evaluate network new state: (u, ux) (mo 1]. 1)), (.)),,atpl V= 11,t(++ 1)))
10: end for
Stop gradient: ={0: .] SG(=[O: N])
12: Random sub-sampling: (Ir, Fi-1, Ii-1, 4.11)
13: /+ Use either EM or Heun integration +/
14: Propagate backward SDE batched points: Y/+1 14 Ay
Use PINNs
instead for FS-PINNS */
Compute self-consistency loss: Let Ei-D (#++1 y4t1)
Compute terminal loss: Cg (u + di²² llue V.
Cases
PINNs
FS-PINNs EM-BSDE (NR)
EM-BSDE
Heun-BSDE
Full Algorithm (ef. Algorithm Results
1OOD HJB
.1260 + .0107
.0737 + .0110
:214 .0452
3626 = .0113
.0493 = .0109
100D BSB
1.5066 + .2349
.0497 + .0031
.1855 .0078
.3735 .0470
.0535 .0113
100D BZ.
3. 3.1259 + 11807
3.5619 = .2716
10D BZ.
3.8566 .0310
.031 + .0011
.1309 .0311
.1903 .0066
.0228 .0016
Batched Algorithm (cf. Algorithm 3)
Results
100D HJB
1362 # .0276
.1888 + .0774
.5214 .0452
.3811 .0884
.0573 = .0106
100D BSB
3.0488 + 1.5625
.0851 + .0027
.1855 .0078
.3668 = .0244
.0472 = .0076
100D BZ
5.4502 .1351
5.7330 .2422
1.7973 .1108
10D BZ
3.8495 + .562
.0270 + .007
.1309 .0311
.1933 = .0022
.0236 .0031
Table 1: Summary of RL2 errors averaged over three different initialization random seeds,
standard satisfactory solution are denoted with -. The first set of results
deviation. Settings that failed to converge
correspond to the full algorithm (see Appendix G.4 detailed description), whereas the second set of results correspond to the batched algorithm (cf. Algorithm 1).
of 100 dimensions using three different initialization seeds for training. The results are reported in
Reference
EM-BSDE Heun-BSDE
PINNs
FS-PINNS EM-BSDE
4.8
FS PINNs
10-2
BSDE
EM-B5DE (NR)
uik. 8
4.7
Error
4.6
10-3
10
0.0
0.2
0.6
0.B
0.0
0.2
0.4
0.6
1.0
Time, tt
Time.
(a) Plot of learned solution different models on the 100D HJB problem.
(b) the RL2 errors across time [0, 1] the
100D HJB case.
Figure 2: A plot of the 100D HIB reference and learned solutions for each model and the associated RL2 errors.

---

Table 1, which shows that for nearly all the cases, the Heun-BSDD method outperforms EM-BSDE methods (i.e., lower RL2 error) as predicted by our analysis in Section 4. Furthermore, Figure illustrates the performance of all methods across time / ε [0, 1] on the 100D HJB case, which also highlights the low RI.2 error of the Beun-BBDDE method. The one exception to the trend is the 100D) BZ. case, where all methods failed to converge to a high-quality solution. We hypothesize due to the high dimensionality of the problem involving fully-coupled SDEs, the optimization landscape for all methods is too complex to recover high-fidelity solutions. To evaluate this hypothesis, we further reduce the dimensionality of the BZ problem to IOD, which restores the relative performance of all methods (cf. Table 1, last row). Another key observation from Table 1 is that FS FS-PINNs and Heun-BSI BSDE perform similarly across all cases, showing that parity between the BSDE and PINNs is restored through Heun integration. Finally, we note that the performance of PINNs is quite poor in comparison to ES-PINNs, which illustrates the relative impact of the sampling distribution j for PINNs methods (ef. (37-391).
To show that the gap between EM and Heun performance cannot be closed with finer discretization meshes, we re-run the 10D BSB example at varying discretization sizes. The results are reported in Figure 3, which show that the EM-BSDE methods only experience minimal improvement with smaller discretization size compared with the Heun-BSDE method. Figure 3 corroborates the findings in Section 4, which shows that the one-step EM-BSDE loss contains a bias term of the same order as the residual error which is not present in the one-step BOuu-BSDDE loss.
Computational considerations. Although Heun-BSDE outperforms EM-BSDE, it comes at a computational cost. In Table 2, we report the average training time for each method on single NVIDIA A100 GPU; on average Heun- BSDE is approximately 6x slower than the EM-BSDE method for the batched algorithm. There are two ma- jor factors to this overhead. First, for the specific ellip- tic/parabolic PDEs we consider in (). EM-BSDE does not require the computationally expensive Hessian com- putation v²u(zz f). which PINNs. FS-PINNs, and Hcun-BSDE all do require. However, this does not necessarily hold true for all PDEs (c.g., Appendix G.6). Second, the Hcun integration requires approximately double the compute of the EM integration this is clearly reflected in the overhead between FS-PINNS and BSDE.
Method
Full
Batched
PINNs
Ix
1x
FS-PINNS
2.64x
1.14x
EM-BSDE
2.83x
0.34x
36.37x
2.03x
Table 2: table average training time overhead relative to PINNs for both the full and batched algorithm runs.
In Figure 4, we show the cuntime-normalian RL2 performance demonstrating that while EM-BSDE shows strong convergence at first, its performance does not improve with more compute. Conversely, both FS-PINNs and Heun-BSDE achieve similar RL2 performance at equal runtimes.
6.2 Multi-Step Self-Consistence Losses
We next consider multi-step celf-consis BSDE losses [5] in order to evaluate the mathematical analysis conducted in Section 5, Specifically, we evaluate both the multi-step LEMA, (0) (cf. (3.7))
PINNS
10: 7
10-1
FS PIN
EM-BSDE
Heun-BSDE
EM-85DE EM-BSDE (NR)
RL2
RL2
10-2_
Heun-BSDC
10
10
0.04D 0.020
0,010
0.005
10
20
50
60
At
Runtime (min)
Figure 3: RL2 performance for 10D BSB at discretiza- tion step-sizes - for N (25, 50, 100, 200).
Figure 4: A plot of the RL2 performance versus run- time for the 100D HIB problem.

---

and (Heume.7(0) (cf. (4.8)) for varying values of skip-length k. For multi-step losses, we also need determine where both EM-BSDE and EM-BSDE (NR) will "reset" the value of Y, ug (X,,,). Note that there are many degrees of freedom here in the multi-step formulation, so we simply pick one choice as a representative choice. For EM-BSDE, to integrate between ln and ((n+1)k- On the other hand, for EM-BSDE (NR), we directly use the value of y from (6.1). We also vary EM-BSDE and EM-BSDE (NR) with N, the number of discretization steps for the interval [O, 1]. varying between E (50, 200} as well. We conduct this experiment on the 10D BSB setting, with the results reported in Figure 5.
set nkk us(( (Xnne and use (6.1)
10"
10a
EM-BSDE
EM-BSDE (N-200)
EM (NR)
BSDE (NR, N=200)
10
EM-BSDE (N-200)
-0:1
EM BSDE (NR, N=200)
Heun BSDE
RU2
10-2
10-2
10-3
10-3
PINNS Baseline
FS-NINGs Baséline
Heun BSDE Bascline
10-4
10-4
10
20
30
40
50
0
50
100
150
200
Skip-Length,
Skip-Length,
(a) Plot of RL2 performance as a function of the skip- length k. with the number of discretization steps also varying in N (50, 200).
(b) Plot of RL2 performance as a function of the skip- length & for EM-BSDE and EM-BSDE (NR), with number of discrctization steps set to N 200.
Figure 5: A plot of RL2 performance of each model on the 10D BSB case at various skip lengths.
Figure 5 shows that while the Heun-BS -SDD performance decreases as the skip-length l increases, the performance of both EM-BSDE methods initially improves with skip-length before then degrading, demonstrating the trade-off between minimizing the bias term and decreasing quality of the self- consistency loss identified in Section 5. Furthermore, this trade-off is present for EM-BSDE at both N = 50 and N 200, illustrating again that these issues for EM are not mitigated with finer discretization step-sizes. Finally, although the EM-BSDE method overall improves its performance by tuning the skip-length & (consistent with the findings from [5]), the Heun-BBDDE model at k still outperforms the best multi-step EM-BSDE model by a significant margin.
7 Conclusion and Future Work
We conducted a systematic study of discretization strategies for BSDE- basee loss formulations used solve digeenimensio PDEs. By comparing the commonly used Maruuama scheme with stochastic Heun integration, we demonstrated that the choice of discretization can significantly impact the accuracy of BSDE-based methods. Our theoretical analysis showed that EM discretization introduces non-trivial bias to the single-step self-consisteng the step-size decreases. On the other hand, we show that this bias issue is not present when utilizing Heun discretization. Finally, the empirical results confirmed that the Heun scheme consistently outperforms EM in solution accuracy and performs competitively with PINNs.
BSDE loss which does not vanish as
Our work underscores the importance of stochastic integrator choice in BSDE solvers and suggest that higher-order schemes though more computationally intensive can offer substantial gains in performance. In future works, we aim to reduce Beun-BSDE's computational costs through methods such as Hutchinson trace estimation [40], reversible Heun [41], and adaptive time-stepping. Furthermore, while this current work focuses on understanding and restoring performance parity to solve
between BSDE and PINNs methods, future work will utilize the advantages of problems such as high-dimensional stochastic control problems in model-free settings (cf. [6]).
10

---

References
[1] Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Physics informed deep learn- ing (part i): Data-driven solutions of nonlinear partial differential equations. arXiv preprint 17X1::7111.106611 2017.
[2] Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Physics informed deep learning (part ii): Data-driven discovery of nonlinear partial differential equations. arXiv preprint 2017.
[3] Weinan E, Jiequn Han, and Arnulf Jentzen. Deep learning- numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations. Communications in Mathematics and Statistics, 5(4))349-380) -8000 2017.
[4] Maziar Raissi. Forkwrd-backward
backward stochastic neural networks: deep learning of high- dimensional partial differential equations. Peter Carr Gedenkschrift: Research Advances in Mathematical Finance, pages 637-655. World Scientific, 2024.
[5] Nikolas Nüsken and Lorenz Richter. Interpolating between bsdes and pinns: Deep learning for elliptic and parabolic boundary value problems. Journal of Machine Learning, 2(1):31-64, 2023.
[6] Yutian Wang and Yuan-Hua Ni. Deep bsde-ml learning and its application to model-free optimal control, 2022.
[7] Sifan Wang, Shyam Sankaran, Hanwen Wang, and Paris Perdikaris. An expert's guide to training physics-ifformed neural networks. arXiv preprint arXiv:2308.08455 2023.
[8] Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Physics-informed neural networks: deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2019.
[9] Weinan E and Bing Yu. The deep ritz method: deep learning- based numerical algorithm for solving variational problems. Communications Mathematics and Statistics, 6(1):1-12, 2018.
[10] Justin Sirignano and Konstantinos Spiliopoulos. Dgm: A deep learning algorithm for solving partial differential equations. Journal of Computational Physics, 375: 2018.
[11] Aditi Krishnapriyan, Amir Gholami, Shandian Zhe, Robert Kirby, and Michael W Mahoney. Characterizing possible failure modes in physics-informed informed neural networks. In Advances in Neural Information Processing Systems, volume 34, pages 26548-26560. Curran Associates, Inc., 2021.
[12] Pratik Rathore, Weimu Lei, Zachary Frangella, Lu Lu, and Madeleine Udell. Challenges in training pinns: a loss landscape perspective. In Proceedings of the 41st International Conference on Machine Learning, ICML'24. 24. JMLR.org, 2024.
[13] Sifan Wang, Xinling Yu, and Paris Perdikaris. When and why pinns fail train: A neural tangent kernel perspective. Journal of Computational Physics, 449: 10768, 2022.
[14] Tao Wang, Bo Zhao, Sicun Gao, and Rose Yu. Understanding the difficulty of solving Cauchy problems with PINNs. In Proceedings of the 6th Annual Learning for Dynamies & Control Conference, volume 242 of Proceedings of Machine Learning Research, pages 453-465. PMLR, 15-17 Jul 2024.
[15] Yi-Yueh Chuang and Lorena Barba. Predictive limitations of physics-informed neural networks in vortex shedding. arXiv preprint arXiv:2306.00220, 2023.
[16] Chenxi Wu, Min Zhu, Qinyang Tan, Yadhu Kartha, and Lu Lu. A comprehensive study of non- adaptive and residual-base adaptive sampling for physics-ifformed neural networks. Computer Methods Applied Mechanics and Engineering, 403:11567, 2023.
[17] Haixu Wu, Huakun Luo, Yuezhou Ma, Jianmin Wang, and Mingsheng Long. Ropinn: Region optimized informed neural networks. arXiv preprint arXiv:2405. 245..24005.4369 2024.
11

---

[18] Chuwei Wang. Shanda Li, Di He, and Liwei Wang. Is L2 physics informed loss always suitable for training physics informed neural network? In Advances in Neural Information Processing Systems, volume 35, pages
Curran Associates, Inc., 2022.
[19] Zhiwei Gao, Liang Yan, and Tao Zhou. Failure-informed adaptive sampling for pinns. SIAM Journal on Scientific Computing, 45(4) 2023.
[20] Rambod Mojgani, Maciej Balajewicz, and Pedram Hassanzadeh. Lagrangian pinns: A causality- conforming solution to failure modes of physics- informed neural networks. arXiv preprint arXiv:2205.02900 2022.
[21] Jiequn Han, Arnulf Jentzen, and Weinan E. Solving high-dimensional partial differential equations using deep learning, Proceedings of the National Academy of Sciences, 15((4))85005- 8510, 2018.
[22] Côme Huré, Huyên Pham, and Xavier Warin. Deep backward schemes for high- high-dimensional nonlinear pdes. Mathematics of Computation, 89(324):15471157, 2020.
[23] Christian Beck, Weinan E, and Arnulf Jentzen. Machine learning approximation algorithms for high-dimensional fully nonlinear partial differential equations and second-order backward stochastic differential equations. Journal of Nonlinear Science, 29: 19:136916,,, 2019.
[24] Akihiko Takahashi, Yoshifumi Tsuchida, and Toshihiro Yamada. A new efficient approximation scheme for solving high-dimensional semilinear pdes: Control variate method for deep bsde solver. Journal of Computational Physics, 454: 10956, 2022.
[25] Kristoffer Andersson, Adam Andersson, and C. W. Oosterlee. Convergence of a robust deep fbsde method for stochastic control. SIAM Journal on Scientific Computing, 45(():A226-A255, 2023.
[26] Zebang Shen, Zhenfu Wang, Satyen Kale, Alejandro Ribeiro, Amin Karbasi, and Hamed Hassani. Self-consistency of the fokker planck equation. In Proceedings of Thirty Fifth Conference on Learning Theory, volume 178 of Proceedings of Machine Learning Research, pages $17-841. PMLR, 02-05 Jul 2022.
[27] Nicholas M Boffi and Eric Hunden-Eijnden Probability flow solution of the fokker fokker-planck equation. Machine Learning: Science and Technology, 4(3):035012, 2023.
[28] Lingxiao Li, Samuel Hurault, and Justin M Solomon. Self-consistent velocity matching of probability flows. In Advances in Neural Information Processing Systems, volume 36, pages 77038-77057 Curran Associates, Inc., 2023.
[29] Fran-François Chassagneux, Junchao Chen, and Noufel Frikha. Deep runge- kutta schemes for bsdes. arXiv preprint arXiv:2212.14372. 2022.
[30] LEE Lagaris, A. Likas, and D.I. Fotiadis. Artificial neural networks for solving ordinary and partial differential equations. IEEE Transactions on Neural Networks, 955):887-1000, 1998.
[31] Aditya Singh, Zeyuan Feng, and Somil Bansal. Exact imposition of safety boundary conditions in neural reachable tubes. arXiv preprint arXiv:2404.00814 2024.
[32] Quentin Chan-Wai-Nee Wai-Nam, Joseph Mikael, and Xavier Warin. Machine learning for semi linear pdes. Journal of Scientific Computing, 79(3):1667-1712, 2019.
[33] Lorene Kapllani and Long Teng. Deep learning algorithms for solving high-dimensiooal nonlinear backward stochastic differential equations. Discrete and Continuous Dynamical Systems Series B (DCDS-B), 29(4): 2024.
[34] Peter F. Kloeden and Eckhard Platen. Numerical Solution of Stochastic Differential Equations. Springer, 1992.
[35] W. Rümelin. Numerical treatment of stochastic differential equations. SIAM Journal on Numerical Analysis, 1943)6604-113 1982.
12

---

[36] Christian Bender and Jianfeng Zhang. Time discretization and markovian iteration for coupled fbsdes. The Annals of Applied Probability, 18(1):143-177, 2008.
[37] Mohammad Amin Nabian, Rini Jasmine Gladstone, and Hadi Meidani. Efficient training of physics-informed neural networks via importance sampling. Computer-Aided Civil and Infrastructure Engineering. 36(8):962-97)) 2021.
[38] Arka Daw. Jie Bu, Sifan Wang. Paris Perdikaris, and Anuj Karpatne. Mitigating propagation failures in thysics-informee neural networks using stain-resample-relesss (R3) sampling. In Proceedings of the 40th International Conference on Machine Learning. volume 202 of Proceedings of Machine Learning Research, pages 7264-7302. PMLR, 23-29 Jul 2023.
[39] Zhengqi Zhang, Jing Li, and Bin Liu. Annealed adaptive importance sampling method in pinns for solving high dimensional partial differential equations. Journal of Computational Physics, 521:11 13561, 2025.
[40] Zheyuan Hu, Zekun Shi, George Em Karniadakis, and Kenji Kawaguchi. Hutchinson trace estimation for high-dimensional and high-order physies-informed neural networks. Computer Methods in Applied Mechanics and Engineering, 424: 1 16883, 2024.
[41] Patrick Kidger, James Foster, Xuechen (Chen) Li, and Terry Lyons Efficient and accurate gradients for neural sdes. In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wort- man Vaughan, editors, Advances in Neural Information Processing Systems, volume 34, pages 18747-18761. Curran Associates, Inc., 2021.
[42] Riu Naito and Toshihiro Yamada. An acceleration scheme for deep learning-based bsde solver using weak expansions. International Journal of Financial Engineering. 07(02):2050012, 2020.
[43] Huyen Pham, Xavier Warin, and Maximilien Germain. Neural networks-basee backward scheme for fully nonlinear pdes, 2021.
[44] Christian Beck, Sebastian Becker, Patrick Cheridito, Arnulf Jentzen, and Ariel Neufeld. Deep splitting method for parabolic pdes. SIAM Journal on Scientific Computing, 43(5):A3135 A3154, 2021.
[45] Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces with applications to pdes, Journal of Machine Learning Research, 24(89):1 -9,, 2023.
[46] Lorenz Richter, Leon Sallandt, and Nikolas Nüsken. From continuous-ting formulations to discretization schemes: tensor trains and robust regression for bsdes and parabolic pdes. arXiv preprint arXivv:307.15599 2023.
[47] Jiequn Han and Jihao Long. Convergence of the deep bsde method for coupled fbsdes. Proba- biliry, Uncertainty and Quantitative Risk, 5(1), 2020.
[48] Terry J. Lyons. Differential equations driven by rough signals. Revista Matemática Iberoameri- cana, 14(2):215-310, 1998.
[49] Jan R. Magnus. The moments of products of quadratic forms in normal variables. Statistica Neerlandica, 32(4):201-210, 1978.
[50] Tuomas Hytönen, Jan Van Neerven, Mark Veraar, and Lutz Weis. Analysis in Banach Spaces: Volume I: Martingales and Littlewood-Paie Paley Theory. Springer, 2016.
[51] Dean S. Clark. Short proof of a discrete gronwall inequality. Discrete Applied Mathematics, 16(3):279-281, 1987.
[52] Prajit Ramachandran, Barret Zoph, and Quoe V Le. Searching for activation functions. arXiv 2017.
preprint arxiv:
[53] Matthew Tancik, Pratul P. Srinivasan, Ben Mildenhall, Sara Fridovich- Keil, Nithin Raghavan, Utkarsh Singhal, Ravi Ramamoorthi, Jonathan Barron, and Ren Ng. Fourier features let networks learn high frequency functions in low dimensional domains. In Proceedings of the 34th International Conference on Neural Information Processing Systems, NIPS '20 Red Hook, NY, USA, 2020. Curran Associates Inc.
13

---

[54] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Mi -Mine,, and Qiao Zhang. JAX: composable transformations of Python+NummP programs, 2018.
[55] Haoyu Han and Heng Yang. On the nonsmooth geometry and neural approximation of the pendulum swing-up. up. arXiv preprint arXiv:2312.17677,
optimal value function 2024.
14

---

NeurIPS Paper Checklist
1. Claims
Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Question:
Answer: [Yes]
Justification: We make two main claims in the abstract and introduction: (i) the performance gap between existing BSDE and PINNS methods explained by an irreducible bias intro- duced through the use of EM integration for BSDE, and (ii) reformulating Stratonovich BSDE loss and integrating it via Heun removes the irreducible bias and restores parity between the two methods. Both claims (i) and (ii) are well supported theoretically (Section and Section 5) and experimentally (Section 6). Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results
be expected to generalize other settings.
is fine to include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: Please see Appendix A for a discussion of the limitations of our work. Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed not penalize honesty concerning limitations.
3. Theory assumptions and proofs
15

---

Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
Justification: Both Lemma 4.1 and Lemma 4.3 have stated their full set of assumptions.
and have accompanying proofs in Appendix D.2 and Appendix D.3. Furthermore, for the theoretical claims made in Section the corresponding derivatives are given in Appendix
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: Both Section 6 and Appendix G contain the necessary details to reproduce the experiments conducted in this paper, including architectural choices, learning rate schedules, batch sizes, benchmark PDE formulas and analytical solutions, and compute platforms used. We have also included our code github repository with a README. md file.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily the architecture clearly and fully.
model architecture, the paper should describe
(c) the contribution a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
16

---

(d)
We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, described in supplemental material?
Answer: [Yes]
Justification: For the submission, we have included an anonymized version of the code in the supplementary material, including a README. md file. For the camera ready version of the paper, we have included a link to our github repository bectps://githi com/ sungje e (ergjjeparr/heee park/heunbsde). Our code on github is released under the Apache 2.0 License.
Guidelines:
The answer NA means that paper does not include experiments requiring code. Please see the NeurlPS code and data submission guidelines (https://nips.. cc/ helliiccgudes///// httgg///oooeeeuummaaaaa
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to reproduce the results. See the NeurlPS code and data submission guidelines (https: /tips.cc/puu wttp:cccccpbccccccccccccccccclct wwthc//gugee//gee heweeidee//ooeemmieeeeeee soiorbiiissioooooooiiooo for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing much information possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: At the beginning of Section 6, we include detailed description of the all methods being evaluated, including the loss functions and how the collocation//raaecttory data is generated. The benchmark PDEs are described in detail in Appendix G.1, which includes the analytical solutions we use to evaluate model quality. We define within the use to assess model fidelity, and describe the
main text the error metric RL2 (6.2) which number of test trajectories we collect for computing RL2, Finally, as noted previous, details about the training pipeline are described in Appendix G.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental
material.
7. Experiment statistical significance
17

---

Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: All models are trained over different training initialization seeds. In Table 1, Table 3, and Table 2, we report the average RI22/rntime + one standard deviation. For the figures, we only report the average in order to avoid cluttering the diagrams, since the standard deviation is modest relative to the average (cf. Table 1).
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: In Appendix G, we describe our compute environment (single node NVIDIA A100 GPU in our internal cluster).
Guidelines:
The answer NA means that the paper does include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported the paper (e.g., preliminary or failed experiments that NeurlPS Code of Ethics https: hetp::/nneui.ppeeee wwwpurips.cccppuulo cwtppblic/eeaeeeeeeeeecec http://iwc/hhiiisiieeeeeeeeeeeee
didn't 'tttt make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the
Answer: [Yes]
Justification:
We confirm that our work conforms with the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
18

---

If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: Our work falls into the category of foundational research; we do not believe it poses any specific and/or unique potential harms.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative socictal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deploymeett of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied partieular applications, let alone deployments. However, if there is direct path any negative applications, the authors should point it out. For example, it is legitimate to point that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, scraped datasets)?
Question:
Answer: [NA]
Justification:
The assets resulting from our work do not pose high risk for misuse.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
19

---

Answer: [NA]
Justification: Our code does not use any existing assets beyond standard open-source machine learning frameworks, which we give proper attribute to.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset. For scraped data from a particular source (e.g.. website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
this information is not available online, the authors are encouraged to reach out to the asset's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification:
We provide documentation via README, md file for the code we release.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/mooel as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Question:
Answer: [NA]
Justification: Our paper does not involve crowdsourcing research with human subjects.
Guidelines:
The answer NA means that the paper does involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
20

---

Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: Our paper does not involve crowdsoureing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with
human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods this research? Note that the LLM used only for writing, editing, or formatting purposes and does not impact the core methodology. scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA]
Justification: LLMs do not constitute important, original, or non-standard component of our work.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non non-standard components.
Please refer to our LLM policy (https: http///eerpp../ee cc/Confferences/20222 for what should or should not be described.
21

---

A Limitations
We provide concise, bulleted list of the limitations present in our work.
(a) Computational overhead of BBDD: As discussed in Section 6.1 (ef. Table 2), the computa- tional overhead of the Heun-BSDE method compared with the existing EM-BSDE methods is non-trivial. Furthermore, we also found the Heun integrator to be more susceptible to floating point imprecision (cf. Table 3), and hence we run our main experiments in float64, which further adds to the computation time. We mitigate these issues using batched computation and random sub-sampling (see Appendix G.4) which helped significantly reduce computation time, but there still remains a computational penalty. While we believe more sophisticated techniques (c.g., randomized trace estimation and more numerically stable Heun integrators discussed in Section 6.1) can help to reduce the computational overhead, we have not yet verified this hypothesis experimentally our current work.
(b) Performance relative to PINNs: Our results in Section 6.1 show that the proposed BE method restores parity with the FS-PINNs method in terms of the RL2 error. While this is significant improvement over EM-BSDE, Eleun-BSDD does not yet provide significant (e.g.. orders of magnitude) performance improvement over the best PINNs method (cf. Table 1); further work is needed to determine whether or not such an improvement is possible. Hence, the current advantage of Heun- BeunE over PINNs is with its model-free capabilities (ef. Appendix C).
(c) Limitations of theoretical analysis: While our theoretical analysis in Section 4 and Section is fairly predictive of practice (ef. Section 6), our analysis is not without its own set of limitations. One limitation is that we only analyze the two extremes of horizon length: the one-step case (Section 4) and the full-horizon N-step case (Section 5) the intermediate regimes (i.e., skip- lengths & satisfying < < N) are studied empirically (cf. Section 6.2). Another limitation is that we do not consider fully-coupled FS-BSDE: in our analysis (e.g., the Bender & Zhang (BZ) PDE, cf. Appendix G. where the forwards SDE (3.3) is allowed to depend on Y,.
B Further Discussion on Related Works
While our work focuses the impact that the widely used EM integration scheme has on the performance of BSDEEbassd solvers, several complementary enhancements to the BSDE loss have been proposed in literature [6, 24, 42-44]. Many of these improvements with our proposed Heun-B BSDE framework. For example. the -eeun method could be adapted to utilize a control variate [24], applied in operator splitting settings [44], and extended to fully non-linear PDEs [43], enabling direct comparison against their EM-BSDE baseline.
naturally compatible
Furthermore, in addition to PINNs and EM--ased BSDE methods discussed in the paper, there are various other deep-learning methods for solving PDEs such as Deep Ritz [9], Neural Opera- tors [45], and tensor trains [46], in addition to various theoretical analyses developed for Itô- dô-based BSDE formulations [47]. We leave extending these approaches and analyses to Stratonovich-base formulations as future research directions.
C Details Regarding Model-Free BSDE Formulation
Suppose that the drift term fax 1) from the forwards SDE (3.3) is unknown. Instead, suppose that our the computational model takes as input realization of ( [o and returns FSDE trajectories (Xx) or more practically a sub-sampled trajectory { Under this computational model, PINNs methods cannot be used to solve the PDE (3.1) if the other terms g. h, and are all known-since the residual term R[u) cannot be computed without knowledge of the drift term f. However, in this settings, BSDE methods including the proposed Heun-BSDE method can still be used. This is because BSDE losses only require access to the FSDE trajectory (X,) and the Brownian motion (B,) used generate it (cf. Equation (3.5)).
A specific example where the proposed computational model is realistic comes from model-free optimal control, and first described in [6]. Consider the following deterministic continuous-time control-affine fully-actuated system:
X, f(x,t) ((X,, X2, ER, rank()(z,t)) (z,t).
(C.I)
22

---

We assume that we do not know the drift term f(,,t,, but we are able to select control inputs Uz and obtain the resulting trajectories (X,): this setting is often called the model-free setting in optimal control and reinforcement learning. In this framework, by setting the input U, to be a nominal input U. injected with excitation noise injected, i.e., U, U. + "Noise,", ", we can select the realization (B,) of Brownian noise and observe the trajectories of forward SDE of the forms:
dx, Lf(X,,t)+ ((X,, e)U, + t)),,) dBB,
(C.2)
where the o indicates the SDE is to be interpreted in terms either Itô or Stratonovich, depending on the context. To rigorously define "Noise," to establish a connection between (C.1) and (C.2) is technical, requiring the use of e.g. rough path theory [48]. We will take more practical approach inspired from [6] and observe how injecting Gaussian noise into discretizations of (C.1) induces stochastic discretizations of (C.2).
Concretely, we proceed as follows. We work with constant step-size integrators, and define integration times I. kr, N, for step-size 7 Given a nominal control input U. which an open-loop (i.e., only time-dependent) signal,' we form our control input U, by injecting Gaussian noise as follows:
U. = +
Wk ~ N(O, Ia).
(C.3)
For what follows we define Ug == Us, and Üg := U. for k e N. We will assume that our dynamics f. & are continuous in both (zr, t), in addition to the nominal signal U, being continuous in t. However, since our signal U, is discontinuous in due to the addition resulting vector field F(x,1) := f(x,() ((x,()( is discontinuous on l. Hence, we will assume that the integration strategies will, order to generate the (k + 1)-th iterate given the k-th iterate, only evaluate the vector field F(x, t) on the half-open interval [tt, fx+1). In particular,
the Gaussian noise, the
will interpret F(x, tx) using the right limit F(x, tf) lim,-eg F(x,t) t) and F(x, tx-1) using the left limit F(x, (in) = lime- F(x,().
Euler seheme. Consider the standard forward Euler scheme used to discretize (C.1) with constant step-size T:
X Xx +[S( (2) ((8),4))4
= Xx + +T[ (x) + 4(A,, 4x ÜnJ + VF4 (&e, (x)wk,
which corresponds to the standard Euler-Maruyama discretization of the Itō variant of (C.2) with
g = Φ.
Heun scheme. Now consider the Heun scheme used to discretize (C.I), again with constant
step-size 7:
X Xx T[S(X, (Xx (x) (Xx,
tw (Xx-1, X Xx T [f tx) 208,, (R()))) + 2
2
where Úū.1 Ūr+t + WF. Using the shorthand (xx, (x,1) $(x, l)Ur, we see that:
X = Xn TF FTTFF(A (8,, tx) + ( Wx)
V [ ((tx),tt+)) Wk.
X -x++ [F((),,t) tx) F ([
2
which corresponds to the stochastic Heun discretization of the Stratonovich variant of (C.2), again with g 2 +.
D Proofs of Main Results
D.1 Auxiliary results
We first recall a standard formula for the variance of Gaussian quadratic forms.
A similar argument can also be made for dependet policies.
23

---

Proposition D.1. Let Q be a d x symmetric matrix and w ~ N(O, Ia). Then,
- =um 21
Proof. We have that Ew (tx(Q)) - w Qw) Ew Ew (w)) - From [cf. 49, Lemma 2.2] we obtain the identity E.d (w Qw)² = Qu))==tt((((( 2 tr()) which concludes the proof.
L
Proposition D.2. Let (X1) for & b denote -a-uue stochastic process, and let Rd x [a, b] + R be an L-Lipschitz function its domain. Suppose that for all Lı, la e [a, (] we have X X1 & MIn 12]. Then,
EJ
S. ((r((,,) - r(Xa, a)d & (b- (b- a)
(D.I)
Furthermore, suppose that for some O & 1, we have N == (b- a)/T N.. Define in for n E N... Then we have that the left-endpoint Riemann sum satisfies:
b b-a L. B(r(Xx,
Blr (n)] L(1 (11))
Σ
(D.2)
N 2=
Proof. First, have:
J. ( -r(de L ["
[ XXa
E
dl
Ja
t-a a
B||X, Xa xalct+ L S( (t - at
< L
( T, X Xalldb+ + L(b a))
2
2L
(b (b- ab-a)3/2 3/2 L(b - a)
3
2
where (a) follows by Jensen's inequality and the second moment assumption on (Xt)s:
L J. E| - -
& S. VEIX, Xal & LM S. [((t- (t--))) 2LM
(b (b-a) (b-a)
This establishes (D.1). We now turn to (D.2). We write:
S. t)
N-1
N-1 t]]d-- Σ B[r(Xen,
Σ
B[r(X,,, tn)]
B[r((t,))-rr - -r()))r(Xee. fm Jdt)
B[(r(X,,t) - r(Xe,,, tn llt
NI. 11 + 72]
1/2 +1) (b- (b-a) 1)
where (a) is from (D.1) and (b) is from the assumption that & 1. This establishes (D.2).
Proposition D.3. Let (X1) for & b denote an Ra-valued stochastic process, and let gd x [[a,b] IR be B-bounded and L-Lipschitz function on its domain. Suppose that for all t2 E [[a b] we have [E|[Xe, Xxall & Mlt1 tal. with M > 1. Then, b- <
E. (f. (((XXx))a)) - -(b 16-aj -(b-a)El(="X < !21 (bb (6/2))/2 5/2
(D.3)
24

---

Proof. We first decompose:
B (s. ((x,)))dt)
= EE (b-a))(Xa:a)+ S. J. (((x, a)){t) ?
= =b (Xa, (B
(r((, -x(aa - a)]dt
+((6 - a) r o) f( [r(Xe,)) -
-r(Xa. -x(xx,)))x(Xx
Next, by Jensen's inequality,
(1((xx)))) rxxxxxx a(]dt
< (b- J. (((((X)))))(XXxoo)))
< (b- a) S(( (E[I)X, - xxa1]] 1]
(t )
ab- [( (M a) (t a)))
(0) & (bb
where in (a) we use the assumptions that M > 1 and < 1. On the other hand, by another application of Jensen's inequality.
2(b - a) R r(xa. ((x(xx.)-
< < - = [ [[(r(x,,))) - r((a, adl
2BL | (B ()/2 (b - a)$1
where (a) uses Proposition D.2, specifically (D.1), and (b) uses the assumptions that M > 1 and b < 1. The claim now follows.
Proposition D.4. Consider the Itō SDE (X1) defined by dx, f(d,, g(Xe, t)dB, with
sup max max(1 B. (x.t)eR<x(a.b)
For any to, t1 E [[.
E X|, -xaall & 2B² [(t1- (
+//t, - tol]
Proof. Assume wlog that t > to. We first decompose:
E1, - X.oll 2 = = 1S.0 f(Xe,t)dt+ S.h 9(X,,))
< 2E
S." "
t(xx.
+ 2E
S. x(Xx,t)
Next, we have:
I.S.a* s(Xx..)ae/ < S.th II(Xx, & (4-44) (L1 - E I s(x,, 2
< <Ba( 1.)2.
On the other hand, by Itō isometry.
*S,," s(Xx,))dB. = -Saa 1) & B²1 (t)).
The result now follows
25

---

Proposition D.5. Consider the Stratonovich SDE ( defined by dx; = toddB. Suppose that for t & [a.b]] b], the map x + g(z,t) is C. and that
m Σ Oag (a, ()1(1))) 10(2.))) .100 < <B.
sup ((a.
max
f( f(x,t)
+
2 k=1
Then for any 1o.11 1 b.
E B| - XXa $ 2B2 [(t to)² 2 + +1t - -tol]
Proof. Consider I)X,, is pathwise equivalent to (Xi),, i.e, (Xi(w)), = (Xi(w)), for a.e. w, and hence the result follows from Proposition D.4.
+ Ek= Oug" (x, ()(()))))))))))))))))))))))))( The Itō SDE dX, J(X,,
D.2 Proofs of Lemma 4.1 and Theorem 4.2
We first restate and prove Lemma 4.1.
Lemma 4.1. Suppose that f.g are bounded and ug is C We have that
1 (x[uo](x, -r (r, V +)]] 0 T LEM. xr, (0, where the (-) hides factors depending on d, the bounds on f.g. and l lea..
(4.2)
Proof. We first introduce two pieces of notation: OC)) and O (.). The former O(.) hides constants that depend arbitrarily on the dimension the bounds on and g, and [uellca, whereas the latter O" (-) in addition also hides constants that depend polynomially on ||w]. The latter polynomial dependence important when we take expectations of powers of O (-) terms, since E| w|l is finite for any finite p E N.
Setting À == (ieer T, IRd I and writing uo. second-order Taylor expansion yields:
u(bttt; 7) - 1(xx,))= = Du(z, 1) 2
0(12)
Du(a,1)ā (Vu(I, 1), Eur - x)) + dru(x,(),,
A DZu(a, ()} (( - x) V22 ((isr
2r ( Vu( 1), Ees - x) 70
Plugging in - a- we obtain:
f((,,t -M(z - - (z, (olct))))V( (791))
=T =T( tr (H - + g(x,t) ) ((z, V
+ + +
where in the last equality we used the definition of the PDE residual from (3.1). Hence, (2,t)
Em -((x)))))))) (.t.::----........--.//////////-
= E. (R[ - - (
+ -0(7//3
=-² (M(U((z)))) + + L o
0(7312)
= =>" (
1 (n((((H,,,)) ( +
0(731))
where the final equality follows from Proposition D.I. The claim now follows.
26

---

Next, we use Lemma 4,1, along with order 1/2 strong convergence for EM integration (Appendix F) show the following result.
Theorem 4.2. Suppose that Sig, ha CO.1, uo C and & 1. We have that:
LEM.- (0) - J B [((RAuo(Xx).))) + where the O() hides constants that depend on d, T, and the Hölder norms of J. ho, and ug-
1 J (H(,,) v'uo(X,, (((())))))
dt+ O(T)12), (4.3)
Proof.
To start, we have that:
N-1
LEM, (0)
Σ Ex, [T (0,,,ttt
N
(a)
Ex, (R[[u))(Xn) bn))² + (H V'uo(Xn (n))²] 0(+212)
N
2
where (a) comes from Lemma 4.1 which holds since (i) g e C implies fig are bounded, and (ïi) E C Our next observation is that the map Rg is also Lipschitz continuous over the domain Id x by our assumptions f. h [o] e C and ug C Let us call this Lipschitz constant LRo- which depends only on the norms [ Igllo., | lo),,
/ |ullca.. Continuing from above,
E x [X,,n, tn)
(Ro bn
( B(B,), Ro (t,,, [n]] )1 E(B.), (Xn,tt) RX t )=]
where in (a) we consider the process {X,,) } as being defined over (AW,), {Br... B,,,i In in place of the process (VTw,), yn in (3.6), which is distributionally equivalent, in (b) we consider the forward SDE (X,), from (3.3) as being coupled with the process {X,),, via the same realization of both Brownian motion (Bt), and Xo Xo vo. Hence we have:
[ [Ro (n)] E(D,), (X [n] < EIB.,,
Ro (o (X,,, (I]
< Ln E.B.). X.]
< LRo V BIB... 1I3.
Now, by definition, since the functions f.g CO.,, then the pair (Jsg) is EM-regular (Definition F.1). From Theorem F.2, we have that (8,), is strong order 1/2 convergent towards (X,)., and hence
[( xx,, [|]]] < B(B,), (eex).(mxxx [m..,.
X.,, |2] < C²T, where the constant C' does not depend on but can depend on d, T, and the norms co. and ||g|oo... Consequently.
LEM, (0) = N n=0
Σ [Ro 02)
(D.4)
Our last step is to approximate the sum N NN11 in (D.4) with the integral. To do this, we will use Proposition D.2. We already have Rs is Lipschitz over gd II.I. Furthermore, since f. C0011 they are both bounded over the domain, and hence Proposition D.4 shows that X| -X & ((1)lt tol for any tu, & e I. Therefore by Proposition D.2 and the assumption J <
T L' t))
1
1 N-1
Σ [.))((X,,,,))]]]0
(a) LEM, (0) + O071131)))))
N n=0
where (a) from (D.4). The result now follows.
D.3 Proofs Lemma 4.3 and Theorem 4.4
We now restate and prove Lemma 4.3.
27

---

Lemma 4.3. Suppose that f. 9. and he are all in C and us is in C We have that
x-2 Houn,
(0,x, (0,
(Rl (x, +
(4.9)
where the O() hides factors depending on d and the Hölder norms of J. 9. ha, and
Proof. Similar to the proof of Lemma 4.1, we let O(+) hide constants that depend on and the Hölder norms Il Igl. [hellc,, and | Iaaa, and O* additionally hides constants that depend
polynomially on ||w
Our first step is to check that hi C under our assumptions. Recalling that (z, t) ha(x, t) tr(HH(x () (z, this is ensured if hg, H e C and ug C which holds since e C implies H с C Furthermore, l O(1). For what follows, we drop the dependency in the notation on 0.
We next Taylor expand 8s+7 up to order т terms. To do this, we observe that by our assumptions on g,g, S, u, the functions F, are both in C*.,1, Hence,
=1)+ - + -21))+++++++++ 0 F( = +
By similar argument,
Hence,
1 -
20st - 21 =
(D.5)
A straightforward computation yields:
Dg(a, [) (rgp::/ee..Diluuucooooooaaaaaaa
and therefore:
} (IGG t) Do(a, () w g(z,t) (Tz,t)Tv²u(aa) t)
(((,)t)
Substituting the above into expression (D.5) for 2 2t.
f(x,t)
2007 2
(((E,1)- - (u,(H(,,t) () u(, ((),
Dg(a, (xx
(u( D) [g((,, 1) w (t, (')((x))TTvu(tee) (V
VF| [v(vvta, () (T((,,)))))))))) D +7
g(I,t)w
(D.6)
Setting A == (āver - 2, 7) Rdd+ we have:
u(t t+ 7) u(x, Du(x, t)A L A (
2
Du(a, 1)ā (Vu(z, 1), 2, x)) 0,u(x,1)T.
A 1 (),()(=)))))))))")
2 ((eerr - x) T 422 u (1447 x)
2
2 (0, Vu 1) Seer 2) 720 ()):
from which we conclude,
u(Kxur,l+ ( () - u(x,1)
= =VVu(z, f). z) d,u(z,t) (fa 0(?2))
28

---

(v(( Do (o(,, d,((x,t)+ w g(x,t) (((,,)))) (( (
On the other hand, from (D.6).
Uct - y1 = (xx,t) - (((((2,))V
2 [Vu(a.t) Dg(z, w (0) (ee,tt)((((( ()(x, [)]
VTVu(a,t) f)r,t)w g(z, (T
Hence,
(f Dt+r au(iuttt++ u)x,t)) (fe+r - 31)
(Vu(z,)),, + 0,u(2, (MNo::/----((DDDD/(.--)))))))))- 1
+[ul(2,t)T++0% (-
To conclude,
End Buc(u((err ++ - === Eny (R( +0()))))) + 0(7913).
Theorem 4.4. Suppose that J. 9, and hg are all Cl.,, up e C and We huve that
LHenn r (0) = - J. [XH[ d(
1
0(7
(4.10)
where the O() hides factors depending on d, T, and the Hölder norms of J, g. hg, and ug.
Proof. The proof follows the structure of Theorem 4.2 closely. We start with:
L (0) = N n=
Σ
ExplT & (0, (n) fn)
(a)
Σ
Ex; (H[[ua]( tnn)) 1(2
N
n=0
where (a) follows from Lemma 4.3. Next, define Ro(x,t) (R (xx t))² and observe that Rs is Lipschitz over Rd x I due to assumptions on g. heg ug. Hence, have the decomposition:
xx; [Ro bn) (a)
[ tn)]
E(B,) Ro X; [n]] (R,,), Ro(X;, tn,)]
where in (a) and (b) we take same steps Theorem 4.2: (a) considers the process {X;},, as being defined over (AW))a (Brownian increments) in place of (VTwn), (4.7), and (b) couples the SDE (X;), together with (X;),, via the same Brownian motion (Bi); and initial condition Xi X; = øo-
Hence, we have:
Ex: (R [.,]] E(B)) [R,, +t)] 0(1)
[Ix; (())) 8 1.].
Since f.ge C1.1 then they are by definition Heun-regular (Definition F.3), and hence by Theorem F.4, [I {= < E B,B,) [maxx.e(......)] XX;, -x; & CºT, where does not depend on but depends on d, T. and the Hölder norms on f. Therefore we have:
1 N
L Heun. (() =
(D.7)
Now we can finish up using the same ending as Theorem 4.2. The only thing we need do differently is to control the second moment EIIX;; X; 1 Since f.g I Proposition D.5 yields that

---

B|X; - -X,, 0(1]lt tol. From this inequality and the Lipschitz continuity of Rs over the
domain Rd I. by Proposition D.2 and the assumption &
T L (x 1
Σ (B),, (x,],,,, R=0
Lounn (0) (
where (a) is from (D.7). The result now follows.
E Analysis of Multi-Step BSDE Losses
We now present the derivations supporting the analysis in Section 5. For what follows, we define the follow FS-PINNs loss:
1 J." ())
LES-PINNa Ezo-po, Bi T
(E.1)
E.1
BSDE loss and Euler-Maruyama discretization
Proposition E.1. Suppose that ug E C We have that:
LBSDE, (0) & LFS-PINNS(O).
(E.2)
Proof. To start, we abbreviate Ro(x, = (a t). By application of Itó's Lemma, we have:
un (r; T) no (X., L. (x,, t)dt J." (x,,(X,, 5(X,,,)) = S. Ro(X,,
which immediately yields the following identity:
1 LESDE, (0) 1E Be T [' (d((x
(E.3)
Thus, the BSDE loss equal to averaging the square of the accumulation of the residual error Rs along the forward SDE trajectories. Hence by Jensen's 's inequality, the BSDE loss is dominated by:
LASDE, & Euroopp T S.' (Roo(X.)))2al (IN.
Note that while the relationship in (E.2) is pointed out in [5, Section 5.2.3], the implications of this inequality are not discussed further their work.
Proposition E.2. Suppose that f,g, hg CO.1, up C2.1, and & 1. Then,
LISDE, (0) LYS-PINN. (f) + 0(71/2),
(E.4)
where the O(-) hides constants depending on the Hölder norms of f. g. he, and us.
Proof. Again, we abbreviate Rg (x, = (octt)=[[uap)) [[ue](xxt)) t). By our assumptions on f.g. ha E C and Ug C we have that Rs 0. with [R O(1). Also, since f.g C by Proposition D.4 we also have that E| X x2, < (())lto - t1l for fo. f1 I. Hence by Proposition D.3,
(
E
t)dt) = - PR tt)] 5(75))).
(E.5)
Furthermore, since Rg E C we can readily check that Rg is Lipschitz on its domain as well, and hence Proposition D.2 yield:
T L. T ())
Σ (ta + 127
(E.6)
N 0=0
30

---

Therefore,
L = Ezo BB Bi
1
N-1
Σ
(L,m
Ro(Xt
N7
N-1
-Ezo-po Bi NT
Σ
(x),
[using (E.5)
N-1
Ezo-po Bi N
Σ
R (,)) +
n=0
1
J
Ezo-po.Bi B T
0
[using (F.6)
= LINN (0) + 0(7212).
Proposition E.3. Suppose the assumptions of Proposition E.2 hold. Then,
LBSDE,T & LBSDE, (0) 0(7)
(E.7)
where the 0()) hides constants depending on the Hôlder norms of f. 9. ha, and ug.
Proof. Follows immediately from Proposition E.1 and Proposition E.2.
Proposition E.4. Suppose that f.9, he CO.1, ug C..,, and T & 1. We have that: L (0) LBSDE, 0(7
(E.8)
where the 0() hides constants depending on the d, T, and the Hölder norms of J, g. hg, and ug.
Proof. We consider the joint forward and backward SDE (cf. (3.3) and (3.4)) with z = (X,, Y
X 5(x;,() dt d
((X,,t) dB [g
(x 0))
xo
(x. s(X,,t)
Y
=F)(2
=Go(2
Given assumptions on f.g. ha, ue. have that both Fø, Gg CO.1, Hence, the pair (Fa, Go) EM-regular (Definition F.1). Therefore, by Theorem F.2, the EM-discretization {(X,, Y()),, (ef. (3.6)), coupled with (Z through Brownian inerements (AW,), satisfies order 1/2 strong convergence to (Zς);:
max max[]|X, Xt,, 12 x2/1e) ne(O.....N)
1/2
Now, define i N == us (RN, (N, T) 8 and Up up Xr; T) - Y.
E| E[ (V (UN 10---771
- 2 VT)
Hence by Cauchy-Sehward
E [[ & E + -1
Since ug E C the function is |ue c... -ipscchiz and therefore:
E [ Vr) < (uo N. T) - uo ( -
< 2 up 8... [[ Xx Xr 2 yf)
< 20 uu (a.1)T.
Furthermore us is also -uollc..-boundee and hence:
E 2 < 2 up 2 22E[[1?11].
31

---

Next, since Fg. Gg are both bounded, then by Proposition D.4, we have that [[ O(1). Putting these bounds together yields:
I
El ((++)
(F.9)
since 7 < 1. To finish the proof, we observe
LFMx.-(0) Bro-rar 7
V
Exg-
T Vr
[coupling with Brownian increments]
Exg-po T 2 + ((7
[using (E.9))
LBSDE, (
Proposition E.S. Suppose that f,g, hg CO.1, ug C².1, and T < 1. We have that:
LEM. (0) LBSDD,, Bias(() 0(71/3)
(E.
2T JT¹ B[Lr((((X,,0) vo (X,, ())
1
Bias 0)
Here, O() hides factors depending on d, T, and the Hölder norms of f. g. he, and ug.
Proof. We have the following:
LEM, (a)
J.¹ T T
E ( tr((H((),,)) ((uo(X,,)))3)) (1"p((XX.))--vooxX.)))) dt + 0(T
LES-PINN (0) Bias ( + (27))
LBSDE,- (0) + (ia (
where (a) holds from Theorem 4.2, and (b) holds from Proposition E.2.
E,2 Stratonovich BSDE and Heun discretization
We first define the Stratonovich variant of the FS-PINNs loss (E.1):
Le.s.PINNs(()) Ex B. T ['" (R]uol(X),
Proposition E.6. Suppose that us C We have that:
Ls. LS.BSDE, & (INN
Proof. We mimic the arguments in Proposition E.1. Abbreviating Ro(x,t) R[uo]<x, t) and using the Stratonovich chain rule, we have:
[ ho(,), in ()'xxx)))''
un - -x..))
('vNo(X2,1) 9(x;, dB.
and hence the following identity which parallels (E.3) holds:
S. R Ro(X;; T
1
Ls BSDD,T(O)) Eoo-ao B,
Now we simply apply Jensen's inequality to conclude:
1 S
2
Ls.BSDE.T()) Ezg-Ho.
T Ro(Xi,,)dt))
Exo- T [." (Ro( (SPS.PINNa())..
32

---

Proposition E.7. Suppose that f. CO.1, g C uo C and T & 1. Then.
LS.BSDE,- (0) LINNS. 0(T))2)
(E.11)
where the 0() hides constants depending on the Hölder norms of f. 9. hg, and ug.
Proof. The proof is nearly identical to the proof of Proposition E.2, but with Ls BSDEE, (0) taking the place of LBSDE,-(()) and rs Pes..PINNN(C (0) taking the place of (ISSPINNs(f)) The only notable
difference is need to establish the condition Xi - X2, 112 < ((1) -1l for 40,4 e I. By our assumption that f C001 and g C*I.1, we have that both (2x,) + Ex= (x, t) and g(x, 1) are bounded, and hence the condition E X% Xi, ((1) holds by Proposition D.5.
Proposition E.8. Suppose the assumptions of Proposition E.7 hold. Then,
LS.NSDEE (0) & Ls LS.NDDE,, O(7212),
where the 0() hides constants depending on the Hölder norms of J. 9, ha, and ug. Proof. Follows immediately from Proposition E.6 and Proposition E.7.
Proposition E.9. Suppose that fr ho CO.11 g C
C and r & 1. We have that:
LHeung (0) LS.EDEET) ((7
(E.12)
where the 0() hides constants depending on the d, T, and the Hölder norms of f. g. hg. and ug.
Proof. We consider the joint forward/backward Stratonovich SDE Z, (X;, Y,,0) from (4.6) the form dz, Fo z t)dt+ Go(Z;" .t) dB, We first show that the pair (Fo, Go) is Heun Hegnlregular (cf. Definition F.3.. A sufficient condition is that (a) Fg E C and (b) Go E C1.1 For condition (a), it is equivalent both f and (2, () (a 1) tr( H (a, uuguu ()) are in C..1; the former by assumption, and the latter holds since ho C 9 C.i,, and us C by assumption. Now for condition (b), it is equivalent both and Vue t) g(x, 1) are in C The former is again by assumption, and the latter holds since ug C and g Ct.l Hence by Theorem F.4, we have that Heun discretization { (2:-)) In from (4.7), coupled with the SDE (Z;"), through Brownian increments (AW, )n, satisfies order 1/2 strong convergence to the SDE (Z7 )1, i.e.
1/2
IF
max (Ix;; X2 ne(D.....NN)
max
12 1..0 vr,
< CC
where does not depend on r. The remainder of the proof proceeds nearly identically to Propo- sition E.4, with the only difference being that in order to argue B[ O(1), we utilize that Fg 00.1 and Gg 1 to invoke Proposition D.5.
Showing that Heun, (0) LS-BSDE, (0) + 0(r)/2). We start from (4.10) and following nearly identical arguments as in the derivation of (0.
Proposition E.10. Suppose that f. 9, and hg are all in C C and r & 1. We have:
LHemn (0) LS.BEDE, (0) 2(2))))
where the O()) hides factors depending on d, T, and the Hölder norms of J, g, he, and ug- Proof.
We proceed similarly to Proposition E.5:
THomn, (0)) (a)
T L. E(R[uo](X;; 0
LL..S.PINN. (0) (12
( LS.BSDE. (0) ((T
where (a) holds from Theorem 4.4, and (b) holds from Proposition E.7.
33

---

F Strong Convergence of Euler-Maruyama and Heun Integration
I= = [o, 7] denote a time interval, and consider functions a IR x H IR and b: R IR ddxm
which define the following SDE:
dx. a (X,, 1)de + b(X,, dB,, Xo ~ Do-
(F.1)
where (Bi) is m-dimensional Brownian motion and Dp is an arbitrary distribution over IR
with bounded second moments, i.e., E Xoll² < co. Here, the pair (a,b) is used instead of (f.g)
to avoid confusion with the forward SDE (3.3), and the notation denotes that the SDE (F.1) is either to be interpreted as an Itō or Stratonovich SDE. We write 6a E x I + R' for E [m] so that = (b'....,b')) .....,,"),, ) i.c., 6º (t, is the k-th column of the matrix b(t,z). We consider discretization time т (0, such that
[T/T] N.. We denote the timesteps (tn In-o and Brownian increments (AW, ) W-1 as ln := nr and AW, = B B,,. We will also often utilize the shorthand notation a, (x) == a(x, (n), by (x) = 6(x, t..) and 6y (1) == 6a (x,,t) for n (O.....N),
In this section, we review standard results regarding convergence of basic stochastic integration schemes (Euler-Maruyama for Itô, stochastic Heun for Stratonovich) for the SDE (F.1).
F.1 Euler-Maruyama Convergence
The Euler-Maruyama scheme for integrating the SDE (F.1) interpreted as an Itô SDE is the following discrete-time process:
X an (Xx)) (X,) Xo Xo-
(F.2
The order 1/2 strong convergence of the Maruyaaa process (F.2) to the Itô SDE (F.1) is thoroughly documented in the literature. Concretely, we will state result from [34]. First, we define the necessary regularity condition on the drift and diffusion terms
Definition F.1 (EM-regularity). The pair (a,b) with Rd x I Rd and g Rd I Relmm is EM-regular if there exists finite K1, K2, K3 such for all z, Rd and s,t EI.
|a(x, t)1 + [1((,,))lopp < Kı (1 |x]),
la(x,t) - a alv.t)ll [b(2) b(y, < K2Il --ll- yl.
(a(r, - a(x, [|b(,, - loop < K3 + | 1/2/2.
By definition the Hölder class C we have that the pair (a satisfies a, C then (a. b) is EM-regular, although Definition 1 is weaker assumption. With this notation of regularity in place, we have the following order 1/2 strong convergence result.
Theorem F22 ([34, Theorem 10.2.2]). Suppose the pair (a,b) is EM-regular (ef. Definition F.I). Then the ltô SDE (X+): defined in (F.1) and the Fuler-Maruyama discretization (Xn)a defined in (F.2) satisfy the following bound:
1/2
|X - -0xxxp1 (O.......
xx12) ECVT,
max
where the constant C does not depend on T.
F.2
Stochastic Heun Convergence
The stochastic Heun discretization of the Stratonovich SDE (X,); defined in (F.1) is the discrete-ting time process with Xo Xo and:
Yn-1 Xn, + am (Xi)+ by, (K))AWn. Xn-1 X,, + ブ an (X.) am. (8 )] r
(F.3a)
1
2 [b [b,. + bn-1 (Pxe)) AW..
(F.3b)
Analogous to the order 1/2 strong convergence results in Appendix F.1 for the EM discretization of the Itō SDE, we also have similar result that holds for the Heun discretization (F3) of the
34

---

Stratonovich SDE (F.1). While we consider such result to be a folklore result, we were unable to find a specific theorem statement in the literature listing out a precise set of sufficient conditions on (a,b) for strong convergence to hold.* Thus, the rest of this sub-section provides a result and mostly self-contained proof that builds on top of EM results stated in Appendix F.I.
We first start with a sufficient regularity condition, which adds few extra assumptions to the EM-regular definition (Definition F.1).
Definition F.3 (Heun-regularity The pair (a, b) with IR x I+ TRd and TR x I + Te is Heun-regular for every t E I and k e [m]], the map + (t, is C* (R.), and there exists finite K2 i [5] such that for all x. y IRd and 8, I I::
|a(x, t)l [|((,,)) op $ K1,
1a(2,, () - aav, 1b(x, b( () op & W
Ia(x, a(1, + [b( b(x, up < Ks(1+ + [x] I[s t|
110 (x, 0, (y. 1) lop < <aa|l=-10 14ll--11 yl,
lop < K5 + |>| [s t//1/2
We note that from the definition of the Hölder classes C and C that if a E C and be C².1, then the pair ba is Heun-regula The following result is the main convergence result for Heun. Theorem E4, Suppose that the pair (a, b) is Heun-regular (cf. Definition F.3). Then the Stratonovich SDE (X,), defined in (F.1) and the stochustic Heun discretization { (Xn)n defined in (F.3 satisfy:
1/2
|Xe, - -xx11) &CVT,
E
max
(F.4)
RE(O.....NN)
where the constant C' does not depend on T.
F22.1 Proof of Theorem E4
Our proof of Theorem F.4 is based on the following reduction. By defining:
1 m Σ k=1
a(x,t) : a(x,t)
az (z, (x,
the Itō SDE:
dx, (xX,, t)dt b(X,, t)dB, X = Xo.
(F.5)
defines an identical process as the Stratonovich SDE (X;); from (F.1), i.e., (X,(w)) = (X,(w)), for almost every l,w. Furthermore, we can consider an Fuler-Maruyama discretization of (F.5):
Å :..=,,, ān (&x)) b, (& eo = Xo.
(F.6)
As the Itô SDE (F.5) and Stratonovich SDE (F.1) identical, then also have that if the pair (a, b) EM-regular (cf. Definition F.1), then the EM discretization (F.6) is order 1/2 strongly convergent to the Stratonovich SDE (F.1). Hence, this reduces the problem to comparing the two discrete processes (X from (F.3) and {8,), from (F6). In particular, we can show that:
1/2
E
max 1&,, 1,112 ne (O......,.)
< CVTT
where again the two processes are coupled under the same Brownian motion (B:), and initial condition Xo. then by triangle inequality we have the desired result Theorem F.4. One advantage of this proof strategy is that we only need to study the evolution of two discrete-time processes, which we can do with purely elementary (discrete-time) martingale techniques, avoiding the need for any stochastic calculus. Indeed, the main tools we utilize are the following two results.
"The closest statement we were able to find in the literature is [41, Theorem D. 12]. which shows order 1/2 strong convergence of the reversible Heun method, which is modified version of the stochastic Heun method that algebraically reversible.
35

---

Proposition F,5 (Doob's 's maximal inequality prector-valuedd cf. [50, Theorem 3.2.21). Let (Xn) neR. denote martingale taking values normed vector space X with norm Il-IIx. We have that for any p E (0, oo) and n N,:
1/8
E X xc[x]
(x
1P
p p-1
Proposition F.6 (Discrete Gronwall inequality, ct. [5I]). Let fand and {,Bn) naN be non-negative sequences satisfying for some a > 0:
In E
Σ k=0
BR
ne N.
Then we have:
n-1 In & exee
Σ BE
nEN. N.
Here, we interpret Σx-o to indicate zero.
Our first step shows that Heun-regularity of the pair (a, b) implies EM-regularity of the pair (a. b).
Proposition F.7. If the pair ba,b) is Heun-regular (ef. Definition F.3), then the pair (a, b) is FM- regular (cf. Definition F.I..
Proof. We let K max
Ki We first check the growth condition on [|(xx)(]:
1 m
|ā(z,t)| = ((xx
Σ
a
(z,
2 k-
(a)
1 m
(b) (z, t)
(x, (z, ② K
m
K2m
K +
Σ
1
x
t)b
< K Σ kml
2
k=1
2
2
where (a) holds since |a(x, 1)]l < K, (b) holds since 116 (x, 1) < 1b(x, 1)pp & K, and (c) holds
since (g t) (v.t)l < K | yll implies that 10.g* x (x, t) op < K.
Next, we check the Lipschitz condition over x on ā(x, I):
1 m
|iā(,, t) - a(y,,)] < la(x,t) 1) alv. 1)1 +
Σ
la (x, (x, (x), (v.
2
(0)
K| I yl + 2
1 (z, (,b, 8 ()]l
m
Σ
[0,b (x, t) - (g, (b)
(8)
K r yl + K2m
K2m
x - y
1 - yll (K + K2 K2m) 12 vl.
2
2
where (a) uses the Lipschitz condition ||a 1) aly, f)l -I and (b) uses (0,
0. (y, t) lop & K||æ vl. [=*( K, 10,b (I, (lop & K, and (x, - 8a (v, &
K| - yl.
Finally, we check the Hölder 1/2 condition over on a(z,():
1
|ā(z, () - -á(x,t) & [a(z 8) a(2,t)) +
Σ
0.b
(a)
m
K(I+ [1|+|x/)|s--- - t| +
Σ
1 (x, (z,
8' (x,1)]
2 -=
1
Σ
(1,,5((r,s)))) O 0. (z.
k=1
36

---

② K(1 [(1+|=|)]s - t/
K2m
(1+ |x|||s s- - t|+ K2m
2 (1+ |<11|s - t/12
2
(9) (K K V | ||s 11/12
where in (a) we use |a(z, 0 1)
[ || t|, in (b) we use [b*(2, -b &
K(1 + |x Dss 4/. |00 (2, 1) up < K, | (z, -0, & K(1 + | - 1/212 and
(x, K, and in (c) use |s t| = |s < vT t/1/
Since the growth, Lipschitz, and Hölder 1/2 conditions on b(x, () are immediate from the Heun regularity assumptions, this concludes the claim.
L
The following result shows that the discrete-time processes (F.3) and (F.6) strongly convergent.
Lemma F.8. Suppose that the pair (a, b) is Heun-regular (cf. Definition F.3) and that & 1. Then, we have that the updates (F3) and (F.6) satisfy:
1&n En11 1/2
) <CVT,
E
max THE (O.......
where the constant C does not depend on
Proof. To start the proof, we recall that b nT, N N = [T/T] which is assumed to be positive integer, and AW, := B Bt..- We will equivalently write AW, VTw,, where {w,) n=o are i.i.d. N(O, mm) random vectors. In order to index the coordinates of both AW, and wn. we use the notation AW; and wh for i E [m]. refer the i-th coordinate of the vector. We let K == ++ Xool + max;e[s] K; to denote a bound on all the parameters from both Heun regularity and the second moment of the initial condition Xo- As with an, bn, we also define ā, (z) == a(x, t.) for n [N]. Furthermore,
will drop the hat notation reduce notational clutter and and write Xn, Yn Xn, Yn. and similarly X, X.: since we are not dealing with the SDEs (F.1) and (F.5), there no risk of confusion with this notation.
To avoid keeping track of explicit dependence on constants besides т, for a set of parameters Q we Co to denote finite constant that that depends only on the parameters listed in Q, and Sq b to denote & Cqb. For example, Ck.m is a constant that depends only on (K, m) (its dependency may be arbitrary however). We also let S denote & Ch where C is universal positive constant.
With the aforementioned notation in place, we have the following discrete-time update rules for (F.6):
X = Xn an (x(X))++++ (X.)
and for (F:3):
Yn-1 Xn + an (Xn)T + (X AWW,
1
1 [b,, (Xn) + buni (Yar)) AW..
Xn-1 Xn - [an (Xa) anti (Y ]]] 2
Let us define the filtration F,, := oww, 1n-1) for n E N, with Fo the trivial o-algebra. We observe key property of both (N,)nn and (X,), is that both Xn, and X,, are F. measaaaaabl Hence by the tower property we have for expressions A, B, ((w,) = [( FrJ] and E (w,,) ) FF.] E,,,, F... B w,,) ]. property we make heavy use of in our calculations. Another simple inequality we make use is that for any q e N, and any set of vectors Ul...., Ug.
Σ
1
which follows by triangle inequality and Cauchy-Schward
Heun update decomposition.
To relate the Heun and EM updates, we write the Heun update as: X Xx, + a. (X,, + b, (XW, En,
(F.7)
37

---

where En contains the residual terms:
1 En 2 ((n-1(Yn=1)) an (X,)]
r+ 2 -bx
d - (X,)
=.
We further decompose F., follows. We first write:
bn-1 (Y=+)) - bn (X,,) (bn=1(Yn+1) by (Yn+1)) (Yn())
b (X)))
Next, we use the Heun-regularity to expand the RHS above:
8 (Yn-1) (X,) 02 (Y - Xxa) + R
(X,) (x)T + bn AW,) R
where the remainder term R satisfies | & KI Y - Xx122 Hence,
(bn (bp(Y bn (xx))))
m
Σ (Yn+1) km1
JObK Σ R Aww,
Now the middle term further decomposes as:
Σ
Σ 0, (X,
AW
k 12=1
Σ 0, (X (X,,) (A AW 7 (k1=k2))
k kg-1
Combining these decompositions together,
1
E
(0n-1(Yn+1) -V. (X,,) Σ
0 (X,)
2
F=1
1 10n-1(Yn - bn (Yn+1)] AW,
(Yn(1) by - bx, - Σ dzb (Xn)) (X.))
T
2
2
2 i=1
1
m
m Σ
[bn-1 (Y - by, (Y AW,
Σ
(X.))
R AW
2
k=1
(X,,)
T m
m
(nn+1(Yn=1) - bn (Yn+1) AW,
Σ
0,b (X,,) an (X,)
+ 2 k-1
Σ R AW n
2
2
k-
=
E..3
Σ
bzb (X (Xa) (AW") AW 7 (()
Thus, (F.7) becomes:
4 X = Xn, (xx)++++++ by (X,)),W, E, Σ
E..
(F.8
which serves as the starting point for what follows.
38

---

Second moment bounds on error terms. Our next step is to bound the second moment of all the error terms separately. Here we make heavy use of the Heun-regularity conditions.
Bound on [E] F.y ? We write:
(0n-1(Yn=1) an (X,,)) = (1a+t1 (Y+t)) an+1 (an-1 (Xx))
1axt1 (x a (Xx,) + -
(xn)) (Xx,l
K|| X,1 K(1 IVVF
< K ((aan (Xx + |x,
<
Hence,
El (Y (xx,)l (r²(" xx
+m) + KK (1+ E|, 13) Sx. (nnI| 17T
which implies:
7
E 1
E amt (Y (0,n)))-,,,,)) (X.))1 Skm (1+
4
Bound on [E] Ex 12: We have:
E 1 1
[[[ 1 bn (Y
4
$ 7 " | 13) wn 12]
" |x,|| 2 Ian (Nx,)) 127²+
1b, (xn)) ww/l [2
K ( (1 |x, |2, + 13)1 wn 131
Skm (1 (1+ X 1?)
Bound E|| E 12
We have:
E8
b.2
2
Σ d b (X,,) (X.,)
73 B b 18p aa
-
3
(xxx))
Sk,m -3.
Bound E 3 12 Recall that the residual R satisfies RR < K| X.I13. We have:
E 5.33
E
Σ R A
S Tm
[[[ 12 || k
S T 3
Σ
[[[ n+t Xx.1 (1
k=1
$ k 2 m Σ
m
E [(((la, 16 (ow (oet
Bound on [E] E...
We have:
ba (X.) (Xn) Ww 1l (00t /(
pb.4
T
E
m Σ
Sk.m 72.
n
4
39

---

Second moment bounds on the Heun process.
We now use (F.8) to write for any E [N]]
E x 2 E
Xa
(x(X))+++
Σ (XX))WW;
E
1=0
n-1
x <BI + +
Σ (X,)
+E Σ
(X)))
+ E
Σ
E
i=0
(F.9)
We focus bounding these second moments, using the fact that for n [N], we have Tn & ~N = x[T/r] & T. First, we have |a:(X:)) & K + Kzm and hence
2
E
Σ a. (x) 1-0
Sn' (x)) Sum (nT)² SK.m.T 1.
4=0
Next, since (b,(X;)AW; forms martingale difference sequence (MDS).
s-1
E
Σ
(X,)
Σ Blb,(x;)1 Skm SK,mT 1.
Next, we have the following bound using the second moment computations for the error terms:
n-1
Σ
E
+ EE
Σ
E
n
Σ
(E
+
13)
Skm
n
Σ
(1+E|X
1=
100
1=
1=0
Kkm
( + TeT
Σ 1=1
E|
S K,m,
T+
Σ
1-0
n=
1-
n-
E
Σ
E
b.2
E
h.3
Σ
Skm (nT) Skm,
Furthermore, since {L
}
an MDS,
E
Σ
E (NT)T SKmTT
Combining these bounds together in (F.9), we obtain:
X SKm.T 1
2
A=0
By the discrete Gronwall lemma (Proposition F.6), we have:
xI < CK,, exp((r"" Ck,,,.)) & exp k.m,T) SKmTT 1.
Hence, we have shown that:
E|X,,|| SKm,T 1.
max ne pe(.......N? N)
Final result. start with:
Our goal now is to estimate B[A for A, == maxxe[m]lldsl,, where o; 8, X, We
8n-1 Ant1 Xn+1
- + (X)) (Xx) (X,) -b, (X,,) AW, E,,
Hence for n E [N].
n-1
|8, 12 1==
Σ
(ax(X)))axxx)))
-
Σ E
1=0
1=0
1=0
40

---

Σ (X.) ((x))]]
[ -
Σ
SK 7 1-
16,112
Σ [4,
+ Σ
m-1
|
Σ
12
Therefore,
A max Sv 2 SK T ke[n]
n-1 Σ
||0:
2 +
Σ
| + Σ 1
3
im0
[b:(X;) - by (x:) AW:
Σ E
max k
+ max
1-0
ke(
1=0
SK 72
Σ f=0
|
E1
n
max ke [n]
AW
max ke[n)
Σ 1-0
F
(F.10)
Now, using nearly identical arguments in the second moment calculation of the error terms, in addition to the uniform bound on ||, for n .N......, we have:
1 + E ) SK.MTT.
E n
On the other hand, since both [[b,(X,) - b. sequences, using Doob's maximal inequality (Proposition F.5),
(A.)
J and (E("") } are both martingale difference
E max
[o,(X.) bx.)) AW.
E max ke(n]
i
b,1
ke
Σ
[ - bx()))] AW;
Σ
E
1-0
1-0
E - AW
-x.)))
Σ
E1 D
|
-
1=
SKm,T
118 SKm,T
E| B[ +
Hence, taking expectation in (F.10) and combining the previous second moment estimates, have:
E 2 SK.,T
E[ +T.
From the discrete Gronwall inequality (Proposition F.6).
[[A?] & CK.m.TT exp(nr Tm,,) SKm.T T,
which completes the proof.
We can now complete the proof of Theorem F44
Proof of Theorem E44. We have:
X,, _, X
(0) X,,
-x.) -X,,- +X) -aa);
where (a) holds since the Stratonovich SDE (Xt) (F.1) and the Itô SDE (X,) (F.5) are identical for a.e. (w, t). Hence we have:
|x 12, 2
[_ - -2,1122
1&G, -x,
max
max
max (O...NN
no (o........))
ne ne(0.... (O.....)
41

---

Next, since the pair (a. is Heun Heun-regular by assumption, then the pair (ā. is EM-regular by Proposition F.7. Hence by Theorem F.2, we have the EM discretization {Xn)n is order 1/2 strongly convergent to the SDE (Xt),, (. ie.. E () I Xn, 1 & CT. Furthermore, by Lemma F.8 we have E (0.... [|x,, X,l1 2 < C'T as well. Note in both cases, C, C" do not depend and hence the proof is complete.
G Experimental Details
We use each algorithm train 8-layer neural network with 64 neurons per layer and swish
activation [52] to model the solution u(x, 1) of a PDE. The boundary condition is enforced by adding
boundary condition penalty Ez-zr [uu (z. T) )(x))]+ +E[]]]]+[[Vaaaa [o T) Vo(x) [2] involving both
the zero-th and first-order values of Φ [4], where the distribution jí' is taken over each method's approximation of the distribution of XT. Additionally, following state-of-the-art PINN architectures practices [7], Fourier embeddings [53] with a 256 embedding dimension and skip connections on odd layers are used. We use a trajectory batch size of 64, translating 64 realizations of the underlying Brownian motions. Additionally, we utilize sub-sampling batch size of 1024 for the batched algorithm runs. We use the Adam optimizer with multi-step learning rate schedule [4] of 10 10-4 and 10-6 at 50k, 75k, and 100k iterations, respectively. All models are trained on a single NVIDIA A100 GPU node in our internal cluster, using the j jax library [54].
G.1 PDE Test Cases
Hamilton- Jacioonn---oo------ Beeli-n (HJB) Equation. First, we consider the following Hamilton-Jacobi. Bellman (HJB) equation studied in [4]:
uu(( (1, () -T[VVV"'(,, [[ux,, E le [0,
For the terminal condition u(x, 7) = <6x) In (.5(1+ [x|²)), the analytical solution given as u(z,t) t) in (E [e|e ( (-9(2) v2BT-))
(G.1)
The HJB PDE related the forward-backward stochastic differential equation the form:
dX, add:; t [0,T],
dY |2,l Z, dB, t [O,T).
where 7' 1, a= == V2, Xo = 0, and Yr ((X)). Additionally, the Stratonovich SDE is given as:
dX; = ood t [0,0
dy 2 2
[T de oZ! o dB.,
2
In our experiments, in order to compute the analytical solution (G.1), we approximate it using 1 Monte-Carlo samples.
Baack-Scholes-Bareee Barenblatt (BSB) Equation. Next, we consider the 100D Sew--Shhooll---rrnbhhie (BSB) equation from [4] of the form
0,((2, - Tr diag( t)] ( DER. t [0. T].
where x is understood to be coordinate-wise. and diag(v) ) is a diagonal matrix with diag(u); v, Given the terminal condition u(x, T) ((x) | the explicit solution to this PDE is
u(x,t) t) exp ((r o (())
The BSB PDE is related to the following FBSDE
dX, = a diag (i) e [0,T],
dY; r 1o=r((i (Y, 2, x,) dt + oft diag(X,)dlB, te (0,T))
42

---

where T = 11,0= .4. r= .05, Xo =
and Yr = ((Tr). The equivalent
Stratonovich SDE given as:
8² dx, 2
X,dl+ diag X,) BB
Z X, Tr [r[diag( (x)ddddi))v(aaa))) V dt 2 diag( dBr, d² 2 dy = r rY, X.)
Fully-Coupled FBSDE. Finally, we consider FBSDE with coupled forward and backwarwds dynamics adapted from Bender & Zhang (BZ) [36]:
dx, YYdB,, tE [0,
7 at+ ZJ te (O,T)).
dY; = -ry++ + e 3r(T-() 2
sin(X;;)) 3-n
where X,, denotes the j-th coordinate of X, Rd Due the dependence of the forward process (X₂) on (Y)), this set of coupled FBSDE does not fit into the mathematical formulation set forth in still apply the BSDE methods described at the beginning of Section 6 by initializing Yo = uo (x, 0) and jointly integrating (X,, Y)). We set 7 = 1, r .I, a = .3, D = 1, Xo (x/2, x/2...., x/2). and (XXr) = o(Xr) = D Σÿ= sin(X,,r). The above FBSDE induced from the following PDE
(3.3) and (3.4). Nevertheless,
0ea(2, - ()) +r(,))))))(x))) 1 -x((T-()) 0²
D d sin sin(x;)
ru(xx 2
D
with the analytical solution u(x,t) e-((T-) D D sin(2;)) Additionally, the equivalent Stratonovich SDE is given as:
0²
dX,- ZiY, oY, o
1 - 3r((T))) 2
d
0² (7 Tr Vu, dt 2 dB
dY, -ry, +
Σ sin sin(X
2
100D HJB
100D BSB
Method
float32
float64
float32
float64
PINNs
0.1281 = .0136
0.1281 + .0171
2.9648 = .8652
1.5066 + 2349
FS-PINNs
0.0838 .0170
0.0702 = .0044
0.0602 .0150
0.0497 .0031
EM-BSDE
0.3776 + .0365
0.3820 + .0219
0.2451 # .0160
0.3735 # .4700
EM-BSDE (NR)
0.4459 .0410
0.4676 .0153
0.1648 .0143
0.1855 .0088
Heun-BSDE
0.0675 = .0053
0.0529 + .0029
0.4587 # .0261
0.0535 + .0113
Table 3: float32 vs float64 Performance in 100D HJB/BSB
G.2 Sensitivity to Floating Point Precision
In BSDE-based losses, floating point errois can accumulate through out integration of the SDEs. leading to poor performance on the trained model. As seen in Table 3, the floating point error is especially apparent the Heun loss on the 100D BSB case where the performance improved by a factor of 10 between float32 and float64. In addition, performance improvements were observed for the PINNs and FS-PINNs models as well. models performed slightly worse at f float64, which may be attributed to the bias term present in its loss. Overall, floating point sensitivity is more apparent in the BSB problem than the HJB problem. We attribute this to the non-trivial forward trajectory in the BSB problem. We leave more numerically stable implementations of the Heun solver in float32, such as PDE non-dimensioonalizatoon [7] and the use of the reversible Heun solver from [41], future work.
the model is
also noted that the EM-BSDE
43

---

G.3 Dimensionality Study
Additionally, we demonstrate scalability of the algorithms by re-running the HJB problem at various dimensions and plotting the RL2 enror for each method. Figure 6 shows EM-BSDE underperforming both ES-PIS and Heun-BSDE across all dimensions tested. Additionally, trajectory-based methods scale more effectively to high-dimensional problems with PINNs-based methods showing slight advantage lower dimensions.
PINNS
10
10-1
FS-PINNS
EM-BSDE
EM-BSDE
Heun-BSDE
Heun-BSDE
RL2 Score
RL2
10-2
10-2
10-3
10
100
200 Input Dimension,dir
300
400
500
50
100
150
200
Runtime (min)
Figure plot of RL2 performance for the HJB problem at various dimensions den {2, 10,50, 50, 100, 500)
Figure plot of RL2 performance vs runtime fixed iteration steps
G.4 Improved Algorithm Schemes
We describe the algorithmic differences between the full roll-out algorithm and the batched sub- sampling variation for the general BSDE loss; a similar algorithm is used FS-PINNS, replacing seff-regularization loss with the PINNs loss.
Algorithm 2 Full BSDE Loss Algorithm
Input: Ncural network ûg (z, t), parameters 6, terminal function é, time step At, trajectory length N.
Output: Self-consistnccy loss L.. and terminal loss L.
1: Sample initial state: (a:[0], [{]])=
(.a,0), with
2: Evaluate network at initial state: (u, un) ( (d(xx(]],[()) V.ūo (c())))))
3: Initialize self-consistency losses: Estep
4: for i ,,..,,..11
1 do
5: Sample Brownian noise: N(O, La)
6: /+ Forward
rollout +/ =[1]= Ar
7: Propagate forward state: æ[i
8:
if no resetting then
9:
/* Use either EM or Heun integration (NR) +/
10:
Propagate backward state: yli 1] y|i] Ay
12:
/* Use either EM or Heun integration */ Propagate backward state: yli 1] + Ay
13:
14:
end if
Propagate time: t[i 1]= t[1] At
Evaluate network at new state: (u, uz)) = (u(x(i + 1],4 1]), V)(((( 1],4[i+ 1)))
17:
Record local residual loss: Estep [1] (u y 1])²
18:
end for
Compute self-consistency loss: C.r 5N Estep
20:
Compute terminal loss: Lo (u Φ) |lu, 2
21:
return (Cer, Lo)
As shown in Algorithms 2 and 3, the new variation of the loss separates the forward and backward propagation which enables random sub-sampling the loss evaluation. full sampling (B N), the batched algorithm recovers loss and gradient values consistent with the original algorithm with
44

---

Algorithm 3 Batched, Sub-sampling BSDE Loss Algorithm (Full description of Algorithm 1)
Input: Neural etwork uo( t), parameters 6, terminal function é, time step At, trajectory length N, evaluation batch B.
Output: Self-consistency loss Ler and terminal loss L
Sample initial state: (z[0], [[]]] (Io. 0), with zio μ
2: Sample Brownian noise: { N - -] N(O, Ia)
3: Evaluate network at initial state: (u, u.)
V.
4: /+ Forward SDE rollout +/
5: for 0...., N
6: /* Use either EM or Heun integration */
7: Propagate forward state: xf2 1] x[i] Ar
8: Propagate time: +[i+1] t[1] At
9: if coupled then
Evaluate network new state: u.) (ūo(x[i+ 1],4[i+ 1)), (u(e 1],4[1+ 1))) 11: end
12: end for
13: Stop gradient: [
Separate states: (1,, N1]0::-11,,
Same permutations */
16: Random sub-sampling: (Iz I4+1, 44, fi-1)
Zit (0: B-
17: Evaluate network at points: (ux w:.) (fo (i), (zs,
18: /+ Use either EM Heun Integration */
Compute backward state at batched point: y4-1 = ue + Ay
20: Evaluate network at +1 points: ui+1 ûs (1++1,41-1)
21: /* Use PINNs loss instead for FS-PINNs */
22: Compute self-coosistency loss: Lor EEED' (M++1 -1+1)²
23: Evaluate network at T: (u, = (N, (
24: Compute terminal loss: L. = (u+ o)² ||M V
25: return (Laz, Lo)
proper scaling. Additionally, the stop gradient on the forward SDE has negligible effects on model performance but improve convergence fixes optimization only the backward SDE PDE.
Note that Algorithms 2 and 3 are simplified for readability and excludes some algorithmic details such trajectory batching and loss weighting.
G.5 Behavior Policy Rollouts for HJB Optimal Control
Suppose our control system is deterministic control-affine system: += f(z)+((++)) + (z)u. For a positive definite R, the HJB equation for stagewise cost c(1)) 2 |ul?, and terminal cost cr is:
a,v (V.V. c- 1 vv 00 c(r)))=c((x),
and the optimal control induced by is πy (z, = (x
Now, suppose we have any rollout policy (Ix t), and we consider Itô stochastic rollouts of the form:
dxx" [f((x)) + g(X," xx)))((()))))) t)]dt + odB,
Now, the optimal value function (* 1) satisfies the following SDE:
d² u 1))
dV (X7,))
0.v avv (f(x7) (x()))))+(()))(x))))))))))))(
+o(VV" (X).t).
d² lg ll 1 - - c+ (g VV') o tv(v'v") dt "VV'',,
noting that the last expression is evaluated at (XF,t) t). Hence, the forward/backward Itô SDEs for given value function V and behavior policy и are:
xx; [f(x7) ((X7)=((7 t) odB;,
(G.2a)
45

---

aYp" 1 vv - (gT, vV)+ o² 2 t((p'v) + +
(G.2b)
noting again that the expressions in (G.2b) are evaluated at . t). Similarly, define the forward/hackward Stratonovich SDEs t) adB,
ax;
(G.3a)
dr," V = - 1 lg - (gm, vV) dt + ovv od
(G.3b)
Note that the Itō BSDE (G.2b) requires explicit Hessian computation v²V while the Stratonovich BSDE (G.3b) does not. This holds for all first-order PDEs, such as in deterministic HJB problems.
These forward/backward SDEs can be used in conjunction with the induced policy Ty from the current value function V. Some care must be taken though when setting up the BSDE losses. In particular, since both the forward SDE trajectory (X,"), ) and the (X, .t) terms which appear the backward SDEs depend implicitly on V, stop-gradient operators should be placed so that gradients not pook-propagated through these values, which can destabilize training.
6.6 Pendulum Swing Up Experiment
In addition to the results above, we include a simple pendulum swing-up optimal control experiment inspired by [55]. Given the pendulum equations of motion,
mit (bŌ - mgl sin 0 - w)]
f(z,u)= +
we define optimal control objective:
[' (((() +(x(T)) u(t)
('(o)) min
where:
(x,u)) +(x)++++++)+)))))))++)) ru² ((x) Q sin + 91 (cos 0 - 11) 9202
with a 92, r O. Observe that this problem setup exactly fits the setup in Appendix G.5, and hence both the forward/backward Itō and Stratonovich SDEs in (G.2) and (G.3) directly apply.
G.6.1
Pendulum Results
Metric
PINNs
FS-PINNs
EM-BSDE
Cost
53.17
46.59
46.42
46.43
PDE Error
2.77
3.38
78.94
18.6
Table 4: Accumulated cost and average PDE emror for the pendulum swing-up up problem.
The results of the pendulum swing-up case are outlined in Table 4. We use the specific constants b = 0.1,1 = 1, 9.8, 91 10, 92 1,r in our experiment. It is observed that while the methods remain similar, the lower PDE error on FS-PINNs and Heun-BSDE signify better learned solutions. Furthermore, in Figure that Heun-BSDE generally has the lowest PDE error with high errors only at the discontinuities.
accumulated cost between the three trajectory-
observe
46

---

PINNS
FS-PINNS
100
-6
101
-75
-5.0
-2.5
2.5
5.0
7.5
-8
-7.5
-5.0
-2.5
0.0
2.5
5.0
7.5
(log
EM-BSDE
Heun-BSDE
W
10t
10<1
-7.
-5.0
-2.5
0.0
2.5
5.0
7.5
-75
-5.0
-25
0.0
2.5
5.0
7.5
Figure 8: PDB error 0 for the pendulum swing up case.
47


Reviews:
[{"review_id": "vjuFMJAyl2", "reviewer": "Reviewer_LXC6", "summary": "This paper studies backward stochastic differential equation (BSDE) based deep learning methods for solving high-dimensional partial differential equations (PDEs). These methods provide an alternative to the popular PINNs and have an advantage that the underlying equation do not need to be known, but their performance has been lagging behind PINNs. The authors identify the reason of this performance gap as the Euler-Maruyama (EM) integrator typically used for BSDE methods, which is based on ito SDE formulation and is shown by the authors to induce significant discretization bias. The authors propose to employ Stratonovich SDE interpretation and use Heun integration scheme, which eliminates the raised issue. The authors empirically demonstrate that the fixed BSDE algorithm is competitive to PINNs in solving three PDE benchmark problems.", "questions": "I have no particular questions, but would like to hear the authors' opinion on W1.", "limitations": "The authors discussed limitations in Sections 6, 7 and Appendix A.", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths\n\nS1. The paper addresses an important problem of engineering alternative approaches to PINNs for PDE solving which have qualitatively different characteristics (e.g., in the required information for training).\n\nS2. The authors theoretically show that BSDE methods with EM integrator suffers from a bias term which can be problematic even in the limit of infinite training data, and demonstrate it through an informative example (Theorem 4.1; Figure 1). The algorithmic fix based on Stratonovich BSDEs and Heun integration are well-motivated and analyzed to avoid the discretization error (Theorem 4.2; Figure 1). A potential fix based on multi-step losses is further analyzed in Section 5.\n\nS3. The experimental results in Section 6 support the theoretical claims made in Sections 4 and 5. The baselines are properly chosen, the results are presented clearly, and limitations are transparently discussed.\n\nWeaknesses\n\nMy concerns are mostly about the the limitations already discussed in the paper by the authors.\n\nW1. As the authors discussed in Section 7, the fact that the proposed method requires computing the Hessian leads to a large computational overhead, including the requirement to use float64 precision. While this is properly acknowledged and discussed (including Appendix E.2), it would be informative to, e.g., have a controlled EM-BSDE baseline which is allowed to use a similar amount of computation for integration, if it is possible.\n\nMinor comment\n- Line 35: involving -> involves?", "was_revised": true, "final_justification": "The authors have clearly addressed the original main concern on computational cost through a careful algorithmic optimization.", "rebuttal": "We thank the reviewer for their thoughtful and constructive review. We greatly appreciate the recognition of the importance of studying alternatives to PINNs approaches for solving PDEs, in addition to the well-motivated approach that we take in this work.\n\n**Response to W1 (Computational Costs):**\n\n*(Note that similar questions were also brought up by Reviewers Fvez and sGAy, so we have duplicated our response to make the rebuttals self-contained).*\n\nWe agree that the computational costs of the Heun integration presented in the original draft are key limitations for practical implementation. Hence, we spent a bit of time investigating how we could improve our implementation. **We are happy to report that we have improved our implementation—for all methods—and the result is that the computational overhead of the Heun implementation is much more reasonable (e.g., <2x overhead over PINNs/FS-PINNs)**. Due to the given rebuttal time constraints, we were unfortunately only able to re-run our new implementation for all methods on the 100D HJB case, but we expect the computational benefits from the implementation changes we made (which we describe below) to carry over to all PDE cases. We will of course re-run all methods over all PDE cases for the revised draft.\n\nWe now describe the implementation changes we made. In the original Heun-BSDE implementation, the forward and backward SDEs rollouts were performed simultaneously, evaluating $L_{Heun_{k,\\tau}}(\\theta)$ at each step. As a result, the loss gradient calculation required back-propagation through the entire rollout; we had to resort to gradient checkpointing when combined with the Hessian calculation of $u$ in order to avoid running out of memory. In our new implementation, we first perform a forward SDE rollout to generate trajectories, $\\mathcal{X}=\\lbrace (x^i_t,t)\\mid t \\in [0,T],i=1,\\dots,N \\rbrace$, then calculate the loss in one pass by batching across $\\mathcal{X}$. Additionally, we further optimize and calculate the loss on $\\mathcal{B}\\subset \\mathcal{X}$, which importantly decouples the batch size from trajectory length. Both EM-BSDE and FS-PINNs methods can be (and were) optimized in the same manner, with the exception of the no-resetting method (EM-BSDE (NR)) (see Lines 229-232) which requires additional considerations due to its propagation of $Y_t$. In the following table, we report the average performance and overhead across three runs for each method with a batch size of 256:\n\n| **100D HJB** | **RL2** | **Runtime (s)** | **Overhead** |\n|-|-|-|-|\n| **PINNs** | 0.1349 ± .0368 | 2170 ± 9 | 1x |\n| **FS-PINNs** | 0.0611 ± .0137 | 2227 ± 4 | 1.03x |\n| **EM-BSDE** | 0.3549 ± .0119 | 384 ± 2 | .18x |\n| **Heun-BSDE** | 0.0464 ± .0050 | 3846 ± 78 | 1.77x |\n\nThese results align closely with the RL2 errors reported in Table 1, albeit with **significantly lower** runtimes (cf. Table 4). We also note some key details:\n- The updated trajectory rollout schemes for FS-PINNs and EM-BSDE are more efficient with nearly zero overhead. This is reflected in the runtime of FS-PINNs vs. PINNs.\n- The Heun-BSDE scheme requires approximately double the compute for Heun-based trajectory rollouts and loss calculation versus FS-PINNs, which is reflected in its overhead.\n- For the specific elliptic/parabolic PDE we consider in (3.1), EM-BSDE does not require the computationally expensive Hessian computation of $u$ in the backwards SDE, while Heun-BSDE (and (FS)-PINNs) does, leading to significantly faster runtimes. This does not necessarily hold true for all PDEs however. In general, $h[u]$ can also depend on $\\nabla^2 u$ as noted in e.g., [Pham et al., 19] referenced by Reviewer DQqj, allowing for fully nonlinear second order PDEs. Alternatively, when solving first-order PDEs (e.g. deterministic HJB problems), the Stratonovich formulation $$ \\mathrm{d}X_t=f(X_t,t)\\mathrm{d}t+\\sigma \\circ \\mathrm{d} B_t $$ $$ \\mathrm{d}Y_t=[\\partial_tu(X_t,t)+\\nabla_xu(X_t,t)\\cdot f(X_t,t)]\\mathrm{d} t + Z_t^T\\circ \\mathrm{d}B_t $$ does not introduce a Hessian term, \nwhile the Ito formulation $$\\mathrm{d}X_t=f(X_t,t)\\mathrm{d}t+\\sigma \\mathrm{d} B_t$$\n$$\\mathrm{d}Y_t=\\left[\\partial_tu(X_t,t)+\\nabla_xu(X_t,t)\\cdot f(X_t,t)+\\frac{\\sigma^2}{2}\\mathrm{Tr}\\left[\\nabla^2_{xx}u(X_t,t)\\right]\\right]\\mathrm{d} t + Z_t^T \\mathrm{d}B_t$$ requires the Hessian term for calculation. In these cases, we actually expect Heun-BSDE to either match or be faster than EM-BSDE.\n\nIn addition, we ran the original EM-BSDE solver to equal compute time as the original Heun-BSDE method on the 100D HJB case, and as predicted by our theoretical analysis, we found that EM-BSDE does not improve significantly with more iterations (see table).\n\n| **Iterations** | **Run Time** | **RL2** |\n|-|-|-|\n| 100,000 | .61 hours | 0.3626 ± .0113 |\n| 2,000,000 | 14.5 hours | 0.3394 ± .0093 |\n\nWhile EM-BSDE converges faster than all other methods, the **RL2 error quickly plateaus**, leading to little improvement despite further optimization.\n\nOverall, the updated results align with the paper’s original claims while **demonstrating the potential for practical application**. In the final revision, we will update all relevant results with the updated implementation discussed above and further include our insights on Hessian computations. Additionally, we will include runtime-normalized performance plots (e.g. error vs wall-clock time) with each method for added clarity and completeness."}, {"review_id": "IGBzk1Y9eH", "reviewer": "Reviewer_Fvez", "summary": "Backward SDE-based deep learning methods offer a powerful framework for solving high-dimensional PDEs. Most existing approaches rely on the Euler–Maruyama (EM) scheme for time integration. This paper revisits that choice carefully and points out a discretization bias introduced by the EM scheme. To address this, the authors propose a Stratonovich-based BSDE formulation paired with a stochastic Heun integrator, which eliminates the bias. The paper presents detailed theoretical analysis of both schemes and supports the findings with three high-dimensional PDE examples, demonstrating notable improvements with the proposed method.", "questions": "Beyond the points raised above, here are several additional comments:\n\n1. Literature discussion contains several inaccuracies:\n\n    - The claim in the abstract that \"existing BSDE-based solvers have empirically been shown to underperform relative to PINNs in the literature\" is quite misleading. In the main text, this is attributed to Ref [6], but the diffusion loss proposed in that work itself is based on a BSDE formulation and shows notably better accuracy than both PINNs and the original Deep BSDE method. So this statement undercuts the actual strengths shown by BSDE-based solvers in the high-dimensional PDE setting relative to PINNs; see also the recent development in arXiv:2409.08526.\n\n    - In the introduction, the authors cite Refs [1,2]—the original PINN papers—as “scalable alternatives” meant to address the curse of dimensionality. This characterization is inaccurate. Those early PINN works (and many follow-ups) primarily focused on low-dimensional problems and were not designed with high-dimensional PDEs in mind. In contrast, BSDE-based methods were motivated from the start by high-dimensional settings, as discussed in Refs [4,22]. To the best of reviewer's knowledge, interest in applying PINNs to high-dimensional PDEs only developed more recently; see, e.g., arXiv:2307.12306 and references therein.\n\n    - Ref [3], the Deep Ritz method, is not a PINN-type approach; it is a distinct variational method and should not be grouped together with PINNs.\n\n1. Table 2 presents computational time differences across methods. It would be helpful to also compare performance under equal computation time rather than fixed iteration counts, as that would give a more meaningful picture of practical efficiency. Again, including diffusion loss from Ref [6] as a baseline is helpful.\n\n1. There is a substantial body of theoretical work on Itô-based BSDE formulations that supports the analysis of deep learning–based methods (e.g., arXiv:1811.01165, arXiv:2307.15496). It would be helpful if the authors could comment on whether a similar theoretical framework could be developed for the Stratonovich-based formulation. Such discussion could strengthen the theoretical appeal of the proposed method.", "limitations": "As acknowledged, the overhead due to second-order derivatives is a main limitation.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 4, "strengths_and_weaknesses": "Strengths:\n\n1. The paper challenges a commonly used but under-examined component—time integration—in BSDE-based PDE solvers and highlights the significant role that the choice of integration scheme plays in solution quality. It is of great originality.\n\n2. The authors conduct a careful and thorough numerical analysis of discretization errors associated with both schemes. The implications are consistent with prior observations in the literature and are well supported by numerical experiments.\n\nWeaknesses:\n\n1. As the authors themselves note, the method proposed in Ref [6] is closely related but is not included as a baseline. This seems like a notable omission.\n\n1. The additional computational cost introduced by the second-order derivatives required in the Stratonovich formulation appears to be a key limitation of the method. Further comparisons that account for this overhead would provide a clearer assessment of the method’s practical value.", "was_revised": true, "final_justification": "The paper is technically sound and is likely to have a significant impact on the actively studied field of solving PDEs and related problems through an SDE formulation.", "rebuttal": "We thank the reviewer for their thoughtful and constructive review. We greatly appreciate the recognition of both the originality of our work and its careful and thorough numerical analysis.\n\n**Response to W1/Q2 (Comparisons to [6]):**\n\nThank you for raising this point. While we do not specifically use the term “diffusion loss” from [6] in the baselines of our experiments, the EM-BSDE implementation of the multi-step self consistency losses $L_{BSDE,H}(\\theta)$ evaluated in Section 6.2 is precisely the method described in [6]; we discuss this in both Footnote 2 and Lines 205-207. In our revision, we will improve the clarity of the writing in the paper describing the connection between the two losses. Given this connection, Figure 4 demonstrates exactly the comparison that the reviewer is inquiring about: tuning the skip length does improve performance consistent with the empirical results in [6], but ultimately does not solve the bias issue in EM-BSDE (cf. Lines 303-305). Again, we will make this connection more explicitly clear in our revision.\n\n**Response to W2/Q2 (Computational Costs):**\n\n*(Note that similar questions were also brought up by Reviewers LXC6 and sGAy, so we have duplicated our response to make the rebuttals self-contained).*\n\nWe agree that the computational costs of the Heun integration presented in the original draft are key limitations for practical implementation. Hence, we spent a bit of time investigating how we could improve our implementation. **We are happy to report that we have improved our implementation—for all methods—and the result is that the computational overhead of the Heun implementation is much more reasonable (e.g., <2x overhead over PINNs/FS-PINNs)**. Due to the given rebuttal time constraints, we were unfortunately only able to re-run our new implementation for all methods on the 100D HJB case, but we expect the computational benefits from the implementation changes we made (which we describe below) to carry over to all PDE cases. We will of course re-run all methods over all PDE cases for the revised draft.\n\nWe now describe the implementation changes we made. In the original Heun-BSDE implementation, the forward and backward SDEs rollouts were performed simultaneously, evaluating $L_{Heun_{k,\\tau}}(\\theta)$ at each step. As a result, the loss gradient calculation required back-propagation through the entire rollout; we had to resort to gradient checkpointing when combined with the Hessian calculation of $u$ in order to avoid running out of memory. In our new implementation, we first perform a forward SDE rollout to generate trajectories, $\\mathcal{X}=\\lbrace (x^i_t,t)\\mid t \\in [0,T],i=1,\\dots,N \\rbrace$, then calculate the loss in one pass by batching across $\\mathcal{X}$. Additionally, we further optimize and calculate the loss on $\\mathcal{B}\\subset \\mathcal{X}$, which importantly decouples the batch size from trajectory length. Both EM-BSDE and FS-PINNs methods can be (and were) optimized in the same manner, with the exception of the no-resetting method (EM-BSDE (NR)) (see Lines 229-232) which requires additional considerations due to its propagation of $Y_t$. In the following table, we report the average performance and overhead across three runs for each method with a batch size of 256:\n\n| **100D HJB** | **RL2** | **Runtime (s)** | **Overhead** |\n|-|-|-|-|\n| **PINNs** | 0.1349 ± .0368 | 2170 ± 9 | 1x |\n| **FS-PINNs** | 0.0611 ± .0137 | 2227 ± 4 | 1.03x |\n| **EM-BSDE** | 0.3549 ± .0119 | 384 ± 2 | .18x |\n| **Heun-BSDE** | 0.0464 ± .0050 | 3846 ± 78 | 1.77x |\n\nThese results align closely with the RL2 errors reported in Table 1, albeit with **significantly lower** runtimes (cf. Table 4). We also note some key details:\n- The updated trajectory rollout schemes for FS-PINNs and EM-BSDE are more efficient with nearly zero overhead. This is reflected in the runtime of FS-PINNs vs. PINNs.\n- The Heun-BSDE scheme requires approximately double the compute for Heun-based trajectory rollouts and loss calculation versus FS-PINNs, which is reflected in its overhead.\n- For the specific elliptic/parabolic PDE we consider in (3.1), EM-BSDE does not require the computationally expensive Hessian computation of $u$ in the backwards SDE, while Heun-BSDE (and (FS)-PINNs) does, leading to significantly faster runtimes. This does not necessarily hold true for all PDEs however. In general, $h[u]$ can also depend on $\\nabla^2 u$ as noted in e.g., [Pham et al., 19] referenced by Reviewer DQqj, allowing for fully nonlinear second order PDEs. Alternatively, when solving first-order PDEs (e.g. deterministic HJB problems), the Stratonovich formulation $$ \\mathrm{d}X_t=f(X_t,t)\\mathrm{d}t+\\sigma \\circ \\mathrm{d} B_t $$ $$ \\mathrm{d}Y_t=[\\partial_tu(X_t,t)+\\nabla_xu(X_t,t)\\cdot f(X_t,t)]\\mathrm{d} t + Z_t^T\\circ \\mathrm{d}B_t $$ does not introduce a Hessian term, \nwhile the Ito formulation $$\\mathrm{d}X_t=f(X_t,t)\\mathrm{d}t+\\sigma \\mathrm{d} B_t$$\n$$\\mathrm{d}Y_t=\\left[\\partial_tu(X_t,t)+\\nabla_xu(X_t,t)\\cdot f(X_t,t)+\\frac{\\sigma^2}{2}\\mathrm{Tr}\\left[\\nabla^2_{xx}u(X_t,t)\\right]\\right]\\mathrm{d} t + Z_t^T \\mathrm{d}B_t$$ requires the Hessian term for calculation. In these cases, we actually expect Heun-BSDE to either match or be faster than EM-BSDE.\n\nIn addition, we ran the original EM-BSDE solver to equal compute time as the original Heun-BSDE method on the 100D HJB case, and as predicted by our theoretical analysis, we found that EM-BSDE does not improve significantly with more iterations (see table).\n\n| **Iterations** | **Run Time** | **RL2** |\n|-|-|-|\n| 100,000 | .61 hours | 0.3626 ± .0113 |\n| 2,000,000 | 14.5 hours | 0.3394 ± .0093 |\n\nWhile EM-BSDE converges faster than all other methods, the **RL2 error quickly plateaus**, leading to little improvement despite further optimization.\n\nOverall, the updated results align with the paper’s original claims while **demonstrating the potential for practical application**. In the final revision, we will update all relevant results with the updated implementation discussed above and further include our insights on Hessian computations. Additionally, we will include runtime-normalized performance plots (e.g. error vs wall-clock time) with each method for added clarity and completeness.\n\n**Response to Q1 (Literature review):** \n\nWe thank the reviewer for the helpful clarifications regarding literature discussions. We agree that the abstract and introduction could be more precise in the positioning of the PINNs and BSDE literature. Specifically,\n- We acknowledge that as noted by the reviewer, [6] demonstrates that BSDE solvers can exceed the performance of standard PINNs solvers in high dimensional settings, with the use of e.g., the proposed diffusion loss. However, as shown in [6], standard (EM-)BSDE methods do indeed struggle compared with PINNs; our statement about existing (EM-)BSDE solvers in the abstract was intended to refer to this phenomenon. It is this gap in performance, in addition to trying to answer why the diffusion loss empirically improves performance, that motivated the present analysis. We apologize for the confusion, and in the final version we will revise the abstract and introduction to clarify this point.\n- We agree that the characterization of early PINNs as “scalable alternatives\" for high-dimensional PDEs should be clarified. As pointed out by the reviewer, the original PINNs works [1,2] cited primarily focused on low-dimensional PDEs. However, we do believe that the motivation for using PINNs like methods for solving high-dimensional PDEs did arise quite early in the development. For example, the initial work on Deep Galerkin Methods (DGMs) from [Sirignano and Spiliopoulos, 2018] was motivated by looking for meshfree algorithms for high-dimensional PDEs. Furthermore, as another example, practitioners from the control theory community started to use PINNs methods to solve Hamilton-Jacobi reachability PDEs [Bansal and Tomlin, 2020] not long after the DGM work, with the explicit motivation to avoid the curse of dimensionality from classic grid-based PDE solvers. Nevertheless, we do agree that the work of [Hu et al., 2023] does seem to be the first work that explicitly studies the extent to which PINNs methods can be scaled to high dimension (e.g., 100k dimensions) through the use of algorithmic modifications to the training. We will certainly update our related work to reflect this chronology more accurately.\n- We agree that [3] is a distinct method from PINNs and should not be grouped under the PINNs umbrella– thanks for catching this mistake. We will correct this in our revision.\n\n**Response to Q3 (Stratonovich-based BSDE theory for deep-learning based methods):**\n\nThis is an interesting question, and we believe the answer to be yes. We answer the question in the context of the self-consistency losses we consider, but similar answers should also hold for the recursive formulations (where separate models are estimated at every timestep, see Line 116 for the distinction). Our result Theorem (4.2) combined with Eq. (4.8) is essentially the starting point for such a Stratonovich theory, by showing that the error of the $L_{Heun,\\tau}(\\theta)$ loss scales as the PDE residual term $R[u_\\theta]$ (integrated along trajectories of $X_t$) plus an $O(\\tau^{1/2})$ term that arises from the discretization. Hence, by choosing the discretization small enough, the difference between the Heun loss and the PDE residual integral can be made small enough. What is left to derive is that having small PDE residual error implies control of e.g., the true backwards SDE paths $(Y_t)$ and approximate path $(\\hat{Y}^\\theta_{nk})$ generated by the model $u_\\theta$, or other measures of quantify solution quality that is useful for downstream applications. This part, while technical in nature, should at a high level follow similar arguments as the cited papers for the Ito setting."}, {"review_id": "kUbSnh3NX3", "reviewer": "Reviewer_sGAy", "summary": "The authors introduce a new integrator based on Stratonovich viewpoint of SDEs to reduce biases induced with Euler based solvers. The authors consider this in the case of solving backwards SDEs for use in semi-linear parabolic PDEs with applications in control. They illustrate theoretically and numerically the advantages of the Stratonovich based BSDE integrator.", "questions": "It’s surprising to see how well FS-PINNs works in the experiments, do the authors have an idea as to why this is the case?\n\nDo the authors have any suggestions for how they can improve upon the computation speed of the proposed method?\n\nUnder a similar computational budget, how do the different methods compare? E.g. for a time constrained problem, which would be the correct method to use?\n\nI’m a bit curious as to why PINNs are much faster, my initial impression is that computing derivatives manually especially in high dimensions would be slower? Is this mainly because these are vector hessian products?\n\nI think a study on dimensionality vs performance could also be quite a bit helpful. \n\nFrom table 3, is there any intuition why the 32 bit performance makes the Euler Maruyama baselines perform better than 64 bit?", "limitations": "Yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\nThe proposed integrator has theoretical and empirical advantages over existing Euler Maruyama and PINNs based approaches for solving high dimensional semi linear parabolic PDEs.\n\nThe problem is studied from a numerical analysis perspective, which provides useful results for other fields beyond machine learning. \n\nWeaknesses:\n\nThe main weakness is the induced computation time from the proposed method. \n\nThe empirical results could probably use a bit more analysis. For example, in one case the vanilla EM method performs better which is quite surprising given that there’s only four examples of empirical studies conducted on the methods. \n\nThis applies to a limited set of PDEs and perhaps some more examples should be motivated within the text.", "was_revised": true, "final_justification": "While I still believe that the authors should do more numerical studies and should adequately describe what performance improvements were done, I think the paper could be an interesting contribution to the field (while echoing some of the other reviewers' comments that more motivation should be added within the manuscript to relate this class of PDEs to more well known concepts within the field).", "rebuttal": "We thank the reviewer for their thoughtful and constructive review. We greatly appreciate the recognition of the usefulness of our numerical analysis perspective, with impact possibly beyond the ML community.\n\n**Response to W1/Q2/Q3 (Computational Costs):**\n\n*(Note that similar questions were also brought up by Reviewers LXC6 and Fvez, so we have duplicated our response to make the rebuttals self-contained).*\n\nWe agree that the computational costs of the Heun integration presented in the original draft are key limitations for practical implementation. Hence, we spent a bit of time investigating how we could improve our implementation. **We are happy to report that we have improved our implementation—for all methods—and the result is that the computational overhead of the Heun implementation is much more reasonable**. Due to the given rebuttal time constraints, we were unfortunately only able to re-run our new implementation for all methods on the 100D HJB case, but we expect the computational benefits from the implementation changes we made (which we describe below) to carry over to all PDE cases. We will of course re-run all methods over all PDE cases for the revised draft.\n\nIn the original Heun-BSDE implementation, the forward and backward SDEs rollouts were performed simultaneously, evaluating $L_{Heun_{k,\\tau}}(\\theta)$ at each step. As a result, the loss gradient calculation required back-propagation through the entire rollout; we had to resort to gradient checkpointing when combined with the Hessian calculation of $u$ in order to avoid running out of memory. In our new implementation, we first perform a forward SDE rollout to generate trajectories, $\\mathcal{X}=\\lbrace(x^i_t,t)\\mid t\\in[0,T],i=1,\\dots,N\\rbrace$, then calculate the loss in one pass by batching across $\\mathcal{X}$. Additionally, we further optimize and calculate the loss on $\\mathcal{B}\\subset\\mathcal{X}$, which importantly decouples the batch size from trajectory length. Both EM-BSDE and FS-PINNs methods can be (and were) optimized in the same manner, with the exception of the no-resetting method (EM-BSDE (NR)) (see Lines 229-232) which requires additional considerations due to its propagation of $Y_t$. In the following table, we report the average performance and overhead across three runs for each method with a batch size of 256:\n\n|**100D HJB**|**RL2**|**Runtime (s)**|**Overhead**|\n|-|-|-|-|\n|**PINNs**|0.1349 ± .0368|2170 ± 9|1x|\n|**FS-PINNs**|0.0611 ± .0137|2227 ± 4|1.03x|\n|**EM-BSDE**|0.3549 ± .0119|384 ± 2|.18x|\n|**Heun-BSDE**|0.0464 ± .0050|3846 ± 78|1.77x|\n\nThese results align closely with the RL2 errors reported in Table 1, albeit with **significantly lower** runtimes (cf. Table 4). We also note some key details:\n- The updated trajectory rollout schemes for FS-PINNs and EM-BSDE are more efficient with nearly zero overhead. This is reflected in the runtime of FS-PINNs vs. PINNs.\n- The Heun-BSDE scheme requires approximately double the compute for Heun-based trajectory rollouts and loss calculation versus FS-PINNs, which is reflected in its overhead.\n- For the specific elliptic/parabolic PDE we consider in (3.1), EM-BSDE does not require the computationally expensive Hessian computation of $u$ in the backwards SDE, while Heun-BSDE (and (FS)-PINNs) does, leading to significantly faster runtimes. This does not necessarily hold true for all PDEs however. In general, $h[u]$ can also depend on $\\nabla^2u$ as noted in e.g., [Pham et al., 19] referenced by Reviewer DQqj, allowing for fully nonlinear second order PDEs. Alternatively, when solving first-order PDEs, the Stratonovich formulation $$\\mathrm{d}X_t=f(X_t,t)\\mathrm{d}t+\\sigma\\circ\\mathrm{d}B_t$$ $$\\mathrm{d}Y_t=[\\partial_tu(X_t,t)+\\nabla_xu(X_t,t)\\cdot f(X_t,t)]\\mathrm{d}t+Z_t^T\\circ\\mathrm{d}B_t$$ does not introduce a Hessian term, \nwhile the Ito formulation $$\\mathrm{d}X_t=f(X_t,t)\\mathrm{d}t+\\sigma\\mathrm{d}B_t$$\n$$\\mathrm{d}Y_t=\\left[\\partial_tu(X_t,t)+\\nabla_xu(X_t,t)\\cdot f(X_t,t)+\\frac{\\sigma^2}{2}\\mathrm{Tr}\\left[\\nabla^2_{xx}u(X_t,t)\\right]\\right]\\mathrm{d}t+Z_t^T\\mathrm{d}B_t$$ requires the Hessian term for calculation. In these cases, we actually expect Heun-BSDE to either match or be faster than EM-BSDE.\n\nIn addition, we ran the original EM-BSDE solver to equal compute time as the original Heun-BSDE method on the 100D HJB case, and as predicted by our theoretical analysis, we found that EM-BSDE does not improve significantly with more iterations (see table).\n\n|**Iterations**|**run time**|**RL2**|\n|-|-|-|\n|100,000|.61 hours|0.3626 ± .0113|\n|2,000,000|14.5 hours|0.3394 ± .0093|\n\nWhile EM-BSDE converges faster than all other methods, the **RL2 error quickly plateaus**, leading to little improvement despite further optimization.\n\nIn the final revision, we will update all relevant results with the updated implementation discussed above and further include our insights on Hessian computations. Additionally, we agree that a study on dimensionality and run-time normalized results would provide valuable insights for the paper. Therefore, we plan to include results from both studies as plots for added clarity and completeness.\n\n**Response to W2/W3 (Empirical Analysis and Generality of PDEs considered):**\n\nFirst, we hypothesize that EM-BSDE vs Heun-BSDE performance on the 100D BZ case in Table 1 is due to the optimization challenges which arise from the high-dimensional, coupled nature of the problem, which may be confounding performance results (cf. Lines 252-256). This is evidenced by the lack of convergence for PINNs, FS-PINNs, and EM-BSDE (NR), in addition to recovering the expected ordering at lower dimensions (10D BZ). \n\n*(Note that the following response’s question was also brought up by Reviewer DQqj, so we have duplicated our response to make the rebuttals self-contained).*\n\nAdditionally, the non-linear PDEs we consider (cf. Eq. (3.1)) are the standard kind of PDEs considered for BSDE methods (see e.g., references [4-6]). There are many types of applications, including in quantitative finance (e.g., multi-asset options pricing via Black-Scholes, local/stochastic volatility extensions from the work of Duprie and Heston), optimal control (e.g., Hamilton-Jacobi-Bellman equations for finding value functions, Hamilton-Jacobi-Issacs for reach-avoid problems), reaction-diffusion in chemistry (e.g., Allen-Cahn and Fisher-KPP equations), committer/exit-time probabilities for transition path analysis (e.g. [R1]), and so on. We will include a paragraph about these applications to our revised draft.\n\nWe also note that another reason we consider the specific elliptic/parabolic form (3.1) given in the paper (i.e., $h[u]$ does not depend on $\\nabla^2u$) is to favor the standard EM-BSDE formulation, since the form (3.1) means that the EM-BSDE approach does not need to compute Hessians of $u$ in the backwards SDE, whereas the Heun-BSDE does.\n\nThe empirical examples tested in the paper were deliberately chosen for their simplicity and existence of analytical solutions in order to focus on isolating the integration bias introduced by EM-BSDE. For additional examples of PDEs, we are currently working to apply these methods to more complex non-linear optimal control problems (e.g., a high-dimensional task involving collision avoidance for many Dubin’s cars), and will report on our results in the final paper.\n\n[R1]: E, Vanden‑Eijnden (2006), Towards a Theory of Transition Paths. J. Stat. Phys.\n\n**Response to Q1 (FS-PINNs Performance):**\n\nThe FS-PINNs performance is in some sense not surprising, since the RL2 metric we evaluate performance over uses the forward SDE trajectories to compute the collocation points (cf. Lines 239-240, Eq. (6.2)). Hence, the FS-PINNs sampling distribution in training matches the performance metric better versus the PINNs sampling distribution in training, which uses a normal distribution fit to match the spatial distribution of the forward SDE trajectories. This explains the difference in performance.\n\nWe note that our choice of RL2 metric is deliberate—using the forward SDE as the sampling distribution is a much more meaningful measure than using e.g., uniform or isotropic Gaussian over a high-dimensional spatial domain. Furthermore, for problems such as solving Hamilton-Jacobi-Bellman (HJB) equations, the forward SDE trajectories actually have semantic meaning for the problem (e.g., either open-loop or closed-loop trajectory rollouts for HJB), which is substantially more useful for downstream tasks than the typical uniform distributions used in standard PINNs approaches. \n\n**Response to Q4 (PINNs Speed):**\n\nWith the exception of EM-BSDE, all models require Hessian computations in the tested cases. The perceived speed advantage of PINNs stems from the other methods requiring trajectory rollouts for sampling, which impacted speed in our original implementation. We address this issue in our new implementation described in (Computational Costs).\n\n**Response to Q6 (f32 vs f64):**\n\nWe observed this as well and believe that one plausible explanation is that float32 works as a form of implicit regularization due to the stochasticity introduced by lower precision (but importantly, the benefits are not enough to overcome the issues with EM-BSDE vs Heun-BSDE, cf. Table 3). Equally likely, there could also be subtle differences in how either jax, the Cuda layer, or the underlying GPU hardware handles f32 vs f64 which causes the discrepancy—these types of numerical issues are unfortunately not uncommon in modern NN training pipelines. To further investigate, we ran a similar test for both EM-BSDE variants on the 100D HJB benchmark which showed minimal differences between f32 and f64 (see table below, cf. Table 3). Interestingly however, there does still remain a persistent gap between f32 and f64, which we leave to future investigation. \n\n|**100D HJB**|**`float32`**|**`float64`**|\n|-|-|-|\n|EM-BSDE|0.3776 ± .0365|0.3820 ± .0219|\n|EM-BSDE (NR)|0.4459 ± .0410|0.4676 ± .0153|"}, {"review_id": "81ZgWBhGwz", "reviewer": "Reviewer_DQqj", "summary": "The article is concerned with deep learning based methods for PDEs and more specifically, the deep BSDE method based on a stochastic reformulation of the PDE. Authors study specific loss function, so-called self-consistent loss, with Euler-Maruyama scheme for discretization of underlying stochastic processes. Paper analyses mathematically the limiting behaviour of loss function as the discretization step-size tends to zero and proves that a bias arises. Then, authors introduce an alternative discretization based on Stratonovich formulation of underlying stochastic differential equations and prove that in this case no bias arises.", "questions": "-\tCan you spell out formally the assumptions that guarantee that (4.3) holds? Moreover, what are the conditions for strong convergence of the EM scheme? Can you provide a reference where this statement can be found?\n-\tCan you provide a summary of all the assumptions needed to make the overall conclusion of the paper valid?\n-\tCan you specify what precise assumptions are needed to make lines 903, 905, 909 valid?\n-\tCan you provide more background on the considered PDEs? Why are these relevant? In which applications does each of them arise? \n-\tCan you discuss how your method compares to state-of-the art methods based on stochastic representations?", "limitations": "Yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "### Strengths\n\n-\tThe paper raises an important point: PINNs and deepBSDE type methods have not been systematically compared in the literature\n- The paper is mostly well-written and presented\n-\tThe observation that the Euler scheme might create a bias for the deep BSDE method is interesting and new\n-\tThe numerical experiments appear to be sound and indicate that the new method performs well, at the expense of higher computational cost\n- The method appears to be new and innovative\n\n\n### Weaknesses\n\n**Imprecise mathematical statements:** My main concern with the paper is its lack of mathematical rigour. Broadly speaking, the paper is written in a way that indicates that probably there are conditions under which the conclusion holds true. However, in the present form of the paper it is impossible to see clearly what these conditions are. To provide a few examples:\n\n-\tTo derive (4.3), which is crucial for the main contribution of the paper, the authors do not spell out formally any assumptions, but just mention *basic regularity assumptions on $R$, $H$, and $u_\\theta$* and *strong convergence of the EM scheme*. What precisely are these conditions? What are the conditions for strong convergence of the EM scheme? Can you provide a reference where this statement can be found?\n-\tAlso in the proofs, there are several statements that are too informal to be checked properly. For example, in line 903: *assuming enough regularity on R*. Again, can you specify what are the precise assumptions here? Similarly, line 905: *where the penultimate equality holds from the standard approximation properties of a left Riemann sum to the Riemann integral*. Can you be more precise here? How do you get the desired decay with $\\tau$? \nMoreover, some assumptions are just swept under the carpet, see line 909: *We assume that sufficient regularity for the order 1/2 strong convergence of EM jointly on …* This seems to be a strong assumption, that is not mentioned in the main text. Here I feel like the presentation would need to be significantly improved to clarify what are the assumptions under which the claimed statements hold.  \n\n**Background on PDEs limited:** The paper does not provide context on the considered class of PDEs. For example, deep BSDE methods apply to specific types of parabolic PDEs (so-called Kolmogorov PDEs). Why are the considered PDEs relevant? In what applications do they arise? Which important PDEs are of this type?\n\n**Benchmarks not reflecting state-of-the-art** The overall motivation is based on numerical experiments from reference [6], which is from 2021 and does not reflect the state-of-the-art. In [6] it appears that PINNs outperform BSDE methods, but in the meantime, several other papers improving upon the deep BSDE method have been employed. To justify the strong conclusions made in the paper, some of these methods would need to be considered as well. Examples:\n\n- Beck, Christian and Becker, Sebastian and Cheridito, Patrick and Jentzen, Arnulf and Neufeld, Ariel. Deep Splitting Method for Parabolic PDEs. SIAM Journal on Scientific Computing.\n- Riu Naito, Toshihiro Yamada. An acceleration scheme for deep learning-based BSDE solver using weak expansions.  International Journal of Financial Engineering 2020 07:02\n- Akihiko Takahashi a, Yoshifumi Tsuchida, Toshihiro Yamada . A new efficient approximation scheme for solving high-dimensional semilinear PDEs: Control variate method for Deep BSDE solver. Journal of Computational Physics, 2022. \n- Jean-François Chassagneux, Junchao Chen, Noufel Frikha. Deep Runge-Kutta schemes for BSDEs, arXiv Preprint 2212.14372. \n- Huyên Pham, Xavier Warin & Maximilien Germain. Neural networks-based backward scheme for fully nonlinear PDEs, SN Partial Differential Equations and Applications.", "was_revised": true, "final_justification": "Overall, the paper is well-written and mathematically rigorous. It provides an analysis of the deep BSDE method and shows how the Euler scheme might create a bias for the deep BSDE method. The numerical experiments indicate that the new method performs well. \nIn the rebuttal, the authors have addressed my concerns about mathematical rigour and have suggested a new theorem to summarize the main findings. This will greatly clarify the presentation. Thus, in my view the only weakness is the limited benchmark comparisons. Experiments could have been strengthened by comparing to further benchmark methods. However, even without these the paper is already good, as the main finding is the convergence issue of the Euler scheme.", "rebuttal": "We thank the reviewer for their thoughtful and constructive review. We greatly appreciate the recognition of the importance of comparing PINNs and BSDE methods, and the novelty of our approach. \n\n**Response to W1/Q1/Q2/Q3 (mathematical statements and assumptions):**\n\nThanks for bringing this up. We would like to clarify a few points. First, our main theorems (Thm. 4.1 and 4.2) are fully rigorous and include all the necessary assumptions in the theorem statement. On the other hand, the discussions that follow the formal theorem statements (including Eq. (4.3) and Line 903, 905) do omit exact technical assumptions, as we opted for clarity/brevity of exposition—we wanted to convey the message of *how* our main theorems could be interpreted without getting bogged down into too many technical details which we felt distracted from the main message. Nevertheless, given the feedback, we will add these details back into the draft.\n\nTo answer the reviewer’s questions, and to give a sense of how we plan to edit the paper (recall we are not allowed to upload a new draft), we’ll derive Eq. (4.3) more carefully below. This will not only illustrate the kinds of assumptions we make, but also fill in some of the analysis details (e.g., how we go from Riemann sum to integral) that the reviewer is asking about. The details presented here will be added to our revision.\n\nThe starting point for Eq. (4.3) is from the proof of Theorem 4.1 (Appendix C.1), where we show that:\n$$\n\\left|\\ell_{\\mathrm{EM},\\tau}(\\theta, x, t) - \\tau^2\\left(  (R[u_\\theta](x, t))^2 + \\frac{1}{2} \\mathrm{tr}( (H(x, t) \\cdot \\nabla^2 u_\\theta(x, t))^2  )  \\right)\\right| \\leq L_{d,u_\\theta} \\tau^{5/2},\n$$\nwhere $L_{d,u_\\theta}$ is a finite constant that depends on the dimension $d$, in addition to the Lipschitz constant of $D^2 u_\\theta$. Note that the latter is bounded, since we assume that $u_\\theta \\in C^{2,1}$ ($C^{2,1}$ is defined in Line 136).\n\nHence at this point, we have that:\n\n$$\\left|L_{\\mathrm{EM},\\tau}(\\theta) - \\frac{1}{N\\tau^2} \\sum_{n=0}^{N-1} \\mathbb{E}_{\\hat{X}_n}\\left[   \\tau^2\\left( (R[u_\\theta]( \\hat{X}_n, t_n))^2 + \\frac{1}{2}   \\mathrm{tr}(  (H(\\hat{X}_n, t_n) \\nabla^2 u(\\hat{X}_n, t_n))^2  )  \\right) \\right]\\right| \\leq  L_1\\tau^{1/2},$$\n\nwhere $L_1 = L_{d,u_\\theta}$ (this was necessary to get openreview to render the previous equation). \n\nOur next step is to argue that the map $(x, t) \\mapsto F_\\theta(x, t) := (R[u_\\theta](x, t))^2 + \\frac{1}{2} \\mathrm{tr}( (H(x, t) \\cdot \\nabla^2 u_\\theta(x, t))^2  )$ is Lipschitz for some constant $L_{F_\\theta}$. This is what the role of “basic regularity assumptions on $R[u_\\theta]$, $H$, and $u_{\\theta}$” (Line 150) plays; there is not a unique set of assumptions that render the map $F_\\theta(x, t)$ Lipschitz. Nevertheless, we specify one sufficient condition now: (a) $H, f, h[u_\\theta]$ are all bounded, and (b) $u_\\theta$ and all its derivatives up to 2nd order are Lipschitz and bounded. This is certainly not the weakest assumption; one could instead assume growth conditions of all the relevant functions and their derivatives as $(x, t)$ to go infinity, and then combine this with a high probability bound on the trajectories $X_t$ on $t \\in [0, T]$. In the interest of simplicity for this rebuttal, we opt for our simpler stated assumption.\n\nWith the Lipschitz condition on $F_\\theta$ in place, we can now apply order $1/2$ strong convergence of EM. Note that reference [34] (cited on Line 147)---specifically Theorem 10.2.2 of [34]---gives conditions on the drift $f$ and diffusion $g$ of the forward SDE for strong convergence to occur; these conditions are roughly Lipschitz $f, g$ and some growth conditions of $f, g$ as a function of $x$. We then combine strong convergence with Lipschitz $F_\\theta$ to argue: \n$$|\\mathbb{E}[ F(\\hat{X}_n, t_n) - F(X(t_n), t_n) ]| \\leq L_F \\mathbb{E}\\| \\hat{X}_n - X(t_n) \\| \\leq L_F C \\tau^{1/2}.$$\n\n(we had to drop the subscripts on $\\theta$ and write $X_{t_n}$ (from the paper) as $X(t_n)$ to get the equation to render). Note that if in the previous step we derive a Lipschitz constant $L_F(x)$ which depends on $x$, the definition of strong convergence involving second moments becomes important as it allows us to decouple the Lipschitz constant from the bound on $\\| \\hat{X}_n - X(t_n) \\|$ via Cauchy-Schwarz. Combining this together with the previous bound, we have\n\n$$\\left| L_{\\mathrm{EM},\\tau}(\\theta) - \\frac{1}{N}\\sum_{n=0}^{N-1} \\mathbb{E}[ F(X(t_n), t_n))  ]\\right| \\leq L_1 \\tau^{1/2} + L_F C \\tau^{1/2}.$$ \n\nOur final step is to treat the second term on the LHS as a Riemann sum approximation to a Riemann integral. From standard approximation properties of the left Riemann sum, we have\n$$\\left| \\frac{1}{T} \\int_0^T \\mathbb{E}[F(X(t), t)] dt - \\frac{1}{N}\\sum_{n=0}^{N-1} \\mathbb{E}[ F(X(t_n), t_n))] \\right| \\leq \\frac{L_F T}{2N} = \\frac{L_F \\tau}{2}.$$\n\nThis verifies the claimed Eq. (4.3):\n$$\\left| L_{\\mathrm{EM},\\tau}(\\theta) - \\frac{1}{T}\\int_0^T \\mathbb{E}[ F(X(t), t))  ] dt \\right| \\leq L_1 \\tau^{1/2} + L_F C \\tau^{1/2} + L_F \\tau/2 = O(\\tau^{1/2}).$$\n\nWe hope this clarifies the conditions under which Eq. (4.3) hold, and makes precise what is assumed in Lines 903, 905, 909. We note that Eq. (4.8) is also derived using similar arguments; we will also add its derivation to the paper.\n\n**Response to W2/Q4 (PDE background):** \n\n*(Note that this question was also brought up by Reviewer sGAy, so we have duplicated our response in both rebuttals to make the rebuttals self-contained).*\n\nThe non-linear PDEs we consider (cf. Eq. (3.1)) are the standard kind of PDEs considered for BSDE methods (see e.g., references [4-6]). There are many types of applications, including in quantitative finance (e.g., multi-asset options pricing via Black-Scholes, local/stochastic volatility extensions from the work of Duprie and Heston), optimal control (e.g., Hamilton-Jacobi-Bellman equations for finding value functions, Hamilton-Jacobi-Issacs for reach-avoid problems), reaction-diffusion in chemistry (e.g., Allen-Cahn and Fisher-KPP equations), committer/exit-time probabilities for transition path analysis (e.g. [R1]), and so on. We will include a paragraph about these applications to our revised draft.\n\nWe also note that another reason we consider the specific elliptic/parabolic form (3.1) given in the paper (i.e., $h[u]$ does not depend on $\\nabla^2 u$) is to favor the standard EM-BSDE formulation, since the form (3.1) means that the EM-BSDE approach does not need to compute Hessians of $u$ in the backwards SDE, whereas the Heun-BSDE does. In general, $h[u]$ can also depend on $\\nabla^2 u$ as noted in e.g., [Pham et al., 19] referenced by the reviewer, allowing for fully nonlinear second order PDEs.\n\n[R1]: E, Vanden‑Eijnden (2006), Towards a Theory of Transition Paths. J. Stat. Phys.\n\n**Response to W3/Q5 (other benchmarks/references):**\n\nThanks for pointing us to these references, and we will further incorporate them into our related work. For ease of reference, we will refer to these papers as [C1] to [C5]. In regards to direct comparison, our work focuses on the impact that the most widely used integration scheme in practice (Euler-Maruyama) has on the performance of BSDEs. Almost all the references above are focused on enhancements to the BSDE loss, and are also compatible with our proposed Heun-BSDE scheme. For example, we could modify our Heun-BSDE method to also utilize a control variate [C3] (i.e., split $u = u_1 + u_2$ each with their own properties), or to use Heun-BSDE in operator splitting settings [C1], or to apply to fully non-linear PDEs [C5], comparing these methods to their EM-BSDE + technique baseline. While this is an interesting direction which we will certainly incorporate in our revision, we believe the merits of our work can be judged on the basis of our current experiments which are focused on a direct comparison of EM vs. Heun, without extra confounding factors.\n\nWe next specifically discuss higher order Itô integration schemes and the paper [C4]. The next step up in complexity from a standard Euler scheme is a trapezoidal integration scheme; Heun integration is the explicit form of this. Direct application of Heun to integrate an SDE converges to the Stratonovich solution, which forms the basis of the study in our work. However, this is not the only approach. There are e.g., implicit Crank-Nicholson (CN) schemes and also higher order explicit Runge-Kutta (RK) approaches, as discussed in [C4]. We will certainly consider an empirical comparison in our future draft. However, we note several points: \n\n- First, the Heun-BSDE + Stratonovich formulation presented has several algorithmic advantages compared to other approaches: (a) it is an *explicit* method, and hence does not require solving implicit equations at every step, (b) it is very simple to implement, unlike e.g., higher order RK methods for Itô which require non-trivial implementations for handling diffusion terms that are not commutative---which is essentially always the case for backwards SDEs where the diffusion term depends on $\\nabla u$. Thus, part of our contribution is also to advocate for the advantages of reinterpreting the SDEs as Stratonovich, in addition to using more sophisticated integrators. \n\n- Second, the schemes presented in [C4] are not directly comparable as described, since they are designed for the recursive BSDE loss formulations, where a separate network is trained at every discretization point. On the other hand, in this paper we focus on self-consistency losses (see the discussion starting on Line 116 for the difference), which are more practical as only one network is trained. The integration schemes can be applied in the self-consistency setting as well (which we will investigate in the revision), but this is not what is specifically considered in [C4]."}]

Author rebuttal:

Meta-review: This paper investigates a critical issue in BSDE-based deep learning methods for solving PDEs: a performance gap relative to PINNs caused by a discretization bias from the standard Euler-Maruyama (EM) integrator. Through both compelling theoretical analysis and empirical evidence, the authors demonstrate that this bias is a major drawback. To resolve it, they introduce a novel approach that uses the Stratonovich SDE interpretation and the Heun integration scheme, which successfully eliminates the bias. While initial concerns about high computational overhead were noted, the authors' rebuttal effectively addresses this weakness, showing substantial efficiency improvements that make their fixed algorithm highly competitive. Overall, the work is a significant contribution, backed by rigorous theory and solid experiments, and presents a promising, viable alternative to PINNs.
Decision: Accept (spotlight)