Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
Authors: Srinath Dama
Kevin Course
Prasanth B. Nair
Venue: ICML 2025 poster (2025)
Keywords: spatio-temporal modeling
time-series modeling
time-shift operator
Khatri-Rao neural operator
neural operator
operator learning
time-series forecasting
Primary area: general_machine_learning->sequential_network_and_time_series_modeling
TL;DR: We present a new operator-theoretic paradigm for temporal and spatio-temporal forecasting problems by learning a continuous time-shift operator.

Abstract: We present an operator-theoretic framework for temporal and spatio-temporal forecasting based on learning a *continuous time-shift operator*. Our operator learning paradigm offers a continuous relaxation of the discrete lag factor used in traditional autoregressive models, enabling the history of a system up to a given time to be mapped to its future values. We parametrize the time-shift operator using Khatri-Rao neural operators (KRNOs), a novel architecture based on non-stationary integral transforms with nearly linear computational scaling. Our framework naturally handles irregularly sampled observations and enables forecasting at super-resolution in both space and time. Extensive numerical studies across diverse temporal and spatio-temporal benchmarks demonstrate that our approach achieves state-of-the-art or competitive performance with leading methods.

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
Srinath Dama
Kevin Course'
Prasanth B. Nair
Abstract
data
We present an operator-theorent framework for temporal and spatio-tempora forecasting based learning a confinuous time-shifi operator. Our operator learning paradigm offers continuous relaxation of the discrete lag factor used in tradi- tional autoregressive models, enabling the history of a system up to a given time to be mapped to its future values. We parametrize the time-shift oper- ator using Khatri-Res Rao neural operators (KRNOs), novel architecture based non-stationary inte- gral transforms with nearly linear computational scaling. Our framework naturally handles irregu- larly sampled observations and enables forecast- ing at super-resolution both space and time. Extensive numerical studies across diverse tempo- ral and spatio-temporal benchmarks demonstrate that our approach achieves state-of-the-ar or com- petitive performance with leading methods.
forecast
time
Figure 1. The top row shows low-resolution test data. the bottom row plot high- fighrreoolution forecast. By parametrizing the time- shift operator with a Khatri- Rao neural operator we can forecast in super-resolution both space and time.
models in natural language processing tasks ( (Vaswani et al., 2017) and computer vision applications (Dosovitskiy et al., 2020), this class of models is gaining popularity in time- series forecasting (Chen et al., 2021; Nie et al., 2023; Zhou et al., 2022; Wu et al., 2023; Liu et al., 2022; 2024; Gru- ver et al., 2024; Zhang et al., 2024). In the world of spatio-temporal forecasting, Gaussian processes (Hameli- jnck et al., 2021), deep operator networks (DeepONets) (Lu et al., 2021), and neural operators (Li et al., 2020b;c;d) have emerged as cornerstones of the literature.
1. Introduction
Time series forecasting fundamental problem in ma- chine learning and statisties with applications to a broad spectrum of problems encountered in all branches of sci- ence, engineering, and finance (Roberts et al., 2013; Milani et al., 2017; Siami-Namini & Namin, 2018). At high-level, time-series problems are concerned with forecasting the fu- ture values of quantities of interest given past observations of the same or correlated quantities.
A major challenge with all autoregressive-style models is that observations are required to be provided a constant frequency at both training and inference time. This require- ment introduces number of challenges in practice. First, when observations are not provided at regular intervals, it is common practice to create a hierarchy of approximations that can negatively impact performance for reasons unre- lated the capacity of the model. Second, in online setting, this requirement necessitates creating pipeline for imputing any missing data points (due sensor error or system latency) before predictions can be made. While neu- ral ordinary differential equations based methods (Chen et al., 2018; Rubanova et al., 2019) have shown tremendous promise for learning from irregularly spaced observations, they are challenging scale and train for large-scale tem- poral and spatio-temporal datasets.
The majority of methods for time-series forecasting largely fall into the categories of autoregressive moving average models and their variants (Box & Jenkins, 1976; Girard, 2004), and deep autoregressive models with memory (El- man, 1990; Hochreiter & Schmidhuber, 1997; Salinas et al., 2020). With the tremendous success of transformer-basee
'Inttiute for Aerospace Studies, University of Toronto, ON, Canada. Correspondence to: Srinath Dama <srinath.dama @ mail.utoronto.ca,, Kevin Course <kevin.course @ www.l.lttoooooc..... Prasanth B. Nair <pras- auth.nair&e utar
In the present work, we propose casting time-series forecast- ing problems as a supervised learning problem of the contin- uous time-shift operator. In contrast to standard autoregres- sive models based on discrete- time (or discrete space-time) representation of the dynamics, the continuous time-shift
Proceedings of the 42nd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
operator maps the entire, continuous history of the dynamics over past time-window into its future values over a subse- quent time-window, Our theortort approach can be viewed as a continuous relaxation of the discrete lag factor in autoregressive models. This offers several practical ad- vantages such as the ability to learn directly from irregularly sampled observations and to forecast at super-resolution in both space and time while retaining the stability of training neural operators; see Figure I.
property: At: A A where ln 4 <
12 < t.. and (2) continuity property: EC > such that 1| 1h
1p 21 - A æ2 (a(,,))) S Cl| - وy >=)1==(())._).)>>>
for all 21, 22 E L ( 4] I):: see Appendix A for a proof.
The proposed time-shift operator formalism offers several compelling advantages for data-driven learning. First, offers robustness through leveraging a richer representation of the system's 's history, potentially mitigating uncertainty in any single observation. Second, it naturally accommodates irregularly sampled observations by treating the history of system dynamics as a function. In other words, since the time-shift operator is continuous in time, it can be learned from irregularly sampled observations (similar to neural ODEs (Chen et al., 2018) but without requiring adjoint sen- sitivity calculations), significant advantage in many prac- tical applications. Third, the continuous-ting formulation enables super-resolution forecasting in both space and time. Moreover, since the operator depends on lp and l,. this rep- resentation enables the study of the dynamics of complex systems over different time scales. In time-series forecast- ing contexts, lp and l are treated as hyperparameters that can be tuned via cross-validation or using techniques such as hypergradient descent (Chandra al., 2022).
To address the complexities of learning the time-shift oper- ator for temporal and spatio-temporal dynamical systems, we propose Khatri-Rao neural operators (KRNOs). KRNOs are a new architecture for operator learning based on non- stationary integral transforms that provides exceptional model flexibility compared to methods based on station- ary kernels (Li et al., 2020b; Kovachki al., 2023; Rah- man et al., 2023), while achieving almost linear scaling. We demonstrate the efficacy of the proposed approach on a suite of challenging test cases, including shallow water simulation (Kissas et al., 2022), a climate modeling prob- lem (Kissas et al., 2022), a set of challenging irregularly sampled time-series benchmarks, the Darts datasets (Herzen al., 2022), and the M4 dataset (Makridakis et al., 2020). In total, we evaluate KRNOs on 39 different test cases and compare performance against numerous modern approaches for temporal and spatio-temporal forecasting to demonstrate its strong generalization capabilities.
The notion of shift operators has been widely studied in functional analysis; see, for example, Marchenko (2006). Recent theoretical work (Zhen et al., 2022b;a) leveraged time-shift operators while studying the relationship between the spectra of the autocorrelation function and the infinite- dimensional Koopman operator (Koopman, 1931b) govern- ing the evolution of observables. However, to the best of our knowledge, the idea of developing operator- sperator-theoretic framework to directly learn the continuous time-shift op- erator from time-series and spatio-temporal observations remains unexplored.
2. Method
We first introduce the continuous time-shift operator for tem- poral and spatio-temporal dynamical systems. Following this, we propose Khatri-Rao neural operators for learning the time-shift operator.
2.1. The continuous time-shift operator
In the present work, we propose to parametrize the time- shift operator using a neural operator. To motivate this, consider the special case where At," (.pp f]: TR*) + (² Ly]; IR*) is a Hilbert-Sehmidt operator (Retherford, 1993). Then there exists a kernel x : [0, [) 10, [0, R satisfying the condition
Consider an ordinary differential equation (ODE) a(t) = ((=())),(() = zo. with Lipschitz continuous F : IR" -> IR" on the time interval J.. While the classical flow map propagates a single state to future times, learning such a map can be challenging in real-world applications involv- ing sparse, noisy measurements. We therefore propose alternative paradigm: the continuous time-shift operator At,",, a causal, continuous-timg time operator that maps the his- tory of z over [tp, t) to its future values over (t, tyl, where 0 < tp < t, S T, i.e., a propagator of the form
=(T) = -A ()() (5(T) (.tt=((t(.))) ()))dd)
(2)
Vr E (1, E]], where the dependence of the kernel on (t. ty) is not explicitly indicated for simplicity of no- tation. It is worth noting that even though the preceding continuous convolution integral representation holds under restrictive assumptions on the dynamics, this representation motivates approximating the time-shift operator for general nonlinear dynamical systems using deep neural operators, which involve a nested composition of integral transforms and point-wise nonlinearities.
x(7) = (14tt=)(( VT E (1,tg)-
(1)
The existence of A (Lt, R²) + [Lt follows from the Picard-Lindelöf theorem and noting that <(T) = 2(t) + Sr F(s(s))ds, T E (t, The time-
shift operator satisfies two key properties: (1) semigroup
2

---

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
We can similarly define the spatio-temporal time-shift oper- ator for scalar field f2 x [0, T] + R, where S2 C Rd-1 (d > 1) denotes a bounded Lipschitz domain. Using the non- overlapping intervltsssa defined previously, the spatio- temporal time-shift operator can be defined as
Table I. Comparison of Khatri-Raa Neural Operator (KRNO), Graph Neural Operator (GNO) (Li et al., 2020c), Multipole Graph Neural Operator (MGNO) (Li et al., 2020d), and Fourier Neural Operator (FNO) (Li et al., 2020b), for computing kernel inte- gral transforms, as compiled by (Kovachki et al., 2023). Here N' << N is constant used to control the variance of the integral transform approximation. Ours is the only approach which allows for exact, non-stationary kernel evaluations while achieving almost linear computational cost.
u(x, 7) = (A u)(, E I TE (t,ty)
(3)
Under the assumption that u lies in the separable Hilbert space U(2 [0, R),, and assuming the spatio-temporal time-shift operator is Hilbert-Schmide operator mapping from U(2 x [tp, t); to U(S x (t,tyl:R), we have the following integral representation
Method GNO MGNO FNO KRNO
Time O(NN')
Non-stationary
Exact kernel'
O(N) ON N) X
((NI+1/4)
V
u(I, T) = (dttu)(a x TEQ, TE (c,4g)
= La.. [[(((a,r)) u((o,T), (y, s)dyds, 5,
these advantages. We will show later that KRNOs provide state-of-the-an art performance across number of benchmarks while inheriting the benefits of neural operators such as he- ing discretizationninvaaian and enabling super-resolution in forecasts (Li et al., 2020b).
(4)
where i: : x2 [0, x 22 x [0, T] + R is a square-integrable kernel. The preceding representation in terms of an integral transform motivates the application of deep neural oper- ators approximate the time-shift operator of complex spatio-temporal dynamical systems. Furthermore, universal approximation results for neural operators (Kovachki et al., 2023) guarantee that, under suitable regularity assumptions, the time-shift operator can be well approximated.
Neural operators Neural operators (NOs) are an expres- sive class of models for approximating maps between func- tion spaces. In contrast to standard multi-layer perceptrons, which are defined by an alternating series of affine maps and nonlinear activations, NOs are defined by an alternating se- ries of linear, kernel integral transforms and nonlinear activa- tions. For simplicity of exposition, consider integral trans- form layer (Li et al., 2020b; Kovachki et al., 2023) that maps the input spatio-temporal vector field Uy 12 x [0, r] + R" to 02+1 Ω x [0, r] + IR,, defined below
worth noting that Li et al. (2020b) considered an op- erator learning test problem where wwoddimesiionl flow field over the time-interval [0, 10] is mapped the interval (10, 7] (for a fixed 7). This was tackled using an autoregres- sive Fourier neural operator (FNO) model and a 3D FNO model, with the latter model making predictions over the entire spatio-temporal domain of interest. The time-shift op- erator learning formalism presented here allows us to view the test-case involving FNO-3D in (Li et al., 2020b) as a special case of the general setting considered here with stationary-kerne' based neural operator parametrization of the time-shift operator and fixed values of (lp. 'f). In the next section, we present a new architecture for parametriz- ing the time-shift operator that enables over an order of magnitude reduction the number of parameters compared to FNO, while achieving superior accuracy.
Dt41(t, 2) (t,))
1 l (. Dueee o''dt ds' Wer x) +
(5)
where K R x xxxxxXxi + RAxP is a matrix-valued kernel, W E RT*P is a weight matrix, and b R29 is bias vector. It is also common to prepend and append the preceding layer by a series of point-wise lifting and projection layers (Kovachki et al., 2023). Note that in (5), the inputs and outputs are assumed to be defined over the same opatio-tempore domain for simplicity - we will later consider the general case when the input and output domains are different.
2.2. Khatri-Rao neural operators (KRNOs)
We now introduce KRNOs, a new operator learning architec- ture based on non-stationary integral transforms, to approxi- mate the time-shift operator of temporal and spatio-temporal dynamical systems. KRNOs offer expressive parametriza- tion of operators using non-stationary integral transform layers which (i) do not require any approximation of the kernel and (ii) scale almost linearly in the number of quadra- ture nodes. As far as we are aware, ours is the only ap- proach for parametrizing neural operators which combines
Rather than computing the integral transforms exactly, NOs propagate evaluations of the intermediate functions at set of quadrature nodes through the network. Let X = (141, ), (12, r (tN, xN E Nxd where
Wee use the term "Exact Kernel" to indicate that the only source of approximation error stems from the quadrature scheme; the non- stationary kernel itself is evaluated without further approximation.
3

---

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
4 I and 2y E Rd-1, N.....N denote quadra- ture nodes in the d-dimensional spatio-temporal domain (N = n², where n is the number of quadrature nodes per dimension), and let w E R denote the vector of quadrature weights. As we will show, while computing the point-wise transformation defined by the weights W and b scales as (, the primary computational bottleneck in comput- ing the output from a kernel integral transform layer arises from approximating the integral over the domain of the input function.
context of Gaussian process (GP) regression build new positive definite kernels and scale GP regression on prod- uct grids (Saatçi, 2011; Wilson et al., 2014). For example, the squared exponential kernel, the Matérn class of kernels, and the spectral mixture product kernel all decompose as a product of the form in Equation (7).
Proposition 1. the quadrature nodes lie on a product grid, X = x(1) x ((d1)),,,,,, where i E R" and xf8) E R denote the quadrature nodes along the time dimension and the ith dimension of the spatial coordinate respectively, and the kernel function has component wise product struc- ture of the form given in Equation (7), then the kernel func- tion evaluated at the quadrature nodes inherits the Khatri- Rao product structure,
In the discussion on computational complexity that follows we omit writing the dependence on and p since these are architectural considerations and the required value for N will be dependent on the complexity of the input func- tion. Letting (x) Xx be the /th layer evaluated at the quadrature nodes, the kernel integral transform can be approximated as
xXx = (1, * * (x(t-))))
ad(-11)) (x(t-)))))))))))
(8)
where 200 (:,-) Rgaxpn is a block-partitioned matrix where block jk is the jkth output from the component kernel x( evaluated on the outer product of the quadrature nodes along the 1 dimension.
La mx(x M,,, (t', xt Nue(t', (t) d'')d'' dx'
RC xX, X) (ec(ddig((w))dd(((
(6)
where vec: NxP + R NP creates a vector from a matrix by stacking columns, diag. RN + N converts a vector into a diagonal matrix, H5 xxX, x'')) E RN represents the kernel evaluated between all the quadrature nodes X in the output domain and single node {t','') x²) } in the input domain. Meanwhile, x) E R.NgxNP is the kernel evaluated between all the quadrature nodes X in both the output and input domains.
A proof for Proposition can be found in Appendix B. A practical consequence of this result is that the computational complexity associated with computing kernel integral trans- forms can be reduced from O(N2) to ((N1+1/d) with O(N2/d + N) memory; see Appendix D for detailed discussion. These complexity estimates have a linear depen- dence on pq, which are architectural parameters common to all NO approaches. Table 1 provides a comparison of KRNOs to other NOs in the literature. To reiterate what was mentioned previously, ours is the only approach which achieves almost linear cost while enabling non-stationary in- tegral transforms without having to approximate the kernel function.
Naively evaluating this kernel integral transform scales as ((2) which is prohibitively expensive for even a modest number of quadrature nodes. In light of this computational challenge, number of approaches have been developed including Graph Neural Operators (Li et al., 2020c), Multi- pole Graph Neural Operators (Li et al., 2020d), and Fourier Neural Operators (FNOs) (Li et al., 2020b; Rahman et al., 2023). As we will show, our approach is the only one which scales almost linearly in the number of quadrature nodes while enabling non- statstttinaay integral transforms with exact kernel evaluations.
In Appendix C, we present a generalization of Proposition I showing that the same computational advantage also holds when the input and output domains are different with dif- ferent resolution quadrature grids. This generalized result enables the design of non-stationary integral transform lay- ers that can be viewed continuous analog of upsampling and downsampling techniques commonly used in convolu- tional neural networks. In addition, this generalization can be leveraged to significantly reduce computational cost and memory requirements in applications with high spatial or temporal resolutions by using lower-resolution quadrature grids within the intermediate kernel integral layers. We provide detailed numerical studies comparing the training and inference costs of KRNO and FNO across different spatial resolutions in Tables 9,10, as well as Figure 7 in Appendix G.
Khatri-Rao product structure In order to achieve almost linear scaling in N without having to approximate the kernel function, we assume that the kernel function decomposes as an element-wise product as follows:
(((t, x), (t',')) = x (t,t)o
(o
lo]k-)) (7)
where x1)) : R x R + R for i = l,...,d, @ denotes the element-wise product, and [z], E R indicates the th element of x. While this assumption may appear limit- ing at first glance, it has been applied extensively in the
In this work, we parametrize each component kernel, x00 RR. x RR. - Raxp by neural network; see Appendix
4

---

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
Table 2. Performance comparison of different NO methods on Darey-flow and hyper-elastic problems. Results with ()', (ff are from Lu et al. (2022) and Li et al. (2023), respectively.
performance statisties of the proposed approach across all test cases are presented in Appendix Table 8. The code- base used to generate the results is available at https: / //gittu..... con/ ht////sinnttaaaaa Shi ft ingTime.
Method
L relative emror
Dlrey-flow
Elasticity
3.1. Spatio forecasting problems
FNO
1.19 = 0.05%!
5.08%t
DeepONet
1.36 + 0. 0.12%
9.65%f
For spatio-temporal problems, consider shallow water simulation (Kissas et al., 2022), and climate modeling dataset (Kissas et al., 2022). For evaluation, we use the L relative error metric, L relative error = ||u(,,t) u(-, ())) /l ( where u(-.t) and a(, () are true and predicted spatial fields at time
KRNO (ours)
0.96 + 0.04%
4.66 + 0.09%
for details. To illustrate the efficacy of KRNO, we first applied our method to the two-dimensional Darcy-flow and hyper-elastic benchmark problems from Lu et al. (2022) and Li et al. (2023), respectively. In these two problems, the goal is to learn mapping between input and output fields over a two-dimensional spatial domain. Figure 11 in Appendix illustrates the predictions from KRNO for these two problems; see Appendix H.5 for additional details. It can be noted from Table 2 that KRNO provides improved performance over FNO (Li et al., 2020b) and DeepONet (Lu et al., 2021) on both problems. We will later benchmark the performance of KRNO on learning the time-shift operator across a variety of challenging datasets to demonstrate that our approach often provides competitive performance when compared with SOTA methods.
Shallow water example Here, the objective is learn spatio-temporal operator that is capable of predicting three field variables (fluid column height p. velocity in the Ti dirdirecto u, and velocity in the r2-direction v) over a future prediction window (t, using historical data from a look-back window [4. i.e., ((2 x 4, R+) + U(Q x (t. ty], I)), where 52 := (0, 1) (0, denotes the spatial domain. It is worth noting that the spatio-temporal forecasting problem statement considered here is signifi- cantly more challenging than the standard benchmark in prior work, which involves mapping the initial condition to the solution at a single future time (Kissas et al., 2022).
The dataset used for this problem is taken from Kissas et al. (2022), which includes simulated data generated on 32 x 32 spatial grid over the time window (0, 1) and collected at every 0.01 seconds. The training and testing datasets each consist of 1000 simulations with different ini- tial conditions. Both the look-back and prediction window period are set to 0.05 seconds. For evaluation on testing data, we use the three field variables from the first 0.05 seconds window to recursively predict their evolution until 0.6 seconds. As baseline method, we consider FNO-3D model (Kovachki et al., 2023) and attention based neural operator LOCA (Kissas et al., 2022) to approximate the time-shift operator alongside the proposed KRNO method.
3. Numerical studies
In addition to the two spatial modeling problems (Darcy- flow and hyper-clasticity), we evaluated the performance of the proposed method on a suite of challenging temporal and spatio-tempora forecasting problems. These datasets included two spatio-temporal datasets, and 21 diverse time- series datasets (5 multivariate irregular time-series datasets, 8 univariate time-series from the Darts collection, 6 datasets corresponding to different seasonalities from the M4 compe- tition, one multivariate time series corresponding to trading prices of 14 cryptocurrencies, and a bivariate time series containing the positions of NBA basketball players). In total, this amounts to 39 test cases or error metries. Across these cases, we compare against numerous modern approaches for temporal and spatio-temporal forecasting problems.
Table 3 and Figure 2 compares the relative L3 error (aver- aged across 1000 test simulations) for the three field vari- ables when training is conducted for 100 epochs. The results indicate that KRNO delivers superior performance relative to FNO-3D and LOCA. In addition, it is worth noting that KRNO only uses 6% of the parameters required by FNO-3D (see Table 3). Predictions from KRNO for a test simulation are shown in Figure 3. Additional numerical results for this test case can be found in Appendix H.6. Climate modeling example In this experiment, we consider nettool time-shift
For all the problems, we first convert the datasets into nec- essary input and output sequence pairs for learning the time-shift operator; see Appendix F for details. The de- fault KRNO architecture used in all experiments has 3 kernel integral layers with 20 channels each, lifting and projection layers are parametrized by MLPs with one hid- den layer containing 128 hidden units, and the kernels in the integral layers are parametrized by MLPs with 3 hid- den layers. The component-we wise kernel function used in KRNO is parameterized by a neural network with three hid- den layers (see Appendix E). Additional results and experi- ment setup details are provided in Appendix H. Aggregated
the problem of approximating operator that maps the surface air temperature and surface air pressure, i.c., u(2 x [t,, t]: R²²) + U(2 x (t, (t,t where Ω = 9-0 90] x [o. 360] denotes the spa-
5

---

Shifting Time: Time-series Forecasting with Khatri-Ran Neural Operators
Table Comparison of the average relative Lt errors the shal- low water problem for the three field variables.
Method
Parameters
L relative enror
p
u
FNO-3D LOCA KRNO
16,62,60 895
0.00211 0.00314
0.02606
0.02637
94, 17.77,22 220 146, 159
0.15221
0.14999
0.00145
0.01497
0.01459
(ww)
pel 12 erno
KRNO LOGA
18-0
1/
(see)
Top figure shows the predictions (Ø, û, é) for the three field variables along with the true fields (p, u, e) as a function of time for a test simulation using KRNO trained for 100 epochs. The bottom figure shows the error bars representing the L relative errors for three field variables across the 1000 test simulations, with the shaded region indicating 1:1 standard deviation. Additional error plots for the three fields are shown in Figure 13 in the Appendix.
Figure 3. Shallow water problem:
10"2
(50c)
Figure 2. Comparison of the average relative L errors func- tion of time for the three field variables (across the 1000 test simulations) obtained using KRNO, FNO-3D and LOCA models.
and Player Trajectory datasets. ?
tial domain defined in terms of latitude and longi- tude. The dataset is taken from Kissas et al. (2022) which is based on the Physical Sciences Labora- tory meteorological data (Kalnay et al., 1996); see https:: //ppl: noaa. gov ncep. reana yysis. surface. html. The training data consists of daily temperature and pressure from 2000 to 2005 (1825 days) over a 72 x 72 spatial grid. The test data contains observations from the years 2005 to 2010 on the same grid. The KRNO operator is trained on temperature and pressure data from a 7-day look-back window, with matching 7-day prediction window. For the evaluation on testing data, used data from the last week of the previous year and recursively predicted the temperature and pressure fields for the whole year. This is repeated for each year in the testing set. Representative predictions for pressure and temperature are shown in Figure 4 along with the corresponding relative 1,2 error.. can be seen that the proposed time-shift operator learning approach performs remarkably well for this dataset.
Distribution shift presents key challenge in temporal fore- casting (Kouw & Loog, 2018; Wang et al., 2021; Kuznetsov & Mohri, 2020). A common practice to tackle this challenge is to use preprocessing strategies, which involve removing known trends and seasonality from the data. To handle distri- bution shifts in some datasets, we use the reversible instance normalization (ReVIN) approach (Kim al., 2021) nor- malize each input sequence and denormalize the model output predictions.
a
3.2.1. IRREGULARLY SAMPLED TIME-SERIES
We use datasets from diverse domains, including healthcare (MIMIC), climate science (USHCN), biomechanics (Hu- man Activity), and physical simulation (MuJoCo); see Ap- pendix 1.0.1 for details. Notably, the MIMIC, USHCN and Human Activity datasets exhibit significant temporal irregu- larities along with missing values for the states. To handle missing data, we augment KRNOs with one-dimensional convolutional network (CNN) feature extractor module, whose parameters are jointly learned; additional details on the KRNO architecture used are provided in Appendix H.O.1 for details.
3.2. Temporal forecasting problems
We evaluate the performance of our approach on both regu- larly and irregularly sampled temporal forecasting problems. For irregular time-series, use the MuJoCo, MIMIC-III, USHCN, and Human Activity benchmarks, along with synthetic 2D spiral dataset (details in Appendix H.I). For regularly sampled time-series, we use datasets from the Darts (Herzen et al., 2022) and M4 (Makridakis et al., 2020) collections, in addition to the Crypto (Ticchi al., 2021)
Table 4 compares the performance of KRNO with latest SOTA methods such as T-PATCHGNN (Zhang et al., 2024). CRU (Schirmer et al., 2022), Neural Flow (Biloš et al., 2021), and Latent-ODEs (Rubanova et al., 2019). The re- sults highlight the strong performance of KRNO. Specifi- cally, KRNO achieves the best MSE on both the MIMIC
6

---

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
Table Comparison of KRNO performance with baseline models on irregular time-series datasets. Test MSE and MAF (mean + std) are shown for each dataset. Top three best results are highlighted in bold, underline, and italic, respectively. Results for baseline models are sourced from Zhang et al. (2024), Table
Algorithm
MIMIC
USHCN
Human Activity
M
MAEx10
MSEX
MAEX10-1
MSEx10-3
MAEx10-2
DLincar
4.90 0.00
16.29 0.05
6,21 = 0.00
3.88 = 0.02
4.03 + 0.01
4.21 = 0.01
TimesNet
5.88 0.08
13.62 =002=00 0.07
5.58 + 0.05
3.60=
0.02
3,56 = 0.02
PutchTST
3.78 0.03
5.75
=.02
4.14
4.80 0.09
Crossforner
2.65 0.10
9.26
3.27 5.04
4,29
4.89 1 0.17
Graph Wavenet
2.93 0.09
1..5
3.16 0.04
0.09
3.40 0.05
MTCNN
2.71=
9.55
5.39
3.34 0.02
0.03
3.53 0.03
StemGNN
1.73 0.02
0.09 0.11
3,40 0.09
0.37
6
CrossGNN
2,95
10.82
5.66 0.04
3.53 0.05
0.10
3.48 0.08
FOURIER/GNN
2.55
10.22 0.08
5.82 + 0.06
3.62 0.07
2.99 0.02
3.42 0.02
GRU-D
1.6 0.03
7.53 + 0.09
5.54 + 0.38
3.40 = 0.28
2.94 0.05
3.51 = 0.06
ScFT
187 0.01
7.84:0.00
5.80 + 0.19
3.70 0.11
12.20 - 0.17
8.43 = 0.05
RainDrop
1.99 0.03
0.07
5.78 + 0.22
3.67 = 0.17
14,92
9.45 = 0.05
Warpformer
1.73 0.04
7.58
5.25 + 0.05
3.23 = 0.05
2.79 1.79+0.00 0.04
3.39 + 0.03
mTANN
1.85 0.06
+.73
5.33 + 0.05
3.26 + 0.10
3.22 +022+0000 0.07
3.81 + 0.07
Latent-ODE
1.89 0.19
0:51:
5.62 1 0.03
3.60 0.12
3.34 0.11
3.94 0.12
CRU
1.97 0.02
7.93
6.09 0.17
3.54 0.18
6.97 0.78
6.30 0.47
Neural Flow
1.87
8.03 0.19
5.35 + 0.05
3.25 0.05
4.05 0.13
4.36
T-PatchGNN
1.69
7.22 + 0.09
5.00 + 0.04
3.08 0.04
2.66 0.03
3.15 0.02
KRNO
1.57 + 0.02
7.43 0.06
4.95 + 0.08
3.06
0.08
2.85 0.03
3.46 0.02
Prodicted p KPs
Tre KPs
Al/M
data. Numerical studies on a synthetic 2D spiral trajectory prediction problem are presented in Appendix H.1.
day 1.00
MuJoCo The MuJoCo dataset used in this study is based simulations of the 'Hopper' model using the DeepMind Control Suite (Tassa et al., 2018). Each trajectory consists of 14-dimensional state vector sampled at 100 evenly-spaced time points. These trajectories are generated from random initial conditions, allowing the deterministic dynamics of the Hopper model to evolve. We used 100 such trajectories from the dataset provided by Rubanova et al. (2019).
200
Proticted in K
Teue 2' la K
day 120
5
The forecasting objective was to predict the values for the next 10 time steps using observations from the preceding 50 time steps. To simulate irregular sampling, we adopted the preprocessing strategy of Zhang et al. (2024), randomly re- moving 30%, 50%, 70% of samples from each trajectory. Table 5 compares KRNO's performance with state-of- the- art (SOTA) methods such as Neural SDEs (Li et al., 2020a; Oh et al., 2024). It can be seen that KRNO outperforms Neural SDES in all dropout settings. It is worth noting that KRNO achieves an error reduction of 46% on regular data and 38% error reduction on the irregularly sampled versions.
100 150 200 220
200
32D
Figure 4. Climate modeling problem: Top two figures show the predicted surface pressure and temperature fields using KRNO model along with the true fields for single day in forecasted year. Bottom figure shows the error bars representing the L relative errors for the five years in test data, with the shaded region indicating 11 standard deviation.
and USHCN datasets, and also secures the best MAE on USHCN. On the MIMIC dataset, KRNO the second- best performing model in terms of MAE, closely following T- PatchGNN. For the Human Activity dataset, T-PatchGNNN shows the strongest performance, with KRNO ranking third in MSE. Overall, KRNO consistently places among the top- performing methods, demonstrating its robustness and effec- tiveness for forecasting with complex, irregularly sampled
3.2.2. REGULARLY SAMPLED TIME SERIES
Darts benchmarks We consider 8 univariate time-series datasets from Darts (Herzen et al., 2022). We compare the performance of the proposed time-shift operator with mod- els such as ARIMA (Box & Jenkins, 1976), neural networks- based models (TCN (Lea et al., 2016), N-BEATS (Oreshkin et al., 2020), N-HITS (Challu al., 2023)). SM-GP (Wilson

---

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
Table Comparison of forecasting performance on the MuloCo dataset under different percentages of missing observations. Results for baseline methods are sourced from Oh et al. (2024), Table 12.
Test MSE
Methods
Regular
30% dropped
50% dropped
70% dropped
GRU-At
0.223 + 0.020
0.198 + 0.036
0.193 0.015
0.196 0.028
GRU-D
0.578 + 0.042
0.608 + 0.032
0.587 0.039
0.579 0.052
GRU-ODE
0.856 + 0.016
0.857 + 0.015
0.832 + 0.015
0.861 + 0.015
ODE-RNN
0.328 + 0.225
0.274 + 0.213
0.237 0.110
0.267 + 0.217
Latent-ODE
0.029 = 0.011
0.056 + 0.001
0.055 0.004
0.058 = 0.003
Augmented-ODE
0.055 = 0.004
0.056 + 0.004
0.057 0.005
0.057 = 0.005
ACE-NODE
0.039 + 0,003
0.053 + 0.007
0.053 0.005
0.052 = 0.006
NCDE
0.028 = 0.002
0.027 + 0.000
0.027
0.026 0.001
ANCDE
0.026 = 0.001
0.025 + 0.001
0.025 0.001
0.024 = 0.001
EXIT
0.026 + 0.000
0.025 + 0.004
0.000
0.026 + 0.001
LEAP
0.022 + 0.002
0.022 + 0.001
0.022 0.002
0.022 + 0.001
Neural SDE
0.028 + 0.004
0.029 + 0.001
0.029 0.001
0.027 = 0.000
Neural LSDE
0.013 + 0.000
0.014 0.001
0.114
0.013 = 0.001
Neural LNSDE
0.012 + 0.001
0.014 + 0.001
0.014 +
0.014 0.000
Neural GSDE
0.013 + 0.001
0.013 + 0.001
0.013 0.000
0.014 = 0.000
KRNO
0.007 1 0.002
0.008 + 0.002
0.0114 0.004
0.0115 0.002
& Adams, 2013) and LLMTIME (Gruver et al., 2024). We used normalized mean absolute error (NMAEX22) as the evaluation metrie. Figure 5 shows the testing errors from KRNO in comparison to other baseline methods presented in (Gruver et al., 2024). The time-shift operator is among the top 3 performing methods on 5 out of 8 datasets in the Darts collection: see Table in Appendix for details.
Table 6. Comparison of sMAPE from KRNO method with other
baseline methods for M4. Results with (.) were taken from Wang et al. (2023).
Method
Quarterly
Weekly
Daily
Montero et al. (2020)
9.733
7.625*
3.097t
Smyl (2020)
9.679
7.817t
3.170
Nbeats-I+G
9.212
KNF
10.008
7,254'
2.990'
Darts
KRNO (ours)
10.503
6,934
3.086
Geometric Mean of NMAE
0.10
models are evaluated using the symmetric mean absolute percentage error (sMAPE) metric used in the M4 compe- tition. A comparison of KRNO performance on M4 data is presented in Table (for full results see Table 13 in Ap- pendix). We observe that KRNO among the top two methods on datasets such as M4-Weekly and M4-Daily, where seasonality trends are not present (Wang et al., 2023).
0.05
0.00
N-BEATS
SM-GP
TCN
N-HTSS
L-aMA-2
GPT-3
ARIMA
TS-KRNO
Crypto and Player Trajectory datasets The Crypto (Tic- chi et al., 2021) dataset is multivariate time series con- taining eight features corresponding to trading prices of 14 cryptocurrencies. The objective is to forecast the returns for all 14 cryptocurrencies. The Player Trajectory dataset is a bivariate time series containing the positions of NBA basketball players. The goal here forecast the positions of the players. For both datasets, we utilized the same train- ing, validation, and test data used by Wang et al. (2023). We employed weighted RMSE (Ticchi et al., 2021) for the Crypto data and RMSE for the Player Trajectory data for evaluation. Our method is compared with KNF and other latest methods such as FedFormer (Zhou et al., 2022). Long Expressive Memory (LEM) (Rusch et al., 2022), Variational
Figure 5. Comparison of geometric mean of normalized MAE on univariate Darts datasets for various methods.
M4 benchmarks The M4 dataset (Makridakis et al., 2020) a collection of 100,000 univariate time series from diverse domains such as finance and demographics. This collection comprises six datasets corresponding to different seasonal- ities, varying from hourly to yearly. On this challenging dataset, the top winning methods in the M4 competition, Smyl (2020) and Montero et al. (2020), Koopman Neural Forecaster (KNF) (Wang et al., 2023), and Nbeats-I+G (Ore- shkin et al., 2020) are considered as baselines. All the
8

---

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
Beam Search (VBS) (Li et al., 2021), Multilayer Percep- tron (MLP) (Faloutsos et al., 2018), and Vector ARIMA (VARIMA) (Stock & Watson, 2001). Table 7 compares KRNO with these baseline methods. KRNO is the second- best method after KNF on Crypto and Player Trajectory datasets (for full results see Table 15 in Appendix H.4).
used the FNO learn map for the solutions of the Navier- Stokes equations from the time interval [0, 10) to (10, 50) and discussed the ability of this strategy to provide super- resolution in both space and time. This specific application can be viewed as an instance of our time-shift operator for- malism with fixed time windows. Our framework is more general, employing learnable non-stationary kernels (via KRNOs) to parametrize the time-shift operator and treating the time window length as hyperparameters. Moreover, as demonstrated in our numerical studies (Section 3), KRNOs are considerably more parameter-efficient than FNOs while often achieving superior accuracy.
Table Comparison of RMSE from KRNO method with other baseline methods on Crypto and Player Trajectory datasets.
Crypto
Basketball (RMSE)
(Weighted RMSE 10 3)
VARIMA MLP TE
8.76
1.26
7.85
1..32
Transformer- models have gained significant popu- larity in time-series forecasting (Zhang et al., 2024; Nie et al., 2023; Chen et al., 2021; Zhou et al., 2022: Wu et al., 2023; Liu et al., 2022; 2024). It is worth noting that most transformer architectures for time-series forecasting assume regularly sampled observations. A recent exception is T-- PatchGNN (Zhang et al., 2024), designed for irregularly sampled data, against which we compare in Section 3.
7.00
1.48
RF+TB
7.84
2.40=mmm
FedFormer
7.462004
1.29
LEM
7.02+0004
1.42+0.02
VBS
19. 19.52=0.00
5.60
KNF
6,91+t01
1.16=0.0
KRNO
6,95
1.25
4. Related work
5. Concluding remarks
The continuous time-shift operator introduced in this work
is distinct from the Koopman operator (Koopman, 1931a), an infnite-dimmesiinnll linear operator acting on a space of observables, which has been extensively studied in dynam- ical systems and machine learning (Mezić, 2021; Brunton et al., 2022; Wang et al., 2023; Liu et al., 2023). The time- shift operator corresponding to a set of sufficiently smooth observables can be viewed as a continuous extension of the Koopman operator (Zhen et al., 2022b). This theoretical connection deserves further study.
We introduced novel operator theorretc framework for tem- poral and spatioltemporal forecasting based on learning a continuous time-shift operator. This paradigm offers prin- cipled and robust approach to map the continuous history of a system to its future evolution, providing a continuous relaxation of traditional discrete-lag autoregressive models. To parameterize this operator, we proposed Khatri-Rao Neu- ral Operators (KRNOs), new architecture that leverages non-stationary integral transforms. KRNOs achieve nearly linear computational scaling with the number of quadrature points while allowing for exact kernel evaluations, a unique combination of efficiency and expressiveness.
Similar neural ODE based methods (Chen et al., 2018; Rubanova et al., 2019; Kidger et al., 2020) and neural SDEs (Li et al., 2020a; Oh et al., 2024), our framework formulates time-series forecasting in a continuous setting. Our approach offers significant computational advantages during training, as it obviates the need for adjoint sensi- tivity methods to compute gradients of the loss function (similar to simulator-free approaches for learning SDEs, see, for example, Course & Nair (2023)). Furthermore, unlike typical NODE and SDE formulations that focus on tempo- ral dynamics, our operator- theooreti framework naturally extends to spatio-temporal problems, inherently enabling super-resolution in both space and time.
The practical advantages of our approach include the ability to handle irregularly sampled observations, perform super- resolution forecasting in both space and time, and achieve significant parameter efficiency. Extensive numerical stud- ies across 39 diverse test cases, encompassing physical simulations and challenging real-world time-series, demon- strate the efficacy and scalability of our method. KRNOs consistently achieve state-of-the-a or highly competitive performance, outperforming leading methods on several benchmarks, particularly for irregularly sampled data and complex spatio-temporal dynamics.
Neural operators, such as the Fourier Neural Operator (FNO) (Li al., 2020b) and its U-Net inspired variant U-shaped Neural Operator (U-NO) (Rahman et al., 2023), have ad- vanced the learning of mappings between functions defined over spatio-temporal domains. Liu-Schiaffini et al. (2023) proposed
We envision several exciting research directions stemming from this work: further thcoretical analysis of the time-shift operator itself, including its connections to other operator- theoretic constructs; development of advanced KRNO vari- ants or alternative parameterizations; and the application of this paradigm to new domains.
recurrent neural operator (RNO) ar- chitecture to address forecasting non-stations non-stationary dynamics on function spaces. As previously noted, Li et al. (2020b)

---

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
Acknowledgements
Chen, R. T. Q., Rubanova, Y., Bettencourt, J., and Duve- naud, D. Neural Ordinary Differential Equations. In Advances in neural information processing systems, pp. 6571-6588 2018.
This work is supported by a Natural Sciences and Engineer- ing Research Council of Canada Discovery Grant,
Impact Statement
Course, and Nair, B. State estimation of a physical system with unknown governing equations. Nature, 622 (7982):261220) October 2023.
Work presented in this paper is mostly theoretical and com- putational in nature. At this stage, we don't t see any potential societal impacts.
De Brouwer, E., Simm, J., Arany, A., and Moreau, Gru-ode-bayes: Continuous modeling sporadically- observed time series. Advances neural information processing systems, 32, 2019a.
References
Ames, W. F. and Pachpatte, B. Inequalities for differential and integral equations, volume 197. Elsevier, 1997.
De Brouwer, E., Simm, J., Arany, A., and Moreau, Y. Gru-ode-bayes: Continuous modeling of sporadically- observed time series. Advances in neural information processing systems, 32, 2019b.
Ba. L., Kiros, R., and Hinton, G. Layer normalization. 1XXiv:607.064500 2016.
Biloš, M., Sommer, J., Rangapuram, S. S., Januschowski, T., and Günnemann, S. Neural flows: Efficient alternative to neural odes. Advances in neural information processing systems, 34:21325-213377 2021.
Dosovitskiy, A., Beyer, I.., Kolesnikov, A., Weissenborn, D., Zhai, X., Unterthiner, T., Dehghani, M., Minderer, M., Heigold, G., Gelly, S., et al. An image is worth 16x16 words: Transformers for image recognition at scale. In International Conference on Learning Representations,
Box, G. E. and Jenkins, G. M. Time series analysis forecasting and control. Holden-Day series in time series analysis. Holden -ay, San Francisco, rev. ed. edition, 1976. ISBN 0-8162-1100 11043..
2020.
Dupont, E., Doucet, A., and Teh, Y. W. Augmented neu- ral odes. Advances in neural information processing systems, 32, 2019.
Brunton, S. L., Budišić, M., Kaiser, E., and Kutz, J. N. Modern koopman theory for dynamical systems. SIAM Review, 64(2):229-35) -340, 2022.
Elfwing. S., Uchibe, E., and Doya, K. Sigmoid-weighte linear units for neural network function approximation in reinforcement learning. Neural networks, 107:3-11, 2018.
Cao, D., Wang, Y., Duan, J., Zhang, C., Zhu, X., Huang, C., Tong, Y., Xu, B., Bai, J.. Tong, et al. Spectral tem- poral graph neural network for multivariate time-series forecasting, Advances in neural information processing systems, 33:17766-17778, 2020.
Elman, L. Finding structure in time. Cognitive Science, 14(2):179-111, April 1990.
Faloutsos, C., Gasthaus, J., Januschowski, T., and Wang, Y. Forecasting big time series: old and new. Proceedings of the VLDB Endowment, 11(12):210-2211 2018.
Challu, C., Olivares, K. G., Oreshkin, B. N., Ramirez, F. G., Canseco, M. M., and Dubrawski, A. Nhits: Neu- ral hierarchical interpolation for time series forecasting. In Proceedings of the AAAI Conference Artificial Intelligence, 2023.
Girard, A. Approximate methods for propagation of un- certainty with Gaussian process models. PhD Thesis, University of Glasgow, 2004.
Chandra, K., Xie, A., Ragan-Kelley, J., and Meijer, E. Gra- The ultimate optimizer. In Advances in Neural Information Processing Systems, 2022.
Gruver, N., Finzi, M., Qiu, S., and Wilson, A. G. Large language models zero-shot time series forecasters. Advances in Neural Information Processing Systems, 36, 2024.
dient descent:
Che, Z., Purushotham, S., Cho, K., Sontag, D., and Liu, Y. Recurrent neural networks for multivariate time series with missing values. Scientific reports, 8(1):6085, 2018.
Hamelijnck, O., Wilkinson, W., Loppi, N., Solin, A., and Damoulas, T. Spatio-temporal variational Gaussian processes. Advances in Neural Information Processing Systems, 34:23621 2021.
Chen, M., Peng, H., Fu, J., and Ling, H. Autoformer: Searching transformers for visual recognition. In Proceedings of the IEEE/CVF international conference 2021.
Herzen, J., Lässig, F., Piazzetta, S. G., Neuer, T., Tafti, L., Raille, G., Pottelbergh, V., Pasieka, M., Skrodzki,
on computer vision, pp.
10

---

Shifting Time:
Time-series Forecasting with Khatri-Rao Neural Operators
A., Huguenin, N., Dumonal, M., Košcisz, J., Bader, D., Gusset, F,, Benheddi, M., Williamson, C.,, Kosinski, M., Petrik, M., and Grosch, G. Darts: User-friendly modern machine learning for time series. Journal of Machine Learning Research, 23 1-6 2022.
Kim, T., Kim, J., Tae, Y., Park, C., Choi, J.-H., and Choo, J. Reversible instance normalization for accurate time-series forecasting against distribution shift. In International
Conference on Learning Representations, 2021.
Kissas, G., Seidman, J. H.,, Guilhoto, L. F., Preciado, V. M., Pappas, G. J., and Perdikaris, P. Learning operators with coupled attention. Journal of Machine Learning Research, 23(215):1633, 2022.
Hochreiter, S. and Schmidhuber, J. Long short-term memory. Neural computation, 98):1735-1780) 1780, 1997. Publisher: MIT press.
Hom, M., Moor, M., Bock, C., Rieck, B., and Borgwardt, K. Set functions for time series. In International Conference on Machine Learning, pp. 4353-4363, PMLR, 2020.
Koopman, B. O. Hamiltonian systems and transformation in hilbert space. Proceedings of the National Academy of Sciences, 17(5):15-311, 1931a.
Huang, Q., Shen, L., Zhang, R., Ding, S., Wang, B., Zhou, Z., and Wang, Y. Crossgnn: Confronting noisy multivari- ate time series via cross interaction refinement. Advances in Neural Information Processing Systems, 36:46885- 46902, 2023.
Koopman, B. O. Hamiltonian systems and transformation in Hilbert space. Proceedings of the National Academy of Sciences, 17(5):315-31,, 1931b. Publisher: National Acad Sciences.
Kouw, W. M. and Loog, M. An introduction to do- main adaptation and transfer learning. arXiv preprint arXiv: 181..1 1806, 2018.
Jhin, S. Y., Jo. M., Kong, T., Jeon, J., and Park, N. Ace- node: Attentive co-evolving neural ordinary differential equations. In Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining, pp. 736-745, 2021.
Kovachki, N., Li, Z., Liu, B., Azizzadenesheli, K., Bhat- tacharya, K., Stuart, A., and Anandkumar, A. Neural oper- ator: Learning maps between function spaces with appli- cations to pdes. Journal of Machine Learning Research, 24(89):1-97,, 2023.
Jhin, S. Y., Lee, J., Jo, M.. Kook, S., Jeon, J., Hyeong, J.. Kim, J., and Park, N. Exit: Extrapolation and neural controlled differential equa- tions for seriss classification and forecasting. In Proceedings of the ACM Web Conference 2022, pp. $102-3112 12, 2022.
interpolation-base
Kuznetsov, V. and Mohri, M. Discrepancy- theory and algorithms for forecasting non-stationary time series. Annals of Mathematics and Artificial Intelligence, 88(4): 367-399, 2020.
Jhin, S. Y., Jo, M., Kook, S., and Park, N. Learnable path in neural controlled differential equations. In Proceedings of the AAAI Conference on Artificial Intelligence, 2023.
Lea, C., Vidal, R., Reiter, A., and Hager, G. D. Tempo- ral convolutional networks: A unified approach to ac- tion segmentation. In Computer Vision-ECcv 2016 Workshops: Amsterdam, The Netherlands, October 8-10 and 15-16, 2016, Proceedings, Part III 14, pp. 47-54. Springer, 2016.
Jhin, S. Y., Shin, H., Kim, S., Hong, S., Jo, M., Park, S., Park, N., Lee, S., Maeng, H., and Jeon, S. Attentive neu- ral controlled differential equations for time-series clas- sification and forecasting. Knowledge and Information Systems, (6(3):1885-1111 2024.
Li, A., Boyd, A., Smyth, P., and Mandt, S. Detecting and adapting to irregular distribution shifts in bayesian on- line learning. Advances in neural information processing systems, 34:681 18:6816-60 6-6828, 2021.
Johnson, A. E., Pollard, T. J., Shen, L., Lehman, L.-w. H., Feng, M., Ghassemi, M., Moody, B., Szolovits, P., An- thony Celi, L., and Mark, R. G. Mimic-iii, a freely ac- cessible critical care database. Scientific data, 3(1):1-9, 2016.
Li, X., Wong. T.-K. L., Chen, R. T. Q., and Duvenaud, D. Scalable Gradients for Stochastic Differential Equations. In Chiappa, S. and Calandra, R. (eds.). Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, volume 108 of Proceedings of Machine Learning Research, pp. 3870-3882. PMLR, August 2020a.
Kalnay, E., Kanamitsu, M., Kistler, R., Collins, W., Deaven, D., Gandin, L.. Iredell, M., Saha, S., White, G., Woollen, J., et al. The NCEP/NCAR 40-year reanalysis project. Bulletin of the American Meteorological Society, 77(3): 437 471, 1996.
Li, Z., Kovachki, N., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A., and Anandkumar, A. Fourier neural operator for parametric partial differen- tial equations. The International Conference on Learning Representations (ICLR 2021), 2020b.
Kidger, P., Morrill, J., Foster, J., and Lyons, T. Neural controlled differential equations for irregular time series. Advances in neural information processing systems, 33: 6696-6707, 2020.

---

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
Li, Z., Kovachki, N., Azizzadenesheli, K., Liu, B., Bhat- tacharya, K., Stuart, A., and Anandkumar, A. Neural operator: Graph kernel network for partial differential equations. arXiv preprint arXiv:2003.0324 03485. 2020c.
Makridakis, S., Spiliotis, E., and Assimakopoulos, V. The m4 competition: 100,000 time series and 61 forecasting methods. International Journal of Forecasting. 36(1):54- 74, 2020.
Li, Z., Kovachki, N., Azizzadenesheli, K., Liu, B., Stuart, A., Bhattacharya, K., and Anandkumar, A. Multipole Graph Neural Operator for Parametric Partial Differential Equations. In Larochelle, H.., Ranzato, M., Hadsell, R., Balcan, M. F., and Lin, H. (eds.), Advances in Neural Information Processing Systems, volume 33, pp. 6755- 6766. Curran Associates, Inc., 2020d.
Marchenko, V. The generalized shift, transformation oper- ators, and inverse problems. In Mathematical events of the twenticth century, pp. 145-162. Springer, 2006.
Menne, M., Williams Jr, C.. Vose, R., and Files, D. Long- term daily climate records from stations across the con- tiguous united states, 2015.
Mezić, I. Koopman operator.geometry. and learning of dynamical systems. NOTICES OF THE AMERICAN MATHEMATICAL SOCIETY, 6807):008711105 2021.
Li, Z. Huang. D. Z., Liu, B., and Anandkumar, A. Fourier neural operator with learned deformations for pdes general geometries. Journal of Machine Learning Research, 24(38)):1-26, 2023.
Milani, E., Andrade, M., and Diniz, C. Generalized normal ARMA model applied to the areas of economy, hydrol- ogy, and public policy. Communications in Statistics Simulation and Computation, 46(7)5511-55835 2017.
Liu, Y., Wu, H., Wang, J., and Long, M. Non-stationary transformers: Exploring the stationarity in time series forecasting. Advances in Neural Information Processing Systems, 35:9881 18933, 2022.
Montero, P., Athanasopoulos, G., Hyndman, R. J., and Ta- lagala, T. S. Fforma: Feature- based forecast model aver- aging. International Journal of Forecasting, 36(1):86-92, 2020.
Liu, Y., Li. C., Wang, J., and Long, M. Koopa: Learn- ing non-stationary time series dynamics with koopman predictors. In seieth Conference on Neural Information Processing Systems, 2023.
Nie, Y., Nguyen, N. H., Sinthong, P., and Kalagnanam, J. A time series is worth 64 words: Long-term forecasting with transformers. In The Eleventh International Conference on Learning Representations, 2023.
Liu, Y., Hu, T., Zhang, H., Wu, H., Wang, S., Ma, L., and Long, M. iTrrasffrmer:: Inverted Transformers Are Effective for Time Series Forecasting. In The Twelfth International Conference on Learning Representations, 2024.
Oh, Y., Lim, D.,, and Kim, Stable neural stochastic dif- ferential equations in analyzing irregular time series data, In The Twelfth International Conference Learning Representations, 2024.
ScuiSchiaffini, M., Singer, C. E., Kovachki, N., Schneider, T., Azizzadenesheli, K., and Anandkumar, A. Tipping point forecasting in non-stationary dynamics on function spaces. arXiv preprint XiXi:2338.00794 2023.
Oreshkin, B. N., Carpov, D., Chapados, N., and Bengio, Y. N-beats: Neural basis expansion analysis for interpretable time series forecasting. In International Conference on Learning Representations (ICLR), 2020.
Loshchilov, I. and Hutter, F. Decoupled weight decay regularization. In International Conference on Learning Representations, 2019.
Rahman, M. A., Ross, Z. E., and Azizzadenesheli, K. U- no: U-shaped neural operators. Transactions on Machine Learning Research, 2023.
Lu, L., Jin, P., Pang, G., Zhang, Z., and Karniadakis, G. E. Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nature Machine Intelligence, 3(3):218-229) March 2021.
Retherford, J. R. Hilbert space: compact operators and the trace theorem. Cambridge University Press, 1993.
Roberts, S., Oshome, M., Ebden, M., Reece, S., Gibson, N., and Aigrain, S. Gaussian processes for time-series mod- elling. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 371(1984):201 10550, 2013.
Lu, L., Meng. X., Cai, S., Mao, Z., Goswami, S., Zhang, Z., and Karniadakis, G. E. comprehensive and fair comparison of two neural operators (with practical tensions) based on FAIR data. Computer Methods in Applied Mechanics and Engineering, 393: 14778, 2022.
Rubanova, Y., Chen, T. Q., and Duvenaud, D. K. La- tent Ordinary Differential Equations for Irregularly- Sampled Time Series. In Advances in Neural Information Processing Systems, pp. 5321-5331, 2019.
Maas, A. L., Hannun, A. Y., Ng, A. Y., et al. Rectifier nonlinearities improve neural network acoustic models. In International Conference on Machine Learning, 2013.
12

---

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
Rusch, T. K., Mishra, S., Erichson, N. B., and Mahoney, M. W. Long expressive memory for sequence model- ing. In The Tenth International Conference on Learning Representations (ICLR 2022), 2022.
Wang, R., Dong, Y., Arik, S., and Yu, R. Koopman neural forecaster for time series with temporal distribution shifts. In International Conference on Learning Representations (ICLR), 2023.
Saatçi, Y. Scalable inference for structured Gaussian process models. PhD Thesis, University of Cambridge, December 2011.
Wilson, A. and Adams, R. Gaussian process kernels for pattern discovery and extrapolation. In International conference on machine learning, pp. 1067-1075. PMLR, 2013.
Salinas, D., Flunkert, V., Gasthaus, J., and Januschowski, T. DeepAR: Probabilistic forecasting with autoregressive recurrent networks. International journal of forecasting, 36(3):1181-1111 2020. Publisher: Elsevier.
Wilson, A. G., Gilboa, E., Nehorai, A., and Cunning- ham, J. P. Fast kernel learning for multidimen- sional pattern extrapolation. In Proceedings of the 27th International Conference on Neural Information Processing Systems- 2. pp. 3626-3634, 2014.
Schirmer, M., Eltayeb. M., Lessmann, S., and Rudolph, M. Modeling irregular time series with continuous recurrent units. In International conference on machine learning, pp. 19388-19405, PMLR, 2022.
Wu, H., Hu, T., Liu, Y., Zhou, H., Wang. J., and Long, M. Timesnet: Temporal 2d-variation modeling for gen- eral time series analysis. In The Eleventh International Conference on Learning Representations, 2023.
Shukla, S. N. and Marlin, B. Multi-time attention net- works for irregularly sampled time series. International Conference on Learning Representations, 2021.
Wu, Z.. Pan, S., Long, G., Jiang, J., and Zhang, C. Graph wavenet for deep The 28th International Joint Conference on Artificial Intelligence (UCCA). International Joint Conferences on Artificial Intelligence Organization, 2019.
graph modeling. In
Siami-Namini, S. and Namin, A. S. Forecasting economics and financial time series: ARIMA vs. LSTM. arXiv preprint arXiv:1803.03360 2018.
Smyl, S. A hybrid method of exponential smoothing and recurrent neural networks for time series forecasting. International Journal of Forecasting, 36(1):75 85, 2020.
Wu, Z., Pan, S., Long, G., Jiang, J., Chang. X., and Zhang, C. Connecting the dots: Multivariate time series forecasting with graph neural networks. In Proceedings of the 26th ACM SIGKDD international conference on knowledge discovery & data mining, pp. 753-763, 2020.
Stock, J. H. and Watson, M. W. Vector autoregressions. Journal of Economic perspectives, 15(4):101-115, 2001.
Tassa, Y., Doron, Y., Muldal, A., Erez, T,, Li, Y., Casas, D. d. L., Budden, D., Abdolmaleki, A., Merel, J., Lefrancq, A., al. Deepmind control suite. arXiv preprint 1rX.001801.00690, 2018
Yi, K., Zhang, Q., Fan, W., He, H., Hu, L.. Wang, P., An, N., Cao, L., and Niu, Z. Fouriergnn: Rethinking multivariate time series forecasting from a pure graph perspective. Advances in neural information processing systems, 36: 69638-69660) 2023.
Ticchi, A., Scherer, A., McIntyre, C., Brito, C. S. N., Snow, D., Develra, dstern, Colless, J., Gar- vey, K., Maggie, Ortiz, M. P. Lynch, R., and Dane, S. G-research crypto forecasting, 2021. URL https: /aap::///agglleee com/ wwmpetitiooioooooooooooo g research -rrpptttoc-o-ttooooooo forecasting.
Zeng, A., Chen, M., Zhang, L., and Xu, Q. Are transformers effective for time series forecasting? In Proceedings of the AAAI conference on artificial intelligence, 2023.
Zhang, J., Zheng, S., Cao, W., Bian, J., and Li, J. Warp- former: A multi-scale modeling approach for irregu- lar clinical time series. In Proceedings of the 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, pp. 3273-3285, 2023.
Van Loan, F. The ubiquitous Kronecker product. Journal of computational and applied mathematics, 123(1-2):85 100, 2000.
Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, E., and Polosukhin, I. At- tention is all you need. Advances in neural information processing systems, 30, 2017.
Zhang, W., Yin, C., Liu, H., Zhou, X., and Xiong, H. Irreg- ular multivariate time series forecasting: A transformable patching graph neural networks approach. In Forty-first International Conference on Machine Learning, 2024.
Wang, R., Maddix, D., Faloutsos, C., Wang. Y., and Yu, R. Bridging moysics-based and data-driven modeling for learning dynamical systems. In Learning for dynamics and control, pp. 385-398. PMLR, 2021.
Zhang, X., Zeman, M., Tsiligkaridis, T., and Zitnik, M. Graph-guided network for irregularly sampled multivari- ate time series. In International Conference on Learning Representations, 2022.
13

---

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
Zhang, Y. and Yan, Crossformer:
Transformer utilizing
cross-dimension dependeney for multivariate time series forecasting. In The eleventh international conference learning representations, 2023.
Zhen, Y., Chapron, B., and Memin, E. Bridging koop- operator and time-series auto-correlation based Hilbert-Schmide operator. Stochastic Transport in Upper Ocean Dynamics, pp. 1 1-20, 2022a.
Zhen, Y., Chapron, B., Mémin, and Peng, L. Eigen- values of autocovariance matrix: A practical method to identify the Koopman eigenfrequencies. Phys. Rev. E,
105, 2022b.
Zhou, T., Ma. Z., Wen, Q. Wang. X., Sun, L., and Jin, R. Fedformer: Frequency enhanced decomposed trans- former for long-term series forecasting. In International Conference on Machine Learning, 2022.
14

---

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
A. Continuity of the time-shift operator for ODEs
Lemma 1. Let ((T) = (()(T)), where F R" + TR" is Lipschitz continuous and T E Upttz] Then the time-shift operator
A (1pp 1); R*) - ((t, [y]: R.) is continuous in the 1,2 sense, L.e, Ep < < 1g, there exists a constant C >>> 0 such that
||4 21 - A 22 (21112( S Ci|21 - 22 (21.12(()))))))))
where and 22 are trajectories of = Fz) corresponding to different initial conditions.
Proof. Let e(T) - (1 (2( with r E 1pp. 4/). where 21, 22 E ((p ky); to different initial conditions. Then, we have
R2) are two ODE trajectories corresponding
|| 21
4 A,, 22 [(((),)))=)) - || - =2 =2ll<>((,,)) - ||e|l
(9)
Noting that ||él|2 [21(2)) P S Lplle|2,, where LF > 0 is the Lipschitz constant of F, we have
d dr |e(()] = 20(T)² de dr < 2l(((T))l1
de
= 2L T [tp, tyl.
(10)
dr
Applying Grönwall's lemma (Ames & Pachpatte, 1997) to the preceding inequality, we have - J." [((()) S (1e(())l '." 220((T-))dr.
| ((1))))))
(11)
Since the trajectories are continuous over [tp, ty]. it follows from the extreme value theorem that [| is bounded over this interval which in turn implies that & E _ (l4p. Lr]; R²) C (2 ();;)")) RR). In addition, since = is square integrable, we
have = E H ([tp. 1.t) Noting that H' ([tp: f); + C(Itp tyl: RP due to the Sobolev embedding theorem, there exists an embedding constant C1 > 0 such that
|e((p))laa I ||ell S C1 (Ieell>(((t.)))))))))))) + [161123(66.))22m))
C1(1+ [F)eeleell))
(12)
Using (9), (11), (12), we have
| 21 A 12 22(((3(()))))))) S C1Is2
22 22lle>((e.)))>>>
where C'
Cx(++L)) V2LF
exp(22 - 1p)) ct
The continuity of the spatio-temporal time-shift operator can be established for time-dependent PDEs under appropriate regularity assumptions. We leave this for future work.
The notion of the continuous time-shift operator can be extended to forced dynamical systems. Consider forced ODE of the form: =(7) = ((T)), where =(r) E Rn is the state, {(T) R.m an exogenous foreing function (or control input), and F : R. x T - I2" is assumed to be sufficiently regular. Similar to the unforced case, we can define a forced time-shift operator, denoted by of which maps the history of the state z over a look-back window [tp. t) and the forcing function & defined over the interval [tp. tyl (which includes both book-back and prediction periods) to the future values of the state z over the interval (t, ty]. The forced time-shift operator ofy is therefore a map from T,2 ((tp 4] IR²) x [14p Ly]; TT") to T.? ty]: R²), which we can define as follows:
Time-shift operator for forced dynamical systems:
x(7) = (g/."" (2/(s..) Eleent)) (7), VJ E (t.tyl.
(13)
where = 2||6,.0 and Ellt,,,,) denote the restriction of z and & to the intervals [4p 1) and [4p, (y], respectively. When parametrizing the forced time-shift operator, is critical to ensure causality, i.c., the predictions made at time instant 7 are not influenced by the future values of the forcing function. This can be achieved by defining the integration domain of the kemel integral transform with respect to & as [tp. 7j, where E (t,tyl denotes the time-stamp at which a prediction is being made. For example, a linear integral transform representation of the forced time-shift operator would take the form: <(T) Jr, R=(T, s)2(s)ds + [r ()(TT,){(ss))) The notion of forced time-shift operators be similarly extended to spatio-temporal systems, enabling the application of our operator-learingg framework to systems with exogenous inputs.
15

---

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
B. Proof for Proposition 2.1
start by briefly clarifying what meant by the quadrature nodes lying on a product grid. An example a two-dimensional product grid is provided in Figure below.
Figure 6. Example of quadrature rule which lies product grid. Circled in blue the quadrature nodes along the first dimension and circled in red are the quadrature nodes along the second dimension.
Proposition states that
1. the quadrature nodes lie product grid, X = 8 x ax1)) ((d1)) where i,)(() E IR" indicate the quadrature nodes along the time dimension and the 7thh dimension of x respectively such that N nd (for the general case N II 1 n2) and
2. the kernel function has component-wise product structure of the form given in Equation (7) reproduced below for clarity:
<(t,,)), ((,t)) x(,)t)))))) @ ( [=k--)
Then the kernel function evaluated at the quadrature nodes inherits the Khatri-Rao product structure provided in Equation (8) and reproduced below.
<(X, -x * ( ( 100 (x((i))) ,z((-
Here (.()),,)) R is block matrix where block jk is the ju: output from the component kernel (()) evaluated on the outer product of the quadrature nodes along the 7th dimension.
Proof. We start by observing that s(X,X) can be dlock-partitioneed into x p blocks of size N x N.
K1,1
K1.2
K1
K2,1
12,
K2
R
xX,x)=
=
(14)
Fig.1
Kq.2
Kgg.P_
Each of these N N blocks inherits the product structure of Equation (8).
Kj,k 1f 1 (X
1]),
(15)
where ( (X):, 11, 1]) E NNN, is the jh output of the 1th component kernel function evaluated on the tth dimension the quadrature nodes. The th block in the kernel evaluated at the quadrature nodes be written as the Kronecker product as follows (Van Loan, 2000),
KjAA = (1) (T,T) @ (0 (2) ,c(-1)))
(
(16)
16

---

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
where (1 (xl x ()) R each x,,(() E RR indicate the one-dimensional quadrature nodes along the time and 7th j,k dimension respectively. The Khatri-Rao product structure follows from substituting (16) into (14).
The above proof can be generalized to cases where the quadrature nodes distributed a product grid with variable number of nodes along each dimension, i.., N 11 12
C. Generalization of Proposition 2.1
In this section, we generalize Proposition to the case where the input and output are defined over different spatio-temporal domains and different quadrature nodes are used for the input and output functions. We start by rewriting the kernel integral transform layer as a map with the input function vt : See x Ie + I2 and the output function vt41 le41 x Ie+1 + 129 as follows
vet ((+111)1,, K(oe)
S I (0.r) x((te-1,ze11h, Ut (f)d)t + W (or + b,,
where k: IR Oig x R. X le++1 - R a matrix-valued kernel, W E RR** a weight matrix, 2 : le+1 - S2z is a map between the output and input domains, and IR9 is bias vector.
Let X₃ and X4+1 denote the sets of quadrature nodes for the input and output domains, respectively. The quadrature nodes over the domain of the input function are assumed to lie on product grid, i.e., X, = fe x x x x...x, x/ (d-1) where E R" and 2t E RT denote the quadrature nodes along the input time dimension and the 7th dimension of mr. respectively, such that N₂- = nę (for the general case when the number of quadrature nodes along the ht dimension is ne, have N2- (i,)),e)) Similarly, the quadrature nodes over the output domain are assumed to lie product (d-1)
grid, X(+1 = fe+1 x x'41 X x... where fe+1 E [R- and 2++1 E IR denote the quadrature nodes along the output time dimension and the 7h dimension of Tz+1, respectively, such that Ne41 n (for the general case with a different number of quadrature nodes along each dimension Ne+1 II "e+1,,) As before, we will consider kernel with component-wiss product structure of the form given in Equation (7).
Similar to the previous proof, we start by observing that x(Xe-1. X²) can be block-partitioned into x blocks of size N4+1 x No, i.e.
K1,1
11,2
11,.
K2,1
K2.p K2,2
x(Xe+1,
(17)
:
Kg.1
Kgg2
Kig.p_
Each Ne+1 N₂ block inherits the product structure, i.e., Kj,AA
". j.b
(t) (Xetil,, 1], Xe(e, -]), where (Xe+1
. - 11, X - 1]) E IR Ne x is the jkth output of the 2th component kernel function evaluated on the 2th dimension of the quadrature nodes. The jhth block in the kernel evaluated at the quadrature nodes can be written as
Ki,kk = (tet1, & (ol 2
(2) (i-1). ( . 1)
(18)
(1-1) where (f x ()) R 1 xxx, Substituting (18) into (17), we have
m(Ne+1, X,) = (T+1,, *
10 * i=2
(i-1) (x(tt)))
(6 1
(19)
where (r,-) E 19 is block partitiort matrix where block jk is the jh output from the component kernel 1e() evaluated on the outer product of the quadrature nodes along the 7th dimension.
It follows from this result that we retain the original computational complexity of the KRNO operator even in situations where the inputs and outputs are defined over different domains. In addition, this result provides the flexibility of designing memory-efficient multi-resolltion neural operators, where the hidden layers operate on arroluierrsoluution representations of the input function. This generalized KRNO integral transform layer can be viewed as a continuous analog of upsampling and downsampling layers used in convolutional neural networks.
17

---

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
D. A Algorithm for Khatri-Rao structured matrix-vecton products
In this section, we present an algorithm to efficiently compute the matrix-vector product associated with the Khatri-Rao product structured matrix defined in (8), without the need explicitly construct the full matrix of size qN x p.V.
Let A E RoNxpN be block structured matrix of the form,
A1,1
A1
A1.p
A2,1
Az,2
A2.p
A
(20)
:
Ag,1
Ag.2
Aq.z_
where each A₁ 8/=14 and A R Assuming q.P << N, the computational complexity associated with the matrix-vector product u = Au can be reduced from O(N2) to O(N++//a) In addition, the memory requirements are also reduced from (N2)) to O(N2/d +N.. An efficient PyTorch implementation outlining the steps is provided below.
def khatri rac._araa___________
A: _____rrr____________ improd
Tensor, "q nl I V: Float Tensor, "pN batch")
Float Tensor, "qN batch"l:
len (A) size the product grid of kernel components)
9, D, -
shape
pN, bs - V shape
X v.reshape (p, -1, ttransassp -1)
for in ra)ge(d):
Gd s[i]_shape[-1]
prep x. shapel
X x.reshape ( (*bs_prod, Ed, -1)
Alil (-3) 0
x.trannspos -1) reshapehe P, bs, -1)
return x.sum (1) transpose (-2, -1). reshape (-1, bs)
We note that the above algorithm is applicable to Khatri-Rao product structured matrix, as defined in (19), where the inputs and outputs are defined over different spatio-temporal domains (with each domain using different set quadrature nodes).
E. Details on KRNO parametrization
As mentioned in the paper, we parametrize each component-wise kernel, k(0 R R - R2 by neural network. All neural nets use skip connections and layer normalization (Ba et al., 2016). In addition, before passing an input into the component function, we apply an input transformation Φ IR x
+ R.m
((t,t')
2 cos ([t
t' w B)
(21)
where w E R2xm and 3 E I Such input feature transforms were found to be beneficial in prior works (Kissas et al., 2022)
F. Practical aspects of learning the time-shift operator
Consider n-dimensional multivariate discrete time-series dataset {2t (2))T 1- This dataset is first converted into pairs of input and output sequences over two non-overlapping time intervals [4, 1] and (1, where 0 < lp < < <, S T, for various time instances t. The input sequence, denoted by U = {z { includes 2; given at P' time steps within the look-back window [tp, t]. The output sequence, denoted by Uj = )/ contains values of 2, given at H time steps within the prediction window (t, ty]. These pairs of sequences, for different values of t, are used to approximate the continuous time-shift operator 1/.5 using KRNO. For estimating the KRNO parameters, we minimize the loss function,
1 EM, |0 - A,, U/: [E>(().,t)):>>> where M is the number of input-output sequence pairs. As discussed previously, M7
bp and fy are hyperparameters of the time-shift operator which are chosen using cross-validation.
Similar to temporal datasets, we consider spatio-tempores data comprising discrete snapshots of spatial fields (u((x) where x E S2 with G, representing a spatial grid over Q. This dataset is converted into pairs of input and output sequences
18

---

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
of spatial fields, (3(z)) and (x), over time intervals U4, 1) and (l, 4y], where ly < < T. The input sequence Up(x)) = {ur(z)) contains spatial fields corresponding to P time steps within the look-back window [t,, t]. The output sequence (z) (u( contains spatial fields over H time steps within the prediction window (t. ty). These sequence pairs are used to learn the spatio-temporal time-shift operator (3) using KRNO by minimizing the loss function,
M E== 1 fe.1 U
1 M
where M is the number of onput- sequence pairs.
Table 8. Summary of all datasets (along with corresponding test cases and error metrics) considered in our numerical experiments. This leads to a total of 39 different test cases or error metrics. The top five models are listed for each case. It can be seen that KRNO achieves top performance on 17/39 cases and top-3 performance 30/39 cases.
Dataset
Test case
Best
Second
Third
Fourth
Fifth
(or error metric)
Darey flow
6
KRNO
FNO
POD--DeepONet DeepONet
DeepONet
Hyper-elastic
Eg
KRNO
FNO
6p
KRNO
HNO-3I)
LOCA
Shallow water
KRNO
FNO-3D
LOCA
Cup
KRNO
FNO-3D
LOCA
Short traj.
KRNO
Latent-ODE
T-PAICHGNN
2D spiral
Long traj.
KRNO
Latent-ODE
1-PATCHGNN
MSE
KRNO T-PATCHGNN
StemGNN
Warpformer
GRU-D
MIMIC
MAE
T-PATCHGNN
KRNO
GRU-D
Warpformer
StemGNN
MSE
KRNO
PAPCHHGNN
Crossformer
Warpformer
Graph Wavemel
USHCN
MAE
KRNO
T-PATCHGNN
Graph Wavenet
Warpformer
Neural Flow
MSE
T-PATCHGNN
Warpformer
KRNO
Graph Wavenet
GRU-D
Human Activity
MAE
T-PATCHGNN
Warpformer
Graph Wavenet
FOURIER/GNN
KRNO
Regular
KRNO
Neural LSDE
Neural LNSDE
Neural GSDE
LEAP
30%
KRNO
Neural GSDE
Neural LNSDE
Neural GSDE
LEAP
MuJoCo
50%
KRNO
Neural GSDE
Neural LNSDE
Neural GSDE
LEAP
70%
KRNO
Neural LSDE
Neural LNSDE
Neural GSDE
LEAP
AirPassengers
LLaMA-2
ARIMA
KRNO
GPT-3
SM-GP
AusBoer
N-BEATS
LLaMA-2
GPT-3
KRNO
ARIMA
GasRateCO2
SM-GP
KRNO
ARIMA
LLaMA-2
N-BEATS
MoathlyMilk
GPT-3
LLaMA-2
KRNO
SM-GP
N-HITS
Darts
sunspots
KRNO
ARIMA
GPT-3
LLaMA-2
N-HITS
Wine
KRNO
GPT-3
ARIMA
TCN
N-HITS
Wooly
N-HITS
ARIMA
SM-GP
KRNO
GPT-3
HeartRate
TCN
GPT-3
SM-GP
KRNO
N-HiTS
Monthly
KNF
Nbeats-I+G
Smyl
Montero al
KRNO
Weekly
KRNO
KNF
Montero et al
Smyl
Daily
KNF
KRNO
Montero et al
Smyl
M4
Hourly
Smyl
KNF
Montero al
KRNO
Yearly
Nbeats-I4G
Smyl
Montero et al
KNF
KRNO
Quarterly
Nbeats-I+G
Smyl
Montero et al
KNF
KRNO
(1 ~
MLP+RevIN+TB
KNF
KRNO
LEM
VARIMA
(6 10)
KNF
KRNO
MLP+RevIN-TB
LEM
FedFormer
Crypto
(11 ~ 15)
KNF
KRNO
LEM
MLP+RevIN+TB
RF+TB
(1 ~ 15)
KNF
KRNO
MLP+RevIN-TB
LEM
FedFormer
(1 10)
VARIMA
KNF
KRNO
LEM
RIP+
(11 ~ 20)
KNF
VARIMA
FedFormer
KRNO
MLP+RevIN-TB
Player Traj
(21 ~ 30)
KNF
KRNO
FedFormer
VARIMA
LEM
15)
KNF
KRNO
VARIMA
FedFormer
LFM
G. Comparison of training and inference costs of KRNO and FNO
In this section, we compare the computational cost and memory requirements of KRNO and FNO-3D for different spatial resolutions using the shallow water dataset. We utilized the default settings for KRNO described earlier, except for the quadrature grid. For all the spatial resolutions, we employed a quadrature grid of size 32 x x22 x 5 in the hidden KRNO integral transform layers. is important to note that for high- high-resolution spatial data, the memory requirements and computational cost are primarily driven by computations in the first and last integral transform layers, which contain the highest number of
19

---

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
Table Comparison of training memory requirements and computational cost of KRNO and FNO-3D models on shallow water problem for different spatial resolutions (with batch size of 8), Two FNO-3D configurations are compared: one with fixed number of Fourier modes m 12. and another with m + that scales with the spatial resolution S x S.
Spatial resolution
GPU memory (MB)
Time (seconds)
# Quadrature nodes
KRNO
FNO-3D
FNO-3D
KRNO
FNO-3D
FNO-3D
(m 12)
(m
+1)
(m 12)
(m + 1)
32 32
5,120 (N)
1,390
854
1,074
0.0279
0.0216
0.0234
64 64
20,480 (4N)
2,776
1,350
2,772
0.0394
0.0252
0.0347
96 96
46,080 (9NN)
4,884
2,124
5,264
0.0626
0.0405
0.0710
128 x 128
81,920 (16N)
7,608
3,318
8,994
0.0999
0.0704
0.1297
160 x 160
128,000 (25N)
10,040
4,872
13,764
0.1584
0.1114
0.2088
Note: = 32 32 x
5120.
Table 10, Comparison of inference memory requirements and computational cost for KRNO and FNO-3D models on shallow water problem for different spatial resolutions (with batch size of 8). Two FNO-3D configurations compared: with fixed number of Fourier modes m 12, and another with m = 2u + that scales with the spatial resolution x S.
Spatial resolution
GPU memory (MB)
Time (seconds)
# Quadrature nodes
KRNO
FNO-3D
FNO-3D
KRNO
FNO-3D
FNO-3D
S x S
(m = 12)
(m = g 11)
(m = 12)
(m _
+ 1)
32 x 32
5,120 (N)
708
942
1,074
0.0107
0.0065
0.0062
64 x 64
20,480 (4N)
1,314
1,294
2,772
0.0134
0.0069
0.0071
96 x 96
46,080 (9N)
2,366
1,622
5,264
0.0185
0.0108
0.0149
128 X 128
81,920 (16N)
3,796
2,428
8,994
0,0312
0.0210
0.0288
160 X 160
128,000 (25N)
5,644
3,292
11,288
0,0502
0.0315
0.0456
quadrature nodes. For the FNO-3D model, the number of Fourier modes (m) in each spatial dimension has significant impact the computational complexity. For instance, when training on data with spatial resolution 160 x 160, increasing from 12 to 81 increases the memory requirements from 4, 872 MB to 13, 764 MB and the training time per iteration from 0.1114 seconds to 0.2088 seconds. In addition, the parameter count of FNO-3D dramatically increases from 5.5 million to 252 million when m increases from 12 to 81. In contrast, KRNO maintains fixed parameter count of 145, 319, independent of the size of the quadrature grid and the resolution the dataset.
The memory requirements and computational cost during training and inference are shown in Table 9, Table 10, and Figure 7. As discussed earlier, the memory requirements and time complexity of FNO-3D grow dramatically with an increase in the number of Fourier modes. We observe that the memory usage of KRNO is significantly higher during training than at inference. We believe that this can be reduced by further optimizing the implementation of Khatri-Rao matrix-vector products.
We would like to highlight that our current implementation uses a mid-point quadrature scheme for the temporal dimension and a trapezoidal quadrature scheme for spatial dimensions to evaluate the integral transform layers. Additionally, in our current implementation, the quadrature nodes are defined over the spatial mesh of the input function. While this approach reasonable for the problems we are considering, high-resolution spatio-temporal datasets. For such datasets, the memory requirements be significantly reduced by using quadrature nodes that are defined over a lower-resolution spatial grid which is independent of the input's spatial resolution.
not the most suitable (and tends be overly conservative) for
Alternatively, we could design a quadrature scheme that uses a low-resolution subsampling of the input as nodes and generates quadrature weights on-the-ff fly to meet a specified target precision. This would not only enhance accuracy and efficiency for high-resolution observations. We plan to explore this in future work.
datasets but also improve the performance KRNO for irregularly spaced
20

---

Shifting Time: Time-series Forecasting with Khatri- Neural Operators
0.200
KRNO (train)
FNO-3D (train)
KRNO (test)
Thme [secondds
FNO-3D (test)
0.125
0.100
0.075
0.050
0.025
No of points in spatio-temporal grid
Figure 7. Training and inference time per iteration for KRNO and FNO-3D models on the shallow water problem. plotted as function of the number of quadrature nodes. For the FNO-3D model, the number of Fourier modes is set m where S is the spatial resolution. See the discussion in Appendix G for more details.
H. Additional experimental details
This section provides additional details on the experimental setup used to generate the results presented in the main text. The default KRNO architecture used in most of the experiments has 3 kernel integral layers with 20 channels each, lifting and projection layers are parametrized by MLPs with one hidden layer containing 128 hidden units, and the kernels in the integral transform layers are parametrized by MLPs with 3 hidden layers with 32 hidden units each. We used the SiLU (Elfwing et al., 2018) activation function as the nonlinearity. In each test case, we used the same input and output quadrature grids in each kernel integral layer except for the final layer. For problems involving high spatial or temporal resolutions, adopting lower-resolution quadrature grids within the internal kernel integral layers recommended as an effective strategy reduce computational costs. The hyperparameter settings used for Darts, M4, Crypto, and Player Trajectory datasets are summarized in Table 12 and Table 14. The AdamW (Loshchilov & Hutter, 2019) optimizer is used for training all the models. All the computations were carried out on a single Nvidia RTX 4090 with 24GB memory.
In all experiments, we treat the fp and fy as fixed hyperparameters. We would like to highlight here that further work is needed to explore the possibility of training a single model on a dataset containing input/output trajectories for different settings of tp and ty. This would enable the possibility of learning a model that can predict the dynamics at different length scales.
H.O.1. ADDITIONAL DETAILS ON MIMIC, USHCN AND HUMAN ACTIVITY BENCHMARKS
MIMIC The MIMIC dataset (Johnson et al., 2016), publicly available clinical database, contains electronic health records from critical care patients. After applying preprocessing steps outlined by Biloš et al. (2021); Zhang al. (2024), we extracted time series data for 23,457 patients. Each time series encompasses 96 variables recorded during the initial 48 hours post-admission. The goal here to forecast target values for the subsequent 24-hour period using the data from the first 24 hours as observations.
The United States Historical Climatology Network (USHCN) dataset (Menne et al., 2015) is utilized to address the common challenge of missing data in climate research, which often arises from issues like sensor malfunctions or data acquisition errors. This dataset originally encompasses daily measurements of five climate variables (precipitation, snowfall, snow depth, minimum and maximum temperature) collected from numerous meteorological stations across the United States. Following preprocessing methodologies adopted from previous studies (De Brouwer et al., 2019b; Schirmer et al., 2022), specific subset is used, focusing on 1, 114 stations and an observational period spanning four years, from 1996 to 2000. We followed the preprocessing steps described by Zhang et al. (2024), which resulted in 26, 736 multivariate time series. The primary task here is to forecast the climate conditions for the subsequent 24 months using the climate data from the preceding 24 months for each time series.
USHCN
21

---

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
Human Activity The Human Activity dataset contains 3D positional data from four sensors worn on the left ankle, right ankle, belt, and chest. This data, initially consisting of 12 irregularly measured variables from five individuals performing activities like walking and sitting, was preprocessed according to Zhang et al. (2024). This resulted in 5, 400 irregularly- sampled multivariate time series, each spanning 4, 000 milliseconds. For forecasting, we use the initial 3, 000 milliseconds of each irregularly sampled multivariate time-series as observed data to predict sensor positions for the subsequent 1, 000 milliseconds.
KRNO architecture details: The MIMIC, USHCN, and Human Activity benchmarks involve irregularly sampled, multivariate time-series with missing observations. To accommodate the missing observations, we need to augment the KRNO architecture with a learnable transformation that maps the irregularly sampled observations with missing states to a feature vector defined over the union of all observation time-stamps across all variables/states. We experimented with different preprocessing architectures and found that the best performance was achieved by augmenting the KRNO architecture with a simple one-dimensional CNN preprocessing module that is shared across the variables (states). This is similar in spirit to approaches used in T-PATCHGNN (Zhang et al., 2024) and PatchTST (Nie et al., 2023), where all the model weights are shared across all the variables (states). It is worth noting here that in the present work only the CNN preprocessing module weights are shared across the states, while the KRNO architecture does not use weight sharing.
We take the union of all time-stamps corresponding to each variable (state) in a single trajectory, and create a binary mask to indicate the missing observations. The CNN module takes this augmented sequence for each variable (state) along with the corresponding mask (treated two-channel input) and transforms it to single-channel feature vector defined over the union of all time-stamps. The transformed sequence of features extracted from the CNN preprocessing module along with the union of time-stamps are then used as inputs to the KRNO architecture. For all three datasets, we used a CNN preprocessing module composed of three layers, each employing a kernel size of 9 with stride length of 1 and padding of 4. The number of output channels for the three layers is set to 30, 20, and 1, respectively. We used the LeakyReLU (Maas et al., 2013) activation function as the nonlinearity. The KRNO architecture used is the default configuration described in Section H. The AdamW optimizer used for training and the learning rate is set to 10-3, The training, validation and test splits are created using 60%, 20%, and 20% of the data, respectively.
Our current KRNO implementation for these datasets uses batch size of 1 since the length of the transformed sequences for each trajectory is different. Further work is required to explore how this limitation can be overcome. One promising approach involves converting variable-length sequences into fixed-length feature vectors, analogous to the methodology in T-PATCHGNN (Zhang et al., 2024), thereby enabling their use as batched inputs for the KRNO model.
Baseline models On irregular time-series benchmarks, we compare the performance of KRNO with numerous SOTA mod- els, including DLinear (Zeng et al., 2023), TimesNet (Wu et al., 2023), PatchTST (Nie et al., 2023), Crossformer (Zhang & Yan, 2023), Graph Wavenet (Wu et al., 2019), MTCNN (Wu et al., 2020), StemGNN (Cao et al., 2020), CrossGNN (Huang et al., 2023), FOURIER/GNN (Yi et al., 2023), GRU-At, GRU-D (Che al., 2018), SeFT (Horn et al., 2020), Rain- Drop (Zhang et al., 2022), Warpformer (Zhang et al., 2023), mTAND (Shukla & Marlin, 2021), CRU (Schirmer et al., 2022), Neural Flow (Biloš et al., 2021), T-PATCHGNN (Zhang et al., 2024), GRU-ODE (De Brouwer et al., 2019a). ODE-RNN, Latent-ODE (Rubanova et al., 2019), Augmented-ODE (Dupont et al., 2019), ACE-NODE (Jhin et al., 2021). NCDE (Kidger et al., 2020), ANCDE (Jhin et al., 2024), EXIT (Jhin et al., 2022), LEAP (Jhin et al., 2023), Neural SDE (Li et al., 2020a), Neural LSDE (Oh et al., 2024), Neural LNSDE (Oh et al., 2024), Neural GSDE (Oh et al., 2024).
H.1. 2D spiral example
This section presents additional numerical results to evaluate the performance of KRNO on irregularly spaced time-series data obtained for the two-dimensional spiral test case from (Chen et al., 2018) and compare against T-PATCHGNN (Zhang et al., 2024) and Latent-ODE (Rubanova et al., 2019). First, we consider 10 irregularly spaced short training trajectories. which represent high levels of irregularity, alongside one equispaced training trajectory. The equispaced trajectory is generated by sampling the states at 100 evenly spaced time stamps over the short time interval [0, 15.61] seconds. To generate training datasets with a high level of irregularity, 100 new time-stamps are obtained by adding random noise ~ (0, 0.156) to the time stamps of the equispaced trajectory. Subsequently, the training trajectory is obtained by sampling the states at the randomly perturbed time-stamps. This randomization procedure was repeated to create 10 different training datasets where the distribution of the time-stamps is highly irregular, making this test case more challenging. Similarly, we generated 10 irregularly spaced trajectories alongside one equispaced training trajectory over the longer time
22

---

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
interval [0, 65.88] seconds.
We consider T-PATCHGNN (Zhang et al., 2024) and Lat (Rubanova et al., 2019) as baseline methods. We trained all models on both short and long trajectories. For models trained on short trajectories, performance is evaluated on test dataset containing the future states at 18 uniformly spaced time stamps over the interval (15.61, 18.45] seconds. While for models trained on long trajectories, performance is evaluated on a test dataset containing the future states at 18 uniformly spaced time stamps over the interval (65.88, 68.72] seconds. All models are trained by fixing the length of the time window corresponding to the input and the output sequences to be 1.419 seconds. We did not tune the KRNO hyperparameters for this experiment. The KRNO architecture consisted of 128 channels in both the lifting and projection layers, with three kernel integral layers, each containing four channels.
The predictive performance of all models is compared in Table 11. When the training trajectories are equispaced, the performance of both Latent-ODE and T-PATCHGNN models is comparable, with KRNO outperforming them by factor training data, we provide the statistics of test MSE over 10 different randomly
of 10. For the case of irregularlysppacee sampled trajectories generated using the procedure described earlier. Table 11 shows that KRNO significantly outperforms both Latent-ODE and T- T-PATCHGNN models on this extrapolation task.
Sample predictions from models trained on short and long trajectories are shown in Figures 8, 9, respectively. KRNO is the only model that is able to consistently predict the future states on both short and long trajectories. Both Latent-ODE and T-PATCHGNN models revert to predicting the training data. In all cases, on average, Latent-ODE performs better than T-PATCIGNN. As shown in Figure 9, there is significant variance in the predictive performance of T-PATCHIGNN across different irregularly spaced trajectories.
We note that reduction in predictive accuracy when the distribution of the time-stamps is highly irregular is influenced by the time intervals where the sampling frequency is low. For instance, lower sampling frequency in the time interval close the testing time window has the most significant impact on predictive accuracy. In such situations, the hyperparameterrs tp and fy would need to be carefully selected to improve the predictive performance.
Another important factor that impacts the predictive accuracy for trregularly-spaced training observations is the quadrature scheme used to approximate the kernel integral transform layers. Our current implementation uses a trapezoidal quadrature rule which is not ideal when the training data is sampled at a highly irregular frequency. It is expected that by adopting a quadrature scheme that obtains weights on-the-fly fly for irregularly spaced data (while meeting a target precision), the accuracy can be improved further.
Table 1. Comparison of test MSE of models trained one equispaced trajectory and 10 trajectories (mean/std) with irregularity in the distribution of the time-stamps.
Train Trajectory
Test MSE
Length
Type
KRNO
Latent-ODE
T-PATCHGNN
Equispaced
0.081
0.931
1.944
Short
Irregular
0.147 + 0.06
0.929 + 0.12
1.740 + 0.40
Equispaced
0.768
4.939
5.221
Long
Irregular
0.974 0.68
6.828 + 3.91
64.023 + 58.9
H.2. Darts benchmarks
For all the datasets in Darts, we used 60%-20%--00% as a train-validatiion-ees split. We performed a grid search on the Darts datasets using the hyperparameters listed in Table 12 to find the optimal hyperparameters. Model selection was done based on the normalized mean absolute error (NMAE) (22) on the validation set. Since the available training data in the Darts dataset is not sufficient to train a deep network, we conducted weight decay tuning to determine the optimal weight decay value using the selected hyperparameters. This optimal weight decay value was then used to train the final model using both the training and validation data. This final model is used to obtain predictions by forecasting recursively until the end of the testing window, shown in Figure 10. The evaluation metric, normalized MAE (NMAE), is computed as follows
MAE(Y,S 8) 1 Li=1 M
MAE $)
NMAE 8)
(22)
mean
where y and y are the truth and predicted time series.
23

---

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
Table 12. Hyperparameter tuning ranges used for Darts dataset.
Learning rate
Integral layer channels [5, 10, 32]
Hidden
Look-back
Prediction
ReVIN
units in
window length
window length
kernel
[le-3,5e] 5e-3]
[32. 64]
10 to 100
5 to 100
[True, False]
H.3. M4 benchmarks
We utilized the training and test datasets from the M4 competition (Makridakis et al., 2020). For all M4 datasets, the last 10% of the data for each time series in the training data is used as validation data. The testing process involves forecasting for specified time period (testing window length) for each seasonality. The testing window lengths for each seasonality are shown in the parentheses next to the seasonality in Table 13.
Table Comparison sMAPE from KRNO with other baseline methods on datasets. Results with (-) were taken from Wang (2023).
Method
Monthly(18)
Weckly(13)
Daily(14)
Hourly(48)
Yearly(6)
Quarterly(8)
Montero et al. (2020)
12.639
7.625
3.0971
11.506
13.528
9.733
Smyl (2020)
12.126
7.817
3.170t
9.328
13.176
9.679
Nbcats-I+G
12.024
12.924
9.212
KNF (Wang et al., 2023)
11.930
7.254
2.990
11.294
13.800
10.008
KRNO(ours)
13.432
6.934
3.086
11.686
14.302
10.503
H.4. Crypto and Player Trajectory benchmarks
For these two datasets, used the same train-test splits Wang et al. (2023). Similar to the M4 datasets, 10% of the data corresponding to each time series in the training data is used for validation. Following (Wang et al., 2023), the prediction window lengths for the Crypto and Player Trajectory datasets were set to 15 and 30, respectively. Table 15 presents comparison of prediction errors between KRNO and baseline methods.
For the Crypto dataset, Table 15 details the test prediction errors for time steps 1-5, 6-10, and 11-15, along with the total error across all 15 steps. Similarly, for the Player Trajectory dataset, Table 15 reports errors for time steps 1-10, 11 11-20, and 21-30, and the total error over all 30 steps.
Table Hyperparameter ranges used in M4, Crypto, and Player Trajectory datasets.
Learning rate
Integral layer channels
Hidden
Look-back window length
Prediction window
ReVIN
units in
kernel
length
[1e-3, 5e-3]
[16, 32]
[32, 64)
3 to 192
1 to 18
[True, False]
24

---

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
Table 15. Comparison of RMSE from KRNO method with other baseline methods on Crypto and Player Trajectory datasets.
Crypto (Weighted RMSE 10
3)
Basketball Player Trajectory (RMSE)
Model
(1~5)
(6~~10)
(1I~15) Total
(1~10)
(11~200
(21~30)
Total
VARIMA
6.09
8.83
10.74
8.76+000
0.22
0.90=
1.98
1.26=0.0
MLP
6.68
7.95
8.64 8.64
7.85
0.73
1.64
2.77
1.91 +0.32
ML.PRevIN+TB
5.03-1008
7. 7.16
2006
7.012:000
0.37
1.162am
2.25
1.48
RF+TB
6.62
7.99
8.51
7.841:001
0.8610011
2.10
3.48
2.40
FedFormer
5.61
7.50+a0s
8.89
7.46#004
0.43
0.92
1.97
1.29 20.0s
LEM
5.27=0m
7.23
8.23+0005
7. 7.02::004
0.33
1.08
2.18=0.m
1.42
VBS
15.23=000
14.4620m1
26.4940.01
19.52
0.90=n00
2.84
9.24=0000
5.6 5.60=0.00
KNF
5.24
7.03:001
7.63:1001
6.91
0.26:20.1
0.84
1.81:000
1.1620.0
KRNO
5.27
7.07=0.1
7.72+001
6.95 #0.16
0.27
0.93+aa.s
1.94+0007
1.25
H.S. Spatial modeling problems
H.5.1. DARCY PROBLEM
Input
Output (3)
Pred (0)
1 0
1.0
0.006
0.6
0.004
0.002
0.2
0.0
Input
Output (y)
Pred (g)
ly
9
1000
800
40
600
400
20
200
0
Figure The top row presents sample prediction from test set for the Darcy-flow problem, while bottom row illustrates sample prediction for the elasticity problem.
We consider steady-state Darcy flow in a two-dimensional unit square domain S E [0, 1] x [0, 1]. governed by the following PDE
(a(x)Vu(x)) = J, E Ω where u(x)
an
(23)
a(x) is the permeability field, u(x) is the pressure field, and is the source term. The goal here is to learn the pressure field
u(x) as a function of the permeability field axx). In this case, the permeability field is modelled as a(x) exp Fx)), where F(x) denotes a truncated Karhunen-Loéve (KL) expansion with 100 terms. We utilized the data provided (Lu et al., 2022) for training and testing purposes. The training dataset consists 1000 snapshots, and the testing data consists of 200 snapshots. The snapshots corresponding both the input permeability field a(x) and target pressure field u(x) are
25

---

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
given on a uniform grid of size 20 x 20. We compared the relative errors of KRNO with other operator-based methods such DeepONet (Lu et al., 2021) and FNO (Li et al., 2020b). The results are presented in Table 2 and Figure 11.
The KRNO architecture used in the experiment has 3 kernel integral layers with 20 channels each, lifting and projection layers are parametrized by MLPs with one hidden layer containing 20 and 128 hidden units. respectively, and the kernels in the integral transform layers are parametrized by MLPs with hidden layers with 128 hidden units each. We used the AdamW optimizer for training and the learning rate is set to 10-3
H.5.2. HYPER-RLASTE PROBLEM
In this problem, we consider a hyper-elastic material within a unit domain Q E [0, x [0, that contains void of arbitrary shape at the center, as illustrated in the bottom row of Figure 11 (Li et al., 2023). The bottom edges of the unit cell are clamped, and we apply tension traction of = [0, 100] to the top edge. The prior of the void radius is r = 0.2 + 1+exp(/) with F ~ (0, 42 (V + +32) - which satisfies the constraint r E [0.2, 0.4]. The dynamics of an clastic material are represented by the following PDE
0.2
8'u u 012
0
(24)
where is the mass density, u is the displacement vector, and o is the stress tensor. Constitutive models are essential to define the behaviour of the system entirely. These models establish the relationship between the strain e and stress tensors, effectively completing the mathematical framework necessary for solving the system. In this context, the material behaviour is characterized by the incompressible Rivlin-Saunders model, which uses energy density function parameters C1 1.863 10t and C2 9.79 x 10², We utilized the data provided in (Li et al., 2023) for both training and testing purposes. The training dataset consists of 1000 snapshots, and the testing data consists of 200 snapshots generated using finite element solver.
The goal here is to learn the stress field as function of the deformed mesh that corresponds to the void shape. We compared the relative test errors corresponding to predictions of our method with other ML-based operator methods such as FNO (Li et al., 2020b) and DeepONet (Lu et al., 2021). The results are presented in Table 2.
The KRNO architecture used in the experiment has 3 kernel integral layers with 20 channels each, lifting and projection layers are parametrized by MLPs with one hidden layer containing 128 hidden units, and the kernels in the integral transform layers are parametrized by MLPs with 3 hidden layers with 64 hidden units each. We used the AdamW optimizer for training and the learning rate is set to 10
H.6. Shallow water simulation
We provide some additional numerical results for the shallow water test case using the default KRNO architecture detailed in Section H. We followed the procedure described in Kovachki et al. (2023) in our numerical studies using FNO. Table 16 and Figure 12 compares results from different models trained for 200 epochs. The results presented show that FNO- 3D performance improves when trained for 200 epochs. We also applied FNO-2D (Kovachki et al., 2023) to learn an autoregressive model that maps the spatio-tempore field at five time instants to the next time step. However, irrespective of the choice of hyperparameters, we were unable to learn an autoregressive model that provided stable predictions over the testing horizon. Due to this numerical issue, the FNO-2D model is only trained for 100 epochs. Representative predictions from all the models are shown in Figures 14, 15, 16, 17.
H.7. Climate modeling example
The KRNO architecture used in the experiment has 3 kernel integral layers with 20 channels each, lifting and projection layers are parametrized by MLPs with one hidden layer containing 20 and 256 hidden units, respectively, and the kemnlss in the integral transform layers are parametrized by MLPs with 2 hidden layers with 64 hidden units each. We used the Adam W optimizer for training and the learning rate is set to 10-3
26

---

Shifting Time: Time-series Forecasting with Khatri- Neural Operators
Equispaced
Irregular traj 3
Irregular traj 7
true (train)
trap (troin)
true (tralo)
(val)
(val)
(val)
true (tost)
tripel (teat)
true (test)
prod (teet)
powd (teat)
pred (seet)
KRNO
10
10
(train) teue (val)
(train)
(tzais)
(val)
crne (yad)
trie (test)
trie (tewt)
Erue (teat)
(best)
(test)
(teat)
Latent- -DE
teur (train)
tron (train)
true (train)
true (vol)
Stoe (vol)
crue (val)
tiue (test)
(Scat)
(test)
prod (tess)
pred (test)
pred (soat)
T-PATCHGNN
10
Figure 8. Comparison of forecasted predictions of KRNO with other models trained on short trajectories with equispaced and high levels irregularity the time-stamps.
27

---

Shifting Time: Time-series Forecasting with Khatri-Ran Rao Neural Operators
Equispaced
Irregular traj 3
Irregular traj 7
(train) true (val)
(train)
true (val)
true [val]
(ceed)
10
prod (toek)
10
peed (tent)
pred (teat)
KRNO
-5
-5
-1
-20
-20
20
true (traim)
20
troe (train)
true (train)
(val) trile (tead)
(val)
trun (teet)
crue (tawt)
pred [teat)
10
(teet)
Latent- -DD
-5
-5
-6
200
-20
-20
-10
-10
10
20
-10
10
20
(treia)
tsuc (srain)
(train)
trne (val)
15
tree (val)
15
(val)
tede (se28)
prod [teat)
pred (fest)
pred (cest)
T-PATCHGNN
10
-20
-10
10
-10
10
20
-10
10
20
Figure Comparison forecasted predictions KRNO with other models trained on long trajectories with equispaced and high levels of irregularity in the time-stamps.
28

---

Shifting Time: Time-series Forecasting with Khatri-Ran Neural Operators
AirPassengers
AusBeet
GasRateCO2
MonthlyMill
Mgp 1 // /
sunspots
Wine
Wooly
HeartRate
A/0 mipd
Ground Truth
TS-KRNO
Figure 10. KRNO predictions the Darts datasets
KRNO PXG LOCA
10-3
Table 16. Comparison of the average relative T.? errors the shallow water problem for the three field vari- ables when training is conducted for 200 epochs.
10-2
10-2
Method
L relative error
10-2
P
u
0.1
0.2 0-3
0.4
0.5
0.6
FNO-3D
0.000719
0.01951
0.01174
LOCA
0.003091
0.15179
0.14942
KRNO
0.000331
0.01339
0.01406
Figure 12. Comparison of the average relative L errors function of time for the three field variables (across the 1000 test simulations) obtained using KRNO, FNO-3D and LOCA models trained for 200 epochs.
29

---

Shifting Time: Time-series Forecasting with Khatri-Rao Neural Operators
time (s)
0.06
0.12
0.18
0.24
0.30
0.36
0.42
0.48
0.54
0.60
P
1.04
1.55
1.000
1.00
1.0
10000
1
1.050
0.025
1002
00O
100
1
2
1p-
P
0..
0.0
0.
10000
100000
000000
20000
100%0
00000
10000
20000
1000000
0
101000
10000
1,00
200000
0200
0000000
00022200
000000
0
0
00.00
2
i
- 55
0.0
0.8
25
325
0.11
000
000
10.00
1000
0.00
00
1000
3..00
20020
2000
0.033
100000
0.020
0
0.00
0010
0.00
2000
4000
2000
200000
101000
100000
600%
000'0
10000
0000
55000 0
00000 o
$200000
10000:00 0.
100'0
200'0
2000'0
000000
-00000
010000
10000
0
2
2
0.3
0.0
2000
00.
0000
0000
0.1
000
10:0
100
4000
008
0.000
0.023
0.000
1023
1000
1022
00
18
0.00
100
11000
2200'0
00000
00000
-0100 0
$900 '
20000
+0000
10000
20000
*000'0'
000000
*20000
000000 O
0-0007:3
10000
1200
1000000 0
1000000
00000'
000025' 0
N090000 o
0
Figure Shallow water problem: Figure shows the predictions (p, ú, the three field variables along with the fields (p. u, v) function of time for test simulation using KRNO trained for 100 epochs.
30

---

Shifting Time: Time-series Forecasting with Khatri- Rao Neural Operators
p
e
P
0.00
0.15
0.30
0.45
0.60
(sec)
rel. 2 error
0.0
0.1
0.2
0.3
0.4
0.5
0.6
t (sec)
Figure 14. Shallow water problem: Top figure shows the predictions (p, à, 6) for the three field variables along with the true fields (p, v) function of time for simulation using FNO-2D trained for 100 epochs. The bottom figure shows the emror bars representing the 2 relative errors for three field variables across the 1000 test simulations, with the shaded region indicating + standard deviation.
31

---

Shifting Time: Time-series Forecasting with Khatri-Ran Rao Neural Operators
p
0.00
0.15
0.30
0.45
0.60
t (sec)
rel. 2 eгтоr
0.1
2
0.0
0.0
0.1
0.2
0.3
0.4
0.5
0.6
(scc)
Figure Shallow water problem: Top figure shows predictions (ő, for the variables along with the true fields (p. u, v) function of time for test simulation using KRNO trained for 200 epochs. The bottom figure shows the error hars representing the 1,² relative errors for three variables across 1000 test simulations, the shaded region indicating standard deviation.
32

---

Shifting Time: Time-series Forecasting with Khatri- Rao Neural Operators
p
0.00
0.15
0.30
0.45
0.60
(sec)
rel. 2 eгтоι
77
27
0.0
0.0
0.1
0.2
0.3
0.4
0.5
0.6
(scc)
Figure Shallow water problem: Top figure shows the predictions (p, for the variables along with the fields (p, u, v) function of time for test simulation using FNO-3D trained for 200 epochs. The bottom figure shows the error bars representing the I? relative errors for three field variables across the 1000 test simulations, with the shaded region indicating + standard deviation.
33

---

Shifting Time: Time-series Forecasting with Khatri-Ran Rao Neural Operators
p
P
0.00
0.15
0.30
0.45
0.60
(sec)
rel. 2 error
0.0
0.1
0.2
0.3
0.4
0.5
0.6
t (sec)
Figure 17. Shallow water problem: Top figure shows the predictions (p, à, 6) for the three field variables along with the true fields (p, v) function of time for simulation using LOCA model trained for 200 epochs. The bottom figure shows the error bars representing the L relative errors for three field variables across the 1000 test simulations, with the shaded region indicating + standard deviation.
34


Reviews:
[{"review_id": "D9qvTyFWbC", "reviewer": "Reviewer_P25q", "summary": "This paper presents a method for time-series forecasting that treats the task as learning a continuous time-shift operator, approximated via a proposed architecture called Khatri-Rao Neural Operator (KRNO). The operator is modeled as an integral transform with a non-stationary kernel, decomposed via Khatri-Rao product structure. The resulting model enables forecasting with irregularly sampled data, super-resolution, and low parameter count. The method is validated across a large suite of temporal and spatio-temporal benchmarks.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your feedback and comments.\n\n### 1. Evaluation on Irregularly Sampled Datasets\nWe agree that numerical studies on regularly sampled time-series do not allow a clear demonstration of the full capabilities offered by KRNO. We now include additional experiments on challenging irregularly sampled time series benchmark datasets (MIMIC, USHCN, Human Activity, and MuJoCo) and the results are compared against a range of alternative approaches such as T-PatchGNN, NeuralSDE, NeuralCDE, PatchTST, and Latent-ODE. KRNO achieves SOTA performance on three of the four benchmarks. Please see Tables in section *Baselines and Evaluation Criteria* of our response to reviewers **ntzg** and **kNx6**.\nAs suggested, we have updated the Appendix to include more details on the Darcy-flow and hyper-elasticity benchmarks. \n\n### 2.  Model Architecture & Implementation Details \nKRNO's architecture follows the general structure of neural operators (lines 129-150), i.e., lifting layers followed by series of kernel integral layers and a projection layer. The key difference is our use of a non-stationary product structured kernel in the kernel integral transform (Equation 6). The KRNO representation of the continuous time-shift operator can be written as:\n$$\n\\mathcal{A}_{t_p}^{t,t_f}:= \\mathcal{P} \\circ \\mathcal{K}_n \\circ \\ldots \\circ \\mathcal{K}_1 \\circ \\mathcal{L},\n$$\nwhere $\\mathcal{K}_i: \\mathbb{R}^{p} \\to \\mathbb{R}^{q}$ is a kernel integral layer defined in Equation 5, $\\mathcal{L}: \\mathbb{R}^{n} \\to \\mathbb{R}^{c}$, and $\\mathcal{P}: \\mathbb{R}^{c} \\to \\mathbb{R}^{n}$ are the pointwise lifting and projection layers parameterized by neural networks with one hidden layer, and $c$ is the number of channels in integral layers.\n\nTo ensure reproducibility, we have shared the source code and scripts used in our studies at [https://anonymous.4open.science/r/KRNO-1F4F/](https://anonymous.4open.science/r/KRNO-1F4F/).\n\n### 3. Theoretical Consistency\nPlease note that the matrix-valued kernels in KRNO are consistently defined as functions of time and space throughout the paper and in our implementation. The kernel integral layers in Equation 5 use these parametrized space-time kernels to learn mappings between function spaces. As a consequence, the numerical implementation is aligned with the operator-theoretic formulation. \n\n### 4. Ablation Studies\nThank you for suggesting insightful ablation studies. We have conducted ablation studies as suggested to better understand the key components that influence generalization. One of the key parameters is the number of integral layers/channels. Larger values of this parameter enhance model capacity at the expense of increased memory usage for high-resolution data. This trade-off can be effectively managed by adjusting the number of quadrature points in the kernel integral layers. The following ablation study on the Elasticity problem (learning stress fields from void deformation in elastic blocks) confirms this flexibility. Using training data given on a 41×41 grid, we notice that reducing latent grid resolution in the integral layers significantly decreases memory requirements with minimal impact on accuracy. We intend to include additional ablation studies to illustrate the impact of the quadrature rule on generalization performance. \n|Latent Grid Resolution|L2 Rel. error|GPU Memory (MB)|\n|:---:|:---:|:---:|\n|16 x 16 |5.20 ± 0.18 % |1,386 |\n|24 x 24 |5.14 ± 0.17 % |1,774 |\n|32 x 32 |5.14 ± 0.26 % |2,186 |\n|40 x 40 |5.12 ± 0.15 % |2,650 |\n|48 x 48 |5.16 ± 0.16 % |3,174 |\n\n### 5. Handling Measurement Noise and Exogenous Inputs\nThe new irregularly sampled benchmarks mentioned earlier are challenging with missing and noisy observations. Our results show that KRNO effectively handles such data. In the present work, we formulate the time-shift operator in the setting of deterministic ordinary/partial differential equations. This formalism enabled sufficient flexibility to provide strong generalization across temporal and spatio-temporal forecasting problems from diverse domains.\n\nFor example, for the case of time-series datasets, the forced time-shift operator can be defined as $X\\_{(t,t_f]}=\\mathcal{A}\\_{t_p}^{t,t_f} (X_{[t_p,t]}, f_{[t_p,t_f]})$, where $X_{[t_p,t]}$ denotes the state trajectory over the time-interval $[t_p,t]$ while $f_{[t_p,t_f]}$ denotes the forcing function over the input and forecasting time-interval $[t_p,t_f]$. To ensure that the time-shift operator is causal (i.e., the predictions made at time instant $\\tau$ are not influenced by the future values of the forcing function), the integration domain of the kernel integral transform w.r.t. $f$ should be set to $[t_p, \\tau]$, where $\\tau \\in (t, t_f]$. A similar approach can be used to define the time-shift operator for forced spatio-temporal dynamical systems. These extensions would enable our approach to be applied to systems with exogenous effects.", "claims_and_evidence": "The main claim of the authors is that they introduce the continuous time-shift operator. There is a drastic disconnect between this claim and the empirical studies. None of the studies demonstrate the importance of continuous time-shift as concept. The continuos shift in either time or space can be achieved by a simple MLP, mapping time and coordinates along with some additional features into predictions. Why do we need KRNO, what is special about it that other architectures cannot do? This is further emphasized by the fact that large part of experiments are conducted on regularly sampled non-spatial datasets that have so many methods that work extremely well on them, that the authors will not have enough space in a 20-page paper to review, explain and present their results. May I suggest that the authors rather focus on a very concrete problem and drill it down, instead of overwhelming the reader with results and studies that are largely irrelevant from the point of view of the problem that they claim to solve?", "methods_and_evaluation_criteria": "- I am not quite convinced that the use of regularly sampled forecasting datasets, such as Darts, Crypto, Baseball or M4, is relevant, as the proposed theory and operator-based approach seem to be most applicable to irregularly spaced and spatio-temporal data. In my view, these results overload the paper and distract attention from other important topics. For example, there is very little discussion of what the proposed model learns i.e., how the time-shift operator behaves or generalizes across inputs. Similarly, there is very little discussion in the text of the implementation details. How does the inference path of the model look like, exactly? Can you provide equations describing it? Additionally, there are no ablation studies.", "theoretical_claims": "- Apparent theory/practice disconnect: the transition between theory (quadrature approximations, kernel formulation) and implementation (data-driven neural parametrization) feels abrupt (i.e. Line 213). The theory defines kernels as functions of coordinates, but the practical implementation uses learnable functions of data, raising questions about the operator-theoretic validity of the implementation.", "experimental_designs_or_analyses": "- The description of many experiments lacks clarity. For example, for Darcy-flow and hyper-elastic benchmark no details are provided at all. What are the problems, why are they relevant for the evaluation of the proposed algorithm? What are dataset/problem sizes, what are the splits? How is the L2-relative error computed, exactly? Similar remarks apply to most benchmarks used in the paper.\n\n- The relevance of Figure 2 is not clear. It would constitute a much stronger case if it contained a comparison with baseline methods that significantly underperform on the presented anecdotes. The fact that the proposed technique apparently does reasonably well on these cases does not render convincing evidence of the proposed framework being able to solve problems that other approaches fail to provide adequate solutions for.\n\n- The paper contains no ablation studies whatsoever. This is a red flag. What are important model components, key assumptions and their effects on model performance?\n\n- Missing comparisons to neural controlled differential equations (CDEs), which are designed for irregular and continuous-time forecasting, are only discussed briefly and not benchmarked empirically.", "relation_to_broader_scientific_literature": "- Overuse of standard mathematical machinery: Much of the theoretical exposition leans heavily on classical results (e.g., Grönwall’s inequality, semigroup continuity), which could have been cited more concisely. Similarly, Proposition 2.1 is a simple consequence of the Khatri–Rao product structure and it does not need a formal statement in the form of proposition or a proof. Its inclusion, while technically correct, adds bulk without substantive new theory.", "essential_references_not_discussed": "[1] Horn et al., Set Functions for Time Series https://arxiv.org/pdf/1909.12064\nOffers a non-sequential way to model irregular time series using permutation-invariant functions.\n\n[2] Kidger et al., Neural Controlled Differential Equations for Irregular Time Series, https://proceedings.neurips.cc/paper/2020/file/4a5876b450b45371f6cfe5047ac8cd45-Paper.pdf, \nmemory-efficient ODE-based approach to irregular time-series", "other_strengths_and_weaknesses": "NA", "other_comments_or_suggestions": "NA", "questions_for_authors": "- Can you clarify the distinction between kernels being functions of time/space vs. functions of the data? Does this violate the assumptions made in the integral operator formulation?\n\n- Why are neural CDEs or latent ODEs not included in the empirical comparisons, especially given the shared goal of handling irregular sampling?\n\n- Can you add ablation studies showing the effects of important model components on model performance? As an example, could you provide an ablation of the quadrature rule used (midpoint vs. trapezoidal vs. learned)? What are other important components and their effects on the model accuracy/speed/memory?\n\n- Is the method capable of forecasting in the presence of measurement noise or exogenous inputs?", "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": "jhYI4W63rr", "reviewer": "Reviewer_ntzg", "summary": "The paper introduces a novel operator-theoretic framework for time-series forecasting, leveraging the Khatri-Rao Neural Operator (KRNO) to learn continuous time-shift operators. By relaxing the discrete lag factor in autoregressive models, KRNO enables super-resolution forecasting in both space and time while handling irregularly sampled observations. The authors demonstrate KRNO’s scalability and competitiveness across benchmark datasets.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your feedback and comments.\n\n## 1. Theoretical Guarantees and Kernel Structure\n\nThe computational complexity of the kernel integral transform layer scales as O($n^2$), where $n$ is the number of quadrature nodes in the input and output domains. Additional assumptions are required to reduce this complexity. For instance, FNOs assume that the kernel is stationary, thereby enabling the use of FFT. In the present work, we overcome this complexity using product structured non-stationary kernels. \n\nThe only approximation error comes from the quadrature scheme used to compute the integral in equation 6; the non-stationary kernel is evaluated exactly, without additional errors. We have added an ablation study to show the influence of the number of quadrature points on the generalization performance of KRNO (please refer to our response to reviewer **P25q** on ablation studies). A particularly attractive feature of KRNO is that it allows us to learn a flexible non-stationary kernel for each component of the product structured kernel. As evidenced by our numerical studies, this additional flexibility provides performance gains over stationary kernel-based neural operator architectures, while significantly reducing the model parameter count.\n\n\n## 2. Implementation and Applicability\n\nIn our exposition of KRNO, we considered the case of product grids to achieve almost linear scalability. It is worth noting that our approach can be extended to unstructured grids by introducing a learnable function that maps the data from an unstructured grid to a latent product grid as in Geo-FNO [1].\n\nWe have added the missing hyper-parameter details for each benchmark in the Appendix. To ensure reproducibility, we have shared the source code and scripts used in our numerical studies at [https://anonymous.4open.science/r/KRNO-1F4F/](https://anonymous.4open.science/r/KRNO-1F4F/).\n\n[1] Li, Zongyi, et al. \"Fourier neural operator with learned deformations for pdes on general geometries.\" JMLR 24.388 (2023): 1-26.\n\n\n## 3. Baselines and Evaluation Criteria:\n\nWe agree that numerical studies on regularly sampled time-series do not allow a clear demonstration of the full capabilities offered by KRNO. We now include additional experiments on challenging irregularly sampled time series benchmark datasets (MIMIC, USHCN, Human Activity, and MuJoCo) and the results are compared against a range of alternative approaches with SOTA performance such as T-PatchGNN [1], NeuralSDE [2], NeuralCDE, PatchTST, and Latent-ODE. KRNO achieves new SOTA performance on three of the four benchmarks; please see the Table below which provides results for three of the datasets. The results for the MuJoCo benchmark are provided in our response to reviewer **kNx6**. \n\n| **Method** | **MIMIC MSE×10⁻²** | **MIMIC MAE×10⁻²** | **USHCN MSE×10⁻¹** | **USHCN MAE×10⁻¹** | **Human Activity MSE×10⁻³** | **Human Activity MAE×10⁻²** |\n|:---:|:---:|:---:|:---:|:---:|:---:|:---:|\n| TimesNet | 5.88 ± 0.08 | 13.62 ± 0.07 | 5.58 ± 0.05 | 3.60 ± 0.04 | 3.12 ± 0.01 | 3.56 ± 0.02 |\n| PatchTST | 3.78 ± 0.03 | 12.43 ±0.10 | 5.75 ± 0.01 | 3.57 ± 0.02 | 4.29 ± 0.14 | 4.80 ± 0.09 |\n| GRU-D | 1.76 ± 0.03 | 7.53 ± 0.09 | 5.54 ± 0.38 | 3.40 ± 0.28 | 2.94 ± 0.05 | 3.51 ± 0.06 | \n| Warpformer | 1.73 ± 0.04 | 7.58 ± 0.13 | 5.25 ± 0.05 | 3.23 ± 0.05 | *2.79 ± 0.04* | *3.39 ± 0.03* |\n| mTAND | 1.85 ± 0.06 | 7.73 ± 0.13 | 5.33 ± 0.05 | 3.26 ± 0.10 | 3.22 ± 0.07 | 3.81 ± 0.07 |\n| Latent-ODE | 1.89 ± 0.19 | 8.11 ± 0.52 | 5.62 ± 0.03 | 3.60 ± 0.12 | 3.34 ± 0.11 | 3.94 ± 0.12 |\n| T-PatchGNN | *1.69 ± 0.03* | **7.22 ±0.09** | *5.00 ± 0.04* | *3.08 ± 0.04* | **2.66 ± 0.03** | **3.15 ± 0.02** |\n| KRNO | **1.57 ± 0.02** | *7.43 ± 0.06* | **4.95 ± 0.08** | **3.06 ± 0.08** | 2.85 ± 0.03 | 3.46 ± 0.02 |\n\n[1] Zhang, Weijia, et al. \"Irregular multivariate time series forecasting: A transformable patching graph neural networks approach.\" ICML 2024.\n\n[2] Oh, Y., Lim, D., & Kim, S. Stable Neural Stochastic Differential Equations in Analyzing Irregular Time Series Data. ICLR 2024.\n\n\n## 4. Computational Cost Analysis\n\nPlease refer to our first response to reviewer **dfXC** for more clarification on the computational complexity and runtime analysis of KRNO compared to FNO.\n\n## 5. Temporal Dependency\n\nYou have raised a very important point about the choice of $t_p$ and $t_f$. In our numerical studies, we treated them as hyperparameters. However, it would be valuable to explore how these parameters can be adaptively chosen and how this could potentially improve generalization. This point is also closely related to the application of KRNO to model long-range dependencies. We plan to pursue these directions in future work.", "claims_and_evidence": "1. Efficiency: KRNO achieves near-linear computational cost for non-stationary integral transforms, scaling better than FNO and other neural operators.\n\n2. Flexibility: KRNO handles irregular sampling and super-resolution forecasting by parametrizing the time-shift operator as a continuous kernel.\n\n3. Superior Performance: KRNO ranks among the top 3 methods on 21/29 test cases and tops 10/29 datasets.\n\nWeakness:\n1. Overstated Generality: The claim that KRNO “inherits the benefits of neural operators” (e.g., discretization independence) overlooks its reliance on specific kernel structures (Equation 7).\n\n2. Lack of Theoretical Guarantees: While Proposition 2.1 supports computational efficiency, there is no theoretical analysis of approximation error or stability for non-product-kernel scenarios.", "methods_and_evaluation_criteria": "1. KRNO is built on non-stationary integral transforms with component-wise kernels decomposed as Khatri-Rao products (Equation 7).\n\n2. The architecture includes lifting/projection layers and three kernel integral transform layers, with neural networks parametrizing each kernel.\n\nWeakness:\n1. Assumption of Product Grids: KRNO’s near-linear complexity relies on input/output data lying on product grids (e.g., time × latitude × longitude). This limits applicability to irregularly gridded or high-dimensional data.\n\n2. Fixed Hyperparameters: The time-shift operator’s boundaries (t_p, t_f) are treated as hyperparameters, but the paper does not discuss adaptive strategies for varying forecasting horizons.\n\n3. Ignoring Temporal Dependencies: The operator learns mappings over fixed windows ([t_p, t] to (t, t_f]), neglecting long-range dependencies beyond the window size.", "theoretical_claims": "There is no theoretical analysis of how approximation errors in the kernel decomposition (Equation 7) propagate to the overall forecasting error.\nWhile KRNO avoids ODE-based adjoints, the paper does not compare its gradient estimation stability to neural ODEs in noisy or high-dimensional settings.", "experimental_designs_or_analyses": "1. Lack of Hyperparameter Tuning Details: Key parameters (e.g., number of layers, hidden units) are not fully documented, reducing reproducibility.\n\n2. Inconsistent Evaluation Protocols: Some experiments (e.g., M4) use recursive forecasting, while others (e.g., Darts) rely on fixed window sizes, complicating comparisons.\n\n3. Missing Baselines: Notable omissions include recent SOTA methods like TimesNet (Wu et al., 2022) and LogTrans (Li et al., 2023) for time-series forecasting.\n\n4. Computational Cost Analysis: While Table 7 compares GPU memory, training/inference times are only shown for shallow water (Figure 9), lacking scalability analysis for larger datasets (e.g., M4).", "relation_to_broader_scientific_literature": "Connects KRNO to operator-learning frameworks (DeepONet, FNO) and highlights advantages over autoregressive models (e.g., Transformer, N-BEATS).", "essential_references_not_discussed": "No", "other_strengths_and_weaknesses": "See review above", "other_comments_or_suggestions": "1. Clarify the scope of KRNO’s applicability (e.g., product grids vs. arbitrary grids).\n\n2. Add error bounds or stability analysis for non-stationary kernels.\n\n3. Include recent SOTA methods (e.g., TimesNet, LogTrans) in benchmarks.\n\n4. Provide full hyperparameter details and reproducibility protocols.", "questions_for_authors": "See review above", "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": "lLP0z7JlrB", "reviewer": "Reviewer_dfXC", "summary": "The paper introduces a novel operator-theoretic approach for time-series forecasting by learning a continuous time-shift operator. This method provides a more flexible alternative to traditional autoregressive models, which rely on discrete time lags. The authors propose Khatri-Rao Neural Operators (KRNOs) as a new architecture to parametrize non-stationary integral transforms, enabling efficient learning of time-dependent dynamics in both temporal and spatio-temporal forecasting.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your feedback and comments.\n\n\n## 1. Computational Complexity and Runtime Analysis\n\nIn Appendix G, we provide a detailed comparison of computational complexity and runtime between KRNO and FNO-3D using the spatio-temporal shallow water problem. While our initial analysis used the default 12 Fourier modes for FNO-3D across all spatial resolutions, we note that when increasing the number of modes for high-resolution datasets, FNO-3D's memory usage and training time exceed KRNO's requirements. The table below presents an updated runtime analysis comparing KRNO with FNO-3D across different spatial resolutions (S × S) using the maximum number of Fourier modes (S/2+1). These results demonstrate that KRNO is both faster and more memory-efficient than FNO-3D at higher resolutions.\n\n| Spatial resolution (S x S) | Memory (MB) | | | | Time (seconds) | | | |\n|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|\n| | Training | | Testing | | Training | | Testing | |\n| | **KRNO** | FNO-3D| **KRNO** | FNO-3D | **KRNO** | FNO-3D | **KRNO** | FNO-3D |\n| 32 × 32 | 1,390 | 1,074 | 708 | 1,074 | 0.0279 | 0.0234 | 0.0107 | 0.0062 |\n| 64 × 64 | 2,776 | 2,772 | 1,314 | 2,772 | 0.0394 | 0.0347 | 0.0134 | 0.0071 |\n| 96 × 96 | 4,884 | 5,264 | 2,366 | 5,264 | 0.0626 | 0.0710 | 0.0185 | 0.0149 |\n| 128 × 128 | 7,608 | 8,994 | 3,796 | 8,994 | 0.0999 | 0.1297 | 0.0312 | 0.0288 |\n| 160 × 160 | 10,040 | 13,764 | 5,644 | 11,288 | 0.1584 | 0.2088 | 0.0502 | 0.0456 |\n\n## 2. Handling Missing or Noisy Data\n\nWe now include additional experiments on challenging irregularly sampled time series benchmark datasets (MIMIC, USHCN, Human Activity, and MuJoCo), which contain missing and noisy observations. The performance of KRNO is compared against a range of alternative approaches with SOTA performance such as T-PatchGNN [1], NeuralSDE [2], NeuralCDE, PatchTST, and Latent-ODE. \n\nKRNO achieves new SOTA performance on three of the four benchmarks; please see the tables in the *Baselines and Evaluation Criteria* section of our responses to reviewers **ntzg** and **kNx6**. From these studies, we see that KRNO is able to handle missing and noisy data. In our future work, we will extend KRNO to support missing data for the case of missing data in spatio-temporal problems.\n\n[1] Zhang, Weijia, et al. \"Irregular multivariate time series forecasting: A transformable patching graph neural networks approach.\" ICML 2024.\n\n[2] Oh, Y., Lim, D., & Kim, S. Stable Neural Stochastic Differential Equations in Analyzing Irregular Time Series Data. ICLR 2024.\n\n## 3. Generalization to High-Dimensional Spatio-Temporal Problems\n\nNumerical studies on spatial modeling problems, 2D spatio-temporal forecasting problems, and temporal forecasting problems show that the KRNO architecture provides better performance than competing methods for this class of problems such as FNO, DeepONet, and LOCA. We plan to evaluate KRNO on challenging 3D spatio-temporal problems in future work.", "claims_and_evidence": "The paper presents several claims, most of which are backed by clear empirical results.", "methods_and_evaluation_criteria": "The paper evaluates KRNO on 29 different forecasting tasks, covering applications in climate modeling, financial markets, and fluid dynamics. The model is tested against leading baselines, including Fourier Neural Operators (FNO), DeepONet, and traditional autoregressive models", "theoretical_claims": "The proofs for Proposition 2.1 and its generalization appear correct in principle and are logically structured, but they lack some explicit derivations and justifications for key claims (especially computational efficiency). While these omissions do not invalidate the results, addressing them would increase clarity.", "experimental_designs_or_analyses": "The paper evaluates KRNO on several temporal and spatio-temporal datasets, including Darts, M4, and physics-based problems (Darcy flow, hyper-elastic problems). It compares KRNO with FNO, DeepONet, and other neural operators. It reports relative error as the primary evaluation metric. The paper claims that KRNO has almost linear computational complexity. However, there is no explicit runtime comparison with other baselines, making this claim difficult to verify.", "relation_to_broader_scientific_literature": "The continuous time-shift operator proposed in this paper generalizes discrete autoregressive models by modeling entire function trajectories instead of relying on discrete time steps. This idea extends prior work on neural operators like Fourier Neural Operators and DeepONet", "essential_references_not_discussed": "No", "other_strengths_and_weaknesses": "The introduction of Khatri-Rao decompositions enhances computational efficiency compared to FNO and DeepONet, making it feasible for large-scale forecasting. The claim of near-linear complexity is promising, though a more explicit runtime analysis would strengthen this claim. The model is evaluated on 29 forecasting tasks, including climate modeling, physics simulations, showing broad applicability.\n\nThe claim that KRNO achieves near-linear complexity is reasonable based on Khatri-Rao properties, but the paper does not provide wall-clock runtime comparisons with FNO, DeepONet, and transformer-based methods and Scalability analysis for increasing dataset sizes and higher-dimensional problems.", "other_comments_or_suggestions": "N\\A", "questions_for_authors": "1. Can you provide runtime benchmarks comparing KRNO with FNO, DeepONet, and standard transformer models?\n2. How does KRNO handle missing or noisy data?\n3. Does KRNO generalize well to very high-dimensional spatio-temporal problems, such as 3D weather forecasting or turbulence modeling?", "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": "QjaxmeJSwY", "reviewer": "Reviewer_kNx6", "summary": "The authors propose a time series forecasting method, leveraging continuous time-shift operators, which act as continuous analogs of the lag factor of the discrete-time autoregressive models or the upsampling/downsampling layers in CNNs, mapping the history of values up to an observation to a future window. To tackle the operator-learning problem, the authors introduce Khatri-Rao neural operators (KRNO) that define non-stationary integral-based transforms with the almost linear cost for spatio-temporal problems. The proposed method is evaluated on several real-world and synthetic spatiotemporal and standard temporal forecasting datasets, showcasing competitive performance with other popular time series methods. An anonymized code repository is also provided for reproducibility purposes.\n\n## update after rebuttal\n\nAssessing the overall impact of the contribution, including the theoretical results and experimental evaluation, and following the author's rebuttal that showcased incremental performance results on real-world irregularly sampled datasets (or without proven statistical significance), I maintain my initial ratings. During the rebuttal, the authors have not adequately addressed significant questions concerning the positioning of the contribution (W1, W2) or presented a thorough computational (e.g., time cost) analysis between the proposed method and considered baselines for the chosen tasks (which is a significant aspect in case of borderline performance improvements). Performance percentage improvements in terms of MSE are misleading (referring to the 2nd decimal with overlapping stds), and the proposed method is outperformed in Tables 5 and 6, while results during the rebuttal on irregular datasets showcase average improvements ranging from 1% to 7% for 2 out of the 3 datasets.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the detailed review and valuable comments. \n\n## 1. Claims and Evidence:\n>**(Almost) Linear Complexity..**\n\nPlease refer to our first response to reviewer **dfXC** for more clarification on the computational complexity and runtime analysis of KRNO compared to FNO.\n\n>**Continuous Analog..**\n\nUsing the proposed continuous time-shift operator, we learn an operator that maps the history of the dynamics over a past time-window into its future values over a subsequent time-window. This enables us to learn from irregularly sampled observations and to forecast at any given time over a fixed time-window, making it a continuous analog of the lag factor of the discrete-time autoregressive models. It is also worth noting that the continuity constant of the time-shift operator can be bounded in terms of the Lipschitz constants of the integral transform layers. We plan to study this in future work to show robustness to distribution shifts similar to the seminal work of Oh et al. (2024) on stable NSDEs. \n\n\n## 2. Baselines and Evaluation Criteria:\nWe agree that numerical studies on regularly sampled time-series do not allow a clear demonstration of the full capabilities offered by KRNO. We now include additional experiments on challenging irregularly-sampled time series benchmark datasets (MIMIC, USHCN, Human Activity, and MuJoCo) and the results are compared against a range of alternative approaches with SOTA performance such as T-PatchGNN, NeuralSDE, NeuralCDE, PatchTST, and Latent-ODE. KRNO achieves new SOTA performance on three of the four benchmarks. \n\nThe results for the MuJoCo benchmark are provided below where we compare against the methods in the ICLR 2024 paper by Oh et al [1] which you brought to our attention. It can be seen that KRNO achieves top performance on this benchmark on all four cases. Scripts for reproducing the results for the additional four benchmarks are provided in the anonymous repository under directories 'code/scripts/mujoco' and 'code/scripts/irregular_time_series/krno'.\n\nThe results for MIMIC, USHCN, and Human Activity can be found in our response to reviewer **ntzg**. \n\n|Methods|Regular|30% dropped|50% dropped|70% dropped|\n|:---|:---:|:---:|:---:|:---:|\n|GRU-Δt|0.223 ± 0.020|0.198 ± 0.036|0.193 ± 0.015|0.196 ± 0.028|\n|GRU-D|0.578 ± 0.042|0.608 ± 0.032|0.587 ± 0.039|0.579 ± 0.052|\n|GRU-ODE|0.856 ± 0.016|0.857 ± 0.015|0.852 ± 0.015|0.861 ± 0.015|\n|ODE-RNN|0.328 ± 0.225|0.274 ± 0.213|0.237 ± 0.110|0.267 ± 0.217|\n|Latent-ODE|0.029 ± 0.011|0.056 ± 0.001|0.055 ± 0.004|0.058 ± 0.003|\n|Augmented-ODE|0.055 ± 0.004|0.056 ± 0.004|0.057 ± 0.005|0.057 ± 0.005|\n|ACE-NODE|0.039 ± 0.003|0.053 ± 0.007|0.053 ± 0.005|0.052 ± 0.006|\n|NCDE|0.028 ± 0.002|0.027 ± 0.000|0.027 ± 0.001|0.026 ± 0.001|\n|ANCDE|0.026 ± 0.001|0.025 ± 0.001|0.025 ± 0.001|0.024 ± 0.001|\n|EXIT|0.026 ± 0.000|0.025 ± 0.004|0.026 ± 0.000|0.026 ± 0.001|\n|LEAP|0.022 ± 0.002|0.022 ± 0.001|0.022 ± 0.002|0.022 ± 0.001|\n|Neural SDE|0.028 ± 0.004|0.029 ± 0.001|0.029 ± 0.001|0.027 ± 0.000|\n|Neural LSDE|0.013 ± 0.000|0.014 ± 0.001|0.014 ± 0.000|*0.013 ± 0.001*|\n|Neural LNSDE|*0.012 ± 0.001*|0.014 ± 0.001|0.014 ± 0.001|0.014 ± 0.000|\n|Neural GSDE|0.013 ± 0.001|*0.013 ± 0.001*|*0.013 ± 0.000*|0.014 ± 0.000|\n|KRNO|**0.007 ± 0.002**|**0.008 ± 0.002**|**0.011 ± 0.004**|**0.012 ± 0.002**|\n\n\n[1] Oh, Y., Lim, D., & Kim, S. Stable Neural Stochastic Differential Equations in Analyzing Irregular Time Series Data. ICLR 2024.\n\n## 3. Questions and Weaknesses:\n**\\[Q1\\]** Thank you for your suggestion. We have restructured the methods section to highlight our contributions more clearly. We have also updated the related work section by discussing  SOTA methods such as Neural SDEs, T-PatchGNN, NeuralCDE, and Latent-ODE for irregularly sampled time-series forecasting. \n\n**\\[Q2\\].** Please refer to our response in section *Baselines and Evaluation Criteria* regarding comparison with SOTA methods and additional benchmarks on irregularly sampled time series datasets with missing observations. For these additional benchmarks, we have now added the results from multiple runs with different random seeds. It can be seen from the results that KRNO achieves new SOTA performance on three of the four new benchmarks we studied.\n\n**\\[Q3\\]** Neural operators have been primarily successful on data from physical systems governed by PDEs. Our comprehensive experiments demonstrate that KRNO generalizes effectively to both physical systems (MuJoCo) and non-physical datasets (healthcare, climate, human activity) with irregular sampling patterns.\n\n**\\[W6\\]** We have compared the memory usage and runtime analysis of KRNO with Neural GSDE for the MuJoCo dataset. For a batch size of 128, KRNO uses 850MB of memory while Neural GSDE uses 306MB of memory. On the GTX4090 GPU, for this setting, we found that the default KRNO configuration is around 10 times faster than Neural GSDE per iteration.", "claims_and_evidence": "A few problematic claims:\n- *(Almost) Linear Complexity Competitive to Operator-Based Methods:* The almost linear complexity of the method compared to other operator-based learning frameworks is not showcased in the main paper but rather in the appendix, where experimental results show that the proposed KRNO method is the worst in terms of train-test times per iteration for the spatiotemporal shallow water dataset (surpassed in several cases by the FNO method that has O(nlogn) complexity).\n- *Continuous Analog of Discrete-Time Autoregressive Models:* It is not very clear how the proposed KRNO method is a continuous analog of the lag factor of the discrete-time autoregressive models, and for which selection of kernels this holds.", "methods_and_evaluation_criteria": "Proposed methods and evaluation criteria (e.g., datasets) are generally appropriate for the task at hand but could be significantly enhanced to enable a thorough experimental evaluation.\n\n**Baseline Methods:** Particularly for the case of irregularly sampled time series, only two baselines are used. More recent state-of-the-art methods can be found here [2].\n\n**Datasets:** The proposed method is evaluated on two spatiotemporal forecasting datasets, including data for shallow water simulation and climate modeling, and several temporal forecasting datasets, including the M4 competition archive, the Darts archive, and the two additional crypto and player trajectory datasets. The datasets are commonly used in the (at least for the temporal) forecasting community. Yet several common datasets are missing in terms of regularly sampled time series (ETT, Electricity, Traffic etc, summarized here https://github.com/thuml/Time-Series-Library). The authors also generate an irregularly sampled synthetic dataset to validate their method for data with non-equidistant observation times. The synthetic 2D spiral dataset is rather limiting for the irregularly sampled time series. More prominent methods perform forecasting on apriori irregularly-sampled data or randomly downsample regular ones and do extrapolation (see [1,2]). \n\n[1] Rubanova, Y., Chen, R. T., & Duvenaud, D. K. (2019). Latent ordinary differential equations for irregularly-sampled time series. Advances in neural information processing systems, 32.\n\n[2] Oh, Y., Lim, D. Y., & Kim, S. (2024). Stable neural stochastic differential equations in analyzing irregular time series data. arXiv preprint arXiv:2402.14989.", "theoretical_claims": "Proofs for the properties of the time-shift operator and Proposition 2.1 have been checked for correctness. No issues detected.", "experimental_designs_or_analyses": "The authors follow standard experimental designs for spatiotemporal and temporal forecasting (including optimization, metrics, and datasets). It is unclear if they measure the variance of models' performances with multiple runs with different random seeds if the results refer to a single fixed seed.", "relation_to_broader_scientific_literature": "This work is related to the problem of time series forecasting, including spatiotemporal and temporal settings, also extending to irregularly sampled time series that are very common in several engineering and scientific domains. It explores operator-based learning methods proposed for solving complex problems in spatiotemporal data and PDEs as an alternative to the complexity of numerical solvers. The proposed method introduces an efficient time-shift operator for capturing multiple levels of granularities and irregular sampling without the need to approximate the kernel while maintaining linear complexity.", "essential_references_not_discussed": "Based on the review comment above (**\"Methods And Evaluation Criteria\"**), several baseline methods for irregularly-sampled temporal forecasting are not discussed, including the recent SOTAs based on Neural SDEs (several methods are mentioned in [1]). For regular temporal forecasting, a more standard benchmark followed in the community is based on Time Series Library [2,3] (methods such as TimesNet, TimeXer, iTransformer, and PatchTST are first in ranking yet not tested in the paper).\n\n[1] Oh, Y., Lim, D. Y., & Kim, S. (2024). Stable neural stochastic differential equations in analyzing irregular time series data. arXiv preprint arXiv:2402.14989.\n\n[2] Time Series Library https://github.com/thuml/Time-Series-Library\n\n[3] Wang, Y., Wu, H., Dong, J., Liu, Y., Long, M., & Wang, J. (2024). Deep time series models: A comprehensive survey and benchmark. arXiv preprint arXiv:2407.13278.", "other_strengths_and_weaknesses": "Summarized **strengths** of the paper:\n- **[S1]** The authors propose a thoroughly presented operation-based theoretical framework that aims to tackle real-world spatiotemporal problems beyond numerical solvers. \n- **[S2]** Broad tasks, such as temporal forecasting and irregular sampling, are approached.\n- **[S3]** The included visualizations enhance the readability of the work.\n\nSummarized **weaknesses** of the proposed study:\n- **[W1]** The work is poorly presented compared with the neural-operator-based related works in the literature (such as FNO). Such related methods are abstractly mentioned in the introduction, and several details are given in the methods section, but very few are in the related work section. I suggest restructuring the related method's details throughout the text to highlight the proposed method’s contributions. \n- **[W2]** The related work section is misplaced and limited regarding references. Only a few methods (Koopman, FNO) are explained in terms of operator-based learning, and one method is explained for irregularly sampled time series modeling (NeuralODEs). Baseline methods used in experiments should be explained in more detail along with more recent methods (see relevant review section above).\n- **[W3]** Several mixed benchmarks in terms of datasets but with important state-of-the-art methods missing, especially for temporal forecasting and irregularly-sampled temporal modeling (see relevant review section above).\n- **[W4]** The proposed method is, in several cases, outperformed by baselines, particularly for the temporal forecasting datasets. However, the performance significance of results is not justified if studied, e.g., were multiple runs with random seeds followed for all datasets and baselines? In some cases, the results for the baseline are directly taken from the relevant papers.\n- **[W5]** The experiments on irregularly sampled time series are limited to simple synthetically generated data (2D spiral), which raises questions about the method's efficacy in real-world setups with irregular timestamps. \n- **[W6]** Studies on the computational complexity of the proposed KRNO methods compared to SOTAs for temporal forecasting (beyond operator-based methods) are missing. This is essential to support the applicability of the process in practical scenarios with vast time series datasets, where simple and lightweight architectures are, in several cases, more favorable.", "other_comments_or_suggestions": "No comments", "questions_for_authors": "1. **[Q1]:** Could authors improve the presentation and structuring of the related work (based on **[W1], [W2]**)?\n2. **[Q2]:** Could you tackle the issues raised on experimental evaluation in terms of temporal SOTAs, methods applied to irregularly sampled data, and real-world irregularly sampled datasets (based on **[W3], [W4], [W5]**)?\n3.  **[Q3]:** To what extent are neural operators for time series forecasting constrained to data exhibiting clear structures, smooth evolution, and spatial correlations? (for instance, NeuralODEs solve irregular sampling problems for specific physical datasets/are not successfully extended to standard forecasting)", "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:

Meta-review: The paper proposes a continuous in-time learning for operators called Khatri-Rao Neural Operator (KRNO). The operator is modeled as an integral transform with a non-stationary kernel, which has also been studied in the past but, this time, decomposed via the Khatri-Rao product structure. The resulting model enables forecasting with irregularly sampled data and super-resolution.

The authors are encouraged to clarify in their papers that there are earlier works that considered mapping from time interval to time interval in function spaces, making the underlying kernels to be time dependent. The early studies of these sorts are UNO and "Tipping Point Forecasting in Non-Stationary Dynamics on Function Spaces".
Decision: Accept (poster)