Embedding and learning with signatures

Introduction

Sequential or temporal data are arising in many fields of research, due to an increase in storage capacity and to the rise of machine learning techniques. An illustration of this vitality is the recent relaunch of the Time Series Classification repository (Bagnall et al. 2018), with more than a hundred new datasets. Sequential data are characterized by the fact that each sample consists of an ordered array of values. The order need not correspond to time, for example, text documents or DNA sequences have an intrinsic ordering, and are, therefore, considered as sequential. Besides, when time is involved, several values can be recorded simultaneously, giving rise to an ordered array of vectors, which is, in the field of time series, often referred to as multidimensional time series. To name only a few domains, market evolution is described by financial time series, and physiological variables (e.g., electrocardiograms, electroencephalograms) are recorded simultaneously in medicine, yielding multidimensional time series. Finally, smartphone and GPS sensors data, or character recognition problems, present both spatial and temporal aspects. These high-dimensional datasets open up new theoretical and practical challenges, as both algorithms and statistical methods need to be adapted to their sequential nature.

Different communities have addressed this problem. First, time series forecasting has been an active area of research in statistics since the 1950s, resulting in several monographs, such as Hamilton 1994, Box et al. 2015 and Shumway and Stoffer 2017, to which the reader is referred for overviews of the domain. Time series are considered as realizations of various stochastic processes, such as the famous ARIMA models. Much work in this field has been done on parameter estimation and model selection. These models have been developed for univariate time series but have been extended to the multivariate case (Lütkepohl 2005), with the limitation that they become more complicated and harder to fit.

More recently, the field of functional data analysis has extended traditional statistical methods, in particular regression and Principal Component Analysis, to functional inputs. Ramsay and Silverman 2005; Ferraty and Vieu 2006 provide introductions to the area. Kokoszka et al. 2017 give an account of recent advances. In particular, longitudinal functional data analysis is concerned with the analysis of repeated observations, where each observation is a function (Greven et al. 2011; Park and Staicu 2015). The data arising from this setting may be considered as a set of vector-valued functions with correlated coordinates, each function corresponding to one subject and each coordinate corresponding to one specific observation.

Although these various disciplines work with sequential data, their goals usually differ. Typically, time series analysis is concerned with predicting future values of one observed function, whereas (longitudinal) functional data analysis usually collects several functions and is then concerned with the prediction of another response variable. However, all these methods rely on strong assumptions on the regularity of the data and need to be adapted to each specific application. Therefore, modern datasets have highlighted their limitations: a lot of choices, in basis functions or model parameters, need to be handcrafted and are valid only on a small-time range. Moreover, these techniques struggle to model multidimensional series, in particular, to incorporate information about interactions between various dimensions.

On the other side, time series classification has attracted the interest of the data mining community. A broad range of algorithms have been developed, reviewed by Bagnall et al. 2017 in the univariate case. Much attention has been paid to the development of similarity measures adapted to temporal data, a popular baseline being the Dynamic Time Warping metric (Berndt and Clifford 1996), combined with a 1-nearest neighbor algorithm. Bagnall et al. 2017 state that this baseline is beaten only by ensemble strategies, which combine different feature mappings. However, a great limitation of these methods is their complexity, as they have difficulty handling large time series. Recently, deep learning seems to be a promising approach and solves some problems mentioned above. For example, Fawaz et al. 2019 claim that some architectures perform systematically better than previous data mining algorithms. However, deep learning methods are costly in memory and computing power, and often require a lot of training data.

The signature has recently received the attention of the machine learning community and has achieved a series of successful applications. To cite some of them, Yang et al. 2016 have achieved state-of-the-art results for handwriting recognition with a recurrent neural network combined with signature features. Graham 2013 used the same approach for character recognition, and Gyurkó et al. 2014 coupled Lasso with signature features for financial data streams classification. Kormilitzin et al. 2016 investigated its use for the detection of bipolar disorders, and Yang et al. 2017 for human action recognition. For introductions to the signature method in machine learning, the reader is referred to the work of Levin et al. 2013 and to Chevyrev and Kormilitzin 2016.

However, despite many promising empirical successes, a lot of questions remain open, both practical and theoretical. In particular, to compute signatures, it is necessary to embed discretely sampled data points into paths. While authors use different approaches, this embedding is only mentioned in some articles, and rarely discussed. Thus, the purpose of this paper is to take a step forward in understanding how signature features should be constructed for machine learning tasks, with a special focus on the embedding step. The article is organized as follows.

