STRIC: Stacked Residuals of Interpretable Components for Time Series Anomaly Detection
Authors: Luca Zancato
Alessandro Achille
Giovanni Paolini
Alessandro Chiuso
Stefano Soatto
Venue: ICLR 2022 Submitted (2022)
Keywords: Anomaly Detection
Time-Series forecasting
Residual Temporal Convolutional Networks

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-of-the-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.

STRIC: STACKED RESIDUALS OF INTERPRETABLE COMPO- NENTS FOR TIME SERIES ANOMALY DETECTION
Anonymous authors Paper under double-blind review
ABSTRACT
We present a residual-style architecture for interpretable forecasting and anomaly detection in mul- tivariate 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-of- the-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.
1 INTRODUCTION
Time series data is being generated in increasing volumes from industrial, medical, commercial and scientific appli- cations. 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 al., 2219)) While recent developments have focused on deep neural networks (DNNs), simple linear models still outperform DNNs in applications that require robustness (Braci & Wagner 2020) and interpretable failure modes (Geiger 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 composed of three modules: interpretable local predictor, a non-linear predictor and novel non- con-aarmettric anomaly detector.
More specifically, STRIC uses a parametric model implemented by sequence residuals blocks with each layer capturing the prediction residual of previous layers. The first layer models trends, the second layer models quasi- periodicity/seasonality at multiple temporal scales, the third layer a general linear predictor, and the last a general non-linear model in the form of 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 The model is trained end-to- end with prediction loss and we automatically select its complexity using an upper bound of the marginal likelihood which, 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 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 statisties (Oreshkin et al. 2019: Cleveland et al. 1990). and incorporates statistics from other time series 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 corre- sponding to the length of observation and enables anomaly detection in the absence of knowledge about the pre- and distibisribu
4. We test our method on standard anomaly detection benchamrks and show merges both the advantages of simple and interpretable linear models and the flexibility of non-linear ones while discounting their major trawbacks: lack of flexibility of linear models and lack of interpretability and overfitting of non-linear ones.
2 RELATED WORK
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 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. error as the discriminant (Brae) & Wagner| (020), 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:
common to use statistics of the prediction
Geiger et al. 2020: Su et al. 2019
Bashar& & Nayak 2020).
Our method follows 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 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 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 introduced (e.g. TCN). In (p.7.dix(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 (T'ang 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.
we
Our method outperforms both classical statistical methods (Braci& & Wagner, 2020) and deep networks (Munir et al. 2019: Geiger et al. 2020) Su et al.] 2019; Bergman & Hoshen] [220 Bashar& & Nayak 2020) on different anomaly detection benchmarks (Laptev & Amizadeh] 2020; Lavin & Ahmad] (015)) (Section (6) Moreover, we show it can be employed to detect anomalous patterns complex data such as text embeddings of newspaper articles (Figure 4).
3 NOTATION
We denote vectors with lower case and matrices with upper case. In particular y is multi-variate time series (iz)) v((1 E R"; we stack observations from time / to L + k and denote the resulting matrix as y = [v((),, v((( + v(! + k - 1)] E R 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 yt and its evaluation at time t ad as (t))) E R. We refer to {u(s), s > t] as the test/future and to {v(s), s < t} as the reference/past intervals. At time t, sub-sequences containing the np past samples up to time t 7p + 1 are given by Y Y(-n,41 (note that we include the
2

---

past
future
Layur impute
IAxrEo-NNAAs- Rcoe!
Amm
7
n
filturs banks
C E'roos 2i sere Mand
PP
Eight
Features
e Io'moss of Far
90 + 1)
of
R 120,
predictors
Frovl't lxeR
Focul! + 1)
Frudt+ DER
Peanc Joor Paol
Future Predictor a
R
Trend
Seasonal
Linear
TCN
Xyrovo R
Xouay "
XaN
(b) Interpretable blocks structure. Time features are extracted independently for each time series (sce Ap- pendix A. for more details).
(a) STRIC predictor architecture.
present data into the past data), while future samples up to time + ny are Y41 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.
4 TEMPORAL RESIDUAL ARCHITECTURE
Our architecture is depicted in Figure Ta Its basic building blocks are causal convolutions (Bai al. 2018). with fixed-sizs size-szze/ii 1-D kernel with input elements from time and earlier. Rather than initializing the convolutional filter randomly, as commonly done in deep learning, we initialize the weights 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 given time series by means of 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 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 model to work wide variety 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 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 TCN model (Sen et al.) 2019). It takes as input the prediction residual of the linear module and outputs a matrix G(Yi-n,+4) E 1p where 7 is the number of output features extracted by the TCN model. The column G(Y: with j = ,,,,.,, Pl of the non-linear features is computed using data up to time ny j (due the internal structure of a TCN network (Bai et al. 2018)). We build a linear predictor on top of G(Y,nn,+1 )) for each single time series independenty: the predictor for the i-th time series is given by: бrcN(t 1); = at G(Yi-n,+1)0 where a, E RR and b, R Since a; combines features (uniformly in time) relevant features across time indices to build the one-step ahead predictor (see Appendix A.
interpret it as feature selector. While b; aggregates
Note that the third layer of the linear module is 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 (c.g., common trend shapes, periodicities) can be extracted. In Appendix A.]
we describe each component of the model detail.
3

---

