Understanding Any Time Series Classifier with a Subsequence-based Explainer

A time series signal is defined as follows:

A set of one or more time series signals forms a time series:

When \(d=1\), the time series is univariate, while if \(d\gt 1\) the series is multivariate. A Time Series Classification (TSC) dataset is a set of time series with a vector of labels (or classes) attached.

For a dataset \(\mathcal {D}\) containing c classes, \(y_i\) can take c different values. When \(c=2\), \(\mathcal {D}\) is a binary classification dataset, while if \(c\gt 2\), then \(\mathcal {D}\) is a multi-class classification dataset. In order to ensure clarity and consistency in notation, we have adopted a tensor-like notation based on [52]. Lowercase letters are used to denote single scalar observations (e.g., x), while bold lowercase letters are used to denote vectors and individual signals of time series (e.g., \({\bf x}\)). Capital letters are used to denote matrices and time series instances (e.g., X), and Euler script letters are used to denote tensors and time series datasets (e.g., \(\mathcal {X}\)). To indicate a specific observation in a time series dataset, we use the notation \(x_{i,j,k}\), where i denotes the \(i^{th}\) multivariate time series in the dataset, j denotes the \(j^{th}\) time-step, and k denotes the \(k^{th}\) signal of the time series. Although time series do not always form a proper tensor since different channels can have different numbers of observations, we use a unique symbol m to denote the length of the time series for simplicity of notation. When indexes are not relevant, we omit them for better readability. We can now define the TSC problem as:

The resulting TSC function f takes as input a time series X and returns \({\hat{y}}\) according to what f learned, i.e., \({\hat{y}}= f(X)\). In general, \({\hat{y}}\) can either be a discrete label or the probability of X belonging to a specific class. We use \(f(\mathcal {X}) = {\hat{{\bf y}}}\) as a shorthand for \(\lbrace f(X) \;|\; X \in \mathcal {X}\rbrace = {\hat{{\bf y}}}\). Typically, the classifier f can be queried at will.

A common way to build classifiers in the time series domain is to use subsequences. In simple terms, a subsequence is a continuous sample of observations from a time series. Formally:

Subsequence-based classification approaches search for patterns that better discriminate the dataset labels, i.e., they try to find those subsequences that are most dissimilar between instances belonging to different classes. The two most common kinds of subsequences used for this purpose are shapelets [95] and symbolic subsequences [61].

Shapelets . Typically, shapelet-based methods extract a set containing p-most discriminative shapelets by minimizing an information gain-like metric. Once the most discriminative shapelets are found, the Shapelet Transform [62] can be applied in order to transform the time series dataset into a simplified representation. Formally:

As a note, the sliding-window distance is method-dependent, and it is calculated w.r.t. the signal the shapelet was extracted from. In practice, the shapelet transform extracts the p-most discriminative shapelets from a time series dataset and returns a new representation of the data where the attributes represent the distances between each time series and the p shapelets. Hence, any classification algorithm can be used, potentially increasing the accuracy while reducing training time.

Symbolic Subsequences . Given that exhaustive subsequences search is computationally expensive, a common approach is to first transform the time series into a simplified representation by using SAX. The Symbolic Aggregate approXimation (SAX) algorithm [61] transforms time series into sequences of strings. For each time series signal \({\bf x} = \lbrace x_{1}, \dots , x_{m}\rbrace \in \mathbb {R}^{m}\), SAX uses the Piecewise Aggregate Approximation (PAA) [46] to split it into w equally sized intervals and averages the values of each interval. Then, the time series signal is discretized using a finite set of symbols, i.e., an alphabet \(\mathbb {A}\). Formally, \(\text{SAX}({\bf x}) = \tilde{{\bf x}} = \lbrace \tilde{x}_{1}, \dots , \tilde{x}_{w}\rbrace \in \mathbb {A}^{w}\), with \(|\mathbb {A}| \gt 1\) being the number of symbols in the chosen alphabet. This approximation reduces running time while denoising the time series. Once the time series are converted to a symbolic representation, the most discriminative symbolic subsequences can be found using different techniques. Similar to the Shapelet Transform, the dataset can be transformed into a new representation, having as features the extracted subsequences and as values 0 or 1 depending on the absence or presence of subsequences inside the time series. Formally:

For interpretability purposes, the symbolic subsequences can be easily mapped back to the original segments of the time series. We emphasize that we use the same symbol \({\varsigma }\) to denote the symbolic subsequence transform and the shapelet transform, given that they both convert a time series dataset into a tabular representation with subsequences as features.

Given a not interpretable, i.e., black-box, time series classifier b and a time series X classified by b, i.e., \(b(X) = {\hat{y}}\), our aim is to provide an explanation e for the decision \(b(X)={\hat{y}}\). More formally:

To build a complete, human-interpretable explanation in the time series domain, we consider three kinds of explanations: saliency maps, examples, and decision rules.

Saliency Maps. . Saliency maps are explanations that highlight the contribution of each feature for the classification [36]. Formally, for time series:

In practice, the saliency map assigns an importance score to each observation in X depending on its contribution to the classification output.

Examples. . Examples, also called example-based (or instance-based) explanations, use whole time series objects to add interpretability to classification models by providing a comparison of the instance to explain with a salient time series. The two main different types of examples are exemplar and counterexemplar instances [34, 37]. Exemplars, also called prototypes, are time series that exemplify the main characteristics that influence the classifier’s decision. Formally:

The explanation is obtained by comparing the instance X for which we have the decision \(f(X)\) with the exemplar \({X_=}\) that represents it. Representation is usually formalized with a notion of similarity. On the other hand, counterexemplars, also called counterfactual instances, are very similar w.r.t. the instance to explain but are classified differently. The explanation is achieved by comparing the minimal differences in shape that lead to a different classification outcome. Formally:

Minimality usually refers to a distance metric, while plausibility is assessed by checking if the instance is not simply an adversarial example [29] and is semantically coherent with the dataset.

