Explaining Time Series Predictions with Dynamic Masks

How can we explain the predictions of a machine learning model? When the data is structured as a multivariate time series, this question induces additional difficulties such as the necessity for the explanation to embody the time dependency and the large number of inputs. To address these challenges, we propose dynamic masks (Dynamask). This method produces instance-wise importance scores for each feature at each time step by fitting a perturbation mask to the input sequence. In order to incorporate the time dependency of the data, Dynamask studies the effects of dynamic perturbation operators. In order to tackle the large number of inputs, we propose a scheme to make the feature selection parsimonious (to select no more feature than necessary) and legible (a notion that we detail by making a parallel with information theory). With synthetic and real-world data, we demonstrate that the dynamic underpinning of Dynamask, together with its parsimony, offer a neat improvement in the identification of feature importance over time. The modularity of Dynamask makes it ideal as a plug-in to increase the transparency of a wide range of machine learning models in areas such as medicine and finance, where time series are abundant.


Introduction and context
What do we need to trust a machine learning model? If accuracy is necessary, it might not always be sufficient. With the application of machine learning to critical areas such as medicine, finance and the criminal justice system, the black-box nature of modern machine learning models has appeared as a major hindrance to their large scale deployment (Caruana et al., 2015;Lipton, 2016;Ching et al., 2018). With the necessity to address this problem, the field Proceedings of the 38 th International Conference on Machine Learning, PMLR 139, 2021. Copyright 2021 by the author(s). of explainable artificial intelligence (XAI) thrived (Barredo Arrieta et al., 2020;Das & Rad, 2020;Tjoa & Guan, 2020).
Saliency Methods Among the many possibilities to increase the transparency of a machine learning model, we focus here on saliency methods. The purpose of these methods is to highlight the features in an input that are relevant for a model to issue a prediction. We can distinguish them according to the way they interact with a model to produce importance scores.
Gradient-based: These methods use the gradient of the model's prediction with respect to the features to produce importance scores. The premise is the following: if a feature is salient, we expect it to have a big impact on the model's prediction when varied locally. This translates into a big gradient of the model's output with respect to this feature. Popular gradient methods include Integrated Gradient (Sundararajan et al., 2017), DeepLIFT (Shrikumar et al., 2017) and GradSHAP (Lundberg et al., 2018).
Perturbation-based: These methods use the effect of a perturbation in the input on the model's prediction to produce importance scores. The premise is similar to the one used in gradient based method. The key difference lies in the way in which the features are varied. If gradient based methods use local variation of the features (according to the gradient), perturbation-based method use the data itself to produce a variation. A first example is Feature Occlusion (Suresh et al., 2017) that replaces a group of features with a baseline. Another example is Feature Permutation that performs individual permutation of features within a batch.
Attention-based: For some models, the architecture allows to perform simple explainability tasks, such as determining feature saliency. A popular example, building on the success attention mechanisms (Vaswani et al., 2017), is the usage of attention layers to produce importance scores (Choi et al., 2016;Song et al., 2017;Xu et al., 2018;Kwon et al., 2018).
Other: There are some methods that don't clearly fall in one of the above categories. A first example is SHAP (Lundberg & Lee, 2017), which attributes importance scores based on Shapley values. Another popular example is LIME in which the importance scores correspond to weights in a local linear model (Ribeiro et al., 2016a). Finally, some methods such as INVASE (Yoon et al., 2019a) or ASAC (Yoon et al., 2019b) train a selector network to highlight important features.
Time Series Saliency Saliency methods were originally introduced in the context of image classification (Simonyan et al., 2013). Since then, most methods have focused on images and tabular data. Very little attention has been given to time series (Barredo Arrieta et al., 2020). A possible explanation for this is the increasing interest in model agnostic methods (Ribeiro et al., 2016b).
Model agnostic methods are designed to be used with a very wide range of models and data structures. In particular, nothing prevents us from using these methods for Recurrent Neural Networks (RNN) trained to handle multivariate time series. For instance, one could compute the Shapley values induced by this RNN for each input x t,i describing a feature i at time t. In this configuration, all of these inputs are considered as individual features and the time ordering is forgotten. This approach creates a conceptual problem illustrated in Figure 1. Context matters. This graph shows a fictional security price over time. A simplistic model recommends to sell just after a local maximum, to buy just after a local minimum and to wait otherwise. When recommending to sell or to buy, 3 consecutive time steps are required as a context.
In this example, the time dependency is crucial for the model. At each time step, the model provides a decision based on the two previous time steps. This oversimplified example illustrates that the time dependency induces a context, which might be crucial to understand the prediction of some models. In fact, recurrent models are inherently endowed with this notion of context as they rely on memory cells. This might explain the poor performances of model agnostic methods to identify salient features for RNNs (Ismail et al., 2020). These methods are static as the time dependency, and hence the context, is forgotten when treating all the time steps as separate features.
If we are ready to relax the model agnosticism, we can take a look at attention based methods. In methods such as RE-TAIN (Choi et al., 2016), attention weights are interpreted as importance scores. In particular, a high attention weight for a given input indicates that this input is salient for the model to issue a prediction. In this case, the architecture of the model allows us to provide importance scores without altering the dynamic nature of the data. However, it has been shown recently that the attention weights of a model can significantly be changed without inducing any effect on the predictions (Jain & Wallace, 2019). This renders the parallel between attention weights and feature importance rather unclear. This discrepancy also appears in our experiments.
Another challenge induced by time series is the large number of inputs. The number of features is multiplied by the number of time steps that the model uses to issue a prediction. A saliency map can therefore become quickly overwhelming in this setting. To address these challenges, it is necessary to incorporate some treatment of parsimony and legibility in a time series saliency method. By parsimony, we mean that a saliency method should not select more inputs than necessary to explain a given prediction. For instance, a method that identifies 90% of the inputs as salient is not making the interpretation task easier. By legibility, we mean that the analysis of the importance scores by the user should be as simple as possible. Clearly, for long time sequences, analysing feature maps such as (c) and (d) on Figure 3 can quickly become daunting.
To address all of these challenges, we introduce Dynamask. It is a perturbation-based saliency method building on the concept of masks to produce post-hoc explanations for any time series model. Masks have been introduced in the context of image classification (Fong & Vedaldi, 2017;Fong et al., 2019). In this framework, a mask is fitted to each image. The mask highlights the regions of the image that are salient for the black-box classifier to issue its prediction. These masks are obtained by perturbing the pixels of the original image according to the surrounding pixels and study the impact of such perturbations on the black-box prediction. It has been suggested to extend the usage of masks beyond image classification (Phillips et al., 2018;Ho et al., 2019). However, to our knowledge, no available implementation and quantitative comparison with benchmarks exists in the context of multivariate time series. Moreover, few works in the literature explicitly mention explainability in a multivariate time series setting (Siddiqui et al., 2019;Tonekaboni et al., 2020). Both of these considerations motivate our proposition.
Contributions By building on the challenges that have been described, our work is, to our knowledge, the first saliency method to rigorously address the following questions in a time series setting.
(1) How to incorporate the context? In our framework, this is naturally achieved by studying the effect of dynamic perturbations. Concretely, a perturbation is built for each feature at each time by using the value of this feature at adjacent times. This allows us to build meaningful perturbations that carry contexts such as the one illustrated in Figure 1.
(2) How to be parsimonious? A great advantage of using masks is that the notion of parsimony naturally translates into the extremal property. Conceptually, an extremal mask selects the minimal number of inputs allowing to reconstruct the black-box prediction with a given precision.
(3) How to be legible? To make our masks as simple as possible, we encourage them to be almost binary. Concretely, this means that we enforce a polarization between low and high saliency scores. Moreover, to make the notion of legibility quantitative, we propose a parallel with information theory 1 . This allows us to introduce two metrics : the mask information and the mask entropy. As illustrated in Figure 3, the entropy can be used to assess the legibility of a given mask. Moreover, these metrics can also be computed for other saliency methods, hence allowing insightful comparisons.
The paper is structured as follows. In Section 2, we outline the mathematical formalism related to our method. Then, we evaluate our method by comparing it with several benchmarks in Section 3 and conclude in Section 4.

