STRIC: STACKED RESIDUALS OF INTERPRETABLE COMPO-
Year: 2021

Abstract:
We present a residual-style architecture for interpretable forecasting and anomaly detection in multivariate time series. Our architecture is composed of stacked residual blocks designed to separate components of the signal such as trends, seasonality, and linear dynamics. These are followed by a Temporal Convolutional Network (TCN) that can freely model the remaining components and can aggregate global statistics from different time series as context for the local predictions of each time series. The architecture can be trained end-to-end and automatically adapts to the time scale of the signals. After modeling the signals, we use an anomaly detection system based on the classic CUMSUM algorithm and a variational approximation of the f -divergence to detect both isolated point anomalies and change-points in statistics of the signals. Our method outperforms state-ofthe-art robust statistical methods on typical time series benchmarks where deep networks usually underperform. To further illustrate the general applicability of our method, we show that it can be successfully employed on complex data such as text embeddings of newspaper articles.

text: We present a residual-style architecture for interpretable forecasting and anomaly detection in multivariate time series. Our architecture is composed of stacked residual blocks designed to separate components of the signal such as trends, seasonality, and linear dynamics. These are followed by a Temporal Convolutional Network (TCN) that can freely model the remaining components and can aggregate global statistics from different time series as context for the local predictions of each time series. The architecture can be trained end-to-end and automatically adapts to the time scale of the signals. After modeling the signals, we use an anomaly detection system based on the classic CUMSUM algorithm and a variational approximation of the f -divergence to detect both isolated point anomalies and change-points in statistics of the signals. Our method outperforms state-ofthe-art robust statistical methods on typical time series benchmarks where deep networks usually underperform. To further illustrate the general applicability of our method, we show that it can be successfully employed on complex data such as text embeddings of newspaper articles.

heading: 1 INTRODUCTION
text: Time series data is being generated in increasing volumes from industrial, medical, commercial and scientific applications. Such growth is fueling demand for anomaly detection algorithms that are general enough to be applicable across domains, yet reliable enough to operate on real-world time series data (Munir et al., 2019; Geiger et al., 2020; Su et al., 2019). While recent developments have focused on deep neural networks (DNNs), simple linear models still outperform DNNs in applications that require robustness (Braei & Wagner, 2020) and interpretable failure modes (Geiger et al., 2020; Su et al., 2019).
To harvest the flexibility and interpretability of engineered modules while enabling end-to-end differentiable training, we introduce STRIC: Stacked Residuals of Interpretable Components. We follow standard practice and consider a two stage anomaly detection pipeline comprising a model of the normal time series and an anomaly detector based on the prediction residuals (Munir et al., 2019). In particular, STRIC is composed of three modules: An interpretable local predictor, a non-linear predictor and a novel non-parametric anomaly detector.
More specifically, STRIC uses a parametric model implemented by a sequence of residuals blocks with each layer capturing the prediction residual of previous layers. The first layer models trends, the second layer models quasiperiodicity/seasonality at multiple temporal scales, the third layer is a general linear predictor, and the last is a general non-linear model in the form of a Temporal Convolution Network (TCN). While the first three layers are local (i.e. applied to each component of the time series separately), the last integrates global statistics from additional time series (covariates). Thanks to the residual structure the interpretable linear blocks do not reduce the representative power of our non-linear architecture: any component of the input time-series which cannot be modeled by the interpretable blocks is processed deeper into our architecture by the non-linear module (see Section 4). The model is trained end-toend with a prediction loss and we automatically select its complexity using an upper bound of the marginal likelihood which, to the best of our knowledge, has never been applied to TCNs before.
Anomalies are detected by checking for time instants in which the prediction residual is not stationary. To avoid any unrealistic assumption on the prediction residuals distribution, we use a likelihood ratio test that we derive from a variational upper bound of f -divergences and that can be computed directly from the data points.
To summarize, our contributions are:
1. We introduce STRIC, a stacked residual model that explicitly isolates interpretable factors such as slow trends, quasi-periodicity, and linearly predictable statistics (Oreshkin et al., 2019; Cleveland et al., 1990), and incorporates statistics from other time series as context/side information.
2. We introduce a novel regularization that is added to the prediction loss and is used for automatic model complexity selection according to the Empirical Bayes framework (Rasmussen & Williams, 2006).
3. We introduce a non-parametric extesion of the CUMSUM algorithm which entails a tunable parameter corresponding to the length of observation and enables anomaly detection in the absence of knowledge about the pre- and post-distributions.
4. We test our method on standard anomaly detection benchamrks and show it merges both the advantages of simple and interpretable linear models and the flexibility of non-linear ones while discounting their major drawbacks: lack of flexibility of linear models and lack of interpretability and overfitting of non-linear ones.

heading: 2 RELATED WORK
text: A time series is an ordered sequence of data points. We focus on discrete and regularly spaced time indices, and thus we do not include literature specific to asynchronous time processes in our review. Different methods for time series anomaly detection (TSAD) can be taxonomized by their choice of (i) discriminant function, (ii) continuity criterion, and (iii) optimization method to determine the tolerance threshold. It is common to use statistics of the prediction error as the discriminant (Braei & Wagner, 2020), and the likelihood ratio between the distribution of the prediction error before and after a given time instant as the continuity criterion (Yashchin, 1993). Recent methods compute the discriminant using deep neural network architectures and euclidean distance as continuity criterion (Munir et al., 2019; Geiger et al., 2020; Su et al., 2019; Bashar & Nayak, 2020).
Our method follows a similar line but introduces novel elements both in (i) and (ii): (i) the discriminant function is the prediction error residual of a novel regularized stacked residual architecture; (ii) the decision function is based on a novel non-parametric extension of the CUMSUM algorithm. The resulting method, STRIC, has the advantage of separating interpretable components due to trends and seasonality, without reducing the representative power of our architecture. At initialization, the system is approximately equivalent to a multi-scale SARIMA model (Adhikari & Agrawal, 2013), which can be reliably applied out-of-the-box on most time series. However, as more data is acquired, any part of the system can be further fine-tuned in an unsupervised end-to-end fashion.
Munir et al. (2019) argue that anomaly detection can be solved by exploiting a flexible model provided a proper inductive bias is introduced (e.g. TCN). In Appendix A.7.2 we show that TCN alone might overfit simple time series. We therefore take their direction a step further, and, differently from previous works (Bai et al., 2018; Munir et al., 2019; Sen et al., 2019; Geiger et al., 2020), we provide our temporal model with an interpretable structure, similar to Oreshkin et al. (2019). Moreover, unlike previous works on interpretability of DNNs (Tsang et al., 2018; Guen et al., 2020), our architecture explicitly imposes both an inductive bias and a regularization which are designed to expose the user both a STL-like decomposition (Cleveland et al., 1990) and the relevant time scale of the signals. Since TCNs tend to overfit if not properly regularized (Appendix A.7.2), we constrain our TCN’s representational power by enforcing fading memory (Zancato & Chiuso, 2021) while retaining what is needed to predict future values.
Our method outperforms both classical statistical methods (Braei & Wagner, 2020) and deep networks (Munir et al., 2019; Geiger et al., 2020; Su et al., 2019; Bergman & Hoshen, 2020; Bashar & Nayak, 2020) on different anomaly detection benchmarks (Laptev & Amizadeh, 2020; Lavin & Ahmad, 2015) (Section 6). Moreover, we show it can be employed to detect anomalous patterns on complex data such as text embeddings of newspaper articles (Figure 4).

heading: 3 NOTATION
text: We denote vectors with lower case and matrices with upper case. In particular y is multi-variate time series {y(t)}t∈Z, y(t) ∈ Rn; we stack observations from time t to t+ k − 1 and denote the resulting matrix as Y t+k−1t := [y(t), y(t+ 1), ..., y(t+ k − 1)] ∈ Rn×k. The row index refers to the dimension of the time series while the column index refers to the temporal dimension. We denote the i-th component of the time series y as yi and its evaluation at time t ad as yi(t) ∈ R. We refer to {y(s), s > t} as the test/future and to {y(s), s ≤ t} as the reference/past intervals. At time t, sub-sequences containing the np past samples up to time t − np + 1 are given by Y tt−np+1 (note that we include the
Ytt−np+1 t
̂y(t + 1)
t − np + 1 t + 1
past future
n ⋮
Linear X̂LIN XLIN ∈ ℝn×np Seasonal X̂SEAS XSEAS ∈ ℝn×np
̂ySEAS(t + 1) ∈ ℝ ̂yLIN(t + 1) ∈ ℝ ̂yTCN(t + 1) ∈ ℝ
X̂TREND
XTREND ∈ ℝn×np
TCNTrend
̂yTREND(t + 1) ∈ ℝ
(a) STRIC predictor architecture.
Z = z1 z2 ⋮ zn
Z ∈ {XTREND, XSEAS, XLIN}
Φ = φ1 φ2 ⋮ φn
Φ ∈ {𝒦TREND, 𝒦SEAS, 𝒦LIN}
Layer inputs
Layer filters banks
Layer features selectors
A = aT1 aT2 ⋮ aTn A ∈ {𝒜TREND, 𝒜SEAS, 𝒜LIN}
Time Features
Future Predictor
Gi
aTi Gi
aTi Gibi
ℝl×np∈
φ1 * zi φ2 * zi
⋮ φl * zi
zii = 1,...,n For each
B = bT1 bT2 ⋮ bTn
Layer predictors
B ∈ {ℬTREND, ℬSEAS, ℬLIN}
∈
∈
ℝ1×np
ℝ
(b) Interpretable blocks structure. Time features are extracted independently for each time series (see Appendix A.1 for more details).
present data into the past data), while future samples up to time t + nf are Y t+nf t+1 . We will use past data to predict future ones, where the length of past and future intervals is an hyper-parameter that is up to the user to design.

