Title: EXPLORING THE INTERPRETABILITY OF LSTM NEU-RAL NETWORKS OVER MULTI-VARIABLE DATA

Abstract: In learning a predictive model over multivariate time series consisting of target and exogenous variables, the forecasting performance and interpretability of the model are both essential for deployment and uncovering knowledge behind the data. To this end, we propose the interpretable multi-variable LSTM recurrent neural network (IMV-LSTM) capable of providing accurate forecasting as well as both temporal and variable level importance interpretation. In particular, IMV-LSTM is equipped with hidden state matrix and update process, so as to learn variableswise hidden states. On top of it, we develop a mixture attention mechanism and associated summarization methods to quantify the temporal and variable importance in data. Extensive experiments using real datasets demonstrate the prediction performance and interpretability of IMV-LSTM in comparison to a variety of baselines. It also exhibits the prospect as an end-to-end framework for both forecasting and knowledge extraction over multi-variate data.
Under review as a conference paper at ICLR 2019ICLR 2016;; Vinyals et al., 2015; Chorowski et al., 2015; Bahdanau et al., 2014) . Nevertheless, current attention mechanism is mostly applied to hidden states across time steps and capture globally temporal information, thereby failing to uncover fine-grained variable level importance.
To this end, in this paper we build interpretable LSTM models with the aim to achieve a unified framework of forecasting and knowledge extraction. In particular, the contribution is fourfold. First, we propose the interpretable multi-variable LSTM, referred to as IMV-LSTM, with hidden state matrix and updating scheme, such that each element of the hidden state matrix encodes information for a certain input variable. Second, based on these variable-wise hidden states we develop a novel mixture temporal and variable attention mechanism. Third, attention values are further summarized to quantify variable-wise temporal importance and overall variable importance. Lastly, we perform extensive experimental evaluation of IMV-LSTM against statistical, machine learning and neural network baselines to demonstrate the superior prediction performance and interpretability of IMV-LSTM. The idea of IMV-LSTM easily applies to other RNN structures, e.g. GRU and stacked recurrent layers. This will be the future work.
In time series analysis, prediction with exogenous variables can formulated as an auto-regressive exogenous model or prediction modeles defined on exogenous variables. Vanilla RNNs have been used to study it in (Zemouri et al., 2010; Diaconescu, 2008) , where interpretability was not investigated.
Recent research on the interpretability of RNNs is categorized into two groups: attention methods and post-analyzing over trained models. Attention mechanism has gained tremendous popularity (

1. INTRODUCTION:

Our daily life is now surrounded by various types of sensors, ranging from smart phones, video cameras, Internet of things, to robots. The observations yield by such devices over time are naturally organized in time series data (Qin et al., 2017; Yang et al., 2015) . In this paper, we focus on multivariable time series consisting of target and exogenous variables. Each variable corresponds to a monitoring over physical world. A predictive model over such multi-variable data aims to predict the future values of the target series using historical values of target and exogenous series.

In addition to forecasting, the interpretability of prediction models is essential for knowledge extraction, variable selection and so on (Hu et al., 2018; Foerster et al., 2017; Lipton, 2016) . For multi-variable time series in this paper, we focus on two types of importance interpretation. (1) Variable-wise temporal importance: exogenous variables present different temporal influence on the target one (Kirchgässner et al., 2012) . For instance, for the exogenous variable instantaneously correlated to the target one, its historical data at short time lags is expected to high importance values. (2) Overall variable importance: exogenous variables and the auto-regressive part of the target variable differ in predictive power, which reflects different variable importance w.r.t. the prediction of the target (Feng et al., 2018; Riemer et al., 2016) . The ability to unveil such knowledge through predictive models enables to fundamentally understand the relevance of exogenous variables to the target one.

Recently, recurrent neural networks (RNNs), especially long short-term memory (LSTM) (Hochreiter & Schmidhuber, 1997) and the gated recurrent unit (GRU) (Cho et al., 2014) , have been proven to be powerful sequence modeling tools in a variety of tasks such as language modelling, machine translation, health informatics, time series, and speech (Ke et al., 2018; Lin et al., 2017; Lipton et al., 2015; Sutskever et al., 2014; Bahdanau et al., 2014) .

However, current RNNs fall short of the aforementioned interpretability for multi-variable data due to their opaque internal states. Specifically, when fed with the multi-variable observations of the target and exogenous variables, RNNs blindly blend the information of all variables into memory cells and hidden states which are used for prediction. It is intractable to distinguish the contribution of individual variables into the prediction through hidden states (Zhang et al., 2017) . Recently, attention-based neural networks have been proposed to enhance the ability of RNN in selectively using long-term memory and the interpretability (Vaswani et al., 2017; Qin et al., 2017; Choi et al., As for post-analyzing based interpretation, (Murdoch et al., 2018; Murdoch & Szlam, 2017; Arras et al., 2017) extracted temporal importance scores over words or phrases of individual language sequences by decomposing the memory cells of trained LSTM. (Chu et al., 2018) proposed interpretation solutions for piece-wise linear neural networks. In (Wang et al., 2018) , it quantified the importance of each middle layer to the output. (Foerster et al., 2017) introduced input-switched affine transformations into RNNs, which analyzed the contribution of input steps via linear methods. Above work focuses on global temporal importance and does not support variable specific temporal interpretation. Our following proposed IMV-LSTM is a combination of attention and post-analyzing for multi-variable time series, where novel variable-wise hidden states and mixture attention enable fine-grained temporal and variable importance interpretation during the training.

Another line of related research is about tensorization and selectively updating of hidden states in RNNs. (Do et al., 2017; Novikov et al., 2015) proposed to represent hidden states as a matrix. (He et al., 2017) developed tensorized LSTM to enhance the capacity of networks without additional parameters. (Kuchaiev & Ginsburg, 2017; Neil et al., 2016; Koutnik et al., 2014) put forward to partition the hidden layer into separated modules with different updates. The hidden state tensors and update processes in existing works do not maintain variable-wise correspondence, thereby lacking the desirable interpretability.

3. INTERPRETABLE MULTI-VARIABLE LSTM:

Assume we have N -1 exogenous time series and a target series y of length T , where y = [y 1 , • • • , y T ] and y ∈ R T .1 By stacking exogenous time series and target series, we define a multi-variable input series as

X T = {x 1 , • • • , x T }, where x t = [x 1 t , • • • , x N -1 t , y t ].

Both of x n t and y t can be multidimensional vector. x t ∈ R N is the multi-variable input at time step t. Given X T , we aim to learn a non-linear mapping to predict the next values of the target series, namely ŷT +1 = F(X T ). Meanwhile, through trained model F(•), we aim to extract variable-wise temporal importance and overall variable importance w.r.t. the prediction from the data. The following described IMV-LSTM can be easily extended to multi-step ahead prediction via iterative methods as well as vector regression (Fox et al., 2018; Cheng et al., 2006) .

3.1. NETWORK ARCHITECTURE:

In IMV-LSTM we develop hidden state matrix and update scheme, which ensure that each element of the hidden state matrix encapsulates information exclusively from a certain variable of the input. It gives rise to mixture attention on both variable and temporal aspects and fine-grained interpretation described below. To distinguish from the hidden state and gate vectors in a standard LSTM, hidden state and gate matrices in IMV-LSTM are denoted with tildes. Specifically, we define the hidden state matrix at time step

t as ht = [h 1 t , • • • , h N t ]

, where ht ∈ R N ×d , h n t ∈ R d . The overall size of the layer is derived as D = N • d. The element h n t of ht is the hidden state vector specific to n-th input variable. Then, we define the input-to-hidden transition as

U j = [U 1 j , • • • , U N j ]

, where U j ∈ R N ×d×d0 , U n j ∈ R d×d0 and d 0 is the dimension of individual variables at each time step. The hidden-to-hidden transition is defined as:

W j = [W 1 j , • • • , W N j ]

, where W j ∈ R N ×d×d and W n j ∈ R d×d . As standard LSTM neural networks (Hochreiter & Schmidhuber, 1997) , IMV-LSTM has the input i t , forget f t , output gates o t and the memory cells c t in the update process. Given the newly incoming input x t at time t and the hidden state matrix ht-1 , the hidden state update is defined as:

EQUATION

where jt = [ j 1 t , • • • , j N t ] has the same shape as hidden state matrix (i.e. R N ×d ). Each element j n t ∈ R d corresponds to the update of the hidden state w.r.t. input variable n. Term W j ht-1 and U j x t respectively capture the update from the hidden states at the previous step and the new input. The tensor-dot operation is defined as the product of two tensors along the N axis, e.g.,

W j ht-1 = [W 1 j h 1 t-1 , • • • , W N j h N t-1 ] where W n j h n t-1 ∈ R d .

Depending on different update schemes of gates and memory cells, we proposed two realizations of IMV-LSTM, i.e. IMV-Full in Equation set 1 and IMV-Tensor in Equation set 2. In these two sets of equations, vec(•) refers to the vectorization operation, which concatenates columns of a matrix into a vector. ⊕ is the concatenation operation. denotes element-wise multiplication. In this paper, matricization(•) reshapes a vector of R D into a matrix of R N ×d .

EQUATION

Equation set 2: IMV-Tensor IMV-Full: With vectorization in Eq. ( 2) and (3), IMV-Full updates gates and memories using full ht-1 and jt regardless of the variable-wise data in them. By simple replacement of the hidden update in standard LSTM by jt , IMV-Full behaves identically to standard LSTM while enjoying the interpretability shown below.

IMV-Tensor: By applying tensor-dot operations in Eq. ( 5), gates and memory cells are matrices as well, elements of which have the correspondence to input variables as hidden state matrix ht does.

W and U have the same shapes as W j and U j in Eq. ( 1)

In IMV-Full and IMV-Tensor, gates only scale jt and ct-1 and thus retain the variable-wise data organization in ht . Meanwhile, based on tensorized hidden state Eq. ( 1) and gate update Eq. ( 5), IMV-Tensor can also be considered as a set of parallel LSTMs, each of which processes one variable series. The derived hidden states specific to individual variables are aggregated on both temporal and variable level through the attention described below.

3.2. MIXTURE ATTENTION:

After feeding a sequence of {x 1 , • • • , x T } into either IMV-Full or IMV-Tensor, we obtain a sequence of hidden state matrices { h1 , • • • , hT }, where the sequence of hidden states specific to variable n is extracted as

{h n 1 , • • • , h n T }.

In this part, we present the novel mixture attention mechanism in IMV-LSTM based on the following idea. Temporal attention is first applied to the sequence of hidden states corresponding to each variable, so as to obtain the summarized history of each variable. Then by using the history enriched hidden state of each variable, global variable attention is derived. These two steps are assembled into a probabilistic mixture model (Zong et al., 2018; Graves, 2013; Bishop, 1994) , which facilitates the subsequent training, inference, and interpretation process.

In particular, the mixture attention is formulated as:

EQUATION

In Eq. ( 8), we introduce a latent random variable z T +1 into the the density function of y T +1 to govern the generation process.

z T +1 is a discrete variable over the set of values {1, • • • , N } corresponding to N input variables. Mathematically, p(y T +1 | z T +1 = n, h n T ⊕ g n ) characterizes the density of y T +1 conditioned on historical data of variable n, while the prior of z T +1 , i.e. p(z T +1 = n | h 1 T ⊕ g 1 , • • • , h N T ⊕ g N ) controls

to what extent y T +1 is driven by variable n. The context vector g n is computed as the temporal attention weighted sum of hidden states corresponding to variable n, i.e., g n = t α n t h n t where attention weight

α n t = exp ( fn(h n t ) ) k exp ( fn(h n k ) ) . f n (•)

can be a flexible function specific to variable n, e.g., neural networks (Bahdanau et al., 2014) . The p(y

T +1 | z T +1 = n, h n T ⊕ g n ) is a Gaussian distribution parameterized by [ µ n , σ n ] = ϕ n ( h n T ⊕ g n ),

where ϕ n (•) can be a feedforward neural network.

The overall variable attention p(z:

T +1 = n | h 1 T ⊕ g 1 , • • • , h N T ⊕ g N ) is derived by a softmax function over {f( h n T ⊕ g n )} N

, where f(•) can be a feedforward neural network shared by all variables.

3.3. LEARNING, INFERENCE, AND INTERPRETATION:

In the learning phase, the set of parameters in IMV-Full or IMV-Tensor as well as the attention functions is denoted by Θ. Given a set of M training sequences {X T } M and {y T +1 } M , the loss function to optimize is defined based on the negative log likelihood of the mixture attention model in Eq. ( 8) plus regularization terms.

In the inference phase, the prediction of y T +1 is obtained by the weighted sum of means as :

ŷT +1 = n µ n • p(z T +1 = n | h 1 T ⊕ g 1 , • • • , h N T ⊕ g N ).

Regarding interpretation, we first illustrate the burden of deciphering variable and temporal importance from the raw attentions mentioned above. For instance, during the training on PLANT dataset used in the experiment section, we collect variable-wise temporal attention and overall variable attention values w.r.t. each training instance at each epoch. In Therefore, we proposed the following summarization method over temporal and variable attentions of data instances to quantify temporal and variable importance. First, we define the function f agg : R A×B → R B , which maps an input Z ∈ R A×B into an aggregated vector z ∈ R B . In the present paper, we choose the simple normalized summation function as:

EQUATION

It is flexible to choose alternative functions for f agg , e.g. robust statistics based methods, etc.

For temporal importance, we collect the temporal attention w.r.t. n-th variable of each instance m as

α n m = [α n m,1 , • • • , α n m,T -1 ] ∈ R T -1 . Then, the unified temporal importance of variable n is obtained as ᾱn = f agg ([α n 1 , • • • , α n M ] ), t ᾱn t = 1, ᾱn t ∈ [0, 1].

The overall variable importance is formulated by using a novel posterior variable attention. Concretely, the posterior variable attention is derived as:

EQUATION

where q n provides more distinguishable attention distribution of different variables by taking the predictive likelihood into account. For instance, Panel (c) in Fig. (2) demonstrates the histograms of posterior and prior attention (i.e. p(z T +1 = n|X T )) of two example variables in PLANT. Compared with priors, the posterior attentions of more important variable "Temperature" further shift rightward, while the posterior of less important variable "P-temperature" moves towards zero.

The unified variable importance can then be defined as q PLANT: This dataset records the time series of energy production of a photo-voltaic power plant in Italy (Ceci et al., 2017) . Exogenous data consists of 9 weather conditions variables (such as temperature, cloud coverage, etc.). The power production is the target. It provides 20842 sequences split into training (70%), validation (10%) and testing sets (20%).

= f agg ([q 1 , • • • , q M ] ), n qn = 1, qn ∈ [0, 1], where q m = [q 1 m , • • • q N m ] ∈ R N is

SML is a public dataset used for indoor temperature forecasting. Same as (Qin et al., 2017) , the room temperature is taken as the target series and another 16 time series are exogenous series. The data were sampled every minute. The first 3200, the following 400 and the last 537 data points are respectively used for training, validation, and test.

4.2. BASELINES AND EVALUATION SETUP:

The first category of statistics baselines includes:

STRX is the structural time series model with exogenous variables (Scott & Varian, 2014; Radinsky et al., 2012) . It is formulated in terms of unobserved components via the state space method.

ARIMAX is the auto-regressive integrated moving average with regression terms on exogenous variables (Hyndman & Athanasopoulos, 2014). It is a special case of vector auto-regression in this scenario.

The second category of machine learning baselines includes:

RF refers to random forests. It is an ensemble learning method consisting of several decision trees (Liaw et al., 2002; Meek et al., 2002) and was used in time series prediction (Patel et al., 2015) .

XGT refers to the extreme gradient boosting (Chen & Guestrin, 2016) . It is the application of boosting methods to regression trees (Friedman, 2001) .

ENET represents Elastic-Net, which is a regularized regression method combining both L1 and L2 penalties of the lasso and ridge methods (Zou & Hastie, 2005 ) and has been used in time series analysis (Liu et al., 2010; Bai & Ng, 2008) .

The third category of neural network baselines includes:

RETAIN uses RNNs to respectively learn weights on input data, which are then used to perform prediction (Choi et al., 2016) .

DUAL is an encoder-decoder architecture, which uses an encoder LSTM to learn weights and then feeds pre-weighted input data into a decoder LSTM for forecasting (Qin et al., 2017) .

In ARIMAX, the orders of auto-regression and moving-average terms are set via the autocorrelation and partial autocorrelation. For RF and XGT, the hyper-parameter tree depth and the number of iterations are chosen from range [3, 10] and [2, 200] via grid search. For XGT, L2 regularization is added by searching within {0.0001, 0.001, 0.01, 0.1, 1, 10}. As for ENET, the coefficients for L2 and L1 penalties are selected from {0, 0.1, 0.3, 0.5, 0.7, 0.9, 1, 1.5, 2}. For machine learning baselines, multi-variable input sequences are flattened into feature vectors. We consider two metrics to measure the prediction performance. Specifically, RMSE is defined as

RMSE = m (y m -ŷm ) 2 /M . MAE is defined as MAE = m |y m -ŷm |/M .

4.3. PREDICTION PERFORMANCE:

We report the prediction errors in Table 1 , each cell of which presents the average RMSE and MAE with standard errors. Note that IMV-Full and IMV-Tensor are single network structures. Their good prediction performance below verifies the idea that instead of complex network architecture in baselines, simple and proper mixture of well-maintained variable-wise hidden states also improves the prediction performance as well as empowering the interpretability shown below. In particular, IMV-LSTM family outperforms baselines by around 80% at most. IMV-Full performs mostly better than baselines, while IMV-Tensor surpasses IMV-Full on NASDAQ and SML datasets. For instance, the importance value at around time lag 1 to 4 of variable "P-temperature" are obviously higher, which indicates that this variable has instant relation to the target. For variable "wind-speed", this heat map tells that its historical observations after time lag 14 could be negligibly correlated to the target.

Panel (c) shows the convergence of the variable importance of IMV-Tensor on PLANT dataset. The ranking of variable importance is slightly different to that in Panel (a). It is because in IMV-Tensor the gate and memory update scheme makes hidden states w.r.t. different variables evolve independently, thereby leading to different hidden states and attention values to IMV-Full. However, we can still decipher something in common. For instance, variable "wind-speed" stays relatively important in both Panel (a) and (c), i.e. rank 1st and 4th respectively. As for temporal importance, in Panel (d) variable "P-temperature" presents temporal importance pattern similar to that in Panel (b). It is also worth exploring whether we can obtain more consistent variable and temporal importance from IMV-Full and IMV-Tensor in the future work. 

4.5. VARIABLE IMPORTANCE FOR PREDICTION:

In this group of experiments, we evaluate the efficacy of variable importance through the lens of prediction tasks. We focus on IMV-LSTM family and RNN baselines, i.e. DUAL and RETAIN. Specifically, for each approach, we first rank variables respectively according to the variable importance in IMV-LSTM and variable attention in DUAL and RETAIN. Then we rebuild datasets only consisting of top 50% ranked variables for each approach (i.e. high importance or attention values) to retrain each model and obtain the prediction errors in Table 2 . (The full ranking of variables is in the appendix.) Ideally, effective variable importance leads to top variables highly related to the target and thus retrained models have comparable errors in comparison to their counterparts in Table 1 . In particular, IMV-Full and IMV-Tensor present comparable and even lower errors in Table 2 , while DUAL and RETAIN have higher errors mostly. An additional advantage of using top variables is the training efficiency, e.g. the training time of each epoch in IMV-Tensor is reduced from ∼16 sec to ∼11 sec. .51 ± 0.12 0.015 ± 0.001, 0.011 ± 0.002 IMV-Tensor 0.12 ± 0.007, 0.10 ± 0.01 157.64 ± 0.14, 128.86 ± 0.13 0.007 ± 0.0005, 0.006 ± 0.0003

5. CONCLUSION:

In this paper, we propose an interpretable multi-variable LSTM (IMV-LSTM) for time series with exogenous variables. Based on the hidden state matrix and update scheme, we present two realizations i.e. IMV-Full and IMV-Tensor as well as developing mixture temporal and variable attention mechanism. It enables to infer and quantify fine-grained variable-wise temporal importance and overall variable importance w.r.t. the target. Extensive experiments exhibit the superior prediction performance, interpretability and efficacy of variable importance of IMV-LSTM.

6. APPENDIX:



In IMV-Full and IMV-Tensor, the learning phase aims to minimize the loss function as follows:

L(Θ) = - M m=1 log N n=1 N (y T +1,m | ϕ n ( h n T,m ⊕ g n m )) exp (f( h n T,m ⊕ g n m )) k exp (f( h k T,m ⊕ g k m )) + λ Θ 2 (11)

The derivation of the posterior attention is as follows. This equation is based on the Bayesian rule to derive the posterior of variable attention, by taking into account both the predictive likelihood w.r.t. the corresponding variable and prior attention.

EQUATION

where the p(y

T +1 | z T +1 = n, h n T ⊕ g n ) is a Gaussian distribution parameterized by [ µ n , σ n ] = ϕ n ( h n T ⊕ g n ),

where ϕ n (•) can be a feedforward neural network. Regarding the summarization of variable-wise temporal attention and variable attention of each instance, since each attention is a sample of discrete distribution (i.e. defined on temporal steps and input variables), alternative summarization method is to fit a Dirichlet distribution to these attentions and use the estimated expectation as the importance values.

Next, we provide the formal analysis about the complexity of IMV-LSTM through Lemma 6.1 and Lemma 6. For IMV-Full, assume each input variable corresponds to one-dimensional time series. Based on Eq. 1, the hidden update has 2D +D 2 /N trainable parameters. Equation set 1 gives rise to the number of parameters equal to that of the standard LSTM. Therefore, the reduce number of parameters is (N -1)D + (1 -1/N )D • D. As for IMV-Tensor, more parameter reduction stems from that the gate update functions in Equation set 2 make use of the tensor-dot operation as Eq. 1. Likewise, the reduced amount of parameters can be derived. As for IMV-Tensor, since the gate update functions in Equation set 2 make use of the tensor-dot operation as Eq. 1, gate update functions have the same computation complexity as Eq. 1. The overall complexity is O(D 2 /N + D), which is 1/N of the complexity of a standard LSTM.

Basically, Lemma 6.1 and Lemma 6.2 indicate that a high number of input variables leads to a large portion of parameter and computation reduction in IMV-LSTM family. In this part, the prediction testing errors of IMV-LSTM family are visualized.

Fig. 4 shows the Q-Q plot of testing errors on each dataset. Fig. 5 demonstrates the testing error curve during the training phase. Basically, we can observe that IMV-Tensor has lower prediction errors as well as more smooth testing error curves. 

:

Vectors are assumed to be in column form throughout this paper.

Code will be released upon requested.

Figures and Tables:

Figure 1: A toy example of a IMV-LSTM with a two-variable input sequence and the hidden matrix of 4-dimensions per variable. Circles represent one dimensional elements. Purple and blue colors correspond to two variables. Blocks containing rectangles with circles inside represent input data and hidden matrix. Panel (a) exhibits the derivation of hidden update jt . Grey areas represent transition weights. Panel (b) demonstrates the mixture attention process. (best viewed in color)

Figure 2: (a) Histograms of overall variable attention at different epochs. (b) Histograms of temporal attention of variable "P-temperature". (c) Prior and posterior attention histograms of two example variables.

the overall variable attention of instance m.NASDAQ is the dataset from(Qin et al., 2017). It contains 81 major corporations under NASDAQ 100, as exogenous time series. The index value of the NASDAQ 100 is the target series. The frequency of the data collection is minute-by-minute. The first 35,100, the following 2,730 and the last 2,730 data points are respectively used as the training, validation and test sets.

IMV-Full: variable-wise temporal importance at different epochs.

IMV-Tensor: variable importance.

IMV-Tensor: variable-wise temporal importance at different epochs.

Figure 3: Variable and temporal importance interpretation during the training of IMV-Full and IMV-Tensor on PLANT dataset. (Best viewed in color)

2. Lemma 6.1. Given time series of N variables, assume a standard LSTM and IMV-LSTM layer both have size D, i.e. D neurons in the layer. Then, compared to the number of parameters of the standard LSTM, IMV-Full and IMV-Tensor respectively reduce the network complexity by (N -1)D + (1 -1/N )D • D and 4(N -1)D + 4(1 -1/N )D • D number of parameters. Proof. In a standard LSTM of layer size D, trainable parameters lie in the hidden and gate update functions. In total, these update functions have 4D•D+4N •D+4D parameters, where 4D•D+4N •D comes from the transition and 4D corresponds to the bias terms.

For time series of N variables and the recurrent layer of size D, IMV-Full and IMV-Tensor respectively have the computation complexity at each update step as: O(D 2 + N • D) and O(D 2 /N + D). Proof. Assume that D neurons of the recurrent layer in IMV-Full and IMV-Tensor are evenly assigned to N input variables, namely each input variable has d = D/N corresponding neurons. For IMV-Full, based on Eq. 1, the hidden update has computation complexity N • d 2 + N • d, while the gate update process has the complexity D 2 + N • D. Overall, the computation complexity is O(D 2 + N • D), which is identical to the complexity of a standard LSTM.

Figure 4: Q-Q plot of testing errors

Figure 5: Testing error curves over the training phase

Variable-wise temporal importance at different epochs.

Figure 8: IMV-Full on NASDAQ dataset. (Best viewed in color)

Variable-wise temporal importance at different epochs.

Figure 9: IMV-Tensor on NASDAQ dataset. (Best viewed in color)

We implemented IMV-LSTM and neural network baselines with Tensorflow 2 . We used Adam with the mini-batch of 64 instances(Kingma & Ba, 2014). For the size of recurrent and dense layers in the baselines, we conduct grid search over{16, 32, 64, 128, 256, 512}. The size of the IMV-LSTM layer is set by the number of neurons per variable selected from {10, 15, 20, 25}. Dropout is selected in {0.8, 0.5}. Learning rate is searched in {0.0005, 0.001, 0.005, 0.01, 0.05}. L2 regularization is added with the coefficient chosen from {0.0001, 0.001, 0.01, 0.1, 1.0}. We train each approach 10 times and report average performance. For baseline DUAL on NASDAQ and SML datasets, we use the hyper-parameters achieving the best performance inQin et al. (2017). On PLANT dataset, hyper-parameters are searched in above sets of values. The window size for NASDAQ and SML, namely T in Sec. 3, is set to 10 according toQin et al. (2017), while for PLANT it is 20 to test long temporal dependency.

Table content: <table/>

RMSE and MAE with std. errors    In this part, we investigate the interpretability of IMV-Full and IMV-Tensor by collecting the variable and temporal importance values during the training under the best hyper-parameters. As far as we know, experiments in previous work using attention in RNNs do not unveil such fine-grained interpretation over both variable and temporal level. Due to the page limitation, the interpretability results on other datasets are in the appendix section.In Fig.3, Panel (a) shows the overall variable importance values w.r.t. training epochs on PLANT dataset. Specifically, as variable importance converges, the ranking of variables is clearly identified at the end of the training, i.e. variables with high importance values are top ranked. Meanwhile, Panel (b) demonstrates the temporal importance values of each variable at the beginning and ending of the training (i.e. epoch 0 and 75). The lighter the color, the higher the temporal importance value of the corresponding time lag. At epoch 0, the randomly initialized network gives rises to similar temporal importance pattern for most of the variables. At epoch 75, diverse patterns w.r.t different variables are learned.

Table content: <table><tr><td>Dataset</td><td>NASDAQ</td><td>PLANT</td><td>SML</td></tr><tr><td>STRX</td><td>0.41 ± 0.01, 0.35 ± 0.02</td><td>231.43 ± 0.19, 193.23 ± 0.43</td><td>0.039 ± 0.001, 0.033 ± 0.001</td></tr><tr><td>ARIMAX</td><td>0.34 ± 0.02, 0.23 ± 0.03</td><td>225.54 ± 0.23, 193.42 ± 0.41</td><td>0.060 ± 0.002, 0.053 ± 0.002</td></tr><tr><td>RF</td><td>0.31 ± 0.02, 0.27 ± 0.03</td><td>164.23 ± 0.65, 130.90 ± 0.15</td><td>0.045 ± 0.001, 0.032 ± 0.001</td></tr><tr><td>XGT</td><td>0.28 ± 0.01, 0.23 ± 0.02</td><td>164.10 ± 0.54, 131.47 ± 0.21</td><td>0.017 ± 0.001, 0.013 ± 0.001</td></tr><tr><td>ENET</td><td>0.31 ± 0.03, 0.21 ± 0.01</td><td>168.22 ± 0.49, 137.04 ± 0.38</td><td>0.018 ± 0.001, 0.015 ± 0.001</td></tr><tr><td>DUAL</td><td>0.31 ± 0.003, 0.21 ± 0.002</td><td>163.29 ± 0.54, 130.87 ± 0.12</td><td>0.019 ± 0.001, 0.015 ± 0.001</td></tr><tr><td>RETAIN</td><td>0.12 ± 0.07 , 0.11 ± 0.06</td><td>250.69 ± 0.36, 190.11 ± 0.15</td><td>0.048 ± 0.001, 0.037 ± 0.001</td></tr><tr><td>IMV-Full</td><td>0.27 ± 0.01, 0.23 ± 0.01</td><td>157.23 ± 0.16,128.13 ± 0.14</td><td>0.015 ± 0.002, 0.012 ± 0.001</td></tr><tr><td colspan="2">IMV-Tensor 0.09 ± 0.005, 0.07 ± 0.004</td><td>159.90 ± 0.22, 129.43 ± 0.10</td><td>0.009 ± 0.0009, 0.006 ± 0.0005</td></tr><tr><td colspan="2">4.4 INTERPRETATION</td><td/><td/></tr></table>

RMSE and MAE with std. errors under top 50% important variables

Table content: <table><tr><td>Dataset</td><td>NASDAQ</td><td>PLANT</td><td>SML</td></tr><tr><td>DUAL</td><td>0.16 ± 0.08, 0.16 ± 0.05</td><td>171.30 ± 0.17, 154.15 ± 0.20</td><td>0.026 ± 0.002, 0.018 ± 0.002</td></tr><tr><td>RETAIN</td><td>0.17 ± 0.03, 0.15 ± 0.02</td><td>226.38 ± 0.72, 167.90 ± 0.81</td><td>0.060 ± 0.001, 0.044 ± 0.004</td></tr><tr><td>IMV-Full</td><td>0.26 ± 0.01, 0.23 ± 0.02</td><td>162.14 ± 0.10, 128</td><td/></tr></table>

RMSE and MAE with std. errors

Table content: <table><tr><td>Dataset</td><td>NASDAQ</td><td>PLANT</td><td>SML</td></tr><tr><td/><td/><td>Without the target variable</td><td/></tr><tr><td>IMV-Full</td><td>0.32 ± 0.01, 0.29 ± 0.01</td><td>163.47 ± 0.21, 135.53 ± 0.12</td><td>0.019 ± 0.001, 0.015 ± 0.002</td></tr><tr><td>IMV-Tensor</td><td>0.11 ± 0.01, 0.09 ± 0.01</td><td>162.10 ± 0.24, 134.99 ± 0.14</td><td>0.011 ± 0.001, 0.010 ± 0.002</td></tr><tr><td/><td colspan="2">Both target and exogenous variables</td><td/></tr><tr><td>IMV-Full</td><td>0.27 ± 0.01, 0.23 ± 0.01</td><td>157.23 ± 0.16,128.13 ± 0.14</td><td>0.015 ± 0.002, 0.012 ± 0.001</td></tr><tr><td colspan="2">IMV-Tensor 0.09 ± 0.005, 0.07 ± 0.004</td><td>159.90 ± 0.22, 129.43 ± 0.10</td><td>0.009 ± 0.0009, 0.006 ± 0.0005</td></tr></table>

Variable importance w.r.t. epochs.

Table content: <table><tr><td/><td>0.200</td><td/><td/><td>AAL AAPL</td><td>CELG CERN</td><td>FOX FOXA</td><td>MXIM MYL</td><td>SWKS SYMC</td></tr><tr><td>Variable importance</td><td>175</td><td/><td/><td>ADBE ADI ADP ADSK AKAM ALXN AMAT AMGN AMZN ATVI AVGO BBBY</td><td>CMCSA COST CSCO CSX CTRP CTSH DISCA DISH DLTR EA EBAY ESRX</td><td>GILD GOOGL INTC JD KHC LBTYA LBTYK LRCX MAR MAT MCHP MDLZ</td><td>NCLH NFLX NTAP NVDA NXPI PAYX PCAR PYPL QCOM QVCA ROST SBUX</td><td>TMUS TRIP TSCO TSLA TXN VIAB VOD VRTX WBA WDC WFM XLNX</td></tr><tr><td/><td/><td/><td/><td>BIDU</td><td>EXPE</td><td>MSFT</td><td>SIRI</td><td>YHOO</td></tr><tr><td/><td/><td/><td/><td>BIIB</td><td>FAST</td><td>MU</td><td>STX</td><td>Auto-regressive</td></tr><tr><td/><td/><td/><td/><td>CA</td><td>FB</td><td/></tr><tr><td colspan="2">AAL AAPL ADBE ADI ADP ADSK AKAM ALXN AMAT AMGN AMZN ATVI AVGO BBBY BIDU BIIB CA CELG CERN CMCSA COST CSCO CSX CTRP CTSH DISCA DISH DLTR EA EBAY ESRX EXPE FAST FB FOX FOXA GILD GOOGL INTC JD KHC LBTYA LBTYK LRCX MAR MAT MCHP MDLZ MSFT MU MXIM MYL NCLH NFLX NTAP NVDA NXPI PAYX PCAR PYPL QCOM QVCA ROST SBUX SIRI STX SWKS SYMC TMUS TRIP TSCO TSLA TXN VIAB VOD VRTX WBA WDC WFM XLNX YHOO Auto-regressive</td><td>8 (a) 9 7 6 5 4 3 Time step lag</td><td>2</td><td>1</td><td/><td/></tr></table>

Variable importance w.r.t. epochs.

Table content: <table><tr><td/><td/><td/><td>AAL</td><td>CELG</td><td>FOX</td><td>MXIM</td><td>SWKS</td></tr><tr><td/><td/><td/><td>AAPL</td><td>CERN</td><td>FOXA</td><td>MYL</td><td>SYMC</td></tr><tr><td/><td/><td/><td>ADBE</td><td>CMCSA</td><td>GILD</td><td>NCLH</td><td>TMUS</td></tr><tr><td/><td/><td/><td>ADI</td><td>COST</td><td>GOOGL</td><td>NFLX</td><td>TRIP</td></tr><tr><td/><td/><td/><td>ADP</td><td>CSCO</td><td>INTC</td><td>NTAP</td><td>TSCO</td></tr><tr><td/><td/><td/><td>ADSK</td><td>CSX</td><td>JD</td><td>NVDA</td><td>TSLA</td></tr><tr><td/><td/><td/><td>AKAM</td><td>CTRP</td><td>KHC</td><td>NXPI</td><td>TXN</td></tr><tr><td/><td/><td/><td>ALXN</td><td>CTSH</td><td>LBTYA</td><td>PAYX</td><td>VIAB</td></tr><tr><td/><td/><td/><td>AMAT</td><td>DISCA</td><td>LBTYK</td><td>PCAR</td><td>VOD</td></tr><tr><td/><td/><td/><td>AMGN</td><td>DISH</td><td>LRCX</td><td>PYPL</td><td>VRTX</td></tr><tr><td/><td/><td/><td>AMZN</td><td>DLTR</td><td>MAR</td><td>QCOM</td><td>WBA</td></tr><tr><td/><td/><td/><td>ATVI</td><td>EA</td><td>MAT</td><td>QVCA</td><td>WDC</td></tr><tr><td/><td/><td/><td>AVGO</td><td>EBAY</td><td>MCHP</td><td>ROST</td><td>WFM</td></tr><tr><td/><td/><td/><td>BBBY</td><td>ESRX</td><td>MDLZ</td><td>SBUX</td><td>XLNX</td></tr><tr><td/><td/><td/><td>BIDU</td><td>EXPE</td><td>MSFT</td><td>SIRI</td><td>YHOO</td></tr><tr><td/><td/><td/><td>BIIB</td><td>FAST</td><td>MU</td><td>STX</td><td>Auto-regressive</td></tr><tr><td/><td/><td/><td>CA</td><td>FB</td></tr><tr><td>AAL AAPL ADBE ADI ADP ADSK AKAM ALXN AMAT AMGN AMZN ATVI AVGO BBBY BIDU BIIB CA CELG CERN CMCSA COST CSCO CSX CTRP CTSH DISCA DISH DLTR EA EBAY ESRX EXPE FAST FB FOX FOXA GILD GOOGL INTC JD KHC LBTYA LBTYK LRCX MAR MAT MCHP MDLZ MSFT MU MXIM MYL NCLH NFLX NTAP NVDA NXPI PAYX PCAR PYPL QCOM QVCA ROST SBUX SIRI STX SWKS SYMC TMUS TRIP TSCO TSLA TXN VIAB VOD VRTX WBA WDC WFM XLNX YHOO Auto-regressive</td><td>8 (a) 9 7 6 5 4 3 Time step lag</td><td>2</td><td>1</td><td/></tr></table>

Variable importance ranking by IMV-Full and IMV-Tensor on NASDAQ dataset.

Table content: <table><tr><td>Dataset</td><td>Method</td><td>Rank of variables according to importance</td></tr><tr><td/><td>IMV-Full</td><td/></tr><tr><td>NASDAQ</td><td/><td/></tr></table>

Peer review (iclr2019):
 The contributions of this paper are in the field of LSTM, where the authors explore the interpretability of LSTM with multivariate data obtained from various and disparate applications. To this end, the authors endow their approach with tensorized hidden states and an update process in order to learn the hidden states. Furthermore, the authors develop a mixture attention mechanism and a summarization methods to quantify the temporal and variable importance in the data. They validate the forecasting and interpretability performance of their approach with experiments. The parer is interesting, well structured and and clearly written. Also, the addressed topic of interpretability is pertinent. However, I have several concerns. 1. In the related work the authors state that “In time series analysis, prediction with exogenous variables is formulated as an auto-regressive exogenous model ” . This is not always right - it is not imperative to add the auto-regressive terms, this is optional and depends on the way we want to formulate our time series forecasting approach and the known constraints. 2. In section 3 — Interpretable Multi-Variable LSTM, by stacking exogenous time series and target series, the authors implicitly formulate their algorithms in a way to consider auto-regression. And I have several concerns with this for time series forecasting. Because, the past is not always a predictor of the future even - particularly in time series context and in industrial settings. And in the occasions where the past allows to predict the future we do not necessarily need to use LSTM to forecast (the notion of persistency in forecasting is enough). Therefore, the power of LSTM in forecasting would have been convincing if you omit the target series in your multi-variable input. 3. In Network Architecture section the authors develop tensorized hidden state and an update scheme. This idea is interesting, I think it would also be good to know what is the algorithmic complexity of this approach? 4. In section 3.3 the authors state that "In the present paper, we choose the simple normalized summation function eq.(9). " Could the authors justify the reason behind this choice? I am not convinced of the reason behind this, especially the authors mention, right after, that "It is flexible to choose alternative functions for f_{agg}" 5. In the experiment section, concerning the prediction performance the authors present a table showing their results, I believe it would have been more compelling to present the prediction results with graphs showing the normalized cumulative errors, as an example. 6. With regard to the interpretation of the results, the authors show the variable importance as a function of the epoch number, it would be equally important to correlate the same figure with the associated prediction results/normalized cumulative errors as a function of the epoch number - this will allow to assess the importance of the interpretability. I think it would be important to further justify the pertinence of this work in terms of interpretability (the statement in the introduction "the interpretability of prediction models is essential for deployment and knowledge extraction" seems to be limited) for example what does it bring knowing the variance importance as a function of the epoch number. As an example, the Pearson correlation coefficient can help select relevant features to a model, and restrict the number of inputs to the relevant ones - can we draw inspiration from this and explain what the authors are proposing in terms of interpretability... Here the idea is to have a motivation presenting the merits of this work, which I think is missing - particularly with the experiments presented here.