Published as a workshop paper at ICLR 2022

N EURAL O PERATOR WITH R EGULARITY S TRUCTURE
FOR M ODELING DYNAMICS D RIVEN BY SPDE S
Peiyan Hu1∗, Qi Meng2†, Bingguang Chen3 , Shiqi Gong3 , Yue Wang2 , Wei Chen3 ,Rongchan Zhu4 , Zhi-Ming Ma3 , Tie-Yan Liu2
1
University of Chinese Academy of Sciences, 2 Microsoft Research Asia
3
Chinese Academy of Sciences, 4 Bielefeld University

A BSTRACTStochastic partial differential equations (SPDEs) are significant tools for modelingdynamics in many areas including atmospheric sciences and physics. Neural Operators, generations of neural networks with capability of learning maps between
infinite-dimensional spaces, are strong tools for solving parametric PDEs. However, they lack the ability to modeling SPDEs which usually have poor regularity
1
due to the driving noise. As the theory of regularity structure has achieved great
successes in analyzing SPDEs and provides the concept model feature vectorsthat well-approximate SPDEs’ solutions, we propose the Neural Operator with
Regularity Structure (NORS) which incorporates the feature vectors for modelingdynamics driven by SPDEs. We conduct experiments on various of SPDEs including the dynamic Φ41 model and the 2d stochastic Navier-Stokes equation, and theresults demonstrate that the NORS is resolution-invariant, efficient, and achieves
one order of magnitude lower error with a modest amount of data.

1

I NTRODUCTION

Stochastic partial differential equations (SPDEs) are significant tools for modeling dynamics inmany areas including atmospheric sciences (Hasselmann, 1976), physics (Uhlenbeck & Ornstein,1930), biology (Wilkinson, 2018), economics (Barone-Adesi & Whaley, 1987), etc. SPDEs generalize partial differential equations via random force terms and they are used to study statisticalmechanics of the dynamics systems. Examples include stochastic Navier-Stokes equations modeling the statistics of turbulent flows (Buckmaster & Vicol, 2019) in atmospheric science and the Φ4model arising in the stochastic quantisation of quantum field theory (Hairer, 2015). Since SPDEs
relate to many scientific open problems, studying the solution of SPDEs from both mathematicalproving and numerical methods is a hot research direction in both math and physics.
Inspired by recent advances in using AI techniques to accelerate scientific computing, we study
using deep learning method for modeling the solution of SPDEs. There have been deep learningmodels arising for modeling dynamics governed by PDE such as Neural Operators (Kovachki et al.,2021), DeepONet (Lu et al., 2019), which model the map between infinite-dimensional functions
and agree with the case of learning solutions of a family of parametric PDEs. However, SPDEs
usually have poor regularity w.r.t the time variable for function-valued noise and singularity w.r.t
space for space-time white noise that these models do not take into consideration. Thus, they cannot
accurately represent the solution of the SPDE.
To deal with the singularity of SPDEs, we incorporate the regularity structure theory (Hairer, 2014)
with neural operator to model the map between (initial condition, driving noise) to the solution. The
key step is to project the driving noise and initial conditions to the model feature vectors in regularity
structure theory which improves the regularity according to the regularity structure theory.
Our Contributions We introduce the Neural Operator with Regularity Structure (NORS) that extends the Neural Operators. This deep-learning-based method has four advantages as follows: (1)∗This work was done when the first author was visiting Microsoft Research Asia.
Corresponding E-mail: meq@microsoft.com.
1
Roughly speaking, regularity describes the smoothness of a function.

†1

Published as a workshop paper at ICLR 2022

The NORS can solve equations whose initial conditions and driving force change and are inputtedsimultaneously, which is beyond the Neural Operators’ capability because of its requirement of regularity. The detailed theory about regularity structure is provided in Section 3. (2) The NORS utilizes
more information from the equations themselves as we take models as our features, which containthe information of the SPDEs’ differential operators and then leads to lower loss. (3) The NORS
inherits Neural Operators’ zero-shot super-resolution, which means it does not need a large amountof data and is mesh-invariant. We test the NORS on the dynamic Φ41 model, reaction-diffusion equation with linear multiplicative noise and the 2d stochastic Navier-Stokes equation. Using the NORS,
both the testing accuracy and sample complexity are enhanced. Specifically, the error is one orderof magnitude lower than other baselines.