heading: 4 TEMPORAL RESIDUAL ARCHITECTURE
text: Our architecture is depicted in Figure 1a. Its basic building blocks are causal convolutions (Bai et al., 2018), with a fixed-size 1-D kernel with input elements from time t and earlier. Rather than initializing the convolutional filter randomly, as commonly done in deep learning, we initialize the weights so that each layer is biased to attend at different components of the signal, as explained in the following.
Linear module. The first (linear) module is interpretable and captures local statistics of a given time series by means of a cascade of learnable linear filters. Its first layer models and removes slow-varying components in the input data. We initialize the filters to mimic a causal Hodrick Prescott (HP) filter (Ravn & Uhlig, 2002). The second layer models and removes periodic components: it is initialized to have a periodic impulse response. Finally, the third layer implements a linear stationary filter bank.
We treat the impulse responses parameters of the linear filters as trainable parameters. To allow our model to work on a wide variety of time scales, we initialize the trend layer with different HP smoothness degrees, while we initialize the periodic and linear-stationary layers with randomly chosen poles (Farahmand et al., 2017) both on the unit circle and within the unit circle, thus allowing to capture different periodicities.
Non-linear module. The second (non-linear) module aggregates global statistics from different time series using a TCN model (Sen et al., 2019). It takes as input the prediction residual of the linear module and outputs a matrix G(Y tt−np+1) ∈ R
l×np where l is the number of output features extracted by the TCN model. The columnG(Y tt−np+1)j with j = 1, ..., np of the non-linear features is computed using data up to time t−np + j (due to the internal structure of a TCN network (Bai et al., 2018)). We build a linear predictor on top of G(Y tt−np+1) for each single time series independenty: the predictor for the i-th time series is given by: ŷTCN(t+ 1)i := aTi G(Y t t−np+1)bi where ai ∈ R
l and bi ∈ Rnp . Since ai combines features (uniformly in time) we can interpret it as a feature selector. While bi aggregates relevant features across time indices to build the one-step ahead predictor (see Appendix A.1).
Note that the third layer of the linear module is a superset of preceding ones, and the non-linear module is a superset of the whole linear module. While this makes the model redundant, we show that this design, coupled with proper initialization and regularization, improves the reliability and intepretability of the final model. We improve filters optimization by sharing their kernel parameters among different time series so that global information (e.g., common trend shapes, periodicities) can be extracted. In Appendix A.1, we describe each component of the model in detail.

heading: 4.1 AUTOMATIC COMPLEXITY DETERMINATION
text: Consider the TCN-based future predictor be ŷTCN(t + 1) := aTG(Y tt−np+1)b = X̂TCNb where X̂TCN ∈ R 1×np is the output of the TCN block which depends on the the past window Y tt−np+1 of length np (the memory of the predictor). 1 Ideally, np should be large enough to capture the “true” memory of the time series, but should not be too large if not necessary (i.e. bias variance trade-off). In practice however, flexible feature extractors (such as TCNs or plain DNNs) are prone to overfitting. Therefore, some regularization is needed to control model complexity and benefit from having a large memory window. In this section, we introduce a novel regularized loss inspired by Bayesian arguments which allows us to use an architecture with a “large enough” past horizon np (i.e., larger than the true system memory) and automatically select the relevant past to avoid overfitting. Such information is exposed to the user through an interpretable parameter λ that directly measures the relevant time scale of the signal.
Bayesian learning formulation: We model the innovations (optimal prediction errors) as Gaussian, so that y(t+ 1) | Y tt−np+1 ∼ N (F ∗(Y tt−np+1)), η 2) where F ∗ is the optimal predictor of the future values given the past. Note that this modeling assumption does not restrict our framework and is used only to justify the use of the squared loss to learn the regression function of the predictor. In practice, we do not know F ∗ and we approximate it with our parametric model. For ease of exposition, we group all the architecture parameters except b in W (linear filters parameters, TCN kernel parameters etc.) and write the conditional likelihood of the future given the past data of our parametric model as p(Y t+nft+1 | b,W, Y tt−np+1) = ∏nf k=1 p(y(t+ k) | b,W, Y t+k−1 t+k−np).
To make the notation simpler, we shall denote with Yf := Y t+nf t+1 ∈ Rnf the set of future outputs over which the predictor is computed and we shall use Ŷb,W ∈ Rnf as the predictor’s outputs. In a Bayesian framework, the optimal set of parameters can be found maximizing the posterior p(b,W | Yf ) over the model parameters. We model b and W as independent random variables:
p(b,W | Yf ) ∝ p(Yf | b,W )p(b)p(W ) (1) where p(b) is the prior associated to the predictor coefficients and p(W ) is the prior on the remaining parameters. The prior p(b) should encode our belief that the prediction model should not be too complex and should depend only on the most relevant past. We model this by assuming that the components of b have zero mean and exponentially decaying variances: Eb2np−j−1 = κλ
j for j = 0, ..., np − 1, where κ ∈ R+ and λ ∈ (0, 1). The corresponding maximum entropy prior pλ,κ(b) (Cover & Thomas, 1991) under such constraints is log(pλ,κ(b)) ∝ −‖b‖2Λ−1 − log(|Λ|) where Λ ∈ Rnp is a diagonal matrix with elements Λj,j = κλj with j = 0, ..., np. Here, λ represents how fast the output of the predictor “forgets” the past. Therefore, λ regulates the complexity of the predictor: the smaller λ, the lower the complexity.
In practice, λ has to be estimated from the data. One would be tempted to estimate jointly b,W, λ, κ (and possibly η) by minimizing the negative log of the joint posterior (see Appendix A.2.1). Unfortunately, this leads to a degeneracy since the joint negative log posterior goes to −∞ when λ → 0. Indeed, typically the parameters describing the prior (such as λ) are estimated by maximizing the marginal likelihood, i.e., the likelihood of the data once the parameters (b,W ) have been integrated out. Since computing (or even approximating) the marginal likelihood in this setup is prohibitive, we now introduce a variational upper bound to the marginal likelihood which is easier to estimate.
Variational upper bound to the marginal likelihood: The model structure we consider is linear in b and we can therefore stack the predictions of each available time index t to get the following linear predictor on the whole future data available: Ŷb,W = FW b where F ∈ Rnf×np is obtained by stacking X̂TCN(Y ii−np+1) for i = t, ..., t+ nf − 1.
Proposition 4.1. Consider a model on the form: Ŷb,W = FW b (linear in b and possibly non-linear in W ) and its posterior in Equation (1). Assume the prior on the parameters b is given by the maximum entropy prior and W is fixed. Then the following is an upper bound on the marginal likelihood associated to the posterior in Equation (1) with marginalization taken only w.r.t. b:
Ub,W,Λ = 1
η2 ∥∥∥Yf − Ŷb,W∥∥∥2 + b>Λ−1b+ log det(FWΛF>W + η2I). (2) This (proved in Appendix A.2) provides an alternative loss function to the negative log posterior which does not suffer of the degeneracy alluded above while optimizing over b, W , λ and κ.
1For simplicity we consider scalar time series, but our approach easily generalizes to multivariate time series (Appendix A.2.2).
Remark: We use batch normalization (Ioffe & Szegedy, 2015) to normalize the output of FW along its rows so that features have comparable scales, this avoids the TCN network to counter the fading regularization by increasing its output scales (see Appendix A.2.3).

heading: 5 ANOMALY DETECTOR
text: In this section, we present our anomaly detection method based on a variational approximation of the likelihood ratio between two windows of model residuals. Our temporal residual architecture model produces the prediction residual after removing trends, periodicity, and stationary (linear) components, as well as considering global covariates. Such a prediction residual is used to test the hypothesis that the time instant t is anomalous by comparing its statistics before t on temporal windows of length np and nf . The detector is based on the likelihood ratios aggregated sequentially using the classical CUMSUM algorithm (Page, 1954; Yashchin, 1993). CUMSUM, however, requires knowledge of the distributions, which we do not have. Unfortunately, the problem of estimating the densities is hard (Vapnik, 1998) and generally intractable for high-dimensional time series (Liu et al., 2012). We circumvent this problem by directly estimating the likelihood ratio with a variational characterization of f -divergences (Nguyen et al., 2010) which involves solving a convex risk minimization problem in closed form.
In Section 5.1, we summarize the standard material necessary to derive our new estimator and the resulting anomaly test. The overall method is entirely unsupervised, and users can tune the scale parameter (corresponding to the window of observation when computing the likelihood ratios) and the coefficient of CUMSUM (depending on the application and desired operating point in the tradeoff between missed detection and false alarms).

heading: 5.1 LIKELIHOOD RATIOS AND CUMSUM
text: CUMSUM (Page, 1954) is a classical Sequential Probability Ratio Test (SPRT) (Basseville & Nikiforov, 1993; Liu et al., 2012) of the null hypothesis H0 that the data after the given time c comes from the same distribution as before, against the alternative hypothesis Hc that the distribution is different. We denote the distribution before c as pp and the distribution after the anomaly at time c as pf .
If the density functions pp and pf were known (we shall relax this assumption later), the optimal statistic to decide whether a datum y(i) is more likely to come from one or the other is the likelihood ratio s(y(i)). According to the Neyman-Pearson lemma, H0 is accepted if the likelihood ratio s(y(i)) is less than a threshold chosen by the operator, otherwise Hc is chosen. In our case, the competing hypotheses are H0 = “no anomaly has happened” and Hc = “an anomaly happened at time c”. We denote with pH0 and pHc the p.d.f.s under H0 and Hc so that: pH0(Y K 1 ) = pp(Y K 1 ) and pHc(Y c−1 1 ) = pp(Y c−1 1 ), pHc(Y K c | Y c−11 ) = pf (Y Kc | Y c−1 1 ). Therefore the likelihood ratio is:
Ωtc := pHc(Y
t 1 )
pH0(Y t 1 )
= pp(Y
c−1 1 )pf (Y t c | Y c−11 )
pp(Y t1 ) = pf (Y
t c | Y c−11 )
pp(Y tc | Y c−11 ) = t∏ i=c pf (y(i) | Y i−11 ) pp(y(i) | Y i−11 )
(3)
To determine the presence of an anomaly, we can compute the cumulative sum Stc := log Ω t c of the (log) likelihood ratios, which depends on the time c, and estimate c∗ using a maximum likelihood criterion, corresponding to the detection function ht = max1≤c≤t Stc. The first instant at which we can confidently assess the presence of a change point (a.k.a. stopping time) is: cstop = min{t : ht ≥ τ} where τ is a design parameter that modulates the sensitivity of the detector depending on the application. The final estimate ĉ of the true change point c∗ after the detection cstop is simply given by the timestamp c at which the maximum of ht = max1≤c≤t Stc is achieved. In Appendix A.3, we provide an alternative derivation that shows that CUMSUM is a comparison of the test statistic with an adaptive threshold that keeps complete memory of past ratios. The next step is to relax the assumption of known densities, which we bypass in the next section by directly approximating the likelihood ratios to compute the cumulative sum.