In Section 2, a brief exposition of the signature definition and properties is given, along with a survey of different approaches undertaken in the literature to combine signatures with machine learning algorithms. Datasets used throughout the paper are also presented in Section 3.

In Section 4, potential embeddings are reviewed and their predictive performance is compared with an empirical study on 3 real-world datasets. This study indicates that the embedding is a step as crucial as the algorithm choice since it can drastically impact accuracy results. In particular, we find that the lead-lag embedding systematically outperforms other embeddings, consistently over different datasets and algorithms. This finding is reinforced by a simulation study with autoregressive processes in Section 5.

In Section 6 the choice of signature domain is investigated. Signatures can be computed on any sub-interval of the path definition domain, and it is natural to wonder whether some local information is lost when signatures of the whole path are computed. The section ends by showing that, with a good embedding, the signature combined with a simple algorithm, such as a random forest classifier, obtains results comparable to state-of-the-art approaches in different application areas, while remaining a generic approach and computationally simple.

In Section 7 some open questions for future work are discussed.

These empirical results are based on three recent datasets, in different fields of application. One is a univariate sound recording dataset, called Urban Sound (Salamon et al. 2014), whereas the others are multivariate. One has been made available by Google 2017, and consists of drawing trajectories, while the other is made up of 12 channels recorded from smartphone sensors (Malekzadeh et al. 2018). They are each of a different nature and present a variety of lengths, noise levels, and dimensions. In this way, generic and domain-agnostic results are obtained. The code is available at https://github.com/afermanian/embedding_with_signatures.

A first glimpse of the signature method

where the supremum is taken over all finite partitions

The set of bounded variation paths is exactly the set of functions whose first derivatives exist almost everywhere. Being of bounded variation is therefore not a particularly restrictive assumption. It contains, for example, all Lipschitz functions. In particular, if $X$ is continuously differentiable, and $\dot{X}$ denotes its first derivative with respect to $t$, then

where $X=(X^{1},\dots,X^{d})$, and $Y=(Y^{1},\dots,Y^{d})$. When $X$ is continuously differentiable, this integral is equal to the standard Riemann integral, that is,

As an example, assume that $X$ is linear, i.e.,

The formula above is useful since in practice only integrals of linear paths are computed, as discussed later in this subsection. It is now possible to define the signature.

$S^{I}(X)_{[s,t]}$ is then said to be a signature coefficient of order $k$.

The signature of $X$ is the sequence containing all signature coefficients, i.e.,

The signature of $X$ truncated at order $K$, denoted by $S_{K}(X)$, is the sequence containing all signature coefficients of order lower than or equal to $K$, that is

For simplicity, when $[s,t]=$, the interval is omitted in the notations, and, e.g., $S_{K}(X)$ is written instead of $S_{K}(X)_{}$.

From these definitions, it follows that the linear interpolation of a (multivariate) time series observed on a finite time horizon will be of bounded variation, and therefore that its signature is well defined. Note that this procedure of mapping a discrete time series into a continuous path is called an embedding, and linear interpolation is only one embedding among others, which will be studied in Section 4. For example, Brownian motion is not of bounded variation but is instead of finite $p$-variation for any $p>2$. However, its signature can still be defined with Itô or Stratonovitch integrals.

and $K+1$ if $d=1$. Unless otherwise stated, it is assumed that $d\neq 1$, as this is in practice usually the case. Thus, the size of $S_{K}(X)$ increases exponentially with $K$, and polynomially with $d$—some typical values are presented in Table 1.

Similarly, coefficients of order 3 can be written as a tensor of order 3, and so on. Then, $S(X)$ can be seen as an element of the tensor algebra

This structure of the tensor algebra will not be used in the present article but is used to derive properties of the signature (Lyons 1998; Friz and Victoir 2010; Hambly and Lyons 2010).

It should be noted that due to the ordering in the integration domain in (2), the signature coefficients are not symmetric. For example, $S^{(1,2)}(X)$ is not the same as $S^{(2,1)}(X)$. Finally, Chevyrev and Kormilitzin 2016 show how, under certain assumptions, empirical statistical moments can be explicitly recovered from signature coefficients. Typically, the empirical mean can be recovered from signature coefficients of order 1, the variance from coefficients of order 2, and so on. Therefore, the larger the truncation order, the more detailed the information encoded in the signature.

As a toy example, consider the linear path (1) again, and assume for simplicity that $d=2$:

Then, for any $[s,t]\subset$ the signature coefficients of order $1$ are

For any index $I=(i_{1},\dots,i_{k})\subset\{1,2\}^{k}$, it is easily obtained that