Decision Rules . Explainable AI methods increasingly simplify explanations to increase user trust by ensuring they identify cause-effect relations between events [12]. In this sense, rule-based explanations are arguably the most interpretable from a human standpoint as they allow the user to understand the reason behind a classifier’s decision in terms of if-then statements. Formally:

Decision rules are generated by rule-based classifiers or can be inferred by analyzing the splits of decision trees. They can be naively applied to time series by considering observations as features. In this case, a condition is in the form \((j,k) \in [v_{\text{low}}, v_{\text{up}}]\), where \((j,k)\) are the indexes identifying the \(j^{th}\) time-step of the \(k^{th}\) signal, and \(v_{\text{low}}, v_{\text{up}} \in \mathbb {\bar{R}}\) are the lower and upper bounds on the observation value. A time series \(X_i\) is covered by the rule if every condition in p is true. A local rule-based model rm can be used for explaining the black-box prediction for X if it imitates sufficiently well its behavior, i.e., if \(rm(X) = b(X)\) and also for every \(X^{\prime }\) in the neighborhood of X, \(rm(X^{\prime }) = b(X^{\prime })\). The definition of the neighborhood is method-dependent. A rule that directly explains the prediction of a black-box is called a factual rule. In contrast, the rules obtained by minimally removing or adding conditions in the factual rule premises are called counterfactual rules. Counterfactual rules are extremely useful for what-if analysis because they allow understanding of the minimal variation that results in a different classification by the black-box. The most notable rule-based XAI approaches for tabular data are Anchor [79], a model-agnostic method explaining the behavior black-boxes with high precision factual rules, and lore [35], which generates the local neighborhood via a genetic algorithm, trains a decision tree to then extract factual and counterfactual decision rules. The first version of lasts [37] extends lore to univariate time series data, generating rules that explain the decision of the black-box in terms of logical conditions based on time series subsequences. In this case, the conditions in the premise refer to subsequences instead of single time series observations, and the feature values represent the presence/absence of subsequences inside the time series. We build upon this to extend the approach to multivariate time series data.

A standard autoencoder (AE) [40] is a type of neural network trained for learning a representation that reduces the dimensionality from \(m \times d\) to k and captures non-linear relationships. An encoder \({g}: \mathbb {R}^{m \times d} {\rightarrow } \mathbb {R}^{q}\), and a decoder \({h}: \mathbb {R}^{q}{\rightarrow } \mathbb {R}^{m \times d}\) are simultaneously trained with the objective of minimizing the reconstruction loss. Starting from the encoding \({\bf z} = {g}(X)\), the autoencoder tries to reconstruct a representation as close as possible to its original input \(X \simeq {h}({\bf z}) = \hat{X}\). Autoencoders learn to encode their input in a latent representation, which is usually of smaller dimensionality and, therefore, simpler and easier to deal with, w.r.t. the original input. The latent space can also be used to sample synthetic instances that can be decoded into completely new and unseen records. In order to use autoencoders as generators, it is useful to have an encoder with a specific latent distribution from which to sample. In this sense, the two leading solutions used in the literature are Variational Autoencoders (VAEs) [51], which learn the parameters of a latent distribution, usually the mean and standard deviation of a Gaussian, or Adversarial Autoencoders (AAEs) [66], that use a GAN inspired approach [28], adding a discriminator in the network architecture, trained to discriminate the latent distribution, constraining it to the desired form.

In order to support the generation of a good explanation, the autoencoder used in our proposal needs to have some desirable properties:

The first property can be ensured by using VAEs or AAEs. These are superior to traditional AEs as generative models because they allow sampling from a known latent distribution. Usually, a VAE is easier to train w.r.t. an AAE, given the former has to optimize only one loss function, while the latter also has to consider the discriminator [34, 66]. For this reason, we adopt VAEs by default as autoencoders for the proposed approach. The second property can be ensured by checking the sampled instances with a discriminator or by using specific sampling techniques that minimize the distribution mismatch [2] (Section 3.4). The third and fourth properties can be achieved with hyperparameter tuning. In our setting, the reconstruction error is measured in terms of mean squared error between the original instances and the reconstructed ones, \(\text{rec}_\text{mse} = \frac{1}{n} \sum _{i=1}^{n}(X_i-\hat{X_i})^2\). The reconstruction accuracy is the percentage of instances in \(\mathcal {X}\) that are correctly classified by the black-box after being reconstructed by the autoencoder. Formally, \(\text{rec}_\text{acc} = \frac{1}{n} \sum _{i=1}^{n} \mathbb {1}_{b(X_i) = b(\hat{X_i})}\), where \(\mathbb {1}\) is a function which outputs 1 if the condition in the subscript is true and 0 otherwise.

Neighborhood generation, being a central issue in the definition of local surrogates, is an increasingly studied topic in the domain of XAI [37, 64, 78]. Local surrogates mimic the local decision boundary of a black-box for the single instance they are tasked to explain. In order to be able to imitate the black-box, a representative neighborhood of the instance to explain has to be defined. It should be composed of (i) similar, i.e., spatially close, instances having the same label as the instance to explain (prototypes, exemplars), but it also should contain (ii) close instances having a different label (counterfactual instances, counterexemplars, distractors). Post-hoc interpretability approaches often rely on perturbations of the input data that query the black-box to understand how its prediction changes. Two of the most notable examples are lime [78] and shap [64]; however, the perturbation methods used by these models do not ensure the generation of realistic data. Moreover, even if these approaches can be applied to time series data with some modifications, they are not explicitly thought for it [34, 37, 54].

