Year: 2021

Abstract:
Time series analysis is a widespread task in Natural Sciences, Social Sciences, and Engineering. A fundamental problem is finding an expressive yet efficient-tocompute representation of the input time series to use as a starting point to perform arbitrary downstream tasks. In this paper, we build upon recent works that use the Signature of a path as a feature map and investigate a computationally efficient technique to approximate these features based on linear random projections. We present several theoretical results to justify our approach and empirically validate that our random projections can effectively retrieve the underlying Signature of a path. We show the surprising performance of the proposed random features on several tasks, including (1) mapping the controls of stochastic differential equations to the corresponding solutions and (2) using the Randomized Signatures as time series representation for classification tasks.When compared to corresponding truncated Signature approaches, our Randomizes Signatures are more computationally efficient in high dimensions, and often lead to better accuracy and faster training. Besides providing a new tool to extract Signatures and further validating the high level of expressiveness of such features, we believe our results provide interesting conceptual links between several existing research areas, suggesting new intriguing directions for future investigations.

heading: 1 INTRODUCTION
text: Modeling time series is a common task in finance, physics, and engineering. A frequent challenge is finding a transformation mapping a newly observed time series into a target one (seq2seq modeling) or into a label summarizing its salient properties (classification). In the absence of any principled model describing such a mapping, one has to infer it from data. The last few years have witnessed the rise of deep neural networks, which have found successful application to problems involving time series in numerous domains (Fawaz et al., 2019; Gamboa, 2017). Nevertheless, their outstanding performance comes at the price of over-parametrization, data hungriness, and expensive training cost (Werbos, 1990; Brown et al., 2020; Teubert et al., 2021; Neyshabur et al., 2018; Marcus, 2018). Furthermore, even if sufficient data is available, the resulting models learn representations of the input data that are highly specialized to the training task and difficult to adapt in different contexts. In addition, the remarkable performance of these methods is often the result of a substantial engineering effort and is not supported by theoretical results.
Reservoir Computing (RC) (Schrauwen et al., 2007) offers an intriguing alternative strategy to cope with the limitations above, yet retaining the universal approximation properties typical of deep neural networks. In RC, the learning is divided into two phases: first, data is passed through an untrained reservoir which extracts a set of task-independent features; second, a simple and efficient-to-train linear map (the readout map) projects such features into the desired output. The critical point is that the design of the reservoir determines the expressiveness of the features, and several alternatives can be found in the literature (see (Gauthier et al., 2021) and references therein).
A powerful reservoir is offered by the Signature Transform, often simply referred to as Signature, stemming from rough path theory (Ben Hambly, 2010; Friz and Hairer, 2020). The Signature of a path consists in enhancing the path with additional curves, which, in the smooth case, corresponds to iterated integrals of the curve with itself. A profound mathematical result (Levin et al., 2013) guarantees that the solution of a (rough) differential equation can be approximated arbitrarily well
by a linear map of the Signature of the input signals (a.k.a. controls). In the machine learning jargon, the Signature of a path can be interpreted as a feature map extracting all its geometrical properties and thus allowing a simple linear map to approximate any function of it. On the other hand, it is often the case that the reservoir features are very high dimensional, and hence are particularly expensive to calculate and use in downstream tasks. In addition, the high-dimensionality of the Signature reservoir poses additional challenges for modern gradient-based optimizers as convergence rates suffer from a linear dependence in the model dimension (Bottou et al., 2018).
Inspired by the remarkable theoretical properties of the Signature reservoir and motivated to fix its practical pitfalls, our contribution is to showcase the effectiveness of Randomized Signatures (Cuchiero et al., 2021b;a), a recently introduced reservoir of random features. These features provably hold the same geometric properties and approximation power as the Signature, yet are often more efficient to compute and of lower dimensionality. To extract Randomized Signatures, we numerically integrate a set of random linear stochastic differential equations driven by the path itself. We assess the expressiveness of the proposed random features on several tasks, including nonparametric (black-box) system identification problems arising from complex controlled nonlinear dynamical systems as well as time series classification.
The paper is structured as follows: In Section 2, we link the Randomized Signature approach with related research areas in the literature; in Section 3, we present relevant theoretical results motivating our method; in Section 4, we assess the effectiveness of Randomized Signatures on various system identification and time series classification tasks, showing promising results when compared to other signature approaches for reservoir computing. Finally, in Section 5, we discuss our conclusions and possible future developments.

heading: 2 RELATED WORKS
text: Random Features and Reservoir Computing. The idea of extracting features based on random operations is not new and has seen a number of successful applications over the past years. Of particular note, the seminal work of Rahimi and Recht (2008) proposes to accelerate kernel machines by learning random features whose inner product matches that of a target shift-invariant kernel. The trade-off between generalization and computational efficiency of learning with random features has then been rigorously studied by Rudi and Rosasco (2017). A conceptually very similar rationale is introduced by a parallel series of works exploring the topic of Reservoir Computing (Schrauwen et al., 2007). Similarly to us, Echo State Networks (Jaeger, 2003) evolve the input state by a series of fixed random projection (the reservoir) and generate the output by applying a trainable linear projection over the hidden states. However, we make the additional step of linking the random features to the Signature of the input path and, as shown below, the evolution of the features is dictated by a randomly-evolved stochastic differential equation driven by the input path.
Controlled Differential Equations. Our work is also related with a series of recent papers investigating the problems of how to process irregular time series and to condition a model on incoming information through the lens of controlled differential equations (CDEs) (Kidger et al., 2020; Morrill et al., 2021). Interestingly, Kidger et al. (2020) show that the action of a linear layer on the final
output of a trained Neural CDEs – the extension of Neural ODEs (Chen et al., 2018) to CDEs – results in a universal function approximator. Differently from them, our method is less computationally expensive as the only parameters we need to train are those of the final linear readout layer. Fundamentally different to our work, the approach of Morrill et al. (2020) first randomly projects the possibly high dimensional controls into a lower dimensional space and then extracts the Truncated Signature from such compressed input. Our method instead extracts a random compression of the signature of the original input controls. Rough Path Theory. Rough path theory is about describing how possibly highly oscillatory (rough) control path interact with nonlinear systems (Lyons, 2014). The concept of Signature is introduced in this context to provide a powerful description of the input path, removing the need to look at the fine structure of the path itself. Recent years have seen a resurgence of these ideas, which have been revisited from a machine learning perspective (Bonnier et al., 2019; Kidger and Lyons, 2021). Our analysis is strongly influenced by the work of Cuchiero et al. (2021a;b), who, starting from the observation that the Signature is an infinite dimensional reservoir, establishes that its information content can be efficiently compressed by a random projection performed by a dynamical system driven by random vector fields. As we show later, this results in a random reservoir which preserves the properties of the Signature while living in a finite dimensional space.

heading: 3 BACKGROUND
text: We provide the theoretical tools supporting the Randomized Signature approach for feature extraction in times series. A formal discussion can be found in (Friz and Hairer, 2020) and (Cuchiero et al., 2021b).