heading: 5.1.1 LIKELIHOOD RATIO ESTIMATION WITH PEARSON DIVERGENCE
text: The goal of this section is to tackle the problem of estimating the likelihood ratio of two general distributions pp and pf given samples. To do so, we leverage a variational approximation of f -divergences (Nguyen et al., 2010) whose optimal solution is directly connected to the likelihood ratio. For different choices of divergence function, different estimators of the likelihood ratio can be built. We focus on a particular divergence choice, the Pearson divergence, since it provides a closed form estimate of the likelihood ratio (see Appendix A.4).
Proposition 5.1. (Nguyen et al., 2010; Liu et al., 2012) Let φ := pf/pp be the likelihood ratio of the unknown distributions pf and pp. Let F := {fi : fi ∼ pf , i = 1, ..., nf} and H := {hi : hi ∼ pp, i = 1, ..., np} be two sets containing nf and np samples i.i.d. from pf and pp respectively. An empirical estimator φ̂ of the likelihood ratio φ is given by the solution to the following convex optimization problem:
φ̂ = arg min φ
1
2np np∑ i=1 φ(hi) 2 − 1 nf nf∑ i=1 φ(fi) (4)
Proposition 5.2. (Liu et al., 2012; Kanamori et al., 2009) Let φ in Equation (4) be chosen in the family of Reproducing Kernel Hilbert Space (RKHS) functions Φ induced by the kernel k. Let the kernel sections be centered on the set of data Str := {F ,H} and let the kernel matrices evaluated on the data from pf and pp be Kf := K(F ,Str) and Kp := K(H,Str). The optimal regularized empirical likelihood ratio estimator on a new datum e is given by:
φ̂(e) = np nf K(e,Str)
( KTp Kp + npγInp+nf )−1 KTf 1. (5)
Remark: The estimator in Equation (5) is not constrained to be positive. Nonetheless, the positivity constraints can be enforced. In this case, the closed form solution is no longer valid but the problem remains convex.

heading: 5.2 SUBSPACE LIKELIHOOD RATIO ESTIMATION AND CUMSUM
text: In this section, we present our anomaly detector estimator. We test for an anomaly in the data Y t1 by looking at the prediction residuals Et1, which provide a sufficient representation of Y t 1 (see Appendix A.5). We therefore assume we are given the prediction errors Et1 obtained from our time series predictor (the predictor should model the normal behaviour). This guarantees that the sequence Et1 is white in each of its normal subsequences. On the other hand, if the model is applied to a data subsequence which contains the abnormal condition, the residuals are correlated.
We estimate the likelihood ratio of pf and pp on a datum et as φ̂t(et). φ̂t is obtained by applying Equation (5) on the past window of size np + nf . At each time instant t, we compute the necessary kernel matrices as Kf (E t t−nf+1, E t t−np−nf+1) and Kp(E t−nf t−np−nf+1, E t t−np−nf+1).
Remark: At time t, the likelihood ratio is estimated assuming i.i.d. data. This assumption holds if no anomaly happened but does not hold in the abnormal situation since residuals are not i.i.d. In Appendix A.5, we prove that treating correlated variables as uncorrelated provides a lower bound on the actual cumulative sum of likelihood ratios. For a fixed threshold, this means the detector cumulates less and therefore requires more time to reach the threshold.
Finally, we compute the detector function by aggregating the estimated likelihood ratios: Ŝtc := ∑t i=c log φ̂i(ei).
Remark: The choice of the windows length (np and nf ) is fundamental and highly influences the likelihood estimator. Using small windows makes the detector highly sensible to point outliers, while larger windows are better suited to estimate sequential outliers (see Appendix A.5).

heading: 6 EXPERIMENTAL RESULTS
text: In this section, we show STRIC can be successfully applied to detect anomalous behaviours on different anomaly detection benchmarks. In particular, we test our novel residual temporal structure, the automatic complexity regularization and the anomaly detector on the following datasets: Yahoo (Laptev & Amizadeh, 2020), NAB (Lavin & Ahmad, 2015), CO2 Dataset (see Appendix A.6). In addition, to show the general applicability and flexibility of our method, we test STRIC on the challenging task of detecting anomalous events in time series generated from embeddings of articles from the New York Times. See Appendix A.7 for the experimental setup and data normalization.
Anomaly detection: While recent works show deep learning models are not well suited to solve AD on standard anomaly detection benchmarks (Braei & Wagner, 2020), we prove deep models can be effective provided they are used with a proper inductive bias and regularization. In Table 1, we compare STRIC against statistical and deep learning based anomaly detection methods. Our experiments follow the experimental setup and evaluation criteria used in Braei & Wagner (2020) and Munir et al. (2019). Note no other method performs consistently (across different datasets) as good as STRIC. In particular, STRIC achieves the best F1 score on Yahoo A3, in Appendix A.7 we show this is mainly
due to STRIC’s predictor. In fact, most of the time series in Yahoo A3 are characterized by trend components and seasonalities which STRIC’s interpretable predictor can easily model (see Appendix A.7.2). In Appendix A.7, we show some ablation studies on the effects of STRIC’s hyper-parameters on its performance. In particular, we find that STRIC is highly affected by the choice of the length of the windows used to estimate the likelihood ratio, while not being much sensitive to the choice of the memory of the predictor (Appendix A.7.1). Interestingly, STRIC does not achieve the optimal F1/AUC compared to linear models on Yahoo A4. The ability of linear models to outperform nonlinear ones on Yahoo A4 is known in the literature (e.g. in Geiger et al. (2020) any non-linear model is outperformed by AR/MA models of the proper complexity). The main motivation for this is that modern (non-linear) methods tend to overfit on Yahho A4 and therefore generalization is usually low. Instead, thanks to fading regularization and model architecture, STRIC does not exhibit overfitting despite having larger complexity than SOTA linear models used in A4. To conclude, we believe that STRIC merges both the advantages of simple and interpretable linear models and the flexibility of non-linear ones while discounting their major drawbacks: lack of flexibility of linear models and lack of interpretability and overfitting of non-linear ones (see Appendix A.8 for a more in depth discussion).
STRIC interpretable time series decomposition: In Figure 2, we show STRIC’s interpretable decomposition. We report predicted signals (first row), estimated trends (second row) and seasonalities (third row) for different datasets. For all experiments, we plot both training data (first 40% of each time series) and test data. Note the interpretable components of STRIC generalize outside the training data, thus making STRIC work well on non-stationary time series (e.g. where the trend component is non negligible and typical non linear models overfit, see Appendix A.7.2).
Ablation study: We now compare the prediction performance of a general TCN model with our STRIC method in which we remove the interpretable module and the fading regularization one at the time. In Table 2, we report the test RMSE prediction errors and the RMSE generalization gap (i.e. difference between test and training RMSE prediction errors) for different datasets while keeping all the training parameters the same (e.g. training epochs, learning rates etc.) and model parameters (e.g. nb = 100). The addition of the linear interpretable model before the TCN slightly improves the test error. We note this effect is more visible on A2, A3, A4, mainly due to the non-stationary nature of these datasets and the fact that TCNs do not easily approximate trends (Braei & Wagner, 2020) (we further tested this
in Appendix A.7.1). While STRIC generalization is always better than a standard TCN model and STRIC’s ablated components, we note that applying Fading memory regularization alone to a standard TCN does not always improve generalization (but never decreases it): this highlights that the benefits of combining the linear module and the fading regularization together are not a trivial ‘sum of the parts’. Consider for example Yahoo A1: STRIC achieves 0.62 test error, the best ablated model (TCN + Linear) 0.88, while TCN + Fading does not improve over the baseline TCN. A similar observation holds for the CO2 Dataset. Fading regularization might not be beneficial (nor detrimental) for time series containing purely periodic components which correspond to infinite memory systems (systems with unitary fading coefficient). In such cases the interpretable module is essential in removing the periodicities and providing the regularized non-linear module (TCN + Fading) with an easier to model residual signal. We refer to Figure 2 (first column) for a closer look on a typical time series in CO2 dataset, which contains a periodic component that is captured by the seasonal part of the interpretable model. To conclude, our proposed fading regularization has (on average) a beneficial effect in controlling the complexity of a standard TCN model and reduces its generalization gap (≈ 40% reduction). Moreover, coupling fading regularization with the interpretable module guarantees the best generalization.
Automatic complexity selection: In Figure 3, we test the effects of our automatic complexity selection (fading memory regularization) on STRIC. We compare STRIC with a standard TCN model and STRIC without regularization as the memory of the predictor increases. The test error of STRIC is uniformly smaller than a standard TCN (without interpretable blocks nor fading regularization). Adding interpretable blocks to a standard TCN improves generalization for a fixed memory w.r.t. standard TCN but gets worse (overfitting occurs) as soon as the available past data horizon increases. On the other hand, the generalization gap of STRIC does not deteriorate as the memory of the predictor increases (see Appendix A.7.1 for a comparison with other metrics).
Anomaly detection on the New York Times dataset: We qualitatevly test STRIC on a time series consisting of BERT embeddings (Devlin et al., 2019) of New York Times articles (Sandhaus, 2008) from 2000 to 2007. We set np = nf = 30 days, to be able to detect change-point anomalies that altered the normal distribution of news articles for a prolonged period of time. Without any human annotation, STRIC is able to detect major historical events such as the 9/11 attack, the 2004 Indian Ocean tsunami, and U.S. elections (Figure 4). Note that we do not carry out a quantitative analysis of STRIC’s predictions, as we are not aware of any ground truth or metrics for this benchmark, see for example Rayana & Akoglu (2015). Additional details and comparison with a baseline model built on PCA are given in Appendix A.8.1.