A crucial feature of the signature is that it encodes geometric properties of the path. Indeed, coefficients of order 2 correspond to some areas outlined by the path, as shown in Figure 1. For higher orders of truncation, the signature contains information about the joint evolution of tuples of coordinates (Yang et al. 2017). Furthermore, the signature possesses several properties that make it a good statistical summary of paths, as shown in the next four propositions.

This proposition is a consequence of the properties of integrals and bounded variation paths (Friz and Victoir 2010, Proposition 7.10). In other words, the signature of a path is the same up to any reasonable time change. There is, therefore, no information about the path parametrization in signature coefficients. However, when relevant for the application, it is possible to include this information by adding the time parametrization as a coordinate of the path. This procedure plays a decisive role in the construction of time embeddings, which will be thoroughly discussed in Section 4.

A second important property is a condition ensuring the uniqueness of signatures.

If $X$ has at least one monotone coordinate, then $S(X)$ determines $X$ uniquely up to translations.

It should be noticed that having a monotone coordinate is a sufficient condition, but a necessary one can be found in Hambly and Lyons 2010, together with a proof of this proposition. The principal significance of this result is that it provides a practical procedure to guarantee signature uniqueness: it is sufficient to add a monotone coordinate to the path $X$. For example, the time embedding mentioned above will satisfy this condition.

This result does not provide a practical procedure to reconstruct a path from its signature. However, this is an active area of research (Chang et al. 2017; Lyons and Xu 2017; Lyons and Xu 2018). In particular, Lyons and Xu 2017 derive an explicit expression of rectilinear paths, defined in Section 4.1, in terms of their signatures; and Lyons and Xu 2018 construct, from the signature of a ${\mathcal{C}}^{1}$ path, a sequence of piecewise linear approximations converging to the initial path.

The next proposition reveals that the signature linearizes functions of $X$. We refer the reader to Király and Oberhauser 2019 for a proof.

This proposition is a consequence of the Stone-Weierstrass theorem. The classical Weierstrass approximation theorem states that every real-valued continuous function on a closed interval can be uniformly approximated by a polynomial function. Similarly, this theorem states that any real-valued continuous function on a compact subset $D$ of bounded variation paths can be uniformly approximated by a linear form on the signature. Linear forms on the signature can, therefore, be thought of as the equivalent of polynomial functions for paths.

This proposition is an immediate consequence of the linearity property of integrals (Lyons et al. 2007, Theorem 2.9). However, it is essential for the explicit calculation of signatures. Indeed, in practice, $X$ is observed at a finite number of times and becomes by interpolation a continuous piecewise linear path. To compute its signature, it is then sufficient to iterate the following two steps:

Compute with equation (4) the signature of a linear section of the path.

Concatenate it to the other pieces with Chen’s formula (5).

and the signature corresponds exactly to an infinite sequence of the moments of the path.

2 Signature and machine learning

In this notation, $x^{k}_{i,j}$ denotes the $k$th coordinate of the $i$th sample observed at time $t_{j}$.

The assumption that $Y_{i}$ is a real number excludes several situations from our study. For example, the goal of functional longitudinal data analysis is usually the prediction of the next functional profile, which does not fall within our setting. Similarly, prediction of functional responses, which are a topic of interest in functional data analysis, or of multiple time points in time series analysis, are not considered.

As an example, consider the Google dataset Quick, Draw! (Google 2017). It consists of the pen trajectories of millions of drawings, divided into 340 classes. Some examples are shown in Figure 2. In this case, the $y_{i}$ are discrete labels of the drawing’s class, and the $\mathbf{x_{i}}$ are matrices of pen coordinates. In this example, $d=2$ and $p_{i}$ varies for each drawing but is typically in the order of a few dozen points.

The first approach is to compute the signature of $X$ on its whole domain, that is, on $$. In this way, the time-dependent input$X$ is mapped into a time-independent finite set of coefficients, that is then fed into a predictive algorithm, typically a feedforward neural network. Any time-independent additional covariates may be added to this feature set. This strategy is implemented by Yang et al. 2017 for skeleton-based human action recognition. From a sequence of human joints’ positions, the authors construct a high dimensional vector of signature coefficients, which is then the input of a small dense network. Gyurkó et al. 2014; Lyons et al. 2014 also apply this method to financial time series, combining it with Lasso and ordinary least squares regression.

A second family of methods consists in describing the input path by a sequence of signature coefficients. There are several variants of this approach, but the one of Wilson-Nunn et al. 2018 is presented here in detail for its simplicity and representativeness. To create a signature sequence, the time interval $$ is divided into a dyadic partition

