Learning Stochastic Multiscale Models
Authors: Andrew Francesco Ilersich
Prasanth B. Nair
Venue: NeurIPS 2025 poster ()
Keywords: Multiscale Models
Stochastic Differential Equations
Spatio-Temporal Data
Variational Inference
Partial Differential Equations
Primary area: machine_learning_for_sciences
TL;DR: We present a framework for learning stochastic multiscale models of spatio-temporal systems.

Abstract: The physical sciences are replete with dynamical systems that require the resolution of a wide range of length and time scales. This presents significant computational challenges since direct numerical simulation requires discretization at the finest relevant scales, leading to a high-dimensional state space. In this work, we propose an approach to learn stochastic multiscale models in the form of stochastic differential equations directly from observational data. Drawing inspiration from physics-based multiscale modeling approaches, we resolve the macroscale state on a coarse mesh while introducing a microscale latent state to explicitly model unresolved dynamics. We learn the parameters of the multiscale model using a simulator-free amortized variational inference method with a Product of Experts likelihood that enforces scale separation.  We present detailed numerical studies to demonstrate that our learned multiscale models achieve superior predictive accuracy compared to under-resolved direct numerical simulation and closure-type models at equivalent resolution, as well as reduced-order modeling approaches.

Learning Stochastic Multiscale Models
Andrew F. Ilersich
Prasanth B. Nair Institute for Aerospace Studies
Institute for Aerospace Studies
University of Toronto Toronto, ON M3H ST6
University of Toronto Toronto, ON M3H 5T6
andrew www.seihhhhhaala...c..am.mmmaa.m utoronto.ca
prasanth.naireutorooin
Abstract
The physical sciences are replete with dynamical systems that require the resolution of wide range of length and time scales. This presents significant computational challenges since direct numerical simulation requires discretization at the finest relevant scales, leading high-dimensional state space. In this work, we pro- pose an approach to learn stochastic multiscale models in the form of stochastic differential equations directly from observational data. Drawing inspiration from physics- based multiscale modeling approaches, we resolve the macroscale state on a coarse mesh while introducing microscale latent state explicitly model unresolved dynamics. We learn the parameters of the multiscale model using a simulator-fee free amortized variational inference method with a Product of Experts likelihood that enforces scale separation. We present detailed numerical studies demonstrate that our learned multiscale models achieve superior predictive ac- curacy compared to under-resolved direct numerical simulation and closure-type models at equivalent resolution, as well as reduced-order modeling approaches.
Introduction
Multiscale phenomena ubiquitous in nature, and broadly fall under two categories. The first are systems with large-scale and small-scale regimes, each of which may follow different laws but interact with each other at some interface, so simulation requires one model both regimes. Such systems are common in materials science [ and biology (2). The second, which we focus on in this work, are systems governed by partial differential equations (PDEs) that require one to resolve wide range of length and time scales to make accurate predictions. A canonical example is weather prediction. Modeling the atmosphere requires resolving phenomena from eloud microphysies (< Im) up to planetary waves (5000 km) 3 In such scenarios, direct numerical simulation (DNS) of the governing equations on grid that resolves the smallest scale of interest is computationally infeasible.
The computational intractability of fully-resolved simulations necessitates working with coarsened representations. This introduces a fundamental multiscale separation between resolved and sub- grid-scale dynamics. To "close" the system, the effect of the sub-grid-scale dynamics on the resolved dynamics must modeled. This approach, known as closure modeling 4 5 has emerged independently across multiple fields. Large-eddy simulation (LES) in fluid dynamics filters out small-scale phenomena and closes the Navier-Stoker equations with -rbdrri--scale model [6]. Similarly, in coarse- grained molecular dynamics, atoms are grouped into larger "beads" to reduce degrees of freedom Recently, machine learning techniques have been used to infer closure models for climate models coarse-grained molecular dynamics (9, and turbulence 110
We propose novel framework for learning stochastic multiscale dynamics directly from observational data. Our approach introduces a latent variable model that explicitly decomposes the system state into macroscale (C) and microscale (7) components, whose dynamics governed by learned coupled system of stochastic differential equations (SDEs). This structure captures complex inter-scale
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Given:
Filter init. cond. to obtain 1.
2.
Solve dynamies to obtain Macmscale cond.
3.
Combine macroscale microccle states obtain Final prediction
Macroscale init.
spatial
Initial condition
dc
and Microscale init. cond.
and
Microscale final cond.
Figure 1: Overview of our multiscale modeling approach.
interactions and propagates uncertainty naturally; see Figure for graphical overview. Our main contributions are:
A novel probabilistic multiscale formulation: We introduce latent variable framework for stochastic multiscale modeling, which involves learning system of coupled SDEs governing distinct, interacting macroscale and microscale latent states directly from data.
Multiscale likelihood for enforcing scale separation: We enforce explicit scale separation by formulating Product of Experts (PoE) multiscale likelihood that ensures macroscale features dominate predictions while microscale components provide only necessary corrections.
Efficient learning and empirical validation: We enable efficient training of our multiscale stochastic model using an amortized, simulator- free variational inference scheme [12] Numerical studies challenging PDE systems demonstrate that explicitly modeling microscale dynamics yields superior predictive accuracy compared to under-resolved DNS and closure-type models at equivalent resolution, as well as reduced-order modeling approaches.
Related work Reduced-order modeling (ROM) methods reduce the computational cost of sim- ulating complex systems through dimensionality reduction of the state space [14]. Traditional ROM methods employ linear projection [6]] [8]; however, the inherent limitations of linear subspaces for systems with slowly decaying Kolmogorov n-width have led to the development of nonlinear manifold methods 19 20
23 ROM methods struggle with multiscale systems (24), It is worth noting that standard ROM methods learn a low- now-dimensional, global latent representations of the dynamics. In contrast, our approach learns spatially diittillyd representation modeling framework more akin to paradigms such as LES.
21
22
motivating specialized closure models 25
26
coarse grid, positioning it as stochastic coarse-grained
There wide body of literature data-driven closure modeling: 4 5 for review. The Mori- Zwanzig (MZ) formalism provides a powerful theoretical framework for deriving coarse-grained models of dynamical systems [13 27 (28), showing that the effects of unresolved scales manifest as memory Knon-Markovian terms the evolution of the resolved variables. The MZ. framework has motivated various approaches for closure modeling (29 4 30 31 32) More recently, Boral et al. [11] proposed a neural SDE- based approach for turbulence closure modeling. This approach, motivated by ideal LES, enables probabilistic closure modeling but requires running the coarse solver in the loop during both training and inference.
Our approach fundamentally differs from data-driven closure modeling, which typically augments the macroscale dynamics with a learned correction that implicitly accounts for the sub-grid effects. In contrast, explicitly model the coupled dynamics of macroscale and microscale latent states, providing a Markovian representation capable of capturing the non-Markovian memory effects arising from coarse-graining (see Appendix BB. Moreover, approach enables prediction full resolution through a decoder, whereas closure models are inherently limited to coarse-scall predictions.
Recent work in weather forecasting has employed hierarchical graph neural networks operating on multiple mesh resolutions to learn auto-regressive models 33 approach learns continuous-time coupled SDEs governing macroscale and microscale dynamics, performing test-time simulations only on a single coarse grid and using high-resolution grids solely decoding targets.
contrast, our multiscale

---

2 Problem setting
We consider dynamical systems governed by PDEs a spatial domain 9 Rd, where d N is the spatial dimension. We denote the PDE solution by the spatio-temporal vector field u E u(2 (8, 2`; Rd=), where U denotes an appropriate function space (typically a Sobolev space) endowed with a suitable norm.
To learn multiscale model from data, work with discretized representation of u on spatio- temporal grid. We introduce a spatial projection operator, P (4(22; IR*) + Rdd., that maps a vector field to its values on spatial grid with points. Using this, define the fully-resolved state ytt) E IR" by projecting the PDE solution onto a fine spatial grid with n points: y(t) := Pn (u(:;)), ()), where the state dimension is ny nydy- This grid is assumed to be sufficiently fine faithfully represent the dynamics and thus serves as our ground truth. The training dataset y = (tt.:i M comprises noisy observations of the fully-resolve state at na time steps, i.e.,
Y/ yt,,) + E,, E ~Mee 10,o²I), for i ..2....,,
(I)
where é; represents zero-mean i...d. Gaussian noise with covariance o² I, and 4, denotes the observa- tion time-stamp. We assume access to this dataset without knowledge of the underlying governing equations.
Our objective to infer a stochastic, continuous-time time, latent variable model from the data. The stochastic framework is essential for capturing both potential inherent randomness in the underlying physics and uncertainty introduced by modeling approximations and data limitations. Next, we present our approach for learning a multiscale SDE model for the dynamics of y by decomposing it into a macroscale state on a coarse grid and a microscale state representing sub-grid-scale features.
3 Multiscale framework
We introduce our stochastic multiscale modeling framework that addresses the challenge of learning the dynamics of the fully-resolved state y(t) E Rhy from the observational data V. Our main idea is to represent the high-dimennional state y through a lower-dimensional latent state, x(t) E IR"x, that explicitly separates the macroscale state denoted by ((t) E F R and the microscale state denoted by n(t) E R"v The temporal evolution of { and is governed by learned system of coupled SDEs, capturing both deterministic dynamics and inherent stochasticity arising from scale separation and model reduction.
Our framework is composed of three core components with learnable parameters that denote by A: a probabilistic encoder that maps the observations of y to the latent state = = (C, 1)) a system of coupled SDEs governing the dynamics of z, and a probabilistic decoder for mapping back to the observation space to make predictions. The following sections detail each of these components.
Probabilistic scale separation (encoder) Inspired by classical multiscale methods like LES, decompose the fully-resolved state y into a large-scale component ア and a residual small-scale component y = y y. In contrast classical deterministic approaches, we define probabilistic encoder po (2 y) that explicitly models the uncertainty in identifying and compressing scale-separated features from the high-dimensional state. We parameterize the mean of this conditional distribution using the operators defined below that are analogous to traditional multiscale analysis:
1. Smoothing Operator Sg R -> RR serves as low-pas, filter that extracts the large-scale component, y : (y)) for example, fixed Gaussian filter or a learned convolutional network.
2. Residual Operator Sy defines the residual component, ÿ == Spt (y) = y - 7. which contains small-scale features.
3. Restriction Operator (Macroscale Encoder) E/ R2 + : maps the smoothed component,
y. to the low-dimensional macroscale latent space, Ç. with [E[[ y] = E (S)()); for example, coarse sampling or a more complex learned mapping.
While we denote ( as vector for notational simplicity, it represents field variables on a coarse spatial grid. This structure informs the parameterization of the drift and diffusion functions, which can be modeled using spatially-awaa aware architectures such as graph neural networks that are capable of handling both regular and irregular grids.
3

---

4. Microscale Encoder E/ : R= - R"" compresses the residual into the compact microscale state, 7) with B/n y) == Eg (S) (y)); for example, learned nonlinear mapping that efficiently captures complex small-scale structures.
We model the conditional distribution of given the fully-resolved state as Gaussian, i.e., (= 1 = (2 ( (y), Σ.). The mean vā(v) E|2 v - [ (S)))) Eg (S,) (v))²]² is constructed from the outputs of the macroscale and microscale encoders, and the covariance matrix Ei E I= is learned parameter representing encoding uncertainty.
Latent stochastic dynamics We model the evolution of the macroscale state, S. and the microscale state, 77, using a coupled system of Itō SDEs:
d []] [
[ SolC, (())) (0 (()) +
[(c((t)
Legn (t) d3 (t)((t Lxp (t) (.S..
(2)
with initial conditions ((4)) po(C ((oo),, n((o)) po(nl v(to)) defined by the encoder at time to- ie.. (=((()) v(t))). The learned drift functions fo : RR' X R. and go R x R² R govern the evolution of the macroscale and microscale states, respectively respectively[TThse depend on coupling functions Φ T + R" and s R.C - R. that mediate inter-scale interactions. The dispersion block matrices (Icce: Len Tng; Lon) modulate the influence of the (ng + dimmensooll Wiener process 3 with components Be and By. capturing both inherent randomness and uncertainty from unresolved scales.
These drift functions
More compactly, we can write the coupled multiscale latent SDE as
de (e(2)dt Lt(t) 2(to) ~ po(2(to) (
where (0(=)= [fe(S, (())), go(n) (())))) is the concatenated drift, and (t) E R>. is the dispersion matrix composed of the blocks Lec: Lem Lnc; Lyn 3 The coupled SDE in ( defines the prior dynamics on the latent trajectory. The solution to this SDE, initialized from the encoded state at
lo. yields the prior predictive distribution p=(2(1) y(to) for > to.
S
S
Probabilistic state reconstruction (decoder) In order to make predictions, we require a decoder, | 2). to map the latent state dynamics governed by the coupled system of SDEs in 33) to the observation space "R"", Towards this end, we introduce macroscale and microscale decoders below:
y
Eg
DD
9
1. Prolongation Operator (Macroscale Decoder) DG R + Ryy maps the macroscale state to the fine grid.
1
Figure 2: Relationships between multi- scale states and operators. Sg: smooth- ing operator; Eg: macroscale encoder; DG: macroscale decoder; Sp residual operator; E/ microscale encoder; D. microscale decoder.
Microscale Decoder D3- Th -> Tv maps the com- pressed microscale state 17 to the fine grid.
Later in Section 4.1
we formulate PoE likelihood to enforce scale separation that will establish the specific func- tional form of the probabilistic decoder (y z)) The pre- dictive distribution of the learned multiscale model that follows from the proposed multiscale likelihood is presented in Section 4.4 relationship between all operators introduced in this section in Figure 2
We summarize the
Remarks on closure models and implicit-scale models Traditional closure modeling aims to derive equations governing only the resolved (macroscale) dynamics by parameterizing the effect of unresolved (microscale) dynamics purely in terms of resolved state variables [4 3) At first glance,
2 More generally, the drift functions can be parametrized as fo(Ç, ((n). x(t)) and go (n, &(C), x()) (e.g. neural networks or physics-basee based parametrizations) with x(t) denoting time encoding. In our imple- mentation,
neural network architecture for fo that respects the spatial locality of the coarse grid variables Ç: see Appendix[]] [o details.
3
The latent SDE model implicitly assumes Markovian dynamics; we discuss the validity of this assumption
in Appendix B

---

it may appear that closure modeling is special case of our framework obtained by setting n. = 0. which reduces the latent state to the macroscale state (2 = () and reconstruction involves only the prolongation operator, with y predicted solely from Ç. However, the fundamental objective in closure modeling is to ensure that the predicted coarse state (or (C)) matches the filtered state y. nor the original observation y. In contrast, our framework, even when n, = 0, aims to predict the fully- fully-resolved state y. We therefore refer to the n. = 0 case as implicit-scale model distinguish it from conventional approaches. It implicitly represents all scales solely through and its associated decoder. Our multiscale framework (with ny > 0) explicitly models sub-grid -rid scales via 1), aiming for improved accuracy in representing the dynamics of the fully-resolved state
Multiscale variational inference
We now present simulator- free stochastic variational inference (SVI) approach to learn the parame- ters of the multiscale SDE model introduced in the previous section. Our approach involves three key ingredients: multiscale PoE likelihood that enforces scale separation, variational approximations for the latent state 2 and model parameters 0, and a reparametrized evidence lower bound (ELBO) that can be maximized without SDE solver in the optimization loop.
4.1 Multiscale likelihood
To explicitly enforce separation between macroscale and microscale states, we employ a PoE per- spective 35
to define the likelihood of observing y given the latent states Ç and v as:
po(v 2=) = poly (,,)= = Z(C, Pi () p2 - DC(())
(4)
where Z(S, f fpi (1)1))))) p2) - DS( vi)dy' denotes the normalization constant.
The term pi(y) () in 4 can be viewed as the "macroscale expert" that evaluates how well the macroscale state explains the full observation y, while p2 (<) 17) functions the "microscale expert" assessing how well the microscale state explains the residual (C) after accounting for the macroscale prediction. This PoE formulation establishes natural hierarchy where macroscale features are modeled first, with microscale features accounting only for unexplained residuals. The normalization constant reinforces scale separation by penalizing large mieroscale contributions, thereby ensuring that the microscale state does not dominate the overall prediction unless necessary.
For Gaussian experts with covariance matrices 2 and Ey, i.e., P1 (y (5) = Ny (C),), x:)) and P2(y-DI (S) m) = NVy -
Eg), the log- ligllikllhhood simplifies to (omitting constants):
1 log (.) = - (1
||((() D2) 1 I0)())IR)) +
2 log Sh + S
(5)
where r()) = y ----yy (C) is the residual after accounting for macroscale prediction, S/ (ES)
and S = (E%) are precision matrices, and A = S (S + -' We use the notation || = x" Az to denote the squared norm induced by the inner product (2. = x" Ay for
I, R" and symmetric positive-definite (SPD) matrix A E g2xx.
Expanding the terms in 5 and rearranging yields (see Appendix C for details):
log (, = - ( + (r((), (H))) sy - 11D 1 +
log S + S%
(6)
where A' = (S + -1s"" is an SPD matrix. follows from the preceding equation that maximizing the proposed multiscale likeliield produces four key effects:
1. Minimizes the (S) + 5T) -nom of the residual r(C), incentivizing the macroscale state to explain as much of the full observation as possible.
Maximizes the S -n-nneer product between the residual r(() and the decoded microscale state D (n), correlating the microscale prediction with the macroscale residual.
Minimizes the A' no-oomm of the decoded microscale state D7()), serving as an adaptive regularizer that keeps microscale contributions small unless necessary to explain the residual.
5

---

4. Maximizes the precision matrices Sj and Sy. which is equivalent to minimizing V and mn.
In summary, the structure of the proposed multiscale log- liellikelihood ensures that macroscale features dominate the overall predictions while microscale components provide only necessary corrections. This can be viewed as embodying Occam's 's razor at the microscale level and maintaining clear multiscale hierarchy. This regularization effect emerges naturally from the PoE formulation and helps prevent overfitting to high-frequency noise the observations.
4.2 Variational approximation and priors
Our multiscale model requires approximating the posterior distribution of the latent state z, whose prior dynamics are governed by (2). We specify variational approximation for the posterior of the latent state, (p l). whose sample paths follow the linear SDE:
d. (-
(t)()) (t))) L(t)d3, 2(to) (Nte E..)
where (t) E R and (t) E R are parametrized in terms of the learnable variational parameters denoted by y, while Mto IR and Enp E R"*x*= represent the mean and covariance derived from our probabilistic encoder (=(z y) for given initial condition y(to). It is worth noting that this linear SDE is used only during training. At test time, predictions are made by solving the learned nonlinear system of SDEs defined in (2).
For model parameters 9,
either adopt a Bayesian approach seek point estimates. In the former case, we specify a prior distribution p()), then construct a variational distribution gp( to approximate its posterior. Our framework can readily incorporate interpretability constraints - for example, by representing drift functions as linear combinations of basis functions with sparsity- inducing priors their coefficients.
4.3
Evidence lower bound
To train our multiscale model, we maximize the ELBO below, which we derive following 36
/
((1) = Σ
E 1=1 5..0
log pa (y xx] 2 J. 2(1).0
E (10,p(=(0),)))))) di - DKC.. (1x (0) (())
(8)
where 1 ~ (0), 24 ~ (=(= t4) and 2(t) ~ ( t). The first term is the expected log likeliioddood of the dataset )' given the latent states, incorporating the multiscale likelihood from Section 4.1 The second term is the KL divergence between the solutions of the nonlinear SDE in Equation (2) and our
linear variational SDE Z with (z()),)) The third term is the KL divergence between variational
(t)=((t) + (t) (2()))) denoting the drift
residual and C(t) (Lo (L)
distribution over model parameters and their prior. For conciseness, we present the ELBO for single time-series trajectory; the extension multiple-traeectee
datasets is discussed in Appendix D
To efficiently maximize the ELBO without requiring a forward SDE solver the training loop, we leverage the reparameterization trick for SDEs from 37 (details in Appendix E)
4.4 Computing the posterior predictive distribution
After training our multiscale model, generate probabilistic predictions by computing the posterior predictive distribution po (v((t) v(40)) given an initial observation y(lo). This distribution accounts for all sources of uncertainty captured by our model.
The prediction process involves three stages. In the first stage, we map the initial observation to a distribution over possible latent states, i.e., po(z(to) v(to)) (N(2(to) Ma, E,,), where
Ho = [ (Sp())t)))) E (S,) (y(to))) and Σ,, Σz. In the second stage. we propagate this initial distribution through our learned latent SDE, de 70(z(s))ds + (o(s)(((ss) E [to, t). yielding the non-Gaussian distribution (=(=)) | y(to)). Finally, we map the latent distribution the observation space using our PoE likelihood. This step yields a Gaussian conditional distribution for given z, whose mean depends on an adaptively weighted contribution of the decoded microscale state (see Appendix [ for detailed discussion):
po 2) = Nyy (2),), Ey), where (z(2) =D((S) + Sy (T) and E, (S$ (5)+ (9)

---

Table 1: Comparison of error statistics obtained using different models for all test cases.
Error mean + std. dev. on test set
Model type
Wave ID
KdV ID
Cylinder 2D
SWE 2D
(nic 20)
(nc = 20)
(nc 512)
(ni 64)
Coarse DNS
1.381
0.382
0.523 + 0.196
0.130 + 0.040
0.861 0.190
DMD
0.492 + 0.066
0.163 + 0.046
0.016 + 0.006
0.531 0.227
POD-SINDY
0.469 = 0.054
0.057 + 0.028
0.051 = 0.024
0.535 + 0.210
Implicit scale
0.565 = 0.116
0.324 + 0.078
0.063 = 0.005
0.426 + 0.099
Our approach
f, = 1
0.517 + 0.095
0.197 + 0.146
0.038 + 0.004
0.543 + 0.360
fy 2
0.210 = 0.095
0.077 = 0.038
0.015 0.003
0.202 0.084
fy = 3
0.105 = 0.047
0.059 = 0.030
0.011 = 0.001
0.190 + 0.075
fy 4
0.089 = 0.035
0.062 = 0.029
0.011 0.002
0.163 0.056
fy = 5
0.079 + 0.031
0.042 + 0.017
0.009 + 0.001
0.152 + 0.047
The posterior predictive distribution then obtained by marginalization:
(v(((t) [((())) = ((((() ( (,(((()), E.)) pa po(2(() (=(()((((t))))))))((())))))))))) ((t))) dz(().
(10)
Since this integral is analytically intractable, we use Monte Carlo sampling leading to the Gaus- sian mixture approximation: po (v(( ((to)) s (1/N) E N(y((t) (z)()()),, Σ,,), where (t)
denotes the state at time f of the ith sample trajectory, generated by integrating the learned SDE (3) starting from initial condition 2t)) N(r,, E.,). We provide further details the prediction
procedure in Appendix F
Results and discussion
We evaluate our multiscale framework on four systems exhibiting slowly-decaying Kolmogorov n--idth and energy cascade between scales: (1) the 1D advecting wave with varying initial conditions, (2) the ID Korteweg-de Vries (KdV) equation with varying initial conditions, (3) a 2D Von Kármán vortex street, and (4) radial dam break modeled with the shallow water equations (SWE) with varying initial conditions Setup/training details including an extra test case are provided in Appendix
We adopt a hierarchical training strategy for our multiscale models, varying R. E f0,...,5. 5}. The implicit scale model was trained first. Its learned parameters were then used initialize the training of the multiscale model with n, = 1. This process was repeated, with each trained model initializing the next in the sequence, incrementing n, by one. For implicit scale and multiscale models, offline training cost therefore scales linearly with ne;; for the four test cases, the training time in minutes is (1) 46 + 51 (2) + 50n,, (3) + 64n,, and (4) 67 79m.
We compare our multiscale models against four baselines: a coarse DNS, a DMD linear system, reduced-order model formed over a POD basis with latent dynamics learned by SINDy, and an implicit scale model (n,, 0). The number of gridpoints in the coarse DNS, the rank of the DMD model, and the latent state dimension in the POD-SINDY model are chosen be equal to the total latent state dimension of our best-performing multiscale model (i.e., ni + 5). Further details are provided in Appendix
Performance is evaluated using the normalized error metric e(t,) = 1ly, - l/tlll I.. where ys is the observation time and (t,,)) is the mean prediction. All experiments were conducted system with 24 CPU cores, 128GB RAM, and an NVIDIA RTXX090 GPU. Results across all test cases are summarized in Table 1 showing the mean and standard deviation of the error metric each test set. Additional visualizations of the results are provided Appendix[]
Code for reproducing our results is available: https http://ttw...c.aalllaee//tttaaat
7

---

+
-wm - + mmmmom
(A)
(B)
(C)
(D)
Figure 3: Advecting wave (top row) and KdV equation (bottom row): scale separation with ny
shown at & = 0.2 (wave) and 1 (KdV). The columns are (A) observation from dataset, (B)
reconstructed full state, (C) prolonged macroscale state, and (D) decoded microscale state.
One-dimensional advecting wave We test on the one-dimensional advection equation with periodic boundary conditions, i.e,, 14 Clug, where 1, E [0, 1], and [0, The initial condition exp(-(x//w) / where w ~ (400001, 0.02]. We generate training trajectories on a fine spatial grid with spacing 0.001 and time-step 0.001 (ng 1001 and my 1000), which are then corrupted by zero-mean Gaussian noise with standard deviation 0.001. The test set covers / E [0, 0.2]. The training dataset comprises 20 trajectories, with each for validation and testing. For multiscale models, we use a coarse macroscale mesh with only ni = 20 points.
On this classic problem, where linear projection methods struggle, the ability to explicitly model sub-grid features proves crucial. As shown in Table [ our multiscale model with n2 5 achieves an 83% error reduction compared to the best baseline (POD-SINDy) 0.469 to 0.079. The multiscale models outperform all baselines on this Figure 3 illustrates the learned scale separation, showing how the microscale component r) captures sub- -rid features that are lost in the coarse macroscale representation.
with mean error dropping from
problem.
One-dimensional Korteweg-de Vries equation The second test problem is the one- oie-dimenso Korteweg-de de Vries (KdV) equation with periodic boundary conditions: 14 uu, V/llaa 0, where x E [o, 10] and = 0.02. The initial condition a cos(TT) exp(-(x 7 where spatial grid with spacing
(2, 0.01) and 8 (11,0.01) We generate training trajectories 0.01 and time-step 0.001 (n2, = 1001, Ny = 1000), which are then corrupted by Gaussian noise with standard deviation 0.01. The training dataset comprises 10 trajectories, with 5 each for validation and testing. Our multiscale models use coarse macroscale mesh with only nę 20 points.
Table shows that our multiscale models dramatically outperform all baselines. The best multiscale model (n, = 5) achieves 87% error reduction compared to the implicit-scale baseline, with mean error dropping from 0.324 to 0.042. Among reduced-order methods, POD-SINDy performs best but still has 36% higher error than our approach.
Figure demonstrates scale separation in the KdV system. The coarse macroscale mesh cannot resolve the wave structure, making the microscale contribution critical. Interestingly, the macroscale state contains apparent small-scale features beyond what simple smoothing would produce. This arises from our learned convolution kernel, which differs from a standard Gaussian filter as it attempts to approximate babb-grd-sccale features (see Appendix J for details).
Two-dimensional flow over a cylinder Our next test case examines a 2D Von Kármán vortex street at Re = 160 from Günther et al. [38]. This dataset consists of nondimensionalized velocity field
(d,, = 2) developing from zero initial condition. The spatial domain is truncated to [-0.5, 3.5] x [-0.5, 0.5] and sampled on a fine 320 80 grid, resulting in a high-dimensional state (n, = 2 - 80 51200). The cylinder is modeled as region of zero velocity. The time domain spans [0, 15] with 1501 time steps, which we partition into training j0, 13), validation (13, 14), and test (14, 15] intervals. Our multiscale models use coarse 32 8 macroscale grid (nę 2 32 8 = 512).
Table demonstrates that our multiscale models significantly outperform the implicit-scabe scale baseline, with the best model (n., = 5) reducing error by 86%. Figure [ reveals how the implicit-scale model fails to resolve the intricate vortex structures in the cylinder wake, whereas the multiscale model captures these flow patterns accurately. Figure 5 provides further insight into the scale separation
8

---

Observation
Coarse DNS
DMD
POD-SINDY
Implicit scale
Ours (nt 5)
Figure 4: Cylinder flow: prediction comparison at = 15 the test interval. The z-component of velocity is shown here, zoomed in the cylinder.
2
(A)
Figure 5: Cylinder flow: scale separation with Ryy 5, shown at 15, zoomed in the cylinder. The columns are (A) observation from dataset, (B) reconstructed full state, (C) prolonged macroscale state, and (D) decoded microscale state.
mechanism, showing how the microscale state encodes the unresolved vortex structures. Notably, as vortices grow larger toward the right the domain, the microscale contribution diminishes, demonstrating model's ability adaptively allocate representational capacity across scales based on the local physics. The periodic dynamics and relatively fast-decaying Kolmogorov n-width in this problem allows DMD to achieve competitive performance (error 0.016), though our multiscale model still achieves 44% lower error.
Two-dimensiooal shallow water equations Our final test case models a radial dam break using the shallow water equations, from the PDEBench dataset [39]. The spatial domain is [-2.5, 2.5) discretized using a 128 x 128 grid, resulting in a high-dimensional state (n, 16384). The fluid elevation field evolves from initial condition of height within radius ~ (0.33, 0.7] of the origin and height elsewhere. The time domain spans [0, 1] with 101 time steps. We corrupt observations with Gaussian noise (standard deviation 0.01) and partition 1000 trajectories into 900 training. 50 validation, and test samples. Our multiscale models use a coarse 8 x macroscale grid (ng = 64).
Table 1 shows our multiscale model (n,g 5) achieves 64% error reduction compared to the best baseline (implicit scale), with mean error dropping from 0.426 152. The coarse DNS exhibits particularly high error (0.861), while reduced-order methods (DMD and POD-SINDy) achieve inter- mediate performance with error of around 0.53. This problem features slowly-decaying Kolmogorov n-width due to significant advection and the diverse ensemble of initial conditions, explaining the challenges faced by linear manifold methods. Figures and ( illustrate how the microscale representation captures the complex flow patterns that enable this performance improvement.
Observation
Coarse DNS
DMD
POD-SINDy
Implicit scale
Ours (n.,
Figure 6: Shallow water: prediction comparison on example trajectory from test set at f 1.
Computational cost studies Figure 8 shows the error-cost tradeoffs for online prediction across all test cases. We use the product of number of drift function evaluations (NFE) and state dimension as a measure computational cost. Our multiscale models achieve favorable error-cost tradeoffs compared to the baselines. Interestingly, increasing n. from to 5 yields significant accuracy gains without proportionally increasing inference cost. This highlights the computational efficiency of our hierarchical structure, which adaptively engages the microscale model only needed.

---

2
+
(A)
(B)
(C)
(D)
Figure 7: Shallow water: scale separation with n, 5, shown at 1. The columns are (A)
observation from dataset, (B) reconstructed full state, (C) prolonged macroscale state, and (D)
decoded microscale state.
10k
10-1
10-1
10-1
Coarse DNS
Belartow
10-1
10-1
10-1
2x 10-1
POD SINDY
Imp. scale
10A NEE state dim.
107
10t
10F
10
Our
NFE state din.
NFE state dim.
NFE x state dim.
Wave
KdV ID
Cylinder 2D
SWE 2D
Figure 8: Mean relative error for each model type vs. NFE x latent state dimension (just state dimension for DNS), as proxy for online inference cost. The cost for single test trajectory.
Concluding remarks
We have introduced a framework for learning stochastic multiscale models directly from data, where macroscalc and microscale states evolve as distinct but interacting dynamical systems. This framing represents fundamental departure from physics-basee based and data-driven closure modeling paradigms. Instead of parameterizing the effects of unresolved scales as corrective term dependent only on the coarse state, introduce distinct, dynamical latent state for the microscale itself. This allows our model to capture the complex, state-dependent interplay between scales that characterizes many challenging physical systems.
The performance improvements are substantial: across test cases, our multiscale models achieved approximately order of magnitude reduction in error compared to implicit-scale baselines. Our approach also consistently outperformed reduced-order modeling methods (DMD and POD-SINDY), with the performance gap most pronounced on problems with slowly-decaying Kolmogorov n-width. Crucially, unlike global ROM approaches that sacrifice spatial structure, our framework maintains a physically interpretable representation the learned microscale corrections.
coarse grid while achieving superior accuracy through
While our framework shows strong performance on the test problems presented, scaling to more complex systems presents opportunities for methodological advances. While the current global compression of effective, it may not be optimal for problems like 3D turbulence where sub-grid dynamics are dominated by local structures. Future work could explore spatially-aware represen- tations for the microscale state. Furthermore, extending the framework to chaotic systems will require careful consideration of long-term stability in the learned dynamics. For systems with distinct scale hierarchies, could envision incorporating multiple nested coarse grids, creating sequence of latent states corresponding to different scales. While this hierarchical generalization presents non-trivial challenges, it represents a compelling direction for future work.
Looking forward, we believe the potential applications extend to numerous fields where multiscale phenomena present fundamental modeling challenges, including climate modeling [40] astrophysical simulations [II], biological pattern formation 42]. and neural dynamics 44 Ultimately, our work demonstrates path toward building models that combine the flexibility of data-driven methods with physically-gee grounded gryouudlly-gooundedd structure of multiscale formalisms. This synthesis of machine learning and physics is particularly promising for complex systems where first-principles closures remain clusive.
10

---

Acknowledgments
This research is supported by an NSERC Discovery Grant, Queen Elizabeth Graduate Scholarship in Science and Technology, and grant from the Data Sciences Institute at the University of Toronto.
References
[1] Weinan E. Principles of multiscale modeling. Cambridge University Press, 2011.
[2] Grace C. Y. Peng, Mark Alber, Adrian Buganza Tepole, William R. Cannon, Suvranu De. Savador Dura-Bernal, Krishna Garikipati, George Karniadakis, William W. Lytton, Paris Perdikaris, Linda Petzold, and Ellen Kuhl. Multiscale modeling meets machine learning: What can we learn? Archives of Computational Methods in Engineering, 28(3):0017-1037, May 2021.
[3] Hannah Christensen and Laure Zanna. Parametrization in Weather and Climate Models. Oxford University Press, 12 2022.
[4] Shady E. Ahmed, Suraj Pawar, Omer San, Adil Rasheed, Traian Iliescu, and Bernd R. Noack. On closures for reduced order models-a spectrum of first-principle machine-learned avenues. Physics of Fluids, 33(9):091 301, 09 2021.
[5] Benjamin Sanderse, Panos Stinis, Romit Maulik, and Shady E. Ahmed. Scientific machine learning for closure models multiscale problems: review. Foundations of Data Science, 7(1):298-337, 2025.
[6] Zhiyin Yang. Large-eddy simulation: Past, present and the future. Chinese Journal of Aeronau- tics, 28(1): 11-24, 2015.
[7] Soumil Y. Joshi and Sanket A. Deshmukh. A review of advancements in coarse-grained molecular dynamics simulations. Molecular Simulation, 47(10-11):786-866) 2021.
[8] Karl Otness, Laure Zanna, and Joan Bruna. Data-driven multiscale modeling of subgrid parameterizations in climate models. ICLR 2023 Workshop: Tackling Climate Change with Machine Learning, 2023.
[9] Linfeng Zhang, Jiequn Han, Han Wang, Roberto Car, and Weinan E. DeePCG: Construct- ing coarse-grained models via deep neural networks. The Journal of Chemical Physics, 149(3):03410,, 07 2018.
[10] Yunpeng Wang, Zelong Yuan, Chenyue Xie, and Jianchun Wang. A dynamic spatial gradient model for the subgrid closure in large-eddy simulation of turbulence. Physics of Fluids, 33(7):075119, 07 2021.
[11] Anudhyan Boral, Zhong Yi Wan, Leonardo Zepeda-Núneed James Lottes, Qing Wang, Yi Fan Chen, John Roberts Anderson, and Fei Sha. Neural ideal large eddy simulation: Modeling turbulence with neural stochastic differential equations. Advances in Neural Information Processing Systems (NeurIPS) 36, 2023.
[12) Kevin Course and Prasanth Nair. Amortized reparametrization: Efficient and scalable variational inference for latent SDEs. S. Levine, editors, Advances in Neural Information Processing Systems, volume 36, pages 78296-78318. Curran Associates, Inc., 2023.
Oh, T. Naumann, Globerson, K. Saenko, M. Hardt, and
[13] Robert Zwanzig. Nonequilibrium statistical mechanics. Oxford university press, 2001.
[14] Peter Benner, Stefano
Alfio Quarteroni, Gianluigi Rozza, Wil Schilders, and Luís Miguel Silveira. Model order reduction: Volumes 1-3. De Gruyter, 2022.
[15] Christophe Audouze, Florian De Vuyst, and Prasanth B. Nair. Reduced-order modeling of parameterized PDEs using ceaee-ppac-pppaameetee principal component analysis. International Journal for Numerical Methods in Engineering, 10(8))1025-1057, 2009.

---

[16] Christophe Audouze, Florian De Vuyst, and Prasanth B. Nair. Nonintrusive reduced-org order modeling parametrized time-dependent partial differential equations. Numerical Methods for Partial Differential Equations, 29(5):1587-1625 2013.
[17] Benjamin Peherstorfer and Karen Willcox. Data-driven operator inference for nonintrusive projection-based based model reduction. Computer Methods in Applied Mechanics and Engineering, 306: 16:1966 19-211, 2016.
[18] Cheng Huang and Karthik Duraisamy. Predictive reduced order modeling of chaotic multi-scale problems using adaptively sampled projections. Journal of Computational Physics, 491:1 12356, 2023.
[19] Wei W. Xing, Akeel A. Shah, and Prasanth B. Nair. Reduced dimensional Gaussian process emulators of parametrized partial differential equations based on isomap. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 471 471(20140697), 2015.
[20] Wei W. Xing, V. Triantafyllidis, Akeel A. Shah, Prasanth B. Nair, and Nicholas J. Zabaras. Manifold learning for the emulation of spatial fields from computational models. Journal of Compulational Physics, 326:666-690, 2016.
[21] Romit Maulik, Bethany Lusch, and Prasanna Balaprakash. Reduced-order modeling of advection-dominated systems with recurrent neural networks and convolutional autoencoders. Physics of Fluids, (33)0037 2021.
[22] Benjamin Peherstorfer. Breaking the Kolmogorov barrier with nonlinear model reduction. Notices of the American Mathematical Society, 69(5)725--73, 2022.
[23] Kookjin Lee and Kevin T. Carlberg. Model reduction of dynamical systems on nonlinear mani- folds using deep convolutional autoencoders. Journal of Computational Physics, 404:108973, 2020.
[24] Tim Keil and Mario Ohlberger. Model reduction for large scale systems. In Ivan Lirkov and Svetozar Margenov, editors, Large-Scale Scientific Computing, pages 16-28. Springer International Publishing, 2022.
[25] Jared L. Callaham, Jean-Christophe Loiseau, and Steven L. Brunton. Multiscale model reduction for incompressible flows. Journal of Fluid Mechanics, 973:A3, 2023.
[26] Changhong Mou, Birgul Koe, Omer San, Leo G. Rebholz, and Traian Iliescu. Data-driven variational multiscale reduced order models. Computer Methods in Applied Mechanics and Engineering, 373:113470, 2021.
[27] Alexandre J. Chorin, Ole H. Hald, and Raz Kupferman. Optimal prediction and the Mori- Zwanzig representation of irreversible processes. Proceedings of the National Academy of Sciences of the United States of America, 97(7)2268-22973 2000.
[28] Alexandre Chorin and Panos Stinis. Problem reduction, renormalization, and memory. Commu- nications in Applied Mathematics and Computational Science, 1(1):1-27, 2006.
[29] Eric J. Parish and Karthik Duraisamy, Non-Markovian closure models for large eddy simulations using the Mori-Zwanzig formalism. Phys. Rev. Fluids, 2:014604, Jan 2017.
[30] Anirban Gupta and Pierre F. J. Lermusiaux. Generalized neural closure models with inter- pretability, Scientific Reports, 13:10634, 2023.
[31] Kevin K. Lin and Fei Lu. Data-driven model reduction, Wiener projections, and the Koopman- Mori-Zwanzig formalism. Journal of Computational Physics, 424:109864, 2021.
[32] Priyam Gupta, Peter J. Schmid, Denis Sipp, Taraneh Sayadi, and Georgios Rigas. Mori- latent space Koopman closure for nonlinear autoencoder. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 481, 2025.
12

---

[33] Remi Lam, Alvaro Gonzalzz Matthew Willson, Peter Wirnsberger, Meire Fortunato, Ferran Alet, Suman Ravuri, Timo Ewalds. Zach Raton-Rooen, Weihua Hu, Alexander Merose, Stephan Hoyer, George Holland, Oriol Vinyals, Jacklynn Stott, Alexander Pritzel, Shakir Mohamed, and Peter Battaglia. Learning skillful medium-range global weather forecasting. Science, 382(6677): 14106667):11116610 14111442,, 2023.
[34] Joel Oskarsson, Tomas Landelius, Mare Peter Deisenroth, and Fredrik Lindsten. Probabilistic weather forecasting with hierarchical graph neural networks. In A. Globerson, L. Mackey, D. Belgrave, A. Fan, U. Paquet, Tomezak, and C. Zhang, editors, Advances in Neural Information Processing Systems, volume 37, pages 41577-41644 Curran Associates, Inc., 2024.
[35] Geoffrey E. Hinton. Products of experts. 9th International Conference on Artificial Neural Networks: ICANN '99, 1999:1-6, 1999.
[36] Cédric Archambeau, Manfred Opper, Yuan Shen, Dan Cornford, and John Shawe-taylor. Varia- tional inference for diffusion processes. In J. Platt, D. Koller, Y. Singer, and S. Roweis, editors, Advances in Neural Information Processing Systems, volume 20. Curran Associates, Inc., 2007.
[37] Kevin Course and Prasanth B. Nair. State estimation of a physical system with unknown governing equations. Nature, 622:261-267, Oct 2023.
[38] Tobias Günther, Markus Gross, and Holger Theisel. Generic objective vortices for flow visual- ization. ACM Transactions on Graphics (Proc. SIGGRAPH), 36(4):141:1-141:111 2017.
[39] Makoto Takamoto, Timothy Praditia, Raphael I Leiteritz, Daniel MacKinlay, Francesco Alesiani, Dirk Pflüeer, and Mathias Niepert. PDEBENCH: An extensive benchmark for scientific machine learning. Proceedings of the 36th Conference on Neural Information Processing Systems (NeurIPS 2022) Track on Datasets and Benchmarks, 35:1596-1611, 2022.
[40] Sungduk Yu, Zeyuan Hu, Akshay Subramaniam, Walter Hannah, Liran Peng, Jerry Lin, Mo- hamed Aziz Bhouri, Ritwik Gupta, Björn Lütjens, Justus C. Will, Gunnar Behrens, Julius J. M. Busecke, Nora Loose, Charles Stern, Tom Beucler, Bryce Harrop, Helge Heuer, Ben- jamin R Hillman, Andrea Jenney, Nana Liu, Alistair White, Tian Zheng, Zhiming Kuang, Fiaz Ahmed, Elizabeth Barnes, Noah D. Brenowitz, Christopher Bretherton, Veronika Eyring, Savannah Ferretti, Nicholas Lutsko, Pierre Gentine, Stephan Mandt, David Neelin, Rose Yu, Laure Zanna, Nathan M. Urban, Janni Yuval, Ryan Abernathey. Pierre Baldi, Wayne Chuang, Yu Huang, Fernando Iglesias-Suarez Sanket Jantre, Po-Lun Ma, Sara Shamekh, Guang Zhang, and Michael Pritchard. ClimSim-Oulines A large multi-scale dataset and framework for hybrid physics-ml climate emulation. Journal of Machine Learning Research, 26(142):1-85 2025.
[41] Arjen van Elteren, Inti Pelupessy, and Simon Portegies 7wart. Multi-scale and multi-domain computational astrophysics. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 372(2021):20133 200033335, 2014.
[42] Mohit Dalwadi and Philip Pearce. Universal dynamics of biological pattern formation in spatio-temporal morphogen variations. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 479(2271): 2023.
[43] Antonio Jimenez-Marin, Thai Diez, Asier Erramuzpe, Sebastiano Stramaglia, Paolo Bonifazi, and Jesus M. Cortes. Open datasets and code for multi-scale relations on structure, function and neuro-genetics in the human brain. Scientific Data, 1 11(1):256, Feb 2024.
[44] ling Li and Panos Stinis. Mori-Zwanzig reduced models for uncertainty quantification. Journal of Computational Dynamics, 6(1):39-68, 2019.
[45] Simo Särkkä and Amo Solin. Applied stochastic differrntiaal equations, volume 10. Cambridge University Press, 2019.
[46] Igor A. Baratta, Joseph P. Dean, Jørgen S. Dokken, Michal Habera, Jack S. Hale, Chris N. Richardson, Marie E. Rognes, Matthew W. Scroggs, Nathan Sime, and Garth N. Wells. DOLFINx: The next generation FEniCS problem solving environment, December 2023.
13

---

[47] Brian M. de Silva, Kathleen Champion, Markus Quade, Jean- Christophe Loiseau, J. Nathan Kutz, and Steven Brunton. PySINDy: A Python package for the sparse identification of nonlinear dynamics from data. Journal of Open Source Software, 5(49), 2020.
[48] Sara M. Ichinaga, Francesco Andreuzzi, Nicola Demo, Marco Tezzele, Karl Lapo, Gianluigi Rozza, Steven L. Brunton, and J. Nathan Kutz. PyDMD: python package for robust dynamic mode decomposition. Journal of Machine Learning Research, 25(417):1-9, 2024.
[49] Bengt Fornberg. Generation of finite difference formulas arbitrarily spaced grids. Mathemat- ics of Computation, 51:699-706, 1988.
[50] Stéphane Popinet. Free computational fluid dynamics. ClusterWorld, 2(6), 2004.
[51] Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. 3rd International Conference on Learning Representations, 2015.
14

---

A Macroscale drift parametrization
We parametrize the macroscale drift o(n)) to respect the underlying spatial structure of the coarse grid. This is achieved by defining the drift at each grid point as the output of a single, shared learnable function, fg. which takes as imput the state values from a local neighborhood or "stencil" around that point. This approach ensures that the learned dynamics are cunslation-qqivarrinn meaning the same physical rule is applied everywhere on the grid.
For a one-dimensional system (with du = 1), the drift component at grid point is defined as [fo(S, ((n)),, Fo (Cz-q..... S5tg: ((1)) for j = 1=...,, nc:
(11)
where q E N defines the stencil half-width, and fo a neural network with weights that are shared across all positions j. Boundary conditions are handled as appropriate for the specific test case. This structure is learnable, nonlinear generalization of a finite difference scheme. While this stencil-based MLP is highly effective for regular grids, the principle of a shared, local operator can be extended to irregular grids using architectures like graph neural networks.
It is worth noting that if the governing equations are known, fo can be derived using a spatial discretization method, such as finite difference (FD) scheme, on coarse macroscale spatial mesh We do not adopt this approach for two primary reasons.
Firstly, even when an FD scheme is derived directly for the coarse grid spacing (Ar) coarse. from the original PDE, its accuracy is limited by truncation errors inherent to coarse discretizations. Standard low-order FD stencils may provide poor numerical approximation of the original PDE operator on such a grid. Secondly, and more crucially, the effective governing dynamics on a coarse grid often differ significantly in their functional form from the original PDE. The process of coarse-rriniin means that unresolved subgrid-scale physics exert an influence on the resolved macroscale dynamics. This influence can manifest as complex, state-dependent modifications to the original PDE terms or even as entirely new effective terms (often referred to as closure terms or subgrid-scale parametriza- tions). simple FD discretization of the original PDE, regardless of the care taken in choosing stencil coefficients for (Az) coarse. is generally incapable of representing these emergent subgrid-scan scale effects. In contrast, by learning the drift function directly from observed trajectories, our approach allows the model to capture these effective coarse-grid dynamics, including any implicit subgrid-scale contributions, without being restricted to the structure of the original PDE. This flexibility is key to accurately modeling multiscale systems on computationally tractable coarse grid.
Finite Difference
Learned
0.05
0.20
0.00
0.00
i idG-anal RRG
-0.05
0.25
0.1
0.0
-0.1
To demonstrate that our representation enables interpretable models of macroscale dynamics to be learned, we perform a numerical study on the KdV test problem from Section 5 of the main paper. Figure D presents a sensitiv- ity study of the localized drift fo (from Eq. (TT) with q = 2) for implicit-scale model. This model, corresponding to our full framework with mieroscale dimension ny 0 (see "Re- marks on closure models and implicit-scale mod- els" in Section 4 of the main paper), was trained on the KdV dataset. For each subplot in Fig- ure 9 corresponding to an input stencil loca- ... 2} (since q = 2 here), the macroscale state component G is varied over [-1, 1], while the other Ç components within the stencil are held at zero. The resulting learned drift is compared to the drift term derived from standard FD approximation of the KdV equation (details of the FD scheme are in Appendix H.2)
0.5
25
0.0
0.0
-25
-0.5
[ (
0.5
05
0.0
-0.5
0.1
-0.1
0.5
0.0
0.5
1.0
Porturbation in G-
tion l E (-2....
Figure 9: Sensitivity of the localized macroscale drift fo to perturbations of its stencil inputs Ci Gi+2. The drift To is from an implicit- scale model trained on the KdV dataset.
15

---

Since the learned macroscale state ( may have different magnitudes or units than the original physical state due to the learned encoder, the outputs of both the learned drift and the FD-derived drift are rescaled to a common range in Figure to enable direct visual comparison.
It can be seen from Figure 9 that the drift derived using FD scheme, not being optimized for the coarse grid, indeed differs from the learned drift function. The unresolved small-scale features introduce complex effects into the coarse-grid macroscale dynamics (such as apparent nonlinearities or memory effects) which a standard FD scheme cannot capture. These results motivate the present approach where the drift functions are learned directly from data.
B Validity of Markovian latent dynamics
Our stochastic multiscale modeling framework models the latent dynamics 2 = (C,7) with an Itô SDE, which assumes the process is Markovian. However, when coarse grassingainngg complex physical systems, the dynamics of the resolved variables often exhibit memory effects, rendering them non-Markovian. As mentioned in the Introduction, the Mori-Zwanzig formalism [3] provides a theoretical basis for understanding these memory effects in the context of closure modeling.
The Mori- Moni-Zwanzig formalism shows that the exact evolution of of resolved variables y (analogous to our macroscale state () interacting with unresolved variables ÿ can be described by a generalized Langevin equation:
d
{ K K(s)r(t- - s)ds
dt ((1)
MM(()
F(1), Noise
(12)
Markow
Memory
where My is a Markovian term that captures the resolved dynamics, F is noise term arising from the unresolved dynamics associated with y. and the convolution integral involving the memory kernel (s)) captures the non-Markovian feedback from the history of y(() due to its coupling with z(().
To obtain a dimeessimmiiinal Markovian representation suitable for learning SDE models, the memory integral in ([2) must be approximated. One common approach (c.g.. (44)) is to introduce hierarchy of auxiliary variables Wo, W1,- ..... Wm to represent the memory term and its derivatives. For instance, defining (n(()= = f K(s)T(I ds, the generalized Langevin equation can be recast as part of an extended, Markovian system:
dy(t)
MD(!) wo(() + F(!)
(13)
d/.
dw)(t) d/.
K(O)(() + 1 (((s))((-s)dss
(14)
where K' denotes the time derivative of the kernel. This procedure can be iterated by defining w for the integral (14) and accounting for the dynamics of W1, which leads to
d
dt (t) = Mp(()= wo(() + F(t)
d
di (I) ((O)T(( + (I)
(15)
di wi (1) = K'OO(()+ + I (')(((( - s)ds,
This iterative process generates a hicrarchy of auxiliary variables wo, wi...., W. that collectively represent the memory effects. The augmented state (7, Wo, Wi...., Wm) becomes approximately Markovian if this hierarchy is truncated at a level m where the remaining integral term is negligible. This truncation an exact representation if the kernel a polynomial of degree p - 1, in which case the hierarchy terminates after p auxiliary variables.
In our multiscale modeling framework, we not explicitly construct these Mori-Zwanzig auxiliary variables. Instead, we posit that our learned latent SDE system for (ç, 1) addresses the memory challenge two ways. Firstly, the explicit microscale state 7 aims to capture some of the fast- evolving components that would otherwise contribute to the memory kernel K ((t) in model for
16

---

alone. By conditioning the evolution of Ç on (and vice-versa) within a coupled Markovian SDE, the effective memory the Ç dynamics may be significantly reduced. Secondly, the flexibility the universal approximators used for the learned drift and diffusion functions of our SDE, operating on the combined state 2, are sufficienlly flexible to implicitly represent any residual short-term memory effects within the chosen finie-dimennionnl latent space.
In summary, while the Markovian assumption for latent dynamics is an approximation for most physical systems, it can be justified if the learned latent space z is structured and expressive enough to capture the dominant interactions and short-term memory. The success of our model in practice (see Section 5 of the main paper) lends empirical support to this assumption for the systems studied.
C Multiscale likelihood derivation
We summarize below the algebraic manipulations used to derive the multiscale log- lig-likelhhood in Section the case of product of Gaussian experts. Adopting the PoE perspective [35]. we define the following multiscale likelihood:
po ((:))))) Z(S,m) (1)( P2 (D()))(())))))
(16)
For the case of Gaussian experts, we have pı (c) NNw Σg) and pzlv DV(() In)
Ni( - (C)) (n), Eg),i.e.
( (y - (F) DG(1))
pi(v = Z1 exp
(17)
P2(U - (2(() =
15(0 (1v DS(()) - (2(n)) (E)) ((v D)(())- (C(G))))(()))
D2(n))
Z2 exp
(18)
where Z1 (2
are the normalization constants of the two 1/// and Ze (2 n222)//2222)) |3 1/2
1/2
Gaussian PDFs.
We define the precision matrices as S = (2: and Sg = (>:) For compactness in the subsequent derivation of the product, we introduce the temporary notational shorthands: fai DD (C) and P2 = ()) - (n). Using these definitions, the product of the two Gaussian PDFs can be written
()(y(() p2(v DD(C)) (n)
(-1 (y m1) 55 -m1 -/ -22) S7 122)
exp
K
e
(y Syy - (n))
exp
(19)
K
where Sy = S/ + S By K (2m) ""
S.
(S + 122,), C = 1 Sylt - SUT + 2bz S2H2), and
The normalization constant of the product of Gaussians is given by
eC Sexp
c / exp (-1 -- n.,) S, (n)) dy
Z(5, = K
K (27))o/?115/11
(20)
which yields the following expression for the normalized PoE likelihood
po(yl (.m) = ",,n)) 1| exp (-) (y my)) "'y 1n)
(21)
17

---

= -2 + (n),, and using the resulting expression for My = (S + 12), the log-likelihood can be written as
Recalling our definitions
2 (v ny) " Sy(u My) - ny (2 + log [
log po (. =
2 (1r(()) 2 ||((() Di + 1 211D((n)1
2
ny log log(2 + log |S + S
1
(22)
2
2
where A = +S and r({) = y- (C( denotes the residual between the full observation and the macroscale model prediction.
To derive the expression for the log-likelihook in 6 which provides valuable insights into the regularization terms arising from the PoE likelihood, we first simplify the of the second and third terms (22) as follows:
2 [I((() - Di + |1(()) =
2 r()) ((7))))))))
DK(m) sTD + 1 2 D3(0) (D)())
2 = (r())))(((())))))))))))))
1
Pi(n)" ())
(23)
2
Noting that
= == 8/ - A = - -) 56 5% [ (S) + -/)-1s S
= 53 (S8 S (S) + S (S) +
(S) + + Sg
(24)
is an SPD matrix, the last term in 23 can written as - [D (n) l which is pon-iise by definition. Substituting (3) into (2) and combining the first two terms yields the expression for the log-likelihood in (6 that we reproduce here for clarity:
log (. |r(S))] ($$+$)) + + (r(s), ((())) 11
Fy log((.)) + log S + S
2
The term (m) implements an automatic scale separation mechanism by strongly pe- nalizing microscale components in directions already well-represented by the macroscale model and allowing the microscale state more freedom directions poorly captured by the macroscale model. This creates an "adaptive regularization" effect, wherein the regularization strength for each component of the microscale state is automatically calibrated based on the relative confidence in macroscale versus microscale predictions for that component, as encoded in the precision matrices St and Sy This is performed automatically based on the learned precision matrices without requiring manual specification of scale boundaries.
D Multiple trajectory training
The SVI approach presented in Section 4 is formulated for dataset with single time series. In this section, we consider the case of multiple trajectories with potentially different initial conditions. We assume that the multiscale model parameters multiple-trajector
shared across all time series. We denote the dataset as y = } where Vi. is the k-th trajectory (ta. 3ar, and Ner is the number of trajectories. The i-th time step of the k-th trajectory is denoted l,k, and the corresponding observation is yiu- The multiple-traeectrry ELBO is then given by simply including an extra summation in
Ner
Th
((p)
[log pe (yae 2x)
Σ
E
(1 (E)()) dt
#
(1),0
(25)
km
DKL (0) p()))
18

---

where Tj is the final time of the k-th trajectory, 6 ~ (+),), 2ihe (2(2 ttx)), and (t) ~ 9p(2x It).
The log- likeliklioood term evaluated each observation time tu of each trajectory is and summed over all observations. The drift residual term is evaluated over the entire time span of each trajectory.
E Reparametrized, amortized ELBO
In this section, outline how multiscale SVI scheme, presented in Section 4 of the main paper, leverages the reparametrization trick and amortized SVI scheme proposed by Course and for efficient training.
Nair 37 12
key challenge evaluating the ELBO defined in (8) of the main paper that we require a SDE solver to generate sample trajectories of the linear variational SDE 7 The reparametrization trick, proposed in [37], allows us to circumvent the need for an SDE solver in the training loop. The core idea is that the marginal distribution of the variational linear SDE takes the form q,
(() Ep(()), where the dynamics of jip and Eup are governed by the following ODE
(t) = - (1))) + (1), (0) 140,
(26)
(t) (tt)==A((tEE,(()) ()) Lo(t)Lo(t)', Ex(0)= Eo-
The reparametrization trick in (37] leverages (26) to express Ag and by. in terms of pip. Ex and their time derivatives. This, in turn, allows the drift residual to be rewritten as ra,
(1)
((t)(=(t) - (*)(1) - (0()))), where B(t) = V~ ((Ex(()) @ ((()))-(((oo(((T))) @ denotes the Kronecker sum. the operator V R"xn - F stacks the columns of a matrix into vector, and V-1 18n - R=*= does the opposite. Consequently, by parametrizing the variational distribution 2 tt directly via Hyp and Ep. all the ELBO terms can be evaluated by drawing samples from the multivariate Gaussian .N (t), (t). This allows the ELBO to be maximized without an SDE solver in the training loop, thereby decoupling the computational cost of training from the stiffness of the variational SDE.
(t)(
To address the challenge of learning the variational distribution g, (z () (i.e., (4p (1), (x)(() over long time horizons, we adopt the amortization strategy from [12]. This involves partitioning the time series into shorter segments and learning the variational distribution locally for each segment. To illustrate this approach, consider trajectory with observation times {f1,..... tmm fm+1, fn. J.
where t = to and fn. = T. We partition this sequence into Tim = [nn/m] non-overlapping segments (sub-intervals) as follows
{ I 10
(b)
Th
(27)
k=1
(k-1) (for A > 1), my m for k < 7m and Mn... m. The variational distribution is then amortized over these 7mm segments, leading to the reformulated ELBO:
where {1 = =
((Ø)
my
E
log po
(k) 1 (v 12 )
1
E (20,0((())))) di
20).
DKI. (1p (0) p(0)),
(28)
where 0 ~ qo (0), 2/ ~ (2 (k) ) z(t) ~ (z t), and vj (k) is the j-th observation in the k-th segment, corresponding to time t..
Within each segment In, the parameters of the variational distribution, php and Egg, are defined as follows: At the discrete observation times (k) within the segment, (t, (8)) and (tp are given
(k)
by applying the probabilistic encoder (defined in Section 3) to the observation y) (h) To obtain
continuous-time representation for Mep and Σpp (and their derivatives مم Ein)) for [t((), (k)
which is necessary for evaluating the integral in (28), we employ a deep kernel interpolation scheme
that operates on these encoder outputs, following [12 The trajectory partitioning parameter E N a hyperparameter of the training process.
5 Recall that the initial conditions jity Lig are given by the probabilistic encoder p=(= v(to)), as discussed (.ction(422oot the main paper.
Note that when na is not divisible by m. the last interval has fewer than m elements.
19

---

F
Computing the posterior predictive distribution
After the model parameters are estimated using the multiscale variational inference framework, the primary objective during inference is to generate probabilistic predictions over given time horizon. Specifically, given an initial observation of the full state y(ta) E IR at time to, we seek the posterior predictive distribution Po(y(t) (((aa)) for any >> to. This distribution quantifies our prediction for y(I) and accounts for all sources of uncertainty captured by the model, including the initial state decomposition, the stochastic evolution of the latent variables over the interval [lo, 4], and the final state reconstruction. The prediction process involves three fundamental stages that bridge the fully-resolved observation (1) and the latent state = 1c- .7 -,= E IR=
Stage 1: Probabilistic encoding of the initial state The inference process commences by mapping the initial observation y(to) into the latent space defined by the macroscale state Ç E IR and the microscale state 7/ E R"v, Since our model posits a latent representation, the true initial latent state =to)) corresponding to y((o)) is inherently uncertain. The learned probabílistic encoder po(z v) provides the distribution over possible initial latent states at time to:
p=(z(to) v(to)) = (=(((t) | 2,,)),
(29)
where Mio E RR= and Σt, &R...., are determined by the learned encoder applied to y(to):
(Sx((())) Pto E[z((o) ((6))] (S, (v((t)))) E.o Cov((((6))
Egg
Here, Σg represents the learned covariance matrix associated with the encoder, capturing the un- certainty introduced by the initial scale separation (S) and microscale compression (E7). This initial distribution (2((t) Mtay (t,) serves as the starting point for evolving the latent dynamics over the time interval (to, t].
Stage 2: Evolving latent state uncertainty via the learned multiscale SDE The core of the
model lies the learned coupled system of multiscale Itô SDEs governing the evolution of the latent macroscale (() and microscale (1)) states, represented compactly
dz = dz + (s()))((s), s [fo
(30)
where (2) E R" is the learned drift function, (s) E Rm=X*= is the learned dependent diffusion matrix, and / E R= is standard Wiener process. Propagating the initial distribution po(2(tn) y(ta)) through this SDE from time to to t yields the distribution of the latent state at the prediction time i, denoted by (=(((t) y(t))) It is worth noting that even if the initial distribution (o(((tp) y(fa)) is Gaussian, the distribution (=()(t) (ff)) after evolution through the nonlinear SDE (30) is non-Gaussian Its probability density function, formally governed by the Fokker-Planck equation associated with (30). rarely admits closed-form analytical solution. We will show next how this challenge can be addressed using Monte Carlo sampling.
Stage 3: Probabilistic reconstruction using the PoE likelihood The final stage connects the distribution of the latent state time 1, i.e., p=(2(() v(f)),, back to the observable state space yf) E R This mapping is defined by the likelihood model pe (y] 2) learned during training. As described in Section 4. and Appendix [ this likelihood leverages the PoE perspective enforce scale separation. In brief, the combination of the macroscale expert (r {) Ny (C). Σ
and the microscale expert P2 (y- - ( = (N (c) (D (n), 23) results in a single Gaussian conditional distribution for y given a specific latent state 2:
( 2) A m,(2), E,),
(31)
where the conditional mean (2)22) Ry and conditional covariance E, E By X are
Mo (2(2) = Ely = = D3 E, (1) and Σ, = Cov(y) 2) = (5/ + 5
(32)
Here, S/ = (21))1 and 5% (- are the precision matrices (parametrized by 0) associated with the Gaussian macroscale and microscale experts, respectively, DG. R. - R. is the macroscale decoder, and D = nd denotes the microscale decoder.
20

---

The PoE likelihood introduces two key statistical nuances. Firstly, the conditional mean (y(2) in 32 not simply the sum of the decoded macrostate D2 (C) and the decoded microstate (n). The contribution from the microscale decoder D2 (1) is adaptively weighted by the matrix W = SnyS
(S) + - This weighting factor emerges naturally from combining the Gaussian experts in the PoE framework and depends on the learned relative precisions (5 SQ). This ensures that the microscale state primarily contributes to explaining aspects of y where the macroscale expert has low precision. Secondly, the conditional covariance Σy in (32) which is independent of 2, quantifies the inherent uncertainty in reconstructing y even if were known perfectly. This uncertainty stems from the combination and potential disagreement between the macroscale and microscale experts, as encoded in their respective covariance matrices Dig and Ey.
The final predictive distribution (e(y(() (((tt)) is obtained by marginalizing the intermediate latent state z(t), which yields
pe pe(u(t) = / po po(u(() (x)))
dz(t)
f N(((() (z)())))((((())))))))))))))()))) E.p) po (=((()) (u(o)) dz((),
(33)
where in the second step we used the PoE likelihood (31 The preceding integral is the expecta- tion of the conditional likelihood po po(y(t) ((t) with respect to the distribution of the latent state (o()(() y(ta)) obtained by propagating the initial condition provided at to through the latent SDE (30). Unfortunately, this integral is analytically intractable due to the nonlinearity of the conditional mean (z()) and the non-Gaussian nature of the distribution po (=((() ((o)). In the present work, we use Monte Carlo sampling to approximate this integral and characterize the posterior predictive distribution as follows:
1. Sample initial latent states: Draw N > 1 independent samples from the initial latent distribution defined in Eq. (9). i.c., 24 (to) ~ (t y(to)) = (N,,, Ex,); for i = ...
2. Propagate samples via SDE: Numerically integrate the learned SDE (30) forward in time from to to t for each initial sample 2; (ta). This yields N samples of the latent state { which collectively form an empirical approximation of the intractable distribution pa po(2(t) (tto)) at time t.
Calculate conditional means: For each latent sample (1), compute the corresponding conditional mean of the observed state y(t) using (32). i.e., My = (1) = (C)(()) + x.,, S D2 (1)).
4. Approximate predictive distribution: The posterior predictive distribution (e(((t) [ ((to)) is approx- imated by the empirical distribution derived from these samples. Specifically, can be viewed as Gaussian mixture model with N components: po (v(() ((4)) R + Σ.=,, N(v(t) pyya Σ.y), where each component (()(t) (y.. X.y) in this mixture represents the distribution of y() conditioned on a specific realization 24 (t) of the latent state evolution. All components share the same conditional covariance M, (the reconstruction uncertainty from the PoE), but are centered at potentially different means Hg reflecting the uncertainty propagated through the latent dynamics.
Predictive mean and covariance The overall moments of the posterior predictive distribution can be estimated from the Gaussian mixture approximation. The predictive mean is the expected value of (() given the initial observation y(4o), which can be estimated by taking the sample average of the conditional means My, obtained from each propagated latent state sample 2, ()), i.e.
1 N
1 N
Ppreed Ely(t) v(to)]]
Ely(t) ( 2x)
My.s.
(34)
The predictive covariance which is the total covariance matrix of y(t) given v(to), capturing the overall prediction uncertainty at time can similarly be estimated as:
Epprea(t) = Cov(y(t) (to)) (eco)/Cov(y(()) =(t)]] + Cov k(t) E(y(t) e(t))
N
E
1 N
Σ (My.a Pepred - (t))
(35)
N
E, + CompleCov((/)
21

---

Table 2: POD-SINDy hyperparameters.
Hyperparameter
Wave ID (ng 25)
KdV ID
Burgers 2D
Cylinder 2D
SWE 2D
(n, 25)
(n, = 69)
(n, 517)
(n, 69)
Polynomial order Threshold
1 0.03
0.1 0.01
0.1
0.01
The total predictive covariance is given by the sum of two distinct terms: (i) the average reconstruction uncertainty, Σy = (S) + + which is independent of the latent state propagation reflecting the precision of the combined experts, and (ii) the propagated uncertainty, SampleCov( ([My,) }). which arises from the variability of the latent state samples (z)( incorporating both the initial uncertainty encoded in Eto and the stochasticity introduced by the SDE dynamics over (4n, 4). The proposed Monte Carlo sampling procedure provides a practical and statistically principled method for generating predictions and quantifying their uncertainty at any future time >> to. It fully leverages the structure of the learned stochastic multiscale model and the scale-separating PoE likelihood, without requiring further analytical approximations beyond the finite number of samples N.
G Coarse DNS, POD-SINDY, and DMD baselines
In this section, we discuss the setup for three of our baseline methods in each test case: coarse DNS, DOD-SINDy, and DMD. The coarse DNS baseline involves simply coarsening the computational mesh used to evaluate the full-order model that generated each dataset. This is straightforward for the wave, KdV, and Burgers test cases, as we created those test cases, and for the shallow water case, as the model that generated the data is publicly available in the associated code repository (39. As the original solver for the cylinder flow dataset was unavailable, we implemented the full-order model using FEniCS|| based on the description given by Günther et al. (38 this code is available in the git repository of this paper. A cubic interpolation scheme was used to map the coarse grid DNS predictions to the fine mesh for computing error metrics.
For the POD-SINDy baseline, implemented using PySINDy 47), we first projected the training data onto POD basis. The latent dimension was set equal to the total latent state dimension of our model ensure fair comparison of dynamics at the same level of compression. We use polynomial basis with thresholded ceast-squares to learn the latent dynamics. We performed hyperparameter search over the polynomial order E [1 and the sparsity threshold E {0.11, 0.03, 0.01, 0.003, 0.001), }. selecting the combination that minimized validation error. Due to the high-dimensional latent space for the 2D test cases, the number of basis functions becomes impractically large when using quadratic polynomials. We therefore only consider linear latent dynamics for the 2D test problems. The final hyperparameters for each case are reported in Table
For DMD, we use the package PyDMD to learn a discrete-time linear model. The rank of the learned operator is chosen be n2g and the Tikhonov regularization parameter is fixed at 0.01 for all cases.
H Detailed setups for test cases
In this section, we provide detailed setups for the test cases, as well as their respective multiscale model architectures and training procedures.
Igor A. Baratta, Joseph P. Dean, Jørgen S. Dokken, Michal Habera, Jack S. Hale, Chris N. Richardson, Marie E. Rognes, Matthew W. Scroggs, Nathan Sime, and Garth N. Wells. DOLFINx: The next generation FEniCS problem solving environment, December 2023
8 Sara M. Ichinaga, Francesco Andreuzzi, Nicola Demo, Marco Tezzele, Karl Lapo, Gianluigi Rozza, Steven L. Brunton, and J. Nathan Kutz. PyDMD: A python package for robust dynamic mode decomposition. Journal Machine Learning Research, 25(417):1 2024
22

---

H.1 One-dimensional advecting wave
The data for this test case is generated from the ID advection equation defined over the domain E [0, 1] with periodic boundary conditions:
du Ou
= C
(36)
01
Ox.
where c 1. The initial condition is u(2 0) (xp(-(z/w)²), where w (000001,..02)
To generate the training trajectories, we use a spatial grid with My = 1000 points (Ar = 0.001) and approximate the spatial derivative with centered finitediffference scheme resulting in the system of ODEs
duy df
4j+1 u/-1
(37)
2
for j = 1,2....., Myy where "j denotes the solution at the j-th spatial gridpoint. The system of ODEs is integrated over the interval [0, using time step of A/ 0.001 (n, 1001). The final observational data is created by corrupting these ground-truth trajectories with i.i.d. zero-mean Gaussian noise (a 0.001). The training dataset comprises 20 trajectories, with each for validation and testing. The test set covers shorter time domain of t E [0, 0.2]. For our multiscale models, the macroscale state defined on coarse mesh with ny 20 points.
H.2 Korteweg-de Vries equation
The Korteweg-de Vries (KdV) equation in ID is given by
du
Iu
0 u
Ot
M Dx
dr
0,
(38)
where 0.02 and the domain is x
[0, 10) with periodic boundary conditions. The initial
condition is u(z, 0) cos (TTz)expp exp
(2-7 where (2, 0.01) and 0..00.....
We generate the observational data using spatial grid with n, 1000 points (Arr 0.01) and approximate the spatial derivatives using finite-difference scheme leading the system of ODEs
duj dt
Ui+i U új 1j+2 34z+1 + 3u, ύj
(39)
2Ar
3r3
for j ,,..,,,...,, where uy denotes the solution the j-th spatial gridpoint. The system of ODEs is integrated over the time interval t E [0, 1] using a time step of At = 0.001 (n, = 1001). The observations are corrupted by i.i.d. zero-mean Gaussian noise with standard deviation 0.01. The training dataset is formed using 10 samples of initial conditions and their corresponding trajectories, and the validation and test sets are formed using 5 samples each. For our multiscale models, we use coarse mesh with ni = 20 points for the macroscale state.
H.3 Burgers' equation in 2D
As an additional test case not included in the main body for brevity, we evaluate the models on the 2D viscous Burgers' equation:
Ou
Ou
Ou
u
+ u
= VVuu u.
(40)
0t
dx
Or
The domain is [o, 1j2 with zero Dirichlet boundary conditions and viscosity v = 0.005. The initial (a) 00.)) (az -0.3"" where a (N,00.01 and
condition is given by u(x1, #2, 0) = aexp (0.2, 0.0001).
We use a uniform Cartesian spatial grid with mesh spacing Ar = 1/127 in each direction (ny = 128 16384) and approximate the spatial derivatives using finite-differece scheme leading to
Bengt Fornberg. Generation of finite difference formulas arbitrarily spaced grids. Mathematics of Computation, 51:699-7066 706, 1988
23

---

Table Error statisties for each model type in Burgers 2D test case.
Model type
Error mean + std. dev. on test set
Coarse DNS
0.560 + 0.273
DMD
0.097 + 0.039
POD-SINDY
0.057 + 0.045
Implicit scale
0.141 + 0.040
Our approach
By 1
0.108 + 0.056
Thy = 2
0.035 + 0.022
ny = 3
0.037 0.029
fh = 4
0.047 + 0.039
On 5
0.028 + 0.019
the system of ODEs
dus.j dt
Ui+1.j
U4-1.5 +
Us.j+1
Ui
14.5
2A.
2Ar
(41)
U1+1,
2u4, + 41-1d
U4,+1
-2us.j + us--1
+ v
Ar
Ax
for = 1,..., 128, = I,.... 128, where u4y denotes the solution the j-t spatial gridpoint. We generate ground-trun trajectories by integrating the system of ODEs over the time interval E [o. 1] with time step of At 0.001 (n, 1001). The observations are corrupted by Gaussian noise with standard deviation 0.001. The training dataset is formed using 20 samples of initial conditions and their corresponding trajectories, and the validation and test sets are formed using 5 samples each. For our multiscale models, the macroscale state is defined on a coarse x 8 grid.
As shown in Table 3 our multiscale model with n, 5 achieves the lowest error, surpassing the strongest baseline, POD-SINDy, and reducing the mean prediction error by approximately 51% (from 0.057 to 0.028). The improvement over the implicit-scale model is even more pronounced (an 80% reduction), demonstrating the critical role of the explicit microscale state in capturing sub-grid phenomena. Figure T0 visually confirms this advantage; while the best baselines struggle to resolve the sharp wavefront, our multiscale model captures it accurately. As shown in Figure [1 this is because the microscale state successfully encodes these fine-scale features, complementing the smooth macroscale state to yield a superior prediction.
H.4 Cylinder flow in 2D
The data for this test case is simulation of a 2D Von Kármán vortex street at a Reynolds number of Re 160 generated using the Gerris flow solver [50] by Günther et [38] The dataset provides
the nondimensionalized, two-component velocity field. The original data is defined the spatial
domain (-0.5, 7.5] 0.5, 0.5) with 640 x 80 spatial grid. For our study, we truncate the domain to [-0.5,33 3.5 x [- 0.5, 0.5] resulting in 320 x 80 spatial grid.
The data consists of single trajectory with 1501 snapshots over the time interval & [0, 15). starting from a zero-velocity initial condition. We partitioned this trajectory temporally for training, validation, and testing, using the intervals [0, 13), (13, 14), and (14, 15). respectively. No noise was added to this dataset. For our multiscale models, the macroscale state is defined on coarse 32 x 8 grid, resulting latent state dimension of nę 2 32 x 512 for the two velocity components.
H.S Shallow water equations in 2D
This test case uses the radial dam break problem from the PDEBench dataset (39]. which simulates the 2D shallow water equations. The dataset provides the scalar fluid elevation field [-2.5, 2.5
24

---

Observation
Coarse DNS
DMD
POD-SINDy
Implicit scale
Ours (n1, 5)
Figure 10: Burgers' equation in 2D: prediction comparison on test trajectory at t 1.
Observation
Approx. full state
Prolonged macro
Decoded micro
Figure 11: Burgers' equation in 2D: scale separation with ne 5, shown
spatial domain discretized using uniform 128 128 grid. The initial condition is binary field: ((x, 0) = 2 v||= r, otherwise u(z, 0) = 1, where r ~ (0.3,, 0.7].
The dataset comprises 1000 trajectories, each generated from different initial radius Each [0, 1]. We use 900 trajectories for
trajectory consists of 101 snapshots over the time interval training, 50 trajectories for validation, and 50 trajectories for testing. We preprocess the data by subtracting from the state everywhere to center it, and corrupt by adding i.i.d. Gaussian noise with standard deviation of 0.01. For our multiscale models, the macroscale state is defined coarse x 8 spatial grid (nç- 64).
H.6 Neural network architectures
This section details the neural network architectures for the smoothing operator, encoder, decoder, and drift functions. All hidden layers in these architectures use a LeakyReLU activation function.
H.6.1 Encoder architecture
The probabilistic encoder models the conditional distribution of the latent state given the fully- resolved state y as Po(z y) [[ [=ū 2g). For all test cases, we use diagonal parametrization of the covariance matrix En n The mean of the encoder is given by
(Sx(t)) (S) (v))
(g(() E[2 vl
where Eg and Eg the macroscale and microscale encoders, respectively, So is the smoothing operator, and S (y) Sy) is the residual component containing small-scale features. The architectures for these components are as follows:
Smoothing operator Sy: We parametrize the smoothing operator for each field variable using single-channel convolutional layer with a stride 1. Let (n, / denote the downsampling factor between the fine grid with ny points and the coarse grid with n. points, where is the spatial domain dimension. and 4s 1 for KdV, while for the 2D test cases the kernel size is set to (4s 1) x (4s 1). We use circular padding of length 3s for advecting wave, and 2s for all other test cases.
set the kernel size to 6s + for advecting wave
2. Macroscale encoder E²: The macroscale encoder maps the smoothed fully-resolved ob- servation, Sg (y), from the fine grid with ny points to the coarse grid with n& points. We implement this as a strided convolution that shares the same learnable kernel as the smooth- ing operator but uses stride of 8, the downsampling factor defined above.
3. Microscale encoder E The input to the microscale encoder the residual ğ = y So(y). The architecture of the microscale encoder consists of series of convolutional layers that
25

---

progressively downsample the input (using stride of 2), followed by final linear layer to produce the 7i vector. The details the architecture used for each test case are provided in T For the ID test cases, we use kernel size of 25 and circular padding of length 12. For the 2D test cases, we use kernel size of x and circular padding of length 4.
Table 4: Architecture of the microscale state encoder for each test case.
Wave ID
KdV ID
Burgers 2D
Cylinder 2D
SWE 2D
Conv layers
Conv layers
Conv layers
Conv layers
Conv layers
# filters
#
# filters
# filters
#
# filters
#
# filters
1
1
2
1
4
2
4
64
2
64
64
16
3
16
32
3
32
32
4
64
4
64
4
16
4
16
4
16
Linear layers
Linear layers
Linear layers
Linear layers
Linear layers
In
Out
#
In
Out
#
In
Out
#
In
Out
#
In
Out
8000
nl 5
8000
Thy
4096
5
6400
4096
H.6.2 Decoder architecture
The decoder maps the latent state 2 = (S, 17) to the parameters the conditional distribution po(yl z). As shown in (32) the conditional mean is weighted sum of the microscale and macroscale contributions, i.e., : E[y 2| = DS(C) + Ey (n). while the conditional covariance is given by Cov(yl 2| = ( - 5' We use a diagonal parametrization for the precision matrices S and Sy The macroscale and microscale decoder architectures are described below:
1. The macroscale decoder, Dÿ single convolutional transpose layer with stride equal to the downsampling factor s. Its weights are tied to the macroscale encoder, serving as its approximate inverse to prolong the macroscale state to the fine spatial grid.
2. The architecture of the microscale decoder mirrors that of the encoder. We use linear layer to project the microscale state to small spatial feature map, followed by series of convolutional transpose layers to progressively upsample the microscale features back to the fine grid resolution. These convolutional transpose layers use the same kernel size, stride, and padding settings as their counterparts in the microscale encoder. The specific architecture for each test case is detailed in Table 5
H.6.3 Macroscale drift
As described in Appendix A the macroscale drift at each grid point is parameterized by a feedforward neural network, fe. whose weights are shared across all spatial locations to enforce translation equivariance. This network takes as input the state values from a local stencil, leading to the input dimension 7:m d, 1)d + n where d, is the number of field variables and d is the spatial dimension. For all cases, we the stencil half-width q = and use an identity coupling function leading to n. = n,, microscale features used as input. The specific MLP architecture for each test case is given Table
For the Wave, KdV, Burgers, and SWE test cases, the physical systems are defined on simple domains and spatially homogeneous. The ranslation-equivariann MLP therefore well-suited architecture. The cylinder flow problem, however, is spatially non non-homogeneous due to the presence of the cylinder obstacle. To allow the shared-weights network learn position-dependen dynamics, we augment its input with four spatial positional encoding features. The input to the elementwise drift fo is augmented with (cos(2mx1/L1)) sin sin(2 L1), cos(22.zg L2), sin(2 L2)), where
26

---

Table Architecture of the microscale state decoder for each test case.
Wave ID
KdV ID
Burgers 2D
Cylinder 2D
SWE 2D
Linear layers
Linear layers
Linear layers
Linear layers
Linear layers
#
In
Out
#
In
Out
#
In
Out
In
Out
#
In
Out
8000 Thp
1
8000 fhy Reg
4096 Mn
6400
4096
ConvT layers
ConvTT layers
ConvT layers
ConvT layers
ConvT layers
#
# filters
#
filters
#
# filters
# filters
#
# filters
2
64
2
64
2
16
16
16
3
16
3
16
3
32
32
32
4
4
4
64
64
64
5
5
5
1
5
2
5
1
Table 6: Architecture of the macroscale drift function for each test case.
Wave ID
KdV ID
Burgers 2D
Cylinder 2D
SWE 2D
Linear layers
Linear layers
Linear layers
Linear layers
Linear layers
#
In
Out
#
In
Out
#
In
Out
#
In
Out
#
In
Out
1
Flin
128
Nin
128
fin
128
Tin +
128
fin
128
2
128
128
2
128
128
2
128
256
2
128
256
2
128
256
3
128
128
3
128
128
3
256
128
256
128
256
128
4
128
4
128
4
128
4
128
2
4
128
1
L1 and Le denote the length of the spatial domain along the coordinates and 22, respectively. This results in an input dimension of nia + 4 for the cylinder flow case.
Table Architecture of the microscale drift function for each test case.
Wave ID
KdV ID
Burgers 2D
Cylinder 2D
SWE 2D
Linear layers
Linear layers
Linear layers
Linear layers
Linear layers
In
Out
In
Out
#
In
Out
#
In
Out
#
In
Out
128
128
128
1
128
1
128
128
512
128
512
128
256
128
256
2
128
256
512
128
512
128
3
256
128
3
256
128
3
256
128
128
71
128
128
Thy
128
fg
128
fby
H.6.4
Microscale drift
In our numerical studies, the microscale state is chosen to have spatial grid structure. We therefore parametrize the microscale drift function using fully-connected neural network. We use non-autonomous model for the drift leading to the input dimension 74n n, + +n2 + For the ID cases, the coupling term & identity (nç = ni)). For the 2D cases, v is learned linear map that projects the n ne-dimensional macroscale state to R (n n,,). The architecture for each case is shown in Table 7
27

---

H.7 Training
Each model trained with the Adam optimizer 10 and the loss function the negative of the ELBO defined in (8) The batch size is 64. The expectation terms in the ELBO are estimated via a single Monte Carlo sample from the variational distributions. To approximate the temporal integral term in the El BO, we use 64 quadrature points. The advecting wave models were trained for 50 epochs, KdV for 100 epochs, Burgers for 100 epochs, cylinder flow for 1000 epochs, and shallow water for 20 epochs.
Each multiscale model trained in hierarchical manner. We begin by training an implicit-scale model with an initial learning rate of 10-3. We apply an exponential scheduler that decays the learning rate by 10% every 2000 optimization steps. We then initialize the multiscale model with ny 1 using the trained parameters from the implicit-scale model: we apply Xavier normal initialization to all new weights and zero initialization all new biases. The multiscale model then trained with a reduced initial learning rate of 10 as its loss landscape is more sensitive to large step sizes than that of the implicit-scad scale model. For a multiscale model with n, 2, we repeat this process using the parameters from the trained model with n, 1, and so on for larger may.
Although training model for large ny directly is possible, we found that the highly non-convex loss landscape can lead to suboptimal local minima. The hierarchical initialization strategy acts as form of curriculum learning, guiding the optimization through a more stable path. This procedure was crucial for achieving the near-monotont performance improvements with increasing Rhy reported in Table I
I Additional plots for results
This section contains expanded plots of the results presented in Section 3 Many of these figures expand on plots presented the body of the paper and earlier in the appendix. They are summarized in Table [
Table 8: Additional result plots in Appendix
Figure in Appendix Relevant earlier figure
Description
Scale separation for...
Wave ID case
KdV ID case
Burgers 2D case
Cylinder 2D case
SWE 2D case
Multiscale predictions for...
18
N/A
Wave ID case
19
N/A
KdV 1D case
20
Burgers 2D case
Cylinder 2D case
SWE 2D case
Baseline comparisons for...
24
N/A
Wave ID case
25
N/A
KdV ID case
26
10
Burgers 2D case
27
28
Cylinder 2D case
29
SWE 2D case
10 Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. 3rd International Conference Learning Representations, 2015
28

---

L.1 Scale separation plots
==0
0.05
t = 0.1
1= 0.15
== 0.2
(
y(ti)
+ g(t))
Yi
Figure 12: Advecting wave: visualization of scale separation with 20 and 5. The columns yı observation from dataset, y(4)): reconstructed full state, 7(4): prolonged macroscale state, (tx)): decoded microscale state.
m
0.25
mmm
== 0.5
m
1= 0.75
1=1
w
yt,))
y(ti) + g(ti)
Yi
Figure 13: KdV equation: visualization of scale separation with ni 20 and n, 5. All plots common æ and y-axis scales. The columns are y: observation from dataset, y(,)): reconstructed full state, g(t,): prolonged macroscale state, (t,)): decoded microscale state.
29

---

0
0.25
0.5
t = 0.75
t=1
y(ti)
g(ti)
+
g(ti)
Yi
002
0.0
0.2
0.4
0.6
0.8
1.0
Figure Burgers' equation 2D: visualization of scale separation with nı 64
and n, = 5. The columns are yı observation from dataset, y(4)) reconstructed full state, ((.): prolonged macroscale state, J(t,): decoded microscale state.
30

---

[
0.25
= 0.5
= 0.75
==1
Yi
y y(t
=
y(ti)
+
g(ti)
1.0
0.5
0.0
0.5
1.0
1.5
Figure 15: Cylinder flow in 2D:
visualization of scale separation on velocity z-component with nÇ 32 512 and ny The columns are yı observation from dataset, yt,)): reconstructed full state, 7(t)): prolonged macroscale state, J(t;)): decoded microscale state.
0.25
0.5
t 0.75
==1
Yi
y(th)
=
y
+
g(ti)
1.0
0.5
0.0
0.5
1.0
1.5
Figure 16: Cylinder flow in 2D: visualization of scale separation on velocity y-component with ne 32 512 and n, The columns are y: observation from dataset, y(t)): reconstructed full state, 7(4): prolonged macroscale state, T(t+)): decoded microscale state.
31

---

0
1 0.25
t = 0.5
t= 0.75
t=1 1
y(ti)
y(ti)
+
g(ti)
Yi
-1.0
0.0
0.5
1.0
Figure 17: Shallow water equations in 2D: visualization of scale separation with n. 8 64 and n,, = 5. The columns are observation from dataset, y(t,): reconstructed full state, (t.)): prolonged macroscale state, 7(4.): decoded microscale state.
32

---

1.2 Multiscale model prediction plots
0
0.05
0.1
== 0.15
1= 0.2
Figure 18: Advecting wave: multiscale model prediction on trajectory from test set with nÇÇ = 20 and n,y 5. Observations represented with black solid curves and model predictions with red dashed curves.
w
0.25
0.5
0.75
Figure 19: KdV equation: multiscale model prediction trajectory from test set with ng 20 and 7y = 5. Observations are represented with black solid curves and model predictions with red dashed curves.
33

---

t
0.25
1= 0.5
0.75
t=1
Observation
Prediction mean
Absolute error
Prediction std. dev.
0.2
0.0
0.2
0.4
0.6
0.8
1.0
0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175
Figure 20: Burgers' equation in 2D: multiscale model prediction trajectory from test set with
nc = 8 x 64 and n. = 5.
34

---

0
= 0.25
= 0.5
0.75
t=1
Observation
Prediction mean
Absolute emror
Prediction std. dev.
-1.0
-0.5
0.0
0.5
1.0
1.5
0.00
0.02
0.04
0.06
0.08
0.10
0.12
0.14
Figure 21: Cylinder flow in 2D: multiscale model prediction of z-component of velocity on test interval with ng = 2 x 32 x 512 and ny 5.
0
= 0.25
= 0.5
0.75
1
Observation
Prediction mean
Absolute error
Prediction std. dev.
-1.0
-0.5
0.D
0.5
1.0
1.5
0.00
0.02
0.04
0.06
0.08
0.10
0.12 0.14
Figure 22: Cylinder flow in 2D: multiscale model prediction on test interval with ni = 2 32 8 512 and My Velocity y-component is shown.
35

---

t t=0
f = 0.25
0.5
0.75
t=1
Observation
Prediction mean
Absolute error
Prediction std. dev.
-1.0
-0.5
0.0
1.0
0.0
0.1
0.2
0.3
0.4
0.5
Figure 23: Shallow water equations in 2D: multiscale model prediction trajectory from test set with nę = x 8 = 64 and n, =
36

---

1,3
Comparison of predictions
t=0
t=0.5
t=0.1
1=0.15
A
t=0.2
Multiscale
Coarse
DMD
POD-SINDY
Implicit scale
Figure 24: Advecting wave: Comparison of predictions trajectory from test set. Observations represented with black solid curves and model predictions with red dashed curves. For the multiscale model, the predictions corresponding ng 20 and n, 5 shown.
t=0
1.2.55
1=0.5
1=0.75
Multiscale SDE
Coarse DNS
DMD
POD-SINDY
Implicit scale
Figure 25: KdV equation: comparison of predictions on trajectory from test set. Observations are represented with black solid curves and model predictions with red dashed curves. For the multiscale model, the predictions corresponding to vıg = 20 and n, = 5 are shown.
37

---

Observation
Multiscale SDE
Coarse DNS
DMD
POD-SINDy
Implicit scale
0.2
0.0
0.2
0.4
0.6
0.8
1.0
Figure 26: Burgers' equation in 2D: Comparison of predictions on trajectory from test set. The rows correspond to time instances & 0, 0.25, 0.5, 0.75, from top bottom. For the multiscale model, the mean predictions corresponding to nę = 64 and n. 5 are shown.
38

---

Figure 27: Cylinder flow 2D: comparison of predictions e--oompoeet of velocity on trajectory from test set. The rows correspond time instances
Implicit scale
t= 0, 0..5, 0.5, 0.75, 1 from top to bottom. For the multiscale model, the mean predictions corresponding to ng = 512 and ny = 5 are shown.
1.5
POD-SINDY
1.0
0.5
DMD
0.0
Coarse DNS
0.5
1.0
Multiscale SDE
Ohservation
39

---

Figure 28: Cylinder flow 2D: comparison of predictions of y-coompoeet of velocity on trajectory from test set. The rows correspond time instances t= 0, 0.25, 0.5, 0.75, 1 from top to bottom. For the multiscale model, the mean predictions corresponding to ng. = 512 and ny = 5 are shown.
Implicit scale
1.5
POD-SINDY
1.0
0.5
DMD
0.0
Coarse DNS
-.5
1.0
Multiscale SDE
Ohservation
40

---

Observation
Multiscale SDE
Coarse DNS
DMD
POD-SINDy
Implicit scale
-1.0
-0.5
0.0
0.5
1.0
Figure 29: Shallow water equations in 2D: Comparison predictions on trajectory from test set. The correspond time instances i 0, 0.25, 0.5, 0.75, from top to bottom. For the multiscale model, the mean predictions corresponding to nıç 64 and n, = 5 are shown.
41

---

Learning smoothing kernels
In this work, choose use a convolution to obtain the macroscale state from the full state. With an appropriate kernel size and stride, the convolution smooths the full state and sparsely samples on the coarse macroscale mesh. Pro- longing the macroscale state onto the full mesh accomplished by transposing the convolution used to obtain it. Because this operation is lossy, the prolonged macroscale state loses sub-grid- scale features of the full state.
0.1
0.0
50
100
High-rosolution mesh index
The macroscale smoothing kernel learned for the multiscale model (ng 20, n, 5) trained on the KdV dataset. Note that adjacent gridpoints the macroscale mesh are 50 grid- points apart on the high-resolution mesh.
Figure 30:
the kernel of the convolution was simply cho- sen to be Gaussian, then this operation would amount to low-pass filter. We choose however to leave the kernel parametrized and to infer it from data. The optimal kernel for a particular dataset may turn out to be non-Gaussian, and indeed this is what we observed in our numerical studies. The kernel for the multiscale model (ng 20, n, 5) trained the KdV dataset is shown in Figure 33 Although the training process did not converge on Gaussian kernel, it did converge on a smooth curve.
w
So W
D
Observation
Smoothed state
Macroscale state
Prolonged state
Figure 31: The macroscale smoothing, restriction, and prolongation operators visualized on a snapshot from the KdV dataset.
Figure 31 illustrates the macroscale smoothing, restriction, and prolongation. As our example state, we use the snapshot at = from a test trajectory of the KdV dataset. We see how the smoothed state and especially the macroscale state its coarsened mesh lose the small-scale features present the original snapshot. The learned kernel however reintroduces some small-scale features in the prolongation, determined by the training process to maximize the likelihood defined in Section 4.1
Figure 32 shows spec- tral plot of the example snapshot from Figure 31 compared to the prolonged macroscale state and a multiscale reconstruction. As expected, the decoded macroscale state accurately follows the true spectrum above the spatial Nyquist frequency of the macroscale mesh. Below this fre- quency, however, the small- scale features introduced by the deconvolution approxi- mately follow the true spec- trum, although they cannot be directly resolved by the macroscale mesh. This sug- gests that the training process selects non- aon-Gassiian kernel not just smooth, but also to encode information about the sub- grid scales that can be partially reconstructed during decoding.
10a
10a
Magnitude
10
10-1
Observation
Prolonged Macro
10
A0 Acroscale
Multiscale
Nycuist Freq.
10-3
10a
10a
Frequency
Figure 32: Spectral plot of observation, multiscale reconstruction (ng = 20, ny = 5), and prolonged macroscale state using a snap- shot from the KdV dataset.
42

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: In the abstract and introduction, propose method for learning stochastic multiscale models of spatio-temport systems. The body the paper details this method and applies it to several test problems of that type.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results
be expected to generalize other settings.
It is fine to include aspirational goals as motivation long it is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: The limitations are discussed in the conclusion (Section 6). Specifically, we note that the latent SDE could incorporate a stabilizing prior, that the microscale state could be formed in a manner that respects spatial structure, that partial knowledge of the physics ought to incorporated as statistical prior, and that our approach may applied to more general class of problems than we consider in this work.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model vell-specifiiattin, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested few datasets with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
43

---

3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: The main theoretical contributions of our work are the multiscale model framework and the multiscale likelihood model. We provide our assumptions and full derivations throughout the body and appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: The training methods, the model architectures, and the variational inference approach are all described in detail. The code is also provided.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary either make it possible for others replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily the architecture clearly and fully.
model architecture, the paper should describe
(c) the contribution a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
44

---

(d)
We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification:
We link an anonymized repository with code in Section 5.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips... cc/
for more details.
While we encourage the release of code and data, we understand that this might be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
/ //nips. ltliss.cc/pubbiccuuuddeoo
for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permirted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Answer: [Yes]
Justification: The experimental setting is described in Results (Section 5) and detailed in Appendix H. This information is also in our code, which we made available.
Guidelines:
The answer NA means that the paper does include experiments.
The experimental setting should be presented the core the paper to level detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: For each experiment, report the error mean and standard deviation obtained on the test set, showing that the error distribution of deffiiinnll-dimensiinne multiscale model is significantly lower than that produced by our baseline models.
Guidelines:
45

---

The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by enror bars, confi-
dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that out of range (e.g. negative етror rates).
If error bars are reported in tables or plots. The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification:
The compute resources and training time are reported in Results (Section 5).
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary failed experiments that didn't make into the paper).
9. Code ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer: [Yes]
Justification: We have reviewed the NeurlPS Code of Ethics and conformed to it in conduct- ing our research.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive socictal impacts and negative socictal impacts of the work performed?
Answer: [NA]
46

---

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

---

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

---

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


Reviews:
[{"review_id": "f34hnVbpr7", "reviewer": "Reviewer_wKE9", "summary": "This paper presents an interesting approach to modeling stochastic multiscale problems by separating scales through coupled SDEs. The work introduces a thoughtful application of Product of Experts (PoE) likelihood to determine the separation between macro and micro scales. The experimental evaluation demonstrates promising results across several PDE problems including KdV, Burgers, and Kármán vortex, showing improvements over implicit scale methods.", "questions": "See above.", "limitations": "The paper does not seem to fully discuss its limitations.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "## Strengths\nThe paper tackles an important and challenging problem in computational physics and machine learning. The use of variational inference for scale separation is mathematically sound and well-motivated. The theoretical framework connecting stochastic multiscale modeling with neural operators is noteworthy and could inspire future work in this direction.\n\n## Weaknesses and Questions:\n- Why not different dynamcis: The latent dynamics equation (2) treats macro and micro scale states similarly. The work could benefit from exploring differentiated dynamics—for instance, allowing macro scale components to be more global (perhaps through convolution) while keeping micro scale components more local. This architectural choice could better reflect the inherent nature of different scales.\n- Why not more scales: can this framework be generalized to more than two scales?\n- Baseline Comparisons: While the variational inference approach is interesting, the current baseline could be too weak. When using implicit-scale, one can use high latent dimension. Generally speaking the machine learning community has seen success with larger number of model parameters in recent machine learning-based surrogate models. To fully demonstrate the method's capabilities, comparisons with established machine learning baselines such as FNO and PDE-refiner, which are multiscale in nature, would strengthen the evaluation and provide better context for the results.\n- Problem Complexity: Although the paper is motivated by large-scale turbulence problems, the current PDE test cases may not fully showcase the method's advantages. These problems could potentially be resolved using simpler approaches like finite difference solvers or standard ML-based surrogate models. The work would benefit from evaluation on more challenging turbulence problems, such as Navier-Stokes equations with higher Reynolds numbers, where scale separation treatment such as LES becomes truly necessary and beneficial.", "was_revised": true, "final_justification": "The author answers my question. The new experiments details on the runtime are very helpful. Therefore I will increase the rating from 3 to 4.", "rebuttal": "We thank you for your constructive review. We are encouraged that you found our approach to be **\"interesting,\"** our use of variational inference **\"mathematically sound and well-motivated,\"** and our overall theoretical framework **\"noteworthy.\"** Your accurate summary confirms our shared understanding of the paper's core components.\n\nBefore addressing your specific questions, we would like to briefly re-emphasize the three key contributions of our work, which you correctly identified as tackling an **\"important and challenging problem.\"** Our framework introduces:\n\n1. **A Structured Latent Space for Explicit Scale Separation:** Unlike traditional approaches that learn a single, monolithic latent representation, we explicitly decompose the system into a **coarse-grained macro-state (ζ)** that lives on a physical grid and a **compressed micro-state (η)** that captures unresolved sub-grid features. This is a deliberate structural choice designed to enhance physical interpretability.\n\n2. **A Principled Mechanism for Learning Scale Interactions:** We introduce the **Product of Experts (PoE) likelihood** as a novel learning objective in this context. As derived in Appendix C, this objective provides a principled, data-driven mechanism for \"closure.\" It automatically learns how to weigh the contributions of the macro- and micro-scales by adaptively regularizing the micro-scale, ensuring it only contributes when necessary.\n\n3. **A Fully Stochastic, Continuous-Time Dynamical Model:** The evolution of the coupled (ζ,η) states is governed by a learned **Stochastic Differential Equation (SDE)** system. This provides a formal framework for modeling continuous-time dynamics and propagating multiple sources of uncertainty, moving beyond deterministic, discrete-time maps.\n\nYour thoughtful questions about architectural choices, baselines, and problem complexity give us an excellent opportunity to clarify the scope of these contributions and to present new evidence that strengthens our evaluation. We address each of your points in detail below.\n\n\n## 1. On Differentiated Dynamics for Macro and Micro Scales\n\n> \"The latent dynamics equation (2) treats macro and micro scale states similarly. The work could benefit from exploring differentiated dynamics...\"\n\nThis is an excellent point. While Eq (2) states the dynamics in their most general form, our framework does indeed admit—and in fact, uses—different architectural choices for the macro and micro components to reflect their distinct natures.\n\nAs detailed in Appendix A, we impose a static graph structure on the macroscale drift f_θ. This respects the spatial locality of the coarse grid and allows for direct comparison with semi-discretized PDEs (as shown in Fig. 9). In contrast, the microscale drift g_θ operates on a non-spatial, compressed latent vector η and is modeled with a standard MLP. This architectural differentiation is a key feature of our implementation. Your suggestion of using convolutions for the macro-drift is very insightful—it is a similar idea to our graph-based approach and a promising direction for future refinement.\n\n## 2. On Generalizing to More Than Two Scales\n\n> \"Why not more scales: can this framework be generalized to more than two scales?\"\n\nThis is a very insightful question. In our framework, the primary distinction is between features that are resolved by a chosen coarse computational mesh (the macroscale) and those that are unresolved (the microscale). By this definition, any additional feature scale would fall into one of these two categories.\n\nHowever, your question points to a more ambitious goal. One could certainly envision a hierarchical extension of our framework with multiple, nested coarse grids, creating a sequence of latent states. This would be analogous to a multiresolution analysis and is a compelling, though non-trivial, direction for future work. We have added a note to this effect in our conclusion.\n\n## 3. On Baseline Comparisons and Latent Dimension\n\n> \"The current baseline could be too weak. [...] comparisons with established machine learning baselines such as FNO and PDE-refiner [...] would strengthen the evaluation.\"\n\nThis is a crucial point that touches on the positioning of our work. While FNO and PDE-refiner are powerful, they represent a different modeling philosophy. Our work aims to move beyond simply building flexible models with an inductive bias that equips them with the capacity to learn from data with multiscale features. Instead, our goal is to introduce a principled statistical framework that enables the data-driven discovery of models that learn to explicitly separate the scales. Our framework retains physical interpretability by:\n\n- Defining the macro-state via a learned linear filter (the smoothing operator).\n- Imposing a graph structure on the macro-dynamics, enabling direct comparison with the underlying PDE operators (Fig. 9).\n\nFNO and PDE-refiner are highly flexible, \"black-box\" function approximators that do not share these structural priors. Therefore, while a performance benchmark is interesting, our primary scientific comparison is against our own implicit-scale model, which isolates the contribution of our core idea: the explicit microscale state.\n\n**New Numerical Study: High-Dimensional Implicit-Scale Model:** You raise an excellent point about latent dimension. To test whether the advantage of our multiscale model comes from structure or simply more parameters, we conducted a new experiment on the 2D Cylinder Flow case. We trained a new, higher-capacity implicit-scale model with four times the number of macroscale grid points (1024 vs. the original 256). Additionally, we performed a coarse DNS on an unstructured mesh of approximately the same size as our original coarse mesh. The results are summarized below:\n\n| Model                  | # Gridpoints | Error mean | Std dev   |\n| ---------------------- | ------------ | ---------- | --------- |\n| DNS                    | 278          | 0.227      | 0.001     |\n| Implicit scale         | 256          | 0.105      | 0.001     |\n| Implicit scale         | 1024         | 0.063      | 0.001     |\n| **Multiscale (n_η=5)** | **256**      | **0.031**  | **0.001** |\n\nThis result shows that while quadrupling the capacity of the implicit-scale model provides a clear benefit, it is still twice as inaccurate as our much smaller, more structured multiscale model. This strongly supports our central claim that the explicit separation of scales is the key driver of performance, not just raw parameter count. This also shows that the coarse mesh compromises the predictive accuracy of DNS substantially, as it is twice as inaccurate as the equivalent-resolution implicit scale model.\n\n## 4. On the Complexity of the Test Problems\n\n> \"The work would benefit from evaluation on more challenging problems, e.g. turbulence.\"\n\nWe agree that scaling to more complex turbulence flow models is the ultimate goal. However, the primary focus of this paper is to introduce and rigorously validate the theoretical and algorithmic framework. The chosen test cases, while classical, are ideal for this purpose as they exhibit the key multiscale phenomena (wave-breaking, shock formation, vortex shedding) needed to clearly demonstrate and diagnose our scale separation mechanism.\n\nResults obtained for the cylinder flow case demonstrates the ability of our methodology to handle complex, spatio-temporal dynamics. The fact that our method succeeds even when the coarse grid cannot resolve the cylinder itself is a strong testament to its ability to learn and represent critical sub-grid phenomena, suggesting it is well-poised for more complex fluid dynamics problems.\n\n## 5. On Discussing Limitations\n\n> \"The paper does not seem to fully discuss its limitations.\"\n\nThank you for pointing this out. We agree that consolidating these points is important for clarity. We will add a dedicated \"Limitations\" section to the conclusion of the revised manuscript. This section will summarize and elaborate on the following key points:\n\n- **Validity of the Markovian Assumption:** As we argue in detail in Appendix B, our (ζ,η) framework is explicitly designed to be a powerful, learned Markovian approximation of the true, potentially non-Markovian coarse-grained dynamics. However, this remains an approximation, and its validity is system-dependent.\n\n- **Latent SDE Stability:** As noted in our results, we did not use stability-preserving priors on our learned dynamics, which may be necessary for more complex systems.\n\n- **Microscale State Structure:** The macro-state ζ has a clear spatial structure, while the micro-state η is a \"global\" compressed vector. Formulating a more spatially-aware microscale representation is a key avenue for future work.\n\n- **Gaussian Likelihood:** Our use of a Gaussian PoE likelihood, while enabling a concise derivation, may be too restrictive for systems with highly multi-modal behavior (e.g., high-Re turbulence).\n\nWe believe this dedicated section will provide a more transparent and comprehensive overview of the current boundaries of our framework. We are grateful for your feedback, which has helped us to significantly improve our manuscript."}, {"review_id": "QiZdleCv7q", "reviewer": "Reviewer_g7k9", "summary": "The proposed method is concerned with learning stochastic dynamics from observed data. The central idea is to de-compose the observed evolution into a coarse macroscopic one $(\\bar y)$, and a fine micro-scale ($\\tilde y = y - \\bar y$). Once learned, a coupled SDE can be used to evolve both macro- and microstate, including their learned interactions. The method is purely data-driven, which sets it apart from numerous existing closure-type works that model sub-grid effects using neural architecture while relying on classical (linear) model-order reduction methods on the coarse grid (e.g. [1] and the references therein). The performance of the method is compared to a modified version where the subgrid evolution is not explicitly handled.", "questions": "The abstract claim that the learned multi-scale model achieves \"superior predictive accuracy compared to direct numerical simulation (...) at equivalent resolution\". When making claims about comparison to the full order model, I would ask for a comparison between run-time of the FOM (at high/moderate resolution) as well as the online/offline cost of the proposed method.\n\nI believe the paper needs to compare to at least one competing method to make a compelling argument. A classical POD method might already be a good candidate.", "limitations": "The statement that the paper \"fundamentally reimagines how to model multi-scale phenomena\" is, in my opinion, too strong.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "The paper is well written and I found it easy to follow. The selfcontained nature of the paper (together with the appendix) is appreciated. All plots, tables and math are clear. The code used to generate the experiments is provided in its entirety for the review process already. The general idea to leverage multi-scale modeling is compelling, albeit not new.\n\nThe topic of neural emulators and data-driven methods to accelerate the simulation of fluid dynamics problems is clearly of big interest, as evidenced by the numerous works on this topic in recent years (although one should also mention more skeptical voices as in [2]).\n\nThe central weakness of the paper is that I am not convinced that the claims made are truly backed up by the numerical results. The examples shown (KdV equation, Burger's equation, flow past a cylinder) are classical and can be solved to high precision by most non-linear projection methods (e.g. registration methods as in [3]) or, in the case of the flow past a cylinder, classical POD methods [4]. Especially for the latter, I am not convinced that the proposed method outperforms a linear projection method.\n\nA comparison to other multi-scale inspired projects such as GraphCast [5] and [6], which operate on much larger problems, is not given.\n\nReferences:\n\n[1] A. Ivagnes, G. Stabile, G. Rozza: Data-driven Closure Strategies for Parametrized Reduced Order Models via Deep Operator Networks. https://arxiv.org/abs/2505.17305v1 (2025)\n\n[2] N. McGreivy and A. Hakim: Weak baselines and reporting biases lead to overoptimism in machine learning for fluid-related partial differential equations. Nature Machine Intelligence (2024)\n\n[3] T. Taddei: Compositional Maps for Registration in Complex Geometries. SIAM Journal on Scientific Computing (2024)\n\n[4] D. Deep, A. Sahasranaman, S. Senthilkumar: POD analysis of the wake behind a circular cylinder with splitter plate. European Journal of Mechanics - B/Fluids (2022)\n\n[5] R. Lam et al: Learning skillful medium-range global weather forecasting. Science (2023)\n\n[6] Oskarsson et al.: Probabilistic Weather Forecasting with Hierarchical Graph Neural Networks. Neurips (2024)", "was_revised": true, "final_justification": "I find this paper to be well done: It is reproducible, clear, and self-contained.\n\nAt the same time, I believe it over-states its limited contribution. I am not convinced it is outperforming well-established methods on the selected benchmarks. However, this is not necessary if the method is novel and interesting. \n\nUltimately, an application to larger problems is needed to properly judge the benefits of this approach over the full order solvers, which do not struggle with the academic examples tackled in this work.\n\nOverall, I want to acknowledge the quality of the paper more than the benchmark results. I have increased the score, but I hope that the authors give a fair discussion of limitations.", "rebuttal": "We thank you for your detailed and thoughtful review. We are particularly pleased that you found the paper \"well written,\" \"easy to follow,\" and that you appreciated the self-contained nature of our manuscript.\n\nFurthermore, we appreciate your reference to the work of McGreivy and Hakim. We share their perspective on the importance of robust evaluation and the dangers of \"overoptimism\" in this rapidly advancing field. It is in this spirit of transparency and rigor that we provided our full codebase and data generation scripts from the outset. We believe that making our work fully reproducible is a crucial step in allowing the community to rigorously test our claims, verify our results, and build upon our framework for broader benchmarking studies.\n\nWe agree with your summary assessment of our work. Our central idea is to decompose the observed evolution into a coarse macro-state and a fine micro-state, governed by a coupled system of SDEs, and you astutely note that our **\"purely data-driven\"** approach **\"sets it apart from numerous existing closure-type works\"** that rely on classical methods. This understanding forms the perfect basis for our discussion.\n\nYou also make the crucial observation that the **\"general idea to leverage multi-scale modeling is compelling, albeit not new.\"** We wholeheartedly agree. The term \"multiscale\" is indeed widely used in various contexts. Our work's primary contribution is to move beyond simply building flexible models with an inductive bias that equips them with the capacity to learn from data with multiscale features. Instead, our goal is to introduce a principled statistical framework that enables the data-driven discovery of models that learn to explicitly separate the scales. To the best of our knowledge, data-driven learning of coupled stochastic dynamics with explicit scale separation is a novel contribution to the field.\n\nYour critical questions about our evaluation provide a welcome opportunity to present numerical results for additional baseline methods and to clarify the specific claims we make within this new framework. We address each of your points in detail below.\n\n## 1. On the Central Weakness: Positioning, Claims, and Comparison to POD\n\n> \"The central weakness of the paper is that I am not convinced that the claims made are truly backed up by the numerical results. The examples shown [...] can be solved to high precision by [...] classical POD methods [4]. Especially for the latter, I am not convinced that the proposed method outperforms a linear projection method.\"\n\nYou raise an excellent point about benchmarking against Reduced Order Model (ROM) approaches based on POD or nonlinear manifold representations. We agree that the problems we considered can be solved with linear and nonlinear Galerkin ROM methods. We would like to clarify our paper's positioning and the rationale behind the benchmarking studies. Our framework represents a different modeling paradigm from traditional ROMs. A key distinction is that traditional ROMs (like POD or autoencoder-based methods) learn a mapping to a global, non-spatial latent vector. Our approach, in contrast, learns a spatially-distributed representation on a coarse grid (ζ). This makes our method a data-driven stochastic coarse-grained model, a paradigm more akin to Large-Eddy Simulation (LES) in computational physics.\n\nViewed through this lens, our \"implicit-scale\" model (n_η=0) serves as a data-driven closure model, and the central scientific question we investigate is: **Within this coarse-graining framework, does augmenting the model with an explicit microscale state (η) improve full-state reconstruction?** Our ablation study directly answers this question. The results in Table 1, showing a 70% error reduction for the cylinder flow when η is added, provide unambiguous evidence for the value of our explicit multiscale structure.\n\nWe present two additional baseline experiments on the cylinder flow example. We performed a coarse DNS on an unstructured mesh of approximately the same size as our coarse mesh, and we trained an implicit scale model on a finer mesh of 1024 gridpoints compared to the original 256. The results are as follows:\n\n| Model                  | # Gridpoints | Error mean | Std dev   |\n| ---------------------- | ------------ | ---------- | --------- |\n| DNS                    | 278          | 0.227      | 0.001     |\n| Implicit scale         | 256          | 0.105      | 0.001     |\n| Implicit scale         | 1024         | 0.063      | 0.001     |\n| **Multiscale (n_η=5)** | **256**      | **0.031**  | **0.001** |\n\nThis shows that while quadrupling the capacity of the implicit-scale model provides a clear benefit, it is **still twice as inaccurate as our much smaller, more structured multiscale model**. This provides strong evidence that the **explicit separation of scales** is the key driver of performance, not just raw model capacity. This also shows that the coarse mesh compromises the predictive accuracy of DNS substantially, as it is twice as inaccurate as the equivalent-resolution implicit scale model.\n\n## 2. On Comparison to Large-Scale Models (GraphCast)\n\n> \"A comparison to other multi-scale inspired projects such as GraphCast [5] and [6] [...] is not given.\"\n\nThank you for this comment. We see our work as complementary to these landmark models. As a coarse-grained model, our framework's philosophy differs from GraphCast. GraphCast learns a powerful mapping but still explicitly represents and computes on all of its meshes, including the finest one. In our framework, the dynamics are computed entirely on the coarse grid for the macroscale state (ζ) and in the compressed latent space for the microscale state (η). The high-resolution grid is only a target for the final decoding step, not part of the dynamics simulation itself. Our focus is on developing and validating this physically-interpretable structure, which could potentially inform future large-scale models.\n\n## 3. On the \"Superior to DNS\" Claim and Runtime Comparison\n\n> \"The abstract claim that the learned multi-scale model achieves 'superior predictive accuracy compared to direct numerical simulation (...) at equivalent resolution'. I would ask for a comparison between run-time...\"\n\nThis is an excellent point. The claim was intended to compare our model to a traditional numerical solver operating on the same coarse macroscale grid. As our new DNS experiment in the table above shows, a standard solver on such a grid is highly inaccurate.\n\n**To provide the concrete cost comparison you requested, we have quantified the number of function evaluations (NFE) for the KdV 1D case as an example:**\n\nThe offline training cost depends on the microscale dimension, as we train the multiscale model in a hierarchical manner. We first train an implicit scale model (48 minutes), then reusing those parameters and augmenting weight/bias matrices where necessary, we train a multiscale model with n_η=1 (56 minutes), then with n_η=2 (56 minutes), and so on. The training cost in minutes comes out to 48 + 56*n_η.\n\nFor the online cost below, the FOM was solved using the LSODA method from scipy.integrate.solve_ivp, which switches between Adams and BDF based on stiffness detection. This was found to be necessary for KdV. For the implicit scale and multiscale models, Euler-Maruyama was used. All methods had rtol=1e-3, atol=1e-5.\n\n| Online computational cost | Cost in NFE per traj. |\n| ------------------------- | --------------------- |\n| FOM                       | ~39,000               |\n| Implicit scale            | ~2,800                |\n| **Multiscale (n_η=5)**    | **~1,900**            |\n\nThis represents a ~20× reduction in computational cost during online simulation. We will clarify the claim in the abstract and add this runtime analysis to the manuscript.\n\n## 4. On the Strength of the \"Reimagines\" Claim\n\n> \"The statement that the paper 'fundamentally reimagines how to model multi-scale phenomena' is, in my opinion, too strong.\"\n\nWe agree. Upon reflection, this phrasing is overly strong. We will revise this in the conclusion to be more precise, for example: \"we have introduced a new framework for learning stochastic multiscale models from data, offering a distinct, data-driven perspective on the classical problem of closure.\" Thank you for this valuable feedback on tone.\n\n## Closing Remarks\n\nYour critical questions have been invaluable in helping us sharpen our paper's core message and strengthen its evaluation. We hope our detailed response and new results addressed your concerns."}, {"review_id": "1giQnHWkpp", "reviewer": "Reviewer_7ou3", "summary": "In this paper, authors propose to learn stochastic multiscale models in the form of stochastic differential equations (SDEs) from observations to better adapt to those various resolutions in real dynamical systems. By optimizing with a modern forward-solver-free amortized variational inference method, the proposed method achieves superior predictive accuracy compared to previous ones.", "questions": "Please refer to the Weaknesses part.", "limitations": "Yes.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n1. The idea of explicitly modeling multiscale dynamics within systems using coupled SDEs is novel and valuable for dealing with real applications which have uncertain resolutions.\n2. Authors use a forward-solver-free amortized variational inference method to improve efficiency. Instead of using an SDE numerical solver, they adopt a variational approximation and reparameterize the drift and diffusion terms of the SDE to avoid the integration process. \n3. Detailed mathematical proofs and sufficient ablation studies make the claims more convincing.\n\nWeaknesses:\n1. As authors claimed, the latent state is divided into a macroscale representation and a microscale one. However, how to decide the real resolution corresponding to them still needs to be considered. E.g., in Sec.1, authors use the example of atmosphere modeling, if the macro scale corresponds to planetary waves (5000km) and the micro scale corresponds to cloud microphysics (< 1m), what about modeling at 10m and 100m? Is there a need for more scales? \n2. In the experimental part, it seems that all the test sets are synthetic ones, it could be better if the method can be evaluated on real data. Besides, the baseline is only implicit-scale models. According to the remarks given on closure models and implicit-scale models in Sec.3, if the closure models predict state matching true state y instead of $\\bar{y}$, it is equivalent to the implicit-scale model, correct? If not, authors should also include closure models into baselines.", "was_revised": true, "final_justification": "This paper proposes explicit multi-scale dynamic modeling through coupled SDEs, a relatively novel and effective approach. The authors provide detailed derivations and descriptions for easy understanding, and demonstrate their effectiveness through original and supplementary experiments. Although the data are synthesized, they demonstrate to some extent that the proposed method surpasses existing methods. Regarding the choice of hyperparameters (scale), the authors also commit to considering this as a future direction. Overall, the novelty and practicality of the proposed method are sufficient to overlook the issues I mentioned in the Weaknesses section. Taking all factors into consideration, I believe this paper exceeds the acceptance bar of NeurIPS, and I will maintain my original rating of 5.", "rebuttal": "We are grateful for your very positive and encouraging review. Thank you for recognizing the novelty and value of our core ideas, including the explicit multiscale SDE modeling, the efficiency of our solver-free variational inference scheme, and the rigor of our mathematical derivations and ablation studies.\n\nYour questions are insightful and touch upon important practical considerations and the relationship between our framework and traditional closure modeling. We elaborate on these points below.\n\n## 1. On Handling More Than Two Scales and Deciding Resolution\n\n> \"...how to decide the real resolution corresponding to them still needs to be considered. ... Is there a need for more scales?\"\n\nThis is a key question relevant to real-world applications, and you have raised a core design choice of our methodology.\n\nYou are right that we do not independently choose two scales. Instead, the scale separation is defined by a single cutoff scale: the resolution of the coarse computational mesh. As you suggest, this coarse grid resolution is a hyperparameter chosen by the user based on the desired level of resolved detail and available computational resources. In our framework:\n\n- The macro-state ($\\zeta$) represents all features large enough to be resolved on this coarse grid.\n- The micro-state ($\\eta$) is tasked with representing the entirety of the unresolved sub-grid spectrum—all features smaller than this cutoff.\n\nFor a system with a continuous spectrum of scales, this binary split is a pragmatic and powerful first approximation. However, your question about a \"mesoscale\" points to an important extension. One could certainly envision a hierarchical version of our framework with a sequence of nested coarse grids, creating multiple latent states corresponding to different bands of the physical spectrum. This is a compelling, though non-trivial, direction for future work, and we have added a note to our conclusion to this effect.\n\n## 2. On Evaluation Strategy, Baselines, and Real Data\n\n> \"it seems that all the test sets are synthetic ones, it could be better if the method can be evaluated on real data. Besides, the baseline is only implicit-scale models. [...] if the closure models predict state matching true state y instead of $\\bar{y}$, it is equivalent to the implicit-scale model, correct?\"\n\nThese are both great points that allow us to clarify and strengthen our evaluation strategy.\n\n- **Evaluation on Real Data:** We agree that evaluation on real-world observational data is the ultimate goal. For this initial methodological paper, we chose to use synthetic data from well-understood PDEs. This provides a known \"ground truth,\" which allows for the precise and unambiguous error quantification necessary to validate that our proposed mechanisms are working exactly as intended. Now that the framework has been validated, applying it to real observational data is a key next step.\n\n- **Relationship to Closure Models:** We agree with your assessment about the relationship between closure models and our implicit-scale baseline. A traditional closure model's objective is to predict the dynamics of the *filtered* state ȳ. Our implicit-scale model is tasked with the significantly harder and more practical problem of reconstructing the **full, high-resolution state y** from only the coarse latent state ζ. Therefore, our implicit-scale model is a stronger and more relevant baseline for our full-state prediction task. We will clarify this distinction in the revised manuscript.\n\n- **Strengthening the Baseline Evaluation:** While the implicit-scale model is our primary scientific control, we took the opportunity during this rebuttal period to benchmark against other relevant methods to further contextualize our results. We conducted two new sets of experiments:\n\n  1. **Comparison with Coarse-Grid DNS:** We compared our model against a standard numerical solver run directly on the coarse grid. For the KdV equation, a **Coarse DNS** on a 20-point grid yielded a high normalized error of **0.574**, whereas our multiscale model achieves an error of **0.044**. This demonstrates that simply using a coarse grid is insufficient; the model must learn from high-resolution data to correct for truncation errors, which our framework successfully does.\n\n  2. **Implicit-Scale Model with Increased Capacity:** To test whether our multiscale model's advantage comes from its structure or simply a higher number of latent variables, we trained a new, higher-capacity implicit-scale model for the 2D Cylinder Flow case with four times the number of macroscale grid points (1024 vs. the original 256). For comparison, we also performed a coarse DNS on an unstructured mesh of approximately the same size as our coarse mesh. The results are:\n\n| Model                  | # Gridpoints | Error mean | Std dev   |\n| ---------------------- | ------------ | ---------- | --------- |\n| DNS                    | 278          | 0.227      | 0.001     |\n| Implicit scale         | 256          | 0.105      | 0.001     |\n| Implicit scale         | 1024         | 0.063      | 0.001     |\n| **Multiscale (n_η=5)** | **256**      | **0.031**  | **0.001** |\n\nThis result shows that while quadrupling the capacity of the implicit-scale model provides a clear benefit, it is **still twice as inaccurate as our much smaller, more structured multiscale model**. This provides strong evidence that the **explicit separation of scales** is the key driver of performance, not just raw model capacity. This also shows that the coarse mesh compromises the predictive accuracy of DNS substantially, as it is twice as inaccurate as the equivalent-resolution implicit scale model.\n\nWe believe these new results further strengthen our claims and provide a more comprehensive validation of our framework.\n\nThank you again for your thoughtful and supportive review. Your perceptive questions have been instrumental in helping us refine the positioning and clarity of our contributions."}, {"review_id": "dTIMQUvF7A", "reviewer": "Reviewer_xrse", "summary": "The authors introduce a data-driven framework for learning stochastic multiscale models of complex physical systems, comprising three core components: an encoder, a latent stochastic differential equation model, and a decoder. Large-scale features are explicitly separated into macroscale and microscale latent states, with their dynamics modeled via coupled SDEs and reconstructed into the high-dimensional space. By leveraging a Product of Experts likelihood and variational inference, the framework achieves efficient training. The authors demonstrate the accuracy across multiple partial differential equation systems.", "questions": "1. The authors discuss reduced-order modeling (ROM) and closure modeling in lines 51–66. I agree the framework differs from closure models, as they focus on coarsened state predictions, while this approach aim to accelerates full-state simulation via latent dynamics. However, doesn’t this framework fall within the ROM category? It seems the proposed framework is a form of ROM, and you can simply called the baseline \"implicit-scale model\" as ROM. How does the proposed method distinguish itself from traditional ROM approaches?? \n2. Can the framework be applied to PDEs with nonuniform spatial discretization?\n3. The latent SDE is described in eq 2. what is the aim of sec 4.2, where a linear SDE in eq 7 is used to approximate the latent SDE in eq 2? \n4. l.101 can you explain why $\\mathbb{E}[\\bar{y}|y] = \\mathcal{S}_{\\theta}(y)$ and how do you guarantee this? \n5. For the experiment parts, will the stochastic multiscale framework give totally the same result as the implicit-scale model when $n_{\\eta}=0$? Is the framework and training procedure of the multiscale framework and baseline totally the same, where the only difference is $n_{\\eta}$? \n6. is it difficult to train with the multiscale model with ELBO in eq 8? \n7. I’m surprised by the poor performance of the baseline in Figure 7 for the cylinder flow experiment. I assume this is not a very tough task when the dimension of latent macroscale state is 512. Can you show the result of the following baselines on the cylinder flow system as in fig 7: naive coarse-scale simulation of the PDE system on 32 * 8 grids; dynamic mode ecomposition with 512 modes?", "limitations": "yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths**\n1. The paper tackles the critical challenge of modeling multiscale physical systems \n2. The paper explicitly separates the microscale and macroscale latent states, enhancing physical interpretability and model accuracy\n3. visual explanation: fig 4, 6, 8 clearly visualize the macroscale state and microscale state \n\n**Weaknesses**\n1. The idea of a multiscale framework comprising of encoder, latent dynamics, and decoder is generally not novel. [1] [2] [3]\n2. Clarity issues: \n      -  l72-76: the definition of fully-resolved state $y$ is a bit confusing. This part can be improved for better comprehension. \n      -  Brief introductions to Large-Eddy Simulation (LES) and Product of Experts (PoE) would enhance understanding\n      -  The notation $S_{\\theta}, D_{\\theta}, T_p, T_r, \\mathcal{E}_{\\theta}$ is complex. Are $\\theta, p, r$ all trainable? Is there any special reason to use $p, r$ instead of $\\theta$ here?  Adopting symmetric notation for microscale and macroscale operators would enhance clarity\n3. See question 1 \n\n--------\n[1] Chen, Peter Yichen, et al. \"CROM: Continuous reduced-order modeling of PDEs using implicit neural representations.\" arXiv preprint arXiv:2206.02607 (2022).\n\n[2] Boral, Anudhyan, et al. \"Neural ideal large eddy simulation: Modeling turbulence with neural stochastic differential equations.\" Advances in Neural Information Processing Systems 36 (2023): 69270-69283.\n\n[3] Stachenfeld, Kim, et al. \"Learned simulators for turbulence.\" International conference on learning representations. 2021.", "was_revised": true, "final_justification": "This paper is well-motivated, with a clear separation between microscale and macroscale latent states to enhance physical interpretability and model accuracy.\n\nOne limitation of the paper is that all experimental data are synthetic. Applying the proposed method to a large-scale real-world problem would further improve the paper.\n\nThe authors have addressed my key concern regarding the performance comparison between their method and linear ROMs on the existing experiments. They included an additional experiment involving a system with slowly decaying Kolmogorov n-width to demonstrate the superior performance of their framework. The authors have committed to incorporating results for this test case and a new fluid flow example in the final version.\n\nBased on these improvements, I have increased the score from 3 to 4.", "rebuttal": "We are grateful for your thorough and insightful review. We are particularly encouraged that your summary and assessment of strengths align perfectly with the core contributions we aimed to make. As you correctly identified, our work introduces a framework that:\n\n- Tackles the **\"critical challenge of modeling multiscale physical systems.\"**\n\n- **\"Explicitly separates the microscale and macroscale latent states,\"** which, as you note, **\"enhanc[es] physical interpretability and model accuracy.\"**\n\n- Leverages a **Product of Experts (PoE) likelihood** to enforce scale separation and variational inference with a Gauss-Markov reparametrization trick for efficient training.\n\n- Provides clear **\"visual explanation\"** of this scale separation in practice, as demonstrated in Figures 4, 6, 8, which you rightly point out as a key strength.\n\nThank you for your thoughtful questions concerning novelty, positioning relative to ROMs, and the strength of our baselines. These are crucial points that we address below, clarifying the positioning of our work and elaborating on our key contributions.\n\n## 1. On Novelty and Relationship to ROM (Weaknesses & Question 1)\n\n> \"The idea of a multiscale framework comprising of encoder, latent dynamics, and decoder is generally not novel. [1] [2] [3] ... doesn't this framework fall within the ROM category? ... How does the proposed method distinguish itself from traditional ROM approaches??\"\n\nThis is a fundamental question of positioning, and we appreciate the opportunity to clarify our contribution. We agree that our method can be broadly categorized as a learned surrogate model, but we argue that it represents a different paradigm from traditional Reduced Order Models (ROMs), including the cited works.\n\nA key distinction is that traditional ROMs learn a mapping to a global, non-spatial latent vector. Our approach, in contrast, learns a spatially-distributed representation on a coarse grid (ζ). This makes it a data-driven stochastic coarse-grained model, a paradigm more akin to Large-Eddy Simulation (LES) in computational physics. This preserves spatial locality and provides a more interpretable latent space.\n\nOur novelty lies in three specific, physically-motivated contributions built on this paradigm:\n\n- **A Structured Latent Space:** We augment the coarse-grained state $\\zeta$ with an explicit, separate latent state for the unresolved sub-grid scales ($\\eta$). This is distinct from [1, 3] which use a single latent space, and a different formulation from the implicit closure in [2].\n\n- **A Principled Learning Objective for Scale Separation:** We introduce the Product of Experts (PoE) likelihood, a novel objective in this context that automatically enforces a scale hierarchy by adaptively regularizing the microscale contribution (Appendix C).\n\n- **A Fully Stochastic, Continuous-Time Model:** We embed this structure within a coupled Stochastic Differential Equation (SDE) system, providing a formal framework for modeling continuous-time dynamics and uncertainty.\n\nWe will revise the manuscript to clarify this positioning, emphasizing our contribution as a data-driven, stochastic coarse-grained modeling framework.\n\n## 2. On Clarity and Notation\n\n> [Comments on notation and clarity]\n\nWe appreciate you pointing out these areas for improvement. We will implement the following changes in the revised manuscript:\n\n- **On fully-resolved state y:** We will expand our definition to be clearer: \"The fully-resolved state y is the system state projected onto a fine spatial grid, which is sufficiently resolved to enable accurate direct numerical simulation and serves as our ground-truth observation.\"\n\n- **On LES/PoE:** We will add brief introductory sentences for both concepts to improve readability and make the paper more self-contained.\n\n- **On Notation:** We agree the notation can be harmonized. Following your suggestion for symmetry, we will revise the notation to use $\\mathcal{E}^\\zeta$  and $\\mathcal{D}^\\zeta$ for the (non-parametrized) macroscale restriction and prolongation operators, and $\\mathcal{E}^\\eta_\\theta$ and $\\mathcal{D}^\\eta_\\theta$ for the (parametrized) microscale encoder and decoder. This clarifies both their function and their trainable nature. To be precise, $\\theta$ represents the set of all learnable model parameters. We will also clarify that $\\mathbb{E}[ \\bar{y} | y] = S_\\theta(y)$ by definition denotes the mean of the learned smoothing operator.\n\n## 3. On the Variational SDE (Question 3)\n\n> \"what is the aim of sec 4.2, where a linear SDE in eq 7 is used to approximate the latent SDE in eq 2?\"\n\nThis is a crucial point we will clarify. These two SDEs play distinct roles in our variational inference scheme:\n\n- **Eq 2 (Nonlinear SDE):** This is the generative model for the latent dynamics. It represents our scientific goal and is the SDE we simulate at test time to make predictions.\n\n- **Eq 7 (Linear SDE):** This is a tool we use for training. It defines a simpler, tractable variational distribution that approximates the true (intractable) posterior over latent trajectories. Using a linear SDE here enables us to leverage an efficient solver-free reparametrization trick from [33, 42]. This enables significant speedups in training by circumventing the need for an SDE solver in the optimization loop. This stands in contrast to other approaches, such as the work by Boral et al. (NeurIPS 2023) which you mentioned ([11] in our submission) which requires a SDE solver to calculate the ELBO and its gradients. Our approach avoids this bottleneck and the potential numerical instabilities that can arise when the stiffness of the SDE increases during the optimization parameter updates.\n\nWe will revise Section 4.2 to state this \"model vs. tool\" distinction more explicitly.\n\n## 4. On Experimental Details and Baselines (Questions 2, 5, 6, 7)\n\n> [Questions about experimental details]\n\nWe are happy to clarify these important experimental details.\n\n- **Application to Nonuniform Grids?** Yes, absolutely. As described in Appendix A, the macroscale state $\\zeta$ is defined on a graph $G=(V,A)$, which naturally handles non-uniform and unstructured discretizations. This is a key strength of our framework.\n\n- **$n_\\eta=0$ Model Equivalence:** Yes, you are exactly right. The implicit-scale model is the multiscale framework with $n_\\eta =0$. We chose this baseline to create the most direct and controlled ablation study possible, isolating the benefit of the explicit microscale state.\n\n- **Training Difficulty:** The ELBO is indeed non-convex. We found training to be stable by employing a standard curriculum learning strategy where the KL divergence terms in the ELBO are initially down-weighted and gradually annealed to full strength. This ensures the model first learns a good reconstruction before regularizing the latent space.\n\n- **Performance of the Baseline in Figure 7:** Your surprise at the baseline performance is well-founded and highlights the difficulty of the task. As you noted, a 512-dimensional latent state seems large, but for a coarse-grained model, this information is spread over a coarse 32 × 8 grid, which is too sparse to resolve the cylinder's boundary layer. This is why the predictions look pixelated.\n\nTo put this in context, we performed two additional baseline experiments on the cylinder flow data. We performed a coarse DNS on an unstructured mesh of approximately the same size as our coarse mesh, and we trained an implicit scale model on a finer mesh of 1024 gridpoints compared to the original 256. The submission results and new results are given below:\n\n| Model                  | # Gridpoints | Error mean | Std dev   |\n| ---------------------- | ------------ | ---------- | --------- |\n| DNS                    | 278          | 0.227      | 0.001     |\n| Implicit scale         | 256          | 0.105      | 0.001     |\n| Implicit scale         | 1024         | 0.063      | 0.001     |\n| **Multiscale (n_η=5)** | **256**      | **0.031**  | **0.001** |\n\nWe did similar additional studies with Burgers:\n\n| Model                  | # Gridpoints | Error mean | Std dev   |\n| ---------------------- | ------------ | ---------- | --------- |\n| DNS                    | 64           | 0.560      | 0.274     |\n| Implicit scale         | 64           | 0.127      | 0.040     |\n| Implicit scale         | 256          | 0.047      | 0.017     |\n| **Multiscale (n_η=5)** | **64**       | **0.028**  | **0.022** |\n\nThis shows that while quadrupling the capacity of the implicit-scale model provides a clear benefit, it is **still twice as inaccurate as our much smaller, more structured multiscale model**. This provides strong evidence that the **explicit separation of scales** is the key driver of performance, not just raw model capacity. This also shows that the coarse mesh compromises the predictive accuracy of DNS substantially, as it is far less accurate than the equivalent-resolution implicit scale model.\n\nAs described in our response to the first comment above, our method belongs to a different modeling paradigm than DMD. DMD would likely perform well on this cylinder flow example, especially with a generous rank of 512. Our work however aims to move beyond flexible models with the capacity to learn from data with multiscale features. Instead, our goal is to enable the data-driven discovery of models that explicitly separate the scales and maintain physical interpretability. DMD is a highly flexible approximator that does not share the same structural priors. Therefore, while a performance benchmark is interesting, our primary comparison is against an implicit-scale model, which isolates the contribution of our core idea: the explicit microscale state.\n\n## Closing Remarks\n\nWe hope these detailed responses and new results have fully addressed your concerns. We are very grateful for your thorough and insightful review, which has helped us significantly improve our manuscript."}]

Author rebuttal:

Meta-review: The paper outlines an approach to learn stochastic multi-scale models by learning a series of coupled SDEs from observational data. Paper has several interesting technical contributions including a principled, probabilistic scheme for scale separation, which might find broad adoption in the field and in the broader ML community. All reviewers and the AC agreed that the paper was well executed and written in its original form and the reviewers concerns primarily focuses on adding additional baseline benchmarks, which were included during the discussion. The main limitation which remained after the discussion is that the experiment only involve small, toy examples which can easily be treated with standard methods, so it is unclear how scalable and generally applicable the presented methodology will be. We encourage the authors to clearly state these limitations in the camera ready version of the manuscript.
Decision: Accept (poster)