heading: 7 DISCUSSION AND CONCLUSIONS
text: We have shown that our interpretable stacked residual architecture and our unsupervised estimation of the likelihood ratio are well suited to solve AD for multivariate time series data. Unlike purely DNN-based methods (Geiger et al., 2020; Munir et al., 2019; Bashar & Nayak, 2020), STRIC exposes to the user both an interpretable STL-like time series decomposition (Cleveland et al., 1990) and the relevant time scale of the time series. Both the interepretable module and the fading memory regularization are important in building a successful model (see Table 2). In particular, the interpretable module helps STRIC generalize correctly on non-stationary time series on which standard deep models (such as TCNs) may overfit (Braei & Wagner, 2020). Moreover, we showed that our novel fading regularization alone can improve the generalization error of TCNs up to ≈ 40% over standard TCNs, provided the periodic components of the time series are captured by the interpretable module (Section 6). We highlight that the overall computational complexity and memory requirement of our method remains the same as standard TCNs, so that our approach can easily scale to large scale time series datasets.
An anomaly is a time instant: at that time instant, either we receive an isolated observation that is inconsistent with normal operation (outlier measurement), or a discrete change occurs in the mechanism that generates the data (changepoint) that persists beyond that time instant (Geiger et al., 2020; Basseville & Nikiforov, 1993). Our method treats these two phenomena in a unified manner, without the need to differentiate between outliers and setpoint changes, with specialized detectors for each. Once the predictor is built, our method can be used online to detect anomalies soon after occurrence without waiting for the entire data stream to be observed and without requiring any knowledge on the prediction error distribution (nominal and faulty).
Making the non-parametric anomaly detector fully adaptive to the data is an interesting research direction: While our fading window regularizer automatically tunes the predictor’s window length by exploiting the self-supervised nature of the prediction task, methods to automatically tune the detector’s window lengths (an unsupervised problem) are an interesting research direction. Moreover, designing statistically optimal rules to calibrate our detector’s threshold τ depending on the desired operating point in the tradeoff between missed detection and false alarms would further enhance the out-of-the-box robustness of our method.
Reproducibility Statement: Our method has been tested on publicly available datasets: Yahoo (Laptev & Amizadeh, 2020), NAB (Lavin & Ahmad, 2015), CO2 Dataset (see Appendix A.6) and NYT (Sandhaus, 2008). We described the data splitting and the data processing steps in Appendix A.7. We describe the major details regarding the implementation of our novel method in Appendix A.1 while in Appendix A.7 we describe both the model structure and training hyper-parameters we used in the experimental section. We do not include the model structures and hyper-parameters of SOTA methods we used as baselines and refer to related literature (referenced in our work) for the implementation details. To further foster reproducibility we shall make our code available.

heading: A APPENDIX
text: A.1 IMPLEMENTATION
Given two scalar time series x andϕ, we denote their time convolution g as: g(t) := (ϕ∗x)(t) := ∑∞ i=−∞ ϕ(i)x(t−i). We say that g is the causal convolution of x and ϕ if ϕ(t) = 0 for t < 0, so that the output g(t) does not depend on future values of x (w.r.t. current time index t). In the following, we shall give a particular interpretation to the signals x and ϕ: x will be the input signal to a filter parametrized by an impulse response ϕ (kernel of the filter). Note any (causal) convolution is defined by an infinite summation for each given time t. Therefore it is customary, when implementing convolutional filters, to consider a truncated sum of finite length. In practice, this is obtained assuming the filter impulse response is non-zero only in a finite window. The truncation is indeed an approximation of the complete convolution, nonetheless it is possible to prove that the approximation errors incurred due to truncation are guaranteed to be bounded under simple assumptions. To summarize, in the following we shall write g(t) := (ϕ∗x)(t) and mean that the impulse response of the causal filter ϕ is truncated on a window of a given length.

heading: A.1.1 ARCHITECTURE
text: Let Y tt−np+1 ∈ R n×np be the input data to our architecture at time step t (a window of np past time instants). The main blocks of the architecture are defined to encode trend, seasonality, stationary linear and non-linear part. In the following we shall denote each quantity related to a specific layer using either the subscripts {TREND, SEAS, LIN, TCN} or {0, 1, 2, 3}.
We shall denote the input of each block asXk ∈ Rn×np and the output as X̂k ∈ Rn×np for k = 0, 1, 2, 3. The residual architecture we propose is defined by the following: X0 = Y tt−np+1 and Xk = Xk−1− X̂k−1 for k = 1, 2, 3. At each layer we extract lk temporal features from the input Xk. We denote the temporal features extracted from the input of the k-th block as: Gk := Gk(Xk) ∈ Rlk×np . The i-th column of the feature matrix Gk is a feature vector (of size lk) extracted from the input Xk up to time t− np − i. To do so, we use causal convolutions of the input signal Xk with a set of filter banks (Bai et al., 2018).
A.1.2 INTERPRETABLE RESIDUAL TCN ON SCALAR TIME SERIES
Modeling Interpretable blocks: In this section, we shall describe the main design criteria of the linear module. For each interpretable layer (TREND, SEAS, LIN), we convolve the input signal with a filter bank designed to extract specific components of the input.
For example, consider the trend layer, denoting its scalar input time series by x and its output by gTREND. Then gTREND is defined as a multidimensional time series (of dimension lTREND := l0) obtained by stacking l0 time series given by the convolution of x with l0 causal linear filters: ϕTRENDi ∗ x for i = 0, ..., l0 − 1. In other words, gTREND := [ϕTREND1∗x, ϕTREND2∗x, ..., ϕTRENDl1 ∗x]T . We denote the set of linear filters ϕTRENDi for i = 0, ..., l0−1 asKTREND and parametrize each filter in KTREND with its truncated impulse response (i.e. kernel) of length k0 := kTREND. We interpret each time series in gTREND as an approximation of the trend component of x computed with the i-th filter. We design each ϕTRENDi so that each filter extracts the trend of the input signal on different time scales (Ravn & Uhlig, 2002) (i.e., each filter outputs a signal with a different smoothness degree). We estimate the trend of the input signal by recombining the extracted trend components in gTREND with the linear map aTREND. Moreover, we predict the future trend of the input signal (on the next time-stamp) with the linear map bTREND.
We construct the blocks that extract seasonality and linear part in a similar way.
Implementing Interpretable blocks: The input of each layer is given by a window of measurements of length np. We zero-pad the input signal so that the convolution of the input signal with the i-th filter is a signal of length np (note this introduces a spurious transient whose length is the length of the filter kernel). We therefore have the following temporal feature matrices: G0 = GTREND ∈ Rl0×np , G1 = GSEAS ∈ Rl1×np and G2 = GLIN ∈ Rl2×np .
The output of each layer X̂k is an estimate of the trend, seasonal or stationary linear component of the input signal on the past interval of length np, so that we have X̂k ∈ R1×np (same dimension as the input Xk). On the other hand, the linear predictor ŷk computed at each layer is a scalar. Intuitively, X̂k and ŷk should be considered as the best linear approximation of the trend, seasonality or linear part given block’s filter bank in the past and future. Our architecture performs the following computations: X̂k := aTkGk and ŷk := X̂kbk for k = 0, 1, 2 where ai ∈ Rlk and bk ∈ Rnp . Note ak combines features (uniformly in time) so that we can interpret it as a feature selector while bk aggregates relevant features across different time indices to build the one-step ahead predictor. Depending on the time scale of the signals it is possible to choose bk depending on the time index (similarly to the fading memory regularization). We experimented both the choice to make bk canonical “vectors” and dense vectors. We found that choosing bk as canonical vectors, whose non-zero entry is associated to the closest to the present time instat provides good empirical results on most cases.
Non-linear module The non-linear module is based on a standard TCN network. Its input is defined as X3 = Y tt−np+1 − X̂0 − X̂1 − X̂2, which is to be considered as a signal whose linearly predictable component has been removed. The TCN extracts a set of l3 non-linear features G3(X3) ∈ Rl3×np which we combine with linear maps as done for the previous layers. The j-th column of the non-linear features G3 is computed using data up to time t − np + j (due to the internal structure of a TCN network (Bai et al., 2018)). The linear predictor on top of G3 is ŷTCN := a T 3 G3b3, where a3 ∈ Rl3 and b3 ∈ Rnp .
Finally, the output of our time model is given by:
ŷ(t+ 1) = 3∑ k=0 ŷk = 3∑ k=0 X̂kbk = 3∑ k=0 aTkGk(Xk)bk.
A.1.3 INTERPRETABLE RESIDUAL TCN ON MULTI-DIMENSIONAL TIME SERIES
We extend our architecture to multi-dimensional time series according to the following principles: preserve interpretability (first module) and exploit global information to make local predictions (second module).
In this section, the input data to our model is Y tt−np+1 ∈ R n×np (a window of length np from an n-dimensional time series).
Interpretable module: Each time series undergoes the sequence of 3 interpretable blocks independently from other time series: the filter banks are applied to each time series independently. Therefore, each time series is processed by
the same filter banks: KTREND, KSEAS and KLIN. For ease of notation we shall now focus only on the trend layer. Any other layer is obtained by substituting ‘TREND’ with the proper subscript (‘SEAS’ or ‘LIN’).
We denote by GTREND, i ∈ Rl0×np the set of time features extracted by the trend filter bank KTREND from the i-th time series. Each feature matrix is then combined as done in the scalar setting using linear maps, which we now index by the time series index i: aTRENDi and bTRENDi. The rationale behind this choice is that each time series can exploit differently the extracted features. For instance, slow time series might need a different filter than faster ones (chosen using aTRENDi) or might need to look at values further in the past (retrieved using bTRENDi). We stack the combination vectors aTRENDi and bTRENDi into the following matrices: ATREND = [aTREND1, aTREND2, ..., aTRENDn]T ∈ Rn×l0 and BTREND = [bTREND1, bTREND2, ..., bTRENDn]T ∈ Rn×np . Non-linear module: The second (non-linear) module aggregates global statistics from different time series (Sen et al., 2019) using a TCN model. It takes as input the prediction residual of the linear module and outputs a matrix GTREND(Y t t−np+1) ∈ R
l3×np where l3 is the number of output features that are extracted by the TCN model (which is a design parameter). The j-th column of the non-linear features GTREND(Y tt−np+1) is computed using data up to time t − p + j, where p is the “receptive” field of the TCN (p < np). This is due to the internal structure of a TCN network (Bai et al., 2018) which relies on causal convolutions and typically scales as O(2h) where h is the number of TCN hidden layers (the deeper the TCN the longer its receptive field). As done for the time features extracted by the interpretable blocks, we build a linear predictor on top of GTREND(Y tt−np+1) for each single time series independenty: the predictor for the i-th time series is given by: ŷTCN(t+1)i := aTi GTREND(Y t t−np+1)bi where ai ∈ R
l3 and bi ∈ Rnp . We stack the combination vectors aTCNi and bTCNi into the following matrices: ATCN = [aTCN1, aTCN2, ..., aTCNn]T ∈ Rn×l3 and BTCN = [bTCN1, bTCN2, ..., bTCNn]T ∈ Rn×np . Finally, the outputs of the predictor on the i-th time series are given by:
ŷ(t+ 1)i = ∑
k∈{TREND,SEAS,LIN}
ak T i Gkibki + aTCN T i GTCNbTCNi.