To sum up, the signature may be used in various ways, and for different applications. Several points of view coexist, and none of them has shown to be systematically better. In particular, on the one hand, signatures may be used to remove temporal aspects and to reduce the dimension of the problem, whereas, on the other hand, they may do the opposite and increase the dimension of the algorithm’s input. Moreover, they are combined with various learning algorithms and it may be hard to distinguish the properties of the signature from those of the algorithms. Nevertheless, all these methods assume that discrete data points have been embedded into actual continuous paths. As will be seen in Section 4, the choice of path is crucial. Therefore, we describe in the next section the datasets used throughout the article to understand their underlying structure and find suitable embeddings.

Datasets

The datasets used in this article have been chosen to cover a broad range of applications while being recent and challenging in various ways. Moreover, they present a variety of sampling frequencies and dimensions. They illustrate therefore different potential embeddings.

First, the Quick, Draw! dataset (Google 2017), which was already discussed in Section 2.2, and illustrated in Figure 2, is a public Google dataset. It consists of $50$ million drawings, each drawing being a sequence of time-stamped pen stroke trajectories, divided into 340 categories. It takes approximately 7 gigabytes of hard disk space and is, therefore, a particularly large dataset. To compute the signature of every sample, it would thus be necessary to design a specific architecture, which cannot be implemented on a standard laptop computer. However, the goal of this study is not to achieve the best possible performance, but to understand embedding properties. Moreover, the experiments should be easily reproducible without requiring much computational ressources. Therefore, only a subset of the data is used: 68 000 training samples in Sections 4.2 and 6.1, 12 million in Section 6.2.

The embedding

Rectilinear path

Another interpolation method is often used in the literature (Chevyrev and Kormilitzin 2016; Kormilitzin et al. 2016) and referred to as an “axis path" or “rectilinear path”. It is also piecewise linear but each linear section is parallel to an axis. In other words, to move from one point $(x_{j}^{1},x_{j}^{2})$ to another point $(x_{j+1}^{1},x_{j+1}^{2})$, a first linear segment goes from $(x_{j}^{1},x_{j}^{2})$ to $(x_{j+1}^{1},x_{j}^{2})$, parallel to the x-axis, and a second segment from $(x_{j+1}^{1},x_{j}^{2})$ to $(x_{j+1}^{1},x_{j+1}^{2})$, parallel to the y-axis. This path is depicted in Figure 7(b). A crucial aspect of this path is that there exists a simple way to reconstruct it from its signature features (Lyons and Xu 2017). Note that for unidimensional data, such as the Urban Sound dataset, the linear and rectilinear interpolations are identical.

Time path

The third approach builds upon the linear path and enriches it by adding a monotone coordinate. This ensures the uniqueness of the signature, as stated in Proposition 2. It usually corresponds to adding the time parametrization as a coordinate of the path, as is done by Yang et al. 2017. Therefore, if $t\mapsto(X^{1}_{t},X^{2}_{t})$ is the linear path described above, which is piecewise linear, the time embedding is the $3$-dimensional path $t\mapsto(X^{1}_{t},X^{2}_{t},t)$, shown in Figure 7(c).

Lead-lag path

In this definition, $X^{1}$ and $X^{2}$ are a linear interpolation of the sequence

in which the last point is repeated twice, and

Stroke path

For the Quick, Draw! data, extra information about pen jumps is provided. In the context of Arabic handwriting recognition, Wilson-Nunn et al. 2018 have introduced the idea of encoding information about jumps into a new coordinate. In essence, the approach is to use a $3$-dimensional path in which the last dimension corresponds to strokes, in a similar way to the encoding of matrix (8). This procedure can be deployed in various ways and we restrict our attention to three of them.

Finally, for comparison purposes, a strictly monotone coordinate is also considered. It has jumps of $1$ when a new stroke begins, and otherwise grows linearly inside one stroke, such that it has increased by $1$ between the beginning and the end of the stroke. In this definition, the goal is to check whether having a strictly monotone coordinate increases accuracy, while in the two previous versions the stroke coordinate is piecewise constant. This embedding can be seen as a mix between time and stroke paths and could inherit the good properties of both. The resulting path is called “version 3" and shown in Figure 7(f).

To conclude, there exists a broad range of embeddings, living in spaces of various dimensions. They lead to different signature features, which therefore do not have the same statistical properties. The embedding choice will prove to have a significant influence on accuracy.

2 Results