4.1 AUTOMATIC COMPLEXITY DETERMINATION
Consider the TCN-based future predictor be VITCN(L + 1) := a² G(Y -(Y--,,,,+ = Xrexb where XTeN E RIx8 is the output of the TCN block which depends on the the past window Y, of length np (the memory of the predictor) Ideally, ", 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 large memory window. In this section, we introduce a novel regularized loss inspired by Bayesian arguments which allows us to use an architecture with "large enough" past horizon Rp (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 + 1) Yi_,,+1 ~ NNFF" (Y==n--1)); 1f²) where Fx is the optimal predictor of the future values given the past. Note that this modeling assumption does not restrict 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,11 b. W. Y? +1= I plu(t + k)] b, W, (4+--))
framework and is used only to justify the use of the squared loss learn
To make the notation simpler, we shall denote with Yg == v4 E R the set of future outputs over which the predictor computed and we shall use Ve,w E TR" as the predictor's outputs.
In a Bayesian framework, the optimal set of parameters can be found maximizing the posterior p(b, W (y)) over the model parameters. We model b and W as independent random variables:
p(b, W YY)) x p(Y, b, W) (W)())((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 R. for = 0, 0,...,7, np 1, where K E R++ and X E (0, 1). The corresponding maximum entropy prior (b) Cover& & Thomas 1991) under such constraints is (p)))...))) < |b61 log((Al) where A E TR" is a diagonal matrix with elements AJJ = RA' with j = O,... Ry. Here, À represents how fast the output of the predictor "forgets" the past. Therefore, À regulates the complexity of the predictor: the smaller X, the lower the complexity.
variances: ED2 -
In practice, À has to be estimated from the data. One would be tempted to estimate jointly b, W, X, i. (and possibly v)) by minimizing the negative log of the joint posterior (see Appendix A..2.]. Unfortunately, this leads to a degeneracy since the joint negative log posterior goes to -oo when - O. Indeed, typically the parameters describing the prior (such as x) 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.
The model structure we consider is linear in b and we can therefore stack the predictions of each available time index to get the following linear predictor on the whole future data available: Yo,w = Fwb where F g x is obtained by stacking XYc
Variational upper bound to the marginal likelihood:
for i = t,... + n2 1.
Proposition 4.1. Consider a model on the form: Yo,w = Fwb (linear in b and possibly non-linear in W) and its posterior in Equation T Assume the prior on the parameters 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 with marginalization taken only w.r.t. b:
Us.WA m Y, P/,w -bTA-16 log det(Fw + I2)).
(2)
This (proved in Appendix A.2) provides alternative loss function the negative log posterior which does not suffer of the degeneracy alluded above while optimizing over b, W, À and k.
For simplicity we consider scalar time series, but our approach easily generalizes to multivariate time series (Appendix A.2.2
4

---

Remark: We use batch normalization (lotfe & 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 A.2.3
5 ANOMALY DETECTOR
In this section, present our anomaly detection method based on 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, well as considering global covariates. Such prediction residual is used to test the hypothesis that the time instant is anomalous by comparing its statistics before on temporal windows of length "p and 7y. 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. (212). We circumvent this problem by directly estimating the likelihood ratio with variational characterization of -divergenness (Nguyen et al. 2010) which involves solving 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).
5.1 LIKELIHOOD RATIOS AND CUMSUM
CUMSUM PPgee 1954) is classical Sequential Probability Ratio Test (SPRT) (Basseville & Nikiforov 1993: Liu
al. 2012 of the null hypothesis Ho that the data after the given time e comes from the same distribution as before, against the alternative hypothesis He that the distribution is different. We denote the distribution before c as Pp and the distribution after the anomaly at time c as pj-
If the density functions Pp and py were known (we shall relax this assumption later), the optimal statistic to decide whether a datum v(i) is more likely to come from one or the other is the likelihood ratio s(y(i)). According to the Geyman--Peason lemma, Ho is accepted if the likelihood ratio syy()) is less than threshold chosen by the operator, otherwise He is chosen. In our case, the competing hypotheses are Ho = "no anomaly has happened and H. "an anomaly happened at time e", We denote with pıa, and pH. the p.d.f.s under Ho and H. that: ps,, (Y-)) Pp (Y,K) |YT
and PH. (YT 1Yf-1) = (Y, (Y.
ps (Y
K
Y Therefore the likelihood ratio is:
2¹ = PH. (Y1) PHo (X1)
}ps (Y
PY Po
II Pr((()) Y s=c Py (v(1)) ( 1
To determine the presence of an anomaly, we can compute the cumulative sum S = log S2 of the (log) likelihood ratios, which depends on the time c, and estimate c* using maximum likelihood criterion, corresponding to the detection function h, maxi St The first instant at which we can confidently assess the presence of a change point (a.k.a, stopping time) is: Coop = min(l= : hy 2 7 where T is design parameter that modulates the sensitivity of the detector depending is simply given by the timestamp e at which the maximum of h, = maxisest St is achieved. In Appendix A.33 we provide an alternative derivation that shows that CUMSUM is 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.
the application. The final estimate è of the true change point c* after the detection Estop
5.1.1 LIKELIHOOD RATIO ESTIMATION WITH PEARSON DIVERGENCE
The goal of this section is tackle the problem of estimating the likelihood ratio of two general distributions Pp and Ps given samples. To do so, we leverage a variational approximation of /-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 closed form estimate of the likelihood ratio (see Appendix A.4.
5

---

Proposition 5.1. (Nguyen et al. 2010 Liu et al. 2012) Let Φ == Ps/Pp be the likelihood ratio of the unknown distributions pj and Py- Let UJ Ji Pf: =I,.. nz) and H. = =h:= : hy ~ Pp: i -I...,npp be two sets
containing f and ny samples i..dd from pf and Pp respectively. An empirical estimator & of the likelihood ratio & is given by the solution the following conver optimization problem:
ny
1 TU
δ arg min
((h))²
off)
(4)
2nP
Proposition 5.2. (Liu et al. 2012 Kanamori et al. 2009) Let Φ in Equation ④ 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 Pε and Pp be K; == K K(F,S, (F,) and Kp = K (H,Strr) The optimal regularized empirical likelihood ratio estimator on a new datum is given by:
ne K(e, Str) (K, (opdnptn)) ) K7L
((e)
(5)
Tt
Remark: The estimator in Equation 55) 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.
5.2 SUBSPACE LIKELIHOOD RATIO ESTIMATION AND CUMSUM
In this section, we present our anomaly detector estimator. We test for an anomaly in the data Y! by looking at the prediction residuals Ef, which provide a sufficient representation of Y (see Appendix A.5). We therefore assume we are given the prediction errors E obtained from our time series predictor (the predictor should model the normal behaviour). This guarantees that the sequence E is white in each of its normal subsequences. On the other hand, if the model applied to a data subsequence which contains the abnormal condition, the residuals are correlated.
We estimate the likelihood ratio of p5 and Pp on datum e as (et). &e is obtained by applying Equation (5) on the past window of size Ωpp + ng. At each time instant 1, we compute the necessary kernel matrices as +1, Ef Ei-n,-ny+1 and ( -1
E7 Ei-,--yy
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 lower bound on the actual cumulative sum of likellhood ratios. For 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 = Σ log (e.)).
Remark:
The choice of the windows length (n, and n2 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
6 EXPERIMENTAL RESULTS
In this section, we show STRIC can successfully applied to detect anomalous behaviours on different anomaly detection benchmarks. In particular, we test our novel residual temporal structure, the automatic complexity regu- larization 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 embed- dings of articles from the New York Times. See Appendix 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 used with a proper inductive bias and regularization. In Table [] we compare STRIC against statistical and deep learning based anomaly detection methods. Our experiments follow the experimental setup and evaluation criteria used in Braci & 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 Fl score on Yahoo A3, in Appendix A.7 we show this is mainly

---

Table 1: Comparison with SOTA anomaly defectors: We compare STRIC with other anomaly detection methods
(see Appendix A.8)
8) on the experimental setup and the same evaluation metrics proposed in (Braei & Wagner 20200
Munir et a 2019) The baseline models are: MA, ARIMA, LOF (Shen al., 2020p, LSTM (Braei & Wagner
2020 Munir et al. 2019). Wavenet (Braci & Wagner 2020). Yahoo EGADS (Munir et al. (019)) GOAD (Bergman
& Hoshen 2020), OmniAnomaly (Su et al. 2019). Twitter AD (Munir et al. 2019), TanoGAN (Bashar Nayak
2020), TadGAN (Geiger et al. 2020) DeepAR (Flunkert et al. (017)) and DeepAn'T (Munir et al. 2019). STRIC
outperforms most of the other methods based on statistical models and based on DNNs. See Table for the same table obtained by looking the relative performance w.r.t. STRIC.
Fi-score
Yahoo A1
Yahoo A2
Yahoo A3
Yahoo A4
NAB Tweets
NAB Traffic
ARIMA
0.35
0.83
0.81
0.70
0.57
0.57
LSTM
0.44
0.97
0.72
0.59
Yahoo EGADS
0.47
0.58
0.48
0.29
OmniAnomaly
0.47
0.95
0.80
0.64
0.69
0.70
Models
Twitter AD
0.48
0.26
0.31
TanoGAN
0.41
0,86
0.59
0.63
0.54
0.51
TadGAN
0,40
0,87
0,68
0.60
0.61
0.49
DeepAR
0.27
0.93
0.47
0.45
0.54
0.60
DeepAn'
0.46
0.94
0.87
0.68
STRIC (ours)
0.48
0.98
0.89
0.68
0.71
0.73
AUC
Yahoo A1
Yahoo A2
Yahoo A3
Yahoo A4
NAB Tweets
NAB Traffic
MA
0.868
0.994
0.994
0.986
ARIMA
0.873
0.989
0.990
0.971
LOF
0,904
0,901
0.641
0.640
0.491
0.428
Models
Wavenet
0,824
0.761
0.580
0.592
LSTM
0.812
0.735
0.578
0,589
GOAD
0.893
0.921
0,888
0.866
0.572
0.641
DeepAnT
0.898
0.961
0.928
0.860
0,554
0.637
STRIC (ours)
0.931
0.999
0.999
0.935
0.658
0.685
due to STRIC's predictor. fact, most the time series in Yahoo A3 are characterized by trend components and seasonalities which STRIC's interpretable predictor easily model (see Appendix (A.7.2). In Appendix A.77] we show some ablation studies on the effects of STRIC's prper-parameters on its performance. In particular, we find that STRIC is highly affected by the choice of the length of the windows used estimate the likelihood ratio, while not being much sensitive to the choice of the memory of the predictor (Appendix (A.7.I). Interestingly, STRIC does not achieve the optimal F1/AUC compared to linear models Yahoo A4. The ability of linear models to outperform non- linear 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 non-linear ones while discounting their major drawbacks: lack flexibility of linear models and lack of interpretability and overfitting of non-linear ones (see Appendix A.8 for more in depth discussion).
STRIC interpretable time series decomposition: 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 non-stationary time series (e.g. where the trend component non negligible and typical non linear models overfit, see Appendix (ppendix(A.7.2).
Ablation study: We compare the prediction performance of a general TCN model with 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 (c.g. n, 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 (Braci & Wagner 2020) (we further tested this
7

---

CO2 Daraset
Yahoo Dataset
NAB Datasct
STRCK 2IN0
w
100
130
202
20
450
520
050
720
140
750
807
350
452
500
550
500
700
520
350
350
520
500
5.00
700
220
250
700
160
222
735
750
775
220
360
350
4.0
453
550
032
652
700
Threstampe
Tinerburppe
Teraitappe
Figure 2:
We test STRIC time series intepretability on different datasets (columns). In each panel, we show both training data and test data (see colors). First row: STRIC time series predictor (output of non non-linear module). Second row: Trend components extracted by the interpretable blocks. Third row: Seasonal components extracted by the interpretable blocks.
Table 2: Ablation study on the RMSE of prediciton errors: We compare Test error and Generalization Gap (Gap.) of standard TCN model with our STRIC predictor and some variation of it (using the same training hyper-parameters) Standard deviations are given in Table
TCN
TCN + Linear
TCN + Fading
STRIC pred
Test
Gap.
Test
Gap.
Test
Gap.
Test
Gap.
Yahoo A1
0.92
0.82
0.88
0.78
0.92
0.48
0.62
0.19
Yahoo
0.82
0.71
0.35
0.22
0.50
0.30
0.16
Datasets
Yahoo A.
0.43
0.30
0.22
0.06
0.40
0.25
0.22
0.03
Yahoo A4
0.61
0.46
0.35
0.16
0.55
0.38
0.24
0.01
CO2 Dataset
0.62
0.48
0.45
0.30
0.61
0.43
0.41
0.08
NAB Traffic
1.06
1.03
1.00
0.96
0.93
0.31
0.11
NAB Tweets
1.02
0.84
0.98
0.78
0.83
0.36
0.77
0.07
in Appendix A.7.1p. 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 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 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 closer look on typical time series in CO2 dataset, which contains 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 (so 40% reduction). Moreover, coupling fading regularization with the interpretable module guarantees the best generalization.
infinite memory systems (systems with unitary
Automatic complexity selection: In Figure 3 test the effects of automatic complexity selection (fading mem- ory 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 uniformly smaller than standard TCN (without in- terpretable blocks nor fading regularization). Adding interpretable blocks to standard TCN improves generalization for 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.I for a comparison with other metrics).
8

---

Standard TCN
Generalization Gap
STRIC No Fadling
1.3
2000 election
9/11 allack
2.-
0.5
1.00 occupation
election
0.5
Housing
Soddam Husseln's deach
0.4
0.2
wM
Memory
2001
2002
2003
2004
2005
Figure 3: Automatic complexity selec- tion: Fading memory regularization pre- serves generalization gap as the memory of the predictor np increases on NAB
Figure 4: Anomaly score on the New York Times dataset. Our method finds anomalies a complex time series consisting of the BERT embedding of articles from the New York Times. Peaks in the anomaly score correspond to historical events that sensibly changed the content of the news cycle.
Tweets.
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 n = 30 days, to be able detect change-point anomalies that altered the normal distribution of news articles for 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, see for example Rayana & Akoglu (2015). Additional details and comparison with a baseline model built on PCA are given in Appendix A.8.
are not aware of any ground truth or metrics for this benchmark.
7 DISCUSSION AND CONCLUSIONS
We have shown that our interpretable stacked residual architecture and 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. Munir
Munir et al., 2019 Bashar & Nayak 2020), STRIC exposes 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 interpretable module helps STRIC generalize correctly on non-stationary time series on which standard deep models (such as TCNs) may overfit (Braci & Wagner 2020). Moreover, we showed that our novel fading regularization alone can improve the generalization error of TCNs up to se 40% over standard TCNs, provided the periodic components of the time series are captured by the interpretable module (Section ()) We highlight that the overall computational complexity and memory requirement of our method remains the same standard TCNs, so that our approach can easily scale to large scale time series datasets.
important in building a successful model (see Table 2 In particular, the
An anomaly is time instant: at that time instant, either we receive isolated observation that is inconsistent with normal operation (outlier measurement), or discrete change occurs in the mechanism that generates the data (change- point) that persists beyond that time instant (Geiger et al. 2020) Basseville & Nikiforov 1993). Our method treats these two phenomena a unified manner, without the need to differentiate between outliers and setpoint changes, with specialized detectors for each. Once the predictor is built, 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 the prediction error distribution (nominal and faulty).
Making the non con-aarmetric 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 T depending on the desired operating point in the tradeoff between missed detection and false alarms would further robustness of our method.
enhance the

---

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 implemen- tation of our novel method in Appendix A. while in Appendix we describe both the model structure and training sypee-parammters we used in the experimental section. We do not include the model structures and hyper-parametert of SOTA methods we used as baselines and refer to related literature (referenced in our work) for the implementation details. To further foster reproducibility shall make our code available.
REFERENCES
Ratnadip Adhikari and R. K. Agrawal. An introductory study on time series modeling and forecasting. CoRR, abs/1302.6613 2013. URL http: org/ abs 1022. 6613.
Shaojie Bai, J Zico Kolter, and Vladlen Koltun. An empirical evaluation of generic convolutional and recurrent networks for sequence modeling. arXiv preprint arXiv: 2018.
Md Abul Bashar and Richi Nayak. Tanogan: Time series anomaly detection with generative adversarial networks. In 2020 IEEE Symposium Series on Computational Intelligence (SSCI), pp. 1778-1785. IEEE, 2020.
Michèle Basseville and Igor V. Nikiforov. Detection of Abrupt Changes: Theory and Application. Prentice-Hall, Inc., USA, 1993. ISBN 0131267809.
detection Liron Bergman and Yedid Hoshen. Classification- anomaly detection for general data. In International Confer- ence on Learning Representations, 2020. URL httll/tttee///eeeerrreeeeeeeeees
Ane Blázquez-Gareía, Angel Conde, Usue Mori, and Jose A Lozano. A review on outlier/anomall detection in time series data. arXiv preprint arXiv:2002.04236. 2020.
Mohammad Braei and Sebastian Wagner. Anomaly detection in univariate simes::::: survey on the state-of-the-an CoRR, abs/2004.00433, 2020. URL
Robert B. Cleveland, William Cleveland, Jean McRae, and Irma Terpenning. Stl: A seasonal-tres -trend decomposition procedure based on loess (with discussion). Journal of Official Statistics, 6:3-73, 1990.
T. M. Cover and J. A. Thomas. Elements of Information Theory. Series in Telecommunications and Signal Processing. Wiley, 1991.
Jacob Devlin, Ming-Wei Wei Chang, Kenton Lee, and Kristina Toutanova. BERT: pre-training of deep bidirectional trans- formers for language understanding. In Jill Burstein, Christy Doran, and Thamar Solorio (eds.), Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, NAACL-HIT 2019, Minneapolis, MN, USA, June 2-7, 2019, Volume (Long and Short Papers), pp. 4171-41866 Association for Computational Linguisties, 2019. doi: 10.18653//1//19911
URL
https: /ttp///ciiiiii
14423
Amir- Amir-massoud Farahmand, Sepideh Pourazarm, and Daniel Nikovski. Random projection filter bank for time series data. In Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vish- wanathan, and R. Garnett (eds.). Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017. URL http:::/pproceedinge www.nnnur..... htt.:ppp...cppppee////////// Paper.
Valentin Flunkert, David Salinas, and Jan Gasthaus. Deepar: Probabilistic forecasting with autoregressive recurrent networks. CoRR, abs/1 ftss//704.41100 2017. httpp/tppp//aa....////oe///
C4110
Alexander Geiger, Dongyu Liu, Sarah Alnegheimish, Alfredo Cuesta-Infante, and Kalyan Veeramachaneni. Tadgan: Time series anomaly detection using generative adversarial networks. arXiv preprint arXivv2009.07766,, 2020.
Vincent Le Guen, Yuan Yin, Jérémic Dona, Ibrahim Ayed, Emmanuel de Bézenac, Nicolas Thome, and Patrick Gallinari. Augmenting physical models with deep networks for complex dynamics forecasting. arXiv preprint arXiv:2010.04456, 2020.
10

---

Sergey loffe and Christian Szegedy. Batch normalization: Accelerating deep network training by reducing internal 031
covariate shift. CoRR, abs/1502.03167, 2015. /tt.
Takafumi Kanamori, Shohei Hido, and Masashi Sugiyama. A least-squares approach direct importance estimation. Journal of Machine Learning Research, 10(48):1391-1445, 2009. URL http:///alll..orr/paarr/////
kanamori 09a l
Nikolay Laptev and Saced Amizadeh. Yahoo! webscope dataset yultttl-bbbcc-----e----rooomoooa CoRR, 2020. htttpps::/wwbsootttooo sandbox, yahoo. com/ htt.aaaaaooyyp.paaaaaaaaaeeeeaaa
Alexander Lavin and Subutai Ahmad. Evaluating real-time anomaly detection algorithms the numenta anomaly benchmark. CoRR, abs/1510.03336, 2015. URL.sstp://arravv..rrr//bbbb 03336
Song Liu, Makoto Yamada, Nigel Collier, and Masashi Sugiyama. Change-point detection time-series data by relative density-ratio estimation. In Structural, Syntactic, and Statistical Pattern Recognition, pp. 363-372, Berlin, Heidelberg, 2012. Springer Berlin Heidelberg. ISBN 978-3-642-34166-
Mohsin Munir, Shoaib Ahmed Siddiqui, Andreas Dengel, and Sheraz Ahmed. Deepant: A deep learning approach for unsupervised anomaly detection in time series. IEEE Access, 7:1991-2005, 2019, doi: 10.1109/ACCESS.2018 2886457.
XuanLong Nguyen, Martin J. Wainwright, and Michael I. Jordan. Estimating divergence functionals and the likelihood ratio by convex risk minimization. IEEE Transactions on Information Theory, 56(11))8877778-811, Nov 2010. ISSN 1557-9654. doi: 10. 1091100 2010..0080.0.
//dx doo.ooorgggooo ob..orrg/10.1110///11.
110.1110991111 2010 206887000000000
Boris N. Oreshkin, Dmitri Carpov, Nicolas Chapados, and Yoshua Bengio. N-BEATS: neural basis expansion analysis for interpretable time series forecasting, CoRR, abs/1905. 2013,,099,, URL http::ppp//xx...orraaaa/////////
10437
Ewan $ Page. Continuous inspection schemes. Biometrika, 41(1/2):100-115 1954.
CE. Rasmussen and CKI. Williams. Gaussian Processes for Machine Learning. Adaptive Computation and Machine Learning. MIT Press, Cambridge, MA, USA, January 2006.
Morten Ravn and Harald Uhlig. On adjusting the Mndrick-prescott filter for the frequency of observations. The Review of Economics and Statistics, 14:371-375, 02 2002. doi: 10.1162/0034446532311111 1604.
Shebuti Rayana and Leman Akoglu. Less is more: Building selective anomaly ensembles with application to event detection in temporal graphs. In Suresh Venkatasubramanian and Jieping Ye (eds.), Proceedings of the 2015 SIAM International Conference on Data Mining. Vancouver, BC, Canada, April 30 May 2. 2015, pp. 622 630. SIAM, 2015. doi: 10.1137/1 97811119744000.700 CRL.hhttps:/ddo: http://oooorggg////////// 0. 1137/1. 9781 61 974010 7.70
Evan Sandhaus. The new york times annotated corpus 1dc2008119. web download. Linguistic Data Consortium, Philadelphia, 6( 6(22)222 2008
Rajat Sen, Hsiang-Fu Yu. and Inderjit S Dhillon. Think globally. act locally: A deep neural network approach to high-dimensional time series forecasting. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Bue, E. Fox, and R. Garnett (eds.), Advances in Neural Information Processing Systems, vol- ume 32. Curran Associates, Inc., 2019. URLhttps: httptt:::/ppoooeediig http///pepoo/.ss////////p// htt.ooeeedingn.nuuiiiis....e.. httpsscccaaaeee/1///ii.. (ee44ceeeooo////eeeeeeeeeee
Paper.pdf
Lifeng Shen, Zhuocong Li, and James Kwok. Timeseries anomaly detection using temporal hi- erarchical one-class network. In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and
H. Lin (eds.). Advances in Neural Information Processing Systems. volume 33, pp. 13016-13026. Curran Associates, Inc., 2020. URL hhttps://proceedingg neurips cccppaper/222//ELL. 97e401a02022011662455 1w1.4958888220047--PPPe......
Ya Su, Youjian Zhao, Chenhao Niu, Rong Liu, Wei Sun, and Dan Pei. Robust anomaly detection for multivariate time series through stochastic recurrent neural network. In Proceedings of the 25th ACM SIGKDD International
Conference on Knowledge Discovery & Data Mining, KDD '19, pp. Association for Computing Machinery. ISBN 9781450362016. doi: 10.1145532925000333666 URL https: //ddl.oogg/0..115//2222222 3330672
New York, NY, USA. 2019.
11

---

Michael E Tippingg Sparse bayesian learning and the relevance vector machine. Journal of machine learning research. 1(1un):211-244, 2001.
Michael Tsang, Hanpeng Liu, Sanjay Purushotham, Pavankumar Murali, and Yan Liu. Neural interaction transparency (nit): Disentangling learned interactions for improved interpretability. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (eds.), Advances in Neural Information Processing Systems, vol- neurips. httpr/ppps.cc/ppapee///////////
ume 31. Curran Associates, Inc., 2018.
Vladimir N. Vapnik. Statistical Learning Theory. Wiley-Interscience, 1998.
Thomas Wolf, Lysandre Debut, Victor Sanh, Julien Chaumond, Clement Delangue, Anthony Moi, Pierric Cistac, Tim Rault, Rémi Louf, Morgan Funtowicz, Joe Davison, Sam Shleifer. Patrick von Platen, Clara Ma. Yacine Jernite, Julien Plu, Canwen Xu, Teven Le Scao, Sylvain Gugger, Mariama Drame, Quentin Lhoest, and Alexander M. Rush, Transformers: State-of-the-an natural language processing. In Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing: System Demonstrations, pp. 38-45, Online, October 2020. Association for Computational Linguisties. URL htRLhhtps::///t www.:wwwwwww........... http://www.......ooooooooooo////
Emmanuel Yashchin. Performance of cusum control schemes for serially correlated observaatooosssssss Technometrics, 35 (1):37-52, 1993. ISSN 00401706. URL http:/httt//// jstor. /rt.tt/o.org/ttaaaee//22222
Luca Zancato and Alessandro Chiuso. A novel deep neural network architecture for non-linear system identification. Peperstins,,,,inne, 54(7):186-191, 2021. ISSN 2405-8963. doi: https///cb..oog//00000...fftoo.. 08..366.
URL https::/wwy www. sciencedirect com http///ww....cccccee/////////aaa IFAC Symposium on System Identification SYSID 2021.
19th
A APPENDIX
A.I IMPLEMENTATION
Given two scalar time series x and ip, we denote their time convolution g as: g(t) == (xxx)) : Ei (())x(t-i). We say that g is the causal convolution of x and yp if pp(t) for f < 0, so that the output g(t) does not depend on future values of x (w.rtt current time index t). In the following, we shall give a particular interpretation to the signals and yo: x will be the input signal filter parametrized by an impulse response yp (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 finite window. The truncation is indeed 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(l) == (4 (1) and mean that the impulse response of the causal filter is truncated on a window of a given length.
A.1.1 ARCHITECTURE
Let Mi-np+1 E Raxnn be the input data to our architecture at time step (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 specific layer using either the subscripts (TREND, SEAS, LIN, TCN] or (0, 1, 3.
We shall denote the input of each block as Xg E Rnnx, and the output as X₂ E RuxA, for = 0, 1, The residual architecture we propose is defined by the following: Xo = Yi-nyt1 and X = Xx-1 - X for k = 1, 3. At each layer we extract l temporal features from the input X₁. We denote the temporal features extracted from the input of the k-th block as: G₁= == G& (Xx) R4xx,, The i-th column of the feature matrix G₁ is a feature vector (of size 4+) extracted from the input X₁ up to time t ny i. To do so, we use causal convolutions of the input signal Xς with a
set of filter banks (Bai et al. 2018).
12

---

A.1.2 INTERPRETABLE RESIDUAL TCN ON SCALAR TIME SERIES
Modeling Interpretable blocks: In this section, 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 FTREND- Then STREND is defined as a multidimensional time series (of dimension TREND (o) obtained by stacking lo time series given by the convolution of x with to causal linear filters: (TREND; I for i = o,... to - 1. In other words, STREND == (PTRENDD zz (TREND22 *I,, .... PTRENDI, * +x]]". We denote the of linear filters (PTREND; for O,..,, l- as KTREND and parametrize each filter in CTREN with its truncated impulse response (i.e. kernel) of length kg- == FIREND
We interpret each time series in FIREND as an approximation of the trend component of computed with the i-th filter, We design each (IREND; 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 STREND with the linear map @TREND- Moreover, we predict the future trend of the input signal (on the next time-stamp) with the linear map STREND.
We construct the blocks that extract seasonality and linear part in a similar way.
The input of each layer is given by a window of measurements of length ni. 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: Go GTREND E o G1 GSEAS R X and G2 GLIN E R
Implementing Interpretable blocks:
The output of each layer X2 is estimate of the trend, seasonal or stationary linear component of the input signal the past interval of length np, so that we have X Ixln, (same dimension as the input Xx)). On the other hand, the linear predictor dR computed at each layer is a scalar. Intuitively, X, and ÿu. 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: Xx+ := a.Gr, GR and ÿk == X&be for k O, 1, 2 where a; E R and 6ş- E " Note ag combines features (uniformly in time) so that we can interpret it as a feature selector while by. 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 by depending on the time index (similarly to the fading memory regularization). We experimented both the choice to make by canonical "vectors" and dense vectors. We found that choosing by 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 = xi-npt1 - X, - X1 - X2. which is to be considered as a signal whose linearly predictable component has been removed. The TCN extracts set of la non-linear features Gx(Xa) Rgaa, 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 - npp + j (due to the internal structure of a TCN network (Bai et al. 2018)). The linear predictor on top of G3 is ŪICN := am G3b3, where a3 E R's and by E RT>
Finally, the output of our time model is given by:
g( 1)
Uk
Xaba
of GX
km
A.1.3 INTERPRETABLE RESIDUAL TCN ON MULTI- DIMENSIONAL TIME SERIES
We extend our architecture to multi-dimennional time series according to the following principles: preserve inter- pretability (first module) and exploit global information to make local predictions (second module).
In this section, the input data to our model Y Y E Rnxx, (a window of length ny 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
13

---

the same filter banks: KIREND. KSEAS and KLIN For ease of notation we shall now focus only on the trend layer. Any other layer obtained by substituting 'TREND" with the proper subscript ('SEAS' or 'LIN").
We denote by &TREND, i E R the set of time features extracted by the trend filter bank KTREND from the i-th time series. Each feature matrix then combined as done in the scalar setting using linear maps, which we now index by the time series index i: @TREND; and ÓTRENDE. 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 (TRENDs) or might need to look at values further in the past (retrieved using DTRENDS). We stack the combination vectors dTREND: and DIREND: into the following matrices: AIREND = [aTRENDI TTREND2, ..., @TREND"] 1T E R=x and BIREND = [IRENDI, DIREND2; ...., bYRENDn]" R"**,.
The second (non-linear) module aggregates global statistics from different time series (Sen et al. 2 using TCN model. It takes input the prediction residual of the linear module and outputs a matrix GTREND (1,,,++1) E R where 13 is the number of output features that are extracted by the TCN model (which is design parameter). The j-th column of the non-linear features (Y
Non-linear module:
is computed using data up to time f - -- j. where 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 linear predictor on top of GTREND (x,,,++1) for each single time series independenty:
the predictor for the i-th time series is given by: ÚTCN(!+ 1); == af (YR +1)b; where a; E R and b, E "
We stack the combination vectors @rcN; and byeN; into the following matrices: ATCN [#TCN1, #TCN2 ..., @rCNa] ] E guxls and BTCN = [bTCNI, bTCN2, ..., brCNn]" 1 E _R*x,
Finally, the outputs of the predictor on the i-th time series are given by:
0(l + 1):
Σ
ax. T Grahas + OTTN; o GTexbreN:-
AE STREND.S
A.1.4 BLOCK STRUCTURE AND INITIALIZATION
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 's kernel parameters K& E R4*X,, where 44 and N₂ are the block's 's number of filters and block's 's kernel size, respectively (without loss of generality, we assume all filters have the same dimension). Each filter (causal 1 ID-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 multidimen- sional time series, we apply the filter banks 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 lo causal Hodrick Prescott (HP) filters (Ravn & Uhlig 2002
of kernel size No. HP filters are widely used to extract trend components of signals (Ravn & Uhlig, 2002). In general 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 hyper-parameter Aup which defines the regularity of the filtered signal (the higher AHP, the smoother the output signal). We initialize each filter with λнр chosen uniformly in log-scale between 10 and 10". 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 λEр and search for the optimal λμP- Instead, we optimize over the impulse response parameters of the kernel which we do not assume live in any manifold (c.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 l 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.
14

---

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 22 increases, this random filter bank is guaranteed to be 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 purt) 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 a-priori standardization (e.g., trend or seasonality removal).
A.2 AUTOMATIC COMPLEXITY DETERMINATION (FADING MEMORY REGULARIZATION)
In this section, we shall introduce a regularization scheme called fading regularization, constrain TCN representa- tional capabilities.
The output of the TCN model is G(Y G())-n,+1) E R where 1a is the number of output features extracted by the TCN model. The predictor build from TCN features is given by: orCN? T GieN(Y? +-n DreNNi where the predictor bIGNs E R" takes as input a linear combination of the TCN features (weighted by oreN)). The j-th column of the non-linear features G(Y +(Yi-n,+1)) ) is computed using data up time I np + (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 = np. the statistic is the one computed on the closest window of time w.n.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, property which is commonly known as fading memory. Currently, this property not built in the predictor ÓTCNE, 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 specific regularization which we now describe.
A.2.1 FADING MEMORY IN SCALAR TIME SERIES
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 no measures can be written as: Ircx(t + 1)- a GTCN((Ti (1-b Xgb. We shall assume that innovations (optimal prediction errors) are Gaussian, so that y(t + 1) Yi (V +1),,,"))); 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 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 6 in the weight vector W (linear filters parameters KTREND, KSEAS, KLIN, linear module recombination weights ALIN, BIREND, BSEAS, BLIN and TCN kernel parameters and recombination coefficients ArcN etc.).
ATREND, A ASEAS, We write the conditional likelihood of the future given the past data of our parametric model as:
ns
p(Yttr 1b, D,W,Y'S yd +/
II plu(t + k) b, W, yt+k-1) A=
b+k-np
(6)
To make the notation simpler, we shall denote by Y= : yrtt" ER" the set of future outputs over which the predictor
is computed and we shall use Y..w E R as the predictor's outputs. Moreover, we shall drop the dependency on the
conditioning past Y Y +1 (which is present in any conditional distribution). Equation (6) becomes: p(Y, 6, W) =
Ihe p(y(t + k) 6. 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| (y)) x p(Y, 6, W))
(7)
15

---

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 too high and therefore it should only depend on the most relevant past.
The prior does not induce hard constraints. 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:
Remark:
log p(b, W yj) = - log pY,; b, W) - p(b) log p( (W). In particular, the first term pY, 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 E R", have zero mean and exponentially decaying variances:
Eby- = 0 and Eb2 カ-ーーー1 = RA for j = 0,...,,,, Rp
(8)
where A E R+ and x E (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 199T). 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 as:
log (b) x [|6]]2 - log |A|
where E R is diagonal matrix whose elements are A, = RX for j = O..., Up -
The parameter À represents how fast the predictor's output 'forgets' the past: the smaller X, 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, X, (and possibly v)) by minimizing the negative log of the joint posterior:
1
arg min Yy - Po,w + log(n?") log(px.. - llo)(P(W)). 6, b,W. 1/2
(10)
Unfortunately, this leads to a degeneracy since the joint negative log posterior goes to - oo when X -> 0.
Bayesian learning formulation for fading memory regularization:
The parameters describing the prior (such as x) 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 the whole future data: Yo,w F where Fw Ra x and its rows are given by (Y; 1) for
f, t,..., + ny mg-1. 1.
We are now ready to find an upper bound to the marginal likelihood associated to the posterior given by Equation 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:
1 IY,
arg min b
Fbll 2 + +b
(11)
1
with Σ = FAFTAFFA AT TT
Equation guarantees that
Y Fbll
log [og/E) Y,T "y + log E|
16

---

where the right hand side is (proportional to) the negative marginal likelihood with marginalization taken only w.r.t b. Therefore, for fixed W.
1 1
v,
%.wl
+ A log Fl + T
is an upper bound of the marginal likelihood with marginalization over 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 1 y Yo.w + |1b1 + log AFw 1 l(W) (,
(12)
17
Remark: log p(W) defines the regularization applied on the remaining parameters of our architecture. In particular, we induce sparsity by applying L regularization on ATREND. ASEAS, ALIN and AreN. Also, we constrain filters parameters to stay close to initialization by applying L regularization KIREND, KSEAS and KLIN.
A.2.2 FADING MEMORY IN MULTIVARIATE TIME SERIES
In the case of multivariate time series, fading regularization be applied either with a single fading coefficient X 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 X.
A.2.3 FEATURES NORMALIZATION
We avoid the non-identifiability of the product Fy by exploiting batch normalization: we impose that different features have comparable means and scales across time indices i = O.... O.....?, " 1. Non-identifiability occurs due to the product b, features have different scales across time indices (i.e., columns the matrix Fw) the benefit of fading regularization might reduced since it can happen that features associated with small b, have large scale so that the overall contribution of the past does not fade. Hence we use batch normalization to normalize time features. Then affine transformation (with parameters to be optimized) to jointly re-seale all the output blocks before the linear combination with
A.3 ALTERNATIVE CUMSUM DERIVATION AND INTERPRETATION
In this section, we describe an equivalent formulation of the CUMSUM algorithm we derived in the main paper. (y(3)
Before change point, by construction we are under the distribution of the past. Therefore, log p>(8) S O Vy, which in turn means that the cumulative sum S will decrease as t increases (negative drift). After the change, the situation is opposite and the cumulative sum starts to show positive drift, since Py- 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 need know is the difference between the value of the cumulative sum of log-likelihood ratios and its minimum value.
sampling y(1) from the future distribution
The CUMSUM algorithm can be expressed using the following equations: Vr St Mr, where my == minj S). The stopping time is defined as: (stop = min(l 0 2 - = min(t S 2 my + T). With the last equation, 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 threshold m, TT. Note that the adaptive threshold keeps complete memory of the past ratios. The two formulations are equivalent because S m: f4.
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 probability distribution IP and random variable f measurable w.r.t. P. we use f fdl denote the expectation of under P. Given samples r(1),..,, I(n) from P. the
17

---

empirical distribution Pn is given by Pn n/2 Ei= 82(1)- We use f fal,, as a convenient shorthand for the empirical expectation
Consider two probability distributions and Q. with P absolutely continuous w.r.t. Q. Assume moreover that both distributions absolutely continuous with respect to the Lebesgue measure fa, with densities pu and qo. respectively, on some compact domain X' C C...
Variational approximation of the f-divergence: The f --divergence between IP and Q is defined as Nguyen et al. 2010)
Dj(P, Q = / pof
40 Po
du
(13)
where / RR- + R is a convex and lower semi-continuous function. Different choices of f result in variety of divergences that play important roles in various fields (Nguyen et al.) 2010). Equation (13) is usually replaced by the variational lower bound:
/ [odQ - (6))FF)
Dy(P.Q) 2 sup deo
(14)
and equality holds iff the subdifferential 0f(m) function of ʃ.
()) contains an element of Φ. Here f is defined as the convex dual
In the following, we are interested in divergences whose conjugate dual function is smooth (which in turn defines commonly used divergence measures such as KI. 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 O ))) = o for some & E 4.
Remark: The dimensooaii optimization problem in Equation sup Bgo-Epff(0 (o).
14 can be written as Dy Q)
In practice, one can have an estimator of any f-divergence restricted to a functional class & by solving Equation [4 (Nguyen al., 2010). Moreover, when P and Q are not known one can approximate them using their empirical counterparts: a and Q.. Then an empirical estimate of the f-divergence is: D Q sup Eq. - (
An estimate of the likelihood ratio can be directly obtained from the varia- tional approximation of f-divergences. The key observation is the following: equality on Equation (14) achieved iff This tells us that the optimal solution to the variational approximation provides us with an estimator of the composite function O du of the likelihood ratio de As long can invert of, likelihood ratio.
Approximation of the likelihood ratio:
Φ = 0 (o)
Pa
uniquely determine the
In the following, we shall get an empirical estimator of the likelihood ratio in two separate steps. We first solve the following:
On
arg max Es, Φ (()),
(15)
which returns an estimator of Of( ()), ), not the ratio itself. And then we apply the inverse of 0 to Φn. We therefore have a family of estimation methods for the likelihood function by simply ranging over choices of f.
Remark: If is not differentiable, then we cannot invert 0f 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.
1.4.1 LIKEL RATIO ESTIMATION WITH PEARSON DIVERGENCE
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).
18

---

The Pearson, x2, divergence is defined by the following choice: [(1)
(t-1)²
The associated convex dual
function is:
(u
2
v²
ssu
+ u.
vER
2
Therefore the Pearson divergence is characterized by the following:
PE(PIQ) = /
Po
40 Po
- 1) du 2 sup Ro0
2
Rpp² Epp.
(16)
ΦE
Solving the lower bound for the optimal @ provides us an estimator of 0 9 the Pearson divergence, we can apply change of variables which preserves convexity of the variational optimization problem Equation (T6) and provides more straightforward interpretation. Let the new variable be : + with E 2, which in this case is nothing but the inverse function of aʃ. We get
2u 1. For the special case of Po
Po
1
sup EGO EpO² EpO
sup Eq² se2
Epp² 2
(17)
2
2
It is now trivial to see that 2 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 Epp - Rog SE# 2
(18)
to build our 'direct' estimator of the likelihood ratio.
Let the samples from and Q be, respectively, zp(i) with 1, L,... By and zi(() with i L,..., ng. We define the empirical estimator of the likelihood ratio o.:
T2e Σ (())
6 arg min
(1)(
(19)
2mp
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 & E ① where & is a RKHS induced by the kernl k.
We exploit the representer theorem to write general function within as:
(xx) Σ k(x, Fer (i))
where we use 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 â.). For now, we do not specify which data should be used as centers (we can use either data from P,, or from Qn or from both or simply use user specified locations).
Let us define the following kernel matrices: K3 := K(Xp. Xtr) Rape Kg = K(Xgr Xtx) E Rer where Xp (I,(i)), X, ;= (i)} and X,, == (xxe(())
We therefore have:
nn
1 ng
Gn = arg minl
( (x)(1))²
(lxx(1))
DEo
2mp
Blg
ny
arg min
k(ap(1), Fer (i))ag)²
ΣΣ
(i)
a.020
2np
TT ==
1=1
arg min
T
KT
K p
1 1 Kya
a,a20
2mp
nf
19

---

Remark: We impose the recombination coefficients to be non negative since the likelihood ratio is a non negative quantity. The resulting optimization problem 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 quadratic regularization term as done in (guuynn et al. 2010) which lead us to the following regularized optimization problem:
1 T K2 Kpa 1 - T K
+ 2 |all 2
arg min
Or 2mp
nt
whose solution is trivially given by:
à pp n
KT Kp n nytTne))
KF1 "
" H K
(20)
nf
The estimator of the likelihood ratio for an arbitrary location is given by the following:
Re (n(()) K (x, Xtr)a = n2 K (x,, Xtr) K Kp + upplner) ) KJ1
(21)
(x
nt
Remark: In the following we shall exploit RBF kernels which are defined by the length scales o.
A.5 SUBSPACE LIKELIHOOD RATIO ESTIMATION AND CUMSUM
In this section we describe our subspace likelihood ratio estimator and its relation to the CUMSUM algorithm. The (=((())) for each time t. We denote Pp as the CUMSUM algorithm requires to compute the likelihood ratio pn()(()) Y normal density and py the abnormal one (after the anomaly has occurred).
We shall proceed express the conditional probability p(y(() | Y) using our predictor. In particular it is always possible to express the optimal (unknown) one-step ahead predictor as:
1 F* FY; = Fl((()) Y [=1
(22)
which is 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,+ + e(t))
(23)
where e(t) = y(t) kY+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 p(e(t) Y Y(-K+1). Where e(t) is the optimal prediction error for each time / and is therefore indipendent on each time I.
Remark: 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! can be written as:
P(Y.) = II p(u(1) [Y--)) = (((())
normal conditions
(24)
t=c
t=c
p(Y.) II p(v(() (Y= = II [(((()) E
abnormal conditions
(25)
t=c
These two conditions in turn influence the log likelihood ratio test as follows: under Ho I Pole(t) while
p)(e())
under He - I' Pr(e((]E:-1) The main issue here is the numerator under Ha: the distribution of residuals
pa(e(()
changes at each time-stamp (it a conditional distribution) and py (e(i) F(-1) is difficult to approximate (it requires the model of the fault). In the following we show that replacing py (e(i) E with (e(((i) allows us to compute
20

---

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 Py(e(i) - as independent random variables (e()()) for = I,..., I allows us to compute lower bound on the log likelihood log St (i.e. the cumulative sum). We denote the cumulative sum of the log likelihood ratio using independent variables as ()(e((
log 2 Ei=c log pp(e())
Proposition A.2. Assume change happens at time so that He is true and the following log likelihood ratio holds true: log 0 = Eize log pp(e(i)
Pr(e(s) Then it holds log 22 2 loggi"
Proof. By simple algebra we can write:
ps(e(i) | 1 Ps(e(i) E pr Pp(e(())
Vi
Pf(e(i)) pp(e())
Now recall the cumulative sum of the log-likelihood ratios taken under the current data generating mechanism py EE) provides an estimate of the expected value of the log-likelihood ratio. Due to the correlated nature of data E the samples are drawn from a multidimensional distribution of dimension sample from this distribution an entire trajectory from c to t),
We now take the expectation of previous formula w.r.t. the 'true' distribution (E{)::
Py (e(i) Pp(e(i))
E
E (E) 2 = (E log II
i=c
-E) log II pf(e(() 1 ps(e(i))
ps (e(())
+E ())) log
Pr (p(()))
(M) II)(e(t)) KL ( (Ir p(
2 KL. (I())(())) (I)()())
i=c where we used the fact the mutual information is always non negative.
i=c
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 [.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--"----
and E5 == Et Ei- So that given current time we look back at a window of length np + ny. The underlying assumption is that under H, normal data are present in E,, and abnormal ones in Ey. We estimate the likelihood ratio (y(e(t)) each time by assuming both E, and Es data are independent (see Proposition A.2) and cumulate their log as f increases.
Py(e(t))
How do Np and n affect our detector? The choice of the windows length (n, and ny) is fundamental and highly influences the likelihood estimator. Using small windows makes the detector highly sensible to both point and sequen- tial outliers, while larger windows better suited estimate only sequential outliers. We now assume np nn and study how small and large values affect the behaviour of our detector in simple working conditions.
In Figure Jnd Figure 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 7) increase of the cumulative sum is candidate abnormal point.
In Figure Figure 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 n are usually not able to capture point anomalies Figure 7 while using small window sizes allow to detect both (at the expenses of more sensitive detector) Figure g
21

---

3.5
3.0
5
2.5
Values
2.0
Values
cumsum
1.5
2
point outlier
1.0
1
cumsum
0.0
change point
40
50
60
70
80
90
100
40
50
60
70
80
90
100
Time
Time
Cumulative sum (blue) ob- tained with our method synthetic example. We use the cumulative sum of estimated likelihood ratios on data in which a change point is present at 60. We use np nz 20 and kernel length scale=0.2
Figure 5: Change point:
Figure 6: Point anomaly: Cumulative sum (blue) ob- tained with our method in a synthetic example. We use the cumulative sum of estimated likelihood ratios data in which a point outlier is present at & 60. We use ny n2 2 and kernel length scale=2.
cumsum
Data
My Mon.
anamalous points
Values
Values
cumsuu
Data
anomalaus points
50
75
100
125
150
175
200
50
75
100
125
150
175
200
Time
Figure 7: Large np and ny: Cumulative sum (blue) obtained with method in synthetic example. We use the cumulative sum of estimated likelihood ratios on data which contain both change points (t 60 and 200) and point outliers (t 80 and = 160). We use mp n 20 and kernel length scale=1.
Figure Small np and ny: Cumulative sum (blue) obtained with our method in synthetic example. We use the cumulative sum of estimated likelihood ratios on data which contain both change points t = 200) and point outliers (t = 80 and t = 160). We use np = ny = 3 and kernel length scale=5.
60 and
A.6 DATASETS
A.6.1 YAHOO DATASET
Yahoo Webscope dataset (Laptev & Amizadeh, 2020) is publicly available dataset containing 367 real and synthetic time series with point amomalies, contextual anomalies and change points. Each time series contains 1420-1689 time stamps. This dataset is further divided into 4 sub-benchmarks A1 Benchmark, A2 Benchmark, A3 Benchmark and A4 Benchmark. A Benchmaark is based on the real production traffic to some of the Yahoo! properties. The other 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).
22

---

A.6.2 NAB DATASET
NAB (Numenta Anomaly Benchmark) (Lavin & Ahmad 2013) is 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 ct al. 2020) we choose subset of NAB benchmark, in particular we focus on the NAB Traffic and NAB Tweets benchmarks.
A.6.3 CO2 DATASET
We test the prediction and interpretability capabilities of our model on the CO2 dataset from kaggle The main goal here is to predict both trend and periodicity of CO2 emission rates different years. Note this is not Anomaly detection task.
Table 3: Datasets summaries. We report some properties the datasets used (see Geiger et al. 2020 for mode details).
Yahoo
NAB
Kaggle
A4
Traffic
Tweets
CO2
signals
67
100
100
100
7
10
9
Properties
anooaliiss
178
200
835
point
68
33
935
833
0
sequential
110
167
2
14
33
anomalous
1669
466
943
837
1560
15651
points (% tot)
1.8%
0.32%
0.56%
0.5%
9.96%
9.87%
data points
94866
142100
168000
168000
15662
158511
4323
A.6.4 NYT DATASET
The New York Times Annotated Corpus Sandhaus, 2008) contains over 1.8 million articles written and published by the New York Times between January 1987 and June 19, 2007. We pre- pre-processed the lead paragraph of each article with pre-trained BERT model (Devlin al., 2019) from the HuggingFace Transformers library (Wolf al. 2020 and extracted the 768-dimensional hidden state of the [CLS) token (which serves as an article-level embedding). For took the mean of the embeddings of all articles from that
each day between January 1, 2000 and June 19, 2007, day. Finally, we computed a PCA and kept the first 200 principal components (which explain approximately 95% of the variance), thus obtaining 200-dimensional time series spanning 2727 consecutive days. Note that did not use any of the dataset's annotations, contrary to prior work such as Rayana & Akoglu 2015).
A.7 EXPERIMENTAL SETUP
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 de- viation equals 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 more sensitive to outliers, thus it would 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 the experiment,
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) use 40%. Such choice dictated by the particular (non uniform) experimental setup reported in (Braei &
23

---

Wagner 2020, Munir et al. 2019) and has been chosen to produce comparable results with state of the art methods present in literature.
Evaluation metries: We compare different predictors by means of the RMSE (root mean squared error) on the one- step ahead prediction errors. Given a sequence of data Y and the one-step ahead predictions YIN the RMSE is defined as: V Ei [iv(i) ((
As done in (Braei & Wagner 2020) we compare different anomaly detection methods taking into account several metrics. We use Fl-Score which is defined as the armonic mean of Precision and Recall (see (Braei & Wagner 2220; 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 rafe on different threshold values. We now follow (Braei & Wagner 2020) describe how AUC is computed. Let the true positive rate and false positive rate be defined, respectively, as: TPR = T and FPR FP where TP stands for true positive, P for positive, PP for false positive and N 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 I 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 random anomalous point in the time series a higher anomaly score than random normal point.
N
Hardware: We conduct our experiments on the following hardware setup:
Processor: Intel(R) Core(TM) 19-10980XE CPU @ 3.00GHz
RAM: 128 Gb
GPU: Nvidia TITAN 12Gb and Nvidia TITAN RTX 24Gb
Hypee--arameters: All the experiments we carried out are uniform on the optimization hyper-parameters. partic- ular 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: to = 10, 1 100, l2 200
linear filters kernel lengths (No. N1. N2): half predictor's '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 chose N₃- 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 T in comparable experimental setup. In particular, we chose comparable window lengths and architecture sizes (same order of mag- nitude 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
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
When computing the F-score we use the predictions of the CUMSUM detector which we collect as 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.
Anomaly scores:
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 Pr it is large when data does not come from Pp (which we consider the reference distribution).
pp"
24

---

Train Error
Test Error
Generalization Gap
RNSE on Vahco A1
Memory Train Error
Memory Test Error
Memory Generalization Gap
RASE on Yahoo
Wea 501000
Memory Train Error
Memory Test Error
Meroory Generalization Gap
To-aaaTT FEN STRE
PTE Sofwinn
RASE on Tahco
20
Memory
Memory Test Error
M Memory
Train
Generalization Gap
RNSE on Yahoo 0V
5.0
29
36
Memory Train Error
Memory Test Error
M Meroory
Generalization Gap
RNSE on C02 Dataset
Ing
Meroory Train Error
Mernory Test Error
M Mernory
RIIS MAB eeatroffic
Generalization Gap
Memory Train Error
Memory Test Error
M Merorr
RMSE on NAB reatweets
Generalization Gap
Memory
Memory
Mernory
Figure 9: Ablation studies on different datasets: Effects of interpretable blocks and fading regularization model's forecasting as the available window of past data increases (memory). Left Panel: Train error. Center Panel: Test error. Right Panel: Generalization Gap. The test error of STRIC is uniformiy smaller than standard TCN (without interpretable blocks nor fading regularization). Adding interpretable blocks to a standard TCN improves generalization for fixed memory w.r.t. Standard TCN but get worse (overfitting occurs) soon the available past data horizon increase. Fading regularization is effective: STRIC generalization GAP is almost constant w.r.t. past horizon.
25

---

A.7.1 ABLATION STUDY
In Figure 9 we show different metries based on the predictor's RMSE (training, test and generalization gap) as 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 10 different random seeds) the values of train and 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 T we show the effects on different choices of the predictor's memory Rored 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 nlee: the best results are achieve for small windows. On the other hand when Tidei is large the performance drops. This due to the type of anomalies present the Yahoo benchmark: most of the them can be considered to be point anomalies. In fact, window det is chosen.
showed in Appendix A.5 our detector is less sensible to point anomalies when a large
In Table we also report the reconstruction error of the optimal predictor given it's memory 7preed Note small memory in the predictor introduce modelling bias (higher training error) while a large memory does not (thanks to fading regularization). As observed in Appendix AA.5 better predictive models provide the detection module with more discriminative residuals: the downstream detection module achieves better F-scores and AUC.
Table 4: Ablation study on the RMSE of prediciton errors with standard deviation on 10 different seeds: We compare a standard TCN model with our STRIC predictor and some variation of it (using the same train hyper- parameters). The effect of adding a linear interpretable model before a TCN improves generalization error. Fading regularization has beneficial effect in controlling the complexity of the TCN model and reducing the generalization gap.
TCN
TCN- Lincar
TCN + Fading
STRIC pred
Train
Thin
Test
Yalooo A1
0.10 0.06
0.9 0.06
0.10 0.03
0.88 0.08
0.05
0.921:000 0:00
0.02
0.00
Dutosets
Yalhoo
0.02
082:2: 0.02
0.13 0.01
0.02
100:0000 0.01
+
0.330 + 0.01
Yahoo A3
0.13 1: 0.01
0:43 0.01
0.16 2 0.01
0.22 1 0.01
0.15 0.01
0.40 0:011:000
0.19 0.01
0.22 1: 0.01
Yalooo A4
0.15 1 0.01
061 1: 0.01
0.19 : 0.1
0.35 0:1
0.17 117:0000 0.01
0.0
0.23 0.01
0.24 2 0.01
0.00
0.62 0.02
0.02
0.05
0.0
0.01
NAB Triffic
0,03-1 1: 0.01
1.06 0.02
0.04 :: 0.01
1.00 0002
0.62:::00 0.01
0.93 0:011:0011
0.63
0.02
0.74 1: 0.02
NAR Tweets
0.18 1 0.05
1.02 2: 0.05
0.20 - 0.05
0.98 = 0.6
0.47 0.47 0.00
0.83 0.00
0.70 0.01
0.77 1 0.01
Table 5: Sensitivity STRIC to dyper-parameters We compare STRIC on different anomaly detection benchmarks datasets using different hyper-paramen papameearaetrrr:: memory of the predictor Npredd and length of anomaly detector windows
np fz Vddet-
Yahoo
Yahoo A2
Yahoo A3
Yahoo A4
F1
AUC
F1
AUC
F1
AUC
F1
AUC
Agredd 10, 2
0.45
0.89
0.63
0.99
0.87
0.99
0.64
0.89
Models
Vlyral 100, Adk 2
0.48
0.9308
0.98
0.9999
0.99
0.9999
0.68
0.9348
Tipral 10. faki 20
0.10
0.58
0.63
0.99
0.47
0.83
0.37
0.72
Tiyred = 100, Maki 20
0.10
0.55
0.98
0.9999
0.49
0.86
0.35
0.76
Yahoo A1
Yahoo A2
Yahoo A.
Yahoo A4
Train
Test
Train
Test
Train
Test
Train
Test
Models
Tepred 10
0.44
0.62
0.16
0.31
0.22
0.23
0.25
0.26
Tppred 100
0.42
0.61
0.14
0.30
0.19
0.22
0.23
0.24
A.7.2 COMPARISON TCN VS STRIC
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
26

---

compare the prediciton enrors of 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.
Tra dut
Dest data
Trm predictions
taed predictions
200
400
600
800
3000
1290
1400
1600
200
400
600
800
1000
1200
14400
1,600
200 400
300
800
1000
12:00
1400
2600
200 400
500
300
1020
1200
1:400
1500
Timustampa
Timestamps
Figure 10: We compare an off-the-shelf TCN against STRIC (time series predictor) the Yahoo dataset A3 Bench- mark. Note the standard TCN overfits compared to STRIC: the standard TCN does not handle correctly the trend component of the signal (First row). If we consider a time series without trend, the standard TCN model performs better but overfitting is still present. particular the generalization gap (measured using squared reconstruction error) for the two models is: Standard TCN 0.3735 and STRIC 0.0135.
TCN
STRIC
Trash data test: data traih prodiction
3-
bect condiction.
Wyy My My
M
5.0
600
700
800
90.0
1.000
1100
32200
500
600
7:00
8:00
9.00
10:00
11.00
1200
Timaatamms
Timastampa
Figure 11: Zoom on the second row of panels in Figure 10 We show the interface between train and test data both on plain TCN and on our STRIC predictor. plain TCN overfits w.r.t. STRIC also when not trend present.
A.8 STRIC SOTA ANOMALY DETECTORS
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 Table G 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 Both for Fl and AUC we report the following for each comparing SOTA method: AUGnvatat AAUCmm:: 100 (similarly for FI).
1-AUCsmIe
To begin with, STRIC outperforms 'traditional' methods (LOF and One-class SVM) which are considered as baselines models for comparing time series anomaly detectors.
27

---

Table 6: Comparison with SOTA anomaly defectors: We compare STRIC with other anomaly detection methods
the experimental setup and the same evaluation metrics proposed in (Braei & Wagner 2020; Munir et al. 2019).
The baseline models are: MA, ARIMA, LOF (Shen al. 2020), LSTM (Braei & Wagner 2020; Munir et al. 2019);
Wavenet (Braci& & Wagner 2020 Yahoo EGADS (Munir et al. 2019) GOAD (Bergman & Hoshen [2020). Om-
niAnomaly (Su et al. 2019). Twitter AD (Munir et 2019), TanoGAN (Bashar & Nayak, 2020), TadGAN (Geiger
al. 2020), DeepAR (Flunkert et al. (017) and DeepAn'T (Munir et al.
2019) STRIC outperforms most of the other
methods based on statistical models and based on DNNs. Same as Table
here we report the relative improvements
w.r.t. STRIC (the higher the better).
Relative Fi-score over
Yahoo A1
Yahoo A2
Yahoo A3
Yahoo A4
NAB Tweets
NAB Traffic
improvement
STRIC in %
ARIMA
-20
-88
-42
6
-33
-37
LSTM
-7
-33
-60
21
Yahoo EGADS
-1
-78
-54
OmniAnomaly
-1
-60
-45
-11
-6
-10
Twitter AD
0
-98
-85
-53
TanoGAN
-11
-85
-73
-13
-36
-44
Models
TadGAN
-13
-85
-65
-20
-25
-47
DeepAR
-29
-72
-79
-41
-37
-32
DeepAnT
-4
-67
-15
0
STRIC (ours)
0
0
0
0
0
0
Relative
AUC im- over
Yahoo A1
Yahoo A2
Yahoo A3
Yahoo A4
NAB Tweets
NAB Traffic
provement
STRIC in %
MA
-47
-98
-98
379
ARIMA
-5
-99
-99
124
LOF
-28
.99
-81
-32
-44
Wavenet
-60
-99
-99
-84
LSTM
-63
-99
-99
-84
Models
GOAD
-37
-99
-99
-51
-19
-12
DeepAnT
-32
-99
-99
-53
-23
-13
STRIC (ours)
Comparison with other Deep Learning based methods: STRIC outperforms most of the SOTA Deep Learning based methods reported in Table TadGAN, TAnoGAN, DeepAn'T and DeepAR (the last one is SOTA time se- ries 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 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 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.
A.8.1 DETAILS ON THE NYT EXPERIMENT
Figure 4 shows the normalized anomaly score computed by STRIC 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-point 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 the high variability of
28

---

1.0
9/11 attack
Election day
indian Ocean tsunami inauguration
0.8
0.8
0.6
0.
02-
0.0
Aug
04
Dec
Figure 12: A closer look at of the change-poing detected by STRIC. Left: Normalized anomaly score (blue line) and normalized frequency of the "Terrorism" descriptor (orange line) around the 9/1 I attack. Right: Normalized anomaly score (blue line) and normalized frequency of the "Earthquakes" descriptor (orange line) in the second half of 2004 and beginning of 2005. The 2004 U.S. election causes an increase in the anomaly score, but the most significant thange-point occurs after the Indian Ocean tsunami.
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 any human annotations. However, they can help interpreting the
found by STRIC.
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 abrupt increase of the frequency of descriptors like "Farthquakes" (orange line) and "Tsunami".
In Figure 12 (right),
not aware of any ground truth metries (e.g.. in Rayana & Akoglu (2015) a similar qualitative result has been reported on the NYT dataset). We therefore tested STRIC against simple baseline which uses PCA BERT features and threshold to detect anomalies. Despite being simple baseline, this method prooved to be highly applied in practice due to its simplicity (Blázquez-García et al. 2020). The PCA threshold baseline is able 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/periode components. For instance, the anomaly score of the simple baseline contains many false alarms which are related normal weekly periodicity that is not easily modeled by the baseline. This does not affect STRIC's predictions, since normal weekly periodicity directly modeled and identified as normal behaviour.
We note this experiment is qualitative and unfortunately
29


Reviews:
[{"review_id": "aX69n3oALPh", "reviewer": "Reviewer_zX4p", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Strengths: The regularization method introduced in Section 4 facilitates avoiding model overfitting. \n\nWeaknesses:\n\n1. Experimental results are not convincing. The chosen baselines are not state-of-the-art solutions for TS anomaly detection, please refer to [1] for details. The benchmark dataset used for comparison should also be changed. For example, as pointed out in [2], Yahoo dataset is not preferred for comparing AD solutions. \n\n2. TS decomposition itself is a well-studied research problem and it is essential to compare the proposed architecture with existing techniques (e.g., [3]) to obtain the residual. \n\n3. The paper is difficult to follow. Firstly, the overall workflow for anomaly detection is not clear, and it is difficult to see the connections between Section 4 and Section 5. Secondly, the model is not described clearly. For example, what is the non-linear module? Is it the TCN block shown in the figure without a title? Thirdly, there does not exist any statement on how the anomaly detector determines a particular time point or subsequence to be an anomaly. \n\n4. The paper is not focused and swings between prediction, anomaly detection, and even text embeddings. \n\n\nReferences: \n\n[1] Lai, Kwei-Herng, et al. \"Revisiting Time Series Outlier Detection: Definitions and Benchmarks.\" (2021).\n\n[2] Wu, Renjie, and Eamonn Keogh. \"Current time series anomaly detection benchmarks are flawed and are creating the illusion of progress.\" IEEE Transactions on Knowledge and Data Engineering (2021).\n\n[3] Wen, Qingsong, et al. \"RobustSTL: A robust seasonal-trend decomposition algorithm for long time series.\" Proceedings of the AAAI Conference on Artificial Intelligence. Vol. 33. No. 01. 2019.", "summary_of_the_paper": "This paper proposed a residual-style architecture namely STRIC for multivariate time series (TS) forecasting and anomaly detection, by first decomposing TS to trend, seasonality, and irregular components as residual, then extracting features from the residual using TCN, and finally introducing a likelihood ratio estimation method. Experimental results show that STRIC can provide interpretability for TS forecasting and outperform some existing anomaly detection models. ", "summary_of_the_review": "While some of the proposed techniques seem to be interesting, the paper is very difficult to follow and the experimental results are not convincing. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "WLSgUFVNDE", "reviewer": "Reviewer_rnBY", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "·\tStrength:\n\no\tThe motivation of each components are clearly justified.\n\no\tThe empirical result on univariate time series data outperforms most of the baselines.\n\no\tTo the best of my knowledge, the Bayesian learning formulation and CUMSUM algorithm with subspace likelihood ratio estimation are novel.\n\n·\tWeakness:\n\no\tEmpirical evaluation is only conducted on univariate time series data, making it not convincing that it is capable of dealing with multivariate time series data (which is very common in many application scenarios). Add some evaluations on multivariate time series may help on eliminate the concern.\n\no\tThe selected datasets has been shown as flawed benchmark [1] and can be easily solved with few lines of rules. Conducting experiment on these datasets making it questionable about the performance of the proposed algorithm. Either adopting new datasets to evaluate the algorithm or involving the one-line rule as baselines may be helpful to eliminate the concern.\n\no\tThe proposed residual framework assume that the input time series is the addition of trend, seasonality and residual series. However, it is very likely that the three components are not linearly dependent to the resulting time series. Providing more justifications on using addition to model aggregate the information may be helpful to eliminate the concern.\n\no\tThere is no empirical evaluation of the interpretable components. The author only show the result of decomposition without illustrating what makes it helpful for interpreting outlier detection result (which is one of the main contributions of this paper). Adding more justification with outlier annotations on Figure 2 may be helpful for readers to understand the contribution in empirical analysis.\n\n[1] Renjie Wu, Eamon Keogh, Current Time Series Anomaly Detection Benchmarks are Flawed and are Creating the Illusion of Progress, IEEE Transactions on Knowledge and Data Engineering, 2021\n", "summary_of_the_paper": "·\tAn interpretable stacked residual framework is proposed for time series anomaly detection with interpretation on seasonality and trend of the target time series.\n\n·\tA novel regularization term is proposed to optimize network parameters with Bayesian Optimization\n\n·\tAn extension of CUMSUM algorithm is proposed to detect outliers without knowledge of pre- and post-distributions of the data.\n", "summary_of_the_review": "·\tTo summarize, this paper proposed a novel interpretable outlier detection algorithm with stack residual blocks, regularized model complexity selection and novel CUMSUM detection algorithm. The proposed method is technically sound. However, the selecting datasets are solely in univariate setting and are shown to be flawed in previous study, which makes the empirical evaluation less credible. In addition, technical justification and empirical analysis of residual framework is missing, making it hard to evaluate the contribution of interpretable components. I believe some more modifications will certainly improve the quality of paper to meet the ICLR standard. Therefore, I vote for marginally below the acceptance threshold.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "TwPcQh5VCP", "reviewer": "Reviewer_TtBt", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Strengths:\n\n* The network architecture (Figure 1) is carefully designed to incorporate useful ideas from both conventional statistical time series models and DNN models to allow the proposed model to get the best of both worlds, including its interpretability and flexibility.\n\n* Although the regularization method for the TCN is inspired by the automatic relevance determination (ARD) method proposed more than two decades ago for (sparse) Bayesian learning of neural networks, its use here is novel as far as I know.\n\n* The nonparametric extension of the classical CUMSUM algorithm for anomaly detection does not require knowledge of the distributions of the two windows of data points, relaxing a major restriction of the original CUMSUM.\n\nWeaknesses:\n\n* Presentation of the experimental results (Table 1) can be improved. For example, some baseline models (e.g., ARIMA) are in fact univariate time series forecasting models. I suppose they handle each dimension of a multivariate time series separately. It would be clearer in the presentation to separate the models into two groups depending on whether they are really multivariate models. In fact, there also exist conventional multivariate time series forecasting models (e.g., VAR). Shouldn’t some of them be included as well? For clarity, it would help to show the dimensionality of each dataset as well, perhaps using a table to summarize the major characteristics (including dimensionality) of all the datasets used.\n\n* While the nonparametric extension of CUMSUM is good, the estimation error can be high (also with high variance) if the windows are not long enough, yet having too long windows may miss some small, subtle anomalies. This is a dilemma. A more detailed analysis is needed.\n\n* The experimental study is a bit limited. It is not clear whether the superior performance of the proposed method shown in Table 1 also holds for a wider range of real-world multivariate time series datasets.\n\nSome general comments:\n\n* You may consider how the proposed automatic complexity determination scheme is related to the attention mechanism.\n\n* For qualitative evaluation, it would help to demonstrate the efficacy of the automatic complexity determination scheme by visualizing the “relevant past” selected automatically.\n\n* In Figure 1, other than the captions of the sub-figures, there should also be a caption for the figure as well. At least it should be shown as “Figure 1” explicitly.\n\n* There are some minor language errors in the paper. It should be proofread carefully to correct them.", "summary_of_the_paper": "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.\n", "summary_of_the_review": "This is a good work which tries to integrate ideas from both conventional time series forecasting models and more recently DNN models. There are some novel ideas which are interesting and are worth studying further. Nevertheless, as mentioned above, the significance of the proposed model for multivariate time series anomaly detection is unknown due to the limitations of the empirical study.\n", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper studies the important problem of time series anomaly detection using deep neural networks (DNNs). Unlike many other DNN models, it focuses on incorporating in its model architecture interpretable components that are inspired by previous studies based on both conventional statistical methods and more recent DNN models.

While the paper has merits as pointed out by the reviewers (esp. TtBt), a number of concerns have also been raised, including the choice of datasets (e.g., by reviewers rnBY and zX4p). We appreciate the authors’ effort by adding some preliminary results of further experiments, but addressing all the concerns thoroughly will need a lot more work to get a scholarly paper that is more ready for publication. We believe this work has potential to be accepted for publication in a reputable venue if the concerns are thoroughly addressed after substantial revision.
Decision: Reject