Physics-Informed Inference Time Scaling for Solving High-Dimensional Partial Differential Equations
Authors: Zexi Fan
Yan Sun
Shihao Yang
Yiping Lu
Venue: ICLR 2026 Poster (2026)
Keywords: AI for Science
Inference-time Scaling
Deep learning
Curse of dimensionality
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)
TL;DR: We introduce SCaSML, a framework that improves pre-trained PDE solvers at inference time without retraining by deriving and efficiently solving a new PDE that governs the model's error, provably accelerating convergence.

Abstract: Solving high-dimensional partial differential equations (PDEs) is a critical challenge where modern data-driven solvers often lack reliability and rigorous error guarantees. We introduce Simulation-Calibrated Scientific Machine Learning (SCaSML), a framework that systematically improves pre-trained PDE solvers at inference time without any retraining. Our core idea is to derive a new PDE, which we term the Law of Defect, that precisely governs the error of a given surrogate model. Because this defect PDE retains the structure of the original problem, we can solve it efficiently with traditional stochastic simulators, yielding a targeted correction to the initial machine-learned solution. We prove that SCaSML achieves a faster convergence rate, with a final error bounded by the product of the surrogate and simulation errors. On challenging PDEs up to 160 dimensions, SCaSML reduces the error of various surrogate models, including PINNs and Gaussian Processes, by 20-80%. SCaSML provides a principled method to fuse the speed of machine learning with the rigor of numerical simulation, enhancing the trustworthiness of Al for scientific discovery.

Published as a conference paper at ICLR 2026
PHYSICS-INFORMED INFERENCE TIME SCALING FOR PARTIAL DIFFERENTIAL
SOLVING EQUATIONS VIA DEFECT CORRECTION
Zexi Fan
Yan Sun & Shihao Yang
School of Metathetical Science
H. Milton Stewart School of Industrial and Systems Engineering
Peking University
Georgia Institute of Technology Atlanta, GA, 30332 shihao, yangeiingeeeee www.ch.teeheee.... edu
Beijing, China
fanzexi francisestunuuu.m pku. edu, on
Yiping Lu
Department of Industrial Engineering and Management Sciences
Northwestern University
Evanston, IL. 60208 yiping Lugnortsenntttt edu
ABSTRACT
Solving high- -ish-niimnsional partial differential equations (PDEs) is a critical chal- lenge where modern data-driven solvers often lack reliability and rigorous error guarantees. We introduce Simulation-Calibratee Calibrated Scientific Machine Learning (SCaSML), framework that systematically improves pre-trained PDE solvers inference time without any retraining. Our core idea is to use defect correction method that derive a new PDE, which we term the Structural-preettee Law of Defect, that precisely describes the eтror of given surrogate model. Because this defect PDE retains the structure of the original problem, we can solve it efficiently with traditional stochastic simulators, yielding targeted correction to the initial machine-learnel solution. We prove that SCaSML. achieves a faster convergence rate, with final error bounded by the product of the surrogate and simulation errors. On challenging PDEs up to 160 dimensions, SCaSML reduces the eгrоr of various surrogate models, including PINNs and Gaussian Processes, by 20-80%. SCaSML provides principled method to fuse the speed of machine learning with the rigor of numerical simulation, enhancing the trustworthiness of AI for scientific discovery.
INTRODUCTION
Solving high-dimensional partial differential equations (PDEs) is a fundamental challenge across science and engineering. Many critical phenomena are modeled by semi-linear parabolic PDEs whose dimensionality scales with the number of underlying components, challenge often termed the curse of dimensionality. Key examples include the imaginary-time Schrödinger equation in quantum many- body systems, nonlinear Black-Schole equations in finance, and the Hamilton-Jacobs Jacobi-BeBlllaan equation optimal control (Bellman 1954). Traditional numerical methods, such finite element and finite difference schemes, become computationally intractable in high dimensions (Larsson & Thomée, 2003). While stochastic simulation methods can be effective, they often suffer from high variance (Briand & Labart 2014). In response, Scientific Machine Learning (SciML) has emerged as powerful alternative, using neural networks and other data-driven models to approximate PDE solutions (Karniadakis al. 2021] Han al., 2018a| Raissi al. 2017). However, the "black-box" nature of these models can introduce subtle biases, and they often lack the rigorous error guarantees of their traditional counterparts, raising concerns about their reliability for safety-critice applications.
Recent breakthroughs in large language models (LLMs) have shown that allocating additional computational resources inference time dramatically improve output quality, a phenomenon

---