Mathematical formulation
In this section, we formalize the problem of feature importance over time as well as the proposed solution. For the following, it is helpful to have the big picture in mind. Hence, we present the blueprint of Dynamask in Figure 2. The rest of this section makes this construction rigorous.

Preliminaries
Let X ⊂ R d X be an input (or feature) space and Y ⊂ R d Y be an output (or label) space, where d X and d Y are respectively the dimension of the input and the output space. For the sake of concision, we denote by [n 1 : n 2 ] the set of natural numbers between the natural numbers n 1 and n 2 with n 1 < n 2 . We assume that the data is given in terms of time series (x t ) t∈[1:T ] , where the inputs x t ∈ X are indexed by a time parameter t ∈ [1 : T ] with T ∈ N * . We consider the problem of predicting a sequence (y t ) t∈[ty:T ] of outputs y t ∈ Y indexed by the same time parameter, but starting at a time t y ≥ 1, hence allowing to cover sequence to vector predictions (t y = T ) as well as sequence to sequence predictions (t y < T ). For the following, it is convenient to introduce a matrix notation for these time series. In this way, X = (x t,i ) (t,i)∈[1:T ]×[1:d X ] denotes the matrix in R T ×d X whose rows correspond to time steps and whose columns correspond to features. Similarly 2 , we denote Y = (y t,i ) (t,i)∈[ty:T ]×[1:d Y ] for the matrix in R (T +1−ty)×d Y . 1 Many connections exist between explainability and information theory, see for instance (Chen et al., 2018). 2 In the following, we do not make a distinction between R T ×d X and X T as these are isomorphic vector spaces. The same goes for R (T +1−ty )×d Y and Y T +1−ty .
Our task is to explain the prediction Y = f (X) of a black box f that has been pre-trained for the aforementioned prediction task. In other words, our method aims to explain individual predictions of a given black box. More specifically, our purpose is to identify the parts of the input X that are the most relevant for f to produce the prediction Y. To do this, we use the concept of masks introduced in (Fong & Vedaldi, 2017;Fong et al., 2019), but in a different context. Here, we adapt the notion of mask to a dynamic setting. Definition 2.1 (Mask). A mask associated to an input sequence X ∈ R T ×d X and a black box f : X T → Y T +1−ty is a matrix M = (m t,i ) ∈ [0, 1] T ×d X of the same dimension as the input sequence. The element m t,i of this matrix represents the importance of feature i at time t for f to produce the prediction Y = f (X). Remark 2.1. For a given coefficient in the mask matrix, a value close to 1 indicates that the feature is salient, while a value close to 0 indicates the opposite.
Now how can we obtain such a mask given a black box and an input sequence? To answer this question, we should let the mask act on the inputs and measure the effect it has on a black box prediction. As the name suggests, we would like this mask to hide irrelevant inputs contained in X. To make this rigorous, it is useful to spend some time to think about perturbation operators.