heading: A.1.4 BLOCK STRUCTURE AND INITIALIZATION
text: In this section, we shall describe the internal structure and the initialization of each block.
Structure: Each filter is implemented by means of depth-wise causal 1-D convolutions (Bai et al., 2018). We call the tensor containing the k-th block’s kernel parameters Kk ∈ Rlk×Nk , where lk and Nk are the block’s number of filters and block’s kernel size, respectively (without loss of generality, we assume all filters have the same dimension). Each filter (causal 1D-convolution) is parametrized by the values of its impulse response parameters (kernel parameters). When we learn a filter bank, we mean that we optimize over the kernel values for each filter jointly. For multidimensional time series, we apply the filter banks to each time series independently (depth-wise convolution) and improve filter learning by sharing kernel parameters across different time series.
Initialization: The first block (trend) is initialized using l0 causal Hodrick Prescott (HP) filters (Ravn & Uhlig, 2002) of kernel size N0. HP filters are widely used to extract trend components of signals (Ravn & Uhlig, 2002). In general a HP filter is used to obtain from a time series a smoothed curve which is not sensitive to short-term fluctuations and more sensitive to long-term ones (Ravn & Uhlig, 2002). In general, a HP filter is parametrized by a hyper-parameter λHP which defines the regularity of the filtered signal (the higher λHP, the smoother the output signal). We initialize each filter with λHP chosen uniformly in log-scale between 103 and 109. Note the impulse response of these filters decays to zero (i.e., the latest samples from the input time series are the most influential ones). When we learn the optimal set of trend filter banks, we do not consider them parametrized by λHP and search for the optimal λHP. Instead, we optimize over the impulse response parameters of the kernel which we do not assume live in any manifold (e.g., the manifold of HP filters). Since this might lead to optimal filters which are not in the class of HP filters, we impose a regularization which penalizes the distance of the optimal impulse response parameters from their initialization.
The second block (seasonal part) is initialized using l1 periodic kernels which are obtained as linear filters whose poles (i.e., frequencies) are randomly chosen on the unit circle (this guarantees to span a range of different frequencies). Note the impulse responses of these filters do not go to zero (their memory does not fade away). Similarly to the HP filter bank, we do no optimize the filters over frequencies, but rather we optimize them over their impulse response (kernel parameters). This optimization does not preserve the strict periodicity of filters. Therefore, in order to keep the optimal impulse response close to initialization values (purely periodic), we exploit weight regularization by penalizing the distance of the optimal set of kernel values from initialization values.
The third block (stationary linear part) is initialized using l2 randomly chosen linear filters whose poles lie inside the unit circle, as done in (Farahmand et al., 2017). As the number of filters l2 increases, this random filter bank is guaranteed to be a universal approximator of any (stationary) linear system (see (Farahmand et al., 2017) for details).
Remark: This block could approximate any trend and periodic component. However, we assume to have factored out both trend and periodicities in the previous blocks.
The last module (non-linear part) is composed by a randomly initialized TCN model. We employ a TCN model due to its flexibility and capability to model both long-term and short-term non-linear dependencies. As is standard practice, we exploit dilated convolutions to increase the receptive field and make the predictions of the TCN (on the future horizon) depend on the most relevant past (Bai et al., 2018).
Remark: Our architecture provides an interpretable justification of the initialization scheme proposed for TCN in (Sen et al., 2019). In particular our convolutional architecture allows us to handle high-dimensional time series data without a-priori standardization (e.g., trend or seasonality removal).

heading: A.2 AUTOMATIC COMPLEXITY DETERMINATION (FADING MEMORY REGULARIZATION)
text: In this section, we shall introduce a regularization scheme called fading regularization, to constrain TCN representational capabilities.
The output of the TCN model is G(Y tt−np+1) ∈ R l3×np where l3 is the number of output features extracted by the TCN model. The predictor build from TCN features is given by: aTCNTi GTCN(Y t t−np+1)bTCNi, where the predictor bTCNi ∈ Rnp takes as input a linear combination of the TCN features (weighted by aTCNi). The j-th column of the non-linear features G(Y tt−np+1) is computed using data up to time t− np + j (due to causal convolutions used in the internal structure of the TCN network (Bai et al., 2018)). One expects that the influence on the TCN predictor as j increases should increase too (in case j = np, the statistic is the one computed on the closest window of time w.r.t. present time stamp). Clearly, the exact relevance on the output is not known a priori and needs to be estimated. In other words, the predictor should be less sensitive to statistics (features) computed on a far past, a property which is commonly known as fading memory. Currently, this property is not built in the predictor bTCNi, which treats each time instant equally and might overfit while trying to explain the future by looking into far and possibly non-relevant past. In order to constrain model complexity and reduce overfitting, we impose the fading memory property on our predictor by employing a specific regularization which we now describe.

heading: A.2.1 FADING MEMORY IN SCALAR TIME SERIES
text: We now follow the same notation and assumptions used in Section 4.1 which we now repeat for completeness.
We consider a scalar time series so that the TCN-based future predictor given the past np measures can be written as: ŷTCN(t + 1) = aTGTCN(T tt−np+1)b = X̂kb. We shall assume that innovations (optimal prediction errors) are Gaussian, so that y(t + 1) | Y tt−np+1 ∼ N (F ∗(Y tt−np+1)), η 2), where F ∗ is the optimal predictor of the future values given the past. Note that this assumption does not restrict our framework and is used only to justify the use of the squared loss to learn the regression function of the predictor. In practice, we do not know the optimal F ∗ and we approximate it with our parametric model. For ease of exposition, we group all the architecture parameters except b in the weight vector W (linear filters parameters KTREND,KSEAS,KLIN, linear module recombination weights ATREND,ASEAS,ALIN,BTREND,BSEAS,BLIN, and TCN kernel parameters and recombination coefficients ATCN etc.). We write the conditional likelihood of the future given the past data of our parametric model as:
p(Y t+nf t+1 | b,W, Y tt−np+1) = nf∏ k=1 p(y(t+ k) | b,W, Y t+k−1t+k−np) (6)
To make the notation simpler, we shall denote by Yf := Y t+nf t+1 ∈ Rnf the set of future outputs over which the predictor is computed and we shall use Ŷb,W ∈ Rnf as the predictor’s outputs. Moreover, we shall drop the dependency on the conditioning past Y tt−np+1 (which is present in any conditional distribution). Equation (6) becomes: p(Yf | b,W ) =∏nf k=1 p(y(t + k) | b,W ). The optimal set of parameters b∗ and W ∗ in a Bayesian framework is computed by maximizing the posterior on the parameters given the data:
p(b,W | Yf ) ∝ p(Yf | b,W )p(b)p(W ) (7)
where p(b) is the prior on the predictor and p(W ) is the prior on the remaining parameters. We encode in p(b) our prior belief that the complexity of the predictor should not be too high and therefore it should only depend on the most relevant past.
Remark: The prior does not induce hard constraints. It rather biases the optimal predictor coefficients towards the prior belief. This is clear by looking at the negative log-posterior which can be directly interpreted as the loss function to be minimized: − log p(b,W | Yf ) = − log p(Yf | b,W ) − log p(b) − log p(W ). In particular, the first term log p(Yf | b,W ) is the data fitting term (only influenced by the data). Both log p(b) and log p(W ) do not depend on the available data and can be interpreted as regularization terms that bound the complexity of the predictor function.
The main idea is to reduce the sensitivity of the predictor on time instants that are far in the past. We therefore enforce the fading memory assumption on p(b) by assuming that the components of b ∈ Rnp have zero mean and exponentially decaying variances:
Ebj = 0 and Eb2np−j−1 = κλ j for j = 0, ..., np − 1 (8)
where κ ∈ R+ and λ ∈ (0, 1). Note the larger variance (larger scale) is associated to temporal indices close to the present time t.
Remark: To specify the prior, we need a density function p(b) but up to now we only specified constraints on the first and second order moments. We therefore need to constrain the parametric family of prior distributions we consider. Any choice on the class of prior distributions lead to different optimal estimators. Among all the possible choices of prior families we choose the maximum entropy prior (Cover & Thomas, 1991). Under constraints on first and second moment, the maximum entropy family of priors is the exponential family (Cover & Thomas, 1991). In our setting, we can write it as: log pλ,κ(b) ∝ −‖b‖2Λ−1 − log |Λ| (9) where Λ ∈ Rnp×np is a diagonal matrix whose elements are Λj,j = κλj for j = 0, ..., np − 1. The parameter λ represents how fast the predictor’s output ‘forgets’ the past: the smaller λ, the lower the complexity. In practice, we do not have access to this information and indeed we need to estimate λ from the data.
One would be tempted to estimate jointly W, b, λ, κ (and possibly η) by minimizing the negative log of the joint posterior:
arg min b,W,λ,κ
1
η2 ∥∥∥Yf − Ŷb,W∥∥∥2 + log(η2)− log(pλ,κ(B))− log(p(W )). (10) Unfortunately, this leads to a degeneracy since the joint negative log posterior goes to −∞ when λ→ 0. Bayesian learning formulation for fading memory regularization:
The parameters describing the prior (such as λ) are typically estimated by maximizing the marginal likelihood, i.e., the likelihood of the data once the parameters (b,W ) have been integrated out. Unfortunately, the task of computing (or even approximating) the marginal likelihood in this setup is prohibitive and one would need to resort to Monte Carlo sampling techniques. While this is an avenue worth investigating, we preferred to adopt the following variational strategy inspired by the linear setup.
Indeed, the model structure we consider is linear in b and we can therefore stack the predictions of each available time index t to get the following linear predictor on the whole future data: Ŷb,W = FW b where FW ∈ Rnf×np and its rows are given by X̂TCN(Y ii−np+1) for i = t, ..., t+ nf − 1.
We are now ready to find an upper bound to the marginal likelihood associated to the posterior given by Equation (7) with marginalization taken only w.r.t. b. Proposition A.1 (from (Tipping, 2001)). The optimal value of a regularized linear least squares problem with feature matrix F and parameters b is given by the following equation:
arg min b
1
η2 ‖Yf − Fb‖2 + b>Λ−1b = Y >f Σ−1Yf (11)
with Σ := FΛF>AT + η2I .
Equation (11) guarantees that 1
η2 ‖Yf − Fb‖2 + b>Λ−1b+ log |Σ| ≥ Y >f Σ−1Yf + log |Σ|,
where the right hand side is (proportional to) the negative marginal likelihood with marginalization taken only w.r.t. b. Therefore, for fixed a W ,
1
η2 ∥∥∥Yf − Ŷb,W∥∥∥2 + b>Λ−1b+ log |FWΛF>W + η2I| is an upper bound of the marginal likelihood with marginalization over b and does not suffer of the degeneracy alluded at before.
With this considerations in mind, and inserting back the optimization over W , the overall optimization problem we solve is
arg min b,W,λ∈(0,1),κ>0
1
η2 ∥∥∥Yf − Ŷb,W∥∥∥2 + ‖b‖2Λ−1 + log |FWΛFW + η2I|+ log p(W ) (12) Remark: log p(W ) defines the regularization applied on the remaining parameters of our architecture. In particular, we induce sparsity by applying L1 regularization on ATREND, ASEAS, ALIN and ATCN. Also, we constrain filters parameters to stay close to initialization by applying L2 regularization on KTREND, KSEAS and KLIN.