In this subsection, the results of our study on embedding performance are presented. To this end, the first approach described in Section 2.2 is implemented. Starting from the raw data, it is first embedded into a continuous path, then its truncated signature is computed and used as input for a learning algorithm. The embeddings described in the previous section are used. Note that the lead-lag path is taken with lag $1$, but other lags will be discussed in Section 6.2. Each feature is normalized by the absolute value of its maximum so that all input values lie in $$. The findings should be independent of the data and the underlying statistical model so a range of different algorithms is used. Their hyperparameters have been set to their default values, without trying to optimize them for each dataset. Indeed, the goal is not to select the best algorithm or to achieve a particularly good accuracy, but rather to compare the performance of different embeddings. The classification metric to assess prediction quality is the accuracy score. Denoting by$(y_{1},\dots,y_{n_{\text{test}}})$the test set’s labels, and$(\hat{y}_{1},\dots,\hat{y}_{n_{\text{test}}})$ the predicted labels, this score is defined by

The four following algorithms have been used throughout the study.

Following Yang et al. 2017, a dense network with one hidden layer composed of 64 units with linear activation functions is first considered. A softmax output layer and the categorical cross-entropy loss are used, which yields a linear model equivalent to logistic regression. This architecture is a sensible choice, since Proposition 3 states that linear functions of the signature approximate arbitrarily well any continuous function of the input path. The Python library keras (Chollet et al. 2015), with TensorFlow backend, is used. The network is regularized by adding a dropout layer after the input layer, with a rate of 0.5. Optimization is done with stochastic gradient descent with an initial learning rate of 1. It is reduced by $2$ when no improvement is seen on a validation set during 10 consecutive epochs. The maximal number of epochs is set to 200 and the mini-batch size to 128.

Furthermore, the performance of a random forest classifier with 50 trees, implemented in scikit-learn (Pedregosa et al. 2011), is tested. It is a nonlinear very popular method initially proposed by Breiman 2001.

The XGBoost algorithm, introduced by Chen and Guestrin 2016, and implemented in the Python package xgboost, is also used. It is a state-of-the-art gradient boosting technique, building upon the work of Friedman 2001. The maximum number of iterations is set to 100 and early stopping with a patience of $5$ is used to prevent overfitting and speed up training. The maximum depth of a tree is set to 3 and the minimum loss reduction to make a split to 0.5.

Finally, a nearest neighbor classifier is run with a default value of 5 neighbors. This method is known to suffer from the curse of dimensionality, so it is of interest to see how the signature truncation order affects its performance.

For each of the algorithms described above and each dataset of Section 3 (Quick, Draw!, Urban Sound, and Motion Sense), the following steps are repeated:

Split the data into training, validation, and test sets, as described in Table 2.

Compute $S_{k}(X_{i})$, the signature truncated at order $k$, for every sample $i$. This results in training, validation and test sets of the form

Fit the algorithm on the training data. Validation data is used when the algorithm chosen is the linear neural network or XGboost, to adapt the learning rate and to implement early stopping, respectively.

Compute the accuracy, defined by (10), on the test set.

The maximal value considered for the truncation order, denoted by $K$, is fixed so that the number of features are computationally reasonable. The meaning of “reasonable” depends on each dataset, as they have a different number of samples and classes. For Quick, Draw!, we will consider up to $10^{5}$ samples, for Urban Sound $10^{4}$ and for Motion Sense $5\times 10^{5}$.

The results of this procedure are plotted in Figures 9, 10 and 11, which correspond respectively to the Quick, Draw!, Urban Sound and Motion Sense datasets. A first observation is that some embeddings, namely the time and lead-lag, seem consistently better, whatever the algorithm and the data used. It suggests that this performance is due to the intrinsic theoretical properties of signatures and embeddings, not to domain-specific characteristics. It is particularly remarkable as the dimension of input streams is different from one dataset to another.

On the other hand, the best embedding is the lead-lag path (green curve), followed closely by the time path (brown curve). The difference between these two embeddings is again most important for the Urban Sound dataset. For the Quick, Draw! data, stroke paths have intermediate results, better than the linear path but still worse than the time and lead-lag paths. Yet stroke paths are the only embeddings in which new information, about pen jumps, is included. It is surprising how little impact this information seems to have on prediction accuracy. Note that in all of these cases, the uniqueness of the signature is ensured so it cannot explain the performance differences.

Good performance of the lead-lag path has already been noticed in the literature. However, up to our knowledge, there are few theoretical results. Still, Flint et al. 2016 have considered a discretely sampled input path $X$, assumed to be a continuous semimartingale, and have studied convergence results of its associated lead-lag path, called Hoff process, when sampling frequency increases. Thus, a lot of questions remain open concerning the statistical performance of the time and lead-lag embeddings, with, to our knowledge, no theoretical result in classification or regression frameworks.