For time series data, it is often hard to generate a meaningful neighborhood in the manifest space by directly perturbing the original instances because the risk of generating unrealistic adversarial examples is relatively high. This is the main reason why generative models such as VAEs and AAEs are increasingly used in this field [7, 34, 37, 43, 54, 74]. In fact, these models are usually trained to produce a specific prior distribution in the latent space, typically a Gaussian. After the training phase, they can be used as generators by sampling new instances or perturbing existing ones. Perturbation techniques can greatly vary: from the usage of different sampling or searching algorithms [74], to genetic approaches [34, 37], to gradient descent-based methods [7, 43]. Even if the technical aspects of these methods are very different, they all rely on autoencoders to generate a suitable latent space to explore and analyze.

In [93], it is shown that latent space operations commonly used in the literature, such as instance interpolation and vicinity sampling, induce a so-called distribution mismatch between the outputs and the prior distribution the model was trained on. This is a delicate issue, given that decoders and generators are usually trained on fixed priors and thus assume that their inputs will have statistical properties that align with their distributions. For this reason, in [2] is proposed a Gaussian-matched neighborhood generation function which avoids the sampling from latent space locations that are highly unlikely given the prior distribution of the autoencoder, and allows the creation of more consistent, i.e., higher quality synthetic instances. Given a latent space vector \({\bf z}\in \mathbb {R}^{{q}}\) in which each scalar entry is independently sampled from a standard Gaussian distribution, i.e., \(\forall z \in {\bf z},\; z\sim \mathcal {N}(0,1)\), a standard Gaussian (ga) vicinity sampling is defined as \({\bf z}_{{\rm\small GA}} = {\bf z}+\theta {\bf u}\) with \(\theta\) being a scaling factor and \({\bf u}\in \mathbb {R}^{{q}}\) a randomly sampled normal Gaussian vector. The Gaussian-matched (gm) operation is defined as \({\bf z}_{{\rm\small GM}} = {\bf z}_{{\rm\small GA}} / \sqrt {1+\theta ^2}\). In [2] it is shown that this approach is guaranteed to produce samples coming from a standard Gaussian distribution, i.e., samples that are indistinguishable from randomly sampled instances from that distribution. The authors show that this property is independent of the number of latent dimensions, and it is proven to work even in high-dimensional latent spaces. Incorporating this operation into our framework helps to reduce the distribution mismatch and generate a synthetic neighborhood that is more consistent, resulting in samples that seem as if they were drawn from the original training set.

The time series X is passed to the encoder part of the VAE that compresses it into a so-called latent representation \({\bf z}={g}(X)\). The time series in Figure 3 is used as a running example. In our case, X originally has 128 observations, and it is compressed in a bidimensional vector, \({\bf z}=[1.344, -2.005]\), depicted in Figure 4 (left). The latent vector \({\bf z}\) can also be passed to the decoder to reconstruct the original time series. Both X and its reconstructed version, \(\hat{X}={h}({\bf z})\), are depicted in Figure 3 (right). The original time series is much more noisy w.r.t. its reconstructed version. This suggests that the autoencoder is able to capture the most relevant features of the time series for the classification, i.e., its general shape, discarding the random noise.

The second step of lasts is to find the closest counterfactual instance w.r.t. \({\bf z}\). To perform this search, lasts adopts an algorithm that iteratively samples latent instances around \({\bf z}\), decodes them, and checks their label using b. The counterfactual search algorithm (cfs) is reported in Algorithm 2. cfs uses a function \({\it sample}\) to sample a distribution of synthetic instances around \({\bf z}\). The function sample can theoretically be any sampling function, ranging from a pure random approach like in lime [78] to a genetic algorithm maximizing a fitness function like in lore [33]. A good generation function is fundamental to create consistentsynthetic instances by avoiding the aforementioned distribution mismatch [2], i.e., the sampling from locations in the latent space that are highly unlikely given the prior distribution of the autoencoder. We emphasize that we employ the term “consistent” as mentioned in [2], i.e., as a synonym for probable, likely, or, in other words, coherent with the autoencoder’s underlying distribution. Our aim is to convey that the generated samples possess a similar distribution to that of the original training data. The first version of lasts presented in [37] uses the genetic approach of lore; in the version presented in this paper, lasts adopts by default a procedure inspired by the growing sphere algorithm proposed in [55].

cfs depends on a threshold \(\theta\) that decides the amount of space around \({\bf z}\) that the generated neighborhood is going to occupy. This threshold can be the radius for a spherical random uniform, as in the original growing sphere [55], or a scaling factor for a Gaussian distribution, as detailed in Section 3.4. The goal of cfs is to systematically explore the latent space in order to find the instance closest to \({\bf z}\) belonging to a different class, i.e., the closest counterfactual. The counterfactual search is performed by iteratively generating a neighborhood Z around \({\bf z}\), using the aforementioned sample function, and by checking if there is at least one counterfactual among the generated instances. The presence of counterfactuals indicates that the sampled neighborhood is still crossing the decision boundary; therefore, the threshold is halved, and the procedure loops until the sampling function does not generate any counterfactual. We highlight that in this setting, the generation happens in the latent space, while the presence of counterfactuals is checked in the time series domain by decompressing the latent instances through \({h}\). At each step of the iteration, the cfs algorithm stores all counterfactuals generated and, once out of the loop, it selects the closest one to \({\bf z}\), i.e., \({{\bf z}_\ne }\), as the best counterfactual. In the running example, the closest counterfactual instance can be viewed in its latent form in Figure 4 (left).