heading: A.2.2 FADING MEMORY IN MULTIVARIATE TIME SERIES
text: In the case of multivariate time series, fading regularization can be applied either with a single fading coefficient λ for all the time series or with different fading coefficients for each time series. In all the experiments in this paper, we chose to keep one single λ for all the time series. In practice, this choice is sub-optimal and might lead to more overfitting than treating each time series separately: the ‘dominant’ (slower) time series will highly influence the optimal λ.

heading: A.2.3 FEATURES NORMALIZATION
text: We avoid the non-identifiability of the product FW b by exploiting batch normalization: we impose that different features have comparable means and scales across time indices i = 0, ..., np − 1. Non-identifiability occurs due to the product FW b, if features have different scales across time indices (i.e., columns of the matrix FW ) the benefit of fading regularization might reduced since it can happen that features associated with small bi have large scale so that the overall contribution of the past does not fade. Hence we use batch normalization to normalize time features. Then we use an affine transformation (with parameters to be optimized) to jointly re-scale all the output blocks before the linear combination with b.

heading: A.3 ALTERNATIVE CUMSUM DERIVATION AND INTERPRETATION
text: In this section, we describe an equivalent formulation of the CUMSUM algorithm we derived in the main paper. Before a change point, by construction we are under the distribution of the past. Therefore, log pf (y)pp(y) ≤ 0 ∀y, which in turn means that the cumulative sum St1 will decrease as t increases (negative drift). After the change, the situation is opposite and the cumulative sum starts to show a positive drift, since we are sampling y(i) from the future distribution pf . This intuitive behaviour shows that the relevant information to detect a change point can be obtained directly from the cumulative sum (along timestamps). In particular, all we need to know is the difference between the value of the cumulative sum of log-likelihood ratios and its minimum value.
The CUMSUM algorithm can be expressed using the following equations: vt := St1 − mt, where mt := minj,1≤j≤t S t j . The stopping time is defined as: tstop = min{t : vt ≥ τ} = min{t : St1 ≥ mt + τ}. With the last equation, it becomes clear that the CUMSUM detection equation is simply a comparison of the cumulative sum of the log likelihood ratios along time with an adaptive thresholdmt+τ . Note that the adaptive threshold keeps complete memory of the past ratios. The two formulations are equivalent because St1 −mt = ht.
A.4 VARIATIONAL APPROXIMATION OF THE LIKELIHOOD RATIO
In this section, we present some well known facts on f -divergences and their variational characterization. Most of the material and the notation is from (Nguyen et al., 2010). Given a probability distribution P and a random variable f measurable w.r.t. P, we use ∫ fdP to denote the expectation of f under P. Given samples x(1), ..., x(n) from P, the
empirical distribution Pn is given by Pn = 1n ∑n i=1 δx(i). We use ∫ fdPn as a convenient shorthand for the empirical
expectation 1n ∑n i=1 f(x(i)).
Consider two probability distributions P and Q, with P absolutely continuous w.r.t. Q. Assume moreover that both distributions are absolutely continuous with respect to the Lebesgue measure µ, with densities p0 and q0, respectively, on some compact domain X ⊂ Rd. Variational approximation of the f-divergence: The f -divergence between P and Q is defined as (Nguyen et al., 2010)
Df (P,Q) := ∫ p0f ( q0 p0 ) dµ (13)
where f : R → R is a convex and lower semi-continuous function. Different choices of f result in a variety of divergences that play important roles in various fields (Nguyen et al., 2010). Equation (13) is usually replaced by the variational lower bound:
Df (P,Q) ≥ sup φ∈Φ
∫ [φdQ− f∗(φ)dP] (14)
and equality holds iff the subdifferential ∂f( q0p0 ) contains an element of Φ. Here f ∗ is defined as the convex dual function of f .
In the following, we are interested in divergences whose conjugate dual function is smooth (which in turn defines commonly used divergence measures such as KL and Pearson divergence), so that we shall assume that f is convex and differentiable. Under this assumption, the notion of subdifferential is not required and the previous statement reads as: equality holds iff ∂f( q0p0 ) = φ for some φ ∈ Φ.
Remark: The infinite-dimensional optimization problem in Equation (14) can be written as Df (P,Q) = supφ∈Φ EQφ− EPf∗(φ). In practice, one can have an estimator of any f -divergence restricted to a functional class Φ by solving Equation (14) (Nguyen et al., 2010). Moreover, when P and Q are not known one can approximate them using their empirical counterparts: Pn and Qn. Then an empirical estimate of the f -divergence is: D̂f (P,Q) = supφ∈Φ EQnφ−EPnf∗(φ). Approximation of the likelihood ratio: An estimate of the likelihood ratio can be directly obtained from the variational approximation of f -divergences. The key observation is the following: equality on Equation (14) is achieved iff φ = ∂f( q0p0 ). This tells us that the optimal solution to the variational approximation provides us with an estimator of the composite function ∂f( q0p0 ) of the likelihood ratio q0 p0
. As long as we can invert ∂f , we can uniquely determine the likelihood ratio.
In the following, we shall get an empirical estimator of the likelihood ratio in two separate steps. We first solve the following:
φ̂n := arg max φ∈Φ
EQnφ− EPnf∗(φ) (15)
which returns an estimator of ∂f( q0p0 ), not the ratio itself. And then we apply the inverse of ∂f to φ̂n. We therefore have a family of estimation methods for the likelihood function by simply ranging over choices of f .
Remark: If f is not differentiable, then we cannot invert ∂f but we can obtain estimators of other functions of the likelihood ratio. For instance, we can obtain an estimate of the thresholded likelihood ratio by using a convex function whose subgradient is the sign function centered at 1.

heading: A.4.1 LIKELIHOOD RATIO ESTIMATION WITH PEARSON DIVERGENCE
text: In this section, we show how to estimate the likelihood ratio when the Pearson divergence is used. With this choice, many computations simplify and we can write the estimator of the likelihood ratio in closed form. Other choices (such as the Kullback-Leibler divergence) are possible and legitimate, but usually do not lead to closed form expressions (see (Nguyen et al., 2010)).
The Pearson, or χ2, divergence is defined by the following choice: f(t) := (t−1) 2
2 . The associated convex dual function is :
f∗(v) = sup u∈R
{ uv − (u− 1) 2
2
} = v2
2 + v.
Therefore the Pearson divergence is characterized by the following: PE(P||Q) := ∫ p0 ( q0 p0 − 1 )2 dµ ≥ sup φ∈Φ EQφ− 1 2 EPφ2 − EPφ. (16)
Solving the lower bound for the optimal φ provides us an estimator of ∂f( q0p0 ) = q0 p0 − 1. For the special case of the Pearson divergence, we can apply a change of variables which preserves convexity of the variational optimization problem Equation (16) and provides a more straightforward interpretation. Let the new variable be z := φ + 1 with z ∈ Z , which in this case is nothing but the inverse function of ∂f . We get
sup φ∈Φ
EQφ− 1
2 EPφ2 − EPφ = sup z∈Z EQz −
1 2 EPz2 − 1 2 (17)
It is now trivial to see that z is a ‘direct’ approximator of the likelihood ratio (i.e., it does not estimate a composite map of the likelihood ratio). Therefore for simplicity, we shall employ
arg min φ∈Φ
1 2 EPφ2 − EQφ (18)
to build our ‘direct’ estimator of the likelihood ratio.
Let the samples from P and Q be, respectively, xp(i) with i = 1, ..., np and xq(i) with i = 1, ..., nq . We define the empirical estimator of the likelihood ratio φ̂n:
φ̂n = arg min φ∈Φ
1
2np np∑ i=1 φ(xp(i)) 2 − 1 nq nq∑ i=1 φ(xf (i)). (19)
A closed form solution: Up to now we have not defined in which class of functions our approximator φ lives. As done in (Nguyen et al., 2010; Liu et al., 2012), we choose φ ∈ Φ where Φ is a RKHS induced by the kernel k. We exploit the representer theorem to write a general function within Φ as:
φ(x) = ntr∑ i=1 k(x, xtr(i))αi,
where we use ntr data which are the centers of the kernel sections used to approximate the unknown likelihood ratio (how to choose these centers is important and determines the approximation properties of φ̂n). For now, we do not specify which data should be used as centers (we can use either data from Pn or from Qn or from both or simply use user specified locations).
Let us define the following kernel matrices: Kp := K(Xp, Xtr) ∈ Rnp×ntr , Kq := K(Xq, Xtr) ∈ Rnq×ntr , where Xp := {xp(i)}, Xq := {xq(i)} and Xtr := {xtr(i)}. We therefore have:
φ̂n = arg min φ∈Φ
1
2np np∑ i=1 φ(xp(i)) 2 − 1 nq nq∑ i=1 φ(xf (i))
= arg min α,α≥0
1
2np np∑ i=1 ( ntr∑ j=1 k(xp(i), xtr(j))αj) 2 − 1 nf nf∑ i=1 ntr∑ j=1 k(xf (i), xtr(j))αj
= arg min α,α≥0
1
2np αTKTp Kpα−
1
nf 1 TKfα
Remark: We impose the recombination coefficients α to be non negative since the likelihood ratio is a non negative quantity. The resulting optimization problem is a standard convex optimization problem with linear constraints which can be efficiently solved with Newton methods, nonetheless in general it does not admit any closed form solution.
We now relax the positivity constraints so that the optimal solution can be obtained in closed form. Moreover we add a quadratic regularization term as done in (Nguyen et al., 2010) which lead us to the following regularized optimization problem:
arg min α
1
2np αTKTp Kpα−
1
nf 1 TKfα+
γ 2 ‖α‖2Φ
whose solution is trivially given by:
α̂ = np nf
( KTp Kp + npγIntr )−1 KTf 1 :=
np nf H−1KTf 1 (20)
The estimator of the likelihood ratio for an arbitrary location x is given by the following:
pq(x) pp(x) ≈ φ̂n(x) = K(x,Xtr)α̂ = np nf K(x,Xtr)
( KTp Kp + npγIntr )−1 KTf 1 (21)
Remark: In the following we shall exploit RBF kernels which are defined by the length scales σ.