heading: 3.1 RANDOMIZED SIGNATURE OF A PATH
text: Let X : [0, T ] → Rd be a continuous piecewise smooth d-dimensional path X = ( X1, · · · , Xd ) . We will refer to X as the control and to its single components Xi as controls. We denote by {e1, . . . , ed} the canonical basis of Rd.
Definition 1 (Signature) For any t ∈ [0, T ], the Signature of a continuous piecewise smooth path X : [0, T ]→ Rd on [0, t] is the countable collection St := ( 1, S1t , S 2 t , . . . ) ∈ ∏∞ k=0 ( Rd )⊗k
where, for each k ≥ 1, the entries Skt are the iterated integrals defined as
Skt := ∑
(i1,...,ik)∈{1,...,d}k
(∫ 0≤s1≤···≤sk≤t dXi1s1 . . . dX ik sk ) ei1 ⊗ · · · ⊗ eik . (1)
We define the Truncated Signature of X of order M ≥ 0 as
SMt := ( 1, S1t , . . . , S M t ) ∈ M∏ k=0 ( Rd )⊗k =: T M ( Rd ) . (2)
A practical example of computation of Signatures is presented in Appendix A.3.
The definition in the last paragraph suggests that the Signatures of X can be used to approximate any regular enough function of X , for instance the solution to differential equations controlled by X . The following result makes this argument precise in the multidimensional setting.
Theorem 1 (Signature is a Reservoir) Let Vi : Rm → Rm, i = 1, . . . , d be vector fields regular enough s.t. dYt = ∑d i=1 V i (Yt) dX i t , Y0 = y ∈ Rm, admits a unique solution Yt : [0, T ] → Rm. Then, for any smooth test function F : Rm → R and for every M ≥ 0 there is a time-homogeneous linear operator L : T M ( Rd ) → R which depends only on (V1, . . . , Vd, F,M, y) s.t.
F (Yt) = L ( SMt ) +O ( tM+1 ) , t ∈ [0, T ]. (3)
This theorem suggests the first M entries of the Signature of X are sufficient to linearly explain the solution of any differential equation driven by it. In addition to this, it supports the claim that such features are valuable for any downstream usage.

heading: 3.2 RANDOMIZED SIGNATURE OF A PATH AND ITS THEORETICAL GUARANTEES
text: Calculating SMt requires the calculation of dM+1−1 d−1 iterated integrals (Appendix A.3) – which in total quickly becomes computationally expensive. Several computational techniques have been developed to circumvent this problem, see, e.g. (Kidger and Lyons, 2021). The next result provides a practical description of how it is possible to reduce the computational burden without losing too much explanatory power. As such, this results provides the theoretical foundation for our approach.
Theorem 2 (Randomized Signature (Informal)) For any k ∈ N big enough and appropriately chosen random matrices A1, . . . , Ad in Rk×k, random shifts b1, . . . , bd in Rk×1, random starting point z in Rk×1, and any fixed activation function σ, the solution of
dZt = d∑ i=1 σ (AiZt + bi) dX i t , Z0 = z ∈ Rk, t ∈ [0, T ]. (4)
– called the Randomized Signature of X – has comparable approximation power as the Signature itself and maintains its geometric properties.
The formal statement of such a result is the combination of Theorem 8, Theorem 9, and Definition 10 (see Appendix). A complete formal discussion as well as the original statement of this theorem can be found in the paper (Cuchiero et al., 2021a). In a nutshell, the Randomized Signature Z follows a dynamics which provides an efficient and powerful compression of the Signature through a low-dimensional random projection. The expressiveness of such projection is guaranteed by a standard Johnson-Lindenstrauss argument (Cuchiero et al., 2021a).
Computational complexity and dimensionality of Randomized Signatures. The computational complexity of calculating Zt is O(k2d), while its dimensionality is O(k). However, we should then ask the crucial question “how should we compare k with M?”. We show experimentally in Section 4 that, in order to match the approximation guarantees of the Truncated Signature of order M , the number of required Randomized Signatures k is not too big — in particular it is not exponential in M . This confirms that working with Randomized Signatures is often less computationally demanding and results in lower-dimensional features.
We conclude the section remarking the fascinating implications of choosing a random Ais and bis.
Theorem 3 (Density of Randomized Signatures (Informal)) For any sequences of time points 0 ≤ t0 < . . . < tN ≤ T and points zt0 , . . . , ztN ∈ Rk we can find a smooth control X such that the solution Z of Equation 4 at time ti is such that Zti = zti , for i = 0, . . . , N . Additionally, if the control X is a d-dimensional Brownian motion, then the solution of Equation 4 at any point in time t > 0 admits a smooth density with respect to Lebesgue measure.

heading: 4 THE EXPRESSIVE POWER OF RANDOMIZED SIGNATURES
text: To compute the Randomized Signature, we rely on an Euler-Maruyama approximation; details are given in Algorithm 1. To start, we experimentally validate some of the results presented above. In
Algorithm 1 Generate Randomized Signature
Require: X ∈ Rd sampled at 0 = t0 < · · · < tN = T , Randomized Signature dimension k, activation function σ. Initialize: Z0 ∈ Rk, Ai ∈ Rk×k, bi ∈ Rk to have iid standard normal entries for n in {1..N} do Ztn = Ztn−1 + ∑d i=1 σ ( AiZtn−1 + bi ) ( Xitn −X i tn−1
) end for
particular, we show that the Truncated Signature of a control lies in the linear span of its Randomized Signatures, thus validating Theorem 2. As a second proof of concept, we reconstruct a control X as a linear combination of its Randomized Signatures, showcasing the expressiveness of such random features. Finally, we combine Theorem 1 and Theorem 2 and we use Z to learn the dynamics of complex nonlinear systems of differential equations and to perform time series classification. Our result show that the Randomized Signature is a powerful random reservoir which guarantees — both in theory and practice — the same geometrical properties and approximation power as the Signature itself, while often being of lower dimensionality and less computationally expensive both at extraction and at deployment time.
Note on the choice of σ. While the choice of σ does not affect the theoretical results (Cuchiero et al., 2021a;b), we found that a careful choice is needed for optimal results. Inspired by seminal works on the stability of deep linear networks (Glorot and Bengio, 2010) and by the connection to neural ODEs (Chen et al., 2018), we selected σ to be linear, with slope 1
d× √ k . This choice, complemented with additional fine-tuning, leads to the best performance for downstream tasks.

