Time Series Kernels based on Nonlinear Vector AutoRegressive Delay Embeddings
Authors: Giovanni De Felice
John Y Goulermas
Vladimir Gusev
Venue: NeurIPS 2023 poster (2023)
Keywords: Time Series
Kernel methods
NVAR processes
Dynamical systems
Reservoir Computing
TL;DR: We introduce a new kernel for univariate and multivariate time series that compares the linear dynamics extracted from NVAR delay embeddings and overcomes downsides of Reservoir Computing based approaches.

Abstract: Kernel design is a pivotal but challenging aspect of time series analysis, especially in the context of small datasets. In recent years, Reservoir Computing (RC) has emerged as a powerful tool to compare time series based on the underlying dynamics of the generating process rather than the observed data. However, the performance of RC highly depends on the hyperparameter setting, which is hard to interpret and costly to optimize because of the recurrent nature of RC. Here, we present a new kernel for time series based on the recently established equivalence between reservoir dynamics and Nonlinear Vector AutoRegressive (NVAR) processes. The kernel is non-recurrent and depends on a small set of meaningful hyperparameters, for which we suggest an effective heuristic. We demonstrate excellent performance on a wide range of real-world classification tasks, both in terms of accuracy and speed. This further advances the understanding of RC representation learning models and extends the typical use of the NVAR framework to kernel design and representation of real-world time series data.

Time Series Kernels based on Nonlinear Vector AutoRegressive Delay Embeddings
Giovanni De Felice
John Goulermas' Vladimir Gusev Department of Computer Science
University of Liverpool
gdefe@liverpool.er ac.uk gusev@liverpool.co ac. ac.uk
Abstract
Kernel design is pivotal but challenging aspect of time series analysis, espe- cially in the context of small datasets. In recent years, Reservoir Computing (RC) has emerged
powerful tool compare time series based on the underlying dynamics of the generating process rather than the observed data. However, the performance of RC highly depends the hyperparameter setting, which is hard to interpret and costly to optimize because of the recurrent nature of RC. Here, we present new kernel for time series based on the recently established equiva- Ience between reservoir dynamics and Nonlinear Vector AutoRegressive (NVAR) processes. The kernel is non-recurrent and depends on small set of meaningful hyperparameters, for which we suggest an effective heuristic. We demonstrate excellent performance on a wide range of weal-world classification tasks, both in terms of accuracy and speed. This further advances the understanding of RC representation learning models and extends the typical use the NVAR framework to kernel design and representation of real-world time series data.
Introduction
Time series are arguably one of the most important types of data in the modern (Hamilton 2020) and ubiquitous both in scientific research (Strogatz 2018) and practical applications (Zhang et al. 2018 Zeroual et al. 2020). A key element in the design of most machine learning protocols for time series is a quantification of similarity (Ding et al. 2008 Abanda et al. 2019 Echihabi
al. 2020). This especially true for kernel methods (Schölkopf et al. (001), which search for linear solutions after projecting the data into a higher (possibly infinite) dimensional space and represent effective alternative to non-linear models. The performance of kernel methods is heavily impacted by the definition of a positive semi-definite (PSD) similarity function for the type of data hand, i.c., kernel. However, their design is challenging for structured data such as univariate (UTS) and multivariate time series (MTS), which exhibit temporal autocorrelation, inter-attributes (or dimensions) dependencies and possibly variety of temporal distortions and misalignments [Paparrizos et al. 2020).
To overcome this, promising approach is identify and compare the underlying dynamics rather than the raw observed data. Notably, kernels based on Reservoir Computing (RC) stand out this context (Chen al. 2013 Bianchi al. (020). RC-based kernels layer of recurrently connected nodes (the reservoir) to map each time series into a rich dynamical feature space. The similarity is then obtained by comparing individual readouts capturing dynamics. In practice, RC is highly sensitive to a large set of hyperparameters having little interpretability, the optimization of which is challenging for the comprehensive range of values.
randomized and untrained
Deceased May 2022
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Recently, theoretical work by Bollt 2021) has demonstrated that a simple RC can be formally rewritten as non-linear vector autoregressive (NVAR) model. In the NVAR framework, the reservoir is replaced by simple concatenation of the input series with time-delayed copies and nonlinear functionals, such as products (Gauthier et al., 2021). This reduces the number of hyperparameters and has proven to be very effective in chaotic systems forecasting (Shahi et al. (022), as much to earn the title of "Next-Generation Reservoir Computing". However, its applicability and performance for real-world data and outside of forecasting dynamical systems are largely unexplored.
In this work, we investigate whether an NVAR process can replace reservoir for kernel design and if the extracted features are just as effective in comparing real-world time series data. We propose a very efficient and effective NVAR kernel for UTS and MTS data. The main idea is follows (Fig. [] First, we enrich each time series with lagged copies of itself and additional nonlinear terms, forming a high-dimensional NVAR embedding. Then, we extract a linear parameterization of the time evolution in the embedding space, representative of the underlying dynamics, analogously to RC representation learning methods (Bianchi et al. 2020). Finally, these fitted parameters are used to compute the similarity between the time series. Our main contributions branch out into different communities investigating time series, namely, kernel design, RC and NVAR, bringing specific benefits to each of them:
From a kernel design perspective, we provide an NVAR kernel for UTS and MTS and a general parameter setting based on simple heuristics. Based on our experiments range of datasets, the proposed approach matches the state-of of-the-art art (SOTA) in terms of accuracy and is considerably faster than SOTA, representing the best compromise between accuracy and efficiency.
wide
From the perspective of RC, we advance its success for TS representation by incorporating NVAR-made dynamics. The resulting model is more accurate, interpretable, and non- recurrent, which overcomes difficulties in hyperparameter optimization and sheds light on the potential RC performance.
From an NVAR perspective, we extend the framework to real-world time series represen- tation and kernels, well beyond the original work on forecasting synthetic, noise-free, and chaotic dynamical systems (Gauthier et al. 2021).
We establish a connection to Takens' theorem (Takens) 1981) and the field of state space reconstruction (Sauer al., 1991), which theoretically underpin our approach, highlighting a compelling avenue for the use of dynamical systems theory in machine learning.
1.1
Notation
Through the paper, indicate variables as lowercase (z); constants uppercase (X); vectors and UTS as bold lowercase (x); matrices and MTS as bold uppercase X. An index between square parenthesis [n]] indicates the n-th sample of set. For an MTS X, the corresponding lower-case X indicates all dimensions at time stamp I while x" indicates dimension d at all timestamps.
2 Related methods
In this section, we provide a background on the previously proposed kernels for time series data, grouped by category. Despite the abundance of time series similarity measures (ang & Shahabi 2004 Paparrizos & Gravano 2015 Janati et al. 2020), we limit the discussion to PSD metrics.
Lock-step Lock-step approaches view time series static vectors and directly use common distance metrics (Chal 2007). The PSD property is usually obtained by superposing a linear or radial kernel (Schölkopf et al. 2002). This is very efficient, but ignores any temporal structure within the scries and is not applicable for different lengths.
Elastic Elastic measures account for time distortions (e.g. shifts) or different lengths between sequences. Lu et al. (2008) proposed to interpolate the series and treat the problem as a distance between curves. The Global Alignment Kernel (GAK) (Cuturi 201 and KDTW (Marteau & Gibet 2014) instead are built on the well-known Dynamic Time Warping (DTW) (Berndt & Clifford 1994 and compute the similarity from the cost of one-to-mang alignments between the series. The

---

Shift-invarinn kernel (SINK) Paparrizos & Franklin 2019) computes cross- corseloorol similarity in the Fourier domain. However, for the aforementioned methods, the computation of the pairwise similarity is computationally expensive and does not account for the relations between different attributes in the multivariate case.
Model-based Model-based kernels process the series with probabilistic or deterministic model and base the similarity on the extracted information. Such a model can be a single generative model, as in the Fisher kernel (Jaakkola al. 1999) the probability product kernels (Jebara al. 2004), a parameterized family of probability distributions, as in the autoregressive kernel (Cuturi & Doucet 2011). Similarly to the latter, the Time Cluster Kernel (TCK) (Mikalsen et al., 2018) is obtained from an ensemble of Gaussian mixture posteriors, sharing the same parametric form but trained on different subsets of the dataset. Among non-proboa non-probabilistie approaches, the Learned Pattern Similarity (LPS) (Baydogan & Runger 2016) uses an ensemble of regression trees to extract local patterns from each series. The downside of all such methods comes from the specific functional form, which may limit the generalizability of the extracted features. Ensemble methods are also unsuited for datasets with few training samples.
Reservoir-based Reservoir Computing (RC) is a class of recurrent neural networks that keep recurrent connections untrained in order to overcome the high cost of pack-propagation through time (Werbos 1990) and the vulnerability to exploding and vanishing gradient (Pascanu et al. (013)). Its simplest form, the Echo State Network (ESN) (Jaeger 200T), is composed of three layers; an input layer, a hidden layer of connected neurons (reservoir) and readout layer:
input:
u1 W
reservoir:
r. -1- (1-a)ra +af(Arr- +u, +x)
(2)
readout:
Yt Wo tc
where x, IRd, u, gd. and r, E Ed, are, respectively, the input, its upscaling projection and the reservoir states at time t (with d, >> d, and ru 0); 0 S < 1 is the leak parameter; ep an additive noise and f is any non-linear activation function. The reservoir (A E xd.) and the projection (Win E Xdx) matrices are kept untrained, while the readout (Woat E xd.) and the bias term (c) are learned by fitting the reservoir states to the output (y+ E Rdv) via ridge regression.
Despite their randomness, the reservoir states constitute rich pool of heterogeneous features that give complete knowledge of the underlying dynamical system generating the observed time series (Lokse) al. 2017 Bianchi et al., 2020). In fact, as opposed to fitting a specific functional form, solid theoretical ground ensures the quality of the reservoir features (Hart et al. 2020: Gonon et al. 2023 Gonon & Ortega 2021 Hart et al. 202T). Nevertheless, performance is highly sensitive the choice of a large set of hyperparameters controlling the reservoir initialization and update (Lukoševičius)
2012). For instance, input scaling ([' |) and spectral radius (p(A)) assure stable reservoir dynamics (Gallicchio) 2019) and control the degree of injected nonlinearity. This is not easy to judge and requires experienced insight into nonlinear dynamics. In general, the interpretation and setting of ESN parameters is complex task and, to this day, an active area of research (Dong et al. 2022 Steiner et al., 2022 Zhang et al., 2022). Simple supervised optimization, such as cross-var cross-validation (CV), is often prohibitive, the optimization space is large and the recursive nature of the reservoir forces the computation of Eq. [ for the whole length of the series before assessing the performance
of one hyperparameter setting.
Regarding RC-based kernels, Chatzis Demiris (201 originally built radial kernel directly on top of the reservoir states. In Cneh et al. (2013), all series are processed by shared reservoir and compared based on the respective readouts, which are individually trained on the one-step-ahead prediction task (y: x)+1). Finally, Bianchi et al. (2020) observed that more information is retained by using y. = rt+1, and an effective kernel can be obtained by applying a radial function on the vectorization of all the readout weights
3

---

T++1
Wout I
|vec (
X
Rigg
R
Raaninn
Linear readout
Input time series: X
NVAR embedding: RNVAR d
Figure Building of each representation vector the NVAR kernel. The input series is concatenated with Ring. which contains & lags of all input dimensions, and with Raoolin. containing products between dimensions random subset of dimensions is considered. The readout then extracts the linear dynamics (Wout) of the embedding states and the bias term (c) by performing a ridge regression fit along the temporal dimension. These fit parameters vectorized to form the representation.
X and Ring. Among all possible concatenation terms, only a
3
Methodology
This section presents the main building blocks of NVAR models, followed by our proposal of an NVAR-based kernel and general setting for the hyperparameters that govern it.
3.1 The Nonlinear Vector AutoRegressive model
Recently, the work of Bollt (2021) has established a connection between the ESN architecture and Nonlinear Vector AutoRegressive (NVAR) models. Under specific conditions, this takes the form of formal equivalence, where the coefficients of an NVAR process can be expressed in terms A. W and Worr By building
such equivalence, in Gauthier al. (2021), ) the input layer (Eq. I and the reservoir (Eq.22 2 of an ordinary ESN are replaced with the deterministic operations that lead to a feature matrix RVVAR- In detail, this is the result of concatenating three terms: the input series (X), lagged copies of the input, each of which is increasingly spaced by s time-stamps (Rigg), ) and all possible products between lagged and unlagged dimensions, up to polynomial order n (Rmontin):
Flag,t X Xz-2a 2s
X
Inontin. = xx Flggtt
[ [x [lag.]]
[8]
(n times)
(4)
R = X Riag Rontin
where -|. is the column-wise concatenation and all unique monomials.
[@] 1 performs the outer product and concatenates
This procedure is very efficient as it is non-recurrent and the size of R AR is usually smaller than typical reservoir. On top of that, it is determined by small set of integer hyperparameters: k. s and n.. In principle, to sustain the equivalence to RC, an infinite amount of lags (k: - xo) should be considered. Nevertheless, it has been observed that small with exceptionally well in typical dynamical systems forecasting tasks (Gauthier et al. [202]] Shahi et al.) 2022; Gauthier al. (022).. However, applicability and performance of NVAR models beyond forecasting yel remain largely unexplored, particularly in the context of representing and comparing real-world time series data. this regard, they present promising solution to significantly enhance RC representations and kernels.
still perform
3.2
NVAR kernel
We present here the NVAR kernel (NVARk), which integrates the NVAR framework with RC-based kernel architectures. Given pair of time series, X Xij], with dimensionality d, of possibly different lengths, NVARK operates in three steps, which we present in the following paragraphs. Fig. I depicts a graphical reference, while Algorithm [ explicits the algorithmic process.
2 other forms of nonlinearity can be used depending on prior knowledge, but the polynomial functionals can be motivated universal approximators of dynamical systems (Franz & Schölkopf 2006)

---

Algorithm NVAR kernel
input MTS list X|1 N), number of lags k. lag size s, embedding size d,
I: L ((a,p) l E {1...., dg p {1, ..... k}} (dimension a lagged p times)
2: P P {((a,p), ) (b. 9)) a E (1.... dz p.q = {0, ..., k}} (all possible products)
3: samples <- sample (d, dz) elements from L U P at random
4: for in [1: .V] do
5: [ < X[i]
6: for elem in samples do
7:
if elem (a, is in L then
8:
concatenate lag term: R[i] <- [R[]]
|| xi xi-pxxli]]]
9:
end if
10:
if elem ((a. p). (b. 9)) is in P then
11:
concatenate nonlinear term: R R
R
x² xi-pxal))
[x]
x [x]
12: end if
13: end for
14: fit a linear model mapping r.- [ to r.
15: vectorize fit parameters to obtain representation 0|i]
16: end for
17: for all representation pairs: Ki +- RBF(O 0[j))
output K
NVAR embeddings The general idea is to transform the time series using the NVAR map defined by Eq. 4 with n 2. However, for a given choice of the number of lags k; the resulting number of additional dimensions scales with k2 For moderately high da, this can result in high collinearity and raise the curse of dimensionality in the following readout, making it slow and inaccurate. It also limits the possible choices of k and, with it, the ability to capture longer-term memory effects. The method, as it is, then only seems tractable in the context of dynamical systems, where de is usually very contained. To address this issue and make the method applicable to higher dimensional datasets, we propose to concatenate just a random subsample of all possible dimensions:
RNVARla X
rd
dd randomly sampled columns in Riag ] with uniform probability and d, a with (ralam Ia= threshold hyperparameter (notation [t] is omitted for all terms and the same operations are applied to time series (j)). This converts the deterministic NVAR map into a random feature map. We support this choice in the following, by pointing out an unexplored connection of NVAR models to the theory of state space reconstruction (SSR) (Kantz & Schreiber) 2004 Strogatz 2018).
It is generally true that we can interpret the observed data as a realization of more complex underlying dynamics (Ø*), acting on some unknown states (x*): xi+1 6" (xx)). SSR focuses on using the observed data (x) to construct embedded states (s,) for which the dynamics is topologically equivalent to ø², and thus a better representation of the underlying system (Casdagli al. 1991). In this regard, the celebrated Takens' theorem (Takens) 1981) claims that a valid embedding, called delay embedding (Packard et al.) 1980), can be formed by concatenating the input with its lags and, most importantly, guarantees that not all lags are needed, but only sufficient number &: 8; = [xt, st-s; 2... Es. More recently, the work of Deyle generalization of Takens' theorem to multivariate inputs in which, to match E, lags can be spread over different observed dimensions and alternative functions of the input can also be considered. All terms in the NVAR representation, i.c., lags and nonlinear terms, are reminiscent of this formulation, hence our choice of not considering all terms, but posing instead an upper threshold on the number of concatenation terms (Eq. 5).
Sugihara (2011) presented
This also allows for different perspective on our approach. We can interpret R AR la_ as an attempt to create states for which the dynamics approximates the underlying one generating the data. As there is generally a loss of information in the generation process, a comparison based on such states would be more truthful than directly using the observed data. In practice, the presence of noise separates us from obtaining such exact states (Casdagli et al. 199T) and places our method under the field of embedology (Sauer et al. 1991). i.c., the building of delay-observvatin maps with special
5