heading: A.5 SUBSPACE LIKELIHOOD RATIO ESTIMATION AND CUMSUM
text: In this section we describe our subspace likelihood ratio estimator and its relation to the CUMSUM algorithm. The CUMSUM algorithm requires to compute the likelihood ratio pf (y(t)|Y t−1 c )
pp(y(t)|Y t−1c ) for each time t. We denote pp as the
normal density and pf as the abnormal one (after the anomaly has occurred).
We shall proceed to express the conditional probability p(y(t) | Y t−11 ) using our predictor. In particular it is always possible to express the optimal (unknown) one-step ahead predictor as:
ŷt|t−1 = F ∗(Y tt−K+1) := E[y(t) | Y tt−K+1] (22)
which is a deterministic function given the past of the time series (whose length is K). So that the data density distribution can be written in innovation form (based on the optimal prediction error) as:
y(t) = F ∗(Y tt−K+1) + e(t) (23)
where e(t) := y(t) − F ∗(Y tt−K+1) is, by definition, the one step ahead prediction error (or innovation sequence) of y(t) given its past. We therefore have: p(y(t) | Y tt−K+1) = p(e(t) | Y tt−K+1). Where e(t) is the optimal prediction error for each time t and is therefore indipendent on each time t.
Remark: In practice we do not know F ∗ and we use our predictor learnt from normal data as a proxy. This implies the prediction residuals are approximately independent on normal data (the predictor can explain data well), while the prediction residuals are, in general, correlated on abnormal data.
To summarize: under normal conditions the joint distribution of Y tc can be written as:
p(Y tc ) = t∏ i=c p(y(i) | Y i−1c ) = t∏ i=c p(e(i)) normal conditions (24)
p(Y tc ) = t∏ i=c p(y(i) | Y i−1c ) = t∏ i=c p(e(i) | Ei−1c ) abnormal conditions (25)
These two conditions in turn influence the log likelihood ratio test as follows: under H0 =⇒ ∏t i=c pf (e(i)) pp(e(i)) while
under Hc =⇒ ∏t i=c pf (e(i)|Ei−1c ) pp(e(i))
. The main issue here is the numerator under Hc: the distribution of residuals changes at each time-stamp (it is a conditional distribution) and pf (e(i) | Ei−1c ) is difficult to approximate (it requires the model of the fault). In the following we show that replacing pf (e(i) | Ei−1c ) with pf (e(i)) allows us to compute a
lower bound on the cumulative sum. Such an approximation is necessary to estimate the likelihood ratio in abnormal conditions, the main downside of this approximation is that the detector becomes slower (it needs more time to reach the stopping time threshold).
Applying the independent likelihood test in a correlated setting: We now show that treating pf (e(i) | Ei−1c ) as independent random variables pf (e(i)) for i = 1, ..., t allows us to compute a lower bound on the log likelihood log Ωtc (i.e. the cumulative sum). We denote the cumulative sum of the log likelihood ratio using independent variables as log Ω̄tc = ∑t i=c log pf (e(i)) pp(e(i)) Proposition A.2. Assume a change happens at time c so that Hc is true and the following log likelihood ratio holds true: log Ωtc = ∑t i=c log pf (e(i)|Ei−1c ) pp(e(i)) . Then it holds log Ωtc ≥ log Ω̄tc.
Proof. By simple algebra we can write: pf (e(i) | Ei−1c )
pp(e(i)) = pf (e(i) | Ei−1c ) pf (e(i)) pf (e(i)) pp(e(i)) ∀i
Now recall the cumulative sum of the log-likelihood ratios taken under the current data generating mechanism pf (Et1) provides an estimate of the expected value of the log-likelihood ratio. Due to the correlated nature of data Et1 the samples are drawn from a multidimensional distribution of dimension t (a sample from this distribution is an entire trajectory from c to t).
We now take the expectation of previous formula w.r.t. the ‘true’ distribution pf (Et1):
Epf (Etc)Ω t c = Epf (Etc) log t∏ i=c pf (e(i) | Ei−1c ) pp(e(i))
= Epf (Etc) log t∏ i=c pf (e(i) | Ei−1c ) pf (e(i)) + Epf (Etc) log t∏ i=c pf (e(i)) pp(e(i))
= MI ( pf (E t c); t∏ i=c pf (e(i)) ) +KL ( t∏ i=c pf (e(i)) ∣∣∣∣∣∣ t∏ i=c pp(e(i)) )
≥ KL ( t∏ i=c pf (e(i)) ∣∣∣∣∣∣ t∏ i=c pp(e(i)) )
where we used the fact the mutual information is always non negative.
How to approximate pre and post fault distributions: Both pp and pf are not known and their likelihood ratio need to be estimated from available data. From Section 5.1.1 we know how to approximate the likelihood ratio given two set of data without estimating the densities. In our anomaly detection setup we define these two sets as: Ep := E t−nf t−np−nf+1 and Ef := E t t−nf+1. So that given current time t we look back at a window of length np + nf . The underlying assumption is that under Hc normal data are present in Ep and abnormal ones in Ef . We estimate the likelihood ratio pf (e(t))pp(e(t)) at each time t by assuming both Ep and Ef data are independent (see Proposition A.2) and cumulate their log as t increases.
How do np and nf affect our detector? The choice of the windows length (np and nf ) is fundamental and highly influences the likelihood estimator. Using small windows makes the detector highly sensible to both point and sequential outliers, while larger windows are better suited to estimate only sequential outliers. We now assume np = nf and study how small and large values affect the behaviour of our detector in simple working conditions.
In Figure 5 and Figure 6 we compute the cumulative sum of log likelihood ratios estimated from data on equally sized windows. Intuitively any local minimum after a ‘large’ (depending on the threshold τ ) increase of the cumulative sum is a candidate abnormal point.
In Figure 7 and Figure 8 we compare the cumulative sum of estimated likelihood ratios on data in which both sequential and point outliers are present. In particular we highlight that large window sizes np and nf are usually not able to capture point anomalies Figure 7 while using small window sizes allow to detect both (at the expenses of a more sensitive detector) Figure 8.

heading: A.6 DATASETS

heading: A.6.1 YAHOO DATASET
text: Yahoo Webscope dataset (Laptev & Amizadeh, 2020) is a publicly available dataset containing 367 real and synthetic time series with point anomalies, contextual anomalies and change points. Each time series contains 1420-1680 time stamps. This dataset is further divided into 4 sub-benchmarks: A1 Benchmark, A2 Benchmark, A3 Benchmark and A4 Benchmark. A1Benchmark is based on the real production traffic to some of the Yahoo! properties. The other 3 benchmarks are based on synthetic time series. A2 and A3 Benchmarks include point outliers, while the A4Benchmark includes change-point anomalies. All benchmarks have labelled anomalies. We use such information only during evaluation phase (since our method is completely unsupervised).

heading: A.6.2 NAB DATASET
text: NAB (Numenta Anomaly Benchmark) (Lavin & Ahmad, 2015) is a publicly available anomaly detection benchmark. It consists of 58 data streams, each with 1,000 - 22000 instances. This dataset contains streaming data from different domains including read traffic, network utilization, on-line advertisement, and internet traffic. As done in (Geiger et al., 2020) we choose a subset of NAB benchmark, in particular we focus on the NAB Traffic and NAB Tweets benchmarks.

heading: A.6.3 CO2 DATASET
text: We test the prediction and interpretability capabilities of our model on the CO2 dataset from kaggle 2. The main goal here is to predict both trend and periodicity of CO2 emission rates on different years. Note this is not an Anomaly detection task.

heading: A.6.4 NYT DATASET
text: The New York Times Annotated Corpus3 (Sandhaus, 2008) contains over 1.8 million articles written and published by the New York Times between January 1, 1987 and June 19, 2007. We pre-processed the lead paragraph of each article with a pre-trained BERT model (Devlin et al., 2019) from the HuggingFace Transformers library (Wolf et al., 2020) and extracted the 768-dimensional hidden state of the [CLS] token (which serves as an article-level embedding). For each day between January 1, 2000 and June 19, 2007, we took the mean of the embeddings of all articles from that day. Finally, we computed a PCA and kept the first 200 principal components (which explain approximately 95% of the variance), thus obtaining a 200-dimensional time series spanning 2727 consecutive days. Note that we did not use any of the dataset’s annotations, contrary to prior work such as Rayana & Akoglu (2015).