heading: 4.1 LEARNING THE SIGNATURES
text: In this experiment, we show that we can approximate the Truncated Signature of a control using a linear combination of its Randomized Signatures. We fix an equally-spaced partition D = {t0, · · · , tN} of [0, 1] and calculate the truncated Signature of order M of a d-dimensional path X on D, i.e. SMti , for i = 0, · · · , N . For a fixed value of k, we use Algorithm 1 to generate a realization of the Randomized Signature Z of X on D. For sake of clarity, the truncated Signatures on D are reshaped into an ( N × (( d(M+1) − 1 ) /(d− 1)− 1 )) -dimensional matrix called StrueM while the Randomized Signatures on D is represented as an (N × k)-dimensional matrix called ZRandk . Finally, we train a Ridge Regression to find a parameter matrix W ∈ Rk×((d (M+1)−1)/(d−1)−1) mapping the Randomized Signatures ZRandk into the truncated Signatures S true M .
More precisely, we consider a 10-dimensional control X = (t,Wt) where W is a 9-dimensional Brownian motion with independent components, and fix M = 6. For each k ∈ {1, · · · , 200}, we generate ZRandk , perform a simple linear regression and calculate the L
2 approximation error. Results in Figure 3 indicate that the approximation error decreases as k increases. This confirms that, as ensured by Theorem 2, the random features constructed as prescribed by Equation 4 hold the same information content as the Signature and thus maintain its approximation power and geometric properties. Finally, we highlight that, if we were to be satisfied with an approximation error of order 10−4, we could achieve it with roughly with k = 190 Randomized Signatures. To conclude, this means that, instead of calculating (( d(M+1) − 1 ) /(d− 1)− 1 ) = 1111110 integrals per time
step, we could just perform k2d = 36100 calculations per time step — which is 3 times cheaper. More remarkably, this high feature dimensionality (1111110 for Signature) makes any model for downstream tasks exponentially more complex than the that based on Randomized Signature (only 190 features). Modeling choices and details are reported in Appendix B.2.

heading: 4.2 RECONSTRUCTION AND RANDOMIZED SIGNATURE AUTOENCODERS
text: The Randomized Signature of X is rich enough to reconstruct it almost perfectly. In practice, we consider an “autoencoder” architecture, whose encoder is our random feature generator and the decoder is a trained linear projection mapping the Randomized Signatures into the original path. For this experiment, we fix an equally-spaced partition D = {t0, · · · , tN} and consider a single
0 50 100 150 200 k
4
2
0
2
4
6
8
10
Va lu
e
L2 Log-Error
Figure 3: Randomized Signature approximates the Signature as k increases. 0.0 0.2 0.4 0.6 0.8 1.0 Time
0.6
0.8
1.0
1.2
1.4
1.6
1.8
2.0
Va lu
e
Out of Sample
Extracted Xt True Xt
Xt X̂tZt
Figure 4: Performance of the Randomized Signature autoencoder.
realization XTrain of an Ornstein–Uhlenbeck process. Ultimately, we train a Ridge Regressor to map the Randomized Signature Z of the one and only train sampleXTrain intoXTrain itself. Table 6 in Appendix shows the average reconstruction error on NTest = 10000 test samples as we repeat the experiment for different choices of k. Such results are in accordance with Figure 4 which shows an example of the true and generated time series of a test sample. This confirms that the Randomized Signature of a path holds an equivalent information content as the path itself while being more usable for any downstream tasks. Modeling choices and details are reported in Appendix B.2.

heading: 4.3 LEARNING SDE SOLUTIONS WITH A LINEAR MAP ON THE RANDOMIZED SIGNATURE
text: In this series of experiments, we provide evidence of the potential of Randomized Signatures in several downstream tasks. First of all, as shown in Appendix B, we prove that our approach can accurately learn the dynamics of several rough differential equations: we highlight that the algorithm has only access to the input controls and the output response of the system but it is completely agnostic with respect to the law of the dynamics itself. To start with, we focus on the nonlinear dynamics of 1-dimensional stochastic Double Well. Secondly, we study the multidimensional case of the 4-dimensional Ornstein–Uhlenbeck. Then, we investigate the performance on the extremely rough dynamics of Fractional Geometric Brownian Motion even for low values of the Hurst parameter H . Additionally, we obtain promising results in learning the dynamics of an Enzyme-Substrate Reactions (Ingalls, 2013). We also demonstrate that this approach has sustainable performance even when the input control is irregularly sampled as well as Neutral Controlled Differential Equations. Due to space constraints, these results are presented in the Appendix. Here, we show the performance in learning highly nonlinear dynamical systems with a real world application.
Tumor Growth Kinetics. In this experiment, we consider predicting the output of a complex controlled nonlinear system using Randomized Signatures. The Simeoni Tumor model (Simeoni et al., 2004) is widely used to predict the evolution of tumor weight Yt under the concentration of a treatment drug Xt. The system is specified by the following equations
du1t = { λ0u 1 t [ 1 + (λ0/λ1wt) Ψ ]−1/Ψ − k2Xtu1t } dt, du2t = [ k2Xtu 1 t − k1u2t ] dt,
du3t = [ k1 ( u2t − u3t )] dt, du4t = [ k1 ( u3t − u4t )] dt, Yt = u 1 t + u 2 t + u 3 t + u 4 t .
(5)
We fix an equally spaced partitionD of [0, 10] and perform Ridge Regression to map 10000 instances of the Randomized Signature of the controls into the respective integrated solution Yt. We sampleXt from the law of a Scaled Squared Brownian motion 0.25·W 2t whereWt is a 1-dimensional Brownian Motion (BM) and the left side of Figure 5 shows the comparison of the true and the generated time series on a test sample along with the respective path of the control. To test the robustness of our approach, we use controls that do not follow the law used in training. More precisely, the middle of Fig. 5 shows what happens if we use 1{W 2t >1} as input control instead of 0.25 ·W 2 t while its right does the same as we use Wt. These two experiments prove that the algorithm learnt much more than merely predicting the behaviour of the system under ordinary circumstances. The first shows that the algorithm has “understood” that the absence of drug leads to the growth of the tumor while it dies back once the medicine is injected. On the other hand, the second suggests that the algorithm cannot behave properly when the input becomes negative — not allowed by the problem definition (negative drug concentration). Further details are reported in Appendix B.2. Electrochemical Battery Model. In this experiment, we learn the dynamics of the electrochemical battery model proposed in (Daigle and Kulkarni, 2013) which returns the voltage Y as the current
X is drawn from the battery. We use the open source NASA Prognostic Model Package (Teubert et al., 2021) to simulate voltage trajectories given input current control paths. On a fixed equally spaced partition D of [0, 500], we model the input current with step functions taking values 0 or 1 on random sub-interval of [0, 500]. To test the robustness of our approach, we train a Ridge Regression to map NTrain = 1000 instances of k-dimensional Randomized Signature of the controls into the respective solutions to which we add white noise. The left of Figure 6 shows the comparison between the ground truth and our prediction on an in-sample trajectory. In the same figure, we compare the generated time series for several choices of k on a test sample along with the path of the input control. In Figure 6 we show the performance of our method in comparison with a nonlinear neural ARX (NNARX) — which is often used in the modern control theory/system identification literature as a benchmarch (see e.g. (Masti and Bemporad, 2018)), given its strong theoretical guarantees (Schoukens and Ljung, 2019). Table 8 in Appendix reports the performances in terms of MSE on 1000 test samples for our models compared to two versions of NNARX: NNARX 22 which has around 1000 trainable parameters just like our model with k = 1000, and NNARX 1000 which matches our best model in terms of MSE but has around 106 trainable parameters. Details are reported in Appendix B.2.
Comparison with Neural Controlled Differential Equations. In Appendix B, we benchmark the performance of our model with that of NCDEs (Kidger et al., 2020), a recently proposed competitive deep learning technique for controlled dynamical system modelling. We test the two models on the task of learning the dynamics of an Ornstein-Uhlenbeck Process, and showcase the superior performance of our method both in terms of computational complexity and accuracy.