The third step we detail is the generation of the neighborhood around \({{\bf z}_\ne }\) using the function \(\mathit {neighgen}\). In principle, \(\mathit {neighgen}\) can be a different sampling function w.r.t. the sample function used for the search. In Section 5, we experiment with different combinations of the two. The sampling is performed extremely close to the black-box decision boundary; therefore, this synthetic neighborhood contains, by construction, both a set of time series having the same class as \({\bf z}\) and a set of time series having a different class, i.e., exemplars, \({Z_=}\), and counterexemplars, \({Z_\ne }\) (Figure 4, left). The number of distinct counterexemplar classes is influenced by the final threshold, and increases as \(\theta\) increases. Note that the sampling is performed in the latent space, while the labels are retrieved by decoding the latent neighborhood and then applying the black-box function to obtain its predictions. For ease of viewing, we depict all the counterexemplar instances with the red color, independently of their specific class. Z is decoded and classified by the black-box into \({\hat{\mathcal {X}}_=}\) and \({\hat{\mathcal {X}}_\ne }\). These instances represent the example-based explanation, showing how the decision of the black-box changes depending on the shape of the time series (Figure 4, right). Moreover, by computing the absolute difference between the decoded closest counterfactual \({\hat{X}_\ne }= b({{\bf z}_\ne })\), and \({\hat{X}}=b({\bf z})\), we can discover the time points for which a change in the values is more likely to modify the decision of the black-box. Thus, the saliency map is defined as \(\Phi = |{\hat{X}}- {\hat{X}_\ne }\)|.

Figure 4 (left) shows the latent neighborhood Z highlighting in green the exemplar instances labeled as \({\it bell}\), and in red the counterexemplar instances labeled with a different class value, in this case cylinder. In gray is depicted a standard bivariate normal distribution, highlighting the probability of hidden vectors in the latent space; formally \(\mathcal {N}_2(\mathbf {\mu }, \Sigma)\) with \(\mathbf {\mu } = [0, 0]^T\) and \(\Sigma = I_2\) where \(I_2\) is a 2-dimensional identity matrix. For ease of viewing, we provide an explanation for an instance at the edge of the distribution. The corresponding instances in the manifest space can be viewed in the same figure to the right. The synthetic neighborhood sampled in the latent space perfectly summarizes the separation of the different class values, unveiling a local decision boundary that is easy to detect even with a simple classifier. Both \({Z_=}\) and \({Z_\ne }\), forming the neighborhood, densely surround \({{\bf z}_\ne }\), helping to capture the black-box behavior locally around the closest decision boundary to \({\bf z}\). Note that \({\bf z}\) remains at the edge of the generated neighborhood by design because our area of interest is the decision boundary and not the instance to explain itself.

Given the decoded local neighborhood \({\hat{\mathcal {X}}}\) and \({\hat{{\bf y}}}= b({\hat{\mathcal {X}}})\), lasts extracts a set of subsequences S and performs the subsequence transform \({\varsigma ^*}\), encoding time series into a space of presence/absence of subsequences (line 7, Algorithm 1), as detailed in Section 3.1. Subsequences can be of many kinds; in Section 5, we test both shapelet-based and SAX-based subsequences. Figure 5 presents an example of SAX-based subsequences and a heatmap depiction of the transformed dataset. Each column represents a time series in \({\hat{\mathcal {X}}}\). For each row, a colored cell indicates the presence of the corresponding subsequence, while a white cell indicates its absence. The colored cells, green and red, indicate the SAX subsequences contained by the exemplar and counterexemplar time series, respectively. In some cases, we can observe a sort of complementarity. In other cases, a time series can simultaneously contain a combination of subsequences describing exemplars and counterexemplars.

Given T and \({\hat{{\bf y}}}\), lasts trains a subsequence-based decision tree classifier \(\mathit {dt}\) that allows to identify subsequence-based factual and counterfactual rules \({r_=}, {{r_\ne }}\) (lines 8 - 9, Algorithm 1). lasts adopts decision trees because factual and counterfactual decision rules can be naturally derived by following the root-leaf paths [33]. Figure 6 (bottom) reports the explanation rules for our example: \({r_=}= \lbrace s_{327} \in X \rbrace \rightarrow {\it bell}\), \({{r_\ne }} = \lbrace s_{327} \not\in X \rbrace \rightarrow \lnot {\it bell}\). The visual representation of the rules shows the position of the subsequences that must be contained and those that must not be contained at their best alignment with X. Looking at the rules, a user can truly understand the reasons for the classification. In this case, the presence of an increasing pattern in the time series differentiates between a bell instance and instances of other classes. In general, rules can be longer and include an arbitrary number of contained and non-contained subsequences, depending on the complexity of the classification task and the resulting surrogate tree.

In summary, lasts explains the prediction of a black-box b for a time series instance X. The final explanation, \(e = \lbrace {\Phi }, ({\hat{\mathcal {X}}_=}, {\hat{\mathcal {X}}_\ne }), ({r_=}, {{r_\ne }}) \rbrace\), is composed of three parts. First, a saliency map, \({\Phi }\), highlights the most important timesteps of the time series, i.e., those timesteps that, if changed, would bring the prediction toward a different class. Second, a neighborhood composed of exemplar and counterexemplar instances. Exemplars, \({\hat{\mathcal {X}}_=}\), are instances close to X in the latent space, with the same label as X. They are prototypes, showing the main characteristics of a specific class. On the other hand, counterexemplars, \({\hat{\mathcal {X}}_\ne }\), are instances close to X in the latent space but classified by the black-box differently and can help the user understand how the prediction of the black-box changes by changing the shape of the time series consistently. Finally, the factual, \({r_=}\), and counterfactual, \({{r_\ne }}\), rules logically explain the prediction of the black-box, both in direct and contrastive ways, showing the minimum change in contained/not contained subsequences to modify the black-box prediction. The complete explanation for our running example is presented in Figure 6.

In Table 1, all the information about the datasets used for evaluating our approach is reported. The training set \(\mathcal {X}_{\text{train}}\) is used both to train the black-box and the autoencoders. In principle, the dataset used for the black-box training, \(\mathcal {X}_{\text{bb}}\), and the dataset used for the autoencoder training, \(\mathcal {X}_{\text{ae}}\), can be different. However, due to the small dimensionality of some datasets, we set \(\mathcal {X}_{\text{bb}} = \mathcal {X}_{\text{ae}} = \mathcal {X}_{\text{train}}\) for all the datasets, with the exception CBF and CBM which are synthetic, and their instances can be sampled at will. The test set \(\mathcal {X}_{\text{test}}\) is used to benchmark both autoencoders and black-box models, and to sample \(\mathcal {X}_{\text{exp}}\), which is composed of a maximum of 50 instances to explain and evaluate.