heading: A.7 EXPERIMENTAL SETUP
text: In this section, we shall describe the experimental setup we used to test STRIC.
Data preprocessing: Before learning the predictor we standardize each dataset to have zero mean and standard deviation equals to one. As done in (Braei & Wagner, 2020) we note standardization is not equal to normalization, where data are forced to belong to the interval (0, 1). Normalization is more sensitive to outliers, thus it would be inappropriate to normalize our datasets, which contain outliers.
We do not apply any deseasonalizing or detrending pre-processing.
Data splitting: We split each dataset into training and test sets preserving time ordering, so that the first data are used as train set and the following ones are used as test set. The data used to validate the model during optimization are last 10% of the training dataset. Depending on the experiment, we choose a different percentage in splitting train and test. When comparing with (Braei & Wagner, 2020) we used 30% as training data, while when comparing to (Munir et al., 2019) we use 40%. Such a choice is dictated by the particular (non uniform) experimental setup reported in (Braei &
2https://www.kaggle.com/txtrouble/carbon-emissions 3https://catalog.ldc.upenn.edu/LDC2008T19
Wagner, 2020; Munir et al., 2019) and has been chosen to produce comparable results with state of the art methods present in literature.
Evaluation metrics: We compare different predictors by means of the RMSE (root mean squared error) on the onestep ahead prediction errors. Given a sequence of data Y N1 and the one-step ahead predictions Ŷ N 1 the RMSE is
defined as: √
1 N ∑N i=1 ‖y(i)− ŷ(i)‖ 2.
As done in (Braei & Wagner, 2020) we compare different anomaly detection methods taking into account several metrics. We use F1-Score which is defined as the armonic mean of Precision and Recall (see (Braei & Wagner, 2020; Munir et al., 2019)) and another metric that is often used is receiver operating characteristic curve, ROC-Curve, and its associated metric area under the curve (AUC). The AUC is defined as the area under the ROC-Curve. This metric is particularly useful in our anomaly detection setting since it describes with an unique number true positive rate and false positive rate on different threshold values. We now follow (Braei & Wagner, 2020) to describe how AUC is computed. Let the true positive rate and false positive rate be defined, respectively, as: TPR = TPP and FPR = FP N , where TP stands for true positive, P for positive, FP for false positive andN for negative. To copute the ROC-Curve we use different thresholds on our anomaly detection method. We therefore have different pairs of TPR and FPR for each threshold. These values can be plotted on a plot whose x and y axes are, respectively: FPR and TPR. The resulting curve starts at the origin and ends in the point (1,1). The AUC is the area under this curve. In anomaly detection, the AUC expresses the probability that the measured algorithm assigns a random anomalous point in the time series a higher anomaly score than a random normal point.
Hardware: We conduct our experiments on the following hardware setup:
• Processor: Intel(R) Core(TM) i9-10980XE CPU @ 3.00GHz
• RAM: 128 Gb
• GPU: Nvidia TITAN V 12Gb and Nvidia TITAN RTX 24Gb
Hyper-parameters: All the experiments we carried out are uniform on the optimization hyper-parameters. In particular we fixed the maximum number of epochs to 300, the learning rate to 0.001 and batch size to 100. We optimize each model using Adam and early stopping.
We fix STRIC’s first module hyper-parameters as follows:
• number of filter per block: l0 = 10, l1 = 100, l2 = 200
• linear filters kernel lengths (N0, N1, N2): half predictor’s memory
In all experiments we either use a TCN composed of 3 hidden layers with 300 nodes per layer or a TCN with 8 layers and 32 nodes per layer. Moreover we chose N3 = 5 (TCN kernels’ lengths) and relu activation functions (Bai et al., 2018).
Comparison with SOTA methods: We tested our model against other SOTA methods (Table 1) in a comparable experimental setup. In particular, we chose comparable window lengths and architecture sizes (same order of magnitude of the number of parameters) to make the comparison as fair as possible. For the hyper-parameters details of any SOTA method we used we refer the relative cited reference. We point out that while the window length is a critical hyper-parameter for the accuracy of many methods, our architecture is robust w.r.t. choice of window length: thanks to our fading regularization, the user is required only to choose a window length larger than the optimal one and then our automatic complexity selection is guaranteed to find the optimal model complexity given the available data Section 4.1.
Anomaly scores: When computing the F-score we use the predictions of the CUMSUM detector which we collect as a binary vector whose length is the same as the number of available data. Ones are associated to the presence of an anomalous time instants while zeros are associated to normality.
When computing the AUC we need to consider a continuous anomaly score, therefore the zero-one encoded vector from the CUMSUM is not usable. We compute the anomaly scores for each time instant as the estimated likelihood ratios. Since we write the likelihood ratio as pfpp , it is large when data does not come from pp (which we consider the reference distribution).

heading: A.7.1 ABLATION STUDY
text: In Figure 9 we show different metrics based on the predictor’s RMSE (training, test and generalization gap) as a function of the memory of the predictor. We test our fading regularization on a variety of different datasets. In all situations fading regularization helps improving test generalization and preserving the generalization gap (by keeping it constant) as the model complexity increases. All plots show confidence intervals around mean values evaluated on 10 different random seeds.
In Table 4 we extend the results we show in the main paper by adding uncertainties (measured by standard deviations on 10 different random seeds) to the values of train and test RMSE on different ablations of STRIC. Despite the high variability across different datasets STRIC achieves the most consistent results (smaller standard deviations both on training and testing).
Finally, in Table 5 we show the effects on different choices of the predictor’s memory npred and length of the anomaly detectors windows ndet on the detection performance of STRIC. Note both F-score and AUC are highly sensible to the choice of ndet: the best results are achieve for small windows. On the other hand when ndet is large the performance drops. This is due to the type of anomalies present in the Yahoo benchmark: most of the them can be considered to be point anomalies. In fact, as we showed in Appendix A.5, our detector is less sensible to point anomalies when a large window ndet is chosen.
In Table 5 we also report the reconstruction error of the optimal predictor given it’s memory npred. Note small memory in the predictor introduce modelling bias (higher training error) while a large memory does not (thanks to fading regularization). As we observed in Appendix A.5 better predictive models provide the detection module with more discriminative residuals: the downstream detection module achieves better F-scores and AUC.

heading: A.7.2 COMPARISON TCN VS STRIC
text: In this section we show standard non-linear TCN without regularization and proper inductive bias might not generalize on non-stationary time series (e.g. time series with non zero trend component) and TCN architecture. In Figure 10 we
compare the prediciton errors of a standard TCN model against our STRIC on the A3 Yahoo dataset. We train both models using the same optimization hyper-parameters (as described in previous section). Note a plain TCN does not necessarily capture the trend component in the test set.

heading: A.8 STRIC VS SOTA ANOMALY DETECTORS
text: In this section we further expand the discussion on the main differences between STRIC and other SOTA anomaly detectors by commenting results obtained in Table 1 and Table 6. Table 6 is used to highlight the relative performance of STRIC when the peformance are nearly saturated (e.g. Yahoo A2 and A3): we report the relative performance indeces against TRAID for all the models we tested in Table 6. Both for F1 and AUC we report the following for each comparing SOTA method: AUCmethod−AUCSTRIC1−AUCSTRIC · 100 (similarly for F1).
To begin with, STRIC outperforms ‘traditional’ methods (LOF and One-class SVM) which are considered as baselines models for comparing time series anomaly detectors.
Comparison with other Deep Learning based methods: STRIC outperforms most of the SOTA Deep Learning based methods reported in Table 1: TadGAN, TAnoGAN, DeepAnT and DeepAR (the last one is a SOTA time series predictor). Note the relative improvement of STRIC is higher on the Yahoo dataset where statistical models outperforms deep learning based ones. We believe this is due to both fading regularization and the seasonal-trend decomposition performed by STRIC.
Despite the general applicability of GOAD (Bergman & Hoshen, 2020) this method has not been designed to handle time series, but images and tabular data. “Geometric” transformations which have been considered in GOAD and actually have inspired it (rotations, reflections, translations) might not be straightforwardly applied to time series. Nevertheless, while we have not been able to find in the literature any direct and principled extension of this work to the time series domain, we have implemented and compared against (Bergman & Hoshen, 2020) by extending the main design ideas of GOAD to time-series. So that we applied their method on lagged windows extracted from time series (exploiting the same architectures proposed for tabular data case with some minor modifications). We report the results we obtained by running the GOAD’s official code on all our benchmark datasets. Overall, STRIC performs (on average) 70% better than GOAD on the Yahoo dataset and 15% better on the NAB dataset.

heading: A.8.1 DETAILS ON THE NYT EXPERIMENT
text: Figure 4 shows the normalized anomaly score computed by STRIC on the NYT dataset, following the setup described in Appendix A.6.4. Some additional insights can be gained by zooming in around some of the detected change-points. In Figure 12 (left), we see that the anomaly score (blue line) rapidly increases immediately after the 9/11 attack and reaches its peak some days later. Such delay is inherently tied to our choice of time scale, that privileges the detection of prolonged anomalies as opposed to single-day anomalies (which are not meaningful due to the high variability of
the news content). The change-point which occurs the day after the 9/11 attack is reflected by a sudden increase of the relative frequency of article descriptors such as “Terrorism” (orange line). Article descriptors are annotated in the NYT dataset, but they are not given as input to STRIC so that we do not rely on any human annotations. However, they can help interpreting the change-points found by STRIC.
In Figure 12 (right), we can observe that the anomaly score (blue line) is higher in the months around the 2004 U.S. election and immediately after the inauguration day. However, the highest values for the anomaly score occur around the end of 2004, shortly after the Indian Ocean tsunami. Indeed, this is reflected by an abrupt increase of the frequency of descriptors like “Earthquakes” (orange line) and “Tsunami”.
We note this experiment is qualitative and unfortunately we are not aware of any ground truth or metrics (e.g., in Rayana & Akoglu (2015) a similar qualitative result has been reported on the NYT dataset). We therefore tested STRIC against a simple baseline which uses PCA on BERT features and a threshold to detect anomalies. Despite being a simple baseline, this method prooved to be highly applied in practice due to its simplicity (Blázquez-Garcı́a et al., 2020). The PCA + threshold baseline is able to pick up some events (2000 election, 9/11 attack, housing bubble) but is otherwise more noisy than STRIC’s anomaly score. This is likely due to the lack of a modeling of seasonal/periodic components. For instance, the anomaly score of the simple baseline contains many false alarms which are related to normal weekly periodicity that is not easily modeled by the baseline. This does not affect STRIC’s predictions, since normal weekly periodicity is directly modeled and identified as normal behaviour.

Review summaries:
This paper aims to boost the performance of deep neural networks (DNNs) for time series applications by focusing on the characteristics of interpretable forecasting and anomaly detection which are important for real-world time series data. The authors propose an end-to-end trainable DNN architecture which is composed of stacked residual blocks to separate signal components including slow trends, quasi-periodicity, and linear dynamics, followed by a temporal convolutional network (TCN) to model other components. Although previous studies showed that conventional simple linear models often outperformed DNN models on typical time series benchmarks that require robustness and interpretability, this paper shows that the proposed DNN model outperforms state-of-the-art robust statistical methods in some datasets, kind of demonstrating the best of both worlds.