Riesz Neural Operator for Solving Partial Differential Equations
Authors: shouyiliu
Xiaokang Yang
Yuntian Chen
Venue: ICLR 2026 Poster (2026)
Keywords: Neural Operator
Local derivative
Taylor expansion
Nonlinear
PDEs
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)
TL;DR: Riesz Neural Operator blends global spectra and local derivatives via the Riesz transform, capturing non-stationarity and outperforming prior PDE models.

Abstract: Local non-stationarity is pivotal to solving partial differential equations (PDEs). However, in operator learning, the spatially local information inherent in the data is often overlooked. Even when explicitly modeled, it is usually collapsed into local superpositions within the model architecture, preventing full exploitation of local features in physical phenomena. To address this limitation, our paper proposes a novel Riesz Neural Operator (RNO) based on the spectral derivative representation. Since PDEs are fundamentally governed by local derivatives, RNO leverages the Riesz transform, a natural spectral representation of derivatives, to mix global spectral information with local directional variations. This approach allows the RNO to outperform existing operators in complex scenarios that require sensitivity to local detail. Our design bridges the gap between physical interpretability and local dynamics. Experimental results demonstrate that the RNO consistently achieves superior prediction accuracy and generalization performance compared to existing approaches across various benchmark PDE  problems and complex real-world datasets, presenting superior non-linear reconstruction capability in model analysis.