2

R ELATED WORK

There have been several popular deep-learning-based methods for modeling the solution of parametric partial differential equations (Lu et al., 2019; Patel et al., 2021; Kovachki et al., 2021; Li et al.,2020b; Bhattacharya et al., 2020; Nelsen & Stuart, 2021; Li et al., 2020a). For example, the Neural
Operator (Kovachki et al., 2021) and Fourier Neural Operator (FNO) (Li et al., 2020b) are representatives which are mesh-independent model, whose architectures approximate Picard iteration forsolving PDEs. Since the solution of SPDEs is determined by both the initial condition and the force,
capturing the structure of the force (e.g., the space-time white noise) is beyond the capability of thesemodels. To handle the case that SPDEs’ solutions depend simultaneously on the initial condition u0
and the force term ξ, Salvi & Lemercier (2021) introduce the neural stochastic partial differential
equation (Neural SPDE), which parameterizes the kernels according to Duhamal’s fix-point formulafor SPDE whose linear differential operators can generate semigroups. In this paper, we adopt another way which first projects the initial condition and force to a set of models. Since the models
incorporate more prior (including the kernel, the initial condition, and the force) of the SPDE, it isexpected to have a better generalization and lower sample complexity.

3

P RELIMINARY

In this section, we introduce background on the regularity structure theory (Hairer, 2014) of SPDEs.
Consider an SPDE on [0, T ] × D with the following form∂t u − Lu = µ(u, ∂1 u, · · ·, ∂d u) + σ(u, ∂1 u, · · · , ∂d u)ξ,u(0, x) = u0 (x),(1)
where x ∈ D ⊂ Rd , t ∈ [0, T ], L is a linear differential operator, ξ is the space-time white noise,u0 : D → R is the initial condition. Under local Lipschitz condition on µ, σ with respect to suitablenorm, this SPDE has a unique mild solution (Hairer, 2014; Salvi & Lemercier, 2021):tL

Z t

ut = e u0 +(t−s)LeZ tµ(us , ∂1 us , · · · , ∂d us )ds +0e(t−s)L σ(us , ∂1 us , · · · , ∂d us )ξds.

(2)

0

Thus in the field of SPDE, the solution is determined by both the initial condition and the force term,i.e., (u0 , ξ) 7→ u. The design of deep learning models such as Neural Operators does not consider
the solution structure of SPDE, therefore, they cannot well approximate the case (u0 , ξ) 7→ u.
It is then natural to utilize the regularity structure theory to help handle the regularity problem. The
concept model in the regularity structure is a collection of model feature vectors, which are multidimensional signals designed to approximate solutions of SPDEs even with low regularity regimes.
The motivation comes from Picard theorem and Taylor expansion. AccordingR to the representationtof the mild solution in Eqn.(2), we define two linear operators I[f ](t) = 0 e(t−s)L f (s)ds andIc [u0 ](t) = etL u0 for any function f defined on [0, T ] × D to Rd . Picard theorem shows that thefollowing recursive sequence approximates the solution u of equation (1) as n → ∞u0t = Ic [u0 ]t ,un+1
= Ic [u0 ]t + I[µ(un ) + σ(un ξ)]t .
t(3)

Using Taylor expansion, we then have the recursive sequence that can approximate u as m, l, n →∞u0,m,l= Ic [u0 ]t ,tun+1,m,l= Ic [u0 ]t +tm
Xµ(k) (0)k=0k!

I[(un,m,l )k ]t +lX
σ (k) (0)
k=0

2

k!

I[(un,m,l )k ξ]t .

(4)

Published as a workshop paper at ICLR 2022

Figure 1: The architecture of our model and the shape of the data of one sample in the 1d case.

Then, the solution of SPDE can be approximated by weighted sum of the features
I[(un,m,l )k ], I[(un,m,l )k ], l = 0, · · · , k; m = 0, · · · , k, where we call n as the height and m, las the width in the approximation. Motivated by this, Chevyrev et al. (2021) develops tool for feature engineering of SPDEs. By the regularity structure theory, the model feature vectors are obtainedby integrals of functionals of u0 and ξ (as I and Ic are convolution operations), whose regularity isproved to be better due to the polishing effect of integrals (Salvi & Lemercier, 2021). To avoid thenumber of model feature vectors grows exponentially, the height of the features is constrained according to the regularity of the SPDE. Please refer the details about the generation of model feature
vectors and its degree constraints in Appendix A.1