Published as a conference paper at ICLR 2026
known inference-tinee scaling (Snell al. 2024 Wei et al. 2022) This success inspires central research question:
Can we leverage additional computatiio at inference time systematically refine and provably improve a pre-trained surrogate model-allocating more compute to harder PDE states just as LLMs spend more search or planning on harder -ithhut any retraining or fine-tuning?
In this work, we provide an affirmative answer by introducing Simulation Caalbbrrte Scientific Machine Learning (SCaSML.,, novel physies-informee framework for improving SciML solvers at inference time. We focus on broad class of semi-linear parabolic PDEs of the form:
du + Lu + F(u, a Vu) 0, on [0, T) Rd Br
(I)
(u(T, g(y),
on Rd
where Lu:= (1u, Vu) Tr(o (e))(u)o) is second-order linear differential operator. SCaSML operates in two stages. First, a standard SciML solver à (c.g.. PINN (Raissi et al.) 2017)) Gaus- sian Process (Chen et al.) (021). or Tensor Network (Richter al. 2021 is trained find an approximate solution. At inference time, rather than directly accepting ú, we invoke defect- correction method (Bank Weiser 1985| Stetter 1978; Böhmer et al., 1984) derive a gov- erning equation for the approximation error-ils defect as Law of Defect (Figure Crucially, unlike classical grid- based defect correction which is intractable in high dimensions, we show that this new PDE describing the exact defect inherits the semi-lineal structure of the original problem. This structural preservation allows to solve it efficiently using well-established stochastic simulation algorithms based on the Feynman-Kae formula. This simulation step acts as a targeted correction, leveraging additional compute refine the initial surrogate prediction.
û. We term this
the wrucurral-prreee
Our main contributions are summarized as follows:
We propose SCaSML, the first physics-informet inference-time scaling framework that improves pre-trained Scisurrogate model at inference-time, without any retraining or fine-tuning. SCaSML uses defect correction method that corrects a pre- pre-trained surrogate SCiML model by deriving and solving new PDE a branching Monte Carlo Simulation that approximates its error, which we call the Structural weweccurallgreeeeeeeee -preserving Law of De 7 Notably, this characterization of the defect is, to our knowledge, the first derivation that preserves the semi-linear structure essential for high-dimensional Monte Carlo solvers.
We theoretically prove that the final error of SCaSML is bounded by the product of the surro- gate model's is error and the simulation error. Analogous the classical defect-correction litera- ture where each correction step systematically improves the convergence rate- we establish analogous result for our Monte-Carlo defect-correction procedure at the first time. The improved convergence rate is corroborated empirically in Section 3 with further comprehensive findings presented in Appendix G.3
We conduct extensive numerical experiments on challenging nigh-dimensional PDEs (up to 160 dimensions). Our results show that SCaSML significantly reduces approximation errors by 20-80% across various surrogate models with high statistical significance (p < 0.001), demonstrating its flexibility, practical efficacy and potential to mitigate the curse of dimensionality. We also demonstrate that, with inference-time scaling, smaller base PINN can outperform a larger PINN under the same inference-time compute budget by spending its additional computation on targeted refinement rather than parameter count. This enables elastic compute: users trade inference time for accuracy on demand.
METHODOLOGY
The core of our SCaSML framework is the derivation of new PDE that describes the error of pre- trained surrogate model. We term this the St st.ructuaal-presevv.. solving this auxiliary PDE inference time, prediction. To build intuition, we first introduce this concept in the context of linear parabolic equations before extending it to the general semi-linear case.
Law of Defect. By compute a precise correction the surrogate's 's
2

---

Published as a conference paper at ICLR 2026
a) SCaSML Framework Pipeline
Different Themes for solving Dimensional PDts
Physics- Informed Inference Times Scaling via Simulation Calibritaed Scientífic Machine learning
Theme Train Surrogate Model
Step 1. Train Surrogate Model
il boo
Artecc
1/10-30
bitecce
Step 2. Model the Residual via PDE
We show thet the defect v can formulated as the solution
another serni inear Paraboli equation [Ext
Step Simulate Residual
bug
789
Cancharic
b) Derivation of the Structural-preseevie
Law of Defect
Semi-linear PDEs:
Step 1: Surrogate Model
Step 2: Law of Defect
Norlineer term
(1--21 (f
Au + f(u) = differential operator du
Lú + f(i) = = (1)
(2)
A + A(x () + ( ((0) =c
Small termis easy to simulitt
e.g.
surrogatee approximately solve eqqattion
semi-ioear equation that
dr
Figure 1: Overview of the SCaSMI. framework. (a) A pre-trained surrogate model ú provides an initial, approximate solution to the PDE. At inference time, SCaSML calculates the defect ä u à via stochastic simulation and adds it back to the surrogate prediction, yielding more accurate final solution u = â + u. (b) The St
Law of De fect is derived by subtracting the PDE approximately satisfied by the surrogate ú from the original PDE. This process yields a new semi-linear PDE that describes the defect ú, enabling its estimation through simulation.
2.1 WARM UP: THE STRUCTURAL LAW OF DEFECT FOR LINEAR PARABOLIC PDES
Let us begin with high-dimensional linear parabolic PDE, a simpler setting that clarifies our core idea:
8u Gr
+ (1, V Vy + Tr Hess, fr, y), on [0, I') x RRd on Rdd
(T 3) g
Suppose we have pre-trained surrogate model á that approximates the true solution u. This surrogate inevitably imperfect, producing a residual when plugged into the PDE. Our goal is to run defect- correction method inference time. The defect-correction method (Böhmer et al., 1984; Stetter 1978) is classical numerical strategy that improves approximate solution by formulating and solving an equation for its residual-induces error. The first step is to find new equation that describes the defect u(r; y) u(r, y) u(r, y). which represents the true, unknown error. To achieve this, we define this residual as:
c(r, v) = f(r,v) - ( au
+ Vo + TT(o Hess, ú ūo))
(3)
0r
By subtracting the equation for á from we arrive at the following governing law.
Definition 2.1 prrserina Law of Defect for Linear PDEs). The defect û == u â the solution to the linear parabolic PDE:
Oū + (1, Vyū) + Tr Hess, ūo) e(r, v(r,y). (0, Rd. Rdd.
dr (T. g(y) =TT.y).
(4)
3

---

Published as a conference paper at ICLR 2026
This structural-preeervii Law of Defect allows us solve for the error i directly. Since (4) is a linear PDE, its solution can be expressed probabilistically via the Feynman- -ac formula:
M(8, x) = To (o(()) a(T, x;") + S." ctt, Xi )
(5)
where (X," he(aTT] is the stochastic process associated with the PDE's linear operator. This representation allows us estimate the defect ä using Monte Carlo simulation.
Remark 2.2 (Regards Training and Inference Separation). Training corresponds solving the PDE globally on the entire domain, learning a map that approximates the solution everywhere. In contrast, inference-tinn -ieee correction solves the PDE only specific, user-specified state. This separation is natural and parallels standard practices in machine learning: a base model is trained once to answer all queries, while compulationally intensive refinement (like beam search or planning) is invoked at inference time only when high precision is required for a specific input. This separation enables "elastic compute," allowing users to trade inference time for accuracy on demand without incurring the massive fixed cost of retraining the global model.
Intuition for Faster Convergence. Why does this two-step process converge faster? The variance of the Monte Carlo estimator for u in (5) depends on the magnitude of the integrand, which is primarily the surrogate's residual e. A more accurate surrogate (i.e., smaller e) leads lower simulation variance. If the surrogate achieves an error of e(û) ~ m "" using m training points, the residual will be of similar order, and the variance of our Monte Carlo estimator will be of order vn m=-27 By averaging over new Monte Carlo paths at inference time, the final statistical error becomes /m-22 27 pmm22/mm--mm m 1-2 Thus, for a total budget of 2m function evaluations, SCaSML achieves a faster convergence rate than both the surrogate (m-7) and a naive Monte Carlo solver (m-11/2)
Why Use Monte Carlo for Correction? Neural networks and other common surrogates exhibit a spectral bias, preferentially learning low-frequency (smooth) components of the solution first (Rahaman et al. 2019). Consequently, the residual error e is often a While challenging for many function approximators, Monte Carlo methods are perfectly suited for this scenario, as their convergence rate is independent of the integrand's smoothness. This makes Monte Carlo an ideal choice for the correction step, as can efficiently average out the complex error signal left behind by the surrogate model.
irregular function.
2.2 EXTENSION TO SEMI-LINEAR PARABOLIC PDES
We now extend the Monte-Carlo based inference time defect-correction procedure to the general semi-linear PDE in [T. Let û be a surrogate solution. We define its residual with respect to the PDE dynamics and the terminal condition as:
c(r, y) 8r ( = g(v) (T,
aa + Lu F ToTTVVV Vvú),
(6)
By subtracting 6 from the original PDE in [T.. we obtain the governing law for the defect ä = u û. Fact 2.3 (Structural-presevving Law of Defect for Semi-linear PDEs). The defect usr, y) = u(r, y)
a(r, y) is the solution the following semi- lini-r parabolic equation:
+ Lü + F(ü, a T TV Vvū)= 0, on [0. T) x Rdd on Rd.
(7)
( = g(y).
where the modified nonlinear term 1' is given by F(ū VV,) = F(u+ ü TV F(ù, VVū +
Law of Defect 7 retains semi-linear struc- ture. This is the key property that allows us to apply powerful stochastic solvers, such as the Multilevel Picard (MLP) iteration (Hutzenthaler al. 2019) estimate the defect ú and correct the initial surrogate u.
Notably, the Structural-preservee
How does the structural-preservvis correction methods? Classical finite element methods admit a weal-chharcterrzedd asymptotic error expansion Strang et al. (1973), which enables defecctccorection schemes systematically remove the
Law of Defect differ from classical defect-
4

---

Published as a conference paper at ICLR 2026
leading enror term and improve convergence rates (Zienkiewicz & Zhu, 1992a b Bank & Weiser 1985). In contrast, no such asymptotic structure is available for neural networks: NN approximations lack any mesh-refinement hierarchy, their enrors do not exhibit a polynomial expansion with respect to a single resolution parameter, and the optinization-inducee perturbative decomposition. A different family of debiasing techniques in numerical PDEs relies on iterative solvers such as Newton methods (Stetter 1978, Dutt et al. 2000 Xu 1994; Böhmer 1981) and quasi-Newton methods (Jameson al. 1974] Heinrichs] 1996). However, these methods present two fundamental limitations in our setting. First, iterative updates produce only approximate corrections, whereas our law of defect is an exact analytical identity that delivers a closed-form unbiased correction in a single step. Second, embedding iterative methods into a Monte-Carlo or Feynman-Kac framework is highly inefficient: each iteration requires recomputing residuals and Jacobian actions through additional Monte-Carlo estimators, producing nested simulation hierarchy whose convergence rate rapidly deteriorates- from the standard (N-1/2 rate for single Carlo level, to O(N-1/4) for second iteration, ( 1/) for a third, and so on as more levels are introduced. Practitioners are therefore forced balance carly-termination errors against the rapidly declining statistical efficiency of nested Monte-Carle estimates, making these approaches both computationally expensive and numerically unstable.
approximation error provides no
Practical Scenarios In many applications, the quantity of interest is required only at a single state rather than across the full domain. For example, in optimal control and financial pricing (e.g.. nonlinear Black-Scholes (Eskiizmirliler et al., [222]] Santos & Ferreira| (024)), practitioners need the value function and its gradient only at the current state to determine the next action or hedge;
forward simulations can then be used to compute the Bellman eтror and correct the current decision. In rare-event analysis and committor problems in molecular dynamics. neural committor estimators (Khoo et al. 2019 Li et al. 2019 Hua et al. 2024 Lucente et al. 2019) can be refined using a small number targeted simulations initiated from designated configurations. In goal-oriented estimation (Becker & Rannacher 1996; 2001), the objective is often specific functional of the solution rather than the full field. In all such settings, training a surrogate to high global accuracy is computationally wasteful. Our method uses the surrogate for fast initial approximation and then applies targeted Monte Carlo refinement at inference time, allocating computational effort precisely where accuracy is needed.
2.3
SIMULATING STRUCTURAL
LAW OF DEFECT USING MULTILEVEL
PICARD ITERATION
The defect PDE 7 is a semi-linear parabolic equation of the same structural form as 24 with different closed nonlinear term. Under the standard regularity assumptions, the pair a
V)ū") admits the Feynman-Kaan and Elwoutt- representations presented in Hence u is the fixed point of the expectation operator & on Lip((0, T] x Rd gl++d) )):
27
28
{)))),
with & given exactly as in the appendix at (28). Intuitively. the operator + is a Feynman type backward propagator. Given an approximation of the solution at a future time, & maps it back to the present by running a stochastic simulation forward in time and averaging over all resulting trajectories. Concretely, for each initial state x, it computes the expected terminal payoff together with the accumulated contribution of the nonlinearity F along the simulated path. The exact solution u" is therefore characterized the fixed point of this propagation: inserting u² into the simulation leaves it unchanged, i.e., tu* = u*. Standard Picard iteration ü = ((ūx) converges to u under standard regularity assumptions (Yong & Zhou 1999. Theorem 3.4).
Multilevel Picard (MLP) method (E et al.
2021 Hutzenthaler et al. 2020a), uses Multilevel Monte Carlo (MLMC) Giles 2008 2015 to simulate the telescoping formulation E(ū,)] =
EB [ + Σm E ( & (ūū )) The MLMC method exploits a hierarchy of approxi- mations T (ūn), (ū,),, ((.,). ranging from the coarsest to the finest resolution. Crucially, consecutive approximations (T(ū)) and ())ū-1))(() are generated using the same underlying sample path i, which induces a strong positive correlation between them. As a result, the variance of their difference is significantly reduced. Moreover, as the level / increases, the iterates converge linearly ú, ú + O, and the variance of the difference decreases linearly toward zero. As a consequence, the required number of samples - can decrease as increases, meaning very few expensive samples are needed at the finest levels. The majority of the computational cost is thereby shifted to the coarser levels, significantly reducing the overall complexity of the estimation.
5

---

Published as a conference paper at ICLR 2026
Build surrogate solution
Formulate St ructural preeerviin Law of Defect describing
7
Use MLP simulate 7 and at inference correct ti
Figure Flow diagram of SCaSML. We formulate the error u - â of surrogate solution the solution Structural -preserring Law of Defect 77. new semi-linear PDE. At inference time, we approximate u - â via solving Structura penectural--pae -prellrreengnng Law of Defect using Multilevel Picard (MLP) iteration. The generated estimation of helps to calibrate the surrogate solution ü.
Another factor affecting the variance is how the time integral is computed; used two MLP variants to simulate structural-preservv
Law of fecect
Quadrature MLP: (E al. 2021) Simulate the time integrals by the Gauss Legendre quadrature.
Full-history MLP:
Hutzenthaler et al. 2021) Simulate the time integrals by Monte Carlo.
We leave the preliminaries and implementation details of the MLP methods Appendix B.2.1 The overall SCaSML procedure, which involves first training surrogate model and then solving the Strucurra -er-epeeeeeeggggee Law of Defect with MLP methods to correct it, is summarized in Algorithm C
2.4
PROVABLY ACCELERATED CONVERGENCE
We now provide theoretical guarantees for SCaSML, showing that it achieves provably faster convergence rate. For simplicity, we present results for the case μ = 0 and o sIg. Our analysis relies on the assumption that the pre- pre-trained surrogate is reasonably accurate.
Why SCaSML Enjoys Provable Faster Convergence Rate. The Monte Carlo error in MLP methods depends on the scale of the terminal defect ğ and the modified nonlinearity F, which depends on the error of the surrogate model. A more accurate surrogate model yields smaller ğ and F. resulting in reduced variance during inference. If the surrogate achieves an eтrоr of e(ū) m (Assumption 2.4 from m training points, then the variance is O(m -2)). During inference, we
m m = m (Blanchet al. 2023). With total computation cost of 2m function evaluations, SCaSML therefore attains a convergence rate that surpasses both the surrogate method m -Y m=> and the standard MLP / Monte-Carlo estimator m -1/2, Full constant-tracking and rigorous proofs are Appendices and E Assumption 2.4 (Surrogate Model Accuracy). Let the true defect be well- bell-behaved such that suPte(o.r] [| lww.. < 0C. We assume the surrogate error is bounded by a measure e(ū), such that for constants CF,1, CF2 > 0:
average over m additional Monte Carlo paths, which reduces the statistical error
1. L Residual: supr.y le(r, b)1 S CF.1e(0).
2. W1.00 Error: sup, [ 1lw I CF.2e(()).
Proof Sketch. Our main theoretical result stems from the observation that the computational complexity of the MLP solver depends on the Lipschitz constant of the nonlinearity F and the magnitude of the "source terms". They appear multiplicatively because nonlinearities- their Lipschitz hounds-propagatt and magnify variance at every Picard iteration. The "source term" driving the Multilevel Picard simulation for the defect is the residual already reduced by the surrogate. At the same time, We show that the regularity in the law of defect no worse than that of the original PDE, ensuring that the refinement introduces additional smoothness requirements. Combining the previous fact, a more accurate surrogate makes the defect PDE "easier" to solve. This leads to main error bound.
Theorem 2.5 (Global L Error Bound). Under standard regularity assumptions on the PDE coeffi- cients (Assumptions E.2 D.7), the global L² error of the SCaSML estimator using a full-history MLP approximation Ü N.M with N levels and use M Monte Carlo samples 1-th level is bounded by:
sup (f
Ü ÜN,Nr(1, ū((,,z)
E(M,N) (Cr e(ū)),

---

Published as a conference paper at ICLR 2026
Table 1: Comparative performance full-history SCaSML against the surrogate model (SR: PINN or GP) and naive MLP solver. We report norms. Bold values indicate the best performance in each category. SCaSML consistently achieves the lowest error across nearly all settings.
total runtime (s) and relative errors in L L and L
Problem
Time (s)
Relative L Error
1 Error
L' Error
SCaSML
MLP
SCaSML
MLP
SCaSML
MLP
SCaSML
10d
0.45
6.77
13.31
5.20E-02
2.275-01
2.74E-02
2.505-01
9.06F-01
L65K-01
5.395-02
1667-01
1.78E-02
204
0.54
6.75
17.11
9.008-02
2.558-01
4.728-02
4.728-01
1.35E400
3.308-01
9.378-02
2.37E-01
4,528-02
22.44
1.45-011
2.38B-01
9.72E-42
2,04E+00
1.392400
7.99-00
1615-01
2.848:01
1.4E-00
60d
0.26
6.94
37.59
3.136-01
2
1.32--01
3,24E-00
2,05E+00
1,578+00
5355-01
4.07E-01
2.063-01
VR-PINN
203
0.54
6.80
10,59
1,178-02
8.30E-02
4,038-83
3.208-02
2.968:01
2.368-00
5.308-03
3. 3.39E-02
1.29K-03
404
0.29
8.11
14.09
4.065-02
1.048-01
2.928-42
8.431B 02
3.571: 01
7.438-02
2.008
4.36E 02
1.24E-02
38.
3.956-02
1.176-01
2.888-42
3.930-01
7.20E-02
1.945-02
4.82E-02
1.223-02
X0M
3.65
11,78
42,50
6,74R-02
1.19F-01
5.648-02
1.90F-001
3.35F-01
1.808-01
5.218402
4 4.735-02
2.468-02
203
1.74
10.55
61.82
1.472-01
1.90-01
6.23E 02
3.541 01
5.72E:01
2.54E
7.018:42
8.00E
2.48E-02
VT-GP
408
1.78
12.28
61.28
1.81E-01
2.20E-01
8.55E 42
4.00E-01
8.71E-01
3.00E-01
9.19E-02
9.06E-02
3.821-02
57.79
2,40F-01
2
1.28F-01
3.84FF-01
9.50F-01
2-84K-01
1.275-01
9.99F-02
6.11K-02
804
1.69
10.12
60,69
2.66F-01
3,02F-01
1.52F-01
3,611-001
1.914400
2,618-011
1.45--01
1.09-001
7.59R-022
1000
0.42
8.27
21.33
7.97E-02
5.63E-00
5.53E
7.82E--0
1.26E+01
6.828-01
1.406-011
1
8.728-02
1200
0.32
8.52
23.98
9.40F-02
5,50F-00
6.66F-02
9.06F-01
1.277401
8.063-01
1.749-01
1:224401
1.068.01
140d
0.44
27.31
9.876-02
5,578:000
6.84E-02
9:908--01
1.278401
8.96--01
1936-01
1.238401
1.12E:01
1600
03.34
8.09
29.95
1,126-01
5,27E:000
9.94F-02
1,404-00
1,28E401
1.30E.00
2
1 1.236401
1.79K-01
1000
0.32
7.59
58.51
1.41F-02
8.90F-02
1.11F-02
9.58F-00
6.37F-01
8.588-402
1.879-02
9.74F-02
1.388-02
DR
120d
0.33
7.16
68.28
1.118-02
9.138402
1.03E-02
7.508-002
5.74E-01
6.50E-02
1.298-022
9:976:02
1.29K.02
1400
0.42
7.75
79.99
3.226-02
8:97E-02
3.00E:42
1.82E-01
8.56E 01
1.72E-01
4,00
9.77E 02
3.75E-02
160d
0.37
7.22
86.77
3.45K-02
9,008-02
3.228:42
2,088-01 01
8.02E:01
1.982-01
4.3089002
9.75E:02
4.002:02
where u ú= (ū, vV is the true defect and its gradient, and E represents the error term of the underlying MLP solver, which depends M and N but is independent of the surrogate.
Theorem 2.5 shows that the final error is the product of the MLP simulation error and the surrogate model error. This synergistic relationship implies that the computational cost to reach a global L error of æ is reduced from -(++)) for a naive MLP solver to (d (2++) e(ū)2++) for SCaSML (see Corollary F.9 in Appendix). This means the cost of our correction step decreases as the quality of the initial surrogate improves. This directly leads to an improved scaling law.
Corollary 2.6 (Improved Scaling Law). Under the assumptions of Theorem 2.5 suppose the surrogate model's ss error scales as e(û) = O(m-)) with m training points. By allocating an additional m samples for the inference-tine time simulation, the total error of the SCaSML procedure improves from O(m77) to O(m (
NUMERICAL RESULTS
We now empirically validate the SCaSML framework across a suite of challenging high- nigh-dimmnsional PDEs. In each experiment, first train baseline surrogate model û (either Physics Informed Neural Network or Gaussian Process) to obtain an approximate solution. Then, at inference time, apply full-history Multilevel Picard (MLP) solver the St sucturr purreeea--preeeeeee --reese rving Law of fect (Fact 2.3 to compute correction term u. The final SCaSML solution is the sum USCASML = à + u.
The primary goal of these experiments is to demonstrate the value added by the correction step. Thus, our key comparison is between the baseline surrogate model (SR) and the final corrected solver (SCaSML). We also include the naive MLP solver for reference, to show that the hybrid approach succeeds where pure simulation often fails. Our implementation leverages JAX (Bradbury al.. 2018) and DeepXDE (Lu et al., 2021 for efficient, parallelized computation.
As shown in Figure 3.. SCaSMI. consistently tightens the error distribution compared to the base surrogate. Refer to Appendix G.6) for detailed pointwise error maps. Figure 33 demonstrates SCaSML's effectivv inference-time scaling: as more computational resources (i.e., Monte Carlo samples) are allocated, the accuracy of the solution progressively improves. A comprehensive comparison of eтror metrics and timings is provided in Table [] and the empirical validation of our theoretical scaling law is shown in Figure 4 More experiments, including statistical significance tests (p & 0.001, Appendix G.4) and fixed-budget efficiency comparisons (Appendix G.7). are shown in the Appendix G
7

---

Published as a conference paper at ICLR 2026
3.1 LINEAR CONVECTION CIVEECTION-DIFUSSIO EQUATION
Problem Formulation. We investigate linear convection-diffusii equation given by ar u(r, y) V. y)) + Auu(r, y) = O, (r: y) E (O,T) x R with the terminal condition u(T, y) 9 +T, E Rd. This PDE admits the explicit solution u(r; y) > 2i=12 34 +
Experimental Setup. The problem is solved over the hypercube [0, 0.5] x [0. 10.5. for dimensions {10, 20, 30, 60) with Dirichlet boundary conditions enforced by the PINN loss. We deploy Physics-Imi Inysies-Informed Neural Network (PINN) with 5 hidden layers, 50 neurons cach, and tanh activation function. Training uses the Adam optimizer (learning rate 7 x 10 = 0.9, 32 = 0.99) for 10t iterations. At each iteration, the network is trained on 2.5 x 10 interior and 102 boundary collocation points. For the inference step, we use 2-level simulation with M = 10 as the basis of Monte Carlo samples at each level in Multilevel Picard iteration for the tabulated results and M E {10.... 16] for the scaling study. A clipping threshold of 0.5(d + 1) is applied to the solution and gradients for both the naive MLP and SCaSML.
Results. As reported in Table 1 (LCD), SCaSML achieves a reduction in the relative L error from 20% to 56.9% compared to the baseline PINN surrogate. Moreover, SCaSML exhibits robust inference scaling, with performance improving as more inference time is allocated (see Figure 10
3.2 VISCOUS BURGERS EQUATION
Problem Formulation. Next, we consider a viscous Burgers equation from (Hutennhalerr et al. %
2019), a standard benchmark for nonlinear PDEs: du dr
+
<- (
1, vuu
+ 2 Ayu
oou Σí (00 V u = O, with terminal condition (T, y) 1Txpp(T+> 25) exp(r+ Σí. +inm mc)
exp v.)
The exact solution
is u(r 3) 1+oxp(r+ Ei 21)
Experimental Setup. We solve the PDE on [0, 0.5 x [-0.5,0.5] for dimensions d (20, 40, 60, 80) with op = 22. We test SCaSML. with two types of surrogates. The PINN was trained for 10a iterations using the Adam optimizer (learning rate x 10-1, 333 = 0.9, and /3 = 0.99), utilizing 2,500 interior, 100 boundary, and 160 terminal condition sample points. A Gaussian Process (GP) regression surrogate was trained over 20 iterations via Newton's 's method, using 1,000 interior and 200 boundary points. For the 2-1evel MLP and SCaSML solvers with the basis of Monte Carlo samples M = 10, we clipping thresholds of 1.0 and 0.01, respectively, to handle the nonlinearity.
Results. SCaSML demonstrates strong performance with both surrogate types. For the PINN surrogate (VB-PINN), it reduces the relative L error by 16.2% 66.1%. For the GP surrogate (VB-GP), the reduction is even more pronounced, ranging from 42.7% to 57.5% (Table []. This highlights SCaSML's versatility as plug-and and--lay corrector for different SciML models.
3.3 DIGN-DIMENSIONAL
EQUATION
Problem Formulation. To showcase SCaSML on problems central to control theory, we tackle a high-dimensional Hamilton-Jacoke Belllllllan (HJB) equation arising from linear-quadratii-Gaant Gaaussian
(LQG) control problem (Han et al. 2018b). The HJB equation is given by du
Or Ayu || V u/ = 0,
1+ 1+ (e.. -w++ " ++2.. where C₁ and 02 are
with terminal condition u(T, y) = log(
independent random draws from interval [0.5, 1.5). The reference solution is computed via u(r, y)
log eeeeeeeppuulT y /2Wr-,)) with sufficiently large sample sizes(e.g. 100d).
Experimental Setup. Following (Hu et al. 2024), we use complex, non-trivial terminal condition and evaluate the problem in very high dimensions, d E {100, 120, 140, 160). The PINN surrogate is trained for 2.5 x 10x iterations on the domain [0, 0.5] x [dd where [d is the unit ball in &d, with 100 interior and 1,000 boundary points per iteration. We use the Adam optimizer with a learning rate of 10 1, B1 = 0.9, and B3 = 0.99. For inference steps, we set total level n = 2 and the basis of Monte Carlo samples each level M 10, where is the total level and M is sample used level I for 0 E < n. To stabilize the simulation for this strongly nonlinear problem, we use a clipping threshold of 10 for the naive MLP and a much smaller threshold of 0.1 for SCaSML,, reflecting
8

---

Published as a conference paper at ICLR 2026
(a) Error Distribution of Test Points
134 10.
Loc Jic
(b) Inference-Time Scaling SCaSML with Increasing Computational Budget
improvmerttso
Improveert
Insrowement
Insrovenere
70.0
30
Improseentt
62.5
10"
10
10"
10"
10'
10'
Evaluation Numbers
N umbbrs
Evaluation Numbers
Evaluation Numbers
Linear Convection-Diffinn
Burgers
Linear Quadratic Goussian
Diffusion
(c) Overall Performance Comparison Across Metrics and Runtimes
Metheds Systems
Relative Error Computation Time
Comp compuration Time Cs. 10g saal)
31--1
24--3
25402
(allaivve (log
1449-10
44:00
34e-11
336-22
3/462
10-
15400
35e-11
16442
124.00
25e-11
25x02
2/400
33w+01
12x+01
VS
PWN
VE-CP
LOG
0A
12
646-21
8.1642
SR.
Time
HLP
Tima
SCISKL
Time
HedE:
Figure 3: Efficiency and performance of the SCaSML methodology. (a) Violin plots showing the distribution of pointwise errors. SCaSML consistently reduces the mean error and tightens the distribution compared the surrogate (SR) model. (b) Inference-time scaling. As the number of inference- eererennee-imeeeeee simulation samples increases, SCaSML's error steadily decreases, demonstrating effective use of additional compute. (c) Summary of performance. The left panel shows that SCaSML (blue stars) consistently achieves lower T.2 error than both the surrogate (SR) and naive MLP methods across all problems. The right panel (heatmap) confirms that SCaSMI. also dominates in J, T and T,! error metrics.
the smaller magnitude of the defect. To accelerate computations, we use Hutchinson's method to stochastically estimate the Laplacian and divergence terms, sampling d/4 dimensions at each step (Hutchinson] 1989: Girard 1989 Shi et al. 2025).
Results. In this challenging high- nigh-ddmensional setting (LQG), the naive MLP solver fails entirely, producing large errors. In contrast, SCaSML successfully refines the PINN solution, reducing the relative 2 error by 11.7% 30.8% and achieving the lowest error across all metrics (Table
3.4 DIFFUSION-EAATTIR EQUATION WITH AN OSCILLATING SOLUTION
Problem Formulation. Finally, we consider a diffusion-reaction system designed to have a highly
oscillatory solution (Gobet & Turkedjiev 2017 Han al. 2018b), making particularly difficult
for standard neural network surrogates Ör A + min 1, (u u² 2} O, where u" (r. v)
1.6 + sin (0.1 Ed_ W:) exp( 0.01d(r-1) the exact solution.
2
Experimental Setup. We solve the problem for dimensions d E {100, 120, 140, 160) the domain [0. 1] x Bd. The PINN surrogate is trained for 2.5 x 10 iterations with 1,000 interior and 1,000 boundary points, using the Adam optimizer with learning rate of 10-3 3,-3,0 ,1 0.9. and B2 = 0.99. For inference steps, we set total level n = 2 and the basis of Monte Carlo samples at each level M = 10. The MLP and SCaSML. solvers use clipping thresholds of 10 and 0.01, respectively.

---

Published as a conference paper at ICLR 2026
Due to the solution's oscillatory nature, we found that the Hutchinson estimator for the Laplacian introduced instability; therefore, we computed the full Laplacian in this experiment.
Results. Even though the PINN surrogate is already quite accurate for this problem, SCaSML, is still able to provide consistent refinment. As shown in Table (DR), SCaSML further reduces the relative I error by 6.6% 10.9%, demonstrating its capability to improve even well well-performing surrogates on complex, high-freqqu figefrreyuency problems.
(a) Illustration of SCaSML's Improved Scaling Law
Surmgett
Faynman Path Simulation
Simubeion-Ciibratee Scientifc Machine Leaningg
Methods
Triining time
Inferrec time
Stating
Law
Error: O,0r7)
Error: 0,07
Error:
(h) Empirical Verification on the Viscous Burgers Equation
Scatll (ym0d3))
Tadatbr
18
(a) d = 20
(b) d = 40
(c) (c)d = 60
(d) d = 80
Figure Empirical verification of the improved scaling law for SCaSML. (a) Conceptual diagram.
The final S error is product of the surrogate error and the simulation error (Theorem 2.5)
By balancing the computational budget between training the surrogate and performing inference- simulation, SCaSML. achieves a faster overall convergence rate (Corollary 2.6). (b) Numerical results. We plot the L error versus the number of collocation points (m) on a log-log scale for GP surrogate and SCaSML. The slope of the line corresponds to the convergence rate 7. SCaSML consistently exhibits a steeper slope than the base surrogate, empirically confirming its accelerated convergence.
4 CONCLUSION AND DISCUSSION
We introduced SCaSML, the first physics-informed inference time scaling framework that integrates surrogate models with Monte- Carlo numerical simulations for solving high- -ish-dimensional PDE. By introducing St strcttrll waw.uraa-rreeeiiiiie
Law of Defect, we use the output of a pre-trained SciMI. solver as an efficient starting point for inference-time corrections. Our theory and experiments show this hybrid approach achieves faster convergence and reduces errors by up to 80% in complex high -igh-dimensional PDEs. SCaSML represents new approach in hybrid scientific computing. Unlike previous work that used machine learning for discovering numerical schemes (Long et al., 2018) as preconditioners (Hsieh et al. 2019), framework uses the machine learning model as a control variate in stochastic immulations to reduce the variance of Monte Carlo simulation. The surrogate handles the low-frequency part, allowing the simulation focus on the small high-frequency residual, and enhances computational efficiency by addressing model bias at inference time. This establishes an elastic compute paradigm, allowing users to trade inference time for accuracy on garnand-achievinng
gains that are often computationally intractable through further training alone. SCaSML is the first inference-time scaling algorithm that enhances the learned surrogate solution during inference without requiring fine-tuning or retraining.
REFERENCES
Randolph E Bank and Alan Weiser. Some a posteriori error estimators for elliptic partial differential equations. Mathematics of computation, 241170)2233-31, 1985.
10

---

Published as a conference paper at ICLR 2026
Pau Batlle, Yifan Chen, Bamdad Hosseini, Houman Owhadi, and Andrew M Stuart. Error analysis of kernel/gp methods for nonlinear and parametric pdes. arXiv preprint 22X.:22305.04962,, 05, 2023.
Roland Becker and Rolf Rannacher. A Afeddbback approach to error control finite element methods: Basic analysis and examples. IWR, 1996.
Roland Becker and Rolf Rannacher.
An optimal control approach to a posteriori error estimation in finite element methods. Acta numerica, 10:1-102, 2001.
Richard Bellman. Dynamic programming and new formalism in the calculus of variations. Pro- 1954. doi: 10. 100.3/pas..00.42333
ceedings of the National Academy of Sciences,
URL https: /tl:://wwwwww /www///www.../w///////////ww pras. htrg///i/aab/// 11t.:/1....00/000
40 231
Jose Blanchet, Haoxuan Chen, Yiping Lu, and Lexing Ying. When can regression-adjusted control variate help? rare events, sobolev embedding and minimax optimality. Advances in Neural Information Processing Systems, 36:36566-36517 2023.
K Böhmer, PW Hemker, and HJ Stetter. The defect correction approach. In Defect Correction Methods: Theory and Applications, pp. 1-32. Springer, 1984.
Klaus Böhmer. Discrete newton methods and iterated defect corrections. Numerische Mathematik,
37(2):167-192, 1981.
James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne and Qiao Zhang. JAX: composable transformations of Python+NumPy programs, 2018. URL
http:/ httt:://gttub.. www. jax ml/jaa
P. Briand and C. Labart. Simulation of BSDEs by Wiener chaos expansion. The Annals of Applied Probability, 24(3):1129-111 2014.
Yifan Chen, Bamdad Hosseini, Houman Owhadi, and Andrew M Stuart. Solving and learning nonlinear pdes with gaussian processes. Journal of Computational Physics, 447:1 10668, 2021.
Yifan Chen, Houman Owhadi, and Florian Schäfer. Sparse factoiiiati for solving nonlinear pdes via gaussian processes, 2024. URL http::://ax..../////atss//
01294
Giuseppe Da Prato and Jerzy Zabezyk. Differentiability of the feynman-kun semigroup and control application. Atti della Accademia Nazionale dei Lincei. Classe di Scienze Fisiche Matematiche Naturali. Rendiconti Lincei. Matematica Applicazioni, 8(3):183-188, 1997.
Alok Dutt, Leslie Greengard, and Vladimir Rokhlin. Spectral deferred correction methods for ordinary differential equations. BIT Numerical Mathematics, 40(2):241-26) 2000.
Weinan E, Martin Hutzenthaler, Arnulf Jentzen, and Thomas Kruse. Multilevel picard iterations
for solving smooth semilinear parabolic heat equations. Partial Differential Equations and Ap- plications, 2(6), November 2021. ISSN 2662-2971. doi: 10.1007744205--01-0008-- URL
http:://dx: doi 10 10000 0210002111-00 021 001-00089--3
K David Elworthy and Xue-Mei Li. Formulae for the derivatives of heat semigroups. Journal of Functional Analysis, 125(1 2:22-288, 1994.
Saadet Eskizzmirllerr Korhan Günel, and Refet Polat. On the solution of the black-scholes equation using feed-forward neural networks. Computational Economics, 58(3)9915-941, 2021.
Michael Giles. Multilevel monte carlo path simulation. Operations research, 56(3):607-617, 2008.
Michael B Giles. Multilevel monte carlo methods. Acta numerica, 24:259-32 1385, 2015.
Antoine Girard. A fast 'monte-carlo cross-validatton procedure for large least squares problems with noisy data. Numerische Mathematik, 56:1-23, 1989.
E. Gobet and P. Turkedjiev. Adaptive importance sampling in least-squares monte carlo algorithms for wark stochastic differential equations. Stochastic Processes and their Applications, 127(4): 1171 1203 2017. ISSN 0304-4149, doi: http:://oo.ogg//....oo//ppppp..0 www. sci enced rect. com/ wttp://cc..aeeeee....e.eeee////a
URL https

---

Published as a conference paper at ICLR 2026
Jiequn Han, Arnulf Jentzen, and Weinan E. Solving high-dimensional partial differential equations using deep learning. Proceedings of the National Academy of Sciences, 115(34):8505-88 00108/85558-111,,,00 8510, 2018a.
Jiequn Han, Arnulf Jentzen, and Weinan E. Solving high-dimensional partial differential equations using deep learning. Proceedings of the National Academy of Sciences, 115(34):8505 8510, 2018b. doi: 10. 10.1073/pnas. http3///aas.111111 1718942115. URL https:/ /ttp:///wwwwww 1073/pnas 171 8942115
/www pnas org/ /000/////000 httg//////b/////////////
Wilhelm Heinrichs. Defect correction for convection-dominate flow. SIAM Journal on Scientific Computing, 17(5):1022-1099 1996.
Jun-Ting Hsich, Shengjia Zhao, Stephan Eismann, Lucia Mirabella, and Stefano Ermon. Learning neural pde solvers with convergence guarantees. arXiv preprint S05i::9900.0200000 06, 2019.
Zheyuan Hu, Khemraj Shukla, George Em Karniadakis, and Kenji Kawaguchi. Tackling the curse of dimensionality with physics-informed neural networks. Neural Networks, 176: 106369, August 2024. ISSN 0893-6080. doi: httff./nuunt.2022. 1016/] neunet. 2024 106369
URL http://dx doi.org/10
Xinru Hua, Rasool Ahmad, Jose Blanchet, and Wei Cai. Accelerated sampling of rare events using neural network bias potential. arXiv preprint arXiv:2401.0660 06936,, 2024.
Michael F Hutchinson. A stochastic estimator of the trace the influence matrix for laplacian smoothing splines. Communications in Statistics-Simulatin Simulation and Computation, 1856 1989.
Martin Hutzenthaler and Thomas Kruse. Multilevel picard approximations of high-dimensional semilinear parabolic differential equations with gradient-dependent nonlinearities. SIAM Journal on Numerical Analysis, 58(2):929 961, January 2020. ISSN 1095-7170. doi: 1
URL http: / www.dooi.oooooooooo..oo 10. 11t.:..11//1111111111111.1111111 1137/77M117011
Martin Hutzenthaler, Arnulf Jentzen, Thomas Kruse, et al. On multilevel picard numerical approxima- tions for high-dimensional nonlinear parabolic partial differential equations and high-dimensional nonlinear backward stochastic differential equations. Journal of Scientific Computing, 79(3):
1534-1571, 2019.
Martin Hutzenthaler, Arnulf Jentzen, Thomas Kruse, Tuan Anh Nguyen, and Philippe von Wurstem- berger. Overcoming the curse of dimensionality in the numerical approximation of semilin- ear parabolic partial differential equations. Proceedings of the Royal Society A: Mathemati- cal, Physical and Engineering Sciences, 476(2244), December 2020a. ISSN 1471-2946 doi: 10.1008/rspa.2019.0 9.060.. http://d: doi www.do.........o w0t.://pa...ap.../. 2019. 0630
Martin Hutzenthaler, Arnulf Jentzen, Thomas Kruse, and Tuan Anh Nguyen. Multilevel picard approximations for high-dimensional semilinear second-order pdes with lipschitz nonlinearities. 2020b.
Martin Hutzenthaler, Arnulf Jentzen, and Thomas Kruse. Overcoming the of dimensionality in the numerical approximation of parabolic partial differential equations with gradient-dependent nonlinearities. Foundations of Computational Mathematics, 22(4))905 9666),99--96,, July 2021. ISSN 1615-3383. doi: 10. 10.10077/10208-021--0000
URL http: //tp:://dx..o wtw.d//.oooooo//0/0.ooooo
s10208- 122208-011000 021-09514-3
Antony Jameson al. Iterative solution of transonic flows over airfoils and wings, including flows at mach 1. Communications pure and applied mathematics, 27(3):283-309, 1974.
George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informed machine learning. Nature Reviews Physics, 3(6):422-440, 2021.
Yuehaw Khoo, Jianfeng Lu, and Lexing Ying- Solving for high-dimensional committor functions using artificial neural networks. Research the Mathematical Sciences, 6(1):1, 2019.
Stig Larsson and Vidar Thomée, Partial differential equations with numerical methods, volume 45. Springer, 2003.
12

---

Published as a conference paper at ICLR 2026
Qianxiao Li, Bo Lin, and Weiqing Ren. Computing committor functions for the study of rare events using deep learning. The Journal of Chemical Physics, 151(5), 2019.
Zichao Long, Yiping Lu, Xianzhong Ma, and Bin Dong. Pde-net: Learning pdes from data. In International conference on machine learning, pp. 3208-3216. PMLR, 2018
Lu Lu, Xuhui Meng, Zhiping Mao, and George Em Karniadakis. Deepxde: A deep learning library for solving differential equations. SIAM Review, 63(1):208-22)) 228, January 2021. ISSN 1095-7200. doi: 10.1 10. 9mm 1274067. BRLhttp::/Gee doi org/ 113./11371111 9111219MM22 274067
Dario Lucente, Stefan Duffner, Corentin Herbert, Joran Rolland, and Freddy Bouchet. Machine learning of committor functions for predicting high impact climate events. arXiv preprint arXiv:1910.11766, 2019.
Nasim Rahaman, Aristide Baratin, Devansh Arpit, Felix Draxler, Min Lin, Fred Hamprecht, Yoshua Bengio, and Aaron Courville. On the spectral bias of neural networks. In International conference on machine learning, pp- 5301-5310. PMLR, 2019.
Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Physics informed deep learning (part i): Data-driven solutions of nonlinear partial differential equations, 2017. URL https
//arxiv. htt.v/orgaabbb//111 10561
Research Group on Stochastic Analysis, University of Duisburg-Essen. Multilevel picard approxi- stwchastissscchttiisitsc.......s scche
mation, 2025. URL https:://// www. aun--uue... Accessed: March 26, 2025.
analysis/mee mlp
Lorenz Richter, Leon Sallandt, and Nikolas Nüsken. Solving digendimensiona parabolic pdes using the tensor train format. In International Conference on Machine Learning, pp. 8998-9009. PMLR, 2021.
Daniel de Souza Santos and Tiago Alessandro Espinola Ferreira. Neural network learning of black- scholes equation for option pricing. arXiv preprint 2405v:2405.05770, 2024.
Sebastian Becker Sebastian Becker, Ramon Braunwarth Ramon Braunwarth, Martin Hutzenthaler Martin Hutzenthaler, Arnulf Jentzen Arnulf Jentzen, and Philippe von Wurstemberger Philippe von Wurstemberger. Numerical simulations for full history recursive multilevel picard approximations for systems high-dimensional partial differential equations. Communications in Computational doi: 1e000//cccp.02222200111111 URL http: www.ooo.ooo.ooooooooooooooooooo. 10. 42080.442220000 wwc.//ccp.......... OA 2020-0130
Physics, 28(5):2109 2138, Jamuary 2020. ISSN 1815- 2406.
Zekun Shi, Zheyuan Hu, Min Lin, and Kenji Kawaguchi. Stochastic taylor derivative estimator:
Efficient amortization for arbitrary differential operators, 2025. URL fittps://axivv.or/
foss/2412.c0m 00088
Charlie Snell, Jaehoon Lee, Kelvin Xu, and Aviral Kumar. Scaling Ilm test-time compute optimally can more effective than scaling model parameters, 2024.
Hans Stetter. The defect correction principle and discretization methods. Numerische Mathematik, 29(4):425 444,,,424443, 1978.
Gilbert Strang, George J Fix, al. An analysis of the finite element method, volume 212. Prentice-hall, 1973.
Michael Unser. A unifying representer theorem for inverse problems and machine learning. Founda- tions of Computational Mathematics, 21(4):941 2021.
Jason Wei, Xuezhi Wang, Dale Schuurmans, Maarten Bosma, Fei Xia, Ed Chi, Quoc V Le, Denny Zhou, et al. Chain-of-thought prompting elicits reasoning in large language models, 2022.
E Weinan, Jiequn Han, and Arnulf Jentzen. Algorithms for solving high dimensional pdes: from nonlinear monte carlo to machine learning. Nonlinearity, 35(1):278, 2021.
Jinchao Xu. A novel two-grid method for semilinear elliptic equations. SIAM Journal on Scientific Computing, 15(1):231-237, 1994.
13

---

Published as a conference paper at ICLR 2026
Shihao Yang, Samuel WK Wong, and SC Kou. Inference of dynamic systems from noisy and sparse data via manifll--oonstraneed gaussian processes. Proceedings of the National Academy of Sciences, 118(15):c2020397118 2021.
2021.
Jiongmin Yong and Xun Yu Zhou. Stochastic controls: Hamiltonian systems and hjb equations. 1999.
URL BtLt::/s://app http:::://ppl.smaaa...... www.nticsccooooe.e........ Org/CorpusID: 1180022771872
Olgierd Cecil Zienkiewicz and Jian Zhong Zhu. The superconvergent patch recovery and posteriori error estimates. part The recovery technique. International Journal for Numerical Methods in Engineering, 33(7):1331-136: 1992a.
Olgierd Cecil Zienkiewicz and Jian Zhong Zhu. The superconvergent patch recovery and posteriori error estimates. part 2: Error estimates and adaptivity. International Journal for Numerical Methods in Engineering, 33(7):1365-1382 1992b.
14

---

Published as a conference paper at ICLR 2026
APPENDIX OUTLINE
This appendix provides supplementary materials to support the main paper. We include a centralized notation glossary, detailed background on the methods used, the complete algorithm, full proofs of our theoretical results, and extensive additional experimental validation including statistical significance tests and computational budget analyses.
The appendix is organized follows:
Appendix A: Notations. We provide centralized glossary defining the mathematical
symbols and operators used throughout the paper and appendices.
Appendix B: Preliminaries. We provide background on the core technical components of the SCaSML framework.
Surrogate Models for PDEs (B.7) : We detail the architectures used: Physics-Informed Neural Networks (PINNs) and Gaussian Processes (GPs).
Multilevel Picard (MLP) Iterations (B2);: We overview the quadrature and full-history MLP methods which form the basis of our inference-ting time correction.
Appendix C: Algorithm. We present the complete SCaSML algorithm in detailed pseu- docode, including practical implementation details like outlier thresholding and Hutchinson's 's estimator for nigh-dimensional Laplacians.
Appendix D: Proof Settings. We establish the common probability space definitions and specific regularity assumptions the surrogate models required for our theoretical analysis.
Appendix E: Proof for Full-History MLP. We provide the theoretical analysis for the full-history MLP variant using Monte Carlo time integration.
Globaal I.2 Error Bound (s E.3.2) ; We derive the error bound for the full-history case.
Improved Scaling Law ( E.3.3
We provide the proof for the accelerated asymptotic
convergence rate of O(m 1/2).
Appendix F: Proof for Quadrature MLP. We present the theoretical analysis for SCaSML
using the quadrature MLP solver.
Global T.2 Error Bound ($
F2.1
We derive the error bound showing dependence on
the surrogate's accuracy e(u).
Commputttional Complexity the SCaSML framework.
F2.2 We prove the reduction in complexity afforded by
Appendix G: Auxiliary Experimental Results. We include comprehensive additional
experiments to validate robustness, statistical sigmificance, and efficiency.
Violin Plots for Error Distribution ($ C.7. Visualizations of the full error distribution for all test cases.
Inference Time Scaling Curves G.2 and G.3 Plots demonstrating monotonic error
reduction with increased inference compute.
Statistical Analysis of L Errors GS.4) Detailed tables reporting means, standard deviations, 95% confidence intervals, and p-values from paired t-tests.
Reelatiee L Error Improvement G.5 Visualization of the percentage error reduction
across all dimensions.
Pointwise Error Reduction Analysis G.6); Scatter plots confirming that SCaSMI. systematically reduces error on the vast majority of individual test points.
Performance Comparison Under Fixed Computational Budgets ($G.7): Pareto effi- ciency analysis comparing SCaSML to baselines when total wall-clock time (Training + Inference) is held constant.
Performance Comparison:
Large PINN vs. SCaSML Correction ($..). A Pareto efficiency analysis comparing SCaSML to PINN with increasing scales the same computing budget.
15

---

Published as a conference paper at ICLR 2026
A
NOTATION
This section establishes the rigorous mathematical framework, including probability spaces, function classes, and norms used throughout the theoretical analysis. We strictly distinguish between spatial functional norms and probabilistic norms to ensure clarity in the convergence analysis.
A.1 GENERAL CONVENTIONS AND GEOMETRY
Tet (0, oc) be a fixed terminal time. We define the spatiotemporal domain fr [0,T] x Rd, where d E N denotes the spatial dimension.
We adopt the unified coordinate convention (1, x)) E Ωtr throughout this appendix. The notation (r, v) is reserved strictly for integration variables within time integrals.
B(RR) denotes the Borel o-algebra on IR
(xr y) denotes the standard Euclidean inner product for æ; y E IRd, and |a]| == V(a,, x)
denotes the Euclidean norm.
A.2 NORMS AND FUNCTION SPACES
Measurable Functions:
Let M(A,B) denote the set of all measurable functions mapping from measurable space A to B.
Spatial Spaces and Norms:
- For any function Φ IRd RR* we define the uniform norm o == supzer [d(x))]-
( T] x Rd) denotes the space of functions that are once continuously differen-
tiable in time and twice continuously differentiable in space.
W7k.00 (Rd) denotes the Sobolev space of functions with essentially bounded weak
derivatives up to order k, equipped with the norm [[ W := Σ l|
.
Probabilistic Spaces and Norms:
Let F. P) be a complete probability space equipped with filtration (F (Fi) satisfying the usual conditions.
For E [1, oo). (P2; R²) denotes the Lebesgue space random variables X Ω I* with finite p-th moment. We explicitly define the probabilistic norm:
| == (E (LXP)))
A.3 PDE FORMULATION AND STOCHASTIC PROCESSES
The Operator: Let μ [0, Rd - Rd and a [0, Rd Raxd We define the second-order linear differential operator L acting on ó E C12 as:
Ca(1, ()),,,)), V ol m)) + Tr ()) Het, 2
The SDE: For any (t, x)) Ωr, let X (X!- )xc(t,TT (X:"")<((()) be the unique strong solution to the stochastic differential equation (SDE):
8 = J. u(r, (tr, + J. o(r, Xi", 8 E [t,
where W is standard d-dimensional Brownian motion under TP.
A.4 DEFECT FORMULATION (TНЕ SCASML ОBJЕCT)
Surrogate and Defect: Let E ( be the surrogate solution. We define the defect pointwise att,, x) u(t, x) x).
The Residual: We define the PDE residual Mr IR as:
0û (t. x) + Lū(t, + F(ū(t, a), TV Vaū(t,z))
ctt, x) = Ot
16

---

Published as a conference paper at ICLR 2026
Modified Nonlinearity F: Let F Ω2 x R R + R be the modified driver defined by: F(u,2,t,a) 2, t,),) F(ū(t, 0, T (V (t, x) + z)) F(i(t, æ), V (tt x)) + <(t,
Crucially, observe that F(0a+1. x) = ctt, æ). This identity bridges the deterministic approximation error and the stochastic driver.
Terminal Defect: 5(x)) == g - x).
B PRELIMINARY
In this section, provide the necessary background the main building blocks of the SCaSML. framework. First, we detail the surrogate models Sais---h-ssiiionnnniin Neural Networks and Gaussian Processes-used generate the initial approximate solution û. Second, we review the Multilevel Picard (MLP) iteration method, the numerical solver we employ at inference time to solve the Structural presservng Law of Defect.
B.1 SURROGATE MODELS FOR PDES
In our experiments, we employ two surrogate models to solve high-dimensional PDEs:a Physics- Informed Neural Network (PINN) and Gaussian Process (GP) regression model. Both models implemented in JAX (Bradbury et al. 2018) and DeepXDE (Lu et al. 202I) to leverage efficient parallelization and runtime performance. Furthermore, Hutchinson's estimator technique as delineated in (Shi et al., (025) is incorporated during the training process to substantially decrease GPU memory consumption, applicable to both the training and inference stages of Physics-Informed Neural Networks (PINN), as well as the inference phase of Gaussian Processes.
B.1.1 PHYSICS-INFORMIE NEURAL NETWORK (PINN)
Neural Networks (PINNs) are designed to approximate solutions of PDEs by embedding physical laws into the learning process. In our framework, the neural network á(t, æ) with parameters & approximates the true solution (t, x) of the given PDE. The training loss is constructed weighted sum of several components, each designed enforce key aspects of the problem's constraints.
The first component is the PDE loss, which ensures that the network output adheres to the governing differential equation. This is achieved by penalizing deviations from the expected behavior defined by the differential operator, evaluated at a set interior collocation points {(4, {x))) The PDE loss is defined
k=
Σ CPDE (0) S
1
S
o² (Lx, mK) + oh Agū(lu, F F(ū, VV,
(10)
Dr
In order to satisfy the prescribed boundary conditions, the model employs Dirichlet boundary loss. This term minimizes the difference between the network output and the given boundary values l((xz)
S2 selected boundary points {(le, 2t& )xx) k- and is expressed
1
(f)
[().. - h(xx)11
(11)
S2
Moreover, the initial conditions of the problem are enforced by an initial loss component. This ensures that the solution time = 0 matches the known initial data g(xx)) for the points ((0,
Ss
1
Cinitial (0)
[ū(0,
qxxx))
(12)
Sa
The overall training objective is then formulated as a combination these losses, with each term scaled by its corresponding weighting coefficient:
(()) an L(0) + 0:2 Lox + 0z Linitial
(13)
This formulation ensures that the PINN not only fits the observed data but also rigorously respects the underlying physical laws, boundary conditions, and initial conditions governing the PDE.
17

---

Published as a conference paper at ICLR 2026
B.1.2 GAUSSIAN PROCESSES
In this section, review the Gaussian Process (GP) framework developed in (Chen al. 2021 Yang et al.
2021; Chen et al., (024) to solve nonlinear PDEs. Consider solving a semi-linear parabolic
PDE
Ou 01 (1, a)) - 7 (u((, æ) An Φ), div, u(t, u)); u 2x = g(z), =
V((,)) E [0, x R"; Yz E R
(14)
where is nonlinear function of the solution and its derivatives, and specifies the terminal condition.
The GP Framework Consider one already sample Mn interior points and Maa boundary points, denoted as Xin { ...,x. x Mas C [0, [] x I2 and Xst {xd ...; ............ x
Mw C {T} R Then,
we assign an unknown GP prior to the unknown function u with mean 0 and covariance function
K: : ((0, TJ x RR)) ((0, T] x IR:)) -> IR, the method aims compute the maximum a posterior estimator of the GP given the sampled PDE data, which leads to the following optimization problem
minimize well
|w|||
s.t.
00 (x = +(((x()))) Aw uxx div, u(xx)), for m = 1,..., Mn= xx] g(x)), for 1,.... Mhi-
(15)
Here, || ll is the Reproducing Kernel Hilbert Space(RKHS) norm corresponding to the ker- nel/covariance function K. Regarding consistency, once K is sufficiently regular, the above solution will converge to the exact solution of the PDE when Mn. Mbi + oc; see (Batlle et al. 2023 Theorem 1.2).
We denote the measurement functions by
2 (u) u >> 8xm u, S m Mn:
2 (u) : Sey o U, S m S Mba-
Ou
3 (u) u > 8xm u S m S Mmn
8 (u) : 6xz s m S Mm (16)
2 (u) : u + 8xm o div, u, < m S Min:
where 6x is the Dirac delta function centered at x. These functions belong to LT the dual space of U, for sufficiently regular kernel functions. We further the shorthand notation o² 23 di de for Mi dimensional vectors and 2² for Ma dimensional vectors as finite dimensional representation for corresponding features. We use [:. to denote the primal-dual pairing, such that for U and ot E u Vi it holds that [u, of = f uxx)o'n (x))xxx For instance, for & we have [u, d.. Ou (.m). Based on the defined notation, we can rewrite the MAP problem 15
((()) (x)dx= DE
as
minimize well
|||||
s.t.
2m (u), 2m Om
(1) (1)
(3)
(3)
(5) (u), o (u), 2m o (u), 1,.... Mn.
(4)
(1) pm (u), m 1,..., Mba-
Am
(4) 7(=m (3)
1, Min:
2m
2m
(2) g(x m 1, Mbd-
2m
(17)
Finite Dimensioanl Representation via Representer Theorem According to Representer The- orem (Chen al. 2021: Unser 2021) show that although the original MAP problem (15) is an infiniteddimensional optimization problem, the minimizer enjoys a finite-dimensional structure
u' (x) K(x, d)a
(18)
where (xx, Φ) the (4M,m + Maa) dimensional vector with entries [((x, x') x')xx (x')dx' (here the integral notation shall be interpreted as the primal-dual pairing as above), i.e.
K (x, Φ) = (((xxx) K (X, Ag ((x,
KK(x.. xin) div. K(x, E R +M
(19)
18

---

Miss Published as a conference paper at ICLR 2026
and E MhM is the unknown coeficients. Based the finite dimensional representation (88 we know
100T =222 2(() ,2( =36)TT] K(d,o)a.
(20)
where 2(1) = (ol( ... (u)] E RM, = (
Mw (u) E
MMi 2(3) [oi(u), (u) [,, (w)] RMn,
[oi(u), (z(w), Ma (w)
RMM, 2(6) = [o)(())." ... S (u)] E RMM, and (()
u) is the kernel matrix as the (4M Maa (4M + Ma) matrix with entries f [x, dm entries of Precisely ((Ø, の) can be written down explicitly as:
where om denotes the
K ((Xins X(n)
K (Xini Noo)
A ((Xis X/n)
K (Xins Xin)
div K K(Xin: Xyx)
81 K (Nhd, Nod)
K (Nnd, xmm)
K (Nhd. Nea)
Agt K ((Nad, x(a)
div K ((rad X
07
Ag K (Nin: x(x)
Dg K(Nin (Nin; x,o)
Ag a K (Niin xnn)
Az
K (Nin: %/n)
Ag div K K(Nin= X(a)
K (Nin X(n)
8 K (Nin: Xaa)
Agr (Niny x(n)
0 8
K (Nin: x(m)
á div K K(Nin: X
01
div. K (Xin; x(m)
div.. K (Xii, %ba)
div, Agr K (Xm:: X
div. ät K(Xi, x(n)
div div. K K(Xim X.n). (21)
Here we adopt the convention that if the variable inside function is set, it means that this function is
applied to every element in this set; the output will be a vector or a matrix, e.g. K ((im: X/n) = exp
I 11
S m, S Min E M MMaXMM in the Gaussian kernel of our numerical experiment, where
2(avaja
is the variance of the equation. Thus the finite dimensional representation (18) can be rewritten in terms of the function (derative) values
ul u((x k(x, (((x,))K ( (1)
(22)
where 2 = (00 (2)
(4) T Mr
20)
Plug the the solution to the following finite-dimensional quadratic optimization optimization problem with nonlinear constraints
representation
22 the original MAP problem (17) we have that 2f is
min 2 (()
ZERAM a MW
subject to
(23)
(4)
(3) (5)
2m
=m
(2) m
= g m
m 1..... Mod-
Solving the Optimization Formulation
To develop efficient optimization algorithms for (23)
observing that the constraints in terms of the other variables, (Chen et al. 2021 2024) reformulate the optimization problem as an unconstrained problem
2(44) (2211)
2m 2m
(3)
(5)
and
2 (xx)) express 2(1) and am
min 200) 203), ER
(2(1) g(xba)) (22)) 2(0)) 2z))) 20 :g0l1/1))((xlaa))) g(xaa): (8); T(2(1,, 2(3)) (8)))
We apply Sparse Cholesky decomposition to the positive-definite (K(o, Φ)) + n/J) as LLT turn, of [K (d), 6) m// -b (LL*) (L- (L-1)) = -1---1 Hence, the loss function is defined as J [zx(1, (00)),(0)) 1 [||2 Optimization is carried out via Newton method in 20 iterations. We initialize (1) 23) 2(3) E TR M. following NO0. 10 10 In each iteration, the
gradient JJ and Hessian V are computed via automatic differentiation. and the Newton direction A= is obtained by solving (v²J Al) A= -VJ, where À 10 is regularization parameter. Then, update J at Newton direction with step size a = 1. Early stopping is triggered when the gradient norm falls below 10-5 Finally, to apply the representer theorem in 22 the algorithm solves the linear system (K ((, Φ) +nI)w zt to obtain the weight vector w and the final PDE solution
given ut (x) K K(x, d)w!.
19

---

Published as a conference paper at ICLR 2026
B.2 QUADRATURE MULTILEVEL PICARD ITERATIONS AND FULL-HISTORY MULTILEVEL PICARD ITERATIONS
Multilevel Picard Iteration (MIP) method (Hutzenthaler al. 2019 is simulation- solver
which solves semilinear parabolic PDEs (Hutzenthaler al. 2019 Han et al. 2018a Weinan et al. 2021.. represented as the following.
Dr t + (M, vy - Tr(o Hess a) F(u y) 0, on (0,T) x RR
2 1 (T, v) ,Wy, IR.
(24)
where T > 0, N. Rd R, [0, RR 1 [0, IR RR. Additionally,
let a be a regular function mapping [0, T) x Rd to a real x d invertible matrix.
The MLP method reformulates the PDE into fixed-poiis problem using the Feynman- Kac formula to represent the solution as the expected value of stochastic process's functional. A Picard scheme iteratively solves this fixed-point problem. The MLP method employs multilevel Monte Carlo approach( Giles 2008). blending coarse and fine discretizations and allocating more samples to deeper iterations to control variance. This strategy ensures computational costs increase moderately with accuracy. According to Feynman- Kac and Bismut-Elworthy
formula(Elworthy & Li 1994
Da Prato & Zabezyk, 1997/. the solution u (u, a Vyu) of semilinear parabolic PDE 24 satisfies the fixed-point equation T(u = x where #:Lip([o, Ti x Rdd, R (+d - Lip((0. T
Rd I +d) ) is defined as
(s)))) =B ((x7*)) (1, [.(65-1)] J.
dW..)
T
+ L E [((v(Xx)))) (1, [[,005,111
t-s S [o(r, (::1)- D dw..) ) dl.
(25)
Here x, and D are defined as
X x [ p(r, (r,x +
Σ
( oy(r, Xx"))
j=1
(26)
D = Inaxd ( (0) (/.
X.")) D dr + Σ 2. X ) dwp.
f=1
where We [0, x Rd is a standard (Fe)xe addatte Brownian motion.
The Feynman-Kac formula gives
20 (s, z) "[ + [ EF (t, (ott (V, Xt,
(27)
Note that o V y 0 appeared on the right-hand side in the fixed point iteration, which necessitates a new representation formula of it be simultaneous with 27 And that Bismut-Eiie Slwoutt-wworthy-Le formula(Elworthy & Li 1994: Da Prato & Zabczyk 1997). which gives
9(x [s)
T-s S'
[o(r, xxx)) D aw,
=E
+ L E F(u [o(t (V uu (t, x² (28)
[0]
/' (o(r, D aw, de,
t-s
20

---

Published as a conference paper at ICLR 2026
Concatenating the solution as u (u, o V yU), can define the iteration operator Φ:Lip([0, T] R R ld) Lip 1 x Rd R as the following
1, -o(s,,)]] J." [olr, DD. dw.
(f(v)) =B o(x7*)
(1,
T-s
+ $ E F(v(t,x(+))) (1, 1.(06s.011 t-s
[ (o(r,x;=)--Dp dvV, )] dt,
(29)
and
28 yield
u ())
(30)
The Multilevel Picard iteration considers simulating the Picard iteartion ug(s,
(+(u==))) (s, E N which is guaranteed to converge to a as k - oo for any s E [0, 0, E
Rd Yong & Zhou| 1999 Page 360, Theorem 3.4). Formally, the MLP method uses MLMC Giles
2008 2015) to simulate the following telescope expansion problem derived from the Picard iteration.
un (8. ui (8. + Σ (1, u((s
= ((u))( [m(,))
+ Σ (u)) $ T (1
(g(x), 0a) +E (o(x)1) g(x))) -((=)) (1, xT)
T [ [o(r, (o(r; D dW, }]
Σ E (F(U(((,))))) 1 (x())(u
1, -ا( S(" (o0r,
T dw.)) dt.
(31)
One can either estimate these integrations with the quadrature method(quadrature MLP (E et al. )02I)) or the Monte- Monte-Carlo metlod MLP (Hutzenthaler et al., (020b)), detialed intruction is shown in demonstrated in [3.2] A comprehensive summary of MLP variants can be found at (Research Group Stochastic Analysis, University of Duisburg-Essee
2025
B.2.1 IMPLEMENTING MULTILEVEL PICARD ITERATIONS
Suppose we are given effective simulators (e.g., Muler-Maruyama or Milstein) parameterized by yp (e.g. discretization level), which produce the numerical approximations
(s,,, xo (1. [[(00.,)]] t-8 [((o(r,
(1.)) kip 1, X2,,,
dw, (32)
where h: denotes the total level, the current level, and i (which may be negative) indexes the sample path. To implement the Multilevel Picard Iterations, we need numerical approximation to the integral f EF EF(u (t, xx""))dtt Following (E et al. 2021 Hutzenthaler et al. 2021), we examine the following two methodologies, using quadrature rule and Monte Carlo algorthhm approximate the integral J (u, X (.")))
Quadrature MLP this approach (D al. 202I), quadrature rules are employed to approximate the time integrals that appear in the MLP formulation. This quadrature-based technique is motivated by the need to efficiently and accurately resolve time integration errors while maintaining the stability of the multilevel scheme. By leveraging well-established Gauss Legendre quadrature, obtain deterministic and high-order accurate approximation that is well-suited to the recursive structure of the SCaSML. algorithm.
21

---

Published as a conference paper at ICLR 2026
Definition B.1 (Gauss Legendre quadrature). For each N, let (= [[-1, denote the d
distinct roots of the Legendre polynomial + n qullalbl 11, b] + IR by
[x22 1)"), and define the function
2x (b a)c (a+ dx, ifa and 21-(a+b) E {cf...... ,. b-a
(b a)c² - (a+ 6)
(f) =
(a+bb)
0g
otherwise.
(33)
The Gauss-Legendre quadrature serve as a fundamental building block to discretize the time variable in the Picard iteration. With these polynomials, one can approximate the time integrals with high-order accuracy while controlling the eтror propagation in the recursive iterations.
(la) Definition B.2 (Quadrature Multilevel Picard Iteration). Let { f N.M.a
S CM(B((O,TT x
R2) 0 B(R x R²) be a family of measurable functions satisfying, for all I,j E N and (8,z) E
(0,7) x Rd, we start with (0,1) (s, (s, O8+1- For n > 0, we define the quadrature SCaSML iteration
M"
U (n,, = =g9(x)) 0a)
M Σ
- g(xx) Ih.0 I/0001111 (,,
(1 Mn-i
(((U))))
n.M (1, x p((i))
k-1:p (8, x, ())) - In(1) F(U M (t, ((((.)
(k (8,
Mn-l
1=0 LE(x.T)
7 (1.3)
k 1- x)
(34)
The use of quadrature in this context is motivated by its ability to yield systematic error control over the temporal discretization, thereby enhancing the stability and accuracy of the multilevel Picard iteration in the simulation calibrated framework.
iimulationcaalibatted
Full-history MLP The full- full-history MLP scheme (Hutzenthaler et al. 2021) adopts Monte Carlo
approach to approximate the time integral led (u, Xx""))ddd instead deterministic quadrature rules with fixed time grids. This modification considerably simplifies error Hnalysis(Hutzenthalee et al. 2020a) and avoids all temporal discretization error.
In the full-history MLP, we employ dime-sample that guarantees unbiased Monte Carlo ap- proximation of time integrals. Let v : S - (0, be collection of independent and identically distributed random variables with density p satisfying TP ((((.)) < <b) = p(s) ds. Consider numeri-
cally approximating the integral 1(fi 5.t f(r) dr with t E (s,T), we construct importance f(1e(1)) 1 (nti)e where is the the e(R(),)
sampling estimator with sample size N:i(f:s,t) N Di=l
rescaled density on (s, T) defined as e(r,s) density s) on (s, via R + (T -s)
and R is the random sample from the
Definition B.3 (Full-history Multilevel Picard Iteration (Hutzenthaler al. 2020a)). Let {ulli)) R.M IJEZ
{
C CM(BI FF be family of measurable functions satisfying.
for all L,j E N and (s,z) [0, T) Rdd we start with U (0,+t) (8, -O4t1- Then, for >
22

---

Published as a conference paper at ICLR 2026
define the full-history SCaSML iteration
M
Un, = [g(x), +
(5,
- g(x)) - (8,
M
n-1
Mm-!
(35)
Σ
Mn Σ
(P X( (s. R))
e(s, R(
1s (
((U (1)
xalil) (R!!)
-h (s, R R!.)))) In s (5 R...).
n.M
here R. R!..0) is i-th sampled time point after t at level which is defined as as (1.1) = s (T - -s) e(1.))
C ALGORITHM
In this section, we describe the complete procedure of Cimuulatio-Caalibatte Scientific Machine Learning (SCaSML) for solving high-dimensional partial differential equations (T.. The SCaSML framework any space-time point (t, æ) can be summarized as follows:
Step 1:Train Base Surrogate. First, surrogate model i is trained to approximately solve the target PDE ( serving as a preliminary estimate of the true solution.
Step iPhysiccs-Iffomee Inference-Time Scaling via the dtructural-preservinn Law of Defect. Recognizing that the defect u == á satisfies semi-linear parabolic equation, termed the Structural-presett
Law of Defect,
L ū + (p. V
Vy + Tr Hess, ūo) + F V, v))== =0, on [O, Rd
ú(T.y) g(v),
on R2
(36)
one obtains an estimate of (() by employing Multilevel Picard iteration, either through quadrature-based MLP (Definition [.22 or full-history MLP (Definition B.3
Step 3:Final Estimation. The final estimate of the solution then given by u(t, x) R att, x) + x(ft.)).
The entire algorithm is detailed in Algorithm 1
We emphasize that the sample-wise iteration in Algorithm T can be substituted by vectorized opera- tions, thereby enabling the algorithm to be applied concurrently to multiple points. These performance enhancements were implemented using JAX and DeepXDE, resulting in a time reduction by a factor of 5x to 10x.
Additionally, methods such as thresholding (Sebastian Becker et al. 2020) and Hutchinson's 'ss esti- mator (Hutchinson 1989; Shi al. 2025) could also be employed within the principal algorithm. Thresholding (Algorithm 2} mitigates numerical instability by methodically "clipping" the defect estimator U, a critical action when the surrogate model yields outlier values or when unbounded growth may manifest during iterative correction phases. Hutchinson's estimator (Algorithm 3 al- leviates the computational and memory demands of PPDE in I by forming an unbiased estimator that necessitates only a subset of second-order derivatives approximating the Laplacian. This partial evaluation only expedites the simulation process but also minimizes peak memory consumption. thus averting out-of- out-of-memory issues.
D PROOF SETTINGS
In the following sections, we establish the rigorous mathematical framework for analyzing the SCaSML method. We proceed in three steps:
1. Notations and Definitions: We define the probability spaces, norms, and function spaces used throughout the proofs.
23

---

Published as a conference paper at ICLR 2026
Algorithm 1 Simulation-Callbraat Scientific Machine Learning for Solving High-Dimensional
Partial Differential Equation
Require: Level и. sample base M, target point (s, a), a surrogate model é, threshold E, (quadrature order Q for using Qudrature MI.P)
1: Train a base surrogate model ú to approximate the PDE solution.
2: Take MLP LaP_Lawwool Lfw__a Defeectfffffsss,.,ee--e-_s__s z, n, M, Q) (1, 0.)) ū(s,æ) estimation of u(s,
3: function MLP LAW OF DEFECT(s, x. n, M, Q)
4: i((s x) ((s,)), a VV,
5: in then
> Start Inference- Timee Scaling via Simulating the St St.uutural --r.sturvvi-esseiinngg Law of De fect.
6:
(n.M.Q(8,3)) + 02+1
7:
return U (2,22) end if
8:
9:
(n, - (g(z), 0,) for = 1 to M do
10:
11:
Sample Feyman-Ka Path (0000 k
T) and Derivative Process I(0,-1 (s, x,
kip
in
(n, t (, Mr (0(220)00000)))) T)) ((0r) 2
12:
13:
end for for
1= 0 do for i = to 1 do using Quadrature MLP to calibrate then Compute Q quadrature points with corresponding weights q (..).)1(1) by B.
15:
18:
for all quadrature points t [s, T do
Sample Feyman-Kac Path x/1,8 x (s, t) and Derivative Process I I11:3) (s, xx,)
according formula 32
ae k- (),
20:
z - J.P.LLawoo of Defect(t, if > 0 then
21:
22:
Zprer +- MLP Law of Defect(t, X (s, t), Q,
23:
AF + F(z) - F(zper)
else
25:
AF F(z)
26:
end if
Ú (,,,QQ(s, x)) + Ún,,L((s, + ( Ma-1
AF I(())) Is,
27:
28:
end for end if
29:
30:
using Full History MLP to calibrate then Sample time step (4,0) e(x, T)
31:
32:
Sample Feyman-Kac Path
x.l.0) k.p
(s, 1,, R(())) and Derivative Process
(1.0) k.p
(s, x. R!!..))
according to formula 32
33:
z + MLP Law. of DefecttR k
(1.) (1.i)
(1.1)) M,,,M,
34:
if > 0 then
35:
Zprv MLP Lap_La_ooff Defectt
1 (8, 1
(ll.))
36:
AF F(z) ((pper)
37:
else
38:
AF f
39:
end ir U Un,M.Q(S,2) Ün.M.((s.2) + a e(a.) Ix- (s, R((())) (1.
AF ((d)) A
40:
41:
end if end for
42:
43:
end for
44:
Ü Un.M.Q(s,a) - Thresholding(s, Ün,MQ(s,
Threshold outliers using Algorithm
45:
return U (,,,.Q(8.z)
end function
24

---

Published as a conference paper at ICLR 2026
Algorithm 2 Thresholding the outliers (Sebastian Becker et al. 2020)
Require: Threshold e, defect estimator Ú
1: function THRESIOLLDNNG(,, Ú)
2: for 5 = 1 to d do
3:
U then
4:
U f-
5:
end
6:
if U < -
then
7:
U. +
8:
end
9:
end for
10: return Clipped Ü
11:
end function
Algorithm 3 Hutchison's estimator for estimating Laplacian (Shi al. 2025
Require: Sample size K, target function
I: function (TTEKK.))
2: Draw K different indices from I,..., d with equal probability 1/d, denoted ......:J8
3: Compute f
4: Compute estimator HTE + d _ D3,5
K
5:
return Laplacian estimator HTE
6:
end function
2. Problem Setup:
We explicitly state the regularity assumptions the original PDE coeffi- cients and the stochastic basis.
3. Surrogate and Defect Properties: We formally define the surrogate model, the defect PDE, and the transfer of Lipschitz properties from the original problem to the defect problem.
D.1 MATUEMATICCAL FRAMEWORK AND DEFINITIONS
In this section, we rigorously define the measure- theorrtic structures, function spaces, and norms required for the convergence analysis. Our framework aligns with the standard stochastic analysis settings found (uutzenthller al., 2020a
al. 2021).
Definition D.1 (Coordinate System and Vector Norms). Throughout this article, we fix a time horizon E (0, oo) and a spatial dimension d E N. We denote the time-space domain by A := [0, x Idd. We consistently use the coordinate notation (t, with t [O,T] and Rd.
For any vector v (U1,..., va) E Rd, we denote the standard Euclidean norm by |o| = (Ei |na 1/1/22 and the inner product by v- For generic vector E IR. (e.g., neural ner- work parameters), we define the discrete p-norm (p [1, 20)) and oc-norm as:
1/p
n
|
Σ
2
and ||=||>0 = max |=/].
1si5n
Definition D.2 (Measurable Spaces and Functions). We denote by B(Rd) the Borel o-algebra Rd For any two measurable spaces (S1,Fi) and (22, we define M(Si, S2) as the set of all measurable mappings from S1 to Sz: M(S,, (f S1 S2 VA (A
When the o-algebras are clear from context (e.g., Borel for topological spaces), simply write MRR RR).
Definition D.3 (Probability Space and L Norms). Let (St, P) be complete probability space. For any measurable random variable E .M(2, I) and E (0), the IP [P())-norm is defined as:
X | == (E (E ([ = ( [(w)) (PP(()
1/p
25

---

Published as a conference paper at ICLR 2026
For the essential supremum norm is defined as:
|X]| == CC 2 D:P([ > =}.
Definition D.4 (Function Spaces), Let D Rd be an open set. For k E N and p E [o] the Sobolev space WPR((D consists of all functions u E (D) such that for every multi-indes a No with aa] I k the weak derivative D exists and belongs to LP(D). We define the norm for W:0.00 (D) as || [ullwx.> (D) = || (ID"u)l..=(D)
Furthermore, let ((O, T x R denote the space of functions @(t, x) that are once continuously differentiable in and twice continuously differentiable in x. This regularity is required for the classical solution and the surrogate i.
Definition D.5 (Extended Real Arithmetic). To handle singularities complexity analysis. we adopt the standard conventions for the extended real number line E R U {-00, oo). Specifically, define 0, 0 00
0, 0 = and VX0 co. For any a > 0 and E R we set co.
00, 0
xc. 26 = 0, and 0 = O.
D.2 PROBLEM SETUP AND REGULARITY ASSUMPTIONS
We now formalize the specific partial differential equation and the stochastic framework used for our theoretical analysis.
D.2.1 STOCHASTIC BASIS
Let T E (0, oo) be the terminal time and d E N be the spatial dimension. Let (2,F,P, P,, be all stochastic processes are assumed to adapted to the usual Filtration (F, (Ft)uejo,T)
To facilitate the Multilevel Picard (MIP) analysis, we assume the existence of family of independent standard Brownian motions. Specifically, let (W(..) L,j E 2} be a collection of independent d-dimensional standard Brownian motions adapted to (Fe) Here, the index I corresponds to the level in the MLP hierarchy, and j corresponds to the Monte Carlo sample index within that level.
D.2.2 THE TARGET PDE
While the SCaSML framework applies to general semi-linear parabolic PDEs, we perform the theoretical analysis the semi-linear heat equation. This corresponds to the generator L with drift and diffusion a = sI, for constant E R {0).
The classical solution C 2 ((0, Rd R) satisfies the terminal value problem:
du (t, x) + Cu(t, x) + (u(t, x) Vau(t, æ)) = 0, (f, E [0, R",
(37)
01
where Lo = a² 2 Av, and subject to the terminal condition u(T, x)) = g(x)).
The nonlinearity : R x Rd - T2 and terminal condition RRd + IR are assumed to be Borel measurable functions.
D.2.3 REGULARITY ASSUMPTIONS
The MLP method achieves dimension-indeppeddent convergence rates under Lipschitz continuity conditions on the problem data. These conditions ensure bounded variance propagation across Picard iterations (Hutzenthaler et al. 2021).
Assumption D.6 (Lipschitz Continuity of Nonlinearity and Terminal Condition). We assume the following:
1. Nonlinearity: There exists constant T. > O such that for all (V1, 21), (v2, 22) and (t, x) E T)) x Idd; [ t,,,, F(22, 22, x)) S L(11 V2l + || 22ll1).
IR R
(38)
2. Terminal Condition: There exists a constant
2 0 such that for all æ, y E Rdd
(g(x)) g(y)l < K [| y1l1-
(39)
26

---

Published as a conference paper at ICLR 2026
D.3 SURROGATE MODEL AND DEFECT PROPERTIES
In this section, we rigorously define the relationship between the pre-trained surrogate model and the defect (error) we aim estimate. We first state the regularity assumptions on the surrogate, then derive the properties of the Defect PDE.
D.3.1 SURROGATE REGULARITY
To ensure the classical defect PDE is well-defined, we assume the surrogate is sufficiently smooth. Let C²2² T] x Rd I) be a deterministic approximation u. To ensure the defect terminal condition is well-behaved, we require the following:
Assumption D.7 (Lipschitz Continuity the Surrogate Terminal), There exists a constant K [0, 00) such that for all n. Rd:
(ū(T - ÜT, 4)1 & | y111.
(40)
D.3.2 THE STRUCTURAL-PRESSENE LAW OF DEFECT
We define the defect ú [0, x Re + R as the pointwise error:
i(l, x)) = u(t, x) - a(),
(41)
The core of SCaSML. is the observation that satisfies a semi-linear PDE of the same structure as the original. We explicitly define the coefficients of this new PDE below.
Definition D.8 (Modified Nonlinearity and PDE Residual). Lere: [0. T) Rd + R be the PDE residual of the surrogate à defined by:
0ù (t. Cū(t, x)) + (()(t o Vxū(t, e(t, x) = al
We define the modified nonlinearity F. x IRd x [0, x Rd - R. for the defect PDE as follows. For any state I. and gradient-state state E Rd at spacetime point (t,
fv. z, 1, x)) == F()(t,x:) u. T (t, x)) + 2) F(u(t, x)) ), o Vaū(t, x:) + c(t, (42)
We similarly define the defect terminal condition g(x) == g(x)) - x))
Lemma D.9 (Structural- Preservvvn Law of Defect). The defect à is a classical solution to the following semi-linear parabolic PDE:
Oü Øt (t, x)) + Lü(t, x) + F(ü. o T Vaūtt, vt,, (0.T) x Rd ū(T,x) = g(x). (43)
Remark D.10 (Why Law of Defect is Easier to Solve). The complexity of MLP depends the magnitude of source term F (Hutzenthaler et al. 2020b Theorem 3.1). Based the Lipschitz continuity of F and the variance pariance-reduction structure inherent to MLMC, Hutzenthaler et al. (2021) shows that the overall computational complexity of MI.P is governed solely by the value of F at the origin. Substituting
and 0 into Definition D.8 [0 0, ett, The "source term" driving the Multilevel Picard simulation for the defect is the residual e, already reduced by an approximate surrogate. If the surrogate is perfect (e - 0), the driving force vanishes, and the variance of the Monte Carlo estimator approaches zero. In our later theorem, we show that the variance of MLP can be controlled by the magnitude of
D.3.3 REGULARITY ESTIMATIONS
MLP complexity is governed by both the smoothness and magnitude of the source term; these factors enter multiplicatively because nonlinearities -via their Lipschitz variance through each Picard iteration. Remark 1.10 that the magnitude component the law of defect be improved using the surrogate. It remains to show that the regularity appearing in the law of defect is no worse than that of the original PDE, ensuring that the refinement does not introduce additional smoothness requirements.
and amplify
27

---

Published as a conference paper at ICLR 2026
Lemma D.11 (Preservation of Lipschitz Constants). Suppose F satisfies Assumption D.6 with
Lipschitz constants L. Then, the modified nonlinearity F satisfies the same Lipschitz condition with the same constants. Specifically. for any fixed (t,, x). and any vectors (01, 21), (02, 22) E. x Rd:
[8 21, x) [62,, 2g, x, S L( 12 || zelli),
(44)
Furthermore, the defect terminal condition ğ is Lipschitz continuous with constants K K K.
Proof. Let w1 = (81, 21) and w2 = (i2, z2). We define the background state vector of the surrogate as Ü = ()) Vz a)). From Definition D.8 the difference is:
((w1,,,)) F(w2, ax F(O + w1) F(O) + F(O + w2) - +c
F( + wi) - -FF + w2).
Note that the shift terms F(Ú) and the residual e(t, x) cancel out exactly. Thus, the Lipschitz continuity of F (Assumption D.6p transfers directly to F:
[F( + w1) F(U + u22) T.] un) (Ú + w2)]l T.|ws 42ll11-
This confirms that I' inherits the Lipschitz constants L. For the terminal condition, since á (",, T| x IZ)) by Assumption D.7 the map á(T, x) is Lipschitz with constants K. The triangle inequality applied ÿ g -(,,)) then yields K K + K.
E PROOF OF FULL--HISTORY MULTILEVEL PICARD ITERATION
This section establishes the theoretical guarantees for SCaSML when the potuctural-preeeviiie Law of Defect is solved using the Full-History Multi- level Picard (MLP) iteration. In contrast the quadrature method, this approach utilizes Monte Carlo sampling for time integration, which relaxes the regularity requirements on the solution.
For the theoretical analysis, retain the setting of the semi-line linear heat equation where μ 0₂ and 0 = SIa for a constant & E R.
E. PROBABILISTIC SETUP AND TIME SAMPLING
To analyze the Full-History estimator, must extend our stochastic basis support random time stepping.
Definition E.1 (Extended Probability Space and Time Sampling)- Let (2,F,P) be a probability space that supports the following independent families of random variables:
(W0-))aee 1. Brownian Motions: A collection (W(-5) J1ajz. of independent d-dimensional standard Brownian motions.
2. Time Step Samples: A collection (1l-.) J.jez of independent random variables distributed on (0, 1) according to a probability density function p (0. - (0, 20).
For our analysis and experiments, we specifically select the density p(s) (1 aass for parameter a E (0,1). This ensures that the cumulative distribution function FT(b) bl-o facilitating efficient inverse transform sampling.
E.2 RELAXED SURROGATE ASSUMPTIONS
key advantage of the Full-History MLP its robustness. Unlike the quadrature scheme, which incurs time discretization error scaling with high-order time derivatives of the solution, the Monte Carlo time integration is unbiased. Consequently, we can drop the higher-order regularity requirement (Assumption F.I Item 3) imposed in [F
Assumption E.2 (Accuracy of the Surrogate Model for Full-History MLP). Let be the solution to the Defect PDE. We assume supte(o.7] ||ū(t, Ilw < co. There exist constants CF.1, CF.2 0 independent of i such that the surrogate error measure e(û) controls the following:
28

---

Published as a conference paper at ICLR 2026
1. Residual Bound ):
sup (1.))
[c(t - CF c(ū).
2. Defect Bound (W1,0));
sup [āi(, (lll. S Cr.2((i) e
However, the singularity of density P requires a specific moment condition ensure finite variance. Assumption E.3 (Integrabilityy of the Residual). There exists p E with p 2 such that for all t E ()) and E p):
So 8912 p(s)-p S ds sup E [((1, x + +oW, -oW.)l < <c.
B E(..T)
(45)
Remark E.4. In Assumption E.3 we explicitly identified F(02+1) with the residual c. This assump- tion ensures that the surrogate residual does not grow too fast in expectation along Brownian paths, and that the time sampling density puls sufficient probubility mass near 0 to counteract the singularity.
E.3
MAIN RESULTS
We now show that, with an appropriately trained surrogate model, the structurallppreeevvii Law of De fect can be simulated with lower complexity than the original PDE. In particular, the error of the full-history MI.P upper-bounded by the surrogate model's error measure e(ū).
E.3.1 SKETCII PROOF
The computational complexity of the MLP solver depends on the Lipschitz constant of the nonlinearity F and the magnitude of the "source terms". The "source term" driving the Multilevel Picard simulation for the defect is the residual already reduced by the surrogate. At the same time, show that the regularity in the law of defect no worse than that of the original PDE, ensuring that the refinement introduces additional smoothness requirements. Combining the previous fact, more accurate surrogate makes the defect PDE "easier" to solve. This leads to our main error bound.
Improved Source Magnitude. Since the source term F(02+1,4,æ) surrogate's 's residual by Definition D.8 Consequently, as the surrogate improves with additional training data, the variance of the Monte Carlo estimator decreases proportionally. Lemma E.5 formalizes this argument.
e(1, a²) is the
2. Complexity of MIP Hutzenthaler et al. [202I) show that the complexity of MLP depends on both the smoothness and the size the source term; these contributions combine multiplicatively since nonlinearities, controlled by their Lipschitz constants, propagate and amplify variance throughout successive Picard iterations. Then we analyze both the magnitude of the source term in Law of Defect and its regularity.
3. Preservation Regularity. By Lemma D.11 the defect nonlincarity F inherits the Lipschitz constants I. of the original F exactly. Thus, the regularity requirements for the ML.P solver remains unchanged, ensuring that the refinement does not introduce additional smoothness constraints.
4. Error Bound. Based the previous inituition, Theorem E.6 bounds the total --2-error as multiplicative form, combining the classical MLP complexity with the surrogate's approximation error. Thus this can leads faster convergence rate if the surrogate's approximation error consistently improves.
5. Complexity Estimate. Substituting the reduced source magnitude into the standard MLP bound yields multiplicative error reduction. Combined with Theorem E.11 we improve ((ds-(2+5) to ((de(ū)2+4 (+6g))22)))
29

---

Published as a conference paper at ICLR 2026
E.3.2 BOUND ON GLOBAL L ERROR
Our proof still utilizes the insight that the overall L error in the MLP mainly hinges on the Lipschitz continuity of the PDE's terminal and solution, as well as the extent of nonlinearity at the origin. We illustrate that the parameter linked to the structural-s wwwtturraappreeiiieee constrained by the surrogate error. Initially, we present lemma demonstrating how the complexity of MLP can be capped by the error assessment.
Law of Defect is
Lemma E.5 (Complexity Estimation via Surrogate Error for Full-History MLP). Under Assumptions suppose p 2 2. There exists a constant Cp > 0 independent of the surrogate
D.6 D.7 E.2 and
such that for all
2:
C' SuPsE((.T)
|| F(0+11,8,8 + +oW, -oW)ll
V max (T' - 1,33K VM
2002(0)
sup (t
2/M
C' (e || x + oW, - oWi)sl
L 2F (M)
< CCr
2
(46)
where the constant C' is defined as:
1 $ C max { 2T r (0m)))) (1 - -a)) max L3 max Tb Tb,
23
r
p + 1
π
2
Proof. We bound the three terms the left-hand side of equation 46 using the L bounds provided by the surrogate accuracy Assumption E.2 We utilize the fact that for any bounded random variable
2, Z (L) S Z| (0)-
Step 1: Bounding the Terminal Condition. As shown in the proof of Lemma R.4 (Step 2), we have:
K <
[T, (T) sup rE|o,7]
[| (R))
(47)
Applying the Defect Bound from Assumption
E.2 (Item 2):
K Cr.
(48)
Step 2: Bounding the Residual Term. Recall that F(0a+1,t,2)) c(t, x). The second term involves the - norm of this residual evaluated along Brownian paths. Since the residual is essentially bounded in space-time:
F(0441,8, + oW, - o oWall 22
= ||<((s + toW -W oWWill LX _ (9)
L 20)
< sup (e(s, a + oW) (W)) weft
< sup (c(s,)))- wERd
(49)
E.2 [Item 1): Applying the Residual Bound from Assumption
sup || ((0(++,,))))) 20, CF.ie(0).
(50)
<E(t.T)
Step 3: Bounding the Defect Norm. The third term involves the - norm of the defect solution and its gradient. Similarly, we bound the stochastic L norm by the deterministic uniform norm:
[[ūs + ooW - -oW)&
< |lū[s )11
L (0)
< | ([w
(51)
Applying the Defect Bound from Assumption E.2 (Item 2):
sup
||ū(s,
-) (lw & Cp
(52)
sE(t.T)
30

---

Published as a conference paper at ICLR 2026
Substituting the bounds equation 48
equation
50
and equation 52 back into equation 46
VV CC e(ū)
CCF. c(ū)
CCr2e(ū)
LHS
v
+
VM
2vM
2
(
ov ToTT +3
C + 2 Cr22 +
C
CF.1 e(ū).
(53)
VM
2
Since M 2 1, we can simplify the coefficient by defining CF == + ((CF, CF.) This proves the lemma.
The above lemma, together with standard error estimates for the full-history MLP, yields the following result.
Theorem E.6 (Bound of Global L Error). Under assumptions D.6
D.7
E.3
and E.2 suppose
p 2 2, a E (2--1)) 1T) E [0, x E RR, 3 = (1-a]lp-
For
UN.M
x) with level
2p
N and sample base M as defined in Algorithm [ it holds that
sup (eft.....t
max
Ü (t, ii,t,, (((t, n
S E(M,N) (Cre(i)),
(54)
(t.
where B)M,N)= [- [<(PN pNN +1
(20) exp 3
VMN-1
Proof. Under assumptions E.3 and D.6 combined with the integrability argument in (Hutzenthaler et al. 2021, Lemma 3.3), the proof of (Hutzenthaler et al., [02T Proposition 3.5) holds. Setting
n this proposition, for all ş (1,..., 1}, , we have
[ (Üx,M((t,)) - úú(t, x(t..z)) |1.e SE(M,N)
max L, 3,3
(55)
VM
C supsejer) ((0,+1)(s, +oW, - -oWil 1
(56)
2
C (e d+1) |i a + oW, oWo)<ll L
2
(57)
Take suP(t..)ejo:TT XR (axx.e(I....d+D)) for the LHS, and note that the RHS does not depend on 5,
we get
sup (1 ( -(((,))(
max
(58)
Tmxxx - 1,,3K
C supsett.T) F(0
toW, oWi)l 1
<E(M,N)
sup ()
VM
2/M
(59)
C' suppsei.m))c(()....)) |i(ss Lu(s + ooW oWi).ll
1F
(60)
2
Substituting the sup term in 58 by Lemma E.5 immediately yields the stated result.
In practice, a common choice for M is [ J. Plugging it in solver w.r.t. N:
E.7 we get the error order of the
Corollary E.7 (Error Order for M = [ )) Under assumptions p (2--1)) x(1-1)) t [0,T), Rd, B (1-a)(p
D.6
D.7
E.3 and
E.2
suppose
2)
It holds that
2p
(Ü (ÜN.M x(t, .l1aa S exp
[Nlog V1-8 ((-B++((11))))) (Cre(i)).
sup (1
max
sEft....
(61)
31

---

Published as a conference paper at ICLR 2026
Proof. First, rewrite E(M, N) in E.6 exponential form:
PN 2 + 1)
(2C)N-1 exp ( (M
N(M,N)
(62)
VMN-I
N
= exx (NN) + N (2g( //
log M
(63)
2
Note that [N2aN] S N28N and that M = |
log M 2 log( (22N - 1) 2 log log(NNaN)
1 = 2,3 log - for 2 21/23 We can simplify
62 to
N (((N) + log(2C) ///
exp
2 log M
(64)
-
S exp
(((N) + N log(20) + 3N
-23 (2,3 log N - 1)
(65)
2
exp
(N(o(1) + (2g( + + (2 log - 1))
(66)
2
exp
log N(-3+((
(67)
Plugging 64 to the conclusion of E.6 we get the result we want.
Corollary E.8 (Improved Scaling Law for M [N23N], D). Under Assumptions D.6
D.7
E.3 and
E.2
suppose that 2 2, E p-2 21-11): (), E TRdd, and define 3
o (p 2) 2p
2
Assume that the error at (t,z) of the surrogate model decays polynomially with respect to the number of training points; namely, e(ū) O(m "), for some ? > 0. Suppose further that m = (d+ 1)5" N28N Then, for all sufficiently large m, the SCaSML. procedure improves the error bound from O(m=)) to O m ( with same points number.
Proof. what follows, we adopt the notation f(m) ~ ((m) to signify that limm-xx g(m) Since m is continuous and strictly increasing function of N, there exists a unique inverse function N N (m). Taking logarithms, we obtain log m log (d + 1) + N (m) log 5 2,3 N (m) log N (m)
f(m)
which follows immediately that log m ~ 23 N(m) log N (m)
Define
log m 23
log(d + + ((m) log
2/
and set x log N, so that the relation x ex holds. The inverse of this equation is given by the Lambert W function, i.e., x (2). Therefore,
log m 23
log(d+ ((() log
N(m) = = w(() z
(())
W
log
log (g((d1)++a N (m) log
28
Since (2) ~ log - log log we can deduce that NV Tog(log log m
2 log log Equivalently,
lng m
N (m) 2/5 log Ing
log 17a
log log log
2,
In contrast to the surrogate model, which uses all m points to achieve an error of O(m>?), the SCaSML method allocates 5 N29N points for training and d5 N28N points for inference (see Foot-
note E.9 thereby yielding an error bound of the form O NN log N (-3 + o(1)) (5 N=2N)-7)=
0 N aN(1+0(1)) m=>).
Substituting the asymptotic expression for N (m), and noting that
Vd + I exp log N(m))
(mm) (m))
(68)
mm
32

---

Published as a conference paper at ICLR 2026
log 5
N ((m)
=vd+ 1 exp(
log m
(69)
log m
log 5
log 5
1
=Vd + eepp(
log ml
+ ol
(70)
log 5
log log m
1 (
=Vd exp(
)) log m)
(71)
log log m
=Vd+ Im - ttoloooooe = O tol log )
(72)
We obtain the SCaSML enror bound 00 m tol ing in (m Hence, for high-dimensional problems where m > and for any fixed y > 0, we conclude that O (m->-++((1) < O m (=>>) thereby demonstrating that the SCaSMI. procedure attains
=
strictly faster rate of convergence.
Corollary E.9 (Error Order for M = [N29N]] 1). Under Assumptions
E.2 E.3 and D.7 suppose
p (2 TT- E (0, R2, 8 (1-a)1p
It holds that
2p
max xe (l.,,)+1) (ÜN,LNERN (t,,)) a(t, x1),
exp (NlogN(BB+ (Cre(a))
sup
(73)
Specifically, this approximator Ü N. |28N requires most [[ J) points for evaluation, as detailed in (Hutzenthaler et al. 2020aa Lemma 3.6).
E.3.3 BOUND COMPUTATIONAL COMPLEXITY
We now define two indicators to quantify the computational complexity of full-history SCaSML:the number of realization variables (RV) and the number of function evaluations (FE).
Definition E.10 (Computational Complexity of full-history SCaSML). We define the following complexity:
Let (RV n,M n.M Jn,M NN satisfy RV O.M and, for all n, M E.
n-1 Σ 1=
RV S dM
M
(1+ RV +Ix(1) RV
(74)
This quantity captures the number of sealar normal and uniform time realizations required to compute one sample of Ú Un,M(s,))
Let (PEn,M)n.Me. N satisfy FE,.M and, for all n. ENN
FFR,M S M
Σ (1 +1N(1) + Ix(1) FFL-I.M)
(75)
This reflects the number of evaluations of I and required to compute one sample of M,, x).
Theorem E.11 (Computational Complexity of hull-history SCaSML). Under assumptions D.6 D.7
E.3 and E.2 suppose 2 2, a E ( p-2 (2(p-1), 2T) and 3 2
(I-a)
2p E (0, g).
For any
N and taking M [N22N J. we have
NNlog N (- (-85+
RV M.M FEN,M Sexp
35 (888+ o(1))
(2+8)
(d (d+1)CCe((i)2+
(ÜN.((t,)) n (t,
sup
max
(1..) SE se
(76)
33

---

Published as a conference paper at ICLR 2026
Proof. From (Hutzenthaler al. 2020a Lemma 3.6), we derive that
RV N,Ar s d(5M)), FFN, S (5M)
(77)
Suppose the maximum error E. To compensate for the (5.M) term in the complexity by the denominator of Theorem E.6 we multiply the complexity by c²+0 and then divide it, and put everything into the exponent:
RV N,M FEN,M
(78)
<(d + 1)(5M)
(79)
=(d + (0(65))=2+t6--0
(80)
PN +1
2+8 (20) I exp BM #
) )
(Cre(ū))2tt(L (d+ 1)(5M)* (d11)(M))
(81)
VMN-1
uu (t, (Uxt))t))-uu(((..))
(2+6)
sup
max
(82)
(1.=)))0TT)xR+ SE (+1.....++)
6 (N- - 1o M + o[N)
exp
N (2+ 6) log(2C) log5) + (2 20+)MM/222 + log M
2
(83)
(2+5)
(d+ (d+1)(Cp((i)2+4
max (UN.u((tx)-u'" -u
sup s
(
62
(84)
-cxp NN(2 (2 (2(2+
log(2C) (2))0)))))( log5 (2+6)0 M2/28
8 220
15 2N
log M + +0(1)
(85)
N
- (2+8)
sup (d+1) (d+1)(Cpe(i)2+t (taa))
u (t,
max
se
(86)
<exp
(2++))) N (N(22 (2 (2 log((C + log 5 5 (2 + 0)3
M2/28_
8 2+0
log M + +o(1)
(87)
N
2
2N
(
(Ux.((t,))- u (t,x))
(2+6)
(d+ (d+1)(Cpe(i))tt
max (t. Se cE(L....d (
sup
(88)
Note that [ < N and << a thus
(2+518 1/29 (2+ a(1 { Therefore, by
2
N
78
exp IN (N(220))0g(2C)
log
(2 + 8)B8 N
M2/28
6
2+ 8
log M + +o(1))
(89)
2N
12+0)
(d ((Cpp(((())2
(ÜN.N((..2) - (t, (
sup
max
5
(90)
6. Sexp (N(2+0))0 + log + a(1 + 5) -
2 2
} log M + +0(1))
(91)
2N
(2+6)
(d 1)CPp((0))2+
(Üx.M(t, - u (t,
sup (t.m) CE ((( (1m.
max
12
(92)
34

---

Published as a conference paper at ICLR 2026
Since M [ log M 2 log(N28N -1) ) = 23 N
- 1 for N 21/28
we can further reduce 89 to:
exp (N((2)
(2+ (N(22 (og( + log 5 + a(1 + 6 - (
( 2+0
log M + +1((1)))
(93)
2N
(2++)
(d+ (d+1)(Crc(i))=)
(UN.M(t, 1 (1,
sup (t.=)
max
(94)
Sexp
N (2++)))og((20) + log a1 + 8
+
0 2+0
- 1
2+
2 log N + o(1)
2
2
2N
2N
(95)
(2+6)
(d (++)(Cre((ū)2=
(UxxM(11,0)) 11 - (1,
max (0,0))/0011) CE (I.....t+))
sup
53
(96)
|Nlog .V -86 (-86+o(1))
=exp
(97)
(24+8)
(d + (C+e)CCe((t)2 (1.))) cE(I....+1)
sup max
(ūN.M(1,))) n (1,2))
-
(98)
The right-hand side of this expression is elearly decreasing for large enough N, and in turn, fi- nite. Hence, quadrature SCaSML boosts quadrature MLP with complexity ((d=-(2+)) corresponding physics-informed inference solver with complexity ((de(t))+6=-(211))
F PROOF FOR QUADRATURE MULTILEVEL PICARD ITERATION
In this section, we present the proof for the Quadrature Multilevel Picard (MLP) iteration method. For simplicity, consider the case where 0 and sl E IR) in the proof. We first establish the mathematical framework and underlying assumptions, then analyze the convergence properties and computational complexity of proposed the error of SCaSML bounded by the product of MIP error and surrogate error. Likewise, the complexity is bounded by the product of MIP error and surrogate error. Both indicate that surrogate models can substantially reduce computational complexity while maintaining accuracy guarantees.
cnnulaaion-calibrater variant. The result shows that
Since the Struct preeservig Law of Defect also semi-linear heat equation, can the quadraaurefful-hissoryy multilevel Picard iteration obtain an estimation Ü(s x)) of 2). In this section, we study the theoretical properties of SCaSML that using Quadrature Multilevel Picard Iteration solve the Struct www.uuraaaeeeeeiieeeeeee investigate the full-history multilevel Picard iteration in the next section.
Law of Defect and we
F.1 SURROGATE ACCURACY AND INTEGRABILITY ASSUMPTIONS
To derive improved convergence rates for the Quadrature MLP, we must quantify the quality of the pre-trained surrogate à. We introduce scalar eтror measure e(ū) : E [0, 20) which serves as a uniform bound on both the PDE residual and the approximation error of the surrogate.
Assumption F.1 (Accuracy of the Surrogate Model for Quadrature MLP). Assumption needed for quadrature MLP builds directly on Assumption E.2 augmenting with an additional higher-order regularity condition required for the quadrature rule.
3. Higher-Order Regularity: To ensure rapid convergence of the time quadrature rules, assume the defect satisfies the following Gevrey-class regularity bounds:
[11,, TV
Voe)
o² Ae) ku
(t,x)
oF
sup kEN
S C..3 e(ū),
(A:133/4
35

---

Published as a conference paper at ICLR 2026
This condition is required only for the Quadrature MLP variant (Appendix F) and is relaxed for the Full-History variant [Appendix(E].
Assumption F.2 (Quadrature Integrability). To ensure the well-posedness of the Feynman-Kan Kac expectations, we assume połynomial growth bounds. There exists p E N such that for the zero vector Oa+1 E R
(0(z) 1+ |
FF a) 1+|
sup WER
sup 1E(0.T).>ERd
00.
(99)
Remark F.3 (Magnitude of Nonlinearity at Zero). Recall from Definition D.8 thtthat that F0 t, e(t, x). Thus, the second term in Assumption effeetivveyy bounds the growth of the surrogate 's residual. In the standard Picard iteration for the defect it, the first iteration is driven solely by this term. A small "magnitude at zero" implies that the fixed-point iteration starts very close to the true solution (zero), minimizing the Monte Carlo work required.
F.2
MAIN RESULTS
We now present main theoretical results, which characterize both the accuracy and computational complexity of our proposed method. These results demonstrate the substantial efficiency gains achieved by incorporating surrogate models into the multilevel Picard framework.
F.2.1 BOUND ON GLOBAL L ERROR
Follows the sam proof sketch as Section E.3.1 the convergence analysis proceeds in two steps. First, we establish "Bridge Lemma" that bounds the complexity-deeeriinng constants of the Defect PDE (specifically the magnitude of the nonlinearity at zero and the Lipschitz constant of the terminal condition) linearly by the surrogate error e(ū). Second, we substitute these bounds into the standard error estimate for Multilevel Picard iterations to prove that the final error is the product of the simulation error and the surrogate error.
Lemma E4 (The Bridge Lemma: Complexity Estimation via Surrogate Error). Suppose Assumptions F22 hold. Then there exists constant Co > independent of the surrogate ú such
0.6
D.7
that:
[11, V ( D (t,
sup (t.))
F(0)L,, + ssup keN
S ((a).
(k1))//
(100)
Proof. We bound each of the three terms on the left-hand side of equation 100 using the surrogate accuracy assumptions defined in Assumption F.
Step 1: Bounding the Residual Term. Recall from Remark F33 that F(0,+1,t,x)) = x).
Applying the - residual bound from Assumption F.I (Item 1):
sup (..
|5 [,)]
sup (t
[e(1,=)) < Co.1 e(ū).
(101)
Step 2: Bounding the Terminal Lipschitz Constant. Let ea be the standard basis vector at index a
in Idd The L norm of K satisfies:
K$
Σ
|| D (T, S
(102)
Using the Defect Bound from Assumption
E.I (Item 2):
K S sup AE(O.T]
[ū Alwi. & CQ))
(103)
Step 3: Bounding the Higher-Order Regularity Term. The third term is directly controlled by the Higher-Order Regularity condition in Assumption (Item 3):
(11, V @ + o A. 1)
sup
S C())
(104)
(1:)3//4
36

---

Published as a conference paper at ICLR 2026
Summing the bounds from Steps 1-3, we define Cq Cq.1 TVT + CQ.3. This yields the desired inequality.
We now combine this lemma with the general convergence theory of Multilevel Picard iterations to state our main result.
F. and F.2 the error of the SCaSML estimator Ü N,N.N with level N, sample base and quadrature order N(as defined in Algorithm 7 satisfies:
Theorem F.5 (Global I,2 Error Bound). Under Assumptions D.6 D.7
sf(1...,++1) - ( EENN (Cq e(ū)),
sup (..7))
max
(105)
where the convergence factor E(N) is defined as:
VNN-3 + (14(4C)
N-1 + + 72N+1
E(N) 7C 1
VNW
with constant C' (ITT+IIVTT(++)))
Proof. We apply the general error bound for Quadrature MLP from (Hutzenthaler & Kruse 2020) to the specific case of the Defect PDE. (Hutzenthaler & Kruse 2020 Corollary 4.7) provides bound of the form:
(t.x)) sef.... (ÜN.N.N(t, - (ü(t, x), o (((t,
sup
max
(106)
(11, V
0 +
(t, x: x)|l
<B(N) x
sup (x.))
(, + sup kEN
(k1)3/4
(107)
Specifically, the second term is the supremum bounded in Lemma F.4 By substituting the result of
Lemma F.4 directly into the corollary, we replace the generic PDE constants with the term Cq e(û). thereby proving the factorization.
L
Corollary F.6 (Asymptotic Error Decay). Under the assumptions of Theorem factor E(N) satisfies the following asymptotic bound N - oo:
E.5 the convergence
E(N) = exp N- + NlogN+OON))
(108)
Consequently, the error decays super-polynomially with respect to the computational depth N.
Proof. We determine the leading order asymptotic behavior log E(N) by analyzing the two summands in the definition of E(N) separately. Recall:
E(N) NC 2N N- INN (14(4C)N-1 + +1)11+112 T2N+1
VNN-3
VNN
=:i(N)
=:{((N
Step 1: Asymptotic of the First Term (NN). Taking the natural logarithm of (N):
N-3
log (N) = log(7 2- + log(2Ce)
2 log
(109)
1 N log N N logg(oC() + log N (3g(3.
3
(110)
2
Observing that N oo, the term 1 log N dominates linear terms ONN), we have:
logT = N log + O(N).
(111)
37

---

Published as a conference paper at ICLR 2026
Step 2: Asymptotic of the Second Term (N)...
We bound the numerator: 14(4C) +1 15(4C) N-1
for sufficiently large N. Thus:
N
log T2(N)= I log (1644C)~11222)
log N
(112)
2
log (4C) 1) + N log( 4C) + 2N log 22 log N
(113)
N log N + N (log(4C) + 2 2og T) 0(1).
(114)
2
Similar to Step 1, the dominant term is - log N:
log (2 - N log N + O(N).
(115)
Since F(N) (N(N) + T2(N), we have log (NN)
1og(2max (7, T2})
log
max (log Tı, lo, T2}. Substituting equation 111 and equation 115
( Nog N +
log (.N) & log
Mlog N + O(N).
2
Exponentiating both sides yields the claim.
F.2.2 BOUND ON COMPUTATIONAL COMPLEXITY
To fully assess the efficiency of our method, introduce two key metrics that capture different aspects of the computational cost.
analyze its computational complexity. We
Definition E.7 (Computational Complexity of Quadrature SCaSML). We define the following com-
plexity measures:
First, let (RNn.M.Q)n.M. C N satisfy RNo.M.Q = 0 and, for all n, M,Q E
Ele (N+,RN,,MQ (IN,., RN-LLMM)]
RND.M.Q S d
(116)
This number represents the total scalar normal random variable realizations required for computing one sample of (x,,).
Second, let (FF J,,M,QEz C N satisfy FF = 0 and, for all M,Q
n-1
FEn,M.Q S M (1 +Ix(1) 1x(1) FEL-IMM)]]
(117)
This quantity reflects the number of evaluations of f and ý necessary to compute of one sample of (n,,.Q(s.2))
These metries provide comprehensive measure of the computational resources required by our method. The first metric, RNN, Q. accounts for the cost of generating random variables, while the second, M.,, captures the number of function evaluations needed.
Theorem E.8 (Complexity of Quadrature SCaSML), Under assumptions any 8 > 0 and all N E N, we have
D.6
D.7
and F.
for
(4+8)
sup (1. CE (E(...,,+1) (Ux...((xx)) 1 (1,2x) u
RN RNN.NNN + FENNNN
max
8(d + ((C=1)(Coe(i))*++ exp
Noog N N( 5 <c..
+0(1)))
(118)
38

---

Published as a conference paper at ICLR 2026
Proof. From established results in (Hutzenthaler et al. 2020a Lemma 3.6), know that for all N N,
RNN,N,N 8dN2N,
FEN,N,N I 8N2N
(119)
We want to use the O(NN/2) denominator in Theorem F.5 to compensate for the N term in the complexity. Suppose maximum error E, and note that N2N (NN/224 < (NN/2)at0, va > we multiply the complexity by g i.e.
(4+6)
sup (t.))E(O,T) XBd CE SE( (E......411) | -UNNNN((.))) u (1,
(RNN.N.N + FEM.N,))
max
<8(d+ 1)N2N ( (2 + +11 +11 N
2N-1EN
VNN-3
+ (14(8(VT + DVTT((++ +) + 4) N-1 1) (4+4)
(Coe(t))**
VNN
<8(d+1)N2N ((24(T+1))* UN (I++))* VN M) (1+6)
Ny
(Cqe(a))
<8(d + 1)(d+1)[Cce(@)tt++ exp Nlog N( N(- + -
(120)
The right-hand side of this expression clearly decreasing for large enough N, and in turn, fi- nite. Hence, quadrature SCaSML boosts quadrature MLP with complexity ((de=(4+)) to corresponding physics-informed inference solver with complexity O(de(ū)*
This theorem provides comprehensive characterization of the computational complexity of our method. The inclusion of the surrogate model error measure e(ū) the complexity bound demon- strates how the quality of the surrogate model directly influences the computational efficiency of our approach. Specifically, a more accurate surrogate model (smaller e(ū)) leads to a lower computational cost for achieving given level of accuracy.
G AUXILIARY EXPERIMENTS RESULTS
We include supplementary experimental results that further validate claims, including detailed еггог distribution plots (violin plots) and additional inference-ting scaling curves for all PDE test cases.
G.1 VIOLIN PL.OT FOR ERROR DISTRIBUTION
In this section, present violin plots of the absolute error distributions for the base surrogate model. the MLP, and the SCaSML method. We uniformly select the test points. By combining kernel density estimation with boxplot-style summaries, these plots capture both the spread and central tendency of the errors. A violin plot exposes the full distribution-inn density, variability, skewness, and nutliers-offering offering much deeper insight into model performance. The width of each violin at a given error level reflects the density of the observations. The results indicate that SCaSML reduces the largest absolute error, lowers the median and produces more accurate points for majority of equations compared to the surrogate and MLP, demonstrating its robustness across different dimensions and equations.
39

---

Published as a conference paper at ICLR 2026
10
(b) 20
30
(d) d
Figure 5: Violin Plot for comparison of the baseline PINN surrogate (black), MLP (gray). applying qudrature SCaSMI. (teal) to calibrate the PINN surrogate on linear canveciin-diiffsion equation for d = 10, 20, 30, 60.
=
20
=
40
d
60
d 80
=
20
d
40
d
d=
60
=
80
(a) SCaSML using Quadrature MLP
(b) SCaSML using full-history MLP
Figure 6: Violin Plot for comparison of the baseline PINN surrogate (black). MLP (gray). applying qudrature SCaSML (teal) to calibrate the PINN surrogate on viscous Burgers' equation equation for 10, 20, 30, 60.
G.2 INFERENCE TIME SCALING CURVE
In this section, illustrate how SCaSMI. enhances estimation accuracy the number of inference- time collocation points increases, as outlined in 2.1 and [22] Our findings indicate that allocating additional computational resources during inference consistently improves estimation accuracy.
G.3 IMPROVED SCALING LAw OF SCASML ALGORIHTS
In this section, we consider the viscous Burgers equation as an illustrative example to demonstrate the improved convergence of SCaSMI. algorithms, as suggested by Corollary E.8
We implemented physics-infone informed neural network (PINN) with five hidden layers, each containing 50 neurons and employing hyperbolic tangent activation functions. Because the number of training number of iterations the PINN, the control group was trained using the Adam optimizer (learning rate 7 x 10 31 0.9, 82 0.99) over iterations set 400, 2 000, 4 0000, 0000, 8 000, and 10 000 (as illustrated along the r-axis). The dataset comprised 2 500 interior points, 100 boundary points, and 160 initial points uniformly sampled from [0. 0.5] [-0.5, 0.5) ensuring that m > 1. To replicate the conditions of Corollary [E.8, the SCaSML group was trained over iterations set [400/ (d 1)]. [2 2000(d 1)J. [1 0000 /2/ (d+1)] |6 (6000/dd+1]] [8 000/(dd
points, m, is proportional
1)J. and [10 000/ 1)]. In addition, we set the inference level as N log m/ (2/2 log log m).
with 3 1/2. Theoretically, SCaSML exhibits improvement y of 1/4 o(1) relative to the control group.
For the Gaussian process regression surrogate model, training was performed over 20 iterations using Newton's method. Due to the increasing inference parameters with m and the consequent GPU memory constraints, was not possible replicate the conditions of Corollary E.8 exactly for the Gaussian process model. Consequently, both the control and SCaSML groups employed identical training sizes, which theoretically does not alter the asymptopic convergence ratefi.e. the slope). Specifically, the training data consisted of the following pairs of interior and boundary points: (100, 20), (200, 40), (300, 60), (400, 80), (500, 100). (600, 120), (700, 140), (800, 160), (900, 180), and (1 000, 200). with the z-axis representing the total number of training points. Again, the inference level was chosen as N = [log m/ m/(23) log log m)] with 3 1/2, and the SCaSML continues to exhibit an improvement of 2 + o(1) relative to the control group.
40

---

Published as a conference paper at ICLR 2026
=
20
= 40
d == 60
d 80
d = 20
d 40
d 60
= 80
(a) SCaSML using Quadrature MLP
(b) SCaSML using full-history MLP
Figure Violin Plot comparison of the baseline Gaussian Process surrogate (black), MLP (gray), applying qudrature SCaSMI. (teal) to calibrate the Gaussian Process surrogate viscous Burgers' equation equation for d 20, 40, 60, 80.
15
d = 100
d = 120
d 140
d = 160
(a) SCaSML using full-history MLP
Figure 8: Violin Plot for comparison of the baseline PINN surrogate (black), MLP (gray), ap- plying qudrature SCaSML (teal) to calibrate the PINN surrogate on LQG control problem for 100, 120, 140, 160.
d = 100
d 120
d = 140
d = 160
(a) SCaSML using full-history MLP
Figure 9: Violin Plot for comparison of the baseline PINN surrogate (black), MLP (gray), applying qudrature SCaSML (teal) calibrate the PINN surrogate on diffusion reaction equation for d 100, 120, 140, 160.
We observe that, for the PINNs, full-history SCaSML achieves near-monotonic error reduction across resolutions (with d ranging from 20 to 80), outperforming quadrature SCaSML, which displays oscillatory behavior at higher dimensions. The Gaussian process-based SCaSML similarly accelerates convergence during training. In both cases, the trajectories generated by SCaSML are generally shifted downward relative to the base models, underscoring its capacity to enhance accuracy without altering the fundamental training dynamics. These findings underscore SCaSML's robustness in diverse settings, ensuring reliable convergence even in sigh-dimensional or non-monotonic scenarios.
41

---

Published as a conference paper at ICLR 2026
65-
Introvemen: (1)
72.5
Inprovementto (1g
60-
ippivmmennt (5g
30-
prpouementt (S
20-
60-
2
67.5-
55-
50
66.0
53-
65
600
26-
40
25-
32-
(a) d = 10
(b) d = 20
(c) d = 30
(d) d = 60
Figure 10: For the linear convection -dffrusion equation, SCaSML for PINNS reliably enhances performance with increased computational resources. Notably, scaling effects are more pronounced lower dimensions, potentially due to the MLP's convergence rate exhibiting a linear dependency on the dimensionality d.
95-
20.5
22.5
ippoysmmnne (5g
26.3
improsmment
60
Improvement
27.
21.5
71.2
Imorovemennt
16:3
65
27.0
20.
25.3
20.
35.0
29.
10
107
10
111
10'
10
15
Evalustion Numbers
Eva lotion Numbers
Evalustion Numbers
Evaluation Kumbers
(a) d 20
(b) 40
(c) 60
(d) d 80
Figure 11: For the viscous Burgers equation, SCaSML with PINN consistently improves performance as the sample size M increases exponentially.
30.0
Insrovemion (8)
23.5
Inprovement (%)
24-
Imgrovanion (5)
11-
Irprosmeent (N1)
33.0-
20
240
32.0
315-
28
28
285-
305-
26.
10
Evalaction Numbers
Eve Numbers
Eosluntion Numbers
Eveluation Nambers
(a) d = 100
(b) d 120
(c) d = 140
(d) d = 160
Figure 12: For the HJB equation, SCaSML with PINN consistently enhances performance with increases in the exponential base of the sample size M. However, the scaling curve plateaus at
M 14. likely due the relatively small clipping range of SCaSML compared to the solution magnitude. In general, a larger clipping threshold permits more outliers, thereby requiring additional samples to mitigate variance and ultimately enhancing accuracy; this trade-off must be considered in light of available computational resources.
1.
20.
Inroowrrre 150
15
morrowomentt.coo 2533
30
10'
1000
30"
12
30'
101
16
30F
Evalustion Numbers
Evaliation Numbees
Evaluation Numbers
Bialuation Numbers
(a) d 100
(b) d 120
(c) 140
(d) d = 160
Figure 13: For the Diffusion Reaction equation, SCaSML with PINN consistently improves perfor- mance as the exponential base of the sample size M increases.
42

---

Published as a conference paper at ICLR 2026
(a) d 20
(b) d 40
(c) d 60
(d) d = 80
Figure 14: We apply quadrature SCaSML to calibrate PINN surrogate for the d-dimensional viscous Burgers equation. All plots employ logarithmic scales on both axes, and the slope denotes the polynomial convergence rate. Numerical results demonstrate that, when collocation points for testing and inference are increased simultaneously, SCaSMI. achieves a faster scaling law than the base surrogate model.
d 20
(b) db)==3 40
(c) 60
(d) d 80
Figure 15:
We apply full-history SCaSML to calibrate PINN surrogate for the d-dimensional
viscous Burgers equation. Numerical results demonstrate that, when collocation points for testing and inference are increased simultaneously, SCaSML achieves faster scaling law than the base surrogate model.
G.4 STATISTICAL ANALYSIS OF ERROR REDUCTION AND CONFIDENCE INTERVALS
In response to reviewer feedback requesting a rigorous statistical validation of our results, we conducted repeated experiment analysis. Unlike the single-run statistics, this procedure accounts for the randomness inherent in both the training process (e.g., neural network initialization, optimizer noise) and the inference process (Monte Carlo sampling seeds).
G.4.1 EXPERIMENTAL DESIGN AND METHODOLOGY
For each problem configuration, we repeated the entire experiment Nreps 10 times with different random seeds. In each repetition, we performed the following steps:
1. Model Training: A new surrogate model (PINN or GP) was trained from scratch (where applicable). Settings are the same with
2. Inference: The baseline surrogate, the naive MLP solver, and the SCaSMI. framework were evaluated on a fixed test set Ntest = 1200 points.
3. Metric Calculation: We computed the Mean Relative 1,2 Error, Mean 1.7 Error, Mean Squared 1,2 Error, and for each run.
From these Nreps repetitions, we calculated the following statisties:
Mean and Standard Deviation: Computed across the 10 independent runs.
95% Confidence Interval (CI): Calculated for the mean metric as [a, b] = [p
1.96
VNr.pa 4 + 1.96
a
Neepe
Paired t-test: We performed paired t-tests to compare the error distributions of SCaSML against the baselines (GP/PINN and MLP) across the repetitions. The null hypothesis is that the mean difference error is zero.
The tables below present the full results. Note that while SCaSML requires more execution time (as expected for inference time scaling), it achieves statistically significant emror reductions (p < 0.001) across all accuracy metrics and dimensions.
43

---

Published as a conference paper at ICLR 2026
Thinhy
Tainkgg Dan
Thaing S2a
(a) d = 20
(b) d 40
(c) d 60
(d) d = 80
Figure 16: We apply quadrature SCaSML to calibrate Gaussian Process surrogate for the d-dimensional viscous Burgers equation. Numerical results demonstrate that, when collocation points for testing and inference are increased simultaneously, SCaSML achieves faster scaling law than the base surrogate model.
Table 2: Statistical analysis for Viscous Burgers (VB) with GP Surrogate (10 repetitions). Compar- isons are pairwise against SCaSML.
20M
40d
Menicc
Method
Mone
Stdd
95% CI
Stat (vs SCiSML)
Mean # Std
95% C1
Stat (vs SCaSMI.)
GP
1/46e-1 +: 2.803
[1.43e11,
1,50c-1]
1=158, p=8e-17
1.8e-1
2: 4.26-3
[1.76-11 1.91e-1]
1=127 1= 1-73.11 padie14
1,2
MLP SCaSML
1.8te- 4: 3
[1.77e-1, 1.8%
p=le-14
2.27e-1
6,40-3
[2.176
2,288
6.16e:2
2:1ee
[5.80e22 6.50e
8.1e 3.1e-
28,55e-2,9.0e-
GP
6.37ee-2
1.4e-3
[6.00-2
7,20--21
1=131, 16
9.4262
1,30-5
19.28e-2,
9.666-111
1=129, p-5e-66
MLP SCaSML
7.62e-2
1.3c33
7. [7.Me-2,
1-81.7, p-3e-14
9.35e
(9.19e22.9.9922]
1480.6 p4to-14
2.49e 1: 6.8c
[2.39e 2.63e:21
4.01e
1.0e-3
2.420-2
3.10-4
[7.20e-3,
8.18c-31
1-79.8 p-4e-14
1.29622
4.3c4
[1.256-2, 1.296-2]
1-110, p-2o-15
L (sq)
M
1.25=22
+: 5 1o-4
[1.13e-2,
1.28e-2]
1-50.1 p3
1.956-2 =: 1.2o-3
[1
10800-2 2. 2.199]
1-46.9 Se-
SCASMI
1.37e-3 + 9.16
[1.25e-3, 1.53e-3]
3. 3.026-3 2.00-4
[2.75e-3, 3.500-3]
60d
80d
Metric
Method
Meon 1: Std
95% CI
Sui (es SCaSML)
Mean :: Sid
95% C1
Stal (vs SCaSML)
GP
2,34-11 1: 6 2o-3
[2.00e-1, 2.4b-1, 2.40e-1]
1=174 p=ste117 1-369 prole11
2.66-1 = 4.6-3
[3.58c-1 2.72e-11
p=2e1_p=2e-16
Rel 2
MIP SCISML.
2,5ta=1 4 8.56-3
2.64c-1] [2.40e-1,
2.766-1
8.26-1
[2.626-1 2.800-1]
p-9-77
1.23e-1
[1334-1, 1.280-1]
1.5361
[1.48o-1, 1.57e-1]
GP MLP SCaSMI
1.26e-1
+: 2.003
[1.23e-1,
1,30e-1]
18 1=222
1.49e11
2.60-3
[1.46e-1,
1.54e=1] p-1e----4--1
1.0le
3.9e
[9.43e22 1.07e-11
1=26.3, p=Se- 10
109e:1 2.9e:3
[L.04e11 1.13ee1
1=35.6, p-
5.96e-2
1: 1.7e-3
[5.68e-2, 6:20e-2]
7.90e-2 1.6e-3
[7.65e-2, 2.22ee)
2.15e-2
6.8c-4
12.05e-2,
2.25e-21
p=8t-11 t-124,
2.89-22
8.12-4
12800-2,3.05e21
p=118,
1.2 (o4)
MLP
2.50e/2 4: 2.1e:3
21
1-27.1. pade- 10
3.076:2
2.2e-3
2,
p=31.8.p-1e-11
SCaSML
6:00e:3 1: 3.2e-1
35.40e33,6,41
9.19e 1: 3.76-4
[8.98e-3. 1.00e:2]
G.5 RELATIVE 1,2 ERROR IMPROVEMENT
In this section, we provide supplementary plots that visualize the relative improvement in L error achieved by SCaSML over the baseline surrogate model (PINN or GP). The percentage improvement is calculated as:
|rrorsumggnel La - ||rrorscasMt.I122 x 100
Improvement % =
|Errorsurogate 1
Prrc -
Fatee
12
Taming
(a) d 20
(b) 40
d 60
(d) d 80
Figure 17: We apply full-history SCaSMI. to calibrate Gaussian Process surrogate for the
d-dimensional viscous Burgers equation. Numerical results demonstrate that, when collocation points for testing and inference are increased simultancously, SCaSML achieves a faster scaling law than the base surrogate model.
44

---

Published as a conference paper at ICLR 2026
Table 3: Statistical analysis for Viscous Burgers (VB) with PINN Surrogate (10 repetitions).
20d
40d
Metric
Method
Mcan + Std
95% C1
Stat (Ns SCISML)
Mean + Std
95% C1
Stat (vs SCaSML)
PINN
1.25e-2 + 3.00-
11-10-2,
1.31e-2]
1-172, pm
4.11o2 1.30-
[4.23-22, 4.67c-2]
1-135, 1-36-16
Rel 1
MLP
8.3 4: 2.2e-
[8.0e 8.78e-3]
1=108 15
2.2-3 3.36-2
1.6e-1
[L.03=1. 1.10e:1] 3.36e223.33-23
p=2e11
SCsSMI.
4.37e-3 2.04
[3,90-33 4.7%o-3]
PINN
5.83e
1.3c-4
6.
1-84.5 14
2.16e:
4.5e-4
12.0%
1-71,7, pa 14 2.222e2
2]
MLP
3/43-2
1.1e3
[3.29e-2, 3.61e-2]
1=91.55 15
4,38e 9.8e-
[4.55e2244555-2
p-3e-pp-e-11
SCASMI.
1/4e-3 1: 4 8e-5
(1.38e-3, 1,55e-3]
1.4ee2 2: 3.5e-
1.45e-2] [1.36e-2,
PINN
5.63e 6
15,00e:5, 6.23e:51
1-57. p-7e-13
1.73e-4
3.3e
18.19--,8..77--
1496.6, p-
L (sp
MLP
2.50e-3
1:
2.882-3]
1-18.00 12
3.29p
2.4e-4
[3.9te
1-51.6, p=2e-11
SCASML
6.900:6 + 9.4e7
[5.55e-6, 8.30e-6)
4.34e4 2.20-5
|4.03e-4 4.06e-4]
60d 95% C1
80d 95% C1
Metric
Method
Meon + Std
Sta (es SC&SMI)
Mean = Skd
Stn (vs SCISMI)
PINN MLP
4.620-2 + 1.163
[4.43--2,
4.75e-2]
1-203, p-91-118 p-75.0 p=76-14
6.59622 4 2.0c-3
[6.17--2, 6.05e-2]
1-168 5-57-17 1-44.8 p-7-112
Rel L
1.17e-1 + 3.26.3
1.25e-1] [1.14e-1,
1.22611 + 3.70-3
[1115-1, 1 1.26-11]
SCISML
3.53e/2 9.5c-4
[3.36/2, 3.65e-21
5.1e22 1.963
2,590e
PINN
2.24e-2 + 4.4o-4
[2.144-2,
2,29e-2]
1-134, 1-
3.20 + 6.60-4
[3.10-22 3.31e-2]
1-182, p-2--17
Z
MLP
4.79e-2
1.76-3
[4.55e-2, 5.04e-3]
1-56.8, p-Mo-13
[4.74-22 5.12e:21
1-50h0),
p-
SCASMI
1.50e-2
1 4.36-4
[1:41e-2 1,55e-3]
2.46072 2: 6.9c-
[2.346-2, 2.58e-2
PINN
8.40e44
[7.72e-4, 8.71e-41
1=128, p=de In
1.76e
8.6e
[L.64e 1.95ee3)
1=136, p-
2.a (ug)
MLP
5
1: 3.5e-1
[4.99e-3,
6.09e:3]
1244.00
12
6.08e-3
3.5e-4
5.377-3.6.35e--
1=36.9. p-4-11
SC5SMI
4.92-4 + 2.20-5
[4.43e-4, 5, 5.18e-4]
1.236-3 = 7,5e-5
[1.13e-3. 1.41e-3]
Table 4: Statistical analysis for Linear Convection -Difuuoionnu (LCD) (10 repetitions).
10d
20d
Metric
Method
Mean d: Std
95% C1
Stat (ys SCaSML)
Mean di Sid
95% C1
Stat (vs SCaSML)
PINN
4.85e:2 1: 2.2e-3
[4.40e2,
5.23e-2]
L-MLL pate- II I=134, 16
5.00e 2: 3.6e-3
[7.80e-2, 9.00 2]
I=39, 1:39.2 1I
Rel 2
MLP
2.30e-1 +: 5.5
2,36e-1] [2.20e-1,
241e-1 2 3.5e-3
2.4e-1] [2.32-1,
1=129 165e116
SCuSMI.
2.77e-2
8604
[2.60-22 2,00--3]
[4846-2, 5.3662)
PINN
3.01e:2 1: 1.1e3
1=35.0, p=35.0.paae--1 II
8.71e22 2.2e-3
[8.24e 2.955e2
pe
MLP
1.68e-1
1.70-3
[1.65e-1, 1.70e-1]
t=285, 1=
2.38e=1
[2.31e-1,
1=241,
SCaSML
1.77e/2 1: 4.3e
[1.71e:2, 1.84e:21
4.676 1.3e-3
[4.50e 2, 4.83e-22
2.20e:3
1.8e44
[1.98e 2.58e:3]
1=26.7, p=
1.29
7.9e-4
[1.15e 1:4le2)
1=36.7, p=te-
r.? (sp)
MLP
4.3ke22 + 1.0e-3
[4.85e-2, 3.162221
1-154, p-le-16
1.02-1
3.1e-3
19:47e-2,
1.05e-1]
1-102, p-
SCASML
7,20e 4 4: 2.3e
[6.82e 7.54e-4
4,50e 2.8c-4
[4.05 3, 3.99ee33
30d 95% CI
60d 95% CI
Metric
Method
Mean 4: Std
Stat (vs SC&SML)
Mean :: Sid
Stat (vs SCaSML)
PINN
1.57e-1 +: 9:56-3
[1.43e11 1.74e-1]
1-19.7, p-1e-08 1-833 p-de-14
285e-1 8.76-3
[2.74-11 2.06e-1]
1-55.2, pele122 1-106, p-3e-15
Rel 2
MLF
2.42e 1 + s.le3
[2.29 2.54-11
2.46e-1
+ 4.8c-3
[2.38c-1, 1,52ee
SCUSMI.
9.50e-2
4.603
[9.13-22 1.070-1]
1.250-1
[1.216-1, 1.30e-1]
PINN
1.8le-1
1: 6.0e3
[1.72e-1, 1.89e-1]
1=37.0, pate-1l
4.90-1 2: 8.2e-3
[4.77e-1, 1=8344 4.59e11
Z1
MI.
R0e--1
1:
[2.80e-1,2.0 2.95e-1
1:316 19
4
[4.14e-1,4244-]]]
t=219.
po4e-11
SCASMI.
1.05-11
+ 8o-3
[1.01e-1, 1.11e-1]
4.76-3 1.96-1
[1 2.076-1
PINN
6.4e-2 9.2e3
[5.50e-2,
8:25e:2]
[=144 p-
4.13e-1 1.7e-2
[2.87e11.4.29-1]
p=55.7.p==e-11
r.2 (4))
ML.
1.52e-1 1 633-3
[1.43e-],
1.60e-1]
1=81.8
p=1e-14
106e-1 7.76-3
[2.96e11 3.20-11
p-4e-16
SCASML
2.4e:2 4: 2.le:3
[2.20e 2, 2.80e 21
7.97e 3.9e-3
[7.36e 8.54ee
These plots directly visualize the 20-80% error reduction claimed in the main text and demonstrate the effectiveness of our correction framework across all test cases and dimensions. The experimental settings are identical to those used for violin plots in Appendix G.T
A30
x3c
850
PM
500R
800
PM
(a) d = 10
(b) d = 20
(c) d = 30
(d) d = 60
Figure 18: Relative L error improvement (%) of SCaSML (full-history) over the bascline PINN surrogate on the linear comvection-diffusia equation for d = 10, 20, 30, 60.
45

---

Published as a conference paper at ICLR 2026
Table 5: Statistical analysis for Hamilton-Tacobi-ellliee (LQG) (10 repetitions).
100d
120d
Metre
Method
Mean + Std
9596 C1
Stat (Ns SCaSML)
Mean + Std
9298 C1
Stat (Ns SCASML)
PINN
9.05022 + 2,50-
[8.57-22,
9.49 2.49)-2]
1-6 5. p-1o-04 1=927, p=
9:136-2 # 2.0c-3
[8.77--2, 9.43e-2]
1-25 9, p-96-10 1=1073, p33--24
1.2
MLP
5,69++0 1: 2.0e2
[5 Rhceee 5,73c-0]
5,50e+0 +: 1.7e-2
5,53c=0] [5,48e+0,
SCaSMI
7:270-2
[6-48-22 9,85o-2]
[6.176-2. 6.68-22
PINN
1.50e
3,0e-3
[1.46e-1, 1.55e 1]
I-7. p-de 05
1.68e1 3.3e-3
[1,62e11 1,74g II
1-146, p=2e (=2531,p=le-2
L1
1.21e+1
1.3e-2
[1.21c+]. 1.21e-1]
26
1.23e=] 1:4e-2
[1.2]c+], 1.23e=1]
SCaSMI
109e-1 1.8e-2
19.68e-2, 1,59p-1]
1.02e-1 2.3e-3
[9.81e-2, 1.05e-1]
PINN
3,70e
[3.34e
1=62 pa2e 04
4.09p
[3,79 4.39
1-75.8, p=óe- 14
L. (x))
MLP
1.16++2 1: 2.8e-1
[1:06e+2.
1.47e-2]
1=1623, 26
1.08++2 :: 3.0e-1
[1.98e+2,
1.49e-2]
1=1576, pale-
SCaSML
2:42e-2 7.0p-3
11916-2,4336-22
2.02622 1.1e.
11.88e22 2:20e-21
140d
160d 95% CT
Motric
Mahee
Mean + Std
95% CT
Sua (es SC&SMI.
Mean = Skd
Soi (ts SC&SMI.
PINN MLP SCaSML
1.03c-1 + 1.805
1066-]] [101e-11
1-377 p-fe-03
1.10e-1 + 2.4-3
[1.66-1, 1.14-1]
1-21.7, p-40-19
Rel L
5.36e40 = 2.0-2
5.402=0] (5,34c00,
1-823, p-3e-23
5,26040 + 196-2
[5,24c10, 5.29c40]
1-570 p=2o-23
8.526 1.562
17.456-2, 1.15e-11
8.4He 5.0p-3
17.88e-2. 9.496-21
PINN
1.92-1 + 3,533
[1877=1, 1.99c-1]
1-4.3, p-26-0
2.11e11 + 3.30-3
[2.05e-1, 2,176-1
1-260,p-9ta-1 p-
L1
1:23c41 1
1.25e=1] [1,22cs1.
1-806 p-26-23
1:23cs1 1.5e-2
[1,23e+1, 1.25e=1]
1-2170, p=5o-27
SCHSMI
1.41e-1 = 3.8e2
1.36-1 1.0e-2
[1.32e11, 1.61e-1]
PINN
5.55e
[5.35e 2, 5.89
1-34. p=3e 03
6:00e
[6.15e
2]
1-25.3, p=le-09
122 (sg)
MLP
1.50++2 + 3.1e-1
[1.50e+2.
1.51e-2]
1=1485. p=
1.52+22 :: 3.0e-
[1.S1e+2. 1.52e-3]
1=1583, p-86-26
SCASMI
5 1.00-222 = 1.5e-2
[2.90e-2, 6 74c-2]
3.5622 2: 4.76-3
[3.45--2,488 4.8762]
Table 6: Statistical analysis for Diffusion- (DR) (10 repetitions).
100d
120d
Metric
Method
Mean d: Std
95% C1
Stat (vs SC&SML)
Mean di Sid
95% C1
Stat (vs SCaSML)
PINN
9.83e-3 + 2.604
(9.45e-3, 1.022-2]
1-14.3, p-2o-177
1.10e22 + 3.00-4
[1044-2,
1.13e-2]
1-25.9, p-9e-10
Rel 2
MLF
8.62e-2 + 2.0-3
9.05e-2] 18,20-22,
1-92.6, p-le-14
9.01e22 # 1.20-3
[8.78-22 9.17e-21 [9.66-3 1.03e-2]
p-1e-11
SCuSMI.
9.10e-3
[8630-3, 9.61e-3]
1.00022 2.8044
PINN
1,20e-2 # 12-4
[1.14e22,
1.240-2]
1-14.1. p-26-77
1.36-2 2: 4.0c-4
[1.29c-2,
1.41o22]
1-167, p-42-15
Z1
MLP
9.37e-2
2.-
[8.97e-2 9.77e-3]
t=100 15
9.776-2
19.5le-2.100ee
E=142
SCASMI
1.14--2 + 4o-4
[1.07e-2 1,19e-2]
1.24672 2.8o
1.286-2] [1:206-2,
2,48e44
1.3e-5
[2.29e44, 2.68e-11
p=2e-07
3.0e-4
[2.79g 3.29e44
1=24.3, p=2e
Z (og)
MLP
1.9ee-2 +: 1.le3
[1.73e-2,
2.08e:2]
1=53.1, p=le-12
2.09-2
4.9c4
[1.98e-2,
2.16e2 p=
SCASML
2.17e-4 + 1.365
2,35e-,
2.5te4 1.46-5
[2.38c 4, 2.722
140M 95% CI
1604 95% C1
Metric
Method
Menn + Sm
Stat (vs SC&SMI)
Mean = Sad
Sun (vs SCASMI)
PINN
3.23e-2 + 5.4e4
[3.14e-2, 3.34e-2]
1-47.4 p=de-12
3.59e2
8.104
[3.47622, 3.71e22]
1-72.4, p-9--41 p-96-14
Rel 1,2
MLF
8.96e + 2.3e33
9,47e-2] 18.09e-2,
1-50 p-le-14
8,74e2
255-3
9.17e22 [8.28e-2
1-62.11
SCISML
3:00 4.8e-4
[2.92e 3,0%
3.37e 8.20-4
2.
PINN
4.05e-2 + 8.0c-
[3.91e-2, 4.21e-2]
1-3.5 p-9t-11
4.50e-2 + 0.524
[4.36e-2,
4.67e22]
1-35.5, p-S-111
1,1
MLP
9.89e/2
2.3e-3
[9.53e-2, 1.Me-11
p-2e-14
9.61e22 2.30-5
100e-11 19.34--2,
1-62.9
p-3c-13
SCASMI
3.776-2
6.56-
[3.64622 9.88o-2
4.226-2 =: 9.904
[4.06072,4 4.396-]]
PINN
2.67e-3 8.76-5
[2.53e-3,
2.85--3
1-4188 p-1e11
3.3103 1.504
[3.100-3,
3:53e-3
1-6899 p-l--3
L (sq)
ML.
2.06--2
1.003
[1
11.03-2,,.33
2.50e=]]
I=55.5
p=le-12
1.96-2
2 1.le-3
[1.776-2,
1=16.99
1=
SCASMI
2.31o-3 + 7.0p-5
[2.10e.3 2.44e-3]
291o-3 = 1.404
[2.70033 3.10e3]
G.6 POINTWISE ERROR REDUCTION ANALYSIS
To further investigate the robustness of our method, we present scatter plots visualizing the pointwise error difference between the baseline methods (Surrogate and Naive MLP) and our proposed SCaSML. The settings are still the same with Appendix G.1
For a given test point x, we calculate the difference in absolute error:
AErrr((z) [rrrrraaceiine (x)) [rrrrscasMI [=)
In the following figures:
Red points (AError > 0) indicate locations where SCaSML has lower error than the baseline.
Blue points (AError < 0) indicate locations where SCaSMI. has higher error.
We provide comparisons for both baselines: Surrogate vs. SCaSML (showing the correction of the initial model) and Naive MLP vs. SCaSML. (showing the benefit of using the surrogate as control
46

---

Published as a conference paper at ICLR 2026
= 20
d= 40
d=
60
d 80
20
d
40
d
60
80
(a) SCaSML using Quadrature MLP
(b) SCaSML using full-history MLP
Figure 19: Relative 2 error improvement (%) of SCaSML over the baseline PINN surrogate on the viscous Burgers' equation for d 20, 40, 60, 80.
d == 20
d
40
d=
60
80
d = 20 d
40
60
= 80
(a) SCaSMI. using Quadrature MI.P
(b) SCaSMI. using full-history MI.P
Figure 20: Relative L error improvement (%) of SCaSML over the baseline Gaussian Process surrogate on the viscous Burgers' equation for d 20, 40, 60, 80.
Mare PRbiee OHP
12:23
616
146
6004
1206
d = 100
d = 120
= 1410
= 160
(a) SCaSML using full-history MLP
Figure 21: Relative 2 error improvement (%) of SCaSML (full-history) over the baseline PINN surrogate the LQG control problem for d 100, 120, 140, 160.
696
1009
3006.
1006.
100
= 120
140
= 160
(a) SCaSMI. using full-history ML.P
Figure 22: Relative 7,2 error improvement (%) of SCaSML. (full-history) over the baseline PINN surrogate on the diffusion reaction equation for d = 100, 120, 140, 160.
variate). Across all experiments, the dominance of red points confirms that SCaSML systematically improves accuracy locally across the high-dimensional domain.
47

---

Published as a conference paper at ICLR 2026
2012
-128
10
d = 20
d 30
d = 60
(a) Baseline PINN SCaSMI.
Sparid reviiie
spaiia Droeesion
Spatia Beesdee
10
30
(b) Naive MIP vs. SCaSMI.
Figure 23: Pointwise enror differences for the Linear Convection- Soavection-Diffussion equation. SCaSML outperforms both the pre-trained PINN and the naive MLP solver across all dimensions.
-100
THAPOS
Spata
Scatal
Sprte
20
d 40
d d= 60
d 80
Baseline PINN SCaSML
47M
PROPOE CME Poraage ADC
= 20
d =
40
d 60
d = 80
(b) Naive MLP vs. SCaSML
Figure 24: Pointwise error differences for the Viscous Burgers' equation (PINN Surrogate). We observe that SCaSML corrects the PINN's error (top) and significantly outperforms the standalone MLP (bottom).
G.7 PERFORMANCE COMPARISON UNDER FIXED COMPUTATIONAL BUDGETS
central question regarding inference-time time scaling is whether the performance gain is simply a result of increased wall-clock time, if the SCaSML framework utilizes computational resources more efficiently than standard training. To address this, we conducted a Fixed Computational Budget analysis.
48

---

Published as a conference paper at ICLR 2026
THO7 -- D
26
22 25
23
20
40
d 60
d 80
(a) Baseline GP vs. SCaSML (Full-Hiistrry)
Spsta Dircedenn
Sadiad
20
40
d
60
d 80
(b) Naive MLP vs. SCaSML -Fulo-Histor))
Figure 25: Pointwise error differences for the Viscous Burgers' equation (Gaussian Process Surro- gate).
3154
Spatud Devenson
Spatel Demonie
Speiial Dimorston
d = 100
d == 120
d= 140
160
(a) Baseline PINN vs. SCaSML
1906
te 1200-4
Pe 130 Pedlin Hattaal houtte
Volee 11M 1007
CROOES what CoA Speni Dreanakm
d = 100
d
120
d 140
d = 160
(b) Naive MIP vs. SCaSMI.
Figure 26: Pointwise error differences for the LQG control problem. The contrast in the bottom row highlights that the naive MLP fails in high dimensions, whereas SCaSML (stabilized by the surrogate) performs well.
We define a "unit budget" based on baseline number of training iterations (e.g.. 2,000 iterations for a PINN), other settings are still the same with G. We then scale this budget by factors of x1, x2, x16. For each budget level, we compare three allocation strategies:
1. Pure Training (Baseline PINN): The entire time budget is allocated to training the neural network. A budget of x implies training for k x Nhuse iterations.
49

---

Published as a conference paper at ICLR 2026
Cheee'S MPNCC oiAAAnu HMA Spet Dmarain
Soslidd
spita Drvesin
-2 -0 20 21
d = 100
d = 120
d= 140
d = 160
(a) Baseline PINN vs. SCaSML
Theet -n10 GROAR Spdtal Spetal/ Demeraa2
-000
Spatial (etvinaa
d 100
d
120
d 140
d = 160
(b) Naive MLP SCaSML
Figure 27: Pointwise error differences for the -eaction equation.
2. Pure Simulation (Naive MLP):
The entire time budget is allocated to generating Monte Carlo paths for the MLP solver.
3. Hybrid Allocation (SCaSML); This represents proposed strategy. We allocate small fraction of the budget (specifically 1/(d+ 1)) to training a "weak" surrogate, and allocate the remaining majority of the budget inference-time correction via the preseural- Law of Defect.
This setup ensures a fair comparison where all methods consume approximately the same total wall-clock time (Training Time + Inference Time). We performed this analysis the Linear Convection-Diffuue using the full- full-history SCaSML variant.
(LCD) equation (d 10, 20) and the Viscous Burgers (VB) equation (d = 20)
The results are visualized in Figure 28
PRNN
MPP
SCI5VL
SCASVE
Relative 122 Ermor
Relative 12 Erao
Relative
xx10
2x10-1
Computing Budget (xbaseline)
Computing Budget (x (xbaseline)
Computing Budget [x [aaeeline))
(a) LCD (d 10)
LCD (d = 20)
(c) VB-PINN 20)
Figure 28: Error Computational Budget. The x-axis represents the total computational budget multiplier (log scale), and the y-axis represents the Relative L Error (log scale). SCaSML (teal triangles) consistently achieves lower error than the Baseline PINN (black circles) and Naive MLP (gray squares) for the same total cost.
As shown in Figure 28 the SCaSML error curve consistently lies below the PINN training curve. This confirms that allocating marginal compute to inference-time correction yields higher return on investment (ROI) than allocating it to further training.
50

---

Published as a conference paper at ICLR 2026
Specifically, for the Viscous Burgers equation (d 20), observe that training the PINN for significantly longer (moving right on the x-axis) results diminishing returns due to the optimiza- components. In contrast, SCaSML leverages the rigorous convergence rate of the Monte Carlo correction to reduce error rapidly. This empirically validates our theoretical claim that the hybrid ML+MC scaling law (0 (0-1/1))/)) is superior to the pure ML scaling law (0(m-7)) 7)).
tion difficulty of high-frequency
G.8 PERFORMANCE COMPARISON: LARGE PINN VS. SCASML CORRECTION
A critical question in SciML is whether the computational budget is better spent on training a larger, more expressive neural network (increasing model capacity) or on post-hoc inference-tine correction (SCaSML). To address this, we conducted a second Fixed Budget analysis where we scaled the model architecture while keeping the number of training iterations fixed.
We define "unit budget" (B = 1) corresponding to our standard PINN configuration: fully connected network with width Wasse = 50 and depth Dosse = As the budget B increases by factors of x x2, x4, we scale the network architecture to increase its capacity. Specifically, the scaled width Wg and depth DB are defined as:
WB Whose VBJ, DB max Dbases
Dbase + (B)J).
(121)
This scaling strategy ensures that the network's 's parameter count and computational cost per iteration grow with the budget, allowing us the limits of model capacity.
We compare three strategies under these scaling rules using the Linear Convection-Diffuuiion (d = 10, 20) and Viscous Burgers (d 20) equations:
1. Large PINN (Model Scaling): We train the scaled network architecture (WB. DB) for fixed number of iterations (Niter = 2000). The optimizer is Adam with a learning rate of 10-4 and Glorot normal initialization. The increased computational cost arises entirely from the more expensive forward and backward passes of the larger model.
2. SCaSML (Inference Correction): We employ the SCaSML framework where the surrogate backbone utilizes the available budget. Crucially, the method allocates resources to the inference- inference-time Monte Carlo correction (using the full-history MLP solver with basis M = 10 and levels N = 2) rather than relying solely on the surrogate's capacity.
3. Naive MLP (Pure Simulation): The entire time budget is allocated to generating Monte Carlo paths for the MLP solver. serving as a pure simulation baseline.
The results (Figure 29 demonstrate that simply increasing the PINN's capacity yields diminishing returns; the model hits a "data efficiency wall" where additional parameters do not translate to proportionally lower errors for high-frequency defects. In contrast, SCaSML consistently achieves lower error for the same total compute time, proving that inference-time correction is more efficient user of marginal compute than model scaling for these high-dimensional problems.
FRN MUP
13-1
PIN
NLF
NU
SCaSUb
Relavtee 12 Emar
SC35VL
lativee Emar
SC35ML
Relsatie
12-1_
13
1.0x
10x
2.5x
Computing Budget (xbeseline)
Computing Budget (xbeseline)
Computing Budget Ix (xbaseline)
(a) LCD (d = 10)
(b) LCD (d = 20)
(c) VB-PINN (d = 20)
Figure 29: Large PINN vs. SCaSML. Comparison of eттоr rates when the computational budget is used to scale up the PINN architecture (black) versus performing inference- time correction (teal).
51


Reviews:
[{"review_id": "plVi0TxaRl", "reviewer": "Reviewer_RnQk", "summary": "This paper proposes to start from an initial solution obtained by neural operators and then solve the residue as a new PDE solving problem. This paper uses two variants of MCMC to solve the new PDE. Theoreical results are provided to show the proposed method has a faster convergence rate. The paper conducted experiments on multiple datasets and achieved the best accuracy among all the baselines.", "questions": "1. Is there any prior work that uses MCMC, deep learning, or traditional numerical solvers for correcting the initial solutions? If there exists some, the authors should discuss and possibly compare with them.\n2. Can the authors justify the improvement drop in Fig. 3?\n3. Have the authors tried variants of PINN and what are the errors?", "limitations": "", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. By formulating the residual part as another PDE, the paper further improves the PINN and GP based surrogates without training any new networks.\n2. The paper theoretically shows the proposed method has a faster convergence rate.\n3. The proposed method is compared on multiple datasets and achieves the best accuracy among all the baselines.", "weaknesses": "1. The idea of correcting the initial solution is not new. The paper chooses to use the existing MCMC variants (quadrature MLP and full-history MLP) for the correction, which are also not new. \n2. It is not clear to me how the paper stands out from the existing works, where correction is applied to the initial solution for improvement.\n3. The experiments are not convincing for that the baselines (PINN) is not tuned enough. A fixed architecture, a fixed optimization strategy, and a fixed number of points are used. To show the effectiveness of the proposed method, the authors can consider starting with a stronger initial solution. More layers and more collocation points will lead to more runtime, and it is interesting to see how that part scales with the complexity compared with the proposed method. For example, for SR, the runtime is about 1 - 2 seconds, while for the proposed method, the runtime is nearly 10 times. I am curious under the similar budget, what will make the proposed method stand out.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "AgrX3rgL8x", "reviewer": "Reviewer_Kghx", "summary": "The paper proposes a method to estimate and correct the error predicted by a surrogate model when solving a PDE equation. \n\nThe method applies to semi-linear parabolic PDE equations and is based on recent results on Multilevel Picard (MLP) solvers.  \n\nThe paper describes the approach, and in the Annex provides extensive details of the approach. \n\nThe paper evaluates the approach on a few tasks: linear convection-diffusion equation, viscous Burger equation, high-dimensional Hamilton-Jacobi-Bellman equation, and diffusion-reaction equation. \n\nIt would be nice to highlight the computational or sample complexity of the method (there is a paragraph in 2.3, maybe a clarification of the impact on the specific cases in the annex would be nice).\n\nThe paper mentions in the first paragraph the use of the method in various applications:\n\n1. Schroedinger equation in quantum many-body systems, \n2. nonlinear Black–Scholes equations in finance, and \n3. the Hamilton–Jacobi–Bellman equation \n\nwhile there are experiments for the last, it would be nice to understand if the application to the other two cases requires some modification to the approach. \n\nIt is not completely clear what is observed. From eq. (2), the forcing term f(r,y) is used to define dht error. It is not clear what is the connection between the error (or residual) and the defect. We have access to the prediction $\\hat u(r,y)$, what is also observed?", "questions": "Mentioned before:\n\n1. Clarify the applicability to other problems. \n2. Give an intuitive explanation of what the expectation operator $\\Phi$ of Eq. (8) does\n3. possibly a complexity (sample, computation) analysis (probably what is described in line 264 and 1662, a summary of the results would be nice).", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "The paper presents a post-prediction method that improves the accuracy by solving a PDE, that describes the \"defect\", which is then computed by solving the \"law of defect\" or the error PDE.", "weaknesses": "Only clarification.\n\nIt would be nice if the authors would release the code, since it would be difficult to build on top of this work.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "hJOHuFAHsc", "reviewer": "Reviewer_Wcpd", "summary": "The paper discusses inference-time scaling for machine-learning (ML) solvers for high-dimensional partial differential equations (PDEs). Concretely, the idea is to train a physics-informed neural network (PINN) or Gaussian process (or any other PDE solver that yields a continuously differentiable solution), derive the defect correction PDE which describes the error of the approximation, and solve this error PDE with a Feynman-Kac-based solver at inference time to correct the PINN prediction.\nThis chaining of solvers improves the accuracy order: the final error is the product of the individual errors. \nThe resulting algorithm is benchmarked on a range of PDEs in up to 160 dimensions, including a Hamilton--Jacobi--Bellman problem.\n\n\n**Summary of my recommendation:**\nI appreciate the idea of inference-time corrections to PINN predictions, and I definitely consider it promising. However, as is, I recommend rejecting this work because the defect correction approach is not as novel as the submission states (details below) and because the submission's experiments have some gaps that need to be filled before publication. I believe that if a future version of this manuscript embedded itself better into the literature of numerical solvers and if the experimental gaps were filled, it would be stronger, and I would probably have given a different score. But I believe that these changes are beyond the scope of a revision during the rebuttal period, which is why I give a \"reject\" score.", "questions": "The following questions don't affect my score, but I think including their answers in a revision would improve the submission.\n\n- Line 038f suggests that ML-based solvers don't suffer from the curse of dimensionality like traditional solvers do. This statement needs to be supported more strongly.\n- Lines 061f: PINNs are not surrogate models. A surrogate model is an approximate model used to replace an expensive simulation. PINNs are (at their core) PDE solvers. I suggest replacing \"surrogate model\" with \"ML-based PDE solvers\".\n- Experiment in Section 3.1: what's the boundary condition? There are boundary collocation points, but the problem statement does not mention boundary conditions.\n- Line 126: How costly is it to compute the full Laplacian in 160 dimensions? \n- Lines 365f: Defect correction methods are sometimes unreliable if the residual is only approximate. How problematic is it to use Hutchinson's estimator for evaluating Laplacians here? For reference, Lines 425f suggest that inaccurate Laplacians are problematic for the experiment in Section 3.4", "limitations": "", "rating": 2, "confidence": 4, "soundness": 2, "presentation": 1, "contribution": 2, "strengths": "Even though the idea of chaining two (or more) solvers, like the submission does with PINNs/GPs and Feynman--Kac simulations, is not new, its application to training-vs-inference stages in physics-informed ML is novel. Inference-time scaling had a significant impact on LLMs (though I'm not an expert on LLMs), and I wouldn't be surprised if it became standard in future versions of physics-informed ML solvers as well. \nFurthermore, I appreciate the mathematical rigour of the analysis and the fact that the experiments cover a range of problems, including linear and nonlinear ones, as well as truly high-dimensional problems (dim $\\gg$ 100).", "weaknesses": "Despite the submission's strengths, I think three weaknesses dominate, which is why I ultimately recommend rejecting this work: \n\n1. What the paper calls the \"law of total defect\" (Section 2.1, Definition 1) is a defect correction approach, and defect correction methods have been studied in numerical analysis for a long time; for example, see:\n\n    > Stetter, Hans J. \"The defect correction principle and discretization methods.\" Numerische Mathematik 29.4 (1978): 425-443.\n\n    > Dutt, Alok, Leslie Greengard, and Vladimir Rokhlin. \"Spectral deferred correction methods for ordinary differential equations.\" BIT Numerical Mathematics 40.2 (2000): 241-266.\n\n    and many others. The submission suggests that the \"law of total defect\" is a central contribution of this work (e.g. in lines 017f, 063f, or 073f), which is misleading because describing the approximation error of a PDE via the defect equation is a standard approach in numerical analysis. Furthermore, it's known that the approximation order of defect correction approaches is the product of the individual orders, which raises questions about the novelty of the submission's theoretical analysis in Theorem 3. In summary, while the application of defect correction techniques to training vs. inference stages of ML-based PDE solvers is new, the general construction is not. For a future version of this manuscript, I would recommend positioning this work more clearly as a ML-application of existing defect correction techniques (adjusting the claims and the related work accordingly).\n\n2. The distinction between training and inference time of ML-based PDE solvers is a bit arbitrary in the current setup: both the PINN and the Monte-Carlo solver contribute to solving the PDE, and there is no reason to delay the Monte-Carlo solver to inference time (unless I have missed something). I imagine that the separation of training and inference would become clearer for parameter estimation problems, where external data enters the loss function, but parameter estimation is not studied in the current version of the submission. I believe that a discussion of inference-time scaling with ML-based PDE solvers requires a more thorough distinction of training and inference time.\n\n3. The numerical results, while broad, are not as conclusive as they could be. On the one hand, the accuracy differences between the methods are small, but standard deviations are not reported. From the current results, it is unclear whether the improvement is statistically significant. On the other hand, I think that the accuracy improvements would be more convincing if they were normalised by budget in the sense that all three \"competitors\" (PINN/GP, MC, PINN/GP+MC) should use similar compute (either function evaluations or runtime/memory). For example, in Table 1, the baselines are slightly less accurate but much faster, and I don't know whether running the baselines for longer would make them more accurate.\n\n\n4. (Minor) Finally, and I consider this to be a smaller weakness than the three above, I wonder whether the volume of the paper (almost 40 pages) is a bit too high for an ICLR paper, and not necessarily proportional to the amount of contributions of this work. I understand it is common practice to defer proofs and extra results to the appendices, but it is currently almost impossible to read the main paper without extensive reference to the supplement. That said, while I think this point is worth mentioning, I consider it less important than the three above.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "CvA9EnQnLJ", "reviewer": "Reviewer_op81", "summary": "Inference-time scaling is a characteristic phenomenon observed in large language models in which allocating additional computation resources at inference time yields an improvement of output quality of LLMs. The paper is the first attempt to bring the idea into physics-informed emulators to solve high-dimensional PDEs. The key assumption to realize the idea is the semi-linearity of parabolic PDEs, and this makes it possible to derive a PDE, called Law of Defect, that governs the residual between the ground-truth solution and surrogate solution for an original PDE. This finding enables the application of standard solvers such as multilevel Picard solver without requiring any additional retraining and fine-turning, and the residual is also utilized to correct the error of the output of a trained surrogate model, which is the proposed inference pipeline to solve the PDEs. The proposed framework is evaluated on a variety of scenarios with various dimensions that control the difficulty of the task and compared to a couple of standard solvers and plain data-driven surrogate models.", "questions": "- What does training points $m$ mean? How does $m$ training points translate to $m$ additional Monte Carlo paths? \n- How and/or along which metric is the computational budget given for MLP in the main experiments ensured to be comparable to that taken for surrogate models and the SCaSML?\n- How long/how many training epochs/what change in hyperparameters is additionally necessary for MLP and SR models to reach the comparable performance of SCaSML reported in Table 1? While theoretical analysis on the variance and computational cost is given, I believe this observation and/or potential contrast could strengthen the significance of the contribution of the paper. \n- In what practical scenario do the authors expect the proposed framework play a crucial role?", "limitations": "", "rating": 4, "confidence": 2, "soundness": 3, "presentation": 1, "contribution": 2, "strengths": "- **Intriguing approach:** The proposed approach is disparate from the other approaches boosting the performance in the post-training phase, especially that requires additional training. The question raised could serve as a new class of tasks in the scientific discovery. \n\n- **Performance gain:** The improvement in the prediction accuracy over plain surrogate models is significant and impressive.", "weaknesses": "- **Practical implication:** Table1 of the main experiments show that surrogate models already outperform the MLP solver in most of the scenarios. It is not very clear to me what practical scenario necessitates an additional computational budget to boost the performance of the trained data-driven models, especially when the performance of the data-driven model is already superior to that of the solvers.\n\n- **Assumption for experiments:** The hyperparameter configuration for the MLP used in the experiment is not stated exhaustively and it makes it hard to access the fairness of the experiments. \n\n- **Mathematical notations for theoretical results:** I see the authors make a tremendous effort to solidify the theoretical ground of the proposed framework. However many mathematical notations are not introduced satisfactorily and I could not access the correctness of the theoretical result. Followings are some of the points which I saw were necessary to be introduced properly when I went through Appendix C. In addition, some of the symbols and notations are apparently coming from the previous sections. I suggest the section dedicated for the notations should be placed at the beginning of the appendix.\n  - Topology and metric on the extend real numbers.\n  - The definition of $\\mathcal{B}(\\mathbb{R}^{d+1})$.\n  - It is not clear that for which mathematical objects the assumptions 5 and 6 are made.\n  - The definition of the terminal function (I speculate that it is a function describing terminal conditions for PDEs)\n  - $\\breve{F}$ in Definition8. Apparently this is different from one introduced in Fact 1 in the main text. \n  - I still do not understand the assertion in Remark 8.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)