To conclude this section, the take-home message is that using the lead-lag embedding seems to be the best choice, regardless of the data and algorithm used. It does not cost much computationally and can drastically improve prediction accuracy. Moreover, the linear and stroke paths yield surprisingly poor results, despite their frequent use in the literature.

3 Running times

To conclude this study on embeddings, this section presents some results on the computational complexity of the different embeddings and truncation order. In Figure 12, the running times for computing signature features and fitting a random forest are shown as functions of the truncation order for various embeddings. The experiments were run on 32 Intel Xeon E5-4660 cores and parallelized with the multiprocessing Python package. The most expensive embedding is the lead-lag, which is not surprising as it doubles the path dimension. Moreover, increasing the truncation order increases exponentially the running time. For Quick, Draw! the running time is of the order of 10-100 seconds, for Motion Sense of the order of 100 seconds and for Urban Sound around 1000 seconds. This is directly linked to the length of the series: the longer the series, the more expensive it is to compute signatures.

In Figure 13 is presented the number of input features for each combination of embedding and truncation order. This is proportional to the memory needed to run each experiment. As given by equation (3), it is clear that the storage cost increases exponentially with the truncation order, which is the main limitation of the signature method.

Simulation study of autoregressive processes

First, Figure 15 shows the results of the same study on the embeddings performance as in Section 4 for different AR(1) processes. The parameter $\phi_{1}$ is equal to $-0.9$, $-1$ and $0.5$, in order to obtain both stationary and nonstationary models. The performance of the different embeddings is plotted against a range of truncation orders. The metric is the $L_{2}$ error on a test set, which is defined by

with the same notations as (10). Contrary to Figures 9, 10 and 11 which use as metric the accuracy, the smaller $S$ the better the prediction. It is clear in Figure 15 that the time and lead-lag embeddings have the smallest errors, confirming the findings of the previous section on real-world datasets. Moreover, the figures are very similar for stationary or nonstationary series, and different strength of time dependence. This shows the generality of the signature method, which does not require strong assumptions on the law of the underlying process.

A natural question is whether the lag parameter is linked to the dependence in the time series. To tackle this issue, Figure 16 shows a boxplot of the $L_{2}$ error as a function of the lag for 3 different values of $p$. The parameters of model (11) are set to the following values: for $p=1$, $\phi_{1}=-0.9$; for $p=3$, $\phi_{1}=\phi_{2}=0$ and $\phi_{3}=-0.9$; for $p=8$, $\phi_{1}=\dots=\phi_{7}=0$ and $\phi_{8}=-0.9$. For each lag, the best truncation order is selected with a validation set, that is the truncation order is chosen to be the one achieving the lowest error on a validation set. The procedure is then evaluated on another test set. This procedure is iterated 20 times to obtain estimates of the variability of the error. Figure 16 shows that when $p$ increases, the best lag increases. For $p=1$, all lags seem to achieve similar errors, for $p=3$, there is an error jump between a lag of 1 and a lag of 2, and for $p=8$ the error decreases towards the best lag of 6. This is strong evidence of the link between the time dependence and the lag parameter.

Finally, Figures 17 and 18 investigate the link between the prediction error, the sample size and the truncation order of the signature for an AR(3) process with the same parameters as before. We choose a lead-lag embedding with a lag of 2, as suggested by Figure 16. For each sample size and truncation order, the model is fitted and evaluated 20 times. In Figure 17, the truncation order is chosen as the minimizer of the error on a validation set. Then, a boxplot of the errors is plotted as a function of the sample size. Both the error and its variance decrease fast when the sample size increases.

On the other hand, Figure 18 shows error boxplots as a function of the truncation order $K$, for various sample sizes. When the sample size is large enough, typically larger than 200, a bias-variance tradeoff can be observed: the error first decreases with the truncation order until a minimum is reached and then the error and its variance increase fast because the number of covariates is too large compared to the number of observations. It is interesting to see that the best truncation order increases when the sample size increases but stabilizes at $K=4$. Moreover, the error variance is very large for a sample size of 10 but stays reasonably small for larger sample sizes.

Signature domain and performance

As discussed in Section 2.2, several authors do not compute signatures on the whole time interval but use instead a partition of $$. The rationale for this division is to describe the path by a sequence of truncated signatures, rather than one signature computed over the whole domain. Therefore, it is not surprising that this approach is typically used in combination with recurrent neural networks. In this section, it is investigated whether signatures computed on a sub-interval contain some local information not present in signatures of the whole path. To this end, following Wilson-Nunn et al. 2018, a dyadic partition of$$ is considered, and defined by (7):