Published as a conference paper at ICLR 2026
RIESZ NEURAL OPERATOR FOR SOLVING PARTIAL DIFFERENTIAL EQUATIONS
Shouyi Li Xiaokang Yang', Yuntian Chen
Shanghai Jiao Tong University
Eastern Institute of Technology
shouyiliues jtu..du.ccm cn, xxyangesjtu.edu.. edu.uunnddccn ychenêeias. acc.cacccoo
ABSTRACT
Local non-stationarity is pivotal to solving partial differential equations (PDEs). However, in operator learning, the spatially local information inherent in the data is often overlooked. Even when explicitly modeled, it is usually collapsed into local superpositions within the model architecture, preventing full exploitation of local features in physical phenomena. To address this limitation, our paper proposes a novel Riesz Neural Operator (RNO) based on the spectral derivative representation. Since PDEs are fundamentally governed by local derivatives, RNO leverages the Riesz transform, a natural spectral representation of derivatives, to mix global spectral information with local directional variations. This approach allows the RNO to outperform existing operators complex scenarios that require sensitivity to local detail. Our design bridges the gap between physical interpretability and local dynamics. Experimental results demonstrate that the RNO consistently achieves superior prediction accuracy and generalization performance compared to existing approaches across various benchmark PDE problems and complex real- world datasets, presenting superior non-linear reconstruction capability in model analysis.
INTRODUCTION
Accurately simulating physical systems governed by partial differential equations (PDEs) remains foundational across fluid dynamics (Herde et al., 2024), materials science (Li et al., 2024), and climate modeling (Kurth et al., 2023). Neural operators have recently emerged as powerful surrogates for PDE solution operators, offering significant reductions in compute time and broad adaptability. However, it is a complex task to capture and predict the non-linear and local changes due to physical information missing, transform without local emphasis, etc. This work leverages local spectral derivatives to improve the precision PDE solutions.
Existing operator-leamning frameworks, including DeepONet (Lu al., 2021; Gu al., 2025), the Fourier Neural Operator (FNO) (Li et al., 2020a), and the Laplace Neural Operator (LNO) (Cao al., 2024a), leverage specific transforms model PDE dynamics. Suppose we analogize the process of function learning to the performance of musical score. In DeepONet, each partitioned ensemble (trunk function) is assigned fixed section, which plays as required by the score. This design elegantly captures the global melody of conctiionto-functtin mappings. However, it lacks control over melodic variations and requires an excessively large ensemble for complex scores. FNO, by contrast, assigns finer-grained roles with varying intensities to all performers but struggles to capture melodic diversity and local harmonies. LNO augments FNO with a smoother, yielding greater stability the cost of detailed melodies. These intriguing comparisons Figure ((a). Transformer based operator learners similarly apply global attention to imputs. As result, physical locality remains underexploited in operator learning. leading to two major shorteomings: (i) Non-stationary information is largely encoded in local spatial features (Pacejka & Besselink, 1997). Since the temporal dimension is typically embedded implicitly during training, such non- stationarity is not explicitly represented (Barberi & Kruse, 2023). We therefore hypothesize that once spatial locality is captured adequately, non-stationary temporal signals can be learned naturally.
presented in
"oorreponning author.

---

Published as a conference paper at ICLR 2026
smoothor
dioctional boods local weight orientation
JC mi
2 trunk
D
baen
IT
brunch
do
DeepONet
FNO
LNO
RNO
z-rlesz
df((t)
Direction Mixer
dz
Mixer
z-kernel
df(t)
dr.
f(t)
dovecctoo desinttve
x- x-kernel
x-riesz
df((t)
space
dy
(d)
y y-kernel
Pounter Inaiforr Fint--rder derivaivv
Abrupt change along time
y-riesz
Ress transform
No abrupt change information
Convolution
frequency
Eime
Figure 1: (a) Taking musical score processing as an example, we compare DeepONet, FNO, LNO, and RNO. (b) In the Riesz space, the Riesz transform yields direction-selective kernels (e.g.. an x x-kernel for the x-direction). Each kernel computes the spectral derivative along its associated direction. (c) The RNO pipeline is formalized mathematically.
(d) In the Fourier spectrum (left), no localized change appears in the figh-freeuency band where a discontinuity is expected, whereas the first-order derivative (right) cleanly captures the emerging jump over time, demonstrating its local sensitivity.
To substantiate this hypothesis, we seek reconstruct spatial locality more accurately and provide the proof in the A.3. (ii) Discarding spatial locality forfeits fine-scale details, especially the small perturbations typical PDE-governed phenomena(Steans. 1985; Xiong al., 2025). From a spectral perspective, these correspond to high-frequency content; neglecting locality therefore suppresses high-frequency representations and weakens the modeling of strong nonlinearities as local detail naturally corresponds to high-freuenny content in images( maggss(Wang et al., 2020). For further details on local necessity comparison, please refer to the test A.1. In the experiments, phase plays crucial role in the local spectral orientation. Although prior work addresses complex physical systems (Li et al., 2021; Wang et al., 2023), emphasizes local structure (Raonic et al., 2023; Wang et al., 2022) and uses derivative optimize the loss(Cheng al., 2025; Qiu et al., 2024), to our knowledge no aperator-learning method that optimizes non-linear dynamics by mixing local derivatives. Concurrently,
increasing share of advances stems from incorporating established architectures (e.g., Transformer( Wu et al., 2023), mamba(T'wari et al., 2025a), and graph network al., 2025), whereas recent attempts modify the neural operator itself have produced limited gains. In these solver variants that introduce external structures (Liu et al., 2025; Kalimuthu et al., 2025), a key property of classical neural operators is disrupted: their ability to learn the intrinsic continuity of the underlying physical fields (Li et al., 2021). This disruption, in turn, limits the continuity of the model with respect to its inputs. Moreover, introducing external structures can partially compromise the efficiency advantages of classical neural operators.
Inspired by the Taylor expansion- where higher-order terms encode subtle, local variations (Kanwal & Liu, 1989; Bruhn et al., 2005), therefore derivatives are naturally suited to capture locality Figure I(c). Yet operator-learning pipelines dominated by global transforms (e.g., Fourier, Laplace) often underrepresent such information, and the energy from abrupt changes is difficult to capture (Alley et al., 2003). To address these bottlenecks, we propose the Riesz Neural Operator (RNO), which improves neural operators' ability to model complex dynamic variations. The Riesz transform in our RNO acts as the conductor's 'ss gestures, injecting directional derivatives that govern local dynamics, while the direction mixer serves as orchestration, harmonizing different directional
2

---

Published as a conference paper at ICLR 2026
components. Analytically, Taylor expansion separates behavior: the DC (direct component) tracks slow trends, whereas abrupt or small-scale variations appear as higher- terms (Hodaei et al., 2017; Bruhn et al., 2005). Building this perspective, we apply the Riesz transform to extract first-order derivative information (direction and rate of change) (Unser & Van De Ville, 2009; Wadhwa et al., 2014; Tamada & Igarashi, 2017), and reproject these derivatives onto the original spectrum to emphasize salient components. The result is not only faithful reproduction of the global score but also nuanced capture of fine-scale, anisotropic variations much like the layered expressiveness of super symphony. The main contributions of this work are as follows:
Motivated by the Taylor expansion, we encode fine-scale structure using local spectral deriva- tives. Integrating the Riesz transform sharpens sensitivity to dynamics while preserving physical interpretability, improving expressivity for non-linear phenomena.
We are the first to explicitly embed the Riesz transform within operator learning, providing direction-aware operator mixer to unify global basis information and local directional derivatives.
Our method supports a principled composition for reconstructing physical dynamics.
RNO achieves state-of-ther -hetttr-rtt performance in extensive benchmarks, covering various complex Navier-Stokes PDEs and the real-world climate ERA5 dataset. Experiments indicate that its low-level feature mixer affords a distinct advantage for modeling strong nonlinearities.
PRELIMINARIES
Neural Operator. Neural operator learning represents a novel paradigm intersecting scientific computing and machine learning, aimed at learning solution operators for PDEs. Owing to their full differentiability, neural operators permit direct parameter optimization, making them ideal for inverse design and a broad spectrum of inverse problem settings. Initially, neural operators originated from DeepONet(LL et al., ). This paradigm later evolved into the classical FNO(Li et al., 2020a;b), which also performs function-tor mnnctiontooffncttion to-funciion mappings in Fourier space and has spawned several FNO variants(Tran et al., 2021; Li et al., 2023). Recently, study employed the Laplace transform as an alternative to the Fourier transform(Cao et al., 2024a) to address non-periodic characteristics data, significantly enhancing the interpretability of neural operators. This highlights the significance of handling non-stationary data(Tiwari et al., 2025b). Additionally, other studies have concentrated on the capability of models to address complex phenomena(Li al., 2021; Lanthaler et al., 2023; Xiao et al., 2024), strengthening the neural operators' ability to handle nonlinearity.
Meanwhile, due the complex variability of real-world data and the higher standards required for cap- turing detailed variations, local properties within neural operators have also attracted attention (Raonic al., 2023; Tripura& & Chakraborty, 2022). Among these studies, the work on neural operators with localized integral and differential kernels emphasizes the necessity of local and differential properties. However, their differential kernels and local kernels are designed separately and introduced into the model as independent modules. Our approach leverages domain knowledge. integrates lask-specific features with directional information, and thus possesses the capability handle local signals, all while maintaining interpretability.
Riesz Methods. Spatial derivatives are widely used in deep learning. In Physics-Informed Neural Networks (PINNs) (Raissi et al., 2019), temporal/spatial derivatives construct PDE residuals, ensur- ing predictions both fit observations and satisfy governing laws. Derivative signals have likewise been incorporated into LSTM variants to better track state changes (DiPietro & Hager, 2020). By contrast, explicit frequency-donee -dommin operations are less common: most prior work analyzes networks' spectral behavior during training (Xu et al., 2019) rather than designing models centered on spectral directionality.
The Riesz transform maps spatial derivatives into the spectral domain while blending integral and differential effects: the integral term preserves neural operator advantages, and the differential term injects local variation (Tamada & Igarashi, 2017; Wadhwa et al., 2014). Here we use the vector Riesz transform from multidimensional signal processing, which generalizes the 1D Hilbert transform to higher dimensions; the relationship is derived in A.7. Formally, the transform is convolution with Calderón calderón-Zygmund kernel (de Francia et al., 1 1986), and its kernel-based processing is illustrated in Figure (bb). A defining property is strong directional selectivity: the transform decomposes signals
3

---

Published as a conference paper at ICLR 2026
along specific spatial directions, enabling analysis of anisotropic structure (e.g., edges and textures) (Unser et al., 2009) and isolation of instantaneous phase. These properties make the Riesz transform well suited to data with complex, non-stationary variability.
METHODS
3.1 PROBLEM DEFINITION
In physics, is widely recognized that PDEs characterize a function by constraining its local derivatives(Wu al., 2023). That is, the evolution of a physical system is described through relations among its partial derivatives:
((x,,(x), u(z), az u(x)) = 0,
where 0z u(x:) denotes the local derivative of with respect to 24. Physically, the presence of many local differential operators means that PDEs contain wide range of higher-order dynamical quantities, which together govern the rich variety of behaviors observed in these systems. More broadly, accurately capturing small-scale variations is beneficial for physical signals. Since deep learning aims to exploit patterns, inputs with natural and approximately orthogonal dynamic components enhance the model's expressive capacity. In addition, many physical governing equations mirror optical-flow formulations that decompose fields into smooth backgrounds and dynamic perturbations; see A.9 for details. Nevertheless, the implications this locality have been insufficiently exploited. We therefore seek principled representation of this local derivative form. As the natural spectral representation of derivatives, the Riesz transform provides more than an engineering convenience: embedding it into neural operators offers mathematically and physically consistent means of unifying global spectral efficiency with local dynamics. For a further relationship discussion of local derivative and PDEs with much higher degree of mathematical rigor, refer the readers to A.1.
3.2 DERIVATIVE INFORMATION IN PHYSICAL FIELDS
Local dynamics via Taylor expansion. Adopting the Eulerian viewpoint commonly used in optical flow analysis(Bruhn et al., 2005), we describe temporal variation of a one dimensional signal by the first-order Taylor series
df(x) dr
1xx Re f(x) + 7(t)
where [(x, () f(x + 7()) and (L) denotes the instantaneous displacement at time (>> 0 as shown Figure The term r(t) (z) therefore encodes motion magnitude implicitly inside the spatial derivative. Neglecting this factor, as in most vanilla neural operator blocks, effectively acts as low-pass filter and can compromise fidelity in highly dynamic regimes.
Directional derivative. Beyond magnitude, first-order derivatives carry orientation information as described in A.1. The directional derivative of u(x, t) along a unit vector n (f,,...,77, naa] is
d du
Fu-n =
T4.
Or:
Consequently, model predictions must depend not only on scalar changes but also on the direction in which those changes occur.
Spectral derivative. gives
Taking the temporal derivative of f(t) and transferring it to the Riesz space
df(()) 7()
(())k(((())
(4)
dt
where & is the wave component and Ø denotes Riesz transform. In this representation, spatial derivative manifests as a 90º phase shift accompanied by an amplitude scaling proportional to In 2D, L- and L, act projectors onto directions defined by k's angles with the axes, allowing efficient computation of spatial (and higher-orden derivatives. This relation enables the Riesz transform to compute spatial derivatives efficiently and incorporate higher-order quantities.
4

---

Published as a conference paper at ICLR 2026
Dynamics in Riesz space. Let h(q. t) be scalar field with spatial coordinates q [41 ......... R" and time t Rt. Its complex-valeee valued spectrum is defined by the exponential kernel
h(k,t) = [ h(q,
h(a t)
kg dq.
k == [k]] E C",
kn]"
(5)
where = 1. The magnitude ||k
12)1/2 controls spatial scale, while the vector
direction of k encodes the dominant orientation of the corresponding spectral component. According the definition of the directional derivative, we apply the n-D Riesz transform & direction yields
h; [
Vh(q.tt) (k)
Eihe
h(q, k)e-gkadg.
kadg.
(6)
==1
showing that differentiation becomes a multiplicative modulation by the directional factor j hıki / [||]
This highlights& key property of the Riesz framework: normalising by ||k|] preserves energy while amplifying orientation selectivity, enabling precise estimation of directional rate of change without additional scaling heuristics.
3.3
RIESZ NEURAL OPERATOR
Operator learning. To enhance representational capacity, the input is first projected into higher- dimensional space by a shallow connection layer, yielding (t). The transformed q(t) then passes through an integral- -ernell layer and linear layer, followed by activation, produce the high- dimensional output representation,
u(I) (Wa(l) + (c+ WED
Here, D a bounded open set; W a linear transformation; and (c q)(() = Sp ctt, {(s) ds denotes the integral operator with kernel By Green's 's function theory (Li et al., 2020b), the kernel ctt, s) can approximate solution operators for PDEs. As shown in Figure 9, rather than relying on single global transform, couple local and global components to more fully map inputs to outputs. Indeed, RNO preserves the compositional structure characteristic of neural operators. Its overall procedure can be decomposed into three modules,
CaordToRion I EM, OR F @RieneTTCCordd
(8)
where denotes operator composition and F is mapping between spaces. Specifically, BowardToRRum X + Rx lifts the input to the Riesz domain, and Frumeeoooordd Ry - y
reconstructs the final output. This modular factorization is generic for neural operators. In the Riesz domain, F -M #R learns and mixes non-linear, direction- wise mappings, where M is the number of directional components. In practice, components are combined via analytic signals (1D) (Marple, 1999) or monogenic signals (multi-D) (Unser al., 2009), providing a natural encoding of signal characteristics. Directional weights emphasize grientation-depender variation. The enhanced Riesz representations are then mapped back by F FR to yield spatial output with enhanced directional features.
2i=1
Riesz conductor. Analogous to convolution performed in a generic frequency domain, the kernel (c* q) is realised by multiplication in Riesz space. For m-dimensional data the i-th directional Riesz transform of q is
j ki - |k|
(.
e-jka dma
The spectral factor j k behaves scale-normalised first-order derivative, thus emphasising local variations while suppressing magnitude. Using the forward (x) and inverse (08 Riesz transforms, the integral operator
(c (t) = -e- ( (t))
VtED, D,
(10)
where "." denotes element- wise complex multiplication. The model overview of RNO is shown in Figure 2. In the Riesz conductor, the directional numbers not chosen randomly but instead based on the data's 's dimensionality. This approach aligns with the inherent characteristics of physical field data. For example, in avo-dimensional data, two orthogonal directions are used, maximizing
5

---

Published as a conference paper at ICLR 2026
data utilization and optimizing results without introducing redundancy, the proof can be found in A.5. Each direction also multiplied by scaling factor Ç. which controls the contribution of each direction. To account for the possibility that incorporating heterogeneous perturbations may introduce data artifacts, we establish an upper bound on Ç.
Eumi(t)
π
(11)
where Wi denotes the frequency along the i-th direction. The upper bound on follows from the aggregate total variation across directions; see A.6 for the full derivation. The resulting spectral enhancement is illustrated in Figure 2: the Riesz integral kernel exhibits increased energy in high- frequency detail bands, mitigating the high-frequency issues reported during training (Xu et al., 2019).
P Lifting
Riesz Layer
Riesz Layer
Riesz Layer I.
Projection
Original Spectall Analsts
Local Bass vvaredarm
R
Invens Ramz transform
c2
Devvooe
ptDe
Riest Spraard Aaalsss
R cc
weight
Ireer transtorm
riar
Rieo
activatio funtiin
The full RNO preserves the canonical neural-operator scaffold. Bottom left: Riesz layer 1. The input ((1) is processed in parallel by (i) direction local transforms, weighted by Ç, to extract orthogonal local features R; (ii) global transform that stabilizes the representation; (iii) linear map W. The mixed output passes through an inverse Riesz transform to reinforce spatial locality, followed by an activation to produce 7(2). Bottom right: Spectral comparison before and after the Riesz integral operator. Energy in the high-frequency band (red circle denotes spectral mean is evaluated under a fixed reference: larger magnitudes appear outward expansion, smaller inward contraction.
Figure 2: Model Overview of RNO. Top:
Local orchestration. introduce a direction-aware mixer in the Riesz domain that orchestrates global-local features across spectral scales. S @coord aggregates global Riesz-spectrum statistics, while directional branches capture local anisotropic variation; their fusion yields a monogenic-se representation that preserves coherence and fine directional detail, beyond magnitude-only operators. Details of the construction are provided in A.10. Data layout and formulas appear Figure 1. In 1-D, Cartesian channels not meaningful, wavelength conditions mixing to fuse forward/backward components (implementation details are in A.1
How derivatives help. Derivatives enhance the operator's physical representation by resolving local variations that global features cannot capture. As in a Taylor expansion, where derivative terms refine the leading field u(x,t), they provide principled way to encode fine-scale spatial and temporal changes. This yields a natural split: low-order global components model smooth, low-frequency structure, while orientation-awar derivatives concentrate energy at high wavenumbers and in anisotropic regions. In RNO, Riesz-based directional derivatives encode oriented changes in latent spectral space, and a small set mixing weights controls their contribution. During training, gradients emphasize directions that most reduce the loss, the derivative channels supply missing detail while the global modes carry the coarse structure, leading to sharper and more physically meaningful predictions. For further discussion, see Sec. A.9.

---

Published as a conference paper at ICLR 2026
Table 1: Test relative lg error (L) on five PDE benchmarks. Lower better. The best (bold) and second-best (underline) results are highlighted. Promotions indicate relative gain of RNO over the second-bes best.
Relative f (4)
Model
Duffing
Beam
Diffusion
Reaction- Diff.
Brusselator
U-NET(2015)
0.3835
0.0522
0.0212
0.1049
0,2557
GALERKIN(202))
0.4239
0,0971
0.0319
0.1983
0,1883
DEEPONET(202I)
0.6289
0,9790
0.1356
0.7969
0.1910
FNO(2020a)
0.4536
0,0821
0.0229
0.1214
0.1827
WNO(2022)
0,2532
0.0571
0.0179
0.1713
0.1789
F-FNO
0.2782
0.0515
0.0143
0.1003
0.1719
LSM(2023)
0.1699
0,0456
0.0145
0.0996
0,1628
LNO(2024a)
0.3325
0,0452
0.0081
0.1355
0,1858
ONO(2024)
0.3491
0.0519
0.0137
0.0989
0.1545
RNO
0.1663
0.0219
0.0079
0.0899
0.1317
Promotion
2.1%
51.6%
2.5%
10.3%
14.7%
Top 4 ranking methods of PDE benchmarks: RNO (ours), LNO (2024), ONO (2024), LSM (2023).
EXPERIMENTS
Benchmark. RNO is evaluated on five diverse PDE benchmarks (Cao et al., 2024a)-Duffing Beam, Diffusion, Caaction-Diffusiinnn and Brusselator(spanning 1-D 3-D, details in A.18). To assess performance on complex and abrupt dynamics, we further test 2-D Navier Stokes equations at Reynolds numbers 40, 500, and 5000 (Li et al., 2021), where higher Reynolds numbers induce richer rotational structure, probing sensitivity to directional capture. We then compare RNO with classical baselines the real-world ERA5 dataset (Hersbach al., 2020). These tasks constitute canonical PDE application scenarios and reflect practical challenges in solving such equations.
Baseline. We compare RNO with seven well-established neural network models: a classic model, operator learning models, Galerkin Trans-
U-Net(Ronneberger al., 2015); two transformer- former(Cao, 2021) and LSM(Wu et al., 2023); and six classical and powerful neural operator models: DeepOnet(Lu al., 2021), FNO(Li al., 2020a), WNO(Tripura & Chakraborty, 2022), F-FNO(Tran et al., 2023), ONO(Xiao et al., 2024) and LNO(Cao et al., 2024a). Each model excels in its respective domain. The comparison enables comprehensive evaluation of RNO's effectiveness.
Implementation details. To ensure fairness, we evaluation metrics and training epochs across methods and use Adam (Kingma & Ba, 2014) optimizer. Metrics: relative l error for the general PDE and ERA5 benchmarks, and mean-squared error (MSE) for Navier Stokes equations. For more details, please refer to A.16.
4.1
MAIN RESULTS
Performance across diverse PDEs. To assess generality, we evaluate RNO on five canonical PDE systems: Duffing, Beam, Diffusion, Reaction diffusion, and Brusselator, spanning mechanical, diffusive, and reactive dynamics. RNO attains the lowest prediction enror on all benchmarks (Table I); on Beam it reduces error by 51.6% relative to the second-best method, This substantial margin underscores RNO's ability to capture structural response behavior with high precision. Similarly, on the Duffing system, classical non-linear oscillator known for its chaotic dynamics, RNO achieves 2. 2.1% improvement, indicating its capacity to model sensitive dependence on initial conditions and inherent nonlinearity effectively.
On more still outperforms prior approaches. While these systems generally more stable and casier to approximate, RNO delivers 10.3% and 2.5% reductions error, respectively, demonstrating its strength not only in complex non-linear settings but also in capturing fine-scale spatial dynamics
tasks, such as Diffusion and Reaction-DiffussionRRee InffusionnRRec-Difii RNO
7

---

Published as a conference paper at ICLR 2026
Ground Truth
FNO
DeepONet
Figure 3: Comparison of high-frequeng components with the Ground Truth.
Table 2: Mean-squared error (4) results of RNO, LNO, WNO, F-FNO, FNO, and DecepONet models on NS datasets.
Datasets
DeepONet
FNO
F-FNO
WNO
LNO
RNO
Re 40
0.0280
0.0078
0.0065
0.0077
0.0060
0.0049
Re 500
3,4082
1.4251
1.2312
0.9388
1,2117
0.4861
Re 5000
6.2721
2.9314
2.8912
2.5914
2.3139
0.9121
in smoother regimes. Finally, for the Brusselator, prototypical model for nonlinear chemical oscillations and pattern formation, RNO achieves 14.7% reduction in relative error. This result is particularly significant as it highlights RNO's capability to handle stiff systems with intricate feedback mechanisms and spatiotemporal instabilities.
Performance on complex tasks. To evaluate RNO under highly non-linear dynamics, we train on two-dimensional incompressible Navier Stokes snapshots three Reynolds numbers Re {40, 500, 5000). As Re increases, the inertial range broadens and turbulent interactions intensify, rendering prediction progressively harder. Mean-squared error (MSE) for RNO, LNO, WNO, FNO, FNO and DeepONet is reported in Table 2. RNO achieves the lowest MSE across all Re; at Re = 5000 it attains MSE = 0.9121, substantially outperforming DeepONet, FNO, F-FNO, WNO and LNO. These results indicate that RNO's recurrent updates enhance stability flow complexity grows. To assess fine-scale reconstruction, we isolate the high-frequency band of predicted velocity fields (Figure 3). RNO most faithfully recovers the nigh-mavenumber energy spectrum; LNO underestimates the spectral tail, whereas FNO exhibits ringing artefacts. DeepONet, however, produces of noise in the high-frequency components. These findings confirm that RNO captures subtle non-linear fluctuations more accurately than the alternatives.
We further evaluate on ERA5 using hourly 850-hPa geopotential-height fields (2012-2022 for training; 2023 for testing). Further details of dataset can be found 18...7. As shown in Table 3, RNO significantly outper- forms FNO. LNO, WNO, and Deep- ONet. RNO recovers both the inner
Table 3: Relative & error (L) results of RNO, LNO, WNO, FNO, and DeeepONet ERA5.
Dataset DeepONet
RNO 0.0093 0.0085 0.0062 0.0022
FNO WNO
LNO
ERA5 0.0912
low-frequency ring and the outer high-frequency core with smaller bias, as confirmed by PSD (Power spectral density) analysis (Figure 4). Across both synthetic turbulence and real-world weather, RNO preserves accuracy as dynamics intensify and faithfully reconstructs high-frequency content. also transfers to operational-scale data without tuning, demonstrating its suitability for modeling and forecasting strongly non-linear spatiotemporal phenomena.
4.2 MODEL ANALYSIS
Ablation study. To assess each module's contribution to the overall architecture, we evaluate four designs two benchmarks (Duffing and Reaction-Diffusion) global-only (g) without mixer (o), local-only (f) without mixer, global+local (g+/) without mixer, and global+local (g+I) with mixer (w). From Table 4 we draw two clear observations: (i) on both datasets, either the global the local module in isolation underperforms their combination; and (ii) the proposed mixer further facilitates effective information integration within the network, yielding additional gains.
8

---

Published as a conference paper at ICLR 2026
Ground Truth
Ours
LNO
WNO
ENO
DeepONNt
Hege
Hope
bojoys Mapooet
W4-11
width
widh
Widi
Figure 4: ERA5 benchmark at t = 5: Comparisons of 2-D PSD (top), spatial error (bottom right) and PSD along the x-axis for each method (bottom left).
Influence of orthogonalization. In practice, we set the number of directions equal to the data dimensionality exploit the input's natural orthogonality. To isolate its effect, vary the direction set and evaluate on three benchmarks (Beam, Reaction- Difffusion, and the Brusselator). We sweep the number of directions from 1 to 5, defining each by uniform partition of the angular domain. We compare RNO with LNO, FNO, and DeepONet to quantify the role of orthogonal representations; in 2-D, two directions align with the native coordinate axes and are exactly orthogonal. Analogously in 3-D, the three canonical axes are orthogonal.
Beam
Reac-Diff.
A
Brusselator
Direction number
-19.4%
-6.9%
df
4
7
+29.49
A
+35.5%
2
3
6
5
8
-62%
-31%
0% 31%
62%
RNO-2
RNO4
RNO-
RNO-1
RNO-1
LNO
FNO
DespONa
Figure 5: Comparison of relative &g errors (left) and rank (right) for RNO, LNO, FNO, and DeepONet Reaction-Diffussonn
on the Beam, Reaction Diffusion, and Brusselator benchmarks as the direction number varies from 1 to 5. RNO-n indicates that the RNO uses n directions.
As shown in Figure 5, RNO achieves its best performance with exactly two directions and remains competitive(typiclly second-best) when the number of directions is an inte- ger multiple of the data dimension (e.g.. 4). With too few directions (e.g., 1), we observe information loss and occasional instability; with too many (e.g.. 3 or 5), performance does not improve because the additional, non-orthogonal di- rections are redundant. Even under suboptimal counts, RNO outperforms LNO, FNO, and DeepONet on most metrics. This highlights that orthogonality is not just heuristic choice, but principled design element in neural operators.
Table 4: Relative la error (4) results of ablation study.
Datasets
g(o)
I(o)
g+l (o)
g+/ (w)
Duffing
0.2098
0,2262
0.1801
0.1663
Reac-Diff..
0.0953
0.1189
0.0903
0.0899
Influence of activation function. In prior experiments, all methods employed the same sin activa- tion to ensure parity. Prior work identifies activation functions as the primary source of nonlinearity in neural operators (Tran et al., 2023; Lu et al., 2021; Li al., 2020b). To assess RNO's intrinsic capacity for non-linear phenomena, we benchmark it under multiple activations (sin, GELU, ReLU, Leaky ReLU, Sigmoid, Tanh, and none) on two representative tasks: Duffing and Reaction Diffuusiinn
As shown in Table 5, on Duffing RNO attains the lowest error under every activation; ReLU performs best, while sin is slightly suboptimal. Even with no activation, RNO surpasses FNO, LNO, and DeepONet by a large margin, indicating stronger intrinsic capacity for non-linear dynamics. On Reaction diffusion, the pattern differs: with Sigmoid or Tanh, RNO no longer exceeds FNO or LNO. Removing the activation improves RNO, approaching the optimum in Table 5. We attribute this behavior to redundant nonlinearities introduced by Sigmoid/Tanhh which exceed the data's representational needs and destabilize gradient propagation. In contrast, eliminating the activation severely degrades LNO, FNO, and DeepONet, underscoring RNO's superior capacity to internalize nonlinear representations.

---

Published as a conference paper at ICLR 2026
Table 5: Relative & error (1) with different activation function on Duffing and benchmarks. Lower is better. The best (bold) and second-best (underline) results are highlighted.
Relative l (L)
Dataset
Model
sin
gelu
relu
leaky relu sigmoid
tanh
None
RNO
0.1663
0.1962
0.1441
0.1481
0.2089
0.1582
0.3215
LNO(2024a)
0.3325
0.4010
0.4527
0.4828
0.4211
0.3662
0.8121
Duffing
FNO(2020a)
0.4536
0.4725
0,4555
0.4612
0.4913
0,4343
0.9296
Dec2ONet(2021)
0.6289
0.5972
0,4278
0.4717
0.9993
0.7042
0.9787
Promotion
+50.0%
+51.1%
+88.2%
+67.9%
+50.4%
+56.8%
+60.4%
RNO
0.0899
0.1011
0.1027
0.0973
0.1323
0.1482
0.1017
(220(224a)
0.1355
0.1277
0.1068
0.1111
0.1147
0.1069
0.2915
Reac-Diff.
FNO(2020a)
0.1214
0.1443
0.1337
0.1500
0.1322
0.1424
0.3511
2202ONet(2021)
0.7969
0.2742
0.2739
0.2754
1.0000
0.4337
0.6967
Promotion
+25.9%
+20.8%
+3.8%
+ +1..4%
- 15.3%
33..6%
+ 66.11
Efficiency analysis. Beyond prediction accuracy, we also assess the computational efficiency of RNO against other neural operators. Figure 6 shows that RNO offers the best trade-off between errr and training time on both Diffusion and Reaction-Diffusion. On Diffusion, achieves the lowest error with training time comparable to the fastest baselines. On both the most accurate and the fastest. Relative the widely used FNO, RNO reduces the relative & error by about 26%, trains nearly 2x faster per epoch, and uses roughly half as many parameters. Despite compact design, LNO is outperformed by RNO by more than 30% in accuracy and over 5x in training speed, showing that the gains do not rely both less accurate and slower, even with similar or larger parameter counts. Overall, RNO dominates the effcrrcc-efficiency- model-size trade-off, making it attractive for practical operator-learnin applications (see Tables in A.12 for details).
it remains
larger model. DeepONet and WNO are
0.6
LNO
RNO
DeepONet LNO
0.5 DeepONet FNO
0.30
Train time (slepoch)
Train time [stepoch)
0.25
0.07M 0.50M Eesteraaeeeio
WNO
FNO
0.3
WNO
0.15
0.2
RNO
0.10
0.1
0.7M
0.62M
0.0
0.00
2
Performance ranking worst, best)
Performance rankiing
worst, best)
Figure 6: Efficiency comparison all compared models on the benchmarks of Diffusion (left) and Reaction- Diff (right) datasets.
CONCLUSION
We presented the Riesz Neural Operator (RNO), which bridges global spectral modeling with local derivative dynamics through Riesz transform. By jointly modeling local daymics (differentiation) and the operator theory (integration), we provide an integrated formulation that gives an intuitive explanation of how RNO enhances the interpretability. Evaluations diverse PDEs, Ntoker-Stokes flows, and ERA5 data confirm consistent gains and robustness. Beyond accuracy, RNO demonstrates that embedding spectral derivatives into operator learning provides natural and principled path toward more expressive models, suggesting broad potential for advancing scientific machine learning.
10

---

Published as a conference paper at ICLR 2026
ACKNOWLEADGEMENTS
This work financially supported by the National Natural Science Foundation of China (No. 12572266), National Key Research and Development Program (2024YFF1500600), Yongjiang Talent Program of Ningbo (No. 2022A-242-Gi as well as by the High-Performance Computing Centers Eastern Institute of Technology, Ningbo, and Ningbo Institute of Digital Twin.
ETHICS STATEMENT
We adhere to the ICLR Code of Ethics pottps://iclee http:pp://iill..ccppil////////// hctpptblccccoooeeffhiiiiccie This
work involves no human participants, personally identifiable information, or sensitive content.
REPRODUCIBILITY STATEMENT
We provide detailed descriptions of benchmarks in A.18. Experimental settings are summarized in
4. The datasets used in this study are publicly available: the diverse PDE benchmarks at (https: / / thub. com/ q any ngcao/ -eeeaaaaee
Stokes with Re and ERAS at (https: : / cds. cl limmate..... coppernicuseeu// All results are averages over 3 runs with different random seeds. Definitions of the evaluation metrics provided in A.17.
{40, 500, 5000} at (https: // /ggtthbb com//nerralpperrtoo/ httpr//neuaaapprrrteeeee
REFERENCES
Robert A Adams and John JF Fournier. Sobolev spaces, volume 140. Elsevier, 2003.
Richard B Alley, Jochem Marotzke, William D Nordhaus, Jonathan T Overpeck, Dorothy M Peteet, Roger Pielke Jr, Raymond T Pierrehumbert, Peter Rhines, Thomas F Stocker, Lynne Talley, et al. Abrupt climate change. science, 299(5615) 2003.
Luca Barberi and Karsten Kruse. Localized states in active fluids. Physical Review Letters, 131(23): 238401, 2023.
Steven S. Beauchemin and John L. Barron. The computation of optical flow. ACM computing surveys 1995.
(TSUR),,7(33)443
Andrés Bruhn, Joachim Weickert, and Christoph Schnörr. Lucas/kanade meets horn/schunck: Com- bining local and global optic flow methods. International journal of computer vision, 61:211-231, 2005.
Qianying Cao, Somdatta Goswami, and George Em Karniadakis. Laplace neural operator for solving differential equations. Nature Machine Intelligence, 6(6):631-640, 2024a.
Shuhao Cao. Choose a transformer: Fourier or galerkin. Advances in neural information processing systems, 34:24924-24940, 2021.
Shuhao Cao, Francesco Brarda, Ruipeng Li, and Yuanzhe Xi. Spectral-refiner: Accurate fine-tuning of spatiotemporal fourier neural operator for turbulent flows. arXiv preprint arXiv:2405.170 172.1 2024b.
Ze Cheng, Zhuoyu Li, Xiaoqiang Wang, Jianing Huang, Zhizhou Zhang, Zhongkai Hao, and Hang Su. Accelerating pde-constrained optimization by the derivative of neural operators. arXiv preprint arXiv:2506. 13120, 2025.
Vaclav Cizek. Discrete hilbert transform. IEEE Transactions on Audio and Electroaconstics, 18(4): 340-343, 1970.
JoséL. Rubio de Francia, Francisco J Ruiz, and JoséL. Torrea. Calderón-zygmund theory for operator- valued kernels. Advances in Mathematics, 62(1):7-48, 1986.
Patrick Delorme. Formule de plancherel pour les espaces symétriques réductifs, Annals of mathemat- ics, pp. 417-452, 1998.
11

---

Published as a conference paper at ICLR 2026
Robert DiPietro and Gregory D Hager. Deep learning: Runs and Istm. In Handbook of medical image computing and computer assisted intervention, pp. 503-03-555 503-519. Elsevier, 2020.
Jianghang Gu, Ling Wen, Yuntian Chen, and Shiyi Chen. An explainable operator approximation framework under the guideline of green's 's function. Journal of Computational Physics, pp. 142244 2025.
Massimo Guiggiani and Paolo Casalini. Direct computation of cauchy principal value integrals in advanced boundary elements. International Journal for Numerical Methods in Engineering. 24(9): 1711-1720, 1987.
Maximilian Herde, Bogdan Raonic, Tobias Rohner, Roger Käppeli, Roberto Molinaro, Emmanuel de Bézenac, and Siddhartha Mishra. Poseidon: Efficient foundation models for pdes, Advances in Neural Information Processing Systems, 37:72525-72624, 2024.
Hans Hersbach, Bill Bell, Paul Berrisford, Shoji Hirahara, András Horányi, Joaquín Muñoz-Sabater, Julien Nicolas, Carole Peubey, Raluca Radu, Dinand Schepers, et al. The era5 global reanalysis. Quarterly journal of the royal meteorological society, 146(730):1999-2049, 2020.
Hossein Hodaei, Absar Hassan, Steffen Wittek, Hipolito Garcia-Gracia, Ramy El-Ganainy, Demetrios N Christodoulides, and Mercedeh Khajavikhan. Enhanced sensitivity at higher-order exceptional points. Nature, 548(7666):187-19) 2017.
Jean-Lin Journé. Calderón-zygmund operators on product spaces. Rev. Mat. Iberoumericuna, 1(3):
55-91, 1985.
Marimuthu Kalimuthu, David Holzmüller, and Mathias Niepert. Loglo-fno: Efficient learning of local and global features in fourier neural operators. arXiv preprint arXivv2500.042260, 2025.
RP Kanwal and KC Liu. A taylor expansion approach for solving integral equations. International Journal of Mathematical Education in Science and Technology, 20(3)4411-41, 1989.
Diederik P Kingma and Jimmy Ba. Adam: method for stochastic optimization. arXiv preprint 14Xiv:1412.69880 2014.
Thorsten Kurth, Shashank Subramanian, Peter Harrington, Jaideep Pathak, Morteza Mardani, David Hall, Andrea Miele, Karthik Kashinath, and Anima Anandkumar, Fourcastnet: Accelcrating global high-resolution weather forecasting using adaptive fourier neural operators. In Proceedings of the platform for advanced scientific computing conference, pp. 1-11, 2023.
Samuel Lanthaler, Zongyi Li, and Andrew M Stuart. Nonlocality and nonlinearity implies universality in operator learning. arXiv preprint arXiv:2304.132211 2023.
Yann LeCun and Yoshua Bengio. Convolutional networks for images, speech, and time series. The handbook of brain theory and neural networks, 1998.
Qian Li, Yuxiao Hu, Juncheng Liu, and Yuntian Chen. Gensr: Symbolic regression based on equation generative space. In The Fourteenth International Conference on Learning Representations.
Zhihao Li, Haoze Song, Di Xiao, Zhilu Lai, and Wei Wang. Harnessing scale and physics: A
multi-graph neural operator framework for pdes on arbitrary geometries. In Proceedings of the 31tt ACM SIGKDD Conference on Knowledge Discovery and Data Mining V. J, pp. 729-740, 2025.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametrie partial differential equations. arXiv preprint arXiv:2010.08895 2020a.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Graph kernel network for partial differential equations. arXiv preprint 200.:000000008888 2020b.
Zongyi Li, Miguel Schh-cciafiii,, Nikola Kovachki, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Learning dissipative dynamics in chaotic systems. arXiv preprint arXiv: 066.06898, 2021.
12

---

Published as a conference paper at ICLR 2026
Zongyi Li, Nikola Kovachki, Chris Choy, Boyi Li, Jean Kossaifi, Shourya Otta, Mohammad Amin Nabian, Maximilian Stadler, Christian Hundt, Kamyar Azizzadenesheli, et al. eeometry-informed neural operator for large-scale 3d pdes. Advances in Neural Information Processing Systems, 36:
35836-35854, 2023.
Zongyi Li, Hongkai Zheng. Nikola Kovachki, David Jin, Haoxuan Chen, Burigede Liu, Kamyar Azizzadenesheli, and Anima Anandkumar. Physies-informed neural operator for learning partial differential equations. ACM/JMS Journal of Data Science, 1(3):1-27, 2024.
Chaoyu Liu, Davide Murari, Lihao Liu, Yangming Li. Chris Budd, and Garola-Bibiane Schönlieb. Enhancing fourier neural operators with local spatial features. arXiv preprint arXiv:2503.17797, 2025.
Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via deeponet based on the universal approximation theorem of operators. Nature machine intelligence, 3(3):218-22) 22-22,, 2021.
Lawrence Marple. Computing the discrete-ting timeeettmemmmee analytic" signal via fft. IEEE Transactions signal processing, 47(9):2600-2666 1999.
Sophocles J Orfanidis. Introduction to signal processing. Prentice-Hall, Inc., 1995.
Thomas Osler. Leibniz rule for fractional derivatives generalized and application infinite series. SIAM Journal on Applied Mathematics, 18(3):658-674, 1970.
HB Pacejka and IJM Besselink. Magic formula tyre model with transient properties. Vehicle system dynamics, 27(S1): 249, 1997.
Yuan Qiu, Nolan Bridges, and Peng Chen. Derivative-enhanced deep operator network. Advances in Neural Information Processing Systems, 777:20945-209883 2024.
Maziar Raissi, Paris Perdikaris, and George Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686-707, 2019.
Bogdan Raonic, Roberto Molinaro, Tim De Ryck, Tobias Rohner, Francesca Bartolucci. Rima Alaifari, Siddhartha Mishra, and Emmanuel de Bézenac. Convolutional neural operators for robust and accurate learning of pdes. Advances in Neural Information Processing Systems, 36: 77187-77200, 2023.
Olaf Ronneberger, Philipp Fischer, and Thomas Brox. U-net: Convolutional networks for biomedical image segmentation. In Medical image computing and cosputer-assisted Intervention-MIC. 2015: 18th international conference, Munich, Germany, October 5-9, 2015, proceedings, part IIl 18, pp. 234-241-1 241. Springer, 2015.
Jonathan D Smith, Kamyar Azizzadenesheli, and Zachary E Ross. Eikonet: Solving the eikonal equation with deep neural networks. IEEE Transactions on Geoscience and Remote Sensing, 59 (12):10655-1066, 2020.
Samuel D Stearns. Of aldapfive signal processing. 1985.
Atsushi Tamada and Michihiro Igarashi. Revealing chiral cell motility by 3d riesz transform- differential interference contrast microscopy and computational kinematic analysis. Nature Com- munications, 8(1):2194, 2017.
Karn Tiwari, Niladri Dutta, NM Krishnan, et al. Latent mamba operator for partial differential equations. arXiv preprint arXiv:2505.19100 2025a.
Kam Tiwari, NM Krishnan, et al. Cono: Complex neural operator for continous dynamical physical systems. APL Machine Learning. 3(2), 2025b.
Alasdair Tran, Alexander Mathews, Lexing Xie, and Cheng Soon Ong. Factorized fourier neural operators. arXiv preprint 1rXiv:21/1.13802 2021.
13

---

Published as a conference paper at ICLR 2026
Alasdair Tran, Alexander Mathews, Lexing Xie, and Cheng Soon Ong. Factorized fourier neural operators. In The Eleventh International Conference on Learning Representations, 2023. URL
https: /ttpr///eppeeeeeeeeeeeeee .eettw..eee//e.e..eee Eorum?id=:m11IP11---
Tapas Tripura and Souvik Chakraborty. Wavelet neural operator: neural operator for parametric partial differential equations. arXiv preprint arXiv:2205.0211 2022.
Michael Unser and Dimitri Van De Ville. Wavelet steerability and the higher-order riesz transform. IEEE Transactions on Image Processing, 19(3):636 (000)66666622202 2009.
Michael Unser, Daniel Sage, and Dimitri Van De Ville. Multiresolution monogenic signal analysis using the riesz-laplace wavelet transform. IEEE Transactions Image Processing, 18(11): 2402-2418, 2009.
Neal Wadhwa, Michael Rubinstein, Frédo Durand, and William T Freeman. Riesz pyramids for fast phase-based video magnification. 2014 IEEE International Conference on Computational Photography (ICCP), pp. 1-10. IEEE, 2014.
Hanchen Wang, Tianfan Fu, Yuanqi Du, Wenhao Gao, Kexin Huang, Ziming Liu, Payal Chandak, Shengchao Liu, Peter Van Katwyk, Andreea Deac, al. Scientific discovery the age of artificial intelligence. Nature, 620(7972):47-60. 2023.
Haohan Wang, Xindi Wu, Zeyi Huang, and Eric P Xing. High-frequende component helps explain the generalization of convolutional neural networks. In Proceedings of the IEEE/CVE conference computer vision and pattern recognition, pp. 8684-8694, 2020.
Rui Wang, Yihe Dong, Sercan Ö Arik, and Rose Yu. Koopman neural forecaster for time series with temporal distribution shifts. arXiv preprint arXiv:2210.0367,, 2022.
Steven Weinberg, Eikonal method in magnetohydrodynamies. Physical Review, 126(6):1899, 1962.
Haixu Wu, Tengge Hu, Huakun Luo, Jianmin Wang, and Mingsheng Long. Solving high-dimensional pdes with latent spectral models. arXiv preprint arXiv:2301. 2023.
Zipeng Xiao, Zhongkai Hao, Bokai Lin, Zhijie Deng, and Hang Su. Improved operator learning by orthogonal attention. In Proceedings of the 41st International Conference on Machine Learning. 54288-54299, 2024.
Xiong Xiong, Kang Lu, Zhuo Zhang, Zheng Zeng. Sheng Zhou, Rongchun Hu, and Zichen Deng. High-frequency flow field super-resolution via physics-informed hierarchical adaptive fourier feature networks. Physics of Fluids, 37(9), 2025.
Zhi-Qin John Xu, Yaoyu Zhang, Tao Luo, Yanyang Xiao, and Zheng Ma. Frequency principle: Fourier analysis sheds light on deep neural networks. arXiv preprint arXiv: 1901 1901.06333.,, 2019.
14

---

Appendix Contents
Supplementary Material for RNO
A.1
Why local direction is important for PDEs?
16
A.2
First-order approximation
18
A.3
Proposition: plug-in stability for temporal learning from spatial locality
19
A.4
Green function and RNO
19
A.5
Orthogonality
20
A.6
Bound
20
A.7
Relationship between the Hilbert and Riesz transforms
21
A.8
Differences between RNO and FNO + CNN
22
A.9
Directional derivatives and full-spectrum modeling
23
A.10
Implementation of the direction mixer
24
A.11
Directionality in one dimension
24
A.12
Efficiency analysis
25
A.13
Derivative fidelity
25
A.14
High-frequency baselines
26
A.15
Robustness to noisy inputs
26
A.16
Implementation details
A.16.1
Model architecture
27
A.17
Evaluation metrics
29
A.18
Supplementary for benchmarks
29
A.18.1
Duffing oscillator
29
A.18.2
Eulcr-Bernoulii beam
29
A.18.3
Diffusion equation
29
A.18.4
Reaction-diffusion system
30
A.18.5
Brusselator reaction diffusion system
A.18.6
Havier-Stokes equation
A.18.7
ERA5 dataset
32

---

Published as a conference paper at ICLR 2026
A APPENDIX
A.1 WHY LOCAL DIRECTION IS IMPORTANT FOR PDES?
Local derivatives inherently directional, with directionality encoded by the phase gradient Vø. In the high-frequency regime, PDEs couple to the data through their principal symbol evaluated Vá. The Riesz transform provides rotation-equivariant estimates of local phase and orientation, thereby supplying neural networks with precisely the quantities that govern PDE propagation. We next analyze each step in detail.
f(x ((x+)(1))+ f(x) (x)
Taylor expansion:
dr
(xx,1)
f(x)
r(t) of(x)
+
Intensity
First-nder expression
Tresee Trandarre
130
Abrupt change along time
Amplitute
No abrupt change information
Prepary
Original interitty
Tuplor upprosimation
Fourier
Finst derivative
Figure 7: The specific structure Taylor expansion and its difference from Fourier transform.
1. LOCAL DERIVATIVES ARE DIRECTIONAL
Given Rd C, the directional derivative along unit vector is
unu(x) n Vu(x). Write Ae with 2 0. wherever u + 0. Then
8, e (8,A iA 0) = u (0, In +i
(12)
Consequently, the normalized directional derivative decomposes neatly into amplitude and phase components:
0,hnA = RR On u
0 = Im
(13)
Accordingly, local derivatives are inherently orientation dependent, with Vo prescribing that orienta- tion.
DIRECTIONALITY IS ENCODED THE PHASE
From Vu = VVVnn A+
we obtain
Vu
Vu
Im
Vln Re
(14)
u
16

---

Published as a conference paper at ICLR 2026
The vector k(x) = Vx((x is the local wavevector and
Vo(x) [((()]
n(x)
is the local orientation. In the high-frequency regime, where amplitude varies slowly compared with phase,
|uu| = Volll
0n an = iA (0nØ) e
so the derivative is dominated by directional phase change.
3. WHY PDES CARE: PRINCIPAL SYMBOL AND HIGH-FREQUENCY LINK
Consider m-th order linear differential operator
P(x, 0) = Σ (x(x) 0n,
(x, Σ (a(x) (iE)a
|o mm
with principal symbol(Li al.) Ppr- Insert ansatz u(x) A(z))ewo(()) (w 1). By repeated Leibniz(Osler, 1970),
P(Aclwo) el ( Pp Vo) A O
(15)
At leading order, the principal symbol ties the dynamics to the phase gradient Vo. Asymptotic solvability at order wh entails the eikonal equation(Smith et al., 2020),
Ppe(z, Vo(z)) =0.
(16)
Hence, rays propagate along the phase gradient Vo. For the transport operator P = 0 a(x) Vz.
the principal symbol is Ppr(I; T) = i(T + a E). Imposing the eikonal condition (Weinberg, 1962)
Ppr(x: (x. a. V = gives t + a- = 0 with g = V so the spatial propagation direction is aligned with V V.
4. RIESZ TRANSFORM: EXTRACTING PHASE AND ORIENTATION
The d-D Riesz transform Ru (Riu,... Rau)) is defined by R,u(w) - Wj i(w).
|wil
Plane wave benchmark. For u(x) = cos(k-z))
k
Ru(x) sin(k.r:)
|k]
Hence Ru points in the wave direction k/ |k| and is 7/2 out phase.
High frequency. Let u(x) = A(x) cos(w ((x)) with smooth A, ø and w >> 1. Locally near ((x) R ((xo) + Vo(an) (x zo). i.e. a plane wave with k = ~ (a))) Then
Volao)
VA + +1
Ru(zo)= A(xo) sin (wo(w))))
(17)
Therefore the orientation is
Ru(x)
VΦ(z)
n (x)
[
and cotation-invarian local phase can be defined by the monogenic signal (Unser et al., 2009).
p(x) /u(x))2 [u,(2) |
((x) = atan2(I|Ru(z)II, u(x)),
so that locally cooos 0, ||Ru|| psin 0.
17

---

Published as a conference paper at ICLR 2026
5. THE RATIONALITY OF USING RIESZ TRANSFORMATION FOR NEURAL OPERATOR
Aligms with principal symbols. By equation 15 15-equation 16, the leading- order response of the PDE depends Vø through Ppr (2, Vő). Riesz features extract the orientation n R Vo/ Vól and the quadrature phase 8, thereby granting the network direct access to the directional variables relevant the PDE.
Rotation covariance. Since R is rotation-covariant, (u, Ru) provides features that transform predictably under rotations. This property is essential for modeling anisotropic diffusion, advection, and wave propagation without relearning all orientations.
Phase emphasis and contrast robustness. Phase (via remains invariant under multiplicative contrast changes u + cu, while the derivatives in equation 13 accentuate 0, Consequently, learning guided toward geometry and propagation rather than raw amplitude.
Differentiable and lightweight. R is linear with fixed multiplier -ig/ll<] and its discrete realizations are convolutional and hackpropagation-fine frrennddl. can serve either as a fixed front-end as parallel stream complementing learned filters.
INFLUENCE OF PHASE AND ORIENTATION
To highlight how phase governs directionality, we conduct an ablation on the Navier-Stokes dataset. Under identical network architecture, we compare two input parameterizations (cf. Figure 8(a)): (i) complex features that retain both magnitude and phase, A(z)e*((*); and (ii) magnitude-only features with the phase suppressed, A (x))e=0 All other training and evaluation settings are held fixed.
As shown in Figure 8(b)-(d), removing phase yields markedly sparser fine-scale structures and substantial misalignment local feature orientations. rollout regimes (Cao et al., 2024b) (e.g.. NO), such local directional data with more numerous localized structures, resulting in larger long- long-horizon errors, as evidenced by Table 2. Moreover, the Riesz transform precisely enhances the directional sensitivity of phase control.
accumulate over time. The effect is amplified for more complex
Transform
Trensoom
A(z)ezp(>á(x))
A(z)cap(i 0)
transform
transform.
Figure 8: (a) Two feature structures; (b) results using magnitude sranitude-ppase features; (c) results using magnitude-only features; (d) orientation comparison selected details.
A.2 FIRST-ORDER APPROXIMATION.
Assume that the image profile f is twice continuously differentiable and that the displacement y(t) is small (Wadhwa al., 2014). Since
[x f(x+
we expand / at position using Taylor's theorem with remainder:
f(x + 7()) f(x) ++))f((x) 7(1²
("(((,))),
for some {(x,1) () between x x + 7(). [f"(=)] < in the relevant neighborhood, then second-order remainder satisfies
(t)
M
('(((x,)))
< 2
lntt)
Hence, whenever (t(t) is small, the remainder is of order (+)) which is negligible compared to first-order term order (t)). Therefore we obtain the first-order approximation
df(x)
I(x, t) = f fx + ((t) = f(x) + (t)
dx
18

---

Published as a conference paper at ICLR 2026
Intuitively, this means that for sufficiently small displacements the change in intensity is well described by a linear response with slope (z), while higher-orde order (curvature) effects only appear at quadratic (and smaller) scales in y(t).
A.3 PROPOSITION: PLUG-IN STABILITY FOR TEMPORAL LEARNING FROM SPATIAL LOCALITY.
Let u : Ω x [0.T + R solve a local evolution law of the form
t,, F(z, t, (x) f)),
((()) = (u, Vu, u),
(18)
on a bounded domain Q Rd, where F is uniformly locally Lipschitz in its last argument: there exists L > 0 such that
F(z,t, F(z,t,z)) s - 2'11 for all (r, E Ω x [0, 2, E RS.
(19)
Assume learned spatial locality representation r(x, () approximates (z, t) with
(u) ( (1 e, s 00.
(20)
Then for any there exists a (pointwise) network go 1R442+0 +R such that 10((r, Ax, 1, (ui) (=(0))007))) S Le 6.
(21)
Proof (sketch). Define G(a,,,,) = F(x, 1, 2). By the Lipschitz property above,
G(I,t,r) G(z,t, v()))
L ||r (w)]l
Taking LP norms yields the term LE by the assumption on r. Continuity of G on com- pact subsets implies, by universal approximation, the existence of a small MLP go such that sUP(z,C.)EX (1r, z) G(æ,, (, 2)1 < 8 on the relevant compact set K. A triangle inequality then gives the stated bound.
In regimes where
0, alz, t) (n(x, t) ()),
(22)
with orientation t) E gd-1 and scalar speed a(x, t), let yp n-Vu and suppose ||@ -pl|e>
a small regressor ag (x, f, estimates a with |ā - L= - I 17. then
||a|| + || || ā ر S
(23)
Thus, temporal accuracy follows directly from the fidelity of spatial directional derivatives and a scalar modulation fit.
Remark (why the RNO representation satisfies the assumptions). The Riesz-based features Ru (.. R.u), with multipliers R2,(E) - ((E), form pnergy-preserving system: Σd_ R |ull? (Plancherel). Hence Ru provides stable, grid-agnostic surrogate of spatial locality (values, gradients, orientations). Once r ~ (u()) (or r z (u, Ru)) holds to accuracy e, Proposition (P4) shows that the non-stationary temporal signal d.u can be learned by low-capacity pointwise map with error bounded by Le + 5. In short: once spatial locality is captured adequately, temporal non-stationarity becomes a stable plug-in regression problem.
A.4 GREEN FUNCTION AND RNO
We recall that many linear (or linearized) PDEs admit a Green representation (Li et al., 2020b) Lu = + u(x) fg G(x, y) f() dy, so the solution operator G ff + u is an integral operator with kernel G. Neural operators approximate such maps by learning an integral kernel Ky so that
(Nof)(x) / (2, v) (1v) dy Re / G(x, v) f(y) dy.
In Fourier Neural Operators, one further assumes ko 4) = ko (x v) and parameterizes its Fourier transform by multiplier Mc(E); the operator becomes (No.) (E) Mo(E) f(E), which is still
19

---

Published as a conference paper at ICLR 2026
bounded convolution operator with (learedd) Green-tyy type kernel Kg that is independent of the discretization..conne hwwp::/ow...oooeeeeecceee/eeeeee
RNO augments this construction with Riesz-based directional derivatives while preserving the same operator class. The ith Riesz transform can be written as a singular integral R, = r; * f with kernel (z) ==//==/=+1 and Fourier multiplier (5) = m.(E) f(E). where m.(5) - ell bounded and homogeneous of order O. Thus R2 a Galderón-Zyggmnn
operator (Journé, 1985), bounded on 7.2 and on Sobolev spaces (Adams & Fournier, 2003), and does not destroy continuity of the solution operator. In RNO we apply the neural operator to the augmented field g = f+ R.f.
so that
(A ()x)) = [ nz( u) dly)dy = S (o(x, v)f(v) dy + Σ my / no(s, v) (Rtf)(v) dy.
Writing R&f = r; * and exchanging integrals gives
(AV (xx) Go(x, 2) (((), [()) Go(a, z) = 2)+ Σ
kin(x rely dy,
so NBRNO again an integral operator with learned kernel Go of Green-function type. In the Fourier domain, this is simply a new bounded multiplier Mo(E) = M1(E) + Σ, wfm, Hence
the Riesz term acts as directional, high-frequency reweighting that injects local anisotropic variation into the representation, while the overall map remains continuous, resolution-invarinn integral operator in the same functional class the original Green-kernel construction.
A.S ORTHOGONALITY
Inner product. Let n 2 2 and wave-vector k (ki,.... kin) E R and magnitude |[k]] -i k22 For any distinct indices + j consider the L2 inner product
(R,, = [ (2(k) (k) dk.
(24)
Then setting E(k) 12 (an even function every ke) yields
(R, Rí) = [ kiky R" |k||2
E(k)
dk.
Odd-even symmetry argument. Fix Fix j. Under the sign change k; ++ -k; the factor kęky flips sign while both ||k| and E(k) remain unchanged; hence the integrand is odd with respect to the k; axis. Because the integration domain R" is symmetric, the integral vanishes:
RR) = 0,
Vizj
(25)
Spatial-domain equivalence. By Parseval's theorem (Orfanidis, 1995), the frequency-domaii result implies spatial orthogonality: (R,f, RRfx 0 for + j. where R, denotes the real-space Riesz transform in direction i. Thus the vector of Riesz components R = (Ri,..... Raf)) forms an orthonormal frame in L2(R²) Consequently, neural operators that process the Riesz vector component-wise inherit stable, non-redundant latent representations in any spatial dimension d.
A.6 BOUND
First-order approximation. The first-order Taylor expansion provides an approximation to the actual intensity function f(x) and its variations, which is reflected in the Riesz transform within the spatial domain (Wadhwa et al., 2014). For rapidly varying data, the accuracy of the first-order Taylor expansion decreases increases. Large values of (>(() may introduce artifacts. To ensure the approximation remains valid, we derive an upper bound for Ç.
Assuming the first-order Taylor expansion holds, the modified intensity is:
f(x + (y((t) x f(x) + (T()) df(x) dx
(26)
20

---

Published as a conference paper at ICLR 2026
The Riesz transform introduces 90-degree phase shift. If f()) cos(wx)), we have: cos(wx + (wr(t)) = cos(wx) (w>(()) sin(wx).
(27)
Using the trigonometric identity:
cos(wr + (w)(()) = cos(wer) cos((wrr()) sin (war) sin((w>(t)),
(28)
require (wr(t) to be small for cos( (t)) 1 and sin((w>((t) (wv(t). Assuming 5%
error tolerance, the second-order term should be less than 5% of the first-order term, yielding the condition:
A
St(t) 12
(29)
where À = 2st is the spatial wavelength. A larger wavelength allows for larger scaling factor Ç; smaller wavelengths require smaller Ç to avoid artifacts. While this bound on is theoretical, wavelength of the data must be determined based on the signal and network's characteristics.
Multidimensional case. For a multi-dimensional setting, let x = (œ1, .., represent the spatial coordinates. The first-order Taylor expansion of f fax becomes:
0f(x)
f' + +x(() ~ { Σ 7:(1)
+
(30)
Or.
Here, (t) the scaling factor in the i-th direction. The scaling factor Ç still amplifies the signal's variations across all dimensions. The error condition becomes:
Σ wer((t)
(31)
where wi the frequency in the i-th direction. The bound provided here is not an exact limit but rather a theoretical upper bound. Since it depends on the angular frequency of the data, the boundary is not fixed. Accordingly, in the mixer A.10 we replace with learnable parameter while constraining its range.
A.7 RELATIONSHIP BETWEEN THE HILBERT AND RIESZ TRANSFORMS
In one dimension, the Hilbert transform Hf (Cizek, 1970) the singular integral
1 p.v. f(y) dy,
Hz)
π
y
with Fourier representation
[ f(E). where "p.v." denotes the Cauchy principal value (Guiggiani Casalini, 1987). In R", the Riesz transform generalizes this construction: for f S(P").
R,f(x) =
f(5) (z).
== L.....,..
so each R, acts as normalized directional derivative with degree-zero Fourier multiplier. For 1, taking & = yields
(I)(x)=
i(E) (x) F-1 [
E
sgn(E) S(E) (a),
recovering exactly the Hilbert transform: R,f H.f.
Collecting the components Rf = (Rif..... Rrf) gives a vector operator that an L L'-iismetry:
Σ |R 1|
Thus the Riesz transform the canonical multi-dimmensonall analogue of the Hilbert transform, widely used to extract directional and phase information in higher-dimensional signal analysis. Making this connection explicit sharpens our understanding of how Riesz transforms operate on derivative- -iie quantities and underpins the RNO framework.
21

---

Published as a conference paper at ICLR 2026
(a) Classical Neural Network: Directly Map
(b) FNO/LNO: Global Transforrttinnn in Fourier/Laplace
Space
Encoder
Decoder
Companents
Mapped Components
R/L
I
F(x)
Fo(e)
Global Transformation
(c) RNO: Local Transformation Riesz Space
Local block
R1
R₂
Global block
Mixer
Components
Mapped Components
Relobal
F(R((
Global Transformation
Figure comparison between RNO and simplified versions of some classical architectures, highlighting the global-local mixer of RNO.
A.8 DIFFERENCES BETWEEN RNO AND FNO + CNN
Resolution invariance and local-global representation. A key advantage of the Fourier Neural Operator (FNO) is resolution invariance: it approximates a continuous operator T grid-agnostic spectral multiplier. Writing one FNO layer as
W++1 of (F(FF(((C (E) +b):
(32)
the learned multiplier (E) evaluated on the Fourier samples of any discretization. Equivalently. the sampling map S₁, approximately commutes with Tasi
Sp TM Rs T00 Sns
(33)
so training and inference can use different resolutions with minimal drift (apply. then sample s sample, then apply).
By contrast, FNO+CNN (Kalimuthu et al., 2025; Liu et al., 2025) augments equation 32 with spatial convolution (K uI( whose kernel E & is defined in pixel units. Under grid change h-th', the effective continuous kernel rescales Kд (x) =<x/h²), hence Kīz (E) depends on h'. This breaks equation 33: the conv branch is not grid-agnostic, making cross-resolution generalization fragile.
We retain the operator view and expose directional locality spectrally via the Riesz transform:
& RTu(E) = |Kll
((E),
& 0,
(34)
scale-normalized directional derivative. Intuitively, the spectral path preserves continuity across grids, while Riesz channels add anisotropic, have-weverumber detail that complements global low-frequency modes.
Complexity and redundancy. Let N= denote grid points and Ciny Cout channel counts. An FNO layer with |A] retained modes costs
ONN log Cin) + [(AI - Cin Comt)
FFF/FFFT
mode mixing
Adding k CNN (LeCunn & Bengio, 1998) head introduces
O(k2 N Cin
local CoIY
22

---

Published as a conference paper at ICLR 2026
Table Comparison of model effiency on the Reaction-Diffusion dataset.
Model
Params
Train time (s/epoch)
Peak memory (MB)
Relative &2 +
E-FNO
793985
0.166
241.07
0.1131
loglo-FNO
0.142 764225
617.22
0.1009
RNO (ours)
172229
0.052
111.17
0.0899
plus meature-map memory N C. For typical k E [3, 7] and moderate |A], this CNN cost dominates as resolution grows, and features remain resolution-bone bouuodd Moreover, since filter responses Hm (E) are unconstrained, concatenating CNN and spectral features often duplicates energy in overlapping bands, inflating rank unless explicit decorrelation is enforced.
In contrast, adding d Riesz channels equation 34 costs
ON log N Cin) +
(CunnCoun)
shared RI/IRT
clementwise multipliers
light mixer
with d = (Q) E (2,3). The overhead is linear in d and grid-agnostic; no extra k2 N spatial conv introduced, and memory adds only d feature maps plus small mixer. The experimental results are shown in the 6, RNO achieves faster speed and higher accuracy.
Why Riesz features are energy-orthogonal, while CNN features can be redundant. By Plancherel (Delomme, 1998),
Σ R
( JR
Σ E
2
-|u1|
so the vector Riesz transform Ru (R1 Rau) an L isometry: it redistributes the same energy across directions without gain or loss. On band-limite discrete grids, the identity holds up to negligible aliasing, and channels are approximately decorrelated for isotropic content. Intuitively, we split energy by orientation rather than duplicating it.
For CNN head with filters {Km IC" and responses (L)()) the feature energy at E scales as Em Hm (5) (0()) Without spectral constraints, multiple filters can attend the same bands; Gram matrix Gwnn (Kem, Ki)) is non-diagonal, yielding correlated channels and effective rank inflation. Our orthogonality study (Fig. 4.2) and analysis in App. A.5 support this view.
A.9 DIRECTIONAL DERIVATIVES AND FULL- SPECTRPMT MODELING
As discussed in 3.3, derivative features central how RNO couples global and local dynamics. We now provide a concise, implementation-Ieve description that links full-spectrum modeling, rank reduction, and learning dynamics.
Representation. In the spectral formulation, the i-th Riesz channel is obtained by the fixed multiplier mk(k) -i, [( implemented via FFT as discussed in A.16.1. The resulting vector of channels i
{Ri(g)) behaves as a normalized first-order directional derivative. In the continuous limit, the transform is energy preserving and components are orthogonal; on discrete, band-limited grids this property holds approximately. Intuitively, the global integral kernel captures smooth structure [xx, whereas anisotropic variations concentrate in the Riesz channels. Together they provide a compact, non-redundant coverage of the full spectrum.
PDE and optical flow. 2023):
PDEs characterize field by constraining its local derivatives (Wu et al.,
Flx,,, 4, u(a, !), 0,u(x, 1), Vu(a, ().x,,
(35)
where 0z, denotes the derivative of u w.r.t. 2;. While derivatives quantify space- time variations, the induced solution operators may couple information globally or scale-locally. Hence, useful representation should make both global modes and orientation-awee awaee local changes explicit, rather than assuming pure locality.
23

---

Published as a conference paper at ICLR 2026
A concise local view is the optical-flow-tyye type (Beauchemin & Barron, 1995) transport along charac- teristics.
a, + vu(zz)
tu) S(,, + Au(x,t),
(36)
whose brightness-constancy case is S= v == This form cleanly separates: (i) the directional derivatíve v Vu, concentrating localized (high-wavent) wavenhaaveumbrr variations; and (ii) the smoothing terms, which induce global or scale-local couplings (low wavenumbers). Guided by this principle, full-spectrum representation should preserve global spectral efficiency for low frequencies while explicitly exposing directional derivatives to complement the global description and mitigate spectral normalized first-order derivatives in the spectral domain- this purpose without committing to a specific PDE class.
bias. The Riesz
Structured rank reduction. The global pathway performs a low-rank approximation in the Riesz basis, which is well-suited for smooth, slowly varying dynamics. In contrast, localized or anisotropic features require high rank in isotropic bases but admit a sparse, orientation-aware encoding in centiiinal-deivvatii channels. RNO exploits this by replacing large local kemel with small directional mixerA. 10: few weights (Wi...., wa) acting on {R.(g)).
2 = Zglobal
+
wi (g)
This realizes a structured "low-rank directional-sparae decomposition: low-rank global modes for low frequencies, and light, direction-sparse mixer for local high-frequency content.
How directional derivatives are learned. Differentiating L w.r.t. w, yields
ac
< aC
0= R.())).
Ow;
i... the residual is projected onto each directional channel. Because {R: are approximately orthogonal, these projections decouple, and training selectively amplifies the oriented components that most reduce error. Moreover, Riesz features (A.1) are rotation-cove covariant: rotation of the input induces the same rotation of the channel vector, so the learned mixing generalizes across orientations without relearning. A bounded amplitude (A.6) prevents excessive high-frequency gain, ensuring numerical stability.
A.10 IMPLEMENTATION OF THE DIRECTION MIXER
The direction mixer essentially performs global-local integration of physical dynamics in the Riesz space. This differs from mixing in the spatial domain, since data representations in the Riesz space are inherently distinct, and can be regarded as latent feature representation. Consequently, the challenge to design a fusion mechanism that preserves data characteristics while naturally integrating physical properties. In the physical domain, the monogenic signal (Unser et al., 2009), high-dimensional extension of analytic signals, offers such an example, as it naturally fuses orthogonal components across scales. For monogenic signal, its general form in R can be expressed
f(x) Φ 222 Tzj.
(37)
Our implementation follows this formulation with several modifications:
Mixer Rigter= + w. wyRj,
(38)
where w denotes learnable, direction-specific parameters that assign spectral weights the cor- responding Riesz components. For n--dimensional data, the formulation naturally extends to n dimensions. To avoid excessive artifacts according to the derived bounds, the values are constrained within (0, 2) according to bound in A.6. This structure enables natural fusion of global and local information, which is key the superior performance of the overall model.
A.11 DIRECTIONALITY IN ONE DIMENSION
In one spatial dimension, the Riesz transform reduces to the Hilbert transform (Cizek, 1970), the canonical operator for building analytic signals (Marple, 1999):
f(T)
Htf = p.v.
dr.
(39)
t
24

---

Published as a conference paper at ICLR 2026
Table Comparison of model size, training efficiency, memory footprint, and accuracy (relative l error) for existing operator-learning baselines and our RNO on the Reaction-Diis dataset.
Model
Params (M) 4
Train time (s/epoch) 1
Peak memory consumption +
Relative & +
DeepONet
567809
0.263
246.30 MB
0.7969
FNO
310913
0.098
230.67MB
0,1214
WNO
192833
0.112
281.11 MB
0,1713
LNO
64241
0.298
1635.57 MB
0.1355
RNO (ours)
172229
0.052
111.17MB
0.0899
Table 8: Comparison of model size, training efficiency. memory footprint, and accuracy (relative l error) for existing operator-lcarning baselines and our RNO on the Diffusion dataset.
Model
Params (M)
Train time (s/epoch) 1
Peak memory consumption
Relative la
DeepONet
403073
0.285
365.67 MB
0,1356
FNO
527745
0,278
313.36 MB
0,0229
WNO
192833
0.112
281.11 MB
0,0179
LNO
76369
0.512
8805.78 MB MB
0.0081
RNO (ours)
150609
0.462
245.71 MB
0.0079
While single spatial dimension lacks intrinsic notion of orientation, the Hilbert transform nonethe- less decomposes any real-valued signal into two oppositely directed components, distinguished by the sign of their temporal frequencies. Define
(+) =
(40)
f-(t) --000-1)110)))))
(41)
so that positive frequencies correspond to energy (or information) propagating away from the source, while negative frequencies correspond to energy returning toward The spectral imbalance between these two sectors thus reveals the net propagation direction. In fluid dynamics, spectra dominated by positive frequencies typically indicate pressure waves radiating from point source. When solutions of the Vavier-Stokes equations predict strong vortical feedback, analysis should instead focus on the band. In contrast, in purely diffusive processes where concentrations relax from high to low, the positive frequency component typically dominates.
A.12 EFFICIENCY ANALYSIS
Full model efficiency results: on Reaction-Diic (Table 7), RNO achieves the lowest error (0.0899), improving over FNO 1214) and LNO (0.1355) by 26% and 34%, while being the fastest /0.o52s/epoch;; ~ 2 vs. FNO, 5X vs. LNO). It remains compact (1.7x fewer parameters than FNO) and memory-efficiae 13x below LNO). On Diffusion (Table 8). RNO again attains the best error (0.0079) slightly ahead of LNO (0.0081) and well ahead FNO/WNO (65%/56% /56%). with only a modest training time overhead; the model is still small (1.5 10S parameters, 3.5x fewer than FNO) and memory usage comparable to DeepONet and markedly below LNO. Overall, the tables indicate that RNO offers the most favorable accuacy--fliciecy--memee
(about half the peak memory of Fourier-basse baselines and nearly
trade-off among compared operators.
A.13 DERIVATIVE FIDELITY.
To further evaluate RNO's ability recover higher-order quantities, examine the reconstruction of physically meaningful first-order derivatives on the Navier-Stokes benchmark at Re = 5000 by comparing vorticity Ug 0,u, (middle row of Fig. 10) and divergence 5 0_u, + "y (bottom row). These quantities separate rotational and compressive components and are sensitive to high- wavenubbermbrr errors around shear layers and filaments.
In the vorticity panels, RNO recovers thin vortex sheets and interface polarity with high fidelity, exhibiting fewer ringing artefacts than LNO and less oversmoothing than WNO/FNO; DeepONet
25

---

Published as a conference paper at ICLR 2026
Table 9: Comparison of relative & errors between RNO and several high operator baselines (WNO, LNO, loglo-FNO, and E-FNO) on the figh-frrqency-enhaacee benchmarks.
and Brusselator
Model
Reac- Diff.
Brusselator
WNO
0.1713
0.1789
LNO
0.1355
0.1858
loglo-FNO
0.1009
0.1679
E-FNO
0.1131
0.1703
RNO (ours)
0.0899
0.1317
largely misses fine eddies. In the divergence panels, RNO produces near-zero residuals-ccnsistenn with incompressiblity-wwaee strongly curved structures.
the baselines show larger coherent divergence patches, especially
The qualitative ranking in Fig. 10 indicates that the proposed directional-spectral representation improves derivative recovery practice, capturing rotational structure while suppressing non-physical divergence.
Ground Truth
LNO
WNO
FNO
DeepONet
Figure 10: Comparison of high-frequency components(top), vorticitymmiddle), divergence(bottomm) with the Ground Truth.
A. 14 HIGI--RREQUENC BASELINES.
To further demonstrate the advantages of our method over other approaches that enhance local high- frequency components, we added a high-frrequncy baseline comparison on the and Brussclator benchmarks. The baselines include LNO, WNO. E-FNO, and loglo-FNO. Table summarizes the relative & errors on the Reaction-Diffusion and Brusselator benchmarks. Our RNO consistently achieves the lowest error both datasets: it improves and from 0.1679 to 0.1317 on Brusselator. Other
Diffusion
the best baseline (loglo-FNO)
from 0. 1009 to 0.0899 on Reaction-Diffusin variants such as WNO, LNO, and E-FNO are clearly less accurate, indicating that the proposed RNO architecture transfers better across different nonlinear reaction-diffusion systems.
A.15 ROBUSTNESS TO NOISY INPUTS.
To test whether the spectral derivative in RNO amplifies high-frequency noise, we inject additive pero-mean Gaussian noise at relative levels a E (o, 0.05, 0.1, 0.2) into the Reaction-Diffusion inputs test time (Table 10; is the standard deviation normalized by the input scale).
26

---

Published as a conference paper at ICLR 2026
Table 10: Relative &g error on the Reaction- dataset under different noise levels (SNR).
Model
SNR
0.05 SNR
0.1 SNR
0.2 SNR
FNO
0.1214
0.2129
0.3137
0.4976
LNO
0.1355
0.2632
0.6622 0.0934
0.3719
RNO
0.0899
0.0926
0.0958
RNO attains the lowest relative & error for all and degrades only mildly: (08999+0.099 0.0958 (<7% increase). In contrast, FNO and LNO grow substantially with noise, from 0.1214 0.4976 and 0.1355 - 0.6622, respectively (roughly -10% and 494%-389% across a:).
The stability aligns with RNO's direetional-derrvaatvv representation: Riesz multipliers mE)
- ((llll satisfy Σ, ( =
avoiding uncontrolled high- bighnwauenumber gain, while the learned mixer weights gate directional contributions. This yields fine-scale fidelity without sacrificing robustness to realistic measurement noise.
A.16 IMPLEMENTAION DETAILS.
All experiments are repeated three times, implemented in PyTorch and conducted on single NVIDIA RTX 3090 GPU (two RTX 3090 GPUs were used for the ERA5 experiments.). For all methods, the performance at the final epoch is recorded as the final result. For fairness in comparison, we re-ran all baseline methods under unified configuration at the same benchmarks.00 On each dataset, we use the same learning rate and number of training epochs across all methods; all models optimized with Adam. For architectures that follow the classical neural-operator design, we align the number of modes and the model width (see Table 11). U-Net-sty style baselines, such U-Net and LSM, implemented with a symmetric 6-layer decoder-decoder structure. The implementation details of RNO the various PDE benchmarks summarized in Table 11.
LNC-Duffing
RNC-Duffing
180--Dear
Valietion
E.
Opodia
Specne
Fruché
"Spoow""
Figure 11: Training loss of LNO and RNO on Duffing and Beam datasets.
It is noted that this unified configuration leads mild discrepancies relative the settings used in the original papers. For example, LNO is trained for 1000 epochs on Duffing in the original work, whereas 500 epochs here; on Beam, we adopt different choices of modes and width. These differences explain the variation in the reported numbers. For completeness, Figure 11 shows the training loss curves of LNO and RNO on the Duffing and Beam datasets under our setup.
1.1 MODEL ARCHITECTURE
Algorithmic outline: Riesz Neural Operator (RNO). Given f S4 - Cm on grid Q2 C Rd construct lifted features (0) Sr -> IRC via pointwise map Pinit that can depend both f(x)
27

---

Published as a conference paper at ICLR 2026
Table 11: RNO hyperparameters on different PDE benchmarks.
Dataset
Batch size
Epochs
Learning rate
Width
Modes
Duffing
500
0.002
16
Beam
40
500
0.002
4
16, 16
Diffusion
20
500
0.002
16, 16
Reaction-Diff.
40
300
0,002
36
8, 8
Brusselator
300
0.005
8, 8,
and the coordinates Each of the L spectral layers (in the implementation L 1) performs
(E) Fa(h (E),
(42)
&
|klla i= 1...,,,,
R((5)
|Ell2
(43)
0,
otherwise,
(E) = 00 (9) (E)
Rc(E) (1),
(44)
(E) (() R(00
15sl S m,j j 1, , ,d
(45)
2 Norm (Fz (20), w (x) = (1
(46)
(47)
(*+1()) = (x) w: (0 (x).
(48)
Here R. are Riesz multipliers implementing normalized directional derivatives: a ...... a are learnable scalar scale parameters (in the 2D implementation, Ol,, 02 correspond to aca le_x, scale_rl, scale_r2); (()(E) CCxC is complex-valued Riesz multiplier supported; 1 is 1 ... x convolution acting pointwise. A normalization layer (e.g. instance normalization) is denoted by Norm.
(e)
(2)
02
(e)
After L layers, a pointwise readout MLP produces the prediction
r(x) (w
(t) sin t,
(49)
Upred (I) W2r(z),
(50)
which matches the fcl-sin-fc2 block in the implementation.
Algorithmic outline: Spatial-Riesz Neural Operator (S-RNO). Algorithm replaces the Fourier- domain Riesz step with a filter-based approximation. For f Ωr x - Ci,, each layer con- volves current features with bank {a:)/z of fixed FIR filters whose frequency responses ap- proximate - producing directional Riesz components ri = a, h11), These define a monogenic/quaternionicc coefficient from which we extract local phase and amplitude along the distinguished axis S; optional band-pass and amplitude-weighe weighted smoothing stabilize the estimate. A scalar phase gain then rotates the Riesz coefficient, and the real part yields thase-modulated feature that is mixed across channels and fused with local residual branch via the nonlinearity, followed by pointwise readout Upred. In short, S-RNO keeps Riesz-style orientation sensitivity but performs with compact spatial filters and explicit phase handling.
Complexity comparison. Let N |62g] be the number of grid points and Cinfout the channel counts. Per layer, RNO (spectral) incurs
ON log N Cin) + (dNC + ONN Cin Coat)
RI/IRT
Riesz multipliers
pointwise mixer
28

---

Published as a conference paper at ICLR 2026
with memory ( C*) and overhead linear in d (grid-agnostir no spatial Ae factor). S-RNO (filter- based) avoids FFTs but applies d small kernels per channel:
Odd N Cm) + O(N Con).
FIR Riesz approx.
channel mixing
with memory also ( C)). Hence S-RNO is attractive on small grids or very small k; whereas for high resolutions (large N) the RNO is typically cheaper when k >log N. Both remain dimension- agnostic; RNO's cost scales with via elementwise multipliers, and S-RNO via filter applica- tions. Therefore, we adopt frequency-doman implementation while also using frequeney-domain caching to accelerate computation.
A.17 EVALUATION METRICS
We evaluate model performance using two standard metries in work: the relative fg and the mean squared error (MSH). The relative fz error measures the discrepancy normalized by the ground-truth energy and defined as
relative error- ΣM (Y; - Y:)
| - 1I
(51)
Σ (V:)
|
The MSE quantifies the average squared prediction error:
MSE
M (r.
(52)
Here Y; denotes the ground-truth value, Y, the model prediction, and M the number samples/points. For (real-valued) signal u(z) with Fourier transform ?(k)) (k), define its power spectral density (PSD) as
(k) (T))
(53)
which describes how the signal power (variance) is distributed over frequencies
A.18
SUPPLEMENTARY FOR RENCHMARKS
A. 18.1 DUFFING OSCILLATOR
A paradigmatic non-linear, damped, and driven oscillator, the Duffing system is governed by mz(t) + cáx( + k1 + (t)r((t) = f(),
(54)
where xx(t) is the displacement, (t)) and ((t) denote velocity and acceleration, and the constants m, ki, and by represent the mass, viscous damping, linear stiffness, and cubic stiffness, respectively. The term f(t) supplies the external forcing.
A. 18.2 BELER-BERNOULI BEAM
For slender Euler-Bemoulli beam subject to transverse loading, the Euler-Lagrange formalism yields
EI d'y(x,t) Ozt
PA Py((,t) f(x,t).
(55)
dt2
with y(z, () denoting the transverse deflection and f(x, 1) the applied load. Material and geometric properties enter through Young's 's modulus E, second moment of area density p. and cross-sectional area A.
A. 1813 DIFFUSION EQUATION
Pure diffusion in one dimension obeys
D B 1) Or22 Øt
0y(sr1) f(x, t),
(56)
where y(x, t) is the scalar field, D the (constant) diffusion coefficient, and ff f(x, t) a distributed source sink.
29

---

Published as a conference paper at ICLR 2026
Algorithm Riesz Neural Operator in dimensions
Require: Spatial dimension d & {1, 2,3,, grid Q4 fmn n==1 Rd
Require: Input channels Cin. hidden width C, output channels Cout
Require: Number of spectral layers (in the implementation
1: Lift (pointwise embedding):
2: for n to do
3: h(0) (x:.) Piuit (m) E IRC e.g. concatenate f(z)) and coordinates 2n, and apply
a linear layer
4: end for
5: for t= 0 L do
6:
(a) Adaptive Riesz spectral update
7:
Reshape ((9) to H IR
RRXN Nx
8:
H + Fa(Hf)
b d-D RT
9:
For each frequency index &. with wave vector {(x) RR, define
(x))
[[E()]]22 > O, i= 1, 1....,d. d.
(N) =
[((())]l2*
0,
otherwise,
10:
Combine identity and Riesz channels with learnable scales
R((c, k) - 00 (c) +
( Ra(k) (c,,
(in 2D, ag 01, az correspond to scale_x, scale_ rl, scale_r2)
11:
Initialize V with zeros (same shape as R)
for all with Kzl S my for = -....,d d do
13:
(s) + We(n) R(;,))
We(R)) E CCX,, implemented by complex weights
14:
end for
15:
b inverse RT (0 + (-rom(()) e.g. instance normalization; in code you may normalize before/after the
V - F-Fa (V)
16:
Riesz layer
17:
(b) Local residual update and fusion
for
N do
19:
(I,) + ol (1(9 (Ix)) Conv2d (width, width, 1)
1 X ... x convolution, e.g.
20: 1(++1) (xx) 2 (x + (x.)
21: end for
22: end for
23: Readout (pointwise regression):
24: for to N do
25: r(z,) + (xn))
b in the implementation o(t) sin(t), layer fcl + sin > final linear layer fc2
26: Upred (xw) + Wor(zn) E Roatt
27: end for
28: return Upred
A. 1..4 REACTION-DIFFUSIO SYSTEM
Coupling diffusion with non-linear kinetics leads to the generic reaction-diffusion form
8 pyl,t) (x,)txtt) D
dylæ, t) f(x,t),
(57)
dx2
where & sets the quadratic reaction rate while y(r, 1) and [(x 1) retain their usual interpretations.
30

---

Published as a conference paper at ICLR 2026
Algorithm 2 Spatial- Riesz Neural Operator (PRNO) with Elter-based Riesz approximation
Require: Spatial dimension d & 21,,, grid 9z (==) J
Rd
Require: Distinguished ID axis S = {1,.... l) (e.g.. time)
Require: Input channels Cin. width C, output channels Cout
Require: Temporal band (flow, fhigh) at sampling rate fa. amplification
Require: Number of layers
Require: Spatial FIR filters {a,} approximating Riesz components (i.e., DTFT Dj(w) wy/||w/| the band of interest), temporal bandpass filter T on (flows fligh), spatial Gaus- sian blur kernel
1: Input-output map:
2: Given f: S4 x S + S-pC RC, define No + Upred-
3: Lift (pointwise embedding):
4: (I s) < (f s)) RC for all (In, 8)
5: for = 0 do
6:
(a) Filter-based Riesz responses on Pp:
7:
(1n, + (x,,)
8:
R. (In, - (aj 11 (z)) for l,...,d
b d-D conv on Sr.
9:
(00) (In:8) (I( 8), Ri (In, s),... R2y (wn,8))
10:
(b) Quaternionic phase pipeline along S:
11:
/ (In, t g0000(nn s) (0)) (In,
1) (with using 5 - 1 1)
12:
(Ap(In: 8), A(In, 8)) < PHASEDIFFANDAMP(( (xn:
Ou (Zin'' f cumulative sum of Ad((n,)) over
Øy
((ou(an:))
temporal bandpass
15:
(o(-,8) f- AMPWEIGITEDBLURR(//
spatial blur on 22
(a 8) < a
mag Fm + PHASESHIFT (q(0) (Inz 8), Pa (J(()))))
18:
mmng (In, 8) t REALPART (q (Q (I (x::s)) s))
19:
(c) Nonlocal-local fusion:
(0) (In, 8) + (() (Ian
X channel mixing
w ( 8) x- (/00 (1(00
21:
22:
(E++)) (xn;)) + (=(10) (1 () 8) w 1 (xn;;) 8))
23: end for
24: Readout (pointwise regression):
25: upred 3) Dout (/ (In, E CCn
26: return Upred
A. 18.5 BRUSSELATOR REACTION-DIFFUSS SYSTEM
The two-species Brusselator, describing an autocatalytic chemical network, is written
Ou
Puu
Su
Do
+a f(t) -1 +011+++++00 vu?
0t
0x12
0n²
(58)
Ov
02
02
+ bu - vu?
Or.2
2y²
defined for (z. e [0, 112 and t e [0,20]. Here and v are concentration fields, Do and D1 their diffusion coefficients, and a, b the kinetic parameters, with ft) introducing time-dependent feed term.
A. 18. NAVIER-STOKES EQUATION
For an incompressible Newtonian fluid with velocity u(x, t) and pressure p(x, t),
0n (u. V)u - -V vAu+f.
v-u=
31

---

Published as a conference paper at ICLR 2026
where v denotes the kinematic viscosity and f body force. Upon nondimensionalliztion using the characteristic speed U and length L,
0,u + (u++(uu V/-
Af Re
V v-u
the sole control parameter is the Reynolds number Re = UL/V, which quantifies the ratio of inertial viscous effects. Higher Reynolds numbers are accompanied by increasingly complex local structures, which is one of the challenges faced by current neural-netwond methods. Accordingly, using configurations spanning range of Reynolds numbers in this experiment to validate the model's capacity to represent local nonlinearities is both reasonable and reliable. As shown in Figure 12, RNO retains excellent predictive performance despite increasing local complexity.
Figure 12: Comparison between the ground truth and RNO predictions across different Reynolds
numbers. (a) Re = 40, (b) Re 500, (c) Re 5000.
A.18.7 ERA5 DATASET
We use the fifth- sethheeenetttonn ECMWF reanalysis (ERAS) as our data source and focus on the geopo- tential height 500 hPa (7500), standard mid-tropospheric diagnostic of large-scale circulation. The native ERAS Z500 fields are regridded to regular 2.5 nattude-longite 72 x 144 global grid. We consider a continuous period from 2013 to 2023, using the years 2013-2022 for training and 2023
mesh, resulting in
year. From this time span construct more than 20,000 input-target sequences for our experiments. The processed dataset occupies approximately 2.5 GB on disk, of which the 2023 test split accounts for about 0.6 GB.
held-out
To reduce the influence of static geographic biases in the 7500 field and to emphasize dynamically relevant circulation anomalies, work with grid-point anomalies rather than raw absolute values. Specifically, for each grid point () and time t, we define
Z Zz.jtt Z
2013-2019
(59)
2013-2019 is the multi-year climatological mean where Zi.ii denotes the original Z500 value and Z
that location computed over 2013-2019. This anomaly transformation removes the reading-order stationary spatial structure (e.g., latitudinal gradients and orographic imprints), encourages the model to focus on temporal variability and circulation patterns, and avoids information leakage because the climatology computed without using the 2023 test year. Unless otherwise stated, all experiments and evaluation metrics are reported in terms of these anomaly fields, while absolute 7500 values can be recovered by adding back the corresponding climatology.
Figure 13: Solution comparison between the ground truth and predictions across ERA5.
DISCUSSION AND LIMINATION
RNO's key insight a scale-wiss wise repartitioning of the data. This reallocation expands the model's representational degrees of freedom. In neural-operator architectures. it provides principled mechanism to surface local components that would otherwise be suppressed. Although these local components may spatially overlap with global features, such overlap does not degrade performance, as
32

---

Published as a conference paper at ICLR 2026
the original signals at those locations underweighted. RNO explicitly amplifies these neglected local contributions and rebalances the component weights. From a Taylor-expansion perspective, these local terms admit further refinement, which leave to future work.
RNO further models the coupling between the data's wave content and the neural- ceural-operator architecture, yielding more complete justification of the design. Nonetheless, RNO has limitations: in some settings, control over local components lacks sufficient precision that an issue we aim to address in subsequent work.
LLM USAGE STATEMENT: LANGUAGE POLISHING ONLY
We used large language model (LLM) solely for language polishing improve clarity and control manuscript length (e.g., correcting grammar and refining phrasing for greater precision.). The LLM was not used for problem formulation, methodology or experiment design, data processing or analysis. result interpretation, or drawing conclusions, All edits suggested by the LLM were reviewed and approved by the authors.
33


Reviews:
[{"review_id": "BzEuGxM84j", "reviewer": "Reviewer_UKYE", "summary": "The paper proposes the Riesz Neural Operator, a novel architecture designed to enhance the ability of neural operators to model complex dynamic variations. The key insight is to overcome the limitation of current methods that often overlook or collapse local non-stationarity and spatial locality. The authors demonstrate that RNO achieves better prediction accuracy and generalization performance on various PDE benchmarks, including Navier-Stokes and the real-world ERA5 climate dataset.", "questions": "- What is the complexity of your Riesz transform implementation.\n\n- Can you differentiate your method from FNO + CNN?", "limitations": "", "rating": 4, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "- The central idea of explicitly embedding the Riesz transform to extract direction-aware spectral derivatives.\n\n- Results on 2D NS and the ERA5 climate dataset shows improvement on high-frequency modeling.", "weaknesses": "- Lack of complexity analysis: This new neural operator architecture needs be accompanied by an analysis of its computational overhead, especially given that Riesz transform is not a popular thing in terms of both a analysis tool and practical implementations. The paper didn't provide comparison of training time, memory consumption, or parameter count against popular baselines, which are known to have very efficient implementations. If RNO is implemented in terms of FFT, it becomes crucial to ablate the model design (see below).\n\n- Novelty of local features. The local-global split has already been applied to standard FNOs using e.g. local convolutions. Given the relation between Riesz transform and FFT, it is not clear whether the improvement comes from the representation itself, or a carefully tuned model architecture/hyper-parameter.\n\n- Motivation of derivatives. The authors motivate their architecture by emphasizing the importance of derivatives in solving PDEs. While this is a generally safe statement, the characteristics of different PDE classes makes it less convincing. There are many (spatially) localized PDEs, but there are also PDEs that are global, or even scale-local. Motivating an architecture based on a vague view of all PDEs does not convince me. Nevertheless, I acknowledge that spectral bias is a very significant deficiency of neural operators. I suggest the authors to emphasize on how RNO properly models the full spectrum, why different rank-reduction method are useful for global vs local dynamics.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "QBGMbbudj0", "reviewer": "Reviewer_F9gf", "summary": "This paper proposes the Riesz Neural Operator (RNO), which leverages the Riesz transform to enhance neural operator learning for PDEs. By mapping inputs into a Riesz domain that encodes directional and derivative information, RNO integrates both local and global spectral components to model anisotropic and high-frequency behaviors more effectively, introducing direction-aware mixing for richer feature representations. The method is supported by comprehensive experimental validation across diverse PDE settings.", "questions": "Please address the noted weaknesses, particularly by providing experimental evaluations of computational complexity and memory usage compared to baselines. Additionally, testing the model’s stability and performance on noisy data and singular solutions would greatly strengthen the validation of the proposed method’s practical robustness.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "This paper presents a novel integration of the Riesz transform to encode directional and derivative information within a neural operator framework. The design is theoretically grounded and suggests the potential for constructing neural operators that effectively capture both local and global information. The proposed architecture is sound and successfully extends the representational capacity of neural operators, and the paper is overall well written. Furthermore, it provides thorough experimental validation across diverse datasets and baselines. In particular, the studies on NS equations with varying Reynolds numbers and the ablation analyses of the proposed architecture were especially effective in verifying the model’s components and advantages.", "weaknesses": "* The paper emphasizes that using the Riesz transform enables efficient representation of higher-order derivatives. However, employing the Riesz transform is primarily a design strategy that maps inputs into a space where higher-order derivatives are representable—it does not inherently guarantee that the model will learn or reproduce higher-order derivatives more accurately. To substantiate this claim, it would be important to experimentally evaluate how well the model reconstructs derivatives of various orders (e.g., $L^2$ and $L^\\infty$ errors of $\\nabla u$, $\\nabla^2 u$, $\\nabla^3 u$).\n\n* The modular structure and the use of complex-valued or monogenic representations may introduce nontrivial computational and memory overhead. A fair comparison of training time, inference time, peak memory allocation during training, and parameter efficiency  against baselines is needed.\n\n* The proposed model may introduce potential instability in practical scenarios, as Riesz transforms can amplify high-frequency noise and encounter numerical issues near singularities. Therefore, experimental validation of stability under such conditions (i.e., non-smooth solutions and noisy data) is necessary. This is especially important because real-world data often contain noise, and many practical PDEs exhibit singularities; thus, these experiments are crucial to assess the practical robustness of the proposed method.\n\n* Tables 3 and 4 could be better positioned; their current placement does not align well with the corresponding discussion in the text.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "1Y8RCsKxm0", "reviewer": "Reviewer_Pgop", "summary": "The key idea of this paper is to embed the Riesz transform, a spectral tool representing directional derivatives, into the operator learning framework. By doing so, RNO blends global spectral features with local directional variations, addressing a noted gap in existing operator-learning methods. The paper demonstrates that RNO achieves state-of-the-art accuracy on a broad suite of PDE benchmarks and a real-world dataset.", "questions": "See cons", "limitations": "", "rating": 4, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. The use of the Riesz transform to incorporate local spatial derivatives into a neural operator is a fresh and well-motivated idea.\n2. The RNO delivers consistent and significant improvements across a wide range of PDE problems.\n3. The paper is supported with rigorous theoretical justification. (Though the reviewer does not have enough time to justify its correctness)", "weaknesses": "1. I think some of the baselines are higher than those reported in the original paper. For instance, LNO on beam & Duffing. On some datasets, this would affect the ranking of different methods. The reviewer would like some explanation on this. \n2. The experiment on Navier-Stokes is not rigorous enough, given the limited baselines.\n3. The writing of this paper could be hard to follow. For instance, some of the hypotheses proposed in the introduction are not fully justified later.\n4. Baseline, specializing in capturing high-frequency features, is also missing.\n5. The implementation code of the paper is also missing.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ZS0dx9rkrX", "reviewer": "Reviewer_xBjE", "summary": "The paper considers the problem of learning PDE solutions with neural operators, and proposes Riesz neural operator (RNO), which works on the spectral domain to improve spatial locality representation. The authors provide approximation error bounds with the use of first order Taylor expansion, and empirically show performance improvements over a couple baselines on several PDEs.", "questions": "1. PSD is used without definition, is that power spectral density?\n2. I believe the Reisz space in the paper means the Reisz field or Reisz representation in signal processing, instead of the vector lattice in functional analysis. The authors should clarify and be more accurate in wording.\n3. How are the directional derivatives effectively learned during training?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The idea of using Reisz transformation to capture local derivatives is interesting and seems beneficial in learning PDEs with highly changing solutions.\n2. The paper is in general clearly written with definitions of the core Reisz transformation explained in the Appendix.", "weaknesses": "1. There seems to be a lack of comparison (both conceptually and experimentally) with an important baseline WNO [1], which uses wavelet transformation and addresses localization in both time/space and frequency.\n\n2. A detailed training procedure is missing. A full algorithm in the appendix should help.\n\n3. The ERA5 problem seems to serve as an important realistic application of the method, but misses key details such as training/testing data size, solution visualization etc. Providing a public link does not make the paper self consistent.\n\n[1] Tripura, Tapas, and Souvik Chakraborty. \"Wavelet neural operator: a neural operator for parametric partial differential equations.\" arXiv preprint arXiv:2205.02191 (2022).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)