heading: 4.4 COMPARISON WITH TRUNCATED SIGNATURES
text: In this experiment, we learn the dynamics of the Fractional Ornstein–Uhlenbeck process using our Randomized Signature model and benchmark our results against a Truncated Signature model of order M . The computational cost of extracting the Truncated Signature of order M = 3 for d controls is O(d3) while for Randomized Signature it is O(k2d). To keep the computational cost of extracting features equal across models, we select k = d. As a result, the number of features of the Signature-based model is O(d3) while that of our model is O(d). In the spirit of Reservoir Computing, we choose both models to be linear. Hence, the number of trainable parameters is equal to the number of features. This implies that, the number of trainable parameters of our model is orders of magnitude smaller than that of our benchmark (left of Figure 7) — this severely affects the training time (see e.g. (Bottou et al., 2018)). Figure 7 shows that our training time is almost independent on d. Instead, the training complexity of the Signature-model is exponential in d and
the Out-of-Sample performance degenerates as the underlying optimization problem explodes in dimension. This result clearly highlights the data-hungriness of Signature-based models when the number of dimensions is high. To conclude, we emphasize that if d = 80, the storage cost of the Truncated Signature of a path is 200MB while its Randomized Signature is just 32KB. As a consequence, training the former onNTrain = 10000 would take an unpractical amount of time and storage memory while we were able to train ours in 4.8896 seconds using just 320MB and attaining an Out-of-Sample Error on NTest = 10000 equal to 0.075329. Modeling choices and details are reported in Appendix B.2

heading: 4.5 TIME SERIES CLASSIFICATION
text: In time series classification (TSC) the goal is to map an input time series to the corresponding label. Given the complex high-dimensional nature of temporal signals, the vast majority of classification algorithms aims at effectively finding expressive representations of the original time series and then maps them into the final prediction. This last operation is often realized via a relatively simple mapping, e.g. a linear map or a shallow neural network. Deep learning algorithms used for TSC include various types of deep convolutional neural networks which extract a hierarchical representation of the input time series by learning the parameters of their filters via backpropgation. These features are finally mapped into the logits associated with each class via a trainable linear projection. Last but not least, it is key to remark that these networks are typically over-parametrized and their training can be computationally demanding as well as challenging from an engineering perspective. In this section, we investigate how we can actually replace deep neural networks used for feature extraction with relatively cheap random layers which iteratively generate a reservoir. Among other notable works, Rocket (Dempster et al., 2020) arguably proposes the most successful solution.
Rocket. This algorithm prescribes the generation of k 1 unidimensional filters characterized by randomly drawn lengths, weights, and jumps. Given a filter of length lfilter, weights ω ∈ Rlfilter , and jump j, it can be easily shown that the feature Zi resulting from its action on an input path X , from position i in X , evolves according to the following iterative update:
Zik = Z i k−1 +Xi+(k×j) × ωk, k ∈ {0, · · · , lfilter − 1} (6)
where Zi0 = b and b is a random real number. We define Z i to be equal to Zilfilter−1. Depending on the values of j and lfilter and on the length lX of the path X , this update produces a vector Z = (Z1, · · · , ZM ) of M := lX − (l − 1) · j random features, used to extract only two nonlinear features: max(Z) and 1M ∑M i=1 1{Zi>0}. Ultimately, this means that if we use k kernels, we extract 2k features per time series which are then used as input to a Ridge classifier to obtain the final prediction. Rocket results in state-of-the-art performances in TSC both in terms of accuracy and computational efficiency, outperforming several deep learning baselines which often involve order of magnitude more parameters to optimize.
An alternative interpretation. Looking at Eq. 6, it is natural to wonder if it can be replaced with an alternative and possibly more principled update scheme driven by the input path and evolved by random projections. With the intent of aligning Rocket to the framework presented in this paper, we replace Eq. 6 with the update equation in Algorithm 1, thus exploiting the expressive power of Randomized Signatures. As detailed in the Appendix B.2, the current implementation of this
algorithm is equivalent to assuming that the random vector field A is approximately diagonal and autonomous. In other words, the Randomized Signatures are never mixed during their evolution and, differently from the update in Eq. 6, the filter weights do not depend on time. These represent two strong limitations which pave the way to several enhancements in future work. Despite these drastic shortcomings, the resulting algorithm outperforms several deep learning baselines on 128 TSC datasets from the UCR archive, as shown in Fig. 8. The original algorithm using Eq. 6 provides better results than ours but does not enjoy the same theoretical guarantees based on the theory of Randomized Signatures presented in Section 3.

heading: 5 CONCLUSIONS
text: Reservoir Computing and the Concept of Signature. The concept of reservoir computing is based on an simple and intriguing idea: generating a large number of features and then training a linear readout to perform an arbitrary downstream task. This simple approach promises to reduce the number of model parameters yet retaining a high level of expressive power. Inspired by results from rough path theory showing that the Signature Transform enjoys the status of reservoir, we propose a random reservoir, i.e. the Randomized Signature, to approximate the Signature Transform and thus inheriting its geometrical properties and expressive power. Following rigorous theoretical arguments, we do so by evolving a latent differential equation controlled by a path and evolved according to random vector fields. The empirical investigations presented in this paper provide evidence supporting that this strategy is successful and can be deployed on several tasks including non-parametric system identification and TSC. Future research and open questions. From a theoretical standpoint, Eq. 4 is written in terms of autonomous vector fields Ai for i = 1, · · · , d, whose entries are drawn according to Gaussian distributions with zero mean and unit variance, and an activation function σ. It is natural to wonder if choosing time-dependent matrices Ai leads to increased expressive power and to what extent the quality of the resulting features depend on a particular choice of σ. Furthermore, since our experiments suggest that the dimensionality k of the random features play a significant role in terms of expressive power, it is not trivial to identify a rigorous criterion to determine it. A precise characterization of this quantity and a clever way to select its value could potentially lead to further savings in terms of computational cost. From a practical point of view, our method proved to be able to accurately predict the response of different dynamical systems driven by heterogeneous and possibly irregular control paths, even when such controls were sampled from different laws with respect to the one used in training. In light of these results, it would be interesting to explore whether the method of Randomized Signatures could be extended to model the response of systems governed by controlled partial differential equations. In the context of TSC, we established a link with existing works using time series representations based on random features and our approach based on the update in Eq. 4. We implemented a first baseline algorithm combining our Randomized Signatures with the algorithm proposed in (Dempster et al., 2020), resulting in superior performance compared to strong deep learning baselines but sub-optimal results compared to the original algorithm. We believe such a gap can be filled by addressing the limitation of the current version of the algorithm, e.g. the assumptions it makes on the vector fields, as discussed at the end of Section 5. As a final consideration, we believe an interesting research direction being to investigate whether our random feature extraction pipeline could be incorporated into hierarchical models, such as deep neural networks, to form hybrid methods including both random and learned components operating in symbiosis.