where $q$ is the dyadic order. For different values of $q$, signature coefficients are computed on each interval $\big[(j-1)2^{-q},j2^{-q}\big]$ of the dyadic partition. Therefore, for each input path $X_{i}$, a collection of signature vectors is obtained, which is then stacked into one large vector. This vector is then the input of a learning algorithm, and the prediction accuracy curves of different dyadic orders are compared. The time embedding with the linear neural network described in Section 4.2 is used. This process is summarized below.

Split the data into training, validation and test sets.

For $j=1,\dots,2^{q}$, compute the signature truncated at order $k$ on $\big[(j-1)2^{-q},j2^{-q}\big]$, denoted by

where $X_{i}$ is the time embedding of sample $\mathbf{x_{i}}$. Repeat this over all training samples.

Fit a linear neural network with this data as features.

The results of this procedure are shown in Figure 19. First, it is clear that a dyadic order of $0$, which corresponds to computing the signature on the whole interval $$, yields the best results. Indeed, the curve is always above the others for the Quick, Draw! and Motion Sense datasets. This is less obvious for the Urban Sound dataset, as the curve$q=0$is really close to the one corresponding to a dyadic order of$1$. However, it still achieves a better result for most truncation orders$k$. This difference between datasets can be linked to their length: it seems that the longer the series, the better the accuracy of thin dyadic partitions. Indeed, high dyadic orders perform best for the Urban Sound dataset, which has an average of 170 000 sampled time points (see Table 2), whereas it is clear that each new dyadic split decreases accuracy for the Quick, Draw! data, which has an average of 44 sampled points.

In a nutshell, little local information seems to be lost when the signature of the whole path is computed. However, it may be worth considering partitions of the path for long streams.

2 Performance of the signature

The message of previous sections is that the lead-lag embedding is the most appropriate in a learning context and that signatures should be computed over the whole path domain. As a natural continuation, the effect of the algorithm is now examined more closely and the prediction scores are compared to the literature. It turns out that the signature combined with a lead-lag embedding has an excellent representation power, to the extent that it achieves prediction scores close to state-of-the-art methods, without using any domain-specific knowledge.

Before starting the comparison, it is worth pointing out that the lead-lag embedding has a hyperparameter that has not yet be tuned, which is the number of lags. It is now selected with the same approach as in previous sections: for each lag, the test accuracy is plotted against the number of features for various truncation orders. The lag which gives a curve above the others is selected. Figure 20 highlights that curves overlap for the Motion Sense and Urban sound cases, therefore, when there is a doubt on which curve is above, the smallest lag is picked. For the Quick, Draw! and Motion Sense datasets, the best lag is then 1, whereas it is 5 for the Urban Sound dataset.

Finally, the truncation order is selected with a validation set: the truncation order achieving the highest accuracy on a validation set is picked. In general, the truncation order and the number of lags could also be selected with cross-validation.

The Motion Sense and Urban Sound datasets do not require a lot of computational resources, as they have a reasonable size (a few hundred samples for Motion Sense and thousand for Urban Sound—see Table 2) and a small number of classes. On the other hand, the Quick, Draw! recognition task, which is comprised of 340 classes, is more involved and requires an elaborate algorithm as well as a significant number of training samples. Therefore, more data will be used than in the previous sections: $12\,185\,600$ training samples and $87\,040$ validation samples.

For the Quick, Draw! dataset, our results are compared to a Kaggle competition (Kaggle.com 2018). In this competition, 1 316 teams competed for a prize of 25 000$. State-of-the-art deep convolutional networks, such as MobileNet or ResNet, trained with several millions of samples, were among the best competitors. Teams on the podium used ensembles of such networks. These winning methods require a lot of computing resources and are specific to images. The metric used in the competition was the mean average precision, defined as follows. Denoting by$\{y_{1},\dots,y_{n_{\text{test}}}\}$ the test set labels, three ranked predictions are made for each sample, denoted by

where $\hat{y}^{1}_{i}$ is the class with the largest probability, $\hat{y}^{2}_{i}$ the second largest, and so on. Then, the mean average precision at rank 3 is defined by

Mean average precision is computed by the competition platform on 91% of a test set of 112 200 samples. The small neural network used in Section 4.2 is enhanced by using ReLU activation functions and adding three hidden layers with 256 nodes. The network is trained during 300 epochs with an Adam optimizer. At each epoch, it is trained on 609 280 samples randomly selected among the 12 185 600 training samples. The best team obtains a $MAP_{3}$ of 95% whereas this small network combined with signature features truncated at order 6 already achieves 54%. The winners use an ensemble of several dozens of deep neural networks, trained on 49 million samples. This kind of architectures requires considerably more computational capacities than ours.