4

L EARNING SPDE S OLUTION VIA MODEL FEATURE VECTORS

We move on to introducing the Neural Operator with Regularity Structure (abbrev. NORS). For
given SPDE which has the form in Eqn.(1), our goal is to learn its solution solution uT at given time
point T under initial condition u0 which is assumed to be generated by a parametric distribution.
According to Eqn. (3) and (4), the solution depends continuously on the model feature vectors noton the (u0 , ξ). Therefore, NORS first maps (u0 , ξ) to the model feature vectors, and then we use
Fourier Neural Operator (abbrev. FNO) to learn the continuous map from the model feature vectors
to the solution.
To represent the continuous input functions u0 and ξ, we discretize the space-time domain D×[0, T ]with D ⊂ Rd onto the grid OX1 ×· · ·×OXd ×OT . Then we use the values of the continuous function
on the grid points to represent them. For one sample of u0 and ξ, we first get the model feature
vectors M of (u0 , ξ) according to data and the form of the equation. As M = {fi }i=1,··· ,m are setof continuous functions, we also use its value on discrete grids to represent them. By concatenating
all the model feature vectors fi and the grid OX1 × · · · × OXd , we get the inputs w0 .
Then, w0 is fed into the FNO and the forward process is expressed as
v 0 (x) = Pθin (w0 (x));

v i+1 (x) = Fθi (v i (x));

ûT (x) = Qθout (v K (x))(5)

for any x ∈ D, where θin , θout , θi , i = 0, · · · , K − 1 are learnable weights, Pθin :RX1 ···×Xd ×(m+d) → RX1 ···×Xd ×h be an embedding neural network to project the input to thelatent feature space, Fθi : RX1 ···×Xd ×h → RX1 ···×Xd ×h be a Fourier layer (Li et al., 2020b) whichapproximates the iteration in Picard’s iteration, Qθout : RX1 ···×Xd ×h → RX1 ···×Xd ×d be a embedding layer to project the latent feature to the output. Here, Xi is the number of grids on dimension
xi , h is the number of hidden channels, and m (the number of model feature vectors) and d (thedimension of region D) are defined before. Since only the FNO contains trainable weights, defining loss function between ûT (x) and the groundtruth u(x) can guide the optimization to learn theweights of FNO. A demonstration of our model is shown in Figure 1.
3

Published as a workshop paper at ICLR 2022

Table 1: Dynamic Φ41 model. We consider the l2 error of the baselines and our model(n = 2 or 3)
in two settings with training data size N = 1000 or 10000.
Model
FNO
NSPDE
Ours(n = 3)
Ours(n = 2)

5

N = 1000ξ 7→ u (u0 , ξ) 7→ u
0.032
0.030
0.009
0.012
0.00030.00110.00020.0012

N = 10000
ξ 7→ u (u0 , ξ) 7→ u
0.027
0.024
0.006
0.006
0.00020.00040.00020.0004

E XPERIMENT

We compare the NORS with other baselines on some significant equations. What we care about
includes the l2 error on two settings: in the setting (ξ 7→ u), the noise ξ changes while the initialcondition u0 is fixed; in the setting ((u0 , ξ) 7→ u), both ξ and u0 vary across samples. We notethat in the following equations, we only consider periodic boundary conditions, but Dirichlet or
Neumann boundary conditions can also be easily complemented. To save space, the details about
the construction of the model M of each SPDE are put into Appendix A.1. We use 32 hidden
channels and 4 Fourier layers for our NORS in all experiments. We use the Adam optimizer totrain for 500 epochs with an initial learning rate of 0.001 in the first two experiments and and the
2d stochastic Navier-Stokes equation on 64 × 64 grid, and 0.01 for the 2d stochastic Navier-Stokes
equation on 16 × 16 grid (after grid search) that are halved every 100 epochs. We randomly split thedataset into training and test sets by 5:1. The NORS codes are deposited in GitHub at https://github.com/Peiyannn/Neural-Operator-with-Regularity-Structure.git.
5.1

DYNAMIC Φ41 M ODEL