To test different framework alternatives, we train and explain a ResNet [39] (res) implemented in keras according to [24], and a k-Nearest Neighbor [85] (knn) baseline as implemented by scikit-learn. Once the best framework setup is found, we choose as a black-box to explain Rocket [17], as implemented by sktime, using the default parameters for the transform and a RidgeClassifierCV as classification model. The ResNet comprises three residual blocks, each containing three convolutional layers, a global average pooling layer, and a dense layer. The layers inside each residual block have respectively \(64, 128, 256\) filters, of size \(8, 5, 3\). We train res with a batch size of 16, monitoring the loss, with a patience parameter of 50 epochs. As optimizer, we select Adam, with the default keras parameters: \(\mathit {learning\_rate} = 0.001\), \(\beta _1=0.9\), \(\beta _2=0.999\), minimizing the sparse categorical crossentropy. For knn we use the Euclidean distance and the k parameter is selected via grid-search with \(k \in [1, |\mathcal {X}_{\text{bb}}|]\).

In the following, we specify the implementation details adopted for each module of lasts³.

Variational Autoencoder. . As autoencoder, we adopt a VAE for the reasons discussed in Section 3. For simplicity, we use the same network structure for all the datasets, i.e., a convolutional autoencoder composed of 16 layers, 8 for the encoder, and 8 for the decoder. The number of filters per layer is 8, and the kernel size is set to 3. The latent dimension is chosen on a dataset basis by starting with \({q}=2\) and iteratively building autoencoders with an increasing value of \({q}\) until the accuracy for the reconstructed instances increases and reaches 0.90, but also keeping \({q}\le m/2\). As an activation function, we always use ReLU. As optimizer, we select Adam with default keras parameters: \(\mathit {learning\_rate}=0.001\), \(\beta _1=0.9\), \(\beta _2=0.999\), minimizing the Mean Squared Error (MSE) with the Kullback–Leibler divergence regularization term. The VAEs are trained with a batch size equal to 32, for a maximum of 8,000 epochs, monitoring the validation loss with a patience of 500 epochs before halving the learning rate and a patience of 1,250 epochs before early stopping. We measure the performance of the autoencoders through the reconstruction error between the original and reconstructed time series, in terms of Mean Squared Error (the lower, the better), and in terms of accuracy of the classifiers on the reconstructed time series (the higher, the better). Since the encoding operation in a VAE is stochastic, \({\bf z}\) can vary slightly. Therefore, to improve the stability of the framework, X is encoded 1000 times, and the latent representation \({\bf z}\) is chosen by checking the most similar reconstruction \(h({\bf z})\) to X using the Euclidean distance. Given a set of 1,000 encodings of X named \(Z^{\prime }\), \({\bf z} = \arg \min _{{\bf z}^{\prime }\in Z^{\prime }} \text{dist}({h}({\bf z}^{\prime }), X)\).

Neighborhood Generators . For the counterfactual search and neighborhood generation, we experiment with four alternatives for the neighgen and sample functions of Algorithm 1 and Algorithm 2: (i) Gaussian (ga) that samples a normal distribution around \({\bf z}\) scaling each vector by a factor \(\theta\) as detailed in Section 3.4; (ii) Gaussian-matched (gm) that is a Gaussian-based sampling [2] that uses distribution matching transport maps to minimize the problem of distribution mismatch. As in Gaussian, it samples the distribution around \({\bf z}\) scaling each vector by a factor \(\theta\) as detailed in Section 3.4; (iii) Uniform Sphere (us) that samples a uniform sphere distribution around \({\bf z}\) with radius \(\theta\); and (iv) Matched Uniform (mu) that combines a Gaussian-matched search to find the closest counterfactual, and Uniform Sphere sampling to generate the neighborhood. In the first three cases, the counterfactual search and the neighborhood generation use the same function, i.e., \({\it sample}={\it neighgen}\), while the last approach uses two different ones, i.e., \({\it sample}\ne {\it neighgen}\).

For the counterfactual search, we generate \(n_s = 10,\!000\) instances at each iteration, starting with a threshold \(\theta = 2\). Once the closest counterfactual is found, the neighborhood generation function \(\mathit {neighgen}\) is run with neighborhood size equal to \(N = 500\) latent instances. If \(\theta\) represents a radius, we perform the final sampling using as \(\theta\) the distance between \({\bf z}\) and \({{\bf z}_\ne }\) while if \(\theta\) represents a scaling factor, we take the last \(\theta\) used in the counterfactual search. In the last generation step, we impose a balance between instance labels by oversampling the minority class. Moreover, we compare these sampling functions against the genetic approach adopted in [37]. For the genetic approach, we also generate \(N = 500\) latent instances using the same parameters as in [37], namely 10 generations, normalized Euclidean distance as the genetic fitness function, probability of mutation equal to 0.5, probability of crossover equal to 0.7.

Subsequence-based Surrogates . We implement the function \({\varsigma ^*}\) of Algorithm 1 in two ways in order to test two different strategies to retrieve subsequences: SAX-based and shapelet-based. To extract SAX-based subsequences, we use the SAX-SEQL algorithm illustrated in [57], which extracts the p-most discriminative subsequences S for \({\hat{{\bf y}}}\), and then converts time series in a binary-valued matrix T using the symbolic subsequence transform detailed in Section 3.1. Furthermore, to extract shapelets-based subsequences, we adopt the LTS algorithm described in [30] that learns the p-most discriminative shapelets S with respect to \({\hat{{\bf y}}}\) via gradient descent. The time series dataset is then converted to a simplified representation via the shapelet transform detailed in Section 3.1. In order to have the same representation as in the SAX-subsequences alternative, the distances in T are replaced with binary values using a threshold \(\tau\) such that \(\forall t_{i,j} \in T\) if the distance \(t_{i,j} \lt \tau\), \(t_{i,j}\) is replaced with 1 else with a 0. The distance threshold \(\tau\) is chosen via grid search by testing the accuracy of the decision tree surrogate for every decile of the distribution of all the distances in T. To further simplify the classification task of the local surrogate, we binarize the label vector by considering only the predicted class of the instance to explain as 1 and all the others as 0. Therefore, \(\forall {\hat{y}}_i \in {\hat{{\bf y}}}\), if \({\hat{y}}_i = b(X)\), \({\hat{y}}_i\) is replaced with a 1 else with a 0.