---

features, overcoming the sensitivity to noise in the original Takens' formulation. This theoretical framework is more extensively discussed in Appendix
An alternative to Eq. 5 would have been to generate the full NVAR representation followed by a dimensionality reduction module in Bianchi al. (2020). In our preliminary evaluation, this proved to be computationally more expensive and led to poorer performance, possibly because of the dimensionality reduction not preserving the delay embedding structure.
Linear readout Individual linear readouts (Eq. 3 are trained by ridge regression to solve a one- step-ahead prediction task on the same embedding states, i..., using Yt r++11 The parameters of each model are then extracted and vectorized to obtain the representations O and Ob]] This step removes the dependence on the length of the TS, which means that the kernel compare TS of different lengths.
As the readout operates on the input concatenated with its lags and quadratic products, these rep- resentation vectors can be interpreted as encapsulating statistics related to the auto-mutual and mutual information within the dimensions of the input. In contrast, representations from RC methods encompass hardly interpretable reservoir dynamics.
RBF aggregation
Finally, radial function is used to ensure the PSD property and form the kernel:
K (X xbI = exp {
0bi ||9 /2
3.3 Hyperparameter setting
We display here the hyperparameters to construct NVARK and propose heuristic expressions.
Number of lags (k) and lag size (s) The parameters and play pivotal roles in the NVAR frame- work (Eq.4) 4 and, consequently, in approach. These key parameters of clear interpretation: k controls the number of lags and extends memory capacity, i.e, capturing significant dependencies that extend further into the past; controls the spacing between lags, which manages the temporal resolution at which these dependencies are probed. In contrast, key RC parameters are arguably less interpretable.
In Takens' formulation, in the case of infinitely long, noise-free time series, there is no importance bias towards any specific lag s. and their number & is the only relevant factor. This is often determined using the False Nearest Neighbours algorithm (Kennel et al. 1992 Wallot & Mønster 2018). However, when dealing with limited noisy observations, the choice of becomes equally critical as it directly relates to the sampling frequency (Broomhead & King 1986). In particular, using too high a value would force modeling of uncorrelated points while too low would model noise. correct choice of s instead allows for extracting seasonality or other relevant trends as well as the underlying dynamics. In practice, for noisy systems, the choice of both k and is inherently difficult problem due to the absence of clear theoretical guidelines (Tran & Hasegawa, 2019). Despite
be exploited to propose task-dependent heuristies (Small & Tse 2004). often requiring a trade-off between redundancy and the irrelevance of concatenated features.
that, interpretability
In this work, propose two settings. The first is rule of thumb and proceeds as follows. We first apply a trendfiltering method to the time series, such as an & trend (Kim et (0099) (this allows for choice of that captures the temporal variability of the system rather than the superposed noise). We then set the product k $ as the average distance between peaks and hollows (dbp) (Kugiumtzis) 1996). Finally, we balance the contributions by adopting = k /dbp. Alternatively, in supervised setting, both parameters are optimized by CV.
Polynomial order (n) The polynomial order controls the amount of non- non-lincarity in the embedding. As in most previous works, we found no need for further complexity than n 2.
Embedding dimension (dz) In Takens' formulation, an embedding size double the dimensionality of the underlying manifold is sufficient for reconstructing the system dynamics. However, when dealing with real-world data, this is unknown, and a trade-off must be struck between the omission of relevant features and redundancy, which can increase collinearity and computation times. In line with the scope of our work, we choose to give priority to fair comparison with previous RC

---

methods (Bianchi et al. 2020), which identify d, 75 as an optimal dimensionality before the readout. In Appendix C.I we undertake brief exploration determine whether this value is also a reasonable choice for NVAR embeddings. While we did not observe any redundancy in univariate datasets, we did observe a performance drop after different di for some multivariate datasets. The study suggests 70 E d, x 100 as suitable range.
Linear readout regularization (Aridge) For the readout, we adopt the optimal regularization (OCReP) proposed in Cancelliere et al. (2015), i.c., set Aridge to the product between the minimum and maximum singular value of RNVARId,
RBF lengthscale (Trer) We set the lengthscale in Eq. to to the median pairwise distance between all the representations 0|i] in the training set.
3.4 Asymptotic complexity
NVARk is based on the individual representations of time series. Hence, the computationally intensive step performed only once per time series and results in O(N) complexity. This also implies that the computation can be parallelized along N. This provides significant computational advantage over comparison- kernels, which must perform expensive computations for all the relevant pairs of series, i... O(N2).
In terms of scalability with the length of the series (T), the advantage of NVAR over any recursive RC approach is the critical difference in how the hidden states are created. Eq. i not iterative and only performs a fixed number d, dz of concatenation operations, independent of the length of the time series. In contrast, reservoir-based based methods are limited by the expensive recursive update of the reservoir state (Eq. 2. which is O(T))). where d, is the size of the reservoir (which is usually quite large). As evidence, improving the scalability of the reservoir update is an active area of research (Dong et al. 2020). The complexity of NVARK is limited to the only ridge regression of the output layer. This corresponds to solving linear system with T training examples and D) features, where D is at most dz. Such a system can usually be solved efficiently, e.g. by using LU (or Cholesky) factorization in asymptotically (T))).
As for the scalability with the dimensionality of the input series, all operations are bounded by dz which we treat as constant and do not discuss complexity in connection to this parameter. Overall, NVARK exhibits an asymptotic complexity at most (NTT)).
Experimental evaluation
In this section, we present an experimental demonstration of the performance of NVARK. In line with the established literature (Cuturi 2011 Cuturi & Doucet 2011 Baydogan & Runger 2016 Paparrizos et al. 2020). our main evaluation consists of time series classification. In particular, we pair different compued kernel matrices with Support Vector Machine (SVM) classifier (Steinwart & Christmann 2008), a popular kernel classification method that lies between simple linear models and more advanced deep learning architectures. The evaluation is performed over 130 UTS datasets from the UCR archive and 19 MTS datasets from the UEA archive (Bagnall A. & E.) (Appendix II.)) and focuses on accuracy, execution time and scalability. In Appendix C.5 A
present snapshot of kernel Principal Component Analysis (kPCA) for dimensionality reduction and visualization.
Performance of the proposed approach is investigated under two different hyperparameter settings. We indicate with NVARK the setting in which and s are determined by the general proposed rule of thumb of Sec. 3.3 The maximum number of dimensions is set to d, = 75. As for the other hyperparameters, we follow the setting given in Sec. 3.3 We instead indicate with NVARK* the configuration in which & and s optimized in CV loop over a small grid. As the datasets cover wide range of series lengths, the CV grid is slightly adapted consequently (Appendix. B.2).
The comparison is performed against previously proposed kernels: the reservoir-model Echo State Network (mmESN) (Bianchi et al., 2020) as the reservoir equivalent to our approach. We set D = 75 as the dimensionality before the linear readout. This allows us to directly compare the quality of the extracted dynamic features with or without reservoir; the Time Cluster Kernel (TCK) (Mikalsen et al. 2018) as representative of model-based kernels and SOTA for MTS; the Shirt-LNvariant
7

---