We first consider the dynamic Φ41 model with the periodic boundary condition. It takes the form∂t u − ∆u = 3u − u3 + ξ, (t, x) ∈ [0, 0.05] × [0, 1](6)u(t, 0) = u(t, 1), (Periodic BC)
u0 (x) = u(0, x) = x(1 − x) + κη(x),
Pk=10
ak
−1
where ξ is the space-time white noise, η(x) = k=−10 1+|k|kπ(x − 0.5)), with ak ∼2 sin(λN (0, 1) with λ = 2, and κ = 0 or 0.1 corresponding to the initial condition is fixed or not.
For this equation, the differential Roperator L is ∆, according to which the operator I and Ic of thetmodel Mn is given by I[f ](t) = 0 e(t−s)∆ f (s)ds and Ic [u0 ](t) = et∆ u0 , where ∆ is the Laplace
operator on D and f : [0, T ] × D → R.
The result is shown in Table 1. We consider two settings, in both of which our architecture outperforms other benchmarks a lot. Even compared with the lowest error of all baselines, our result isabout a tenth of it in the (u0 , ξ) 7→ u setting, while the result of ξ 7→ u setting is even better. We
also note that our model can perform well with few data and low height.
5.2

R EACTION -D IFFUSION E QUATION WITH L INEAR M ULTIPLICATIVE F ORCING

As the dynamic Φ41 model is a parabolic equation with additive forcing, we then consider a parabolic
equation with multiplicative forcing, which is given by∂t u − ∆u = 3u − u3 + uξ, (t, x) ∈ [0, 0.05] × [0, 1]u(t, 0) = u(t, 1), (Periodic BC)
u0 (x) = u(0, x) = x(1 − x) + κη(x),(7)

where ξ is the space-time white noise, and η(x) is the same as the Φ41 model. As the form of theoperator I and Ic of this equation is the same as the Φ41 model, the model Mn can be constructedsimilarly. Please check the details in Appendix A.1.
4

Published as a workshop paper at ICLR 2022

Table 2: Reaction-Diffusion equation with linear multiplicative forcing. We compare the l2 errorof the baselines and our model(n = 2 or 3) with training data size N = 1000 or 10000.
Model
FNO
NSPDE
Ours(n = 3)
Ours(n = 2)

N = 1000ξ 7→ u (u0 , ξ) 7→ u
0.00360.00630.00160.00620.00060.00050.00060.0006N = 10000
ξ 7→ u (u0 , ξ) 7→ u
0.00350.00370.00120.00260.00050.00030.00050.0003

Table 3: 2d stochastic Navier-Stokes equation. We compare the l2 error of the baselines and ourmodel(n = 2 or 3) in two settings with 1000 training samples. While solving the equation on 64×64
grid, we train the model on 64 × 64 and 16 × 16 grid respectively.
Model
NSPDE
Ours(n = 3)
Ours(n = 2)

64 × 64 grid
ξ 7→ ω (ω0 , ξ) 7→ ω
0.039
0.031
0.00170.00290.00180.002816 × 16 grid
ξ 7→ ω (ω0 , ξ) 7→ ω
0.074
0.063
0.00200.00340.00220.0030

The results in Table 2 show that our model has one order of magnitude lower error in both of thetwo settings. The experiments on the two equations clearly show that the effectiveness of the modelfeature vectors and the worse generalization of FNO on SPDEs.
5.3

2 D S TOCHASTIC NAVIER -S TOKES E QUATIONAs both NSPDE and our model claim mesh-invariance, we evaluate both the l2 error and the meshinvariance property of NSPDE and our model on a 2d Navier-Stokes equation for an incompressible
flow:∂t w − ν∆w = −u · ∇w + f + σξ, (t, x) ∈ [0, 0.05] × [0, 1]2(8)
ω(0, x) = ω0 (x)(9)

where u is the velocity field, ω = ∇ × u is the vorticity, ω0 is the initial vorticity, f is the deterministic force defined as in (Li et al., 2020b), ξ is the random force rescaled by σ = 0.05 defined as in(Salvi & Lemercier, 2021), and the viscosity parameter ν = 10−4 .
Our target is to model the vorticity ω, which is harder to learn compared with the velocity
u. According
to the form of the equation, the operation I and Ic in model Mn is defined as
R t (t−s)ν∆I[f ](t) = 0 e
f (s)ds and Ic [ω0 ](t) = etν∆ ω0 , where ∆ is the Laplace operator definedon the 2d space. While solving the 2d Navier-Stokes equation on 64 × 64 grid, we train the modelon 64 × 64 and 16 × 16 grid respectively to test the mesh-invariant property of our model. We follow
the data generation process in (Salvi & Lemercier, 2021) to generate the ground-truth for training.
As shown in Table 3, the error of our model is one order of magnitude lower in the ξ 7→ ω and(ω0 , ξ) 7→ ω settings. Besides, we solve the equation on 64 × 64 grid, then train on the 64 × 64 gridand 16 × 16 grid. From the results, we verify the resolution-invariance of our model.

6

C ONCLUSION AND F UTURE W ORK

In this work, we introduce NORS as a strong SPDE-solving tool with the zero-shot super-resolution
property. By incorporating the regularity structure, the NORS absorbs both the advantages of Neural
Operators and regularity structure, and makes up for the shortcomings. Not only can the NORS
learn solution operators (u0 , ξ) 7→ u of SPDEs, but also has a much lower error. In the future, as
the NORS requires that the differential operator L is already known, we can extend this method byparameterizing the kernel, which will be able to handle the inverse problem that some part of theequations is unknown.
5

Published as a workshop paper at ICLR 2022

R EFERENCES
Giovanni Barone-Adesi and Robert E Whaley. Efficient Analytic Approximation of American Option Values. Journal of Finance, 42(2):301–320, June 1987. URL https://ideas.repec.
org/a/bla/jfinan/v42y1987i2p301-20.html.
Kaushik Bhattacharya, Bamdad Hosseini, Nikola B Kovachki, and Andrew M Stuart. Model reduction and neural networks for parametric pdes. arXiv preprint arXiv:2005.03180, 2020.
Tristan Buckmaster and Vlad Vicol. Convex integration and phenomenologies in turbulence. 2019.
Ilya Chevyrev, Andris Gerasimovics, and Hendrik Weber. Feature engineering with regularity structures. 2021.
M. Hairer. A theory of regularity structures. Invent. Math, 198(2):269–504, 2014.
Martin Hairer. Regularity structures and the dynamical ϕ43 model. 2015.
K. Hasselmann. Stochastic climate models part i. theory. Tellus, 28(6):473–485, 1976. doi: 10.3402/
tellusa.v28i6.11316. URL https://doi.org/10.3402/tellusa.v28i6.11316.
Nikola B. Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya,
Andrew M. Stuart, and Anima Anandkumar. Neural operator: Learning maps between function
spaces. CoRR, abs/2108.08481, 2021. URL https://arxiv.org/abs/2108.08481.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Andrew Stuart, Kaushik Bhattacharya, and Anima Anandkumar. Multipole graph neural operator for parametric partial differential equations. Advances in Neural Information Processing Systems, 33:6755–6766, 2020a.
Zongyi Li, Nikola B. Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya,
Andrew M. Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. CoRR, abs/2010.08895, 2020b. URL https://arxiv.org/abs/2010.
08895.
Lu Lu, Pengzhan Jin, and George Em Karniadakis. Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators. arXiv
preprint arXiv:1910.03193, 2019.
Nicholas H Nelsen and Andrew M Stuart. The random feature model for input-output maps between
banach spaces. SIAM Journal on Scientific Computing, 43(5):A3212–A3243, 2021.
Ravi G Patel, Nathaniel A Trask, Mitchell A Wood, and Eric C Cyr. A physics-informed operator regression framework for extracting data-driven continuum models. Computer Methods in Applied
Mechanics and Engineering, 373:113500, 2021.
Cristopher Salvi and Maud Lemercier. Neural stochastic partial differential equations. CoRR,
abs/2110.10249, 2021. URL https://arxiv.org/abs/2110.10249.
G. E. Uhlenbeck and L. S. Ornstein. On the theory of the brownian motion. Phys. Rev., 36:823–841,
Sep 1930. doi: 10.1103/PhysRev.36.823. URL https://link.aps.org/doi/10.1103/
PhysRev.36.823.
D.J. Wilkinson. Stochastic Modelling for Systems Biology, Third Edition, volume 1. Taylor &Francis Ltd, 2018.

6