In the following, we analyze the effect of various alternatives in terms of neighborhood generation and subsequence types on the explanation returned by lasts.

We measure the performance in terms of fidelity, aiming to evaluate how good the explanation model is at mimicking the black-box decisions. The fidelity (fid) is simply defined as the accuracy between the prediction of the black-box and that of the explanation model [9]. In our setting, we compare the prediction of the black-box for the explanation dataset, \({\hat{{\bf y}}}= b(\mathcal {X}_\text{exp})\), and \({\hat{{\bf y}}}^{\prime } = \lbrace {\hat{y}}^{\prime } | \forall X \in \mathcal {X}_\text{exp},\; {\hat{y}}^{\prime } = \mathit {dt}({\varsigma ^*}({h}({g}(X))))\rbrace\), where \(\mathit {dt}\) is the local subsequence-based decision tree learned for each X processed by lasts. Formally, the fidelity is the percentage of times in which \({\hat{{\bf y}}}= {\hat{{\bf y}}}^{\prime }\), i.e., \({\it fid} = \frac{1}{n} \sum _{i=1}^{n} \mathbb {1}_{{\hat{{\bf y}}}= {\hat{{\bf y}}}^{\prime }}\).

In [37], we showed that lasts is more faithful than a version of lasts using the same random neighborhood generation adopted by lime [78]. In the following, we compare the version of lasts proposed in [37] using the genetic-based neighborhood generation (named gen-shp) with the alternatives of lasts proposed in this paper by combining the various sampling strategies for the neighborhood generation and types of subsequences. Besides, we show that extracting an explanation from the local neighborhood of a given instance is a winning strategy compared to an approach that builds a single global interpretable surrogate. Thus, we compare lasts against shapelet/SAX-based global decision tree (\(\mathit {gdt}\)) classifiers, namely glo-shp and glo-sax, trained on \(\mathcal {X}_\text{bb}\). In this case, the fidelity is calculated as the accuracy between \({\hat{{\bf y}}}=b(\mathcal {X}_\text{exp})\) and \({\hat{{\bf y}}}^{\prime }=\mathit {gdt}(\mathcal {X}_\text{exp})\).

To help readability, we provide the average ranking (\(\overline{{\it rk}}\)) (lower is better) at the bottom of each table, which is the average of the ranks of each method for all datasets and black-boxes. Table 2 reports the values of the fidelity. We observe very high fidelity for all the approaches. However, the null hypothesis that all methods are equivalent is rejected (\(\mathit {p{-}value} \lt 0.01\)) from the non-parametric Friedman over multiple datasets and black-boxes. We notice that the global approaches glo-shp and glo-sax always have slightly lower values than the local alternatives of lasts. Among them, the variants employing SAX-based \(\mathit {dt}\) generally have higher fidelity than those using shapelet-based \(\mathit {dt}\). Concerning the neighborhood generation options, results are very similar even though mu seem to return slightly better results. We better elaborate on this fact in the following paragraph.

The following test is to estimate the precision and coverage of the factual and counterfactual rules (indicated as \(\mathit {pre}_{r_=}\), \(\mathit {pre}_{{r_\ne }}\), \(\mathit {cov}_{r_=}\), \(\mathit {cov}_{{r_\ne }}\), respectively). The coverage is measured as the relative number of time series that respect the factual/counterfactual rule premises in the neighborhood \({\hat{\mathcal {X}}}\), while the precision is measured as the relative number of covered time series for which the prediction outcome is correct. In addition, we measure also the cohesion and separation of the neighborhoods \({\hat{\mathcal {X}}}\) through the silhouette (\(\mathit {sil}\)) coefficient [85] with respect to the two clusters identified as \({\hat{\mathcal {X}}_=}\) and \({\hat{\mathcal {X}}_\ne }\) [32]. Table 3 illustrates the average precision, coverage, and silhouette scores when applying lasts for explaining the res classifier. For every metric, the higher the value, the better the score. Again, we notice that SAX-based methods seem to perform better than shapelet-based ones. More in detail, mu-sax has always the highest \(\mathit {pre}_{r_=}\), while ga-sax is among the bests for \(\mathit {pre}_{{r_\ne }}\). Factual coverage scores are, in general, quite close, while shapelet-based models usually have better counterfactual coverage. Finally, for the silhouette, the approaches using the uniform sphere neighgen function (mu and us) perform better by quite a margin. In summary, we can state that lasts explanations obtained using SAX and uniform sphere or matched uniform are the best ones. However, in the following, we only consider lasts using SAX-based subsequences and the mu neighborhood generation, unless otherwise specified, because the gm sampling does not guarantee a perfectly centered neighborhood generation around the closest counterexemplar. This problem can be detected by computing the Local Outlier Factor (LOF) [10] for \({\bf z}\), as implemented in sklearn. In particular, a LOF score of \(-1\) indicates an outlier, while a score of 1 indicates an inlier. In our tests LOF_MU = 1 for every dataset and for both black-boxes, while LOF_GM \(\in [-0.55, 0.64]\), indicating a higher degree of outlierness. Moreover, given that the radius of the uniform sampling is, by construction, as small as possible, the problem of distribution mismatch in a Gaussian space is minimized in most cases.