Perturbation operators
The method that we are proposing in this paper is perturbation-based. Concretely, this means that a mask is used to build a perturbation operator. Examples of such operators are given in (Fong et al., 2019) in the context of image classification. Here, we propose a general definition and we explain how to take advantage of the dynamic nature of the data to build meaningful perturbations. By recalling that the mask coefficients indicate the saliency, we expect the perturbation to vanish for features x t,i whose mask coefficient m t,i is close to one. This motivates the following definition. Definition 2.2 (Perturbation operator). A perturbation operator associated to a mask M ∈ [0, 1] T ×d X is a linear operator acting on the input sequence space Π M : R T ×d X → R T ×d X . It needs to fulfil the two following assumptions for any given (t, i) ∈ [1 : T ] × [1 : d X ]: 1. The perturbation for x t,i is dictated by m t,i : where π is differentiable for m ∈ (0, 1) and continuous at m = 0, 1.
2. The action of the perturbation operator is trivial when the mask coefficient is set to one :

Loss
Black Box Figure 2. Diagram for Dynamask. An input matrix X, extracted from a multivariate time series, is fed to a black-box to produce a prediction Y. The objective is to give a saliency score for each component of X. In Dynamask, these saliency scores are stored in a mask M of the same shape as the input X. To detect the salient information in the input X, the mask produces a perturbed version of X via a perturbation operator Π. This perturbed X is fed to the black-box to produce a perturbed prediction YM. The perturbed prediction is compared to the original prediction and the error is backpropagated to adapt the saliency scores contained in the mask.
Moreover, we say that the perturbation is dynamic if for any given feature at any given time, the perturbation is constructed with the values that this feature takes in neighbouring times. More formally, the perturbation is dynamic if there exist a couple (W 1 , Remark 2.2. The first assumption ensures that the perturbations are applied independently for all inputs and that our method is suitable for gradient-based optimization. The second assumption 3 translates the fact that the perturbation should have less effect on salient inputs. Remark 2.3. The parameters W 1 and W 2 control how the perturbation depends on neighbouring times. A non-zero value for W 1 indicates that the perturbation depends on past time steps. A non-zero value for W 2 indicates that the perturbation depends on future time steps. For the perturbation to be dynamic, at least one of these parameters has to be different from zero.
This definition gives the freedom to design perturbation operators adapted to particular contexts. The method that we develop can be used for any such perturbation operator. In our case, we would like to build our masks by taking advantage the dynamic nature of the data into account. This is crucial if we want to capture local time variations of the features, such as a quick increase of the blood-pressure or 3 Together with the continuity of π with respect to the mask coefficients.
an extremal security price from Figure 1. This is achieved by using dynamic perturbation operators. To illustrate, we provide three examples: where π g is a temporal Gaussian blur 4 with Similarly, π m can be interpreted as a fade to moving average perturbation with where W ∈ N is the size of the moving window. Finally, π p is similar to π m with one difference: the former only uses past values of the features to compute the perturbation 5 : Note that the Hadamard product M X used in (Ho et al., 2019) is a particular case of π m with µ t,i = 0 = so that this perturbation is static. In the following, to stress that a mask is obtained by inspecting the effect of a dynamic perturbation operator, we shall refer to it as a dynamic mask (or Dynamask in short).
In terms of complexity, the computation of π g requires O(d X · T 2 ) operations 6 while π m requires O(d X · T · W ) operations 7 . When T is big 8 , it might therefore be more interesting to use π m with W T or a windowed version of π g . With this analysis of perturbation operators, everything is ready to explain how dynamic masks are obtained.

Mask optimization
To design an objective function for our mask, it is helpful to keep in mind what an ideal mask does. From the previous subsection, it is clear that we should compare the black-box predictions for both the unperturbed and the perturbed input. Ideally, the mask will identify a subset of salient features contained in the input that explains the black-box prediction. Since this subset of features is salient, the mask should indicate to the perturbation operator to preserve it. More concretely, a first part of our objective function should keep the shift in the black-box prediction to be small. We call it the error part of our objective function. In practice, the expression for the error part depends on the task done by the black-box. In a regression context, we minimize the squared error between the unperturbed and the perturbed prediction: Similarly, in a classification task, we minimize the crossentropy between the predictions: Now we have to make sure that the mask actually selects salient features and discards the others. By remembering that m t,i = 0 indicates that the feature x t,i is irrelevant for the black-box prediction, this selection translates into imposing sparsity in the mask matrix M. A first approach for this, used in (Fong & Vedaldi, 2017;Ho et al., 2019), is to add a l 1 regularisation term on the coefficients of M. However, it was noted in (Fong et al., 2019) that this produces mask that vary with the regularization coefficient λ in a way that renders comparisons between different λ difficult. To solve this issue, they introduce a new regularization term 6 One must compute a perturbation for each dX · T component of the input and the sums in each perturbation have T terms 7 One must compute a perturbation for each dX · T component of the input and each moving average involves 2W + 1 terms 8 In our experiments, however, we never deal with T > 100 so that both approaches are reasonable. to impose sparsity: where . denotes the vector 2-norm, vecsort is a function that vectorizes M and then sort the elements of the resulting vector in ascending order. The vector r a contains (1 − a) · d X · T zeros followed by a · d X · T ones, where a ∈ [0, 1]. In short, this regularization term encourages the mask to highlight a fraction a of the inputs. For this reason, a can also be referred to as the area of the mask. A first advantage of this approach is that the hyperparameter a can be modulated by the user to highlight a desired fraction of the input. For instance, one can start with a small value of a and slide it to higher values in order to see features gradually appearing in order of importance. Another advantage of this regularization term is that it encourages the mask to be binary, which makes the mask more legible as we will detail in the next section.
Finally, we might want to avoid quick variations in the saliency over time. This could either be a prior belief or a preference of the user with respect to the saliency map. If this is relevant, we can enforce the salient regions to be connected in time with the following loss that penalizes jumps of the saliency over time: For a fixed fraction a, the mask optimization problem can therefore be written as Note that, in this optimization problem, the fraction a is fixed. In some contexts, one might want to find the smallest fraction of input features that allows us to reproduce the black-box prediction with a given precision. Finding this minimal fraction corresponds to the following optimization problem: where ε is the threshold that sets the acceptable precision. The resulting mask M * a * is then called extremal mask. The idea of an extremal mask is extremely interesting by itself: it explains the black-box prediction in terms of a minimal number of salient features. This is precisely the parsimony that we were referring to in the introduction.
Finally, it is worth mentioning that we have here presented a scheme where the mask preserves the features that minimize the error. There exists a variant of this scheme where the mask preserves the features that maximize the error. This other scheme, together with the detailed optimization algorithm used in our implementation, can be found in Appendix B.

Masks and information theory
Once the mask is obtained, it highlights the features that contain crucial information for the black-box to issue a prediction. This motivates a parallel between masks and information theory. To make this rigorous, we notice that the mask admits a natural interpretation in terms of information content. As aforementioned, a value close to 1 for m t,i indicates that the input feature x t,i carries information that is important for the black box f to predict an outcome. It is therefore natural to interpret a mask coefficient as a probability that the associated feature is salient for the black box to issue its prediction. It is tempting to use this analogy to build the counterpart of Shannon information content in order to measure the quantity of information contained in subsequences extracted from the time series. However, this analogy requires a closer analysis. We recall that the Shannon information content of an outcome decreases when the probability of this outcome increases (Shannon, 1948;MacKay, 2003;Cover & Thomas, 2005). In our framework, we would like to adapt this notion so that the information content increases with the mask coefficients (if m t,i gets closer to one, this indicates that x t,i carries more useful information for the black-box to issue its prediction). To solve this discrepancy, we have to use 1−m t,i as the pseudoprobability appearing in our adapted notion of Shannon information content.
Definition 2.3 (Mask information). The mask information associated to a mask M and a subsequence ( Remark 2.4. Conceptually, the information content of a subsequence measures the quantity of useful information it contains for the black-box to issue a prediction. It allows us to associate a saliency score to a group of inputs according to the mask. Remark 2.5. Note that, in theory, the mask information diverges when a mask coefficient is set to one. In practice, this can be avoided by imposing M ∈ (0, 1) T ×d X .
As in traditional information theory, the information content is not entirely informative on its own. For instance, consider two subsequences indexed by, respectively, A, B with |A| = |B| = 10. We assume that the submask extracted from M with A contains 3 coefficients m = 0.9 and 7 coefficients m = 0 so that the information content of A is I M (A) ≈ 6.9. Now we consider that all the coefficient extracted from M with B are equal to 0.5 so that I M (B) ≈ 6.9 and hence I M (A) ≈ I M (B). In this example, A clearly identifies 3 important features while B gives a mixed score for the 10 features. Intuitively, it is pretty clear that the information provided by A is sharper. Unfortunately, the mask information by itself does not allow to distinguish these two subsequences. Hopefully, a natural distinction is given by the counterpart of Shannon entropy.
Definition 2.4 (Mask entropy). The mask entropy associated to a mask M and a subsequence ( Remark 2.6. We stress that the mask entropy is not the Shannon entropy for the subsequence (x t,i ) (t,i)∈A . Indeed, the pseudo-probabilities that we consider are not the probabilities p(x t,i ) for each feature x t,i to occur. In particular, there is no reason to expect that the probability of each input decouple from the others so that it would be wrong to sum the individual contribution of each input separately 9 like we are doing here.
Clearly, it is desirable for our mask to provide explanations with low entropy. This stems from the fact that the entropy is maximized when mask coefficients m t,i are close to 0.5. In this case, given our probabilistic interpretation, the mask coefficient is ambiguous as it does not really indicate whether the feature is salient. This is consistent with our previous Since masks coefficients take various values in practice, masks with higher entropy appear less legible, as illustrated in Figure 3. Therefore, we use the entropy as a measure of the mask's sharpness and legibility in our experiments. In particular, the entropy is minimized for perfectly binary masks M ∈ {0, 1} T ×d X , which are easy to read and contain no ambiguity. Consequently, the regularization term L a (M) that we used in Section 2.3 has the effect of reducing the entropy. Since our adapted notions of information and entropy rely on individual contributions from each feature, they come with natural properties.  Masks with high entropy appear less legible, especially for long sequences.
Proof. See Appendix A.
Together, these properties guarantee that I M and S M define measures for the discrete σ-algebra

Experiments
In this section 10 , we evaluate the quality of our dynamic masks. There are two big difficulties to keep in mind when evaluating the quality of a saliency method. The first one is that the true importance of features is usually unknown with real-world data. The second one is that the the performance of a saliency method in identifying relevant features depends on the black-box and its performances. In order to illustrate these challenges, we propose three different experiments in ascending order of difficulty. In the first experiment, a white box with known feature importance is used so that both difficulties are avoided. In the second experiment, a blackbox trained on a dataset with known feature importance is used and hence only the second difficulty is encountered. In the third experiment, a black-box trained on a real-world clinical dataset is used and hence both difficulties are encountered. For each experiment, more details are given in Appendix C.

Feature importance for a white-box regressor
Experiment In this experiment, we work with a trivial white-box regressor whose predictions only rely on a known 10 Our implementation can be found at https://github. com/JonathanCrabbe/Dynamask.
where A T and A X respectively give the salient times and features: Note that, in this experiment, d Y = 1 so that we can omit the second index for f (X). In our experiment, we consider two scenarios that are known to be challenging for state of the art saliency methods (Ismail et al., 2020). In the first one, the white-box depends on a small portion of salient features |A X | d X . In the second one, the white-box depends on a small portion of salient times |A T | T . We fit a mask to this black-box by using the squared error loss together with the regularization. In our experiment: T = d X = 50. The salient time and features are selected randomly, the input features are generated with an ARMA process. We repeat the experiment 10 times.
Metrics Since the true salient features are known unambiguously in this setup, we have a large set of metrics at hand to evaluate the performance of a saliency method. To measure the proportion of identified features that are indeed salient, we use the area under the precision curve (AUP, higher is better). To measure the portion of salient features that have indeed been identified, we use the area under the recall curve (AUR, higher is better). To measure how much information the saliency method predicts for the salient region, we use the mask information content (I M (A), higher is better). To measure the sharpness of the explanation in this region, we use the mask entropy (S M (A), lower is better).
Benchmarks The fact that we are dealing with a whitebox regressor already disqualifies some methods such as FIT (Tonekaboni et al., 2020) or DeepLIFT. We compare our method with Feature Occlusion (FO), Feature Permutation (FP), Integrated Gradient (IG) and Shapley Value Sampling (SVS) (Castro et al., 2009). For fair comparison, we use these baselines to evaluate the importance of each feature at each given time. As explained in Section 2.4, a mask can be associated to each saliency method.

Feature importance for a black-box classifier
Experiment We reproduce the state experiment from (Tonekaboni et al., 2020) but in a more challenging setting. In this case, the data is generated according to a 2-state hidden Markov model (HMM) whose state at time t ∈ [1 : T ] (T = 200) is denoted s t ∈ {0, 1}. At each time, the input feature vector has three components (d X = 3) and is generated according to the current state via x t ∼ N µ st , Σ st . To each of these input vectors is associated a binary label y t ∈ {0, 1}. This binary label is conditioned by one of the three component of the feature vector, based on the state: The label is then emitted via a Bernoulli distribution y t ∼ Ber(p t ). The experiment proposed in (Tonekaboni et al., 2020) focuses on identifying the salient feature at each time where a state transition occurs. However, it is clear from the above discussion that there is exactly one salient feature at any given time. More precisely, the set of salient indexes is A = { (t, 2 + s t ) | t ∈ [1 : T ] }. We generate 1000 such time series, 800 of them are used to train a RNN black-box classifier f with one hidden layer made of 200 hidden GRU cells 11 . We then fit an extremal mask to the black-box by minimizing the cross entropy error for each test time series. We repeat the experiment 5 times.
Metrics We use the same metrics as before.

Feature importance on clinical data
Experiment We reproduce the MIMIC mortality experiment from (Tonekaboni et al., 2020). We fit a RNN blackbox with 200 hidden GRU cells to predict the mortality of a patient based on 48 hours (T = 48) of patient features (d X = 31). This is a binary classification problem for which the RNN estimates the probability of each class. For each patient, we fit a mask with area a to identify the most important observations x t,i . Since the ground-truth feature saliency is unknown, we must find an alternative way to assess the quality of this selection. If these observations are indeed salient, we expect them to have a big impact on the black-box's prediction if they are replaced. For this reason, we replace the most important observations by the time average value for each associated feature x t,i →x t,i = 1 T T t=1 x t,i . We then compare the prediction for the original input f (X) with the prediction for the input where the most important observations have been replaced f X . A big shift in the prediction indicates that these observations are salient. We repeat this experiment 3 times for various values of a.
Dataset We use the MIMIC-III dataset , that contains the health record of 40, 000 ICU deidentified patients at the Beth Israel Deaconess Medical Center. The selected data and its preprocessing is the same as the one done by (Tonekaboni et al., 2020).
Metrics To estimate the importance of the shift induced by replacing the most important observations, we use the crossentropy between f (X) and f X (CE, higher is better). To evaluate the number of patient whose prediction has been flipped, we compute the accuracy of f X with respect to the initial predictions f (X) (ACC, lower is better).
Benchmarks We use the same benchmarks as before.

Discussion
The results are shown in Figure 4 & 5. The observations selected by Dynamask have the most significant impact when replaced. The high accuracy suggests that the perturbation we use (replacing the most important observations by their time average) rarely produces a counterfactual input.

Conclusion
In this paper, we introduced Dynamask, a saliency method specifically designed for multivariate time series. These masks are endowed with an insightful information theoretic interpretation and offer a neat improvement in terms of performance. Dynamask has immediate applications in medicine and finance, where black-box predictions require more transparency.
For future works, it would be interesting to design more sophisticated consistency tests for saliency methods in a dynamic setting, like the ones that exist in image classification (Adebayo et al., 2018). This could be used to study the advantages or disadvantages of our method in more details.
Another interesting avenue would be to investigate what the dynamic setting can offer to provide richer explanations with some treatment of causality (Moraffah et al., 2020).

A. More details on the mathematical formulation A.1. Proofs
In this subsection section, we prove the propositions form the main paper. Proof. Let us proof all the properties one by one.
Positivity By definition, all coefficients from the mask are normalized: Positivity follows trivially from the properties of the logarithm function The same goes for the entropy Additivity We first note that the proposition follows trivially if the sets are distinct A ∩ B = ∅: (2) We shall now prove the additivity property in general by using these two ingredients. First we note that the set A ∪ B can be written as the disjoint union A [B \ (A ∩ B)]. It follows that where we have used the additivity property for disjoint sets in the first equality and the fact that (A ∩ B) ⊂ B in the second equality. The same reasoning holds for the entropy.
Monotonicity To prove the monotonicity property, it is useful to note that if A ⊂ B, we can use (2) to write where we have used the information positivity to produce the inequality. The same reasoning holds for the entropy.

A.2. Normalized information and entropy
In this subsection, we introduce the normalized counterparts of our information theoretic metrics. It is important to keep in mind that all the available information for issuing a black-box prediction is in [1 : T ] × [1 : d X ]. Therefore, the monotonicity property allows to introduce normalized counterparts of the information and the entropy.
Definition A.1 (Normalized metrics). The normalized mask information associated to a mask M and a subsequence .
The same goes for the related normalized mask entropy .
Remark A.1. By the monotonicity and the positivity properties, it is clear that This gives a natural interpretation of these quantities as being, respectively, the fraction of the total information and entropy contained in the subsequence A according to the mask M.
The normalized version of the metrics allow to measure what percentage of the total mask information/entropy is contained in a given subsequence.

A.3. Definition of a mask for other saliency methods
In this section, we explain how to associate a mask to any saliency method. Suppose that a given method produces a score matrix R ∈ R T ×d X that assigns an importance score r t,i for each element x t,i of the input matrix X. Then, if we normalize the coefficients of the score matrix, we obtain an associated mask: where (1) T ×d X denotes a T × d X matrix with all elements set to 1. This mask can subsequently be used to compute the mask information content and entropy. In our experiments, we use this correspondence to compare our method with popular saliency methods.
B. More details on the implementation B.1. Algorithm

Algorithm 1 Dynamask
Input: input sequence X ∈ R T ×d X , black-box f , perturbation operator Π, mask area a ∈ [0, 1], learning rate η ∈ R + , momentum α ∈ R + , initial size regulator The mask optimization algorithm is presented in Algorithm 1. In the algorithm, we used the notation (0.5) T ×d X for a T × d X matrix with all elements set 13 to 0.5. Similarly, (0) T ·d X ·(1−a) denotes a vector with T · d X · (1 − a) components set to 0 and (1) T ·d X ·a denotes a vector with T · d X · a components set to 1. The symbol ⊕ denotes the direct sum between two vector spaces, which is equivalent to the concatenation in Algorithm 1. The error part of the loss L e depends on the task (regression or classification), as explained in Section 3 of the paper. The momentum and the learning rate are typically set to 1, the number of epoch is typically 1000. We also use the clamp function, which is defined component by component as Finally, we note that the mask size regularization coefficient λ a grows exponentially during the optimization to reach a maximum value of δ · λ 0 at the end of the optimization. In practice, it is initialized to a small value (typically λ 0 = 0.1) and dilated by several order of magnitude during the optimization (typically δ = 1000). In this way, the optimization procedure works in two times. At the beginning, the loss is dominated by the error L e (M) so that the mask increases the mask coefficients of salient features. As the regulation coefficient λ a increases, the regulation term becomes more and more important so that the mask coefficients are attracted to 0 and 1. At the end of the optimization, the mask is almost binary.

B.2. Deletion variant
We notice that Algorithm 1 produces a mask that highlights the features that allow to reproduce the black-box prediction by keeping the error part of the loss L e (M) to be small. However, it is possible to highlight important features in another way. For instance, we could try to find the features that maximizes the prediction shift when perturbed. In this alternative formulation, the mask is obtained by solving the following optimization problem: Note that, in this case, the sign of the error part is flipped in order to maximizes the shift in the prediction. Moreover, the error is now evaluated for 1 − M rather than M. This is because important features are maximally perturbed in this case. In this way, a salient feature x t,i can keep a mask coefficient m t,i close to 1 while being maximally perturbed. The regulator stays the same in this deletion variant, as the mask area still corresponds to the number of mask coefficients set to 1. We use the deletion variant to obtain the masks in the experiment with clinical data in the main paper.

C. More details on the experiments C.1. Metrics
We give the precise definition of each metric that appears in the experiments. Let us start with the metrics that are defined when the true importance is known. Definition C.1 (AUP,AUR). Let Q = (q t,i ) (t,i)∈[1:T ]×[1:d X ] be a matrix in {0, 1} T ×d X whose elements indicate the true saliency of the inputs contained in X ∈ R T ×d X . By definition, q t,i = 1 if the feature x t,i is salient and 0 otherwise. Let M = (m t,i ) (t,i)∈[1:T ]×[1:d X ] be a mask in [0, 1] T ×d X obtained with a saliency method. Let τ ∈ (0, 1) be the detection threshold for m t,i to indicate that the feature x t,i is salient. This allows to convert the mask into an estimatorQ(τ ) = (q t,i (τ )) (t,i)∈[1:T ]×[1: Consider the sets of truly salient indexes and the set of indexes selected by the saliency method We define the precision and recall curves that map each threshold to a precision and recall score: The AUP and AUR scores are the area under these curves Remark C.1. Roughly speaking, we consider the identification of salient features as a binary classification task. Each saliency method can thus be seen as a binary classifier for which we compute the AUP and the AUR. Remark C.2. Integrating over several detection thresholds allows to evaluate a saliency method with several levels of tolerance on what is considered as a salient feature.
In our experiment with MIMIC-III, since the ground true feature importance is unknown, we use the following metrics defined for a binary classification problem 14 . Definition C.2 (CE, ACC). Consider a classifier f that maps the input X to a probability f (X) ∈ [0, 1]. LetX be a perturbed input produced by a saliency method 15 . We define the function that converts a probability into a class class(p) = 0 if p < 0.5 1 else.
14 In our experiment, each input corresponds to a patient. Class 0 indicates that the patient survives and class 1 indicates that the patient dies. 15 In our experiment, we replace the most important features by the time average of the corresponding feature.
To measure the shift in the classifier's prediction caused by the perturbation of the input for several test examples {X k | k ∈ [1 : K]}, we use the binary cross-entropy (or logloss) To measure the number of prediction flipped by the perturbation, we use the accuracy We reproduce our experiment several times to get an average and a standard deviation for all of these metrics.

C.2. Computing infrastructure
All our experiments have been performed on a machine with Intel(R) Core(TM) i5-8600K CPU @ 3.60GHz [6 cores] and Nvidia GeForce RTX 2080 Ti GPU.

C.3. Details on the rare experiment
Data generation Since this experiment relies on a whitebox, we only have to generate the input sequences. As we explain in the main paper, each feature sequence is generated with an ARMA process: with ϕ 1 = 0.25, ϕ 2 = 0.1, ϕ 3 = 0.05 and t ∼ N (0, 1).
We generate one such sequence with t ∈ [1 : 50] for each feature i ∈ [1 : 50] by using the Python statsmodels library.
In the rare feature experiment, 5 features are selected as salient. Their indices are contained in A X and drawn uniformly without replacement from Mask fitting For each time series, we fit a mask by using the temporal Gaussian blur π g as a perturbation operator with σ max = 1 and by using the squared error loss. A mask is fitted for each value of a ∈ (n + 1) · 10 −3 | n ∈ [0 : 49] .
In our experiments, we don't consider a > 0.05. This is because we found experimentally that the error L e (M * a )  Runtime For rare time, finding the best mask takes on average 15.7s. For rare feature, finding the best mask takes on average 20.7s.
Illustrations To illustrate the results of our experiments, we show the saliency masks produced by various methods for the rare feature experiment in Figure 8 & 9 and for the rare time experiment in Figure 10 & 11. For all of these examples, we notice that Dynamask identifies a bigger portion of the truly salient inputs, which illustrates the bigger AUR reported in the main paper.

C.4. Details on the state experiment
Data generation The data generation is governed by a Hidden Markov Model (HMM). The initial distribution vector for this HMM is given by π = (0.5, 0.5) and its transition matrix is C = 0.1 0.9 0.1 0.9 .
At each time, the input feature vector has three components (d X = 3) and is generated according to the current state via x t ∼ N µ st , Σ st with mean vectors depending on To each of these input vectors is associated a binary label y t ∈ {0, 1}. This binary label is conditioned by one of the three component of the feature vector, based on the state: The length of each time series is fixed to 200 (T = 200).
We generate 1000 such time series, 800 are used for model training and 200 for testing.

Model training
We train a RNN with one layer made of 200 GRU cells trained using the Adam optimizer for 80 epochs (lr = 0.001, β 1 = 0.9, β 2 = 0.999 and no weight decay).
Mask fitting For each test time series, we fit a mask by using the temporal Gaussian blur π g as a perturbation operator with σ max = 1. A mask is optimized for each value of a ∈ 0.15 + 2n · 10 −2 | n ∈ [0 : 10] . We keep the extremal mask for a threshold set to ε = 0.9 · L e M = (1) T ×d X .
Runtime Finding the extremal mask for a given input takes 49.8s on average.
Illustrations To illustrate the results of our experiments, we show the saliency masks produced by various methods on Figure 12 & 13. By inspecting these figures, we notice that only Dynamask and Integrated Gradients seem to produce saliency maps where the imprint of the true saliency can be distinguished. One advantage of Dynamask is the contrast put between these salient inputs and the rest. In the case of Integrated Gradients, we see that many irrelevant inputs are assigned an important saliency score, although smaller than the truly salient inputs. This is because gradients are computed individually, without the goal of achieving parsimonious feature selection.
In addition, we have reported the mask entropy for each of the methods. As claimed in the main paper, Dynamask produces mask that have significantly lower entropy. Among the methods that produce masks with high entropy, we notice two trends. Some methods, such as RETAIN, produce masks where a significant portion of the inputs are assigned a mask entropy close the 0.5. As discussed in the main paper, these significance of the saliency scores is limited in this situation, since no clear contrast can be drawn between the saliency of different inputs. On the other hand, some methods like FIT produce masks with many different masks coefficients, which renders the saliency map somewhat fuzzy. In both cases, the high entropy detects these obstructions for legibility.

C.5. Details on the mimic experiment
Data preprocessing The data preprocessing used here is precisely the same as the one described in (Tonekaboni et al., 2020), we summarize it here for completeness. We use the adult ICU admission data from the MIMIC-III dataset   Mask fitting For each test patient, we simply fit a mask with a = 0.1 by maximizing the cross-entropy loss in the deletion variant formulation of Dynamask. We use the fadeto-moving average perturbation π m with W = 48. The hyperparameters for this optimization procedure are η = 1, α = 1, λ 0 = 0.1, λ c = 0, δ = 1000, N = 1000.
Runtime Fitting a mask for a given patient takes 3.58 s on average.
Illustrations To illustrate the results of our experiments, we show the 10% most important features for patients that are predicted to die on Figure 14 & 15 and for patients that are predicted to survive on Figure 16 & 17. In each case, we indicate the cross-entropy between the unperturbed and the perturbed prediction, as defined in Definition C.2. We note that Dynamask identifies the features that create the biggest shift. Qualitatively, we note that Dynamask seems to focus much more on the input that appear at latter time. This is consistent with the observations in (Ismail et al., 2019): these inputs are the most important for the black-box, since it is trained to predict the mortality after the 48 hours and RNNs have short term memory.

C.6. Influence of the perturbation operator
To study the effect of the perturbation operator choice, we have performed the following experiment: in the setup of the state experiment, we optimize 100 masks on distinct examples by using a Gaussian blur perturbation (π g , σ max = 1) and a fade-to-moving average perturbation (π m , W = 3). We do the same with a fade-to-past average perturbation that only uses past values of the features: (π p , W = 6). For each pair of perturbation operators, we compute the average accuracy between the associated masks (i.e. the fraction of inputs where both masks agree). For each method, we report the AUROC for the identification of true salient inputs. The results are reported in Table 4. We observe that π g , π m , π p generally agree and offer similar performances. The black-box predicts that this patient will survive. In each case, the cross entropy (CE) between the unperturbed and the perturbed prediction is reported.