For the Urban Sound dataset, state-of-the-art results are obtained by Ye et al. 2017. The authors combine feature extraction with a mixture of expert models and achieve 77.36 % accuracy, defined by (10). The feature extraction step is specific to sound data and is based on several ingredients, such as whitened spectrogram, dictionary learning, soft-thresholding, recurrence quantification analysis, and so on. These crafting operations make use of a lot of domain-specific knowledge and cannot be extended easily to other applications. On the other hand, it is clear from Figure 10 that a random forest classifier performs well with signature features. Therefore, its hyperparameters are tuned with a lead-lag embedding, a lag of 5, and a signature truncated at order 5. An accuracy of 70 % is obtained with 460 trees with a maximum depth of 30 and in which 500 random features are considered at each split.

Finally, Malekzadeh et al. 2019 tackle the problem of mobile sensor data anonymization. They build a deep neural network architecture that preserves user privacy but still detects the activity performed. The architecture is built on autoencoders combined with a multi-objective loss function. There is a trade-off between activity recognition and privacy but good activity recognition results are achieved. Performance of the classifier of Malekzadeh et al. 2019 is measured with the average $F_{1}$ score, defined as follows. Assume there are $C$ different classes, and denote by $(y_{1},\dots,y_{n_{\text{test}}})$ the test labels, and by $(\hat{y}_{1},\dots,\hat{y}_{n_{\text{test}}})$ the predicted ones. Then, the $F_{1}$ score is defined by

Malekzadeh et al. 2019 report an average $F_{1}$ score above 92%, while the signature truncated at order 3 and combined with a XGBoost classifier achieves a $F_{1}$ score of 93.5%. These two scores are close, but the signature approach is computationally much less demanding.

Despite not being tuned to a specific application, the combination signature + generic algorithm achieves results close to the state-of-the-art in several domains, while requiring few computing resources and no domain-specific knowledge. Indeed, it takes approximately 52 seconds to compute the signature at order 3 of 68 000 Quick, Draw! samples on one core of a laptop, which results in $0.0008$ second per sample. Besides, signature computations can be parallelized, making the approach scalable to big datasets. Lastly, the signature method achieves its best results for the high dimensional Motion Sense dataset, which suggests that it is especially relevant for multidimensional streams.

Conclusion

The signature method is a generic way of creating a feature set for sequential data and has recently caught the machine learning community’s attention. Indeed, it yields results competitive with state-of-the-art methods, while being generic, computationally efficient, and able to handle multidimensional series. One of its appealing properties is that it captures geometric properties of the process underlying the data and does not depend on a specific basis. In this paper, its use in a learning context, and several of its successful applications have been reviewed. The use of signatures relies on representing discretely sampled data as continuous paths, a mechanism called embedding. In the literature, authors use various embeddings, without any systematic comparison. We have compared different common embeddings and concluded that the lead-lag seems to be systematically better, whatever the algorithm or dataset used. Moreover, we have pointed out that the signature of the whole path appears to contain as much information as the signature of subpaths, therefore encoding both global and local properties of the input stream.

Our study is a first step towards understanding how signature features can be used in statistics, and a lot of issues remain open, both practical and theoretical. First, it would be of great interest to understand the theoretical statistical properties of embeddings, in particular, to explain the good performance of the lead-lag path. Moreover, in Section 2.2, we have seen that signature features may be combined with feedforward, recurrent, or convolutional neural networks. For each of these architectures, the point of view on signature features is different: they are considered respectively as a vector, a temporal process, or an image. A more detailed understanding of these representations would be valuable. Finally, it could be worth investigating the robustness of the signature method when the truncation order becomes large. Indeed, Figures 9, 11, and 10 suggest that the signature may be robust to dimension: the accuracy curves do not decrease when the number of features becomes large, even when a nearest neighbor algorithm is used with more than a hundred thousand features. This phenomenon may deserve a more in-depth study.

Acknowledgements

This work was supported by a grant from Région Ile-de-France. I thank Gérard Biau (Sorbonne Université), Benoît Cadre (Université Rennes 2) and Terry Lyons (Oxford University) for stimulating discussions and insightful suggestions; and Patrick Kidger (Oxford University) for his thorough proofreading of the manuscript. I also thank the Editor, the Associate Editor, and two anonymous referees for their careful reading of the paper and constructive comments, which led to a substantial improvement of the article.

References