The time required by lasts for the explanation is mainly affected by the cfs algorithm, and the subsequence transform (\({\varsigma ^*}\)). The cfs algorithm’s most expensive operations are the black-box (b) and decoder (h) functions. They dominate the time complexity of all the sample functions described in Section 4.3, which is only \(O(q)\), where q is the dimensionality of the latent vector \({\bf z}\) [88]. Following Algorithm 2, b and h must be repeated until no counterexemplar in the latent space is found. Therefore, the worst-case scenario can appear when the latent encoding of \({\bf z}\) is extremely close to the decision boundary, and we have to reduce the \(\theta\) parameter a great number of times. This is extremely unlikely, and we found that in most of our tests, starting with \(\theta =2\), only a few iterations are needed to exit the loop, with the worst case not exceeding 20 iterations. In Table 4, we report the average runtime and standard deviation for the two critical phases and for the complete explanation process (all). We notice that on every dataset, the cfs strategy is one order of magnitude faster than the genetic one⁴. The subsequence extraction runtime is comparable only for the Coffee dataset but, on average, is at least four times faster when using SAX instead of shapelets. In conclusion, we notice that the complete explanation process is much faster in the novel version of lasts when using SAX with respect to the one proposed in [37] and the current one when using shapelets.

In summary, mu-sax is the most well-rounded approach: (i) it generates a neighborhood of well-separated exemplar and counterexemplar instances around the decision boundary, which contains \({\bf z}\) as an inlier, (ii) it has excellent performance as a local surrogate, with high fidelity scores, factual and counterfactual rules precision, and (iii) it ensures that the closest counterexemplar is sampled from a “probable” portion of latent space, minimizing the distribution mismatch.

Since it is hard to validate the usefulness of the generated factual and counterfactual instances with an experiment involving humans, inspired by [50], we tested their effectiveness with a memory-based machine learning technique. This experiment gives an objective and indirect estimation of the usefulness of exemplars and counterexemplars by checking how the performance of a simple classifier changes by training it on an increasing number of instances. For each instance \(X \in \mathcal {X}_\text{exp}\) we apply lasts, i.e., we encode X, find its closest counterexemplar, generate its latent neighborhood, and decode it into \({\hat{\mathcal {X}}}\). From this synthetic dataset \({\hat{\mathcal {X}}}\), using the black-box to retrieve the predicted labels, we extract n exemplars, i.e., n instances having the same class as X, and n counterexemplars from each other class, i.e., \(n(c-1)\) instances having a different class w.r.t. X. This extraction is random and without replacement. Then, we use the selected exemplars and counterexemplars to train a 1-Nearest Neighbor (1-NN), and we classify X. For each dataset, we compare the classification performances on X with a 1-NN trained on real time series, n per class, randomly selected from \(\mathcal {X}_\text{exp} \setminus \lbrace X\rbrace\). Specifically, \(n \in \lbrace 1, 2, 4, 8, 16\rbrace\) and the Euclidean distance is used as a distance function. Figure 7 shows the accuracy of a 1-NN using an increasing number of exemplars and counterexemplars as training. The plot aggregates the accuracy for each dataset at each n, showing the estimate of the central tendency and the respective confidence interval. On average, the performance of lasts is higher and more constant, revealing that few exemplars and counterexemplars are a good proxy for recognizing the classification outcome. On the other hand, the accuracy of real increases more steeply on the various datasets with increasing n, meaning that more exemplars and counterexemplars are necessary to distinguish the reasons for the classification. Hence, we can state that exemplars and counterexemplars help in discovering the decision boundary and in highlighting similarities and differences. This experiment shows that time series must be carefully chosen, as lasts does, in order to be suitable for class recognition.

In these experiments, we compare lasts with shap [64] employed to solve the time series black-box outcome explanation problem as in [4]. Similar to the proposal in [71], we adapt shap to time series classifiers by performing an adaptive segmentation [49] of the time series in order to divide it into meaningful intervals that are used as features by shap. As is commonly done with images, the absence of a feature is simulated via linear interpolation by connecting the observations before and after the ablated segment.

A common strategy to validate a saliency-based explanation is to observe how the performance of the black-box changes by adding/removing features in order of importance [76]. Regarding deletion, the intuition is that removing the most important time-steps in a time series will force the black-box to change its decision. On the other hand, the insertion evaluation adopts a complementary approach by starting from an “empty” time series and adding the most important time-steps. From a practical standpoint, the removal of time-steps is approximated by taking the average of the time series values. Time-steps are added/deleted in order of importance one by one, and the black-box prediction is checked at each step. Insertion/deletion metrics are computed as \(\text{AUC} / md\), i.e., the Area Under the Curve of the line obtained by computing the accuracy of the black-box after each insertion/deletion, divided by the total number of time-steps m of every dimension d of the time series. For the insertion benchmark, we want the black-box performance to improve as fast as possible; therefore, a high score is desirable. On the contrary, we expect a sudden drop in performance for deletion benchmarks, resulting in a lower score. Results can be seen in Table 5. For this benchmark, lasts scores very similarly to shap, indicating that both methods are able to indicate important observations of the time series.

Furthermore, we measure the stability of an explainer as its ability to produce a similar explanation for close instances, i.e., given similar instances, their saliency map should also be similar. In practice, given an instance X, we find its closest instance \(X^{\prime }\) in the latent space using the Euclidean distance. Then, we compare their saliency maps using the Mean Absolute Error (MAE). Intuitively, stable explanations should have lower MAE, while unstable explanations should have a higher error. As can be seen in Table 5, lasts outperforms shap in stability for all but one dataset, indicating its ability to give similar explanations for similar instances.