Kernel (SINK) (Paparrizos & Franklin, 2019) as representative of Fourier-based methods and SOTA for UTS. A comparison with the Global Aligmment Kernel (GAK) (Cuturi, 201T) is presented in Appendix C.4 Guidelines in Appendix B contain the hyperparameter settings and links to the code implementations for NVARK and all the aforementioned baselines.
Despite the fact that some datasets in our benchmark can be also investigated using Deep Learning (DL) methods ([smail Fawaz al. 2019), the use case for kernel methods is substantially different, usually focused on small datasets, where the application of DL is not trivial. In particular, kernels can compute the similarity between just two time series right away. Consequently, in the literature, kernels are usually not directly tested against DL. Nevertheless, in our work, we compare to ESN- based representations, which are known to be competitive with variety of DL architectures (Bianchi et al.) 2020 Shahi et al. (022). Thus, ESN can be seen as proxy for more complex architectures.
As preprocessing, we apply rero-padding to match all series to the maximum length in each dataset (to allow comparison with TCK) and follow the common practice of standardizing all series to zero mean and unit variance. For each dataset-method pair, we allow for maximum execution time of 20 hours.
4.1 SVM classification of UTS
We show here SVM classification for univariate datasets. SVM requires the setting of the hyperplane regularization C, which we optimize in a CV loop from logarithmically spaced values in [10-3, ,10%). This is the only free parameter in all methods except NVARK* in which is jointly optimized with and s. For the CV loop, we adopt 10-fold CV with the size of the validation set corresponding to 33% of the train set. SVM accuracy for each dataset obtained by averaging different seeds. Tab. reports the average classification accuracy for all kernels within the presented framework, well as the relative ranking (individual results are reported Appendix Overall,
jointlly
repetitions with
Table 1: SVM classification of UTS datasets, performance metrics across 130 datasets.
NVARK*
NVARK
RMESN
TCK
SINK
AVERAGE SCORE
0.803
0.779
0.757
0.726
0.771
AVERAGE RANK
/
2.14
2.62
2.85
2.20
N FIRST RANKED
/
52
26
23
43
MEDIAN TIME (S) / DATASET
14.27
0.59
40.87
44.24
23.73
NVARK obtains very competitive results, with the highest average accuracy and rank. The case of univariate datasets is comfortable scenario for NVARK to operate in as, in most cases, all possible terms in RVVAR can be considered without exceeding d, and the effect of collincarity is minimized. NVARK and SINK appear indistinguishable under Friedman Nemenyi statistical test, which places NVARK under the SOTA in terms of accuracy. Interestingly, the difference between NVARK and mEESN is significant at the 95% level with p-valuewvARI--mm mMFS 0.013. Results for NVARk* show that the performance of NVARK can be very efficiently improved by a CV optimization of just the two integer embedding hyperparameters.
4.2 Execution time and scalability
In terms of execution time, the scalability properties of NVARK (Sec. 3.4 manifest in exceptional performance. In Fig. (1eft) and Tab. [ (last row), we compare running times to compute the train- train and the train-test kernel matrices and observe that the improvement of NVARK is consistent across datasets and baselines, with speedup factors up to x 100 in the high-T and high-N regimes.
Results for NVARK* can be achieved with a slowdown factor of ~ 25 (parallelized grid search + 1 iteration with optimal parameters) with respect the heuristic- based NVARK. Interestingly, for 113/130 datasets, NVARk* execution time is faster than a single iteration for mEESN.
For study of the scaling, in Fig. 2 (middle, right) we progressively vary the length of the Abnormal- Heartheat dataset (by interpolation) and the number of training samples (by random sampling) the Crop dataset. Results show that NVARK performs exceptionally well in all regimes, while rmESN
8

---

excedded
FIEES TOX SINK
wimiSNN TTCK+
SINK
1000
Abnornal Camooroo1Bearbboo Beartbeat
1000
Exop
180
Time (5)
TTmme
20
40
50
80
100
120
10
100
1000
10k
10
300
1000
10k
UCR Cotaset index
Series length
Training samples
Figure 2: Left: For all univariate datasets. the plot shows the ratio between the execution time to compute the train-train and train-test kernels of different methods and NVARK. Middle: Scalability of kernels with the series length. Right: Scalability of kernels with the training size.
scales poorly with the length of the series and TCK and SINK with the size of the training sample. NVARk, instead, is advantaged by its non-recursive structure and linear complexity in N.
4.3 SVM classification of MTS
We present here SVM classification over the 19 MTS datasets. We not include SINK as there is no extension evaluation provided for MTS. Aside from that, the evaluation is conducted identically to the univariate case. For multivariate datasets, we expect to drift away from the ideal operating
Table 2: SVM classification of MTS datasets. Average across 10 seeds. Best accuracy is in bold.
DATASET
NVARK*
NVARK
RMESN
TCK
SPOKENARABICDUS
0,980
0.980=a.0o
0.948
0.985
JAPANESEVOWELS
0.981 1.0.05
0.975
0.968 004
0.961 10.004
PENDIGITS
0.986
0.983 00000
0.976 001
0.954 10.001
RACKETSPORTS
0..09
0.847 0.000
0.859 4:0 008
0.872 0.0006
LSST
0.566
0.543320.00
0.597 1.59772.mm 20.000
0.474400.00
LIBRAS
0.972
0.900
0.8244..0m
0. 0.793
FINGENMOVEMENTS
0.572
0.581 03.032
0.576 - .22
0.505
NATOPS
0.896 0..06
0.895 +0.019
0.828 +0.001
0.806
CHARACTERTRAL.
0.982++0.00
0.915=uu.000
0.987 =00.002 002
0.974++0.002
ERING
0. 0.898
0.864
0.814 04.00
0.957
BASICMOTIONS
1.000
0.938 027
0.995 011
1.000
ARTICULARY WORDREC.
0.978 006
0.976 00.
0.986 1786.1000 008
00.03
EPILEPSY
0.086
1986600.000
0.97 1:0 005
0.947 0.947 10.005
UWAVEGESTURELIB.
0.946
0,939 #0.011
0.862
0. 0.868+0.007
SELFREGUULATIONSCP)
0.7 0.746+0.017
0.688-0.032
0.769
0.827+4.007
SELFREGULATIONSCPP
0.566
0.567 20.028
0.51 1 .207
0.488
CRICKET
0..00
0.986 = 000
0.982 006
0.958 00.000
STANDWALKJUMP
0..03
0.607 0.009
0.373 o58
1..03
EIGENWORMS
0.944
0.950
0.939
0.594
AVERAGE SCORE
0.858
0.848
0.830
0.809
AVERAGE RANK
1.68
2.05
2.26
N FIRST RANKED
11
3
5
MEDIAN TIME (S) DATASET 19.67
3.87
26.47
184.08
scenario of NVARK. In particular, all the possible dynamic features may be equally effective and, by randomly selecting few, important correlations may be missed or spurious ones captured. Despite that, the experiments Tab. Tab.2 surprisingly good results, even for datasets for which d, is just fraction of all possible concatenation terms. We believe that this finding strengthens the connection to the generalized Takens' theorem. It is also interesting to notice good performance in diverse cases of time series that do not immediately arise from dynamical systems. Exploring this aspect constitutes interesting future works, e.g., time series representation learning and state space reconstruction as possible underpinning.

---

Finally, we observe that the two hyperparameter configurations, NVARK* and NVARK, tend to be ranked next to each other. This supports the quality of the proposed general setting, which allows the kernel to extract meaningful features while preserving great computational efficiency.
4.4 Ablation study
In order to tease out which NVARK component most important, here present an ablation study, We separately consider fitting the linear readout directly on the input time series (no concatenation), sampling only linear terms from Eq. [S] and sampling only non-linear terms. This is an interesting ablation study that, to the best of our knowledge, has not been investigated in previous works. Results are presented in Tab.
Table 3: Ablation study across 130 UCR univariate and 19 multivariate UEA datasets
no concat
linear only
non-linear only
all
balanced
SVM score for univariate datasets
mean std
0.536
0.707
0.764
0.779
0.779
0.010
0.008
0.006
0.006
0.006
SVM score for multivariate datasets
mcan
0.765
0.840
0.840
0.848
0.849
std
0.012
0.008
0.013
0.014
0.012
As general trend, observe that using all terms leads to the best results, followed by non-linear only, linear only, and fitting directly on the input series. Note also that results are in line with our interpretation using the generalized Taken's 's theorem. For univariate datasets, linear only tends to underperform, as the total number of concatenated dimensions is very small (equal to k) and may be not sufficient in reconstructing the underlying manifold. In the multivariate case, the two variants perform very similarly, as the lags of all attributes may already constitute a considerable pool of terms.
To conclude, we consider an additional variant of Eq. 5 in which we first concatenate all the lag terms, and then fill the remaining dimensions, up to dy,, by sampling nonlinear terms. Interestingly. for multivariate datasets, this leads to slightly better performance and a reduction in variance. We believe this might eventually compensate for the imbalance between all possible linear and non-linear terms, while also reducing the chance of adding noise from uncorrelated dimensions.
5 Conclusions
We have proposed a kernel for time series that integrates the NVAR framework into reservoir- based kernel architectures. NVARK compares time series based on the linear dynamics of NVAR embeddings, which are built from concatenating lags and nonlinear functionals the original series. In terms of accuracy, NVARK outperforms the corresponding RC architecture. Computationally, it is exceptionally efficient and based on a few integer hyperparameters, which together allow for further improvement of the results with simple supervised grid-based optimization.
As for future work directions, we believe it would be effective target specific weaknesses of NVARK. In particular, as the input dimensionality of the time series approaches the maximum embedding dimension (d the proposed approach finds little margin for concatenating more lags and non-linear terms. Such cases would require substantially increasing d, and, consequently, place excessive strain on the linear readout, which inherently possesses limited expressive capacity. We then expect the effectiveness of NVARK to decrease. To overcome this, one can consider replacing the random sampling of dimensions with more refined strategies that prioritize the selection of meaningful terms. Alternatively, it would be interesting to explore how different linear layers, e.g. Lasso, would perform in this regime. Alongside these, different unsupervised strategies can be considered learn optimal values for the embedding parameters. In particular, we note an unexplored affinity with signature transforms (Chevyrev & Kormilitzin, (016), which we plan on deepening to infer optimal settings for the maximum dimensionality of the NVAR embedding.
10

---

Acknowledgements
G.D. acknowledges the Beckers Group for funding this research. V.V.G. thanks Leverhulme Trust for support via the Leverhulme Research Centre for Functional Materials Design. We would like to thank Chao Huang, Filippo Maria Bianchi, Davide Bacciu, and Claudio Gallicchio for their fruitful feedback and support for the work. Also, we would like to thank all the reviewers for their invaluable comments,
References
Abanda, A., Mori, U. and Lozano. J. A. A review on distance based time series classification. Data Mining and Knowledge Discovery, 33(2):378 -1238788812211 2019.
Bagnall A., Lines J., V. W. and E., K. The UEA and UCR time series classification repository.
WWW
Eimsserreesclasssficaaiii com
Baydogan, M. G. and Runger, G. Time series representation and similarity based on local autopatterns. Data Mining and Knowledge Discovery, (0(2):47-5500 2016.
Berndt and Clifford. Using dynamic time warping to find patterns time series. KDD workshop. 10 (16):359-370, 1994.
Bianchi, F. M., Scardapane, S., Løkse, S., and Jenssen, Reservoir computing approaches for representation and classification multivariate time series. IEEE transactions on neural networks and learning systems, 32 19/2169-2179, 2020.
Bollt, explaining the surprising success of reservoir computing forccaster chaos? The universal machine learning dynamical system with contrast to VAR and DMD, Chaos: An Interdisciplinary Journal of Nonlinear Science, 31(1):013100 108, 2021.
Broomhead, King, G. Extracting qualitative dynamics experimental data. Physica Nonlinear Phenomena, 20(2-3):22 (2-23))221--2 212-226,, 1986.
Cancelliere, R., Gai, M., Gallinari, P., and Rubini, L.. OCReP: An optimally conditioned regularization for pseudoinversion based neural training. Neural Networks, 71:76-87, 2015.
Casdagli, M. Eubank, S., Farmer, J. D., and Gibson, State space reconstruction in the presence of noise.
Physica D: Nonlinear Phenomena, 51(1-3):52 1991.
Cha. S.-H. Comprehensive survey distance/similarity measures between probability density functions. City, 1 (2):1, 2007.
Chatzis, S. P. and Demiris, Y. Echo state gaussian process. IEEE Transactions on Neural Networks, 22(9): 1435-1445, 2011.
Chen, H., Tang. F., Tino, P., and Yao, X. Model-based kernel for efficient time series analysis. Proceedings of the 19th ACM SIGKDD international conference on Knowledge discovery and data mining, pp. 392 400, 2013.
Chevyrev, I. and Kormilitzin, A. A primer on the signature method machine learning. arXiv preprint 2016.
Cuturi, M. Fast global alignment kemels, Proceedings of the 28th international conference on machine learning (ICML-II), pp. 929-936, 2011.
Cuturi, M. Doucet, Autoregressive kernels for time series, arXiv preprint arXiv:1101.0673, 2011.
Deyle, E. R. and Sugihara, G. Generalized theorems for nonlinear state space reconstruction. Plos one, 6(3): e18295, 2011.
Ding, H. Trajcevski, G., Scheuermann, P., Wang, X., and Keogh, E. Querying and mining of time series data: experimental comparison of representations and distance measures. Proceedings of the VLDB Endowment, I (2):1542-15 1552, 2008.
Dong, J., Ohana, R., Rafayclyan, M., and Krzakala, E. Reservoir computing meets recurrent kernels and structured transforms. Advances in Neural Information Processing Systems, 13:16785-16796. 2020.
Dong, J., Börve, E., Rafayelyan, M. and Unser, M. Asymptotic stability in reservoir computing. In 2022
International Joint Conference Neural Networks (IJCNN). pp. 01-08. IEEE, 2022.
11

---

Echihabi, K., Zoumpatianos. K., Palpanas, T., and Benbrahim, H. The leraaean hydra of data series similarity search: An experimental evaluation of the state of the art. arXiv preprint arXiv:2006. 11554, 2020.
Franz, M. O. and Schölkopf. B. A unifying view of wiener and volterra theory and polynomial kernel regression. Neural compuiation, 18(12):3097331) 3118, 2006.
Gallicchio, C. Chasing the echo state property. European Symposiun on Artificial Neural Networks, Computa- tional Intelligence and Machine Learning, 2019.
Gauthier, D. J., Bollt, E., Griffith, A., and Barbosa, W. A. Next generation reservoir computing. Nature communications, 12(1):1-8, 2021.
Gauthier, D. J., Fischer, I., and Röhm, A. Learning unseen coexisting attractors. Chaos: An Interdiscipliaary Journal of Nonlinear Science, 32(1 1):1 13107, 2022.
Gonon, L. and Ortega, J.PP. Fading memory echo state networks are universal. Neural Networks, 138: 10-13, 2021.
Gonon, I.,, Grigoryeva, I.,, and Oriega, J.-P. Approximation hounds for random neural networks and reservoir systems. The Annals of Applied Probability, 33(1):28-60 (9))2228669906 2023.
Hamilton, J. D. Time series analysis. Princeton university press, 2020.
Hart, A., Hook, J., and Dawes, J. Embedding and approximation theorems for echo state networks. Neurul Networks, 128:234-247, 2020.
Hart, A. G., Hook, J. L., and Dawes, Echo state networks trained by Tikhonov least squares are (μ) approximators of ergodic dynamical systems. Physica D: Nonlinear Phenomena, 421:132882, 2021.
Ismail Fawaz, H., Forestier, Weber, J., Idoumghar, L., and Muller, P.-. Deep learning for time series classification: a review. Data mining and knowledge discovery, 33(4):917-95) 963, 2019.
Jaakkola, S., Dickhans, M., Haussler, D., et al. Using the Fisher kernel method to detcct remote protein homologics. ISMB, 99:149-158, 1999.
Jaeger, The "echo state" approach analysing and training recurrent neural networks-with erratum note. Bonn, Germany: German National Research Center for Information Technology GMD Technical Report, 148 (34):13, 2001.
Janati, H., Cuturi, M., and Gramfort, A. Spatio-temporal alignments: Optimal transport through space and time. International Conference on Artificial Intelligence and Statistics, pp. 1695-1704, 2020.
Jebara, Kondor, R., and Howard, A. Probability product kernels. The Journal of Machine Learning Research. 5:819-844, 2004.
Kantz, H. Schreiber, T. Nonlinear time series analysis, volume 7. Cambridge university press, 2004.
Kennel, M. B., Brown, R., and Abarbanel, H. D. Determining embedding dimension for shase-space reconstruc- tion using geometrical construction. Physical review 45(6):3403, 1992.
Kim, S.-J., Koh, K., Boyd, S., and Gorinevsky, D. & trend filtering. SIAM Review, 51(2):339 2009.
Kugiumtzis, D. State space reconstruction parameters in the analysis of chaotic time series-the role of the time window length. Physica D: Nonlinear Phenomena, 95(1):13 28, 1996.
Løkse, S., Bianchi, F. M., and Jenssen, R. Training ccho ccho state networks with regularization through dimensionality
reduction. Cognitive Computation, 9:364-378, 2017.
Lu, Z., Leen, T. K., Huang, Y., and Frdogmus, D. A reproducing kernel hilbert space framework for pairwise time series distances. In Proceedings of the 25th international conference on Machine learning, pp. 624-631 2008.
Lukoševičius, M. A practical guide to applying echo state networks. In Neural networks: Tricks of the trade, pp. 659 686. Springer, 2012.
Marteau. P.-F. and Gibet. S. On recursive edit distance kernels with application to time series classification. IEEE transactions on neural networks and learning systems, 26(6): 121-1133, 2014.
Mikalsen, K. Ø., Bianchi, F. M., Eeguero-Ruiz C., and Jenssen, R. Time series cluster kemel for learning similarities between multivariate time series with missing data. Pattern Recognition, 76:569-581 2018.
12

---

Packard, N. H. Crutchfield. J. P., Farmer, J. D., and Shaw. R. S. Geometry from time series. Physical review letters, 45(9):712, 1980.
Paparrizos, J. and Franklin, M. Grail: efficient time-series representation learning. Proceedings of the VLDB Endowment, 12(11):1762-177 2019.
Paparrizos, J. and Gravano, L. k-shape: Efficient and accurate clustering of time series. Proceedings of the 2015 ACM SIGMOD international conference management data, pp. 1855-1870, 2015.
Paparrizos, J., Liu, C., Elmore, A. J., and Franklin, M. J. Debunking four long-standing misconceptions of time-series distance measures. Proceedings of the 2020 ACM SIGMOD International Conference on Management of Datu. pp. 1887-1905, 2020.
Pascanu, R., Mikolov, T., and Bengio, On the difficulty of training recurrent neural networks. International conference on machine learning, pp. 1310-1310 1318, 2013.
Sauer, T., Yorke, A., and Casdagli, M. Embedology. Journal statistical Physics, 65(3):57--616, 1991
Schölkopf, B., Herbrich, R., and Smola, A. A generalized representer theorem. International conference on computational learning theory, pp. 416-426, 2001.
Schölkopf, B., Smola, A. J., Bach, F., et al. Learning with kernels: support vector machines, regularization,
optimization, and beyond. MIT press, 2002.
Shahi, S., Fenton, F. H., and Cherry, E. M. Prediction of chaotic time series using recurrent neural networks and reservoir computing techniques: comparative study. Machine Learning with Applications, 8: 8:100300, 2022.
Small, M. and Tse, C. K. Optimal embedding parameters: modelling paradigm. Physica D: Nonlinear Phenomena, 194(3-4):23-226, 2004.
Steiner, P., Jalalvand, A., and Birkholz, P. Cluster-based input weight initialization for echo state networks. IEEE Transactions on Neural Networks and Learning Systems, 2022.
Steinwart, I. and Christmann, A. Support vector machines. Springer Science & Business Media, 2008
Strogatz, S. H. Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering. CRC press, 2018.
Takens, F. Detecting strange attractors in turbulence. Dynamical systems and turbulence, Warwick 1980, pp. 366-381, 1981.
Tavenard, R., Faouzi, J., Vandewiele, G., Divo, F., Androz, G., Holtz, C., Payne, M., Yurchak, R., RuBwurm, M., Kolar, K., and Woods, E. Tsleam, a machine learning toolkit for time series data. Journal of Machine Learning Research, 21((18):1166 2020. URL http::/jmmmm wecpllloorrggggoe wttg:/pppes//v1/20--eoe 001.09991tttmmm..
Tran, Q. H. and Hascgawa, Y. Topological time-series analysis with dclay-variant embedding, Physical Review E, 99(3)))32209, 2019.
Wallot, S. and Mønster, D. Calculation of average mutual information (AMI) and false-neares neighbors (FNN) for the estimation of embedding parameters multidimensional time series Matlab. Frontiers in psychology, 9:1679, 2018.
Werbos, P. J. Backpropagation through time: what does and how to do Proceedings of the IEEE, 78(10): 1550-1560, 1990.
Whitney, H. Differentiable manifolds. Annals of Muthematics, pp. 645-680, 1936.
Yang. K. and Shahabi, C. A ACAA-asse similarity measure for multivariate time series. Proceedings of the 2nd
ACM international workshop on Multimedia databases, pp. 65-74, 2004.
Zeroual, A., Harrou, F., Dairi, A., and Sun, Y. Deep learning methods for forecasting COVID-19 time-series data: A comparative study. Chaos, Solitons & Fractals, 140:1 10121, 2020.
Zhang, Z.. Zhang, K., and Khelifi, A. Multivariate time series analysis in climate and environmental research. Springer, 2018
Zhang, Z., Zhu, Y.. Wang, X., and Yu, W. Optimal echo network parameters based behavioural spaces. Neurocomputing. 503:299-313, 2022.
13

---

Appendix
A Takens' delay embeddings
This section introduces the field of state space reconstruction in dynamical system theory (Strogatz 2018 and deepens the theoretical interpretation of Sec. [.2] As in Deyle & Sugthara (201T), we start by defining some key concepts that might be useful to introduce the reader to the subject:
a manifold is, intuitively, generalized surface M C IR- that locally looks like an Euclidean space. More rigorously, locally diffcomorphie to some IR";
dynamical system Φ M M is diffeomorphism defining, for discrete times, trajectories on manifold M;
an embedding is a 1:1 smooth map J': : M, - M2 between two manifolds. Embeddings preserve topological properties, e.g. resolves all trajectories without crossings, and geometrical invariants, such eigenvalues and fixed points;
an observation function y : M + Rº is any function that assigns a vector to each point on M.
It is generally truc that can imagine observed MTS X RTXD as realization of complex and unknown smooth dynamical system o" acting on states x" on an unknown manifold M";
x = x² (x()
The observed and original values are related by an observation function y" : :M - R that projects the underlying states onto what we observe:
(+
(8)
This projection can drastically reduce the available information, but one can wonder to which extent it is possible recover the original system when only the observed vectors (x,) are accessible. We focus here the state reconstruction problem, which is that of finding a valid embedding f .M - IR" of the original manifold M
S, = f(x))
This is, by using the observed states (x+). creating richer states (s:) for which the dynamics is topologically equivalent to the unknown ø² In fact, if f is an embedding, then smooth dynamics & is induced on the space of reconstructed vectors:
(s.) I 8* of (s:)
(10)
and is equivalent to Φ² If we succeed, instcad of the observed x, we can use the reconstructed vectors as a better representation of the underlying evolving system x".
Whitney's 's thcorem (Whiney 1936) states the general result that every compact m-dimensional manifold admits embedding in R 2 For dynamical systems this means that, set of 2m + observation functions, where
y > R. is sufficient construct valid embedding:
S: = [y) (x.) ), (x(x))..... 92m=41 (x])]
(11)
the case scalar observations (D=1), Takens' theorem (Takens 1981) considers a single observation function (y), i.e. the projection on the only observed component, but acting on multiple lagged versions of the system:
(..ee)))))))) v(0² (x+)). v( (x,)
(12)
[ -e-r: ..... [mm]
The idea behind Eq. 12 is that the information about the unobserved state variahles is contained in the relations between past and future values of the observed time series and 2m + 1 dimensional embedding (often less) sufficient to extract that.
Generalized versions of Takens' theorem have also been formulated. In particular, we find in Deyle & Sugiharr (201T) an extension to cases where lags of multiple observation functions can be used to create the embedding. This is useful result in cases where one cannot consider many lags, but can instead make use of fewer lags spread across multiple available dimensions. In more detail, the allowed set of observation functions any {yx 12m+11 containing terms {y. that cannot be written as a lagged version of any other element in the set fk=1
and terms (v,) } that satisfy y4 = yv (ō (x*)) for choice of b. Both Takens" and Whitne'' 's thcorems can be recovered from this formulation: reduces to the case of Takens' when only one {g-} = y is chosen and all the others belong to the family of (vv): it reduces to the case of Whitney's 's when {u- = {uu } and {yr} = 0.
All terms in the NVAR representation are reminiscent of this formulation. This points to the direction of RNV being closer approximation of the underlying state than the initial data X and also suggests that the final number of dimensions plays more important role than which are included and which not. In practice the dimensionality m the original manifold M* unknown, making the selection the embedding dimensionality a potentially rich field of study. On top of that, the observed data may be corrupted by multiple noise sources, which increases the importance of sampling the series correctly with meticulous choice of the lag size (r s).
14

---

B Experimental evaluation details
B.1 Datasets
Details of the 130 Univariate and 19 Multivariate datasets reported in Tab, Both extensive range data points series lengths. From UEA archive, we have excluded dige-diimnsional datasets NVAR kernel is not appropriate, and 5 datasets which SVM accuracy was < 0.5 for all approaches, which indicates that kernel methods are, in general, not a suitable solution.
which the application
Table Attributes of UTS and MTS datasets. The splitting between the training and testing is provided by the respective archive.
UTS DATASET
Nerain
Ntest
D
N Nelass
130 DATASETS
16 + 8926
20 + 16800
15 + 18530
2 + 60
MTS DATASET
Ntraón
Ntest
T min
T Tmaz
D
N
class
SPOKENARABICDIG.
6599
2199
93
13
10
APPNESEVOWELS
270
370
29
12
9
PENDIGITS
7494
3498
2
10
RACKETSPORTS
151
152
30
30
6
4
LSST
2459
2466
36
6
LIBRAS
180
180
45
45
2
15
FINGERMOV.
316
100
50
50
28
2
NATOPS
180
180
51
51
24
6
CHARACTERTRAJ.
1422
1436
60
182
3
20
ERING
30
270
65
65
4
6
BASICMOTIONS
40
40
100
100
6
4
ARTICULARY
ARTICULARYWORDRICC WORDREC.
275
300
144
144
9
25
EPILEPSY
137
138
206
206
4
UWAVEGESTUREELIB
120
320
315
315
3
8
SELFREGULATIONSCHI
268
293
896
896
6
2
SELFREGULLTIONSCF2
200
180
1152
1152
7
2
CRICKET
108
72
1197
1197
6
12
STANDWALKJUMP
12
15
2500
2500
4
3
EIGENWORMS
128
131
17984
17984
6
5
B.2 Implementations and details for reproducibility
For reproducibility, list here the parameter setting and the code implementation for each method. Experiments were run on 32-core AMD Ry Threadripper CPU.
NVAR kernel (NVARK)
the allowed number of lags l: v /average distance between peaks and hollows, after series have heen & trend-filtered (Kim et al. 2009
lag size s l:
order of the polynomial functionals 2;
maximum dimensionality embedding d. 75;
readout regularization set by the OCReP method (Cancellicre et al. 2015):
RBF lengthscale the median pairwise distance between the representations.
An casy-to-use Python implementation of the NVAR kernel is made publicly available at bhttps://githu com/ com/ http://fffooaarrkkeeeeeee
15

---

Table 5: Grid values for the CV optimization of the NVARK embedding parameters.
DATASET
CV GRID
T< 400
k=11,2,3,4,5,kk s-[12,2333
T2 100
K= (=11223345.022000g05=515.1..
SPOKENARABICDIGITS
k= [=11,2.3,4,k* (=1,,,33,,4
JAPANESEVOWELS
PENDIGITS
RACKETSPORTS
LSST
LIBRAS
FINGERMOVEMENTS
NATOPS
CHARACTERTRANECTORIES
K= K=11.233,,110,20.1 |= (=[1,5,00,s*)
ERING
BASICMOTIONS
EPILEPSY
UWANEGESSTURELIBBARY
SELFREGGULAATIONNSCPP
SELFREGULATIONSCP2
CRICKET
STANDWALKJUMP
EIGENWORMS
optimized NVAR kernel (NVARK*)) The lag size and the allowed number of lags l are optimized in a
CV loop scrolling small grid. Values for the grid are reported in Tab. ke and s" refer to the values obtained from the rule of thumb.
In general, one should allow for larger values for k and in the case of longer series. A higher it allows capturing longer-term dynamics, while larger avoid oversampling the series. In the case multivariate datasets, observed that many lags are hardly required. This is possible because few lags spread across different dimensions can be sufficient to construct the embedding. Also, the pool of features (Ring Roontin]) largely increases with the dimensionality of the dataset and k. which affects the stability of the approach.
reservoir-model ESN (rmESN) For the mESN, we follow the setting given in Bianchi et al. 2020
number of internal reservoir units R 800;
spectral radius = 0.99;
leakage;
percentage of non-zero connectivity within the reservoir states 8 = 0.25;
input scaling 0.15:
noise level in the reservoir update = 0.001;
transient dropped states Narop = 5; number of dimensions after the dimensionality reduction module D 75 (via tenor-PCA, also proposed the same paper);
favoring comparability with NVARK, we set the readout regularization by the OCReP method; favoring comparability with NVARK, we set the RBF lengthscale to the median pairwise distance between the representations.
As implementation, we used the public Python code provided by the authors: https://github com/ atw.rriiooonnnnddcttrrrinnn ceadd-lusteennngeee wwtthngwwthh--hh Reservoir Reservoir-Conputiing Comput-ing
Fil SeteooM///immmsereesee wtw.eeess-clssfffaatttiicoocoeee and
16

---

Time Cluster Kernel (TCK) For TCK, we follow the guidelines given in the original paper and in the provided code. For dataset of N samples and D attributes, these are:
minimum percentage subsample minN 0.8;
min number of attributes for each GMM min D for univariate datasets, min 2 for multivari-
datasets;
max number of attributes for each GMM maa = min(ceil D),
min length of time segments for each GMM minT
max length of time segments for each GMM maxl min(floor * 25);
number of mixture components for each GMM C if < 100, C 40 otherwise;
number of randomizations for each number of components G = =30
number of iterations I 20.
The Matlah implementation TCK that we use is part of the code available hhttps///githu com/ FilippoMB/TCK
Global Alignment Kernel (GAK) We set the multiplicative factor for the bandwidth of the radial function to OGAK = 2. Python code is available from the tslearn package (Tavenard et al. 2020). As reported in the package documentation, the method can lead to numerical issues for long series.
Shift INvariant Kernel (SINK) For SINK, adopted the Python implementation provided by the author /tttpa///githbb com hean/TeeDatunOrrggreaaaa petpaaatuunngggral-pyyt.ia We use r and retain all the energy after courier-transforminng the series.
C Additional experiments
C.1 Feature redundancy
For studying whether previously adopted values for the embedding dimensionality (d,) are also suitable for NVARK, we adopted an experimental approach and studied how many features can be included before experiencing redundancy. We do by plotting, for few univariate and multivariate datasets, the accuracy of SVM classifier against the d, parameter (Fig. While do not observe redundancy in univariate datasets,
Arrowitedd (805)
Scredhnn Mil
GundPoar (777)
J.comeentoowee 12/2)
Libra (MF
HouseTwenty 1619)
-ricket (7P)
wowooooatttebbbaay (717)
Tegers/oveeents [V1I
Ecllppy (23)
SVM accuracy
0.8
SVM acouracy
0.8
0.7
0.7
0.5
0.5
10
5 7 %go
NVARenbeddinng dimensionallty
embedding dimensionslity
Figure For some univariate (left) and multivariate (right) time series datasets, the plots show the accuracy of an SVM classifier against the NVAR embedding dimensionality dr. Embedding parameters are indicated as (k/s) and determined via the proposed rule of thumb. Features redundancy manifests for some multivariate datasets (Cricket, UWWaveGestureLibrary as an accuracy drop.
we do observe drop in performance after different d, for some multivariate datasets. We identify the region 70 d. 100 gond, though non-optimal, setting, which is in line with previous RC studies (Bianchi 2020
17

---

C.2 Sensitivity analysis
In this section, we address the sensitivity of NVAR to the different hyperparameters:
0.
0.0
088
0.M
0.
180
0.24
devedeentl
110
300
17
10
Teur
10
Figure 4: Sensitivity of NVARk to the hidden dimensionality (left), evaluated across 16 multivariate datasets; sensitivity to the regularization parameter (center) and RBF lengthscale (right), evaluated across 20 randomly sampled univariate datasets. Red arrows indicate our parameter's choice.
d,: show in Fig. (left), the average performance for multivariate datasets not considerahly impacted the range d, [70, 150). As expected, we observe significant drop in performance for high values (redundancy).
Arige: show in Fig. [4 (center) little sensitivity in the region Aridge E [5, 50). Results averaged across pool of 20 randomly sampled univariate datasets. The employed OCReP method (Cancelliere ct al. 2015) outperforms the best choice by ~ 2%, as adapts the value to each individual dataset.
yer in Fig. 4 (right), we range different multiplicative factors for the median pairwise distance between all representations and show little sensitivity in neighborhood of our choice "Khm 1.
& and 8: The sensitivity to the embedding parameters strongly depends on the dataset, although generally tends to be high. In response, we introduce our dataset-specific heuristic adapting idcas from precedent study (Kugiumtzis (996). which we find works well in practice. A reasonable idea of the average sensitivity these parameters NVAR, where & and s are chosen by the heuristic, to the performance of NVARK*, where and are optimized. Unfortunately, understanding sensitivity by considering other settings is challenging, as the existing literature predominantly focuses on the noise- free case. Furthermore, there is no strong evidence for the wider applicability of such approaches outside dynamical systems.
be obtained by comparing the performance of
C.3 Uneven lengths
NVARk can process and compare time series of different lengths. We report in Tah. 6 the SVM accuracy scores for three datasets with uneven lengths.
Table SVM classification accuracy for datasets with uneven lengths marked with the [dL] notation. For comparison, here are copied the respective results from Sec. 4.3 after zero padding preprocessing,
DATASET
NVARK"
NVARK
SPOKENARABICDIGITS [DL]
0.971 003
0.971
SPOKENARABICDIGITS
0.980 004
0.980
JAPANESE [DL]
0.973
0.975
JAPANESE VOWELS
0.981 008
0.97510.007
CHARACTERTRAAECTORESS [DL]
0.983 005
0.93340.000
CHARACTERTRAJECTORNE
0.982 000
0.915
C.4 Comparison with Global Alignment Kernel
We compare here the performance of the NVAR kernel against the popular Global Alignment Kernel (GAK) baseline from Cuturi (20IT). This has been moved from the main text as authors of SINK have already demonstrated superior performance to GAK, following similar evaluation. As reported in the tslearn package documentation, the provided implementation
lead to numerical issues when used with relatively long time series. To face this, we evaluate GAK on time series linearly interpolated this threshold. We report in Tab, 7 the average results across 122 univariate datasets. We report in Tab, B results for the multivariate datasets. Datasets interpolated to T=400, before being processed by GAK, are marked with the symbol (*).
400 when their length exceeds
18

---

Table 7: Average performances for NVARK and GAK across 130 univariate datasets.
NVARK*
NVARK
GAK
AVERAGE SCORE N FIRST RANKED
0.803
0.779
0.728
77
56
Table 8: SVM classification of MTS datasets. Average across 10 seeds. Best accuracy is in bold.
DATASET
NVARK
GAK
SPOKENARARICDIC
0.980 =000000000000
time limit
JAPANESEVOWELS
0.975
0.974 002
PENDIGITS
0.983_00000
0.978
RACKETSPORTS
0.847 0.847=0.000 20.007
0.83620.000 20.008
LSST
0.543 =0.00
0.500 000
LIBRAS
0.900
0.783
FIVGERMOVEMENTS
0.581 00032
0.565 035
NATOPS
a..19
0.918 000.000
CHARACTERTRAJ.
0.915 008
0.988 ::0 000
ERING
0.861
0.931 2: 002
BASICMOTIONS
0.938 0.0027
0.995 011
ARTICULARY WORDDRE...
0.976
0.980 000
EPILEPSY
0.986 0.986=0.000
0.845+0.000
UWAVEGESTURELIB.
0.939 =0.011
0.865 002
SELFREGULATIONSCP
0,688
0.902 1.0000
SELFREGULATIONNSCEE
0.567 #0.028
0.526 1:0.01
CRICKET
0.986 00000
0.958 4.0.000
WALKJUM
0.607 =0.000 049
0.333 172
WORMSW
0.950 =0.012
0.48220.005
AVERAGE SCORE
0.848
0.756
N FIRST RANKED
13
6
C.5 Example kPCA visualization
Despite the underlying equivalence between NVAR and simple RC, better performance be obtained by finding better parameter setting, which not trivial for RC. In our approach, better setting easier to find with fewer and more interpretable parameters. Here, we aim to visually highlight the possihle extent differences between NVAR and mESSN i.c., theoretically equivalent methods under diffcrent working points. To do projected the data onto the first principal components using different kernel matrices (kPCA). When inspecting the projections for different datasets, observed few critical differences, predominantly for longer series. One example illustrated Fig. 5 CinCEEGTTo dataset). These appear in form a qualitatively
NVARK
NVARK
MESN
CinCECGTorso
Figure 5: kPCA visualization of CinCECGTorso dataset. Projection onto the first three principal components obtained with different kernels. Different colors correspond to different classes.
better grouping and class separation in the projected space. Tracing back the origin of this improvement, e.g.,
19

---

better working point better resilience to noise, would represent interesting future work. For an alternative visualization, display here the two-dimensional projection of the CinCEGGTorso dataset using different kernels.
D Individual results for SVM classification of UTS
We report here the individual accuracies obtained under the SVM classification framework of univariate time series of Sec. [..] Each accuracy is obtained as mean 10 repetitions with different random seeds, readability, the standard deviation over the runs is not reported. Datasets interpolated to T=400, before being processed by GAK. are marked with the symbol (*).
Table 9: SVM classification of 130 UTS datasets, individual performance. Results are averaged
across 10 different seeds.
NVARK*
NVARK
rmESN
TCK
GAK
SINK
ACSFI
0.835
0.826
0.792
0.646
0.396
0.800
AbnormalHeartbaat
0.707
0.737
0.688
0.733
0.716(*)
0.730
Adiac
0.749
0.709
0.638
0.623
0.598
0.731
X
0.673
0.676
0.608
0.471
0.527(
0.584
AllGesture Wiimote
0.733
0.631
0.693
0.444
0.548(
0.651
WiGestuzeWimoteZ
0.674
0.625
0.591
0.452
(.505(#)
0.622
ArrowHead
0.823
0.834
0.577
0.672
0.794
0.842
BME
0.973
0.927
0.983
0.967
0.953
0.971
Beef
0.847
0.590
0.783
0.610
0.897(#)
0.933
BeetleFly
0.800
0.950
0.880
0.645
0.850(#)
0.850
BirdChicken
0.850
1.000
0.785
0.590
0.700(#)
0.800
CBF
0.979
0.958
0.894
0.968
0.925
0.994
Car
0.833
0.733
0.612
0.693
0.825
0.800
CatsDogs
0.737
0.726
0.709
0.663
0.626(
0.725
Chinatown
0.950
0.962
0.944
0.985
0.973
0.969
ChlorineConcentration
0.868
0.701
0.847
0.568
0.601
0.879
CinCECGTorso
0.914
0.849
0.834
0.850
0.888(*)
0.749
Coffee
0.964
0.964
1.000
0.971
0.968
1.000
Computers
0.776
0.685
0.712
0.658
0.527(+)
0.595
CricketX
0.697
0.645
0.649
0.618
0.681
0.744
Yricket
0.682
0.554
0.586
0.656
0.635
0.723
CricketZ
0.697
0.680
0.656
0.647
0.658
0.761
Crop
0.731
0.703
0.753
0.729
time limit
time limit
DiatomSineReduction
0.968
0.912
0.858
0.911
0.964
0.977
DistaPhallaxxOOttieeeggCCooo
0.763
0.741
0.718
0.750
0.705
0.719
DisalPhallaxXOuiineCCrre..
0.779
0.775
0.733
0.787
0.740
0.764
DistalPhalanXTW
0.705
0.676
0.644
0.654
0.675
0.655
DodgerLoopDay
0.462
0.392
0.448
0.679
0.635
0.575
DodgerLoopGame
0.790
0.776
0.774
0.825
0.871
0.868
DodgerLoop Weekend
0.957
0.920
0.914
0.986
0.986
0.978
ECG200
0.880
0.890
0.877
0.895
0.895
0.846
ECG5000
0.940
0.940
0.940
0.938
0.939
0.943
ECGFiveDays
0.984
0.995
0.878
0.831
0.997
1.000
EOGHorizontalignnal
0.566
0.573
0.538
0.419
0.422
0.512
EOGVerticalSignal
0.479
0.453
0.283
0.323
0.373
0.417
Earthquakes
0.748
0.748
0.748
0.748
0.750 (+)
0.748
Devicess
0.654
0.650
0.652
0.643
time limit
time limit
EthanolLevel
0.579
0.360
0.587
0.382
0.607(*)
0.562
FaceAll
0.772
0.775
0.856
0.749
0.884
0.768
FaceFour
0.784
0.540
0.519
0.925
0.857
0.830
UCCesUCR
0.903
0.811
0.808
0.852
0.890
0.895
Fifty Words
0.562
0.565
0.659
0.676
0.720
0.751
20

---

NVARK
NVARK
imESN
TCK
GAK
SINK
Fish
0.960
0.966
0.811
0.769
0.857(+)
0.906
FordA
0.982
0.928
0.945
0.674
0.910())
0.925
FordB
0.868
0.800
0.825
0.583
0.759(*)
0.777
FreezerReggulrTTriin
0.989
0.956
0.964
0.963
0.985
0.980
FraezerSmallTrain
0.836
0.905
0.822
0.743
0.691
0.700
Fungi
0.656
0.952
0.866
0.930
0.982
1.000
G
0.662
0.591
0.535
0.448
0.535
0.554
GestureMidAArrD2
0.546
0.565
0.452
0.295
0.538
0.483
G
0.238
0.248
0.204
0.177
0.290
0.215
GesturePPablee21
0.773
0.573
0.591
0.739
0.693( *
0.725
Gesture/Pbblle22
0.722
0.551
0.594
0.565
0.513(+)
0.618
GunPoint
0.987
1.000
0.963
0.961
0.960
0.980
GeppointAgeSpaa
0.972
0.994
0.979
0.971
0.946
0.936
0.975
0.997
0.992
0.984
0.996
0.997
GunPointOld Young
1.000
1.000
1.000
1.000
1.000
1.000
Ham
0.717
0.739
0.707
0.730
0.748(*)
0.714
HandOutlines
0.948
0.945
0.866
0.884
0.917(*)
0.914
Haptics
0.525
0.399
0.437
0.444
0.445(
0.445
Herring
0.625
0.634
0.594
0.577
0.594(*)
0.605
House Twenty
0.916
0.923
0.900
0.893
0.785(+)
0.815
InlineSkate
0.543
0.470
0.443
0.310
0.340(#)
0.401
InsectEPGRegular
1.000
1.000
1.000
1.000
1.000())
0.956
InsectEPGGSmal/Trainn
1.000
1.000
1.000
1.000
1.000(*)
0.936
InsectWingbeatSound
0.312
0.372
0.515
0.638
0.632
0.617
ItalyPoweDDemand
0.963
0.952
0.952
0.950
0.959
0.949
LargeKiichenAAplianees
0.755
0.690
0.887
0.552
0.512(*)
0.794
Lightning2
0.738
0.725
0.811
0.715
0.797(*)
0.820
Lightning7
0.630
0.612
0.759
0.681
0.699
0.781
Mallat
0.959
0.936
0.909
0.939
0.924(*)
0.937
Meat
0.950
0.760
0.912
0.952
0.930(+)
0.933
Medicallmages
0.713
0.742
0.726
0.698
0.761
0.740
Medetranne
0.931
0.905
0.950
0.938
0.935
0.933
MiddlePPallanxxuttinneegeee
0.591
0.630
0.578
0.574
0.631
0.631
MiddeePhallnnOuutiieCCoreett
0.866
0.834
0.779
0.786
0.632
0.811
MiddIeePhalaxxTW
0.552
0.558
0.560
0.588
0.580
0.594
MixedShapesReggulrTTrie
0.946
0.936
0.913
0.874
0.903(*)
0.923
MixedShapesSmallTraii
0.913
0.867
0.798
0.817
0.846(*)
0.885
MoteStrain
0.812
0.792
0.694
0.882
0.861
0.870
VonlnvasiveeetalECGGG
0.930
0.924
0.935
0.890
0.922(+)
0.937
NonlnvasvveeetulECCGTTooaaa
0.943
0.935
0.938
0.903
0.942(*)
0.950
OSULeaf
0.955
0.909
0.845
0.475
0.585(*)
0.707
OliveOil
0.867
0.867
0.793
0.853
0.407(*)
0.700
PLAID
0.737
0.749
0.521
0.424
0.396(+)
0.393
MhalangesOutinesCCoree
0.826
0.834
0.784
0.797
0.664
0.798
Phoneme
0.297
0.267
0.335
0.221
0.117(*)
0.193
PickupGesture WiimoteZ
0.600
0.560
0.680
0.550
0.642
0.614
PigsirvayPressuree
0.856
0.702
0.226
0.109
0.122(*)
0.216
PigArPrresuure
0.954
0.962
0.514
0.210
0.248(*)
0.779
PigCVP
0.926
0.885
0.815
0.152
0.173(*)
0.601
Plane
1.000
1.000
0.981
0.990
0.973
0.968
PowerCons
0.883
0.873
0.943
1.000
1.000
0.983
ProximaPPialannOOttiiAAAgggon
0.820
0.849
0.849
0.853
0.840
0.837
MoxiimalPhallnxxQutineCCCCC
0.900
0.900
0.858
0.855
0.778
0.849
ProximaIPhalanxTWW
0.771
0.795
0.763
0.801
0.802
0.797
RefrigerationDevices
0.541
0.523
0.542
0.512
0.419(*)
0.483
Rock
0.644
0.686
0.514
0.546
0.618(*)
0.508
ScreenType
0.552
0.416
0.483
0.402
0.455(*)
0.423
21

---

NVARK*
NVARk
rmESN
TCK
GAK
SINK
Semg2HandGGenderCh2
0.871
0.840
0.829
0.886
0.912(*)
0.909
Semg/HanddMovementhh2
0.717
0.562
0.600
0.756
0.768(
0.829
SemgHHandSuhjecrthh
0.816
0.737
0.711
0.841
0.893(+)
0.913
ShakeGestureWimotez
0.800
0.800
0.842
0.634
0.704
0.680
ShapeletSim
0.950
0.950
0.871
0.537
0.462(#)
0.800
ShapesAll
0.867
0.843
0.782
0.745
0.774(*)
0.834
SmallKKirchenAppllaanes
0.656
0.607
0.706
0.630
0.508(+)
0.684
SmoothSubspacce
0.960
0.945
0.903
0.980
0.943
0.900
SonyAlIORRoosSurficcll
0.945
0.930
0.807
0.936
0.708
0.894
Sony-AlBORobotSurfeee2
0.918
0.928
0.834
0.829
0.816
0.879
Starl.ightCurves
0.976
0.966
0.969
0.937
0.946(*)
0.959
Strawberry
0.951
0.965
0.962
0.966
0.924
0.951
SwedishLeaf
0.922
0.929
0.916
0.907
0.895
0.930
Symbols
0.911
0.945
0.950
0.877
0.897
0.913
SyntheticControl
0.980
0.963
0.909
0.990
0.979
0.994
ToeSegmentation!
0.904
0.970
0.890
0.615
0.635
0.864
ToeSegmentation2
0.862
0.883
0.814
0.778
0.795
0.865
Trace
1.000
1.000
1.000
0.931
0.922
0.960
TwoLeadECG
0.978
1.000
0.871
0.692
0.861
0.978
TwoPatterns
0.988
0.822
0.767
0.997
0.953
0.999
UMD
0.993
0.986
0.856
0.963
0.961
0.932
0.967
0.706
0.727
0.946
0.958(+)
0.965
UWaveGestureL X
0.697
0.698
0.631
0.766
0.772
0.792
UWaveGestureLibrag Y
0.624
0.623
0.630
0.668
0.682
0.689
UWaveGessureLibbaayyy
0.704
0.695
0.621
0.684
0.712
0.735
Wafer
0.993
0.992
0.995
0.993
0.995
0.998
Wine
0.722
0.815
0.794
0.707
0.606
0.667
WordSynonyms
0.434
0.440
0.572
0.596
0.647
0.669
Worms
0.779
0.822
0.808
0.569
0.621(*)
0.670
WormsTwoClas
0.805
0.827
0.805
0.568
0.662
0.717
Yoga
0.831
0.822
0.793
0.795
0.845(*)
0.877
22


Reviews:
[{"review_id": "H2585hAhrK", "reviewer": "Reviewer_RJLc", "summary": "In their paper, the authors propose a new method for deriving a kernel from time series data. They combine ideas from reservoir computing with nonlinear vector autoregressive models, based on recent theoretical work exploring their similarities.\nThe primary idea of the paper is to construct a kernel by modeling a time series' dynamics using a nonlinear vector autoregressive model first. The vectorized fit parameters are then used to constitute the kernel representation of the respective time series.\nThe authors investigate the efficacy of this kernel in SVM classification across a range of unimodal and multimodal time series. They demonstrate state-of-the-art performance when compared to three other approaches.", "questions": "Relationship to Taken's Theorem: The choice of methodology is explained with reference to Taken's delay embedding. However, the connection to Taken's theorem in some of the practical cases isn't entirely clear to me. In the original next-generation RC approach, chaotic attractors are considered, where Taken's Theorem applies, while in your setting, many of the time series used in the classification tasks are not derived from something that can be as straightforwardly modeled via an underlying attractor (e.g., Japanese Vowels). It might be interesting to investigate further for which types of time series your approach is particularly suited and outperforms other methods. It's also worth exploring whether this relates to the suitability of a dynamical systems description for the underlying dynamics.\n\nl.56: The resulting model is more accurate, interpretable, and non-recurrent.\nWhat does interpretable refer to in this context? The single k-PCA example does not make this point strong enough in my mind, and it is not apparent to me how the derived parameters are more interpretable than e.g. the parameters of an RC model.\n\nSmall typos:\n226, implies\n236: only the ridge\n328: Computationally, it is exceptionally", "limitations": "The authors explicitly mention that high-dimensional datasets are the main weakness of their approach. I don't foresee any direct negative societal impact.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Overall, the paper is well-written and easy to follow. The main ideas and contributions are clearly delineated. Further, code is provided with the submission. The method manages to significantly reduce the number of hyperparameters needed to model dynamics compared, for instance, to RC approaches. It also achieves state-of-the-art performance on a range of unimodal and multimodal time series classification tasks, compared to other kernel based methods. It presents an interesting contribution that as far as I can judge is novel, and leverages a recent development in a related field to come up with a new application.", "weaknesses": "Comparison methods: If I understood correctly, no hyperparameters were tuned for the comparison methods in the experiment, while you optimized the delay-embedding parameters for your method. While I appreciate that a thorough hyperparameter tuning for some of the comparison methods can be computationally prohibitive and the smaller amount of hyperparameters is one of the strengths of your approach, it would make the comparisons fairer if some tuning took place, or if you explained in more detail why you didn't perform any tuning.\n\nLine 217: \"We propose the region 75 ≲ ¯dr ≲ 100 as a good, though non-optimal, setting\" The cutoff around 75-100 is justified experimentally in the appendix. The phrase \"non-optimal setting\" is repeated in the appendix. What do you mean by \"non-optimal\" in this context? Why would you choose a non-optimal setting? While I understand that this cutoff is motivated by decreased performance due to redundancy, it seems highly dependent on the dataset, and I'm not sure I understand why this reduces performance for the MTS so drastically. Since this is one of the central components/hyperparameters of your algorithm, some more clarification would be useful here.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We warmly thank reviewer $\\color{red}{\\textbf{RJLc}}$ for the thorough review and positive feedback.\nWe really appreciate the kind words on the clarity of our paper, as well as the novelty and innovative utilization of recent developments from a related field.\nLastly, we are also grateful for spotting a few typos.\nAs for the raised concerns and questions, individual points are addressed below.\n\n\n**W1. No hyperparameters were tuned for the comparison methods**\\\nThis is only partially true.\nWe indeed found fine-tuning baseline hyperparameters to be mostly computationally unfeasible and, considering NVARk's heuristic foundation, might also not yield a fair comparison.\nHowever, to the best of our knowledge, there are no rule-of-thumb guidelines facilitating the choice.\nOur choice of baseline hyperparameters is then based on what is provided in the authors' original papers, with the belief that suit their approaches well and are transferable, especially given that our task and benchmarking archives are similar (or a subset).\nWhen applicable, we consistently apply the same choices for both NVARk and rmESN, i.e. OCReP optimization for $\\lambda_{ridge}$ and heuristic for $\\gamma_{rbf}$.\nAlso, note that our choice for the hidden dimensionality ($\\bar{d}_r = 75$) is optimized for rmESN in the authors' paper (see Supplementary Material for [a]), and is not optimal for NVARk.\nFinally, we would like to emphasize that performance for optimized NVARk (NVARk*) is reported separately.\n\n\n**W2. Why would you choose a non-optimal setting for $\\bar{d}_r$?**\\\nIn line with the scope of our work, the main motivation for using $\\bar{d}_r = 75$ (l.252) is to prioritize a fair comparison with rmESN.\nOur term \"non-optimal\" refers to the lack of fine-tuning of $\\bar{d}_r$.\nThe fact that this value lies within the range $70-100$, where we empirically observe an absence of feature redundancy, made us accept this value.\nWe have made this clearer at ll.213-217, with a recall at l.252.\nAs for the observed performance drop after this range, please find a more extensive discussion in our common response \\#2.\n\n\n**Q1. It’s worth exploring whether performance relates to the suitability of a dynamical systems description**\\\nWe agree; having found good performance in diverse cases of time series that are not immediately arising from dynamical systems, is one of the most interesting outcomes of our work.\nDespite our paper would certainly benefit from inspecting such aspects, we argue that it would be better suited for dedicated future work, e.g. compelling follow-ups on time series representation learning and state space reconstruction as a possible underpinning.\nWe believe these results, as well as future directions, are of great interest to both dynamical systems and time series communities.\nWe thank the reviewer for highlighting this future direction and we are happy to include it in the Conclusions section of our paper.\n\n\n**Q2. What does interpretable refer to in this context?**\\\nWe kindly refer the reviewer to our common response \\#1.\n\n\nWe hope to have clarified all of the reviewer’s concerns and are happy to provide further details.\n\n**References:**\\\n[a] Bianchi et al., Reservoir Computing Approaches for Representation and Classification of Multivariate Time Series, IEEE Transactions on Neural Networks and Learning Systems, 2020.\n"}, {"review_id": "yi8pfKlaaF", "reviewer": "Reviewer_hHQP", "summary": "This work proposes a feature extraction based on a non-linear vector autoregressive model (called NVAR kernel in the paper). The NVAR method constructs a deterministic feature matrix with the original input time series, lagged versions of the time series (parametrized by the lag and spacing parameters) and the products between lagged and unlagged dimensions, up to some polynomial order. The resulting feature map is then subsampled to obtain a random feature map. A linear readout layer performs next-step prediction on the embedding states and the resulting parameters form the feature map for a radial basis function kernel. The authors form a connection to state space reconstruction theory to justify their approach of subsampling by connecting it to state space reconstruction theory. The method is compared to other kernel choices on univariate (109) and multivariate (16) datasets for classification. The authors also demonstrate the scalability of their approach by measuring execution time and provide heuristics for hyperparameter choices. ", "questions": "See Weaknesses. ", "limitations": "Limitations are not discussed in the paper (or I missed it). A short discussion of limitations would be appreciated. ", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Originality\nThe underlying ideas in this paper have been proposed before in the time series space (lags and polynomial features; randomized feature maps), but the combination proposed here is novel for time series classification. The NVAR formalism in time series classification is novel as well (to the best of my knowledge). \n\nQuality\nThe algorithm choices presented in this paper are motivated by referring to previous work and connecting to state space reconstruction theory. While I’m not an expert on this, I appreciate the effort the authors made to justify the randomized subsampling step. The paper evaluates the proposed approach against other kernel choice baselines from different classes (reservoir, model-based, Fourier-based) univariate and multivariate datasets. The authors also analyze the execution time and scalability. I appreciate the efforts from the authors to provide heuristics for hyperparameter choices of their approach and showing the gain of additional cross-validation. I think the paper misses the comparison to a critical baseline, which I will discuss in the Weaknesses section. \n\nClarity\nOverall, the paper was well written but I think the clarity could be improved by introducing some core concepts into the introduction to readers that are not familiar with reservoir computing. For example, the authors state in the introduction that reservoir computing has a “large set of hyperparameters” but do to mention what they mean. I think a bit more detail to introduce the reader would be helpful (although I would consider this largely a stylistic choice as most of this detail comes later in the paper). \n\nSignificance\nTime series classification is an important application of time series research. Especially in the healthcare domain, time series classification methods have to run on devices with limited compute (smartphones or other edge devices). As such, a classification method that does only require minimal compute is significant. However, since this method is specific fairly specific to time series classification, the broader impact is limited. \n", "weaknesses": "I agree the rationale of the authors to select one representative baseline methods per category (reservoir, model based) etc. However, I think the the MINIROCKET method (Dempster et al. KDD 21: https://dl.acm.org/doi/10.1145/3447548.3467231) shares many properties of the proposed method (per time series computation is the most expensive step, suitable for downstream linear classification, low computational cost). Since this method provides high classification accuracy at low computational cost, I think it should be one of the baselines in the paper (it can also deal with uneven sequence length due to its final pooling operation). I would appreciate if the authors consider to compare against the method. At the very least, I would expect that this method is discussed in the related work (even though it is not from the reservoir computing literature). ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We deeply thank reviewer $\\color{blue}{\\textbf{hHQP}}$ for the extensive review and positive feedback.\nWe really appreciate the positive remarks on the general quality and clarity of our paper.\nSimilarly, we value the highlighting of a form of novelty and significance.\nAs for the raised concerns and questions, individual points are addressed below.\n\n\n**W1. A bit more detail to introduce the reader to RC would be helpful**\\\nWe agree that the paper would benefit from a more extensive introduction to RC and its hyperparameters, as their non-interpretability and large number partially motivate our approach.\nOur decision not to delve into extensive discussions originated from the need to adhere to the space limitations.\nAs a middle ground, we now mention a few representative examples around l.120, i.e. spectral radius and input scaling (control non-linearity and relative effect of the current input as opposed to the history), and leaking rate (control speed of reservoir updates).   \nWe also refer the reader to Appendix B for a comprehensive list. \n\n\n**W2. Broader impact is limited**\\\nWhile demonstrating versatility across tasks is valuable, we respectfully disagree that assessing one task suggests limited impact.\nWithin the established literature, typical comparison of kernels uses classification only [a]-[c], but their application notably extends beyond that.\nAs in Mikalsen et al. [d], we also present a snapshot of kPCA. \nMost importantly, we would like to emphasize that the numerical accuracy of TS classification constitutes only a part of our work. \nRather, we anticipate an impact on RC representation learning, for which we demonstrate superior performance with higher efficiency and interpretability.\nAdditionally, we demonstrate success of extending applicability of NVAR to noisy multivariate time series analysis, which we argue is of interest to both communities.\n\n**W3. MINIROCKET should be one of the baselines**\\\nWe thank the reviewer for bringing MINIROCKET to our attention, we have included it in the background section.\nOur choice of representative for SOTA was based on the most extensive evaluation of univariate TS metrics to date [e], from where SINK comes out as the strongest trade-off between accuracy and efficiency.\nWe agree that MINIROCKET, as a representation learning method, can be used in conjunction with an RBF kernel, which would allow comparison.\nIn fact, this line of work certainly warrants further investigation.\nHowever, we are unable to provide proper comparison due to the time constraints.\nOn one hand, one should ensure fairness of the hyperparameter setting, with no bias towards the UCR classification task (note that all heuristics for NVARk, except arguably $\\bar{d}_r$, are task agnostic).\nOn the other, the NVARk processing is of a causal type, while our understanding is that MINIROCKET attends to both the past and future of each timestamp.\nAs such, we believe that a bidirectional variant of NVAR representation would be better for such comparison.\nOn a side note, the focus of our work is more on the expressive power of RC methods rather than TS classification.\n\n\n**W4. Limitations are not discussed**\\\nWe kindly refer the reviewer to ll.330-331, where the main limitation is mentioned in the manuscript.\nPlease also find a further discussion in our common response \\#2.\n\nWe hope to have clarified all of the reviewer’s concerns and are happy to provide further details.\n\n**References:**\\\n[a] Paparrizos et al., Grail: efficient time-series representation learning, VLDB Endowment, 2019.\\\n[b] Cuturi and Doucet, Autoregressive Kernels for Time Series, arXiv, 2011.\\\n[c] Baydogan et al., Time series representation and similarity based on local autopatterns, Data Mining and Knowledge Discovery, 2016.\\\n[d] Mikalsen et al., Time series cluster kernel for learning similarities between multivariate time series with missing data, Pattern Recognition, 2017.\\\n[e] Paparrizos et al., Debunking four long-standing misconceptions of timeseries distance measures, ACM SIGMOD, 2020.\n"}, {"review_id": "yuVFNQA4xb", "reviewer": "Reviewer_8sG2", "summary": "The authors introduce a new time series kernel based on Nonlinear Vector AutoRegressive (NVAR) processes, following recent literature on its equivalence to reservoir dynamics. The kernel operates on time delay embeddings and enables the computation of similarity between time series with different lengths. The proposal is paired with SVM models and evaluated in classification tasks for both univariate and multivariate cases. The obtained accuracy and speed results are promising when compared to other kernel methods and a reservoir computing alternative.", "questions": "- I believe a more thorough discussion on how Taken's theorem holds with random lagged dimension dropping (Eq. (5)) is required. Perhaps some of the content of the supplementary material could be brought to the main text.\n\n- In practice, is there any interplay between the choice of $k$ lagged steps and polynomial order $n$? In which scenarios n > 2 would be useful? \n\n- How the proposed method behave for larger dimensional datasets? It would be interesting to include the dataset dimension in Table 2.\n\n- Lines 315-316: \"NVARk can, in fact, get around both issues with its non-recursive structure and the absence of a training phase.\" -> I suppose the ridge regression step and the paired SVM are not considered here, which should be made clearer.\n\n- Since the tuned NVARk* usually presents even better results, it would be interesting to verify its trade-off by including it in the experiments of Section 4.3 (Execution time and scalability).\n\n- It would be valuable to comment about the use of the introduced NVARk in the context of Gaussian processes models. Maybe a brief discussion on how it would compare with other GP approaches for sequences [1,2,3,4] could be included.\n\nReferences\n\n[1] Frigola-Alcade R. and Rasmussen, C. Integrated pre-processing for Bayesian nonlinear system identification with Gaussian processes. IEEE CDC, 2013.\n\n[2] Mattos, C. et al. Recurrent Gaussian processes. ICLR, 2016.\n\n[3] Al-Shedivat, M. et al. Learning scalable deep kernels with recurrent structure. JMLR, 2017.\n\n[4] Toth, C. and Oberhauser H. Bayesian learning from sequential data using Gaussian processes with signature covariances. ICML, 2020.", "limitations": "The manuscript states sufficiently its limitations.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "The manuscript is well written and organized. The motivation is clear and the tackled problem (building kernels for time series data) is relevant to the community. The discussion of the related methods is comprehensive, although I have missed a nod to Gaussian process models (see 'Questions' section). \n\nThe main proposal seems to merge the contributions by Bianchi et al. (2020) and Bollt (2021), evaluating the NVAR equivalence to RC in the task of time series classification. The simple to understand heuristics that form the NVARk is a relevant strength of the work.\n\nThe main proposal, NVARk, is explained in detail and the provided source code that implements Algorithm 1 is easy to follow. I also praise the provided supplementary material and the qualitative result for KPCA visualization.", "weaknesses": "The performed experiments are quite extensive, but the results presentation can be improved. For instance, Table 1 could also present the average rank of each method and be paired with boxplots of accuracy results across the 122 univariate datasets. Beside the per dataset results, Table 2 should also include the metrics considered in Table 1 (rank, average score), for the sake of homogeneity.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We warmly thank reviewer $\\color{green}{\\textbf{8sG2}}$ for the extensive review and positive feedback.\nWe particularly appreciate the positive remarks on clarity of the presentation and soundness of the work, which we are happy to see reflected in the scores.\nIn addition, we value your comment on the relevancy of kernels for time series analysis.\nAs for the raised concerns and questions, individual points are addressed below.\n\n\n**W1. Results presentation can be improved**\\\nWe agree, our initial presentation choices were mostly dictated by the space constraint. \nIn the new version of the manuscript, we now provide the average rank metric, which reflects what is observed in AVERAGE SCORE: \n\nAVERAGE RANK: NVARk = 2.1, SINK = 2.2, rmESN = 2.6, TCK = 2.8\n\nWe now also provide consistent metrics for univariate and multivariate data, i.e.\n\nAVERAGE SCORE (multivariate): NVARK* = 0.890; NVARk = 0.877; rmESN = 0.857; TCK = 0.842\\\nAVERAGE RANK (multivariate): NVARk = 1.62; rmESN = 2.12; TCK = 2.25\\\nN FIRST RANKED (multivariate): NVARk = 10; rmESN = 2; TCK = 4\n\nAs for the boxplot, we argue that it might confuse the reader as we do not claim a statistically significant difference with respect to SOTA.\n\n\n**Q1. More thorough discussion on how Taken’s theorem holds with random lagged dimension dropping is required**\\\nWe agree that the paper would benefit from a more thorough introduction to state space reconstruction for the unfamiliar reader.\nTakens' theorem establishes an upper limit on the necessary concatenated lag count that is necessary to form a valid equivalent of the generating dynamical system.\nIn Takens' formulation, there is no importance bias towards any specific lag and their number is the only relevant factor.\nIn essence, this rationale underpins our adoption of random lag dropping, as long as 'enough' concatenated dimensions remain.\nIn our next revised version, we will prioritize integrating the main text with material from Appendix A.\n\n\n**Q2. Is there any interplay between the choice of lagged steps and polynomial order?**\\\nThere is, but we could not find any use for it.\nThe polynomial order ($n$) regulates the degree of interactions between dimensions.\nWith $n=2$, NVAR captures mutual correlations between attributes and auto-correlations across lags.\nWith $n=3$, three-way interactions are considered as well, such as dimensions being correlated only if a third one exhibit certain patterns.\nIf these are important for the dataset at hand, then $n=3$ could, in principle, be considered.\nMeanwhile, $k$ influences the model's memory span.\nIn words, more lags would enable the model to explore these interactions further back in time.\nHowever, in practice, relevant three-way interactions tend to be rare and our random sampling often leads to inadvertent spurious correlations (something that would appear even more often with a high $k$).\nAs a result, we have never found any advantage in using $n > 2$.\n\n\n**Q3. How the proposed method behave for larger dimensional datasets?**\\\nWe kindly refer the reviewer to our common response \\#2.\n\n\n**Q4. Clarify \"absence of training phase\"**\\\nThe reviewer is right in saying that ridge regression and SVM are not considered in the mentioned statement.\nAt ll.315-316, we have substituted this terminology with a reminder of NVARk's linear complexity in the number of time series ($N$).\n\n\n**Q5. It would be interesting to verify trade-off of optimized NVARk**\\\nThe reported superior performance for NVARk* can be obtained with a median slowdown factor of $\\sim25$ (parallelized grid search + 1 iteration with best parameters) with respect to the heuristic-based NVARk. \nMost unfavorable cases are high-T datasets.\nAs for asymptotic behavior, this is trivially similar to NVARk scaled by the grid size.\nNVARk* also places favorably with respect to the considered baselines.\nWe report here the median computation times:\nNVARk = $0.5 s$, NVARk* = $13.7 s$, SINK = $24.2 s$, rmESN = $38.3 s$, TCK = $44.9 s$.\nInterestingly, we also observed that NVARk* execution time is faster than 1 iteration of rmESN for 107/122 univariate datasets. \n\n\n**Q6. It would be valuable to comment about the use in the context of GP**\\\nThis is an interesting direction that can be further explored. \nWe find these two directions appealing.\nAs a first example, such MTS kernels can, in general, be used in multi-task GP, where they can be paired with a scalar kernel over time (like RBF) for simultaneous modeling of multiple time series:\n$$ k(t_1, t_2, i, j) =  RBF(t_1, t_2) \\cdot K_{MTS}(i,j) $$\ne.g. spatio-temporal data or medical data for different patients.\nSecondly, and more specifically to NVARk, we recognize a compelling similarity between building our feature map (linear terms and polynomials) and different orders of information in signature transforms.\nAs in a reference suggested by the reviewer, signatures are being employed as a building block in deep kernels to build larger GP models, which might pave the path for interesting future directions.\n\nWe hope to have clarified all of the reviewer’s concerns and are happy to provide further details.\n"}, {"review_id": "81dboNGRml", "reviewer": "Reviewer_Nbgs", "summary": "This paper presents a kernel that can be applied to time series - both\nunivariate and multivariate - that draws inspiration from reservoir\ncomputing.  Practically, it constructs an expanded time series by\nsampling lags and then computing polynomial combinations of the\noriginal and lagged data through time.  Experiments show that the\nkernel, when combined with an SVM, is more accurate and more efficient\nthan a number of relevant competitors.\n", "questions": "(1) How sensitive is the approach to hyperparameter choice, both in\nterms of complexity and accuracy?\n\n(2) Are both lags and products important for accuracy?  Why or why not?\n", "limitations": "There is no such discussion in the paper.\n", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "This paper pulls on an interesting thread of work in reservoir\ncomputing (RC) that makes connections between RC and nonlinear vector\nautoregressive models (NVARs).  RC can be computationally expensive,\nand the NVAR approach simply combines time series with lagged copies\nand nonlinear functionals, like products.  The paper further connects\nthe use of lags to Takens' theorem about reproducing noise-free\nnonlinear dynamics through state space reconstruction using lagged\nobservations.  The discussion about the various connections to this\nwork is interesting.\n\nThe paper is clear and relatively easy to read and understand.  RC and\nNVAR are summarized in a useful way without getting bogged down in\ndetails, but with enough information to follow the presentation.  The\nfigures are also very helpful in understanding the paper.\n\nOne potential weakness of the approach is the relatively large number\nof hyperparameters, including the number and size of the lags, the\norder of the polynomial regression, and the embedding dimension.  The\npaper, though, does a good job presenting and arguing for heuristic\nchoices for these parameters that seem to perform well in practice.\n\nFinally, empirical results suggest that the proposed approach is as\naccurate as the state of the art and more efficient.\n", "weaknesses": "The primary contribution of the paper is a kernel that is accurate and\nfast to compute.  There are claims in the introduction where the\nauthors suggest that the kernel is more interpretable than competitors\nthough there is no discussion of that in the remainder of the paper.\n\nThe approach is conceptually simple, which is good - fit a model to an\naugmented time series (lags and products) and use the model parameters\nas a vector-based representation of the time series.  The\nhyperparameter choices - of which there are many - are made\nheuristically.  Some indication of the sensitivity of the approach to\nvariation in those choices would be good.  That is, are the heuristics\na result of just those values working empirically; or is the method\nhighly sensitive to those choices?\n\nThe complexity analysis would benefit from doing a direct comparison.\nThat section claims that the proposed methods does one expensive\ncomputation (fitting a model) per, compared to other methods that do\nan expensive computation per pair of instances.  What is the total\ncomplexity of the approaches?  If competing approaches have less\nexpensive per-pair cost they may still be more efficient overall.\n\nIt's unclear what the kPCA plots add to the paper.  What dataset is\nrepresented in the figure?  Why did you choose that dataset?  I'm sure\nthere are other datasets for which the separation for rmESN is better\nthan for NVARk.\n\nAn ablation study in which the MTS were augmented with just lags or\nwith just products would help tease out which element is most\nimportant.  Takens says nothing about nonlinear features.\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We warmly thank reviewer $\\color{magenta}{\\textbf{Nbgs}}$ for the extensive review and constructive feedback.\nAmong all, we really appreciate the kind words on the clarity of our paper and figures, as well as showing interest in how we connect our work to different research areas.\nAs for the raised concerns and questions, individual points are addressed below.\n\n\n**W1. There is no discussion of interpretability in the remainder of the paper**\\\nFor a discussion on our interpretability claims, we kindly refer the reviewer to our common response \\#1.\n\n\n**W2. How sensitive is the approach to hyperparameter choice?**\\\nIn addressing this, we split hyperparameters into classes, as we believe a separate discussion is needed for $k$ and $s$. \n\n-  \n   - $\\bar{d}_r$: as we show in Fig. 1 (top) (attached pdf) average performance for multivariate datasets are not considerably impacted in the range [75, 150]. As expected, we observe a significant drop of performance for high values (redundancy).\n   - $\\lambda_{ridge}$: we show in Fig. 1 (bottom-left) little sensitivity in the region $[5, 50]$. Results are averaged across a pool of 20 randomly sampled univariate datasets. The employed OCReP method outperforms the best choice by $\\sim 2\\%$, as it adapts the value to each individual dataset.\n   - $\\gamma_{rbf}$: in Fig. 1 (bottom-right), we range different multiplicative factors and show little sensitivity in a neighborhood of the unity.\n\n- The sensitivity to $k$ and $s$ strongly depends on the dataset, although it generally tends to be high.\nIn response, we introduce our dataset-specific heuristic adapting ideas from a precedent study [a], which we find works well in practice.\nA reasonable idea of the average sensitivity to these parameters can be obtained by comparing the performance of NVAR, where k and s are chosen by heuristic, to the performance of NVARk*, where k and s are optimized.\n    Unfortunately, understanding sensitivity by considering other settings is challenging, as the existing literature predominantly focuses on the noise-free case.\n    Furthermore, there is no strong evidence for the wider applicability of such approaches outside dynamical systems (see **Q1** from reviewer $\\color{red}{\\textbf{RJLc}}$).\n\n\n\n**W3. What is the total complexity of the approaches?**\\\nAs evidenced in Section 3.4, NVARk exhibits an overall complexity of $\\mathcal{O}(NT\\bar{d}^2_r)$. \nSINK and GAK report $\\mathcal{O}(T \\log(T))$ and $\\mathcal{O}(\\min(T_1, T_2))$ respectively, which do not counterbalance their $\\mathcal{O}(N^2)$ nature.\nA separate discussion is needed for TCK, which gets around its cubic complexity by putting an upper threshold on the length of the sampled time segments ($\\mathcal{O}(T_{max}^3)$).\nAs for dimensionality $\\bar{d}_r$, we treat it as a constant and do not discuss complexity in connection to this parameter.\n\n\n**W4. It’s unclear what the kPCA plots add to the paper**\\\nDespite the underlying equivalence between NVAR and simple RC, better performance can be obtained by finding a better parameter setting, which is not trivial for RC.\nIn our approach, a better setting is easier to find with fewer and more interpretable parameters. \nWe mean to use the kPCA example (*CinCECGTorso* dataset) to visually highlight the possible extent of the differences between NVAR and rmESN, i.e. theoretically equivalent methods under different working points.\nAs a side note, it highlights that a kernel can be used for visualization.\n\n\n\n**W5. Are both lags and products important for accuracy?**\\\nThis is an interesting ablation study that has not been investigated in our work and previous works to the best of our knowledge.\nWe provide results for such additional experiments in Tab. 1 (attached pdf).\nWe separately consider the case of fitting directly on the input (no concatenation), sampling linear terms only, and sampling non-linear terms only. \nAs a general trend, we observe that using all terms leads to the best results, followed by \\textit{non-linear} and \\textit{linear}.\nNote also that all configurations perform better than fitting directly on the input series.\nIn the univariate case, the \\textit{linear} variant tends to underperform, as the total number of concatenated dimensions is very small (equal to $k$).\nFor multivariate datasets, the two variants perform very similarly, as all lags of attributes may already include a considerable number of dimensions.\nBoth results are also in line with our interpretation using the generalized Taken's theorem. \nInspired by this study, we have also tested an additional variant in which we add all linear terms, and then fill remaining dimensions (up to $\\bar{d}_r$) by sampling nonlinear ones.\nInterestingly, this leads to even better results and a reduction in variance.\nWe believe this compensates for the imbalance between all possible linear and non-linear terms, while also reducing potential noise from uncorrelated dimensions.\nWe sincerely thank the reviewer for raising this question, which led to additional interesting results and the integration of this comprehensive analysis into the main experimental section.\n\n**W6. There is no discussion of limitations **\\\nWe kindly refer the reviewer to lines 330-331, where the main limitation is mentioned in the manuscript.\nPlease also find a further discussion in our common response \\#2.\n\nWe hope to have clarified all of the reviewer’s concerns and are happy to provide further details.\n\n**References:**\\\n[a] Paparrizos et al., State space reconstruction parameters in the analysis of chaotic time series - the role of the time window length, Physica D, 1996."}, {"review_id": "LBqI7lUblZ", "reviewer": "Reviewer_XEWY", "summary": "This focus of this paper on a new type of kernel for reservoir computing (using nonlinear vector autoregressive delay embeddings) that is faster and more accurate than related kernels. \n\nThe three main contributions\n* Introduction of a new kernel for time-series modeling\n* State of the art results on uni- and multimodal problems on real-world datasets (related work has used synthetic datasets.)  \n* Theoretical connections to Taken's theorem and the field of state space reconstruction. \n\nThere are only a small number of (intuitive) hyper parameters for tuning the model. While this approach does not seem suitable for large datasets (where deep learning-based solutions reign), it appears to be a nice and efficient solution for smaller time-series problems. ", "questions": "At the end of section 2, there is discussion on the lack of stability of techniques for reservoir computing. Perhaps I'm being naive, but this seems obvious given that many of the parameters (e.g., W_in and A) are randomly initialized and not optimized over -- this means that parameters that are updated (e.g., W_out and c) need to compensate for potentially poor choices of W_in and A. I know this isn't something you introduced, but generally speaking why is this a reasonable idea?\n\nThe approach is described as \"a simple concatenation of the input series with time-delayed copies and nonlinear functionals.\" Could this approach simply be implemented using convolutions (with dilations)? Would there be a more straightforward way of writing out this math in this case?", "limitations": "This paper is more theoretical and does not touch on noteworthy societal impacts. ", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Overall the paper is well written. I don't have a deep understanding of reservoir computing but I think they authors do a nice job of connecting different technical areas. If I was working on problems in this area I would be strongly interested in building a deeper understanding of this paper. \n\nThe authors do a nice job of contextualizing the work, especially given that this is deviates from many recent deep learning oriented time-series papers. \n\nThe core idea behind the kernel is elegant and clearly has a large impact on efficiency of the approach.  \n\nI appreciate the discussions on asymptotic complexity and the qualitative notes and visualizations. ", "weaknesses": "It looks like the authors only use 122 UCR unimodal datasets but my understanding is that the most recent UCR benchmark (since 2018) has been using 128 unimodal datasets. Is there a good reason for this discrepancy? Similarly, the UEA corpus should have 30 datasets whereas only 22 are used in this paper. \n\nMinor: it's challenging to interpret the graphs (especially Fig 2) when printed in black and white. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We warmly thank reviewer $\\color{orange}{\\textbf{XEWY}}$ for the thorough review and positive feedback.\nWe really appreciate the favorable comments on clarity of the writing and elegancy of the underlying idea, as well as positive remarks on how we bring together different research areas.\nAs for the raised concerns and questions, individual points are addressed below.\n\n\n**W1. Is there a reason for benchmarking on a subset of datasets?**\\\nAmong the univariate datasets, we have excluded 8 datasets with a high $N_{classes} / N_{train}$, which are then incompatible with our choice of stratified cross-validation for the SVM.\nResults for those datasets can still be obtained by considering a shuffle split rather than stratified, which we report in Tab. 2 (attached pdf).\nAmong multivariate ones, we have excluded 6 high-dimensional datasets for which the application of the NVAR kernel is not appropriate, and 5 datasets for which SVM accuracy was $\\ll 50 \\%$ for all approaches, which indicates that kernel methods are, in general, not a suitable solution.\nThanks to the reviewer's comment, we now include 3 additional multivariate datasets that we did not originally consider.\nCorresponding results are shown in Tab. 2.\nAverage accuracy for the additional univariate and multivariate datasets indicates superior performance of NVARk with respect to the baselines.\nAll results have been integrated into the experimental section of the manuscript.\n\n\n**W2. Figures in black and white**\\\nWe thank the reviewer for pointing this out, we have included a line where, for interpretation of colors, we refer the reader to the web version of the article.\n\n\n**Q1. Why are unoptimized parameters in RC a reasonable idea?**\\\nOne of the greatest advantages of keeping the parameters of the recurrent update unoptimized is to avoid exploding or shrinking the gradient when backpropagating through time.\nThis is a complementary solution to maintaining the memory unaltered with the popular gating mechanisms, e.g. LSTMs and GRUs, but much more efficient.\nFurthermore, despite their randomness, the reservoir dynamics have the potential to amplify the relevant features in the input data, actually reducing the burden on the linear readout.\nHowever, if the choice of the hyperparameters (controlling such dynamics) is poor, an adverse effect is obtained.\n\n\n**Q2. Could the approach be rewritten using dilated convolutions?**\\\nIt is interesting to reason about this.\nConvolutions can be declined into products between the input and a sliding filter or between two delayed signals.\nRegarding the first type, in our approach, there is no weighting of the input dimensions until the final readout. \nDespite this layer certainly involving some terms of the form $w_k x_{t-j}$, non-linear terms are different.\nFor what concerns these non-linear terms, they are products between the input and its lags.\nIt is then definitely possible to think of a rewriting in terms of causal convolutions of the input with itself, despite the advantages not being fully clear to us.\nFor a full rewriting, one can think of merging both considerations. \nHowever, at first sight, this seems far from straightforward.\n\nWe hope to have clarified all of the reviewer’s concerns and are happy to provide further details.\n"}]

Author rebuttal:
We would like to remark here our thanks to all reviewers for their extensive reviews.
We are happy to hear that reviewers acknowledged novelty ($\color{red}{\textbf{RJLc}}$, $\color{blue}{\textbf{hHQP}}$) and relevancy ($\color{blue}{\textbf{hHQP}}$, $\color{green}{\textbf{8sG2}}$) of our work, as well as praising its interdisciplinarity ($\color{red}{\textbf{RJLc}}$, $\color{magenta}{\textbf{Nbgs}}$, $\color{orange}{\textbf{XEWY}}$).
We are also grateful to all 5 reviewers for all sharing their appreciation for the clarity of the writing.

In order to strengthen our paper, incentivized by the reviewers' feedback, we have brought sever minor modifications together with the following significant additions:

- we now compare computational time of NVARk* with respect to NVARk and baseline methods, highlighting performance trade-off ($\color{green}{\textbf{8sG2}}$);
- we now provide an additional study of the sensitivity of our approach to the different hyperparameters, either by providing additional experiments (Fig. 1) or extensive discussion ($\color{magenta}{\textbf{Nbgs}}$); 
- we now present a novel ablation study (Tab. 1), deepening the importance of linear and non-linear components of NVARk ($\color{magenta}{\textbf{Nbgs}}$);
- we now extend our evaluation to 8 more univariate and 3 more multivariate datasets (Tab. 2) ($\color{orange}{\textbf{XEWY}}$).

Please find below responses which we hope help to clarify shared concerns across more reviewers.
Finally, please find attached a one-page pdf providing the additional experiments that are referenced across our responses.
We look forward to the reviewers' individual responses to our rebuttal.

**1. Discussion of the interpretability claims**\
Our claims refer to NVAR's ability to circumvent the inherent randomness of RC methods.
As in ll.184-186, NVARk representation vectors encapsulate the mutual and auto-mutual information within the dimensions of the input.
In contrast, representations from RC methods encompass hardly-interpretable reservoir dynamics. 
Additionally, the key parameters of NVARk are easy to interpret: $s$ controls the spacing between lags and just has to be set high enough not to sample noise; $k$ controls the number of lags and extends memory to find the relevant patterns.
On the other side, for instance, input scaling and spectral radius are much less interpretable.
These are two notable RC parameters controlling, among other aspects, the degree of nonlinearity injected into the reservoir.
How much of this is required by the task is not easy to judge and requires experienced insight into nonlinear dynamics or, most of the time, traditional hyperparameter tuning.
In our view, the absence of heuristics to choose these parameters also supports this.
We have now made this connection between our claims and the above considerations more explicit in Sec. 3.2 of our updated paper.


**2. Limitations and behavior with high-dimensional datasets**\
As the input dimensionality approaches $\bar{d}_r$ (our limit), there is little margin for adding more dimensions.
Here, concatenation of lags and non-linear terms would require substantially increasing $\bar{d}_r$.
This places excessive strain on the linear readout, which inherently possesses limited expressive capacity.
This component struggles to effectively model intricate relationships among numerous features and is increasingly overloaded with spurious correlations.
We then expect the effectiveness of NVARk to decrease.
To overcome this, one can consider replacing random sampling with strategies that prioritize the selection of meaningful terms.
In fact, we believe this is a promising future direction for our work.
Alternatively, it would be interesting to explore how different linear layers, e.g. Lasso, would perform in this regime.
We acknowledge the significance of this aspect in our work, which may not have received sufficient attention in the main text. 
We now integrate the provided discussion into the Conclusions of the manuscript.


Meta-review: Five reviewers were assigned to this paper. Four reviewers support acceptance, and a fifth reviewer, Nbgs, supports rejection. Reviewer Nbgs, who supports rejection, did not reply to the authors' rebuttal or get involved in the post-rebuttal discussion. I read Reviewer Nbgs comments, and I did not find any of these comments as a cause for rejection. The comments are mostly in the sense of requesting further clarifications, which the authors have provided in the rebuttal. Therefore, I accept the paper.
Decision: Accept (poster)