heading: A THEORY DETAILS
text: A.1 BASIC DEFINITIONS
First of all, we define the concept of admissible tensor norms, which we assume to have in this work.
Definition 2 (Admissible Tensor Norms) Let E := Rd and ⊗ be a tensor product such that the tensor powers of E, ( E⊗k : k ≥ 1 ) , are equipped with a family (‖ · ‖E⊗k : k ≥ 1) of norms satisfying:
1. For j, k ∈ N and for all h ∈ E⊗j and l ∈ E⊗k
‖h⊗ l‖E⊗(j+k) ≤ ‖h‖E⊗j‖l‖E⊗k ;
2. For any permutation σ of {1, . . . , k}, ‖l1 ⊗ . . .⊗ lk‖E⊗k = ∥∥lσ(1) ⊗ . . .⊗ lσ(k)∥∥E⊗k ;
3. For any bounded linear functionals f onE⊗j and g onE⊗k, there exists a unique bounded linear functional, denoted as f ⊗ g, on E⊗(j+k) such that for all h ∈ V ⊗j and l ∈ E⊗k
f ⊗ g(h⊗ l) = f(h)g(l).
A family of tensor norms satisfying these conditions is called a family of admissible tensor norms.
Then we define the space in which the Signature of a path lies in, that is the following Tensor Algebra.
Definition 3 (Tensor Algebra T ( Rd ) ) We define the tensor algebra on Rd as
T ( Rd ) := ∞∏ k=0 ( Rd )⊗k
as well as its truncated version of order M ≥ 0 as
T M ( Rd ) := M∏ k=0 ( Rd )⊗k ,
where ( Rd )⊗k is the space of tensors of shape (d, . . . , d) given by Rd ⊗ · · · ⊗ Rd for k times.
Definition 4 (Concatenation Operation) We define the Concatenation Operation ∗ such that for any given couple of continuous piecewise smooth paths X : [0, s] → Rd and Y : [s, t] → Rd, their image through ∗ is the continuous piecewise smooth path X ∗ Y : [0, t]→ Rd defined by
X ∗ Yu := { Xu if u ∈ [0, s] Xs + Yu − Ys if u ∈ [s, t]
Definition 5 (Inverse Operation) We define the Inverse Operation←− such that for any continuous piecewise smooth path X : t → Rd, its image through←− is the continuous piecewise smooth path←− X t := XT−t, for each t ∈ [0, T ].
Definition 6 (Topological Space) A topological space is an ordered pair (X, τ), where X is a set and τ is a collection of subsets of X , satisfying the following axioms:
1. The empty set and X itself belong to τ ;
2. Any arbitrary (finite or infinite) union of members of τ still belongs to τ ;
3. The intersection of any finite number of members of τ still belongs to τ .
The elements of τ are called open sets and the collection τ is called a topology on X .
Definition 7 (Arcwise Connected Topological Space) A topological space (X, τ) is said to be arcwise connected if any two distinct points x, y ∈ X can be joined by an arc, that is a continuous map α : [0, 1]→ X such that α(0) = x and α(1) = y. We say that X is uniquely arcwise connected if for any two distinct points x, y ∈ X , there exists a unique path in X that joins them.
Definition 8 (R-tree) An R-tree is a uniquely arcwise connected metric space, in which the arc between two points is isometric to an interval.
Definition 9 (Tree-Like) A continuous piecewise smooth path x : [0, T ] → Rd is tree-like if there exists an R-tree τ , a continuous map φ : [0, T ] → τ and a map ψ : τ → V such that φ(0) = φ(T ) and x = ψ ◦ φ.
A.2 SUPPORTING RESULTS
Lemma 4 (Johnson-Lindenstrauss Lemma, (Cuchiero et al., 2021a)) Given an M -dimensional Hilbert space (H, 〈·, ·〉H) and Q a N -point subset of H , for any 0 < < 1, for each k ∈ N satisfying the so called Johnson-Lindenstrauss constraint
k ≥ 24 logN 3 2 − 2 3 ,
there exists a linear map f{Q,M,k} : H → Rk that embedsQ into Rk in an almost isometric manner. More specifically, we have that
(1− ) ‖a1 − a2‖2H ≤ ∥∥∥f{Q,M,k} (a1)− f{Q,M,k} (a2)∥∥∥2 ≤ (1 + ) ‖a1 − a2‖2H ,
for each a1,a2 ∈ Q.
Remark 1 When there is no need to specify the dependence of f{Q,M,k} from Q, M , or k, we will omit them up to referring to is simply as f . Finally, we define f∗ : Rk → H to be the adjoint map of f with respect to a fixed inner product 〈·, ·〉 in Rk.
To apply such Lemma in our context, we select M ≥ 0, equip TM ( Rd )
with an inner product such that
〈ei1 ⊗ · · · ⊗ eiM , ej1 ⊗ · · · ⊗ ejM 〉 := δi1j1 · · · δiM jM ,
where {ei1 ⊗ · · · ⊗ eiM }i1,...,iM∈{1,...,d} is the canonical basis of T M ( Rd ) . Therefore, we have
that ( TM ( Rd ) , 〈·, ·〉 ) is an Hilbert space.
Theorem 5 (Existence and Uniqueness of the Signature, Lemma 2.10 (Terry J. Lyons, 2004)) The following controlled differential equation
dSt = d∑ i=1 St ⊗ eidXit , S0 = 1
has a unique solution which, at each t ∈ [0, T ] is the Signature St of X on [0, t].
Notation 1 When there is no ambiguity about X , we will refer to its Signature as S, while we use SX if more paths are involved.
Theorem 6 (Theorem 2.29, (Terry J. Lyons, 2004)) Given a continuous piecewise smooth path X : [0, T ]→ Rd, its Signature is the banal Signature 1 := (1, 0, 0, · · · ) if and only if X is tree-like. In particular SX∗ ←− X = 1.
Now, we provide results to show the relevance of these features. First of all, the following theorem ensures that the Signature of a path encodes its essence and characterizes it completely.
According to Theorem 6, the concatenation X ∗ ←− X of X with its inverse1 ←− X has the same Signature as the constant path, but cannot be reparametrised to be constant. Similarly, if X,Y, Z are nonconstant paths, then SX∗Y ∗ ←− Y ∗Z∗ ←− Z ∗ ←− X = 1, but X ∗ Y ∗ ←− Y ∗ Z ∗ ←− Z ∗ ←− X is not a path of the form γ ∗←−γ for any path γ. The formal definition of tree-like path is given in Definition 10.
Theorem 7 (Characterizing Nature of the Signature, Corollary 1.4 (Ben Hambly, 2010)) Given a couple of continuous piecewise smooth paths X and X̂ , then SX = SX̂ if and only if X ∗ ←− X̂ is tree-like.
This result is actually much stronger as it implies that the solution of any differential equation controlled by X is fully determined by the vector fields and S. In particular, Theorem 1 shows that the solution of any differential equation controlled by X is essentially linear in S.
Adapting Theorem III.8. (Cuchiero et al., 2021b) as done in (Cuchiero et al., 2021a) one obtains:
Theorem 8 (Randomized Signature) For any fixed integer M ≥ 0, any fixed partition D = {t0, · · · , tN} of t such that 0 ≤ t0 < · · · < tN ≤ T , let us consider the Truncated Signature of order M of X at such times, that is Q := { SMt0 , · · · ,S M tN } such that its elements all lie
in ( TM ( Rd ) , 〈·, ·〉 ) . Let us now select 0 < < 1, k ∈ N satisfying the associated JohnsonLindenstrauss constraint, let f be the implied Johnson-Lindenstrauss map and f∗ its adjoint map. Then, the solution of the controlled differential equation in Rk
dZt = d∑ i=1 f (f∗ (Zt) ei) dX i t , Z0 = f(1)
on D, that is {Zt0 , · · · , ZtN }, are called the Randomized Signature of X on D. It holds that each Ztk is the projection of S M tk from ( TM ( Rd ) , 〈·, ·〉 ) onto Rk via f . To conclude, since f∗ (Ztk) is
close to SMtk in the norm of T M ( Rd ) , we conclude that Z preserves the geometric properties and the approximation power of S.
Adapting Theorem III.7 in (Cuchiero et al., 2021b) it can be be shown that (asymptotically) the JL projected vector fields stem from random matrices:
Theorem 9 (Z is a random dynamical system) For k 7−→ ∞, for each i ∈ {1, · · · , d} the linear vector fields f{k} ( f{k}∗(·)ei ) : Rk → Rk are square matrices with asymptotically normally distributed, independent entries.
Leveraging this results, we define the following:
Definition 10 (Localized Randomized Signature) For any random matrices A1, . . . , Ad in Rk×k and shifts b1, . . . , bd in Rk×1 such that maximal non-integrability holds on a random starting point z ∈ Rk, any fixed activation function σ, and d-dimensional control X , the solution of
dZt = d∑ i=1 σ (AiZt + bi) dX i t , Z0 = z, t ∈ [0, T ]. (7)
is called the Localized Randomized Signature ofX and we often refer to it as Randomized Signature.
Theorem 10 (Signature is a Reservoir (Restatement of Theorem 1)) Let Vi : Rm → Rm, i = 1, . . . , d be vector fields regular enough such that dYt = ∑d i=1 V i (Yt) dX i t , Y0 = y ∈ Rm, admits a unique solution Yt : [0, T ]→ Rm. Then, for any smooth test function F : Rm → R and for every M ≥ 0 there is a time-homogeneous linear operator L : T M ( Rd ) → R which depends only on
(V1, . . . , Vd, F,M, y) such that F (Yt) = L ( SMt ) +O ( tM+1 ) , and t ∈ [0, T ].
1←−X t := XT−t
Proof: We provide just a sketch. For every smooth test function F : Rm → R, the formal Taylor expansion along controls together with the fundamental theorem of calculus gives us that
F (Yt) = F (y) + d∑ i=0 ∫ t 0 ViF (Ys) dX i s, t ≥ 0. (8)
For y ∈ Rm, this equation can be inserted into itself leading to a generalized version of Taylor expansion for controlled ordinary differential equations
M∑ k=0 ∑ (i1,...,ik)∈{0,...,d}k Vi1 . . . VikF (y) ∫ 0≤t1≤...≤tk≤t dXi1t1 . . . dX ik tk +RM (F, t)
for M ≥ 0, with the remainder satisfying
RM (F, t) = ∑
(i1,...,iM+1)∈{0,...,d}M+1
∫ 0≤t1≤...≤tM+1≤t Vi1 . . . ViM+1F (Yt1) dX i1 t1 . . . dX iM+1 tM+1 .
A.3 PRACTICAL EXAMPLE OF THE COMPUTATION OF SIGNATURES (ONE DIMENSION).
To enhance intuition, we exemplify the computation of Signatures in the simplest setting possible. Let X : [0, T ] → R, then S1t = ∫ t 0 dXs — which is exactly Xt −X0. To get S2t , we instead have
to compute the following iterated integral: S2t = ∫ t
0 (∫ v 0 dXs ) dXv . Signatures of higher order S j t
are computed in a similar way, by iteratively integrating the path j times. As a practical example, let Xt = t. Then it is easy to see that S j t = tj j! . Now let Y (t) be an analytic function of time for
which we have Y (t) = ∑∞ j=0 Y (j) 0 tj
j! . Taylor’s theorem combined with the previous computation of the Signatures implies that Y can be approximated as a linear map of the Signatures of t. Finally, note that Sjt in this case gets smaller and smaller in magnitude as j increases. This suggests that the truncated Signature can be safely used to approximate Y .
Computational Complexity and dimensionality of Signatures. The computational complexity for computing SMt is O(dM ). Indeed, consider d = 2: S2t is a 2 × 2 matrix with elements ∫ t 0 (∫ v 0 dX1s ) dX1v , ∫ t 0 (∫ v 0 dX1s ) dX2v , ∫ t 0 (∫ v 0 dX2s ) dX1v and ∫ t 0 (∫ v 0 dX2s ) dX2v . For M =
3, the object to compute is instead a 2 × 2 × 2 tensor, containing all integrals of the type∫ t 0 (∫ w 0 (∫ v 0 dXi1s ) dXi2v ) dXi3w for all i1, i2, i3 ∈ {1, 2}. Hence, the complexity — as well as the features dimensionality — scales exponentially in M .

heading: B EXPERIMENTS
text: In this section, we prove that our approach can accurately learn the dynamics of several rough differential equations. First of all, we focus on dynamics of 1-dimensional stochastic Double Well which is characterized by a highly nonlinear dynamics. Secondly, we study the multidimensional case of the 4-dimensional Ornstein–Uhlenbeck. Then, we investigate the performance on the rough dynamics of Fractional Geometric Brownian Motion even for choices of extremely low values of the hurst parameter H . We also obtain promising results in learning the dynamics of an EnzymeSubstrate Reactions. Next, we show how in high dimensions randomized signatures are computationally efficient, and also present one experiment in comparison with neural controlled differential equations (Kidger et al., 2020). To conclude, we show how the randomized signature approach can be used in cases where the time grid is irregularly sampled.
B.1 PRELIMINARY RESULTS: ROUGH DIFFERENTIAL EQUATIONS
1-Dimensional Stochastic Double Well. Let us recall that the dynamics of the 1-Dimensional Stochastic Double Well process is given by
dYt = θYt ( µ− Y 2t ) dt+ σdWt, Y0 = y0 ∈ R, t ∈ [0, 1].
where Wt is a 1-dimensional Brownian motion, and (µ, θ, σ) ∈ R×R+×R+. Let us fix y0 = 1 and (µ = 2, θ = 1, σ = 1), and the partition D of [0, 1] to have N = 101 equally spaced time steps. We train a Ridge Regression with regularization parameter λ = 0.001 to map instances of Randomized Signature of the controls Z into the respective solution Yt. We repeat the experiment on different values of the number NTrain of train samples and dimension k of the Randomized Signature. On the left of Figure 9, we plot an example of the trajectory of Z while, on its right, we plot the comparison of the true and the generated time series on an out-of-sample case. The following table shows the performance in terms of L2 relative error on 10000 test samples:
4-Dimensional Ornstein–Uhlenbeck process. Let us recall that the dynamics of the 4- Dimensional Ornstein–Uhlenbeck process is given by
dYt = (µ−ΘYt) dt+ ΣdWt, Y0 = y0 ∈ R4, t ∈ [0, 1]
where Wt is a 4-dimensional Brownian motion, and (µ,Θ,Σ) ∈ R4 × R4×4 × R4×4. Let us fix y0 = 1, µ = 1, Σ = 14, Θi.j = i/j, the partition D of [0, 1] to have N = 101 equally spaced time steps, and k = 708. Finally, we train a Ridge Regression with λ = 0.001 on NTrain train sample and Figure 10 shows the comparison of out-of-sample generated and true trajectories while Table 2 reports the performance in terms of L2 relative error on 10000 test samples.
1-Dimensional Geometric Fractional Brownian Motion. Let us recall that the dynamics of the 1-Dimensional Geometric Fractional Brownian Motion process given by
dYt = Yt
( µdt+ σdB
(H) t ) , Y0 = y0 ∈ R+,
where µ ∈ R, σ ∈ R \ {0}, and B(H)t is a one-dimensional fractional Brownian motion of hurst parameter H ∈ (0, 1). In this experiment, we fix y0 = 1 and (µ = 1, σ = 1), we let H vary in {0.1, 0.2, 0.3, 0.4, 0.5}, and the partitionD of [0, 1] is made ofN equally spaced times. We highlight that, as per Theorem D.3.2 in (Biagini et al., 2008), when H < 1/4, the fractional Brownian motion is so rough that it is no longer a Rough Path and its Signature is no longer defined. Finally, we train a Ridge Regression with λ = 0.001 on NTrain = 10000 train samples of Z into the logarithm of the respective solution log(Yt) which implies that we need to then apply the exp function to our predictions. The right of Figure 11 shows the comparison of the generated and true trajectories on an out-of-sample trajectory while the following table shows the Relative L2 Error as we vary the number of time steps N and hurst parameter H:
Enzyme-Substrate Reactions. First, we consider a set of highly nonlinear CODEs modeling the interaction between a catalyst (enzyme, with abundance Ct) and a substrate S (with abundance St). Substrate (e.g. lactose) is injected via the control variable Xt, reacts with the enzyme (e.g. lactase)
and the concentration Yt of the product chemical (e.g. glucose) is observed (Ingalls, 2013). dSt = (Ct − k1St (1− Ct)) dt+Xtdt dCt = − (Ct − k1St (1− Ct)) dt− k2Ctdt dYt = k2Ctdt.
Following (Ingalls, 2013) we choose (k1, k2) = (30, 10), set (S0, P0, Y0) = (0, 0, 0) and consider the evolution on 0 ≤ t ≤ 1. We fix the time grid to have N = 101 equally spaced time steps and the control Xt to follow the law of W 2t where Wt is a 1-dimensional Brownian Motion (to ensure positivity). We train a Ridge Regression with regularization parameter λ = 10−3 to map 105 instances of 222-dimensional Randomized Signature Z of the controls into the respective solution Yt. On the left of Figure 12, we plot the comparison of the true and the generated time series on a test sample. As we can see, the model has learnt to correctly map a trajectory of Xt to the respective system response Yt. More surprisingly, the right of such a figure shows that our model is able to predict the correct output even if we stimulate the system with substrate injection that follows a completely different law such as 0.5 · 1{W 2t >0.5}. This suggests that the system was correctly identified out of distribution.
Comparison with Neural Controlled Differential Equations In this experiment, we benchmark the performance of our Randomized Signature model with that of Neural Controlled Differential Equations (NCDEs) (Kidger et al., 2020), a recently proposed competitive deep learning technique for controlled dynamical system modelling. We test the two models on the task of learning the dynamics of the 1-dimensional Ornstein-Uhlenbeck Process. As usual, to train our model, we fix an equally spaced partition D of [0, 1], extract the Randomized Signature of the controls for NTrain samples, and train a Ridge Regression to map them into the respective integrated solutions. The NCDE models instead parametrize the potentials of a latent controlled differential equation with feed-forward neural networks with 1 hidden layer of nnodes. Figure 13 shows the comparison of our model against several version of the NCDEs. Additionally, Table 4 and Table 5 show the number of trainable parameters, average L2 Relative Errors, and training times for all the models forNTrain = 100 and NTrain = 1000 respectively. We observe that our model is the best one from all aspects and that not even a considerable amount of training data is beneficial for NCDEs.
B.2 ADDITIONAL DETAILS
Learning the Signatures. In this experiment, we fix the partition D = {t0, · · · , tN} of [0, 1] to have 201 equally-spaced time steps and the regularization parameter of the Ridge Regression is λ = 0.001.
Reconstruction and Randomized Signature Autoencoders. In this experiment, we fix the partition D = {t0, · · · , tN} of [0, 1] to have 101 equally-spaced time steps, fix k = 111, and consider a single realization XTrain of an Ornstein–Uhlenbeck process of asymptotic mean µ = 2, mean reversion θ = 1, and volatility σ = 1. We train a Ridge Regression with regularization parameter λ = 0.001 to map the Randomized Signature Z of the one and only train sample XTrain into XTrain itself.
Comparison with Vanilla Signature Approach Let us recall that the dynamics of the dDimensional Fractional Ornstein–Uhlenbeck process is given by
dYt = (µ−ΘYt) dt+ ΣdBHt , Y0 = y0 ∈ Rd, t ∈ [0, 1]
where BHt is a d-dimensional Fractional Brownian motion of hurst parameter H ∈ (0, 1), and (µ,Θ,Σ) ∈ Rd ×Rd×d ×Rd×d. Let us fix y0 = 1, µ = 1, Σ = 1d, Θi.j = i/j, the partition D of [0, 1] to have N = 101 equally spaced time steps, and H = 0.3.
1-Dimensional Stochastic Double Well - Irregularly Sampled Time Grid. Let us recall that the dynamics of the 1-Dimensional Stochastic Double Well process is given by
dYt = θYt ( µ− Y 2t ) dt+ σdWt, Y0 = y0 ∈ R, t ∈ [0, 1] (9)
where Wt is a 1-dimensional Brownian motion, and (µ, θ, σ) ∈ R × R+ × R+. For each train and test sample, the partition D of [0, 1] is made of N randomly drawn times. More precisely, D = {0, t1, · · · , tN−1, 1} such that tk = 1/(1− exp(−sk)) and {s1, · · · , sN−1} are N − 2 independent realizations of a uniform distribution U [0, 1] sorted increasingly. As a result, the probability that two samples share the sameD is null. We train a Ridge Regressor withNTrain = 10000 train samples and Figure 11 shows the comparison of an out-of-sample generated and true trajectory. Finally, Table 7 shows the Relative L2 Error on 10000 test samples as we vary the number of time steps N and k and while comparing it to the respective experiment in case the time grid is regularly spaced. As we can see, even though the performance are worse than the regularly sampled setup, this technique proves to be anyway reliable on irregularly sampled regimes. In this experiment, we fix y0 = 1 and (µ = 2, θ = 1, σ = 1) while the regularization parameter of the Ridge Regression is λ = 0.001.
Tumor Growth Kinetics. We set (k1, k2, λ0, λ1,Ψ) = (10, 0.5, 0.9, 0.7, 20), and( u10, u 2 0, u 3 0, u 4 0, Y0 ) = (2, 0, 0, 0, 2). We fix the partition D of [0, 10] to have 101 equally spaced time steps and we again train a Ridge Regression with regularization parameter λ = 0.001 to map 10000 instances of the 222-dimensional Randomized Signature Z of the controls into the respective solution Yt.
Electrochemical Battery Model. In this experiment, we fix the partition D of [0, 500] to have 1001 equally spaced time steps, the regularization parameter of the Ridge Regression is λ = 10−3 while k ∈ {50, 166, 1000}. Regarding the latest experiment, we simply learn to map the Random Signature of the controls into the respective Y to which we added white noise with variance 0.01. Regarding NNARX, for the regressor we follow the results presented in (Masti and Bemporad, 2018), where na = nb = 12, and nk = 1 lead to the best results. Regarding the neural network, we use a feedforward neural network with input dimension na + nb + 1 = 25 and 2 hidden layers each with either 22 (NNARX22) or 1000 (NNARX1000) hidden units. For the loss function, we used the mean square error optimized using Adam with learning rate 0.01 for 100 epochs. Table 8 reports the performances in terms of MSE on 1000 test samples for our models compared to NNARX.
Time Series Classification
DATASET DETAILS We run our TSC experiments on the UCR archive, one of the most popular benchmarks for TSC methods (Dau et al., 2019; Fawaz et al., 2019). It consists of 128 heterogeneous datasets of different sizes comprising possibly irregularly sampled one-dimensional time series of various lengths. Rocket is one of the best performing algorithms on this benchmark.
ROCKET When implementing Rocket, we generate k kernels each characterized by a certain random length lfilter, random weights ω ∈ Rlfilter , and random jump j. Since the lengths and the jumps are randomly sampled, the length M of the vector Z = (Z1, · · · , ZM ) of the random features varies accordingly. Given that the number of possible combinations of the values of lfilter and
j is limited, we can group the kernels based on such combinations. As a result, among each of these groups, the length M is fixed and, additionally, each kernel is applied to the very same sub-timeseries ( Xi, Xi+j, · · · , Xi+(lfilter−1)×j ) of X starting from position i. If we now focus on one of these subgroups and let us assume it contains k̃ kernels of length lfilter and jump j. As a result of the application of each kernel on X , we obtain the following matrix Z0,1 Z1,1 . . . . . . ZM,1 Z0,2 Z1,2 . . . . . . ZM,2 ... ... ... . . .
... Z0,k̃ Z1,k̃ . . . . . . ZM,k̃  where Zi,r is the result of the action of the r-th kernel on the sub-time-series( Xi, Xi+j, · · · , Xi+(lfilter−1)×j ) of X starting from position i. Finally, out of each row, we extract the features described above.
Turning to our approach, we generate k kernels each characterized by a certain length lfilter, weights ω = (ω0, ω1) ∈ R2, and jump j. As before, let us again group them based on all the possible combinations of lfilter and j and focus on only one of such subgroups. As a consequence, we have k̃ kernels each characterized by weights (ωr0, ω r 1), r ∈ {1, · · · , k̃} with fixed length and jump. The application of each of these k̃ kernels following the update rules proposed in Algorithm 1 is equivalent to extracting a 1-dimensional Randomized Feature of each of the sub-time-series( Xi, Xi+j, · · · , Xi+(lfilter−1)×j
) . As a result, we obtain the following matrix Z0,1 Z1,1 . . . . . . ZM,1 Z0,2 Z1,2 . . . . . . ZM,2 ... ... ... . . .
... Z0,k̃ Z1,k̃ . . . . . . ZM,k̃  where Zi,r is the 1-dimensional Randomized Signature of ( X0, Xj, · · · , X(lfilter−1)×j ) via the r-th kernel. Let us we define the random matrix A guiding the evolution of the Randomized Signature as the k̃× k̃ diagonal matrix having the weights ωr0 on the diagonal while the random bias b is a column vector having the ωr1 as entries. Then, the i-th column is the k̃-dimensional Randomized Signature of the sub-time-series ( Xi, Xi+j, · · · , Xi+(lfilter−1)×j ) . Finally, out of each row, we extract the features described above.
Comparison of Computational Complexities and Dimensions The following table reports key figures regarding the computational complexity and dimensionality of both Signature and Randomized Signature.

Review summaries:
This paper proposes a Signature-based method for learning feature representation of time series, denoted as $d$-dimensional path $X$. By definition, Signature is a sequence of $M$ tensors where each tensor is a sum of integral terms defined on dimensions of  $X$. Theoretically Signature is guaranteed to linearly explain the solution of any differential equation driven by $X$ which implies that it encodes well $X$. Computing Signature however has time complexity $O(d^M)$ rendering it impractical for learning purposes. To alleviate the issue, the paper introduces Randomized Signature, which can be computed in $O(k^2 d)$ only where $k$ is the number of Randomized Signatures. Each of such Signature has $O(kd)$ dimensions. As stated by Theorem 2 in the paper, Randomized Signature approximates well the original Signature. Experiments on differential equations show Randomized Signature finds solutions close to ground-truth solutions. Experiments on time series classification also show its potential.