Since it is impossible to know exactly if explanations are correct with real datasets and standard black-box models, we use a synthetic experiment to check if the saliency maps obtained with agnostic explainers match with custom-defined ground truths, i.e., to check if they are correct. To perform this study, inspired by [31], we generate a univariate synthetic dataset having two classes. Time series belonging to each class have a distinct and easily identifiable pattern that unequivocally defines their class, plus some random noise. Then, we build a synthetic classifier that performs the classification by only looking at the presence of these predefined patterns. The time-steps of the time series that are checked by the synthetic model during classification are thus known, and the ground truth is defined as a vector of length m with the value 1 if the pattern is present at a given time-step and 0 otherwise. In other words, the only important points for the time series classification have a value of 1, while the noise has a 0 value. Given a saliency-based explanation returned as a vector of length m, we compare it with the ground truth by first normalizing it in the range \([0, 1]\) and then computing MAE between the ground truth and the normalized saliency vector. Intuitively, if the saliency vector correctly identifies important and irrelevant points, the MAE will tend to 0. To improve the test’s significance, we perform it on six synthetic classifiers that use different random and continuous subsets of the original patterns. A comparison of the saliency maps returned by lasts and shap w.r.t. the ground truth is depicted in Figure 8, while box-plots of the MAE for each dataset can be viewed in Figure 9. In general, lasts performs better, i.e., has a lower median MAE in four of the six tests. Moreover, lasts has a lower standard deviation, indicating that even its worst saliency maps are not that far from the ground truth. In general, we observe that lasts tends to give a more targeted and precise explanation, while shap tends to give importance to more points in the time series, resulting in a very wide interquartile range for the MAE.

To the best of our knowledge, lasts is the only local agnostic time series explainer that outputs rules as an explanation. Thus, we decided to compare the factual and counterfactual rules returned by lasts with those of a global decision tree, as in Section 5.3, and also with anchor [79]. In order to adapt anchor to time series, we consider each observation as a separate feature. As a note, anchor can return only factual rules whose conditions depend on single time series observations. We use precision and coverage metrics to evaluate the goodness of a rule. Moreover, given that simpler explanations are to be preferred, we also measure the length of the returned rules.

The results of precision, coverage and length for the factual rules returned by lasts, glo-sax and anchor are presented in Table 6. Regarding precision, anchor is the clear winner in all but one dataset. This result is not surprising given that anchor, by definition, constructs rules that are guaranteed to have a precision above 0.95. lasts scores slightly lower, with a precision that does not drop under 0.9, while glo-sax is the clear loser in this benchmark. The coverage metric helps to show the whole picture of these benchmarks. lasts performs best, followed by glo-sax and anchor in last place. This indicated that while rules returned by anchor are indeed slightly more precise, they are also less generalizable, i.e., they cover a much lower number of instances. Furthermore, as shown by the average lengths, factual rules returned by anchor and glo-sax are also considerably longer, i.e., more difficult to understand from a human standpoint. As a note, anchor is also extremely inefficient for longer time series, requiring hours of runtime to explain a single instance. For this reason, completing the benchmarks for the ART dataset was impossible. Regarding counterfactual rules, results are presented in Table 6. lasts performs better than a global surrogate in precision and length, tying in coverage. This experiment demonstrates that, while the number of instances covered by the rules is comparable for the two methods, counterfactual rules returned by lasts are more precise, shorter, and thus easier to understand.

This section shows qualitative examples from two real-world datasets, ECG5000 and Libras. Figure 10 presents an explanation of an instance from the ECG5000 dataset. ECG5000 contains \(5,\!000\) heartbeats belonging to five different classes, one corresponding Normal instances, and four corresponding to different kinds of Abnormal heartbeats. The instance we are trying to explain is correctly classified by rocket as Normal (Figure 10 top-left). By looking at the difference between exemplars and counterexemplars, we can clearly see that the main difference between normal and abnormal instances is in the rightmost part of the time series, which is much lower for abnormal time series, and presents an evident V-shape. More specifically, these abnormal series belong all to the Premature Ventricular Contraction class. The saliency map confirms the assessment deduced by the example-based explanation, highlighting only the last observations of the time series. The rules show the other main difference between classes. The factual rule, \({r_=}= \lbrace s_{221} \in X \wedge s_{264} \in X\rbrace \rightarrow {\it Normal}\), shows that normal instances contain subsequence \(s_{264}\), while the counterfactual rule, \({{r_\ne }} = \lbrace s_{221} \in X \wedge s_{264} \not\in X\rbrace \rightarrow {\it Premature Ventricular Contraction}\), shows that abnormal time series have a flatter shape, not containing \(s_{264}\). In general, we do not expect the saliency map and the rules to cover the same exact areas because the saliency map only highlights the parts of the time series that change the most between classes, whereas the rules can be based on subsequences emphasizing even a small shape change from any part of the time series.

In Figure 11, we present an explanation of a multivariate time series from the Libras dataset. This dataset contains instances having two signals each, belonging to 15 classes that correspond to different hand movements. The instance we are trying to explain is labeled as face-up curve, and it is correctly classified by rocket. The example-based part of the explanation shows exemplars and counterexemplars that are extremely similar to the naked eye. This means that even very small changes in the time series shape can result in a change of prediction from the black-box. In this sense, the closest instances belong all to class horizontal wavy. The most salient observations are highlighted in the saliency map and show that the top signal, and in particular the central sigmoidal part, is the most relevant for the change in classification. This is also confirmed by the factual and counterfactual rules, \({r_=}= \lbrace s_{412} \in {\bf x}_0 \wedge s_{1026} \not\in {\bf x}_0 \wedge s_{1855} \in {\bf x}_1 \wedge s_{2184} \in {\bf x}_1 \rbrace \rightarrow {\it face-up curve}\), \({{r_\ne }} = \lbrace s_{412} \in {\bf x}_0 \wedge s_{1026} \in {\bf x}_0 \wedge s_{1855} \in {\bf x}_1 \wedge s_{2184} \in {\bf x}_1 \rbrace \rightarrow {\it horizontal wavy}\). The rules show that the most significant subsequence is \(s_{1026}\), given that its presence/absence results in a change in class. Figure 11 (bottom-right) shows a counterexemplar which contains \(s_{1026}\), i.e., has a lower dip in the central sinusoidal pattern w.r.t. X.