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

1 INTRODUCTION

Time series data is being generated in increasing volumes from industrial, medical, commercial and scientific applications. Such growth is fueling demand for anomaly detection algorithms that are general enough to be applicable
across domains, yet reliable enough to operate on real-world time series data (Munir et al., 2019; Geiger et al., 2020;
Su et al., 2019). While recent developments have focused on deep neural networks (DNNs), simple linear models
still outperform DNNs in applications that require robustness (Braei & Wagner, 2020) and interpretable failure modes
(Geiger et al., 2020; Su et al., 2019).

To harvest the flexibility and interpretability of engineered modules while enabling end-to-end differentiable training,
we introduce STRIC: Stacked Residuals of Interpretable Components. We follow standard practice and consider a two
stage anomaly detection pipeline comprising a model of the normal time series and an anomaly detector based on the
prediction residuals (Munir et al., 2019). In particular, STRIC is composed of three modules: An interpretable local
_predictor, a non-linear predictor and a novel non-parametric anomaly detector._

More specifically, STRIC uses a parametric model implemented by a sequence of residuals blocks with each layer
capturing the prediction residual of previous layers. The first layer models trends, the second layer models quasiperiodicity/seasonality at multiple temporal scales, the third layer is a general linear predictor, and the last is a general
non-linear model in the form of a Temporal Convolution Network (TCN). While the first three layers are local (i.e.
applied to each component of the time series separately), the last integrates global statistics from additional time series
(covariates). Thanks to the residual structure the interpretable linear blocks do not reduce the representative power
of our non-linear architecture: any component of the input time-series which cannot be modeled by the interpretable
blocks is processed deeper into our architecture by the non-linear module (see Section 4). The model is trained end-toend with a prediction loss and we automatically select its complexity using an upper bound of the marginal likelihood
which, to the best of our knowledge, has never been applied to TCNs before.

Anomalies are detected by checking for time instants in which the prediction residual is not stationary. To avoid any
unrealistic assumption on the prediction residuals distribution, we use a likelihood ratio test that we derive from a
variational upper bound of f -divergences and that can be computed directly from the data points.

To summarize, our contributions are:


-----

1. We introduce STRIC, a stacked residual model that explicitly isolates interpretable factors such as slow
trends, quasi-periodicity, and linearly predictable statistics (Oreshkin et al., 2019; Cleveland et al., 1990), and
incorporates statistics from other time series as context/side information.

2. We introduce a novel regularization that is added to the prediction loss and is used for automatic model
complexity selection according to the Empirical Bayes framework (Rasmussen & Williams, 2006).

3. We introduce a non-parametric extesion of the CUMSUM algorithm which entails a tunable parameter corresponding to the length of observation and enables anomaly detection in the absence of knowledge about the
pre- and post-distributions.

4. We test our method on standard anomaly detection benchamrks and show it merges both the advantages of
simple and interpretable linear models and the flexibility of non-linear ones while discounting their major
drawbacks: lack of flexibility of linear models and lack of interpretability and overfitting of non-linear ones.

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 to asynchronous time processes in our review. Different methods for time series
anomaly detection (TSAD) can be taxonomized by their choice of (i) discriminant function, (ii) continuity criterion,
and (iii) optimization method to determine the tolerance threshold. It is common to use statistics of the prediction
error as the discriminant (Braei & Wagner, 2020), and the likelihood ratio between the distribution of the prediction
error before and after a given time instant as the continuity criterion (Yashchin, 1993). Recent methods compute the
discriminant using deep neural network architectures and euclidean distance as continuity criterion (Munir et al., 2019;
Geiger et al., 2020; Su et al., 2019; Bashar & Nayak, 2020).

Our method follows a similar line but introduces novel elements both in (i) and (ii): (i) the discriminant function is
the prediction error residual of a novel regularized stacked residual architecture; (ii) the decision function is based on
a novel non-parametric extension of the CUMSUM algorithm. The resulting method, STRIC, has the advantage of
separating interpretable components due to trends and seasonality, without reducing the representative power of our
architecture. At initialization, the system is approximately equivalent to a multi-scale SARIMA model (Adhikari &
Agrawal, 2013), which can be reliably applied out-of-the-box on most time series. However, as more data is acquired,
any part of the system can be further fine-tuned in an unsupervised end-to-end fashion.

Munir et al. (2019) argue that anomaly detection can be solved by exploiting a flexible model provided a proper
inductive bias is introduced (e.g. TCN). In Appendix A.7.2 we show that TCN alone might overfit simple time series.
We therefore take their direction a step further, and, differently from previous works (Bai et al., 2018; Munir et al.,
2019; Sen et al., 2019; Geiger et al., 2020), we provide our temporal model with an interpretable structure, similar to
Oreshkin et al. (2019). Moreover, unlike previous works on interpretability of DNNs (Tsang et al., 2018; Guen et al.,
2020), our architecture explicitly imposes both an inductive bias and a regularization which are designed to expose
the user both a STL-like decomposition (Cleveland et al., 1990) and the relevant time scale of the signals. Since
TCNs tend to overfit if not properly regularized (Appendix A.7.2), we constrain our TCN’s representational power by
enforcing fading memory (Zancato & Chiuso, 2021) while retaining what is needed to predict future values.

Our method outperforms both classical statistical methods (Braei & Wagner, 2020) and deep networks (Munir et al.,
2019; Geiger et al., 2020; Su et al., 2019; Bergman & Hoshen, 2020; Bashar & Nayak, 2020) on different anomaly
detection benchmarks (Laptev & Amizadeh, 2020; Lavin & Ahmad, 2015) (Section 6). Moreover, we show it can be
employed to detect anomalous patterns on complex data such as text embeddings of newspaper articles (Figure 4).

3 NOTATION

We denote vectors with lower case and matrices with upper case. In particular y is multi-variate time series {y(t)}t∈Z,
_y(t)_ R[n]; we stack observations from time t to t + k 1 and denote the resulting matrix as Yt[t][+][k][−][1] := [y(t), y(t +
_∈_ _−_
1), ..., y(t + k − 1)] ∈ R[n][×][k]. The row index refers to the dimension of the time series while the column index refers
to the temporal dimension. We denote the i-th component of the time series y as yi and its evaluation at time t ad as
_yi(t) ∈_ R. We refer to {y(s), s > t} as the test/future and to {y(s), s ≤ _t} as the reference/past intervals. At time t,_
sub-sequences containing the np past samples up to time t _np + 1 are given by Yt[t]_ _np+1_ [(note that we include the]
_−_ _−_


-----

**past** **future**


**Layer inputs** **For each**

_i = 1,...,n_ _zi_

_Z ∈{XTREND, XSEAS, XLIN}_ _z1_

_Z =_ _z⋮2_ _φ1 * zi_

Φ ∈{𝒦Layer filters banksTREND, 𝒦SEAS, 𝒦LIN} Φ = _φφz⋮n12_ _φφ2l *⋮ * z zii_

**Layer features selectors** _φn_ **Time Features** _Gi_

_A ∈{𝒜TREND, 𝒜SEAS, 𝒜LIN}_ _A =_ _aa⋮12[T][T]_ _ai[T][G]i_ ∈ ℝ[l][×][n][p]

_an[T]_

∈ ℝ[1×][n][p]

**Layer predictors**

_B ∈{ℬTREND, ℬSEAS, ℬLIN}_ _b1[T]_

_B =_ _b⋮2[T]_ **Future Predictor** _ai[T][G]i[b]i_ ∈ ℝ

_bn[T]_


_XSEAS ∈ℝ[n][×][n][p]_


_XLIN ∈ℝ[n][×][n][p]_


(b) Interpretable blocks structure. Time features are
extracted independently for each time series (see Appendix A.1 for more details).


_XTREND ∈ℝ[n][×][n][p]_

|n t −n+ 1 t|Col2|⋮|
|---|---|---|
|p Y t−n+1||t t + 1|
||p y(t + 1) yTREND(t + 1) ∈ℝ ySEAS(t + 1) ∈ℝ yLIN(t + 1) ∈ℝ yTCN(t Trend Seasonal Linear TCN XTREND XSEAS XLIN||


(a) STRIC predictor architecture.


present data into the past data), while future samples up to time t + nf are Yt[t]+1[+][n][f] . 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 1a. Its basic building blocks are causal convolutions (Bai et al., 2018), with
a fixed-size 1-D kernel with input elements from time t and earlier. Rather than initializing the convolutional filter
randomly, as commonly done in deep learning, we initialize the weights so that each layer is biased to attend at
different components of the signal, as explained in the following.

**Linear module. The first (linear) module is interpretable and captures local statistics of a given time series by means**
of a cascade of learnable linear filters. Its first layer models and removes slow-varying components in the input
data. We initialize the filters to mimic a causal Hodrick Prescott (HP) filter (Ravn & Uhlig, 2002). The second layer
models and removes periodic components: it is initialized to have a periodic impulse response. Finally, the third layer
implements a linear stationary filter bank.

We treat the impulse responses parameters of the linear filters as trainable parameters. To allow our model to work on
a wide variety of time scales, we initialize the trend layer with different HP smoothness degrees, while we initialize
the periodic and linear-stationary layers with randomly chosen poles (Farahmand et al., 2017) both on the unit circle
and within the unit circle, thus allowing to capture different periodicities.

**Non-linear module. The second (non-linear) module aggregates global statistics from different time series using**
a TCN model (Sen et al., 2019). It takes as input the prediction residual of the linear module and outputs a matrix
_G(Yt[t]−np+1[)][ ∈]_ [R][l][×][n][p][ where][ l][ is the number of output features extracted by the TCN model. The column][ G][(][Y][ t]t−np+1[)][j]
with j = 1, ..., np of the non-linear features is computed using data up to time t − _np + j (due to the internal structure_
of a TCN network (Bai et al., 2018)). We build a linear predictor on top of G(Yt[t]−np+1[)][ for each single time series]
independenty: the predictor for the i-th time series is given by: ˆyTCN(t + 1)i := a[T]i _[G][(][Y][ t]t_ _np+1[)][b][i]_ [where][ a][i]
_−_ _[∈]_ [R][l][ and]
_brelevant features across time indices to build the one-step ahead predictor (see Appendix A.1).i ∈_ R[n][p] . Since ai combines features (uniformly in time) we can interpret it as a feature selector. While bi aggregates

Note that the third layer of the linear module is a superset of preceding ones, and the non-linear module is a superset
of the whole linear module. While this makes the model redundant, we show that this design, coupled with proper
initialization and regularization, improves the reliability and intepretability of the final model. We improve filters
optimization by sharing their kernel parameters among different time series so that global information (e.g., common
trend shapes, periodicities) can be extracted. In Appendix A.1, we describe each component of the model in detail.


-----

4.1 AUTOMATIC COMPLEXITY DETERMINATION


Consider the TCN-based future predictor be ˆyTCN(t + 1) := a[T] _G(Yt[t]−np+1[)][b][ = ˆ]XTCNb where_ _X[ˆ]TCN ∈_ R[1][×][n][p] is the
output of the TCN block which depends on the the past window Yt[t]−np+1 [of length][ n][p] [(the memory of the predictor).][1]

Ideally, np should be large enough to capture the “true” memory of the time series, but should not be too large if not
necessary (i.e. bias variance trade-off). In practice however, flexible feature extractors (such as TCNs or plain DNNs)
are prone to overfitting. Therefore, some regularization is needed to control model complexity and benefit from having
a large memory window. In this section, we introduce a novel regularized loss inspired by Bayesian arguments which
allows us to use an architecture with a “large enough” past horizon np (i.e., larger than the true system memory)
and automatically select the relevant past to avoid overfitting. Such information is exposed to the user through an
interpretable parameter λ that directly measures the relevant time scale of the signal.

**Bayesian learning formulation: We model the innovations (optimal prediction errors) as Gaussian, so that y(t + 1) |**
modeling assumption does not restrict our framework and is used only to justify the use of the squared loss to learnYt[t]−np+1 _[∼N]_ [(][F][ ∗][(][Y][ t]t−np+1[))][, η][2][)][ where][ F][ ∗] [is the optimal predictor of the future values given the past. Note that this]
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(Yt[t]+1[+][n][f] _b, W, Yt[t]_ _np+1[) =][ Q][n]k=1[f]_ _[p][(][y][(][t][ +][ k][)][ |][ b, W, Y][ t]t+[+]k[k][−]n[1]p_ [)][.]
_|_ _−_ _−_

To make the notation simpler, we shall denote with Yf := Yt[t]+1[+][n][f] R[n][f] the set of future outputs over which the
_∈_
predictor is computed and we shall use _Y[ˆ]b,W ∈_ R[n][f] as the predictor’s outputs.

In a Bayesian framework, the optimal set of parameters can be found maximizing the posterior p(b, W _Yf_ ) over the
_|_
model parameters. We model b and W as independent random variables:
_p(b, W_ _Yf_ ) _p(Yf_ _b, W_ )p(b)p(W ) (1)
_|_ _∝_ _|_
where p(b) is the prior associated to the predictor coefficients and p(W ) is the prior on the remaining parameters. The
prior p(b) should encode our belief that the prediction model should not be too complex and should depend only on the
_most relevant past. We model this by assuming that the components of b have zero mean and exponentially decaying_
variances: Eb[2]np _j_ 1 [=][ κλ][j][ for][ j][ = 0][, ..., n][p][ −] [1][, where][ κ][ ∈] [R][+][ and][ λ][ ∈] [(0][,][ 1)][. The corresponding maximum]
_−_ _−_

entropy prior pλ,κ(b) (Cover & Thomas, 1991) under such constraints is log(pλ,κ(b)) ∝−∥b∥Λ[2] _[−][1][ −]_ [log(][|][Λ][|][)][ where]
Λ ∈ R[n][p] is a diagonal matrix with elements Λj,j = κλ[j] with j = 0, ..., np. Here, λ represents how fast the output
of the predictor “forgets” the past. Therefore, λ regulates the complexity of the predictor: the smaller λ, the lower the
complexity.

In practice, λ has to be estimated from the data. One would be tempted to estimate jointly b, W, λ, κ (and possibly η)
by minimizing the negative log of the joint posterior (see Appendix A.2.1). Unfortunately, this leads to a degeneracy
since the joint negative log posterior goes to −∞ when λ → 0. Indeed, typically the parameters describing the prior
(such as λ) are estimated by maximizing the marginal likelihood, i.e., the likelihood of the data once the parameters
(b, W ) have been integrated out. Since computing (or even approximating) the marginal likelihood in this setup is
prohibitive, we now introduce a variational upper bound to the marginal likelihood which is easier to estimate.

**Variational upper bound to the marginal likelihood: The model structure we consider is linear in b and we can**
therefore stack the predictions of each available time index t to get the following linear predictor on the whole future
data available: _Y[ˆ]b,W = FW b where F ∈_ R[n][f][ ×][n][p] is obtained by stacking _X[ˆ]TCN(Yi[i]−np+1[)][ for][ i][ =][ t, ..., t][ +][ n][f][ −]_ [1][.]

**Proposition 4.1. Consider a model on the form:** _Y[ˆ]b,W = FW b (linear in b and possibly non-linear in W_ _) and its_
_posterior in Equation (1). Assume the prior on the parameters b is given by the maximum entropy prior and W is_
_fixed. Then the following is an upper bound on the marginal likelihood associated to the posterior in Equation (1)_
_with marginalization taken only w.r.t. b:_

2

_Ub,W,Λ = η[1][2]_ _Yf −_ _Y[ˆ]b,W_ + b[⊤]Λ[−][1]b + log det(FW ΛFW[⊤] [+][ η][2][I][)][.] (2)

This (proved in Appendix A.2) provides an alternative loss function to the negative log posterior which does not suffer
of the degeneracy alluded above while optimizing over b, W, λ and κ.

1For simplicity we consider scalar time series, but our approach easily generalizes to multivariate time series (Appendix A.2.2).


-----

**Remark: We use batch normalization (Ioffe & Szegedy, 2015) to normalize the output of FW along its rows so that**
features have comparable scales, this avoids the TCN network to counter the fading regularization by increasing its
output scales (see Appendix A.2.3).

5 ANOMALY DETECTOR

In this section, we present our anomaly detection method based on a variational approximation of the likelihood ratio
between two windows of model residuals. Our temporal residual architecture model produces the prediction residual
after removing trends, periodicity, and stationary (linear) components, as well as considering global covariates. Such a
prediction residual is used to test the hypothesis that the time instant t is anomalous by comparing its statistics before
_t on temporal windows of length np and nf_ . The detector is based on the likelihood ratios aggregated sequentially
using the classical CUMSUM algorithm (Page, 1954; Yashchin, 1993). CUMSUM, however, requires knowledge
of the distributions, which we do not have. Unfortunately, the problem of estimating the densities is hard (Vapnik,
1998) and generally intractable for high-dimensional time series (Liu et al., 2012). We circumvent this problem by
directly estimating the likelihood ratio with a variational characterization of f -divergences (Nguyen et al., 2010) which
involves solving a convex risk minimization problem in closed form.

In Section 5.1, we summarize the standard material necessary to derive our new estimator and the resulting anomaly
test. The overall method is entirely unsupervised, and users can tune the scale parameter (corresponding to the window
of observation when computing the likelihood ratios) and the coefficient of CUMSUM (depending on the application
and desired operating point in the tradeoff between missed detection and false alarms).

5.1 LIKELIHOOD RATIOS AND CUMSUM

CUMSUM (Page, 1954) is a classical Sequential Probability Ratio Test (SPRT) (Basseville & Nikiforov, 1993; Liu
et al., 2012) of the null hypothesis H0 that the data after the given time c comes from the same distribution as before,
against the alternative hypothesis Hc that the distribution is different. We denote the distribution before c as pp and
the distribution after the anomaly at time c as pf .

If the density functions pp and pf were known (we shall relax this assumption later), the optimal statistic to decide
whether a datum y(i) is more likely to come from one or the other is the likelihood ratio s(y(i)). According to the
Neyman-Pearson lemma, H0 is accepted if the likelihood ratio s(y(i)) is less than a threshold chosen by the operator,
otherwise Hc is chosen. In our case, the competing hypotheses are H0 = “no anomaly has happened” and Hc = “an
anomaly happened at time c”. We denote with pH0 and pHc the p.d.f.s under H0 and Hc so that: pH0 (Y1[K][) =][ p][p][(][Y][ K]1 [)]
and pHc (Y1[c][−][1]) = pp(Y1[c][−][1]), pHc (Yc[K] _Y1[c][−][1]) = pf_ (Yc[K] _Y1[c][−][1]). Therefore the likelihood ratio is:_
_|_ _|_

1 [)] 1 )pf (Yc[t] 1 ) _c_ 1 ) _t_ _pf_ (y(i) _Y1[i][−][1])_
Ω[t]c [:=][ p][H][c] [(][Y][ t] _[|][ Y][ c][−][1]_ = _[p][f]_ [(][Y][ t] _[|][ Y][ c][−][1]_ _|_ (3)

_pH0_ (Y1[t][) =][ p][p][(][Y][ c][−][1]pp(Y1[t][)] _pp(Yc[t]_ 1 ) [=] _i=c_ _pp(y(i)_ _Y1[i][−][1])_

_[|][ Y][ c][−][1]_ Y _|_

To determine the presence of an anomaly, we can compute the cumulative sum Sc[t] [:= log Ω]c[t] [of the (log) likelihood]
ratios, which depends on the time c, and estimate c[∗] using a maximum likelihood criterion, corresponding to the
detection function ht = max1 _c_ _t Sc[t][. The first instant at which we can confidently assess the presence of a change]_
_≤_ _≤_
of the detector depending on the application. The final estimatepoint (a.k.a. stopping time) is: cstop = min{t : ht ≥ _τ_ _} where τ ˆc is a design parameter that modulates the sensitivity of the true change point c[∗]_ after the detection cstop
is simply given by the timestamp c at which the maximum of ht = max1≤c≤t Sc[t] [is achieved. In Appendix A.3,]
we provide an alternative derivation that shows that CUMSUM is a comparison of the test statistic with an adaptive
threshold that keeps complete memory of past ratios. The next step is to relax the assumption of known densities,
which we bypass in the next section by directly approximating the likelihood ratios to compute the cumulative sum.

5.1.1 LIKELIHOOD RATIO ESTIMATION WITH PEARSON DIVERGENCE

The goal of this section is to tackle the problem of estimating the likelihood ratio of two general distributions pp and
_pf given samples. To do so, we leverage a variational approximation of f_ -divergences (Nguyen et al., 2010) whose
optimal solution is directly connected to the likelihood ratio. For different choices of divergence function, different
estimators of the likelihood ratio can be built. We focus on a particular divergence choice, the Pearson divergence,
since it provides a closed form estimate of the likelihood ratio (see Appendix A.4).


-----

**Proposition 5.1. (Nguyen et al., 2010; Liu et al., 2012) Let φ := pf** _/pp be the likelihood ratio of the unknown_
_distributions pf and pp. Let F := {fi : fi ∼_ _pf_ _, i = 1, ..., nf_ _} and H := {hi : hi ∼_ _pp, i = 1, ..., np} be two sets_
_containing nf and np samples i.i.d. from pf and pp respectively. An empirical estimator_ _φ[ˆ] of the likelihood ratio φ is_
_given by the solution to the following convex optimization problem:_


_np_

_φ(hi)[2]_
_−_ _n[1]f_
_i=1_

X


_nf_

_φ(fi)_ (4)
_i=1_

X


_φˆ = arg min_


2np


**Proposition 5.2. (Liu et al., 2012; Kanamori et al., 2009) Let φ in Equation (4) be chosen in the family of Reproducing**
_Kernel Hilbert Space (RKHS) functions Φ induced by the kernel k. Let the kernel sections be centered on the set of_
_data Str := {F, H} and let the kernel matrices evaluated on the data from pf and pp be Kf := K(F, Str) and_
_Kp := K(H, Str). The optimal regularized empirical likelihood ratio estimator on a new datum e is given by:_

_φˆ(e) =_ _[n][p]_ _K(e,_ _tr)_ _Kp[T]_ _[K][p]_ [+][ n][p][γI][n]p[+][n]f _−1Kf[T]_ [1][.] (5)

_nf_ _S_
 

**Remark: The estimator in Equation (5) is not constrained to be positive. Nonetheless, the positivity constraints can**
be enforced. In this case, the closed form solution is no longer valid but the problem remains convex.

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 Y1[t] [by looking at the]
prediction residuals E1[t][, which provide a sufficient representation of][ Y][ t]1 [(see Appendix A.5). We therefore assume]
we are given the prediction errors E1[t] [obtained from our time series predictor (the predictor should model the normal]
behaviour). This guarantees that the sequence E1[t] [is white in each of its normal subsequences. On the other hand, if]
the model is applied to a data subsequence which contains the abnormal condition, the residuals are correlated.

We estimate the likelihood ratio of pf and pp on a datum et as _φ[ˆ]t(et)._ _φˆt is obtained by applying Equation (5)_
on the past window of size np + nf . At each time instant t, we compute the necessary kernel matrices as
_Kf_ (Et[t] _nf +1[, E]t[t]_ _np_ _nf +1[)][ and][ K][p][(][E]t[t][−]n[n]p[f]_ _nf +1[, E]t[t]_ _np_ _nf +1[)][.]_
_−_ _−_ _−_ _−_ _−_ _−_ _−_

**Remark: At time t, the likelihood ratio is estimated assuming i.i.d. data. This assumption holds if no anomaly**
happened but does not hold in the abnormal situation since residuals are not i.i.d. In Appendix A.5, we prove that
treating correlated variables as uncorrelated provides a lower bound on the actual cumulative sum of likelihood ratios.
For a fixed threshold, this means the detector cumulates less and therefore requires more time to reach the threshold.

Finally, we compute the detector function by aggregating the estimated likelihood ratios: _S[ˆ]c[t]_ [:=][ P]i[t]=c [log ˆ]φi(ei).

**Remark: The choice of the windows length (np and nf** ) is fundamental and highly influences the likelihood estimator.
Using small windows makes the detector highly sensible to point outliers, while larger windows are better suited to
estimate sequential outliers (see Appendix A.5).

6 EXPERIMENTAL RESULTS

In this section, we show STRIC can be successfully applied to detect anomalous behaviours on different anomaly
detection benchmarks. In particular, we test our novel residual temporal structure, the automatic complexity regularization and the anomaly detector on the following datasets: Yahoo (Laptev & Amizadeh, 2020), NAB (Lavin &
Ahmad, 2015), CO2 Dataset (see Appendix A.6). In addition, to show the general applicability and flexibility of our
method, we test STRIC on the challenging task of detecting anomalous events in time series generated from embeddings of articles from the New York Times. See Appendix A.7 for the experimental setup and data normalization.

**Anomaly detection: While recent works show deep learning models are not well suited to solve AD on standard**
anomaly detection benchmarks (Braei & Wagner, 2020), we prove deep models can be effective provided they are used
with a proper inductive bias and regularization. In Table 1, we compare STRIC against statistical and deep learning
based anomaly detection methods. Our experiments follow the experimental setup and evaluation criteria used in Braei
& Wagner (2020) and Munir et al. (2019). Note no other method performs consistently (across different datasets) as
good as STRIC. In particular, STRIC achieves the best F1 score on Yahoo A3, in Appendix A.7 we show this is mainly


-----

Table 1: Comparison with SOTA anomaly detectors: We compare STRIC with other anomaly detection methods
(see Appendix A.8) on 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 et al., 2020), LSTM (Braei & Wagner,
2020; Munir et al., 2019), Wavenet (Braei & Wagner, 2020), Yahoo EGADS (Munir et al., 2019), 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., 2017) and DeepAnT (Munir et al., 2019) . STRIC
outperforms most of the other methods based on statistical models and based on DNNs. See Table 6 for the same table
obtained by looking at the relative performance w.r.t. STRIC.

**F1-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
Twitter AD **0.48** 0 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
DeepAnT 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
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

**Models**

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. In fact, most of the time series in Yahoo A3 are characterized by trend components and
seasonalities which STRIC’s interpretable predictor can easily model (see Appendix A.7.2). In Appendix A.7, we
show some ablation studies on the effects of STRIC’s hyper-parameters on its performance. In particular, we find that
STRIC is highly affected by the choice of the length of the windows used to estimate the likelihood ratio, while not
being much sensitive to the choice of the memory of the predictor (Appendix A.7.1). Interestingly, STRIC does not
achieve the optimal F1/AUC compared to linear models on Yahoo A4. The ability of linear models to outperform nonlinear ones on Yahoo A4 is known in the literature (e.g. in Geiger et al. (2020) any non-linear model is outperformed
by AR/MA models of the proper complexity). The main motivation for this is that modern (non-linear) methods tend
to overfit on Yahho A4 and therefore generalization is usually low. Instead, thanks to fading regularization and model
architecture, STRIC does not exhibit overfitting despite having larger complexity than SOTA linear models used in
A4. To conclude, we believe that STRIC merges both the advantages of simple and interpretable linear models and the
flexibility of non-linear ones while discounting their major drawbacks: lack of flexibility of linear models and lack of
interpretability and overfitting of non-linear ones (see Appendix A.8 for a more in depth discussion).

**STRIC interpretable time series decomposition: In Figure 2, we show STRIC’s interpretable decomposition. We**
report predicted signals (first row), estimated trends (second row) and seasonalities (third row) for different datasets.
For all experiments, we plot both training data (first 40% of each time series) and test data. Note the interpretable
components of STRIC generalize outside the training data, thus making STRIC work well on non-stationary time
series (e.g. where the trend component is non negligible and typical non linear models overfit, see Appendix A.7.2).

**Ablation study: We now compare the prediction performance of a general TCN model with our STRIC method in**
which we remove the interpretable module and the fading regularization one at the time. In Table 2, we report the test
RMSE prediction errors and the RMSE generalization gap (i.e. difference between test and training RMSE prediction
errors) for different datasets while keeping all the training parameters the same (e.g. training epochs, learning rates
etc.) and model parameters (e.g. nb = 100). The addition of the linear interpretable model before the TCN slightly
improves the test error. We note this effect is more visible on A2, A3, A4, mainly due to the non-stationary nature of
these datasets and the fact that TCNs do not easily approximate trends (Braei & Wagner, 2020) (we further tested this


-----

CO2 Dataset


Yahoo A3 Dataset


NAB Dataset


4 2

train data 2
test data 1

2 train predictions

test predictions 0 0

0 1

STRIC's predictions 2

2

100 150 200 250 300 350 400 450 500 660 680 700 720 740 760 780 800 350 400 450 500 550 600 650 700

4 2

2

1

2

Trend 0 0

0 1

2

2

100 150 200 250 300 350 400 450 500 300 400 500 600 700 800 900 1000 1100 350 400 450 500 550 600 650 700

4 2

2

2 1

0 0

0

Seasonality 1

2 2

2

200 220 240 260 280 300 320 340 650 675 700 725 750 775 800 825 850 300 350 400 450 500 550 600 650 700

Timestamps Timestamps Timestamps


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-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
a standard TCN model with our STRIC predictor and some variation of it (using the same training hyper-parameters).
Standard deviations are given in Table 4.

**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 A2 0.82 0.71 0.35 0.22 0.71 0.50 **0.30** **0.16**
Yahoo A3 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**

**Datasets** 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.74** **0.11**
NAB Tweets 1.02 0.84 0.98 0.78 0.83 0.36 **0.77** **0.07**

in Appendix A.7.1). While STRIC generalization is always better than a standard TCN model and STRIC’s ablated
components, we note that applying Fading memory regularization alone to a standard TCN does not always improve
generalization (but never decreases it): this highlights that the benefits of combining the linear module and the fading
regularization together are not a trivial ‘sum of the parts’. Consider for example Yahoo A1: STRIC achieves 0.62 test
error, the best ablated model (TCN + Linear) 0.88, while TCN + Fading does not improve over the baseline TCN. A
similar observation holds for the CO2 Dataset. Fading regularization might not be beneficial (nor detrimental) for time
series containing purely periodic components which correspond to infinite memory systems (systems with unitary
fading coefficient). In such cases the interpretable module is essential in removing the periodicities and providing
the regularized non-linear module (TCN + Fading) with an easier to model residual signal. We refer to Figure 2 (first
column) for a closer look on a typical time series in CO2 dataset, which contains a periodic component that is captured
by the seasonal part of the interpretable model. To conclude, our proposed fading regularization has (on average) a
beneficial effect in controlling the complexity of a standard TCN model and reduces its generalization gap (≈ 40%
reduction). Moreover, coupling fading regularization with the interpretable module guarantees the best generalization.

**Automatic complexity selection: In Figure 3, we test the effects of our automatic complexity selection (fading 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 is uniformly smaller than a standard TCN (without interpretable blocks nor fading regularization). Adding interpretable blocks to a standard TCN improves generalization
for a fixed memory w.r.t. standard TCN but gets worse (overfitting occurs) as soon as the available past data horizon
increases. On the other hand, the generalization gap of STRIC does not deteriorate as the memory of the predictor
increases (see Appendix A.7.1 for a comparison with other metrics).


-----

1.0

0.8

0.6

0.4

0.2


0 20 Memory40 60 80 100

Standard TCN
STRIC No Fading
STRIC


Figure 3: Automatic complexity selec**tion: Fading memory regularization pre-**
serves generalization gap as the memory
of the predictor np increases on NAB
Tweets.


Figure 4: Anomaly score on the New York Times dataset. Our
method finds anomalies in 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.


**Anomaly detection on the New York Times dataset: We qualitatevly test STRIC on a time series consisting of**
BERT embeddings (Devlin et al., 2019) of New York Times articles (Sandhaus, 2008) from 2000 to 2007. We set
_np = nf = 30 days, to be able to detect change-point anomalies that altered the normal distribution of news articles_
for a prolonged period of time. Without any human annotation, STRIC is able to detect major historical events such
as the 9/11 attack, the 2004 Indian Ocean tsunami, and U.S. elections (Figure 4). Note that we do not carry out a
quantitative analysis of STRIC’s predictions, as we are not aware of any ground truth or metrics for this benchmark,
see for example Rayana & Akoglu (2015). Additional details and comparison with a baseline model built on PCA are
given in Appendix A.8.1.

7 DISCUSSION AND CONCLUSIONS


We have shown that our interpretable stacked residual architecture and our unsupervised estimation of the likelihood
ratio are well suited to solve AD for multivariate time series data. Unlike purely DNN-based methods (Geiger et al.,
2020; Munir et al., 2019; Bashar & Nayak, 2020), STRIC exposes to the user both an interpretable STL-like time series
decomposition (Cleveland et al., 1990) and the relevant time scale of the time series. Both the interepretable module
and the fading memory regularization are important in building a successful model (see Table 2). In particular, the
interpretable module helps STRIC generalize correctly on non-stationary time series on which standard deep models
(such as TCNs) may overfit (Braei & Wagner, 2020). Moreover, we showed that our novel fading regularization alone
can improve the generalization error of TCNs up to ≈ 40% over standard TCNs, provided the periodic components
of the time series are captured by the interpretable module (Section 6). We highlight that the overall computational
complexity and memory requirement of our method remains the same as standard TCNs, so that our approach can
easily scale to large scale time series datasets.

An anomaly is a time instant: at that time instant, either we receive an isolated observation that is inconsistent with
normal operation (outlier measurement), or a discrete change occurs in the mechanism that generates the data (changepoint) that persists beyond that time instant (Geiger et al., 2020; Basseville & Nikiforov, 1993). Our method treats
these two phenomena in a unified manner, without the need to differentiate between outliers and setpoint changes,
with specialized detectors for each. Once the predictor is built, our method can be used online to detect anomalies
soon after occurrence without waiting for the entire data stream to be observed and without requiring any knowledge
on the prediction error distribution (nominal and faulty).

Making the non-parametric anomaly detector fully adaptive to the data is an interesting research direction: While our
fading window regularizer automatically tunes the predictor’s window length by exploiting the self-supervised nature
of the prediction task, methods to automatically tune the detector’s window lengths (an unsupervised problem) are
an interesting research direction. Moreover, designing statistically optimal rules to calibrate our detector’s threshold
_τ depending on the desired operating point in the tradeoff between missed detection and false alarms would further_
enhance the out-of-the-box robustness of our method.


-----

**Reproducibility Statement: Our method has been tested on publicly available datasets: Yahoo (Laptev & Amizadeh,**
2020), NAB (Lavin & Ahmad, 2015), CO2 Dataset (see Appendix A.6) and NYT (Sandhaus, 2008). We described the
data splitting and the data processing steps in Appendix A.7. We describe the major details regarding the implementation of our novel method in Appendix A.1 while in Appendix A.7 we describe both the model structure and training
hyper-parameters we used in the experimental section. We do not include the model structures and hyper-parameters
of SOTA methods we used as baselines and refer to related literature (referenced in our work) for the implementation
details. To further foster reproducibility we shall make our code available.

REFERENCES

Ratnadip Adhikari and R. K. Agrawal. An introductory study on time series modeling and forecasting. _CoRR,_
[abs/1302.6613, 2013. URL http://arxiv.org/abs/1302.6613.](http://arxiv.org/abs/1302.6613)

Shaojie Bai, J Zico Kolter, and Vladlen Koltun. An empirical evaluation of generic convolutional and recurrent
networks for sequence modeling. arXiv preprint arXiv:1803.01271, 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`ele Basseville and Igor V. Nikiforov. Detection of Abrupt Changes: Theory and Application. Prentice-Hall, Inc.,
USA, 1993. ISBN 0131267809.

Liron Bergman and Yedid Hoshen. Classification-based anomaly detection for general data. In International Confer_[ence on Learning Representations, 2020. URL https://openreview.net/forum?id=H1lK_lBtvS.](https://openreview.net/forum?id=H1lK_lBtvS)_

Ane Bl´azquez-Garc´ıa, Angel Conde, Usue Mori, and Jose A Lozano. A review on outlier/anomaly detection in time
series data. arXiv preprint arXiv:2002.04236, 2020.

Mohammad Braei and Sebastian Wagner. Anomaly detection in univariate time-series: A survey on the state-of-the-art.
_[CoRR, abs/2004.00433, 2020. URL https://arxiv.org/abs/2004.00433.](https://arxiv.org/abs/2004.00433)_

Robert B. Cleveland, William S. Cleveland, Jean E. McRae, and Irma Terpenning. Stl: A seasonal-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 Chang, Kenton Lee, and Kristina Toutanova. BERT: pre-training of deep bidirectional transformers 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-HLT 2019, Minneapolis, MN, USA, June 2-7, 2019, Volume 1 (Long and Short_
_Papers), pp. 4171–4186. Association for Computational Linguistics, 2019. doi: 10.18653/v1/n19-1423. URL_
[https://doi.org/10.18653/v1/n19-1423.](https://doi.org/10.18653/v1/n19-1423)

Amir-massoud Farahmand, Sepideh Pourazarm, and Daniel Nikovski. Random projection filter bank
for time series data. In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (eds.), _Advances in Neural Information Processing Systems,_ volume 30.
Curran Associates, Inc., 2017. [URL https://proceedings.neurips.cc/paper/2017/file/](https://proceedings.neurips.cc/paper/2017/file/ca3ec598002d2e7662e2ef4bdd58278b-Paper.pdf)
[ca3ec598002d2e7662e2ef4bdd58278b-Paper.pdf.](https://proceedings.neurips.cc/paper/2017/file/ca3ec598002d2e7662e2ef4bdd58278b-Paper.pdf)

Valentin Flunkert, David Salinas, and Jan Gasthaus. Deepar: Probabilistic forecasting with autoregressive recurrent
[networks. CoRR, abs/1704.04110, 2017. URL http://arxiv.org/abs/1704.04110.](http://arxiv.org/abs/1704.04110)

Alexander Geiger, Dongyu Liu, Sarah Alnegheimish, Alfredo Cuesta-Infante, and Kalyan Veeramachaneni. Tadgan:
Time series anomaly detection using generative adversarial networks. arXiv preprint arXiv:2009.07769, 2020.

Vincent Le Guen, Yuan Yin, J´er´emie Dona, Ibrahim Ayed, Emmanuel de B´ezenac, Nicolas Thome, and Patrick
Gallinari. Augmenting physical models with deep networks for complex dynamics forecasting. arXiv preprint
_arXiv:2010.04456, 2020._

10


-----

Sergey Ioffe and Christian Szegedy. Batch normalization: Accelerating deep network training by reducing internal
[covariate shift. CoRR, abs/1502.03167, 2015. URL http://arxiv.org/abs/1502.03167.](http://arxiv.org/abs/1502.03167)

Takafumi Kanamori, Shohei Hido, and Masashi Sugiyama. A least-squares approach to direct importance estimation.
_[Journal of Machine Learning Research, 10(48):1391–1445, 2009. URL http://jmlr.org/papers/v10/](http://jmlr.org/papers/v10/kanamori09a.html)_
[kanamori09a.html.](http://jmlr.org/papers/v10/kanamori09a.html)

Nikolay Laptev and Saeed Amizadeh. Yahoo! webscope dataset ydata-labeled-time-series-anomalies-v1 0. CoRR,
[2020. URL https://webscope.sandbox.yahoo.com/catalog.php?datatype=s&did=70.](https://webscope.sandbox.yahoo.com/catalog.php?datatype=s&did=70)

Alexander Lavin and Subutai Ahmad. Evaluating real-time anomaly detection algorithms - the numenta anomaly
[benchmark. CoRR, abs/1510.03336, 2015. URL http://arxiv.org/abs/1510.03336.](http://arxiv.org/abs/1510.03336)

Song Liu, Makoto Yamada, Nigel Collier, and Masashi Sugiyama. Change-point detection in 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-3.

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):5847–5861, Nov 2010. ISSN
[1557-9654. doi: 10.1109/tit.2010.2068870. URL http://dx.doi.org/10.1109/TIT.2010.2068870.](http://dx.doi.org/10.1109/TIT.2010.2068870)

Boris N. Oreshkin, Dmitri Carpov, Nicolas Chapados, and Yoshua Bengio. N-BEATS: neural basis expansion analysis
[for interpretable time series forecasting. CoRR, abs/1905.10437, 2019. URL http://arxiv.org/abs/1905.](http://arxiv.org/abs/1905.10437)
[10437.](http://arxiv.org/abs/1905.10437)

Ewan S 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 hodrick-prescott filter for the frequency of observations. The Review
_of Economics and Statistics, 84:371–375, 02 2002. doi: 10.1162/003465302317411604._

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.9781611974010.70. URL https://doi.org/10.1137/1.9781611974010.70.](https://doi.org/10.1137/1.9781611974010.70)

Evan Sandhaus. The new york times annotated corpus ldc2008t19. web download. Linguistic Data Consortium,
_Philadelphia, 6(12):e26752, 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´e-Buc, E. Fox, and R. Garnett (eds.), Advances in Neural Information Processing Systems, vol[ume 32. Curran Associates, Inc., 2019. URL https://proceedings.neurips.cc/paper/2019/file/](https://proceedings.neurips.cc/paper/2019/file/3a0844cee4fcf57de0c71e9ad3035478-Paper.pdf)
[3a0844cee4fcf57de0c71e9ad3035478-Paper.pdf.](https://proceedings.neurips.cc/paper/2019/file/3a0844cee4fcf57de0c71e9ad3035478-Paper.pdf)

Lifeng Shen, Zhuocong Li, and James Kwok. Timeseries anomaly detection using temporal hierarchical 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 https://proceedings.neurips.cc/paper/2020/file/](https://proceedings.neurips.cc/paper/2020/file/97e401a02082021fd24957f852e0e475-Paper.pdf)
[97e401a02082021fd24957f852e0e475-Paper.pdf.](https://proceedings.neurips.cc/paper/2020/file/97e401a02082021fd24957f852e0e475-Paper.pdf)

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. 2828–2837, New York, NY, USA, 2019._
[Association for Computing Machinery. ISBN 9781450362016. doi: 10.1145/3292500.3330672. URL https:](https://doi.org/10.1145/3292500.3330672)
[//doi.org/10.1145/3292500.3330672.](https://doi.org/10.1145/3292500.3330672)

11


-----

Michael E Tipping. Sparse bayesian learning and the relevance vector machine. Journal of machine learning research,
1(Jun):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[ume 31. Curran Associates, Inc., 2018. URL https://proceedings.neurips.cc/paper/2018/file/](https://proceedings.neurips.cc/paper/2018/file/74378afe5e8b20910cf1f939e57f0480-Paper.pdf)
[74378afe5e8b20910cf1f939e57f0480-Paper.pdf.](https://proceedings.neurips.cc/paper/2018/file/74378afe5e8b20910cf1f939e57f0480-Paper.pdf)

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´emi 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-art 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 Linguistics. URL https://www.aclweb.org/anthology/2020.emnlp-demos.6.](https://www.aclweb.org/anthology/2020.emnlp-demos.6)

Emmanuel Yashchin. Performance of cusum control schemes for serially correlated observations. Technometrics, 35
[(1):37–52, 1993. ISSN 00401706. URL http://www.jstor.org/stable/1269288.](http://www.jstor.org/stable/1269288)

Luca Zancato and Alessandro Chiuso. A novel deep neural network architecture for non-linear system identification.
_IFAC-PapersOnLine, 54(7):186–191, 2021. ISSN 2405-8963. doi: https://doi.org/10.1016/j.ifacol.2021.08.356._
[URL https://www.sciencedirect.com/science/article/pii/S2405896321011307.](https://www.sciencedirect.com/science/article/pii/S2405896321011307) 19th
IFAC Symposium on System Identification SYSID 2021.

A APPENDIX

A.1 IMPLEMENTATION

Given two scalar time series x and ϕ, we denote their time convolution g as: g(t) := (ϕ _x)(t) :=_ _i=_
_∗_ _−∞_ _[ϕ][(][i][)][x][(][t][−][i][)][.]_
We say that g is the causal convolution of x and ϕ if ϕ(t) = 0 for t < 0, so that the output g(t) does not depend
on future values of x (w.r.t. current time index t). In the following, we shall give a particular interpretation to the

[P][∞]
signals x and ϕ: x will be the input signal to a filter parametrized by an impulse response ϕ (kernel of the filter). Note
any (causal) convolution is defined by an infinite summation for each given time t. Therefore it is customary, when
implementing convolutional filters, to consider a truncated sum of finite length. In practice, this is obtained assuming
the filter impulse response is non-zero only in a finite window. The truncation is indeed an approximation of the
complete convolution, nonetheless it is possible to prove that the approximation errors incurred due to truncation are
guaranteed to be bounded under simple assumptions. To summarize, in the following we shall write g(t) := (ϕ _∗_ _x)(t)_
and mean that the impulse response of the causal filter ϕ is truncated on a window of a given length.

A.1.1 ARCHITECTURE

main blocks of the architecture are defined to encode trend, seasonality, stationary linear and non-linear part. In theLet Yt[t]−np+1 _[∈]_ [R][n][×][n][p][ be the input data to our architecture at time step][ t][ (a window of][ n][p] [past time instants). The]
following we shall denote each quantity related to a specific layer using either the subscripts {TREND, SEAS, LIN,
TCN} or {0, 1, 2, 3}.

We shall denote the input of each block as Xk ∈ R[n][×][n][p] and the output as _X[ˆ]k ∈_ R[n][×][n][p] for k = 0, 1, 2, 3. The residual
architecture we propose is defined by the following: X0 = Yt[t]−np+1 [and][ X][k] [=][ X][k][−][1] _[−]_ _X[ˆ]k−1 for k = 1, 2, 3. At each_
layer we extract lk temporal features from the input Xk. We denote the temporal features extracted from the input of
the k-th block as: Gk := Gk(Xk) ∈ R[l][k][×][n][p] . The i-th column of the feature matrix Gk is a feature vector (of size lk)
extracted from the inputset of filter banks (Bai et al., 2018). Xk up to time t − _np −_ _i. To do so, we use causal convolutions of the input signal Xk with a_

12


-----

A.1.2 INTERPRETABLE RESIDUAL TCN ON SCALAR TIME SERIES

**Modeling Interpretable blocks: In this section, we shall describe the main design criteria of the linear module. For**
each interpretable layer (TREND, SEAS, LIN), we convolve the input signal with a filter bank designed to extract
specific components of the input.

For example, consider the trend layer, denoting its scalar input time series by x and its output by gTREND. Then
_gTREND is defined as a multidimensional time series (of dimension lTREND := l0) obtained by stacking l0 time series_

[and parametrize each filter ingiven by the convolution ofϕTREND1 _∗x, ϕTREND2_ _∗x, ..., ϕ x K withTRENDTREND l0 causal linear filters:l with its truncated impulse response (i.e. kernel) of length1 ∗x][T]_ . We denote the set of linear filters ϕTRENDi ∗ _x for i = 0 ϕ, ..., lTREND0 −i for1. In other words, i = 0 k, ..., l0 :=0 k−TREND1 g asTREND K._ TREND :=

We interpret each time series in gTREND as an approximation of the trend component of x computed with the i-th filter.
We design each ϕTRENDi so that each filter extracts the trend of the input signal on different time scales (Ravn & Uhlig,
2002) (i.e., each filter outputs a signal with a different smoothness degree). We estimate the trend of the input signal by
recombining the extracted trend components in gTREND with the linear map aTREND. Moreover, we predict the future
trend of the input signal (on the next time-stamp) with the linear map bTREND.

We construct the blocks that extract seasonality and linear part in a similar way.

**Implementing Interpretable blocks: The input of each layer is given by a window of measurements of length np.**
We zero-pad the input signal so that the convolution of the input signal with the i-th filter is a signal of length np (note
this introduces a spurious transient whose length is the length of the filter kernel). We therefore have the following
temporal feature matrices: G0 = GTREND ∈ R[l][0][×][n][p], G1 = GSEAS ∈ R[l][1][×][n][p] and G2 = GLIN ∈ R[l][2][×][n][p] .

The output of each layer _X[ˆ]k is an estimate of the trend, seasonal or stationary linear component of the input signal on_
the past interval of length np, so that we have _X[ˆ]k ∈_ R[1][×][n][p] (same dimension as the input Xk). On the other hand, the
linear predictor ˆyk computed at each layer is a scalar. Intuitively, _X[ˆ]k and ˆyk 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:Note ak combines features (uniformly in time) so that we can interpret it as a feature selector whileX[ˆ]k := a[T]k _[G][k][ and][ ˆ]yk := X[ˆ]kbk for k = 0, 1, 2 where ai ∈_ R[l][k] and bk aggregates bk ∈ R[n][p] .
relevant features across different time indices to build the one-step ahead predictor. Depending on the time scale of
the signals it is possible to choose bk depending on the time index (similarly to the fading memory regularization).
We experimented both the choice to make bk canonical “vectors” and dense vectors. We found that choosing bk as
canonical vectors, whose non-zero entry is associated to the closest to the present time instat provides good empirical
results on most cases.

**Non-linear module The non-linear module is based on a standard TCN network. Its input is defined as X3 =**
_Yt[t]−np+1_ _[−]_ _X[ˆ]0 −_ _X[ˆ]1 −_ _X[ˆ]2, which is to be considered as a signal whose linearly predictable component has been_
removed. The TCN extracts a set of l3 non-linear features G3(X3) ∈ R[l][3][×][n][p] which we combine with linear maps
as done for the previous layers. The j-th column of the non-linear features G3 is computed using data up to time
_t −_ _np + j (due to the internal structure of a TCN network (Bai et al., 2018)). The linear predictor on top of G3 is_
_yˆTCN := a[T]3_ _[G][3][b][3][, where][ a][3]_

_[∈]_ [R][l][3][ and][ b][3] _[∈]_ [R][n][p] [.]

Finally, the output of our time model is given by:


_a[T]k_ _[G][k][(][X][k][)][b][k][.]_
_k=0_

X


_yˆ(t + 1) =_


_yˆk =_
_k=0_

X


_Xˆkbk =_
_k=0_

X


A.1.3 INTERPRETABLE RESIDUAL TCN ON MULTI-DIMENSIONAL TIME SERIES

We extend our architecture to multi-dimensional time series according to the following principles: preserve interpretability (first module) and exploit global information to make local predictions (second module).

In this section, the input data to our model isseries). _Yt[t]−np+1_ _[∈]_ [R][n][×][n][p][ (a window of length][ n][p] [from an][ n][-dimensional time]

**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: KTREND, KSEAS and KLIN. For ease of notation we shall now focus only on the trend layer. Any
other layer is obtained by substituting ‘TREND’ with the proper subscript (‘SEAS’ or ‘LIN’).

We denote by GTREND, i ∈ R[l][0][×][n][p] the set of time features extracted by the trend filter bank KTREND from the i-th
time series. Each feature matrix is then combined as done in the scalar setting using linear maps, which we now index
by the time series index i: aTRENDi and bTRENDi. The rationale behind this choice is that each time series can exploit
differently the extracted features. For instance, slow time series might need a different filter than faster ones (chosen
using aTRENDi) or might need to look at values further in the past (retrieved using bTRENDi). We stack the combination
vectors aTRENDi and bTRENDi into the following matrices: ATREND = [aTREND1, aTREND2, ..., aTRENDn][T] _∈_ R[n][×][l][0] and
_BTREND = [bTREND1, bTREND2, ..., bTRENDn][T]_ _∈_ R[n][×][n][p] .

**Non-linear module: The second (non-linear) module aggregates global statistics from different time series (Sen**
et al., 2019) using a TCN model. It takes as input the prediction residual of the linear module and outputs a matrix
_GTREND(Yt[t]−np+1[)][ ∈]_ [R][l][3][×][n][p][ where][ l][3] [is the number of output features that are extracted by the TCN model (which]
is a design parameter). The j-th column of the non-linear features GTREND(Yt[t]−np+1[)][ is computed using data up to]
time t _p + j, where p is the “receptive” field of the TCN (p < np). This is due to the internal structure of a TCN_
_−_
network (Bai et al., 2018) which relies on causal convolutions and typically scales as O(2[h]) where h is the number of
TCN hidden layers (the deeper the TCN the longer its receptive field). As done for the time features extracted by the
interpretable blocks, we build a linear predictor on top of GTREND(Yt[t]−np+1[)][ for each single time series independenty:]
the predictor for the i-th time series is given by: ˆyTCN(t+1)i := a[T]i _[G][TREND][(][Y][ t]t_ _np+1[)][b][i]_ [where][ a][i]
_−_ _[∈]_ [R][l][3][ and][ b][i] _[∈]_ [R][n][p] [.]
We stack the combination vectors aTCNi and bTCNi into the following matrices: ATCN = [aTCN1, aTCN2, ..., aTCNn][T] _∈_
R[n][×][l][3] and BTCN = [bTCN1, bTCN2, ..., bTCNn][T] _∈_ R[n][×][n][p] .

Finally, the outputs of the predictor on the i-th time series are given by:


_yˆ(t + 1)i =_


_akTi_ _[G][ki][b][ki]_ [+][ a][TCN]Ti _[G][TCN][b][TCN][i][.]_
_k∈{TRENDX,SEAS,LIN}_


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 theand block’s kernel size, respectively (without loss of generality, we assume all filters have the same dimension). Each k-th block’s kernel parameters Kk ∈ R[l][k][×][N][k], where lk and Nk are the block’s number of filters
filter (causal 1D-convolution) is parametrized by the values of its impulse response parameters (kernel parameters).
When we learn a filter bank, we mean that we optimize over the kernel values for each filter jointly. For multidimensional time series, we apply the filter banks to each time series independently (depth-wise convolution) and improve
filter learning by sharing kernel parameters across different time series.

**Initialization: The first block (trend) is initialized using l0 causal Hodrick Prescott (HP) filters (Ravn & Uhlig, 2002)**
of kernel size N0. HP filters are widely used to extract trend components of signals (Ravn & Uhlig, 2002). In general
a HP filter is used to obtain from a time series a smoothed curve which is not sensitive to short-term fluctuations and
more sensitive to long-term ones (Ravn & Uhlig, 2002). In general, a HP filter is parametrized by a hyper-parameter
_λHP which defines the regularity of the filtered signal (the higher λHP, the smoother the output signal). We initialize_
each filter with λHP chosen uniformly in log-scale between 10[3] and 10[9]. Note the impulse response of these filters
decays to zero (i.e., the latest samples from the input time series are the most influential ones). When we learn the
optimal set of trend filter banks, we do not consider them parametrized by λHP and search for the optimal λHP. Instead,
we optimize over the impulse response parameters of the kernel which we do not assume live in any manifold (e.g.,
the manifold of HP filters). Since this might lead to optimal filters which are not in the class of HP filters, we impose
a regularization which penalizes the distance of the optimal impulse response parameters from their initialization.

The second block (seasonal part) is initialized using l1 periodic kernels which are obtained as linear filters whose poles
(i.e., frequencies) are randomly chosen on the unit circle (this guarantees to span a range of different frequencies). Note
the impulse responses of these filters do not go to zero (their memory does not fade away). Similarly to the HP filter
bank, we do no optimize the filters over frequencies, but rather we optimize them over their impulse response (kernel
parameters). This optimization does not preserve the strict periodicity of filters. Therefore, in order to keep the optimal
impulse response close to initialization values (purely periodic), we exploit weight regularization by penalizing the
distance of the optimal set of kernel values from initialization values.

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 l2 increases, this random filter bank is
guaranteed to be a universal approximator of any (stationary) linear system (see (Farahmand et al., 2017) for details).

**Remark: This block could approximate any trend and periodic component. However, we assume to have factored out**
both trend and periodicities in the previous blocks.

The last module (non-linear part) is composed by a randomly initialized TCN model. We employ a TCN model due to
its flexibility and capability to model both long-term and short-term non-linear dependencies. As is standard practice,
we exploit dilated convolutions to increase the receptive field and make the predictions of the TCN (on the future
horizon) depend on the most relevant past (Bai et al., 2018).

**Remark: Our architecture provides an interpretable justification of the initialization scheme proposed for TCN in**
(Sen et al., 2019). In particular our convolutional architecture allows us to handle high-dimensional time series data
without a-priori standardization (e.g., trend or seasonality removal).

A.2 AUTOMATIC COMPLEXITY DETERMINATION (FADING MEMORY REGULARIZATION)

In this section, we shall introduce a regularization scheme called fading regularization, to constrain TCN representational capabilities.

The output of the TCN model is G(Yt[t]−np+1[)][ ∈] [R][l][3][×][n][p][ where][ l][3] [is the number of output features extracted by the]
TCN model. The predictor build from TCN features is given by: aTCN[T]i _[G][TCN][(][Y][ t]t_ _np+1[)][b][TCN][i][, where the predictor]_
_−_
_bnon-linear featuresTCNi ∈_ R[n][p] takes as input a linear combination of the TCN features (weighted by G(Yt[t]−np+1[)][ is computed using data up to time][ t][ −] _[n][p]_ [+][ j][ (due to causal convolutions used in the] aTCNi). The j-th column of the
internal structure of the TCN network (Bai et al., 2018)). One expects that the influence on the TCN predictor as j
increases should increase too (in case j = np, the statistic is the one computed on the closest window of time w.r.t.
present time stamp). Clearly, the exact relevance on the output is not known a priori and needs to be estimated. In
other words, the predictor should be less sensitive to statistics (features) computed on a far past, a property which is
commonly known as fading memory. Currently, this property is not built in the predictor bTCNi, which treats each
time instant equally and might overfit while trying to explain the future by looking into far and possibly non-relevant
past. In order to constrain model complexity and reduce overfitting, we impose the fading memory property on our
predictor by employing a specific regularization which we now describe.

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 np measures can be written
as: ˆyTCN(t + 1) = a[T] _GTCN(Tt[t]−np+1[)][b][ = ˆ]Xkb. We shall assume that innovations (optimal prediction errors) are_
Gaussian, so thatvalues given the past. Note that this assumption does not restrict our framework and is used only to justify the use y(t + 1) | Yt[t]−np+1 _[∼N]_ [(][F][ ∗][(][Y][ t]t−np+1[))][, η][2][)][, where][ F][ ∗] [is the optimal predictor of the future]
of the squared loss to learn the regression function of the predictor. In practice, we do not know the optimal F _[∗]_
and we approximate it with our parametric model. For ease of exposition, we group all the architecture parameters
except b in the weight vector W (linear filters parameters KTREND, KSEAS, KLIN, linear module recombination weights
_ATREND, ASEAS, ALIN, BTREND, BSEAS, BLIN, and TCN kernel parameters and recombination coefficients ATCN etc.)._
We write the conditional likelihood of the future given the past data of our parametric model as:

_nf_

_p(Yt[t]+1[+][n][f]_ _b, W, Yt[t]_ _np+1[) =]_ _p(y(t + k)_ _b, W, Yt[t]+[+]k[k][−]n[1]p_ [)] (6)
_|_ _−_ _|_ _−_

_k=1_

Y

To make the notation simpler, we shall denote by Yf := Yt[t]+1[+][n][f] R[n][f] the set of future outputs over which the predictor
_∈_
is computed and we shall useconditioning pastnf _Yt[t]−np+1_ [(which is present in any conditional distribution). Equation (6) becomes:]Y[ˆ]b,W ∈ R[n][f] as the predictor’s outputs. Moreover, we shall drop the dependency on the[ p][(][Y][f] _[|][ b, W]_ [) =]
_k=1_ _[p][(][y][(][t][ +][ k][)][ |][ b, W]_ [)][. The optimal set of parameters][ b][∗] [and][ W][ ∗] [in a Bayesian framework is computed by]
maximizing the posterior on the parameters given the data:
Q

_p(b, W_ _Yf_ ) _p(Yf_ _b, W_ )p(b)p(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 be too high and therefore it should only depend on the most
_relevant past._

**Remark: The prior does not induce hard constraints. It rather biases the optimal predictor coefficients towards the**
prior belief. This is clear by looking at the negative log-posterior which can be directly interpreted as the loss function
to be minimized: log p(b, W _Yf_ ) = log p(Yf _b, W_ ) log p(b) log p(W ). In particular, the first term
log p(Yf _b, W_ ) is the data fitting term (only influenced by the data). Both − _|_ _−_ _|_ _−_ _− log p(b) and log p(W_ ) do not depend on
the available data and can be interpreted as regularization terms that bound the complexity of the predictor function. |

The main idea is to reduce the sensitivity of the predictor on time instants that are far in the past. We therefore enforce
the fading memory assumption on p(b) by assuming that the components of b ∈ R[n][p] have zero mean and exponentially
decaying variances:
Ebj = 0 and Eb[2]np−j−1 [=][ κλ][j][ for][ j][ = 0][, ..., n][p] _[−]_ [1] (8)

where κ ∈ R[+] and λ ∈ (0, 1). Note the larger variance (larger scale) is associated to temporal indices close to the
present time t.

**Remark: To specify the prior, we need a density function p(b) but up to now we only specified constraints on the first**
and second order moments. We therefore need to constrain the parametric family of prior distributions we consider.
Any choice on the class of prior distributions lead to different optimal estimators. Among all the possible choices of
prior families we choose the maximum entropy prior (Cover & Thomas, 1991). Under constraints on first and second
moment, the maximum entropy family of priors is the exponential family (Cover & Thomas, 1991). In our setting, we
can write it as:
log pλ,κ(b) ∝−∥b∥Λ[2] _[−][1][ −]_ [log][ |][Λ][|] (9)
where Λ ∈ R[n][p][×][n][p] is a diagonal matrix whose elements are Λj,j = κλ[j] for j = 0, ..., np − 1.

The parameter λ represents how fast the predictor’s output ‘forgets’ the past: the smaller λ, the lower the complexity.
In practice, we do not have access to this information and indeed we need to estimate λ from the data.

One would be tempted to estimate jointly W, b, λ, κ (and possibly η) by minimizing the negative log of the joint
posterior:

1 2

arg min _Yf_ _Yb,W_ + log(η[2]) log(pλ,κ(B)) log(p(W )). (10)
_b,W,λ,κ_ _η[2]_ _−_ [ˆ] _−_ _−_

Unfortunately, this leads to a degeneracy since the joint negative log posterior goes to −∞ when λ → 0.

**Bayesian learning formulation for fading memory regularization:**

The parameters describing the prior (such as λ) are typically estimated by maximizing the marginal likelihood, i.e., the
likelihood of the data once the parameters (b, W ) have been integrated out. Unfortunately, the task of computing (or
even approximating) the marginal likelihood in this setup is prohibitive and one would need to resort to Monte Carlo
sampling techniques. While this is an avenue worth investigating, we preferred to adopt the following variational
strategy inspired by the linear setup.

Indeed, the model structure we consider is linear in b and we can therefore stack the predictions of each available time
index t to get the following linear predictor on the whole future data: _Y[ˆ]b,W = FW b where FW ∈_ R[n][f][ ×][n][p] and its rows
are given by _X[ˆ]TCN(Yi[i]−np+1[)][ for][ i][ =][ t, ..., t][ +][ n][f][ −]_ [1][.]

We are now ready to find an upper bound to the marginal likelihood associated to the posterior given by Equation (7)
with marginalization taken only w.r.t. b.
**Proposition A.1 (from (Tipping, 2001)). The optimal value of a regularized linear least squares problem with feature**
_matrix F and parameters b is given by the following equation:_


1

_f_ [Σ][−][1][Y][f] (11)
_η[2][ ∥][Y][f][ −]_ _[Fb][∥][2][ +][ b][⊤][Λ][−][1][b][ =][ Y][ ⊤]_


arg min


_with Σ := F_ ΛF _[⊤]A[T]_ + η[2]I.

Equation (11) guarantees that


1

_f_ [Σ][−][1][Y][f] [+ log][ |][Σ][|][,]
_η[2][ ∥][Y][f][ −]_ _[Fb][∥][2][ +][ b][⊤][Λ][−][1][b][ + log][ |][Σ][| ≥]_ _[Y][ ⊤]_

16


-----

where the right hand side is (proportional to) the negative marginal likelihood with marginalization taken only w.r.t. b.
Therefore, for fixed a W,
1 2

_η[2]_ _Yf −_ _Y[ˆ]b,W_ + b[⊤]Λ[−][1]b + log |FW ΛFW[⊤] [+][ η][2][I][|]

is an upper bound of the marginal likelihood with marginalization over b and does not suffer of the degeneracy alluded
at before.

With this considerations in mind, and inserting back the optimization over W, the overall optimization problem we
solve is

1 2

arg min _Yf_ _Yb,W_ + _b_ Λ[−][1][ + log][ |][F][W][ Λ][F][W][ +][ η][2][I][|][ + log][ p][(][W] [)] (12)
_b,W,λ_ (0,1),κ>0 _η[2]_ _−_ [ˆ] _∥_ _∥[2]_
_∈_

**Remark: log p(W** ) defines the regularization applied on the remaining parameters of our architecture. In particular,
we induce sparsity by applying L[1] regularization on ATREND, ASEAS, ALIN and ATCN. Also, we constrain filters
parameters to stay close to initialization by applying L[2] regularization on KTREND, KSEAS and KLIN.

A.2.2 FADING MEMORY IN MULTIVARIATE TIME SERIES

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

A.2.3 FEATURES NORMALIZATION

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

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.
Before a change point, by construction we are under the distribution of the past. Therefore, log _[p]p[f]p[ (](y[y])[)]_

_[≤]_ [0][ ∀][y][, which]
in turn means that the cumulative sum S1[t] [will decrease as][ t][ increases (negative drift). After the change, the situation is]
opposite and the cumulative sum starts to show a positive drift, since we are sampling y(i) from the future distribution
_pf_ . This intuitive behaviour shows that the relevant information to detect a change point can be obtained directly from
the cumulative sum (along timestamps). In particular, all we need to know is the difference between the value of the
cumulative sum of log-likelihood ratios and its minimum value.

The CUMSUM algorithm can be expressed using the following equations: _vt_ := S1[t] :=
minj,1 _j_ _t Sj[t][. The stopping time is defined as:][ t][stop][ = min][{][t][ :][ v][t][ ≥]_ _[τ]_ _[}][ = min][{][t][ :][ S]1[t]_ _[−]_ _[m][t][, where][+][ τ]_ _[}][. With the][ m][t]_
_≤_ _≤_
last equation, it becomes clear that the CUMSUM detection equation is simply a comparison of the cumulative sum of[≥] _[m][t]_
the log likelihood ratios along time with an adaptive threshold mt +τ . Note that the adaptive threshold keeps complete
memory of the past ratios. The two formulations are equivalent because S1[t] [=][ h][t][.]

_[−]_ _[m][t]_

A.4 VARIATIONAL APPROXIMATION OF THE LIKELIHOOD RATIO

In this section, we present some well known facts on f -divergences and their variational characterization. Most of the
material and the notation is from (Nguyen et al., 2010). Given a probability distribution P and a random variable f
measurable w.r.t. P, we use _fdP to denote the expectation of f under P. Given samples x(1), ..., x(n) from P, the_
R

17


-----

empirical distribution Pn is given by Pn = _n[1]_ _ni=1_ _[δ][x][(][i][)][. We use]_ _fdPn as a convenient shorthand for the empirical_

expectation _n[1]_ _ni=1_ _[f]_ [(][x][(][i][))][.]
P R

Consider two probability distributionsP P and Q, with P absolutely continuous w.r.t. Q. Assume moreover that both
distributions are absolutely continuous with respect to the Lebesgue measure µ, with densities p0 and q0, respectively,
on some compact domain X ⊂ R[d].

**Variational approximation of the f-divergence: The f** -divergence between P and Q is defined as (Nguyen et al.,
2010)

_q0_
_Df_ (P, Q) := _p0f_ _dµ_ (13)

_p0_

Z  

where f : R → R is a convex and lower semi-continuous function. Different choices of f result in a variety of
divergences that play important roles in various fields (Nguyen et al., 2010). Equation (13) is usually replaced by the
variational lower bound:



[φdQ − _f_ _[∗](φ)dP]_ (14)


_Df_ (P, Q) ≥ _φsup∈Φ_


and equality holds iff the subdifferential ∂f ( _p[q][0]0_ [)][ contains an element of][ Φ][. Here][ f][ ∗] [is defined as the convex dual]

function of f .

In the following, we are interested in divergences whose conjugate dual function is smooth (which in turn defines
commonly used divergence measures such as KL and Pearson divergence), so that we shall assume that f is convex
and differentiable. Under this assumption, the notion of subdifferential is not required and the previous statement reads
as: equality holds iff ∂f ( _p[q][0]0_ [) =][ φ][ for some][ φ][ ∈] [Φ][.]

**Remark:** The infinite-dimensional optimization problem in Equation (14) can be written as Df (P, Q) =
supφ∈Φ EQφ − EPf _[∗](φ)._

In practice, one can have an estimator of any f -divergence restricted to a functional class Φ by solving Equation (14)
(Nguyen et al., 2010). Moreover, when P and Q are not known one can approximate them using their empirical
counterparts: Pn and Qn. Then an empirical estimate of the f -divergence is: _D[ˆ]_ _f_ (P, Q) = supφ∈Φ EQn _φ_ _−_ EPn _f_ _[∗](φ)._

**Approximation of the likelihood ratio: 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) is achieved iff
_φ = ∂f_ ( _p[q][0]0_ [)][. This tells us that the optimal solution to the variational approximation provides us with an estimator of]

the composite function ∂f ( _p[q][0]0_ [)][ of the likelihood ratio][ q]p[0]0 [. As long as we can invert][ ∂f] [, we can uniquely determine the]

likelihood ratio.

In the following, we shall get an empirical estimator of the likelihood ratio in two separate steps. We first solve the
following:

_φˆn := arg maxφ∈Φ_ EQn _φ −_ EPn _f_ _[∗](φ)_ (15)

which returns an estimator of ∂f ( _p[q][0]0_ [)][, not the ratio itself. And then we apply the inverse of][ ∂f][ to][ ˆ]φn. We therefore

have a family of estimation methods for the likelihood function by simply ranging over choices of f .

**Remark: If f is not differentiable, then we cannot invert ∂f but we can obtain estimators of other functions of the**
likelihood ratio. For instance, we can obtain an estimate of the thresholded likelihood ratio by using a convex function
whose subgradient is the sign function centered at 1.

A.4.1 LIKELIHOOD 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, or χ[2], divergence is defined by the following choice: f (t) := (t−21)[2] . The associated convex dual

function is :

_f_ (v) = sup _uv_ = _[v][2]_

_[∗]_ _u∈R_ _−_ [(][u][ −]2 [1)][2] 2 [+][ v.]

n o

Therefore the Pearson divergence is characterized by the following:


_q0_ 2
_PE(P_ Q) := _p0_ 1 _dµ_ sup EQφ (16)
_||_ Z  _p0_ _−_  _≥_ _φ∈Φ_ _−_ [1]2 [E][P][φ][2][ −] [E][P][φ.]

Solving the lower bound for the optimalthe Pearson divergence, we can apply a change of variables which preserves convexity of the variational optimization φ provides us an estimator of ∂f ( _p[q][0]0_ [) =] _pq00_ _[−]_ [1][. For the special case of]

problem Equation (16) and provides a more straightforward interpretation. Let the new variable be z := φ + 1 with
_z ∈Z, which in this case is nothing but the inverse function of ∂f_ . We get

sup EQφ EQz (17)
_φ∈Φ_ _−_ [1]2 [E][P][φ][2][ −] [E][P][φ][ = sup]z∈Z _−_ 2[1] [E][P][z][2][ −] 2[1]

It is now trivial to see that z is a ‘direct’ approximator of the likelihood ratio (i.e., it does not estimate a composite
map of the likelihood ratio). Therefore for simplicity, we shall employ


1
(18)
2 [E][P][φ][2][ −] [E][Q][φ]


arg min
_φ∈Φ_

to build our ‘direct’ estimator of the likelihood ratio.


Let the samples from P and Q be, respectively, xp(i) with i = 1, ..., np and xq(i) with i = 1, ..., nq. We define the
empirical estimator of the likelihood ratio _φ[ˆ]n:_


_np_

_φ(xp(i))[2]_
_−_ _n[1]q_
_i=1_

X


_nq_

_φ(xf_ (i)). (19)
_i=1_

X


_φˆn = arg min_
_φ∈Φ_


2np


**A closed form solution: Up to now we have not defined in which class of functions our approximator φ lives. As**
done in (Nguyen et al., 2010; Liu et al., 2012), we choose φ ∈ Φ where Φ is a RKHS induced by the kernel k.

We exploit the representer theorem to write a general function within Φ as:


_ntr_

_k(x, xtr(i))αi,_
_i=1_

X


_φ(x) =_


where we use ntr data which are the centers of the kernel sections used to approximate the unknown likelihood ratio
(how to choose these centers is important and determines the approximation properties of _φ[ˆ]n). For now, we do not_
specify which data should be used as centers (we can use either data from Pn or from Qn or from both or simply use
user specified locations).

Let us define the following kernel matrices: Kp := K(Xp, Xtr) ∈ R[n][p][×][n][tr], Kq := K(Xq, Xtr) ∈ R[n][q][×][n][tr], where
_Xp := {xp(i)}, Xq := {xq(i)} and Xtr := {xtr(i)}._

We therefore have:


_np_

_φ(xp(i))[2]_
_−_ _n[1]q_
_i=1_

X


_nq_

_φ(xf_ (i))
_i=1_

X


_φˆn = arg min_
_φ∈Φ_

= arg min
_α,α≥0_

= arg min
_α,α≥0_


2np

1

2np


_np_

_i=1_

X


_nf_

_i=1_

X


_ntr_

_k(xp(i), xtr(j))αj)[2]_
_−_ _n[1]f_
_j=1_

X


_ntr_

_k(xf_ (i), xtr(j))αj
_j=1_

X


1

_α[T]_ _Kp[T]_ _[K][p][α][ −]_ [1] 1[T] _Kf_ _α_
2np _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 is a standard convex optimization problem with linear constraints which
can be efficiently solved with Newton methods, nonetheless in general it does not admit any closed form solution.

We now relax the positivity constraints so that the optimal solution can be obtained in closed form. Moreover we add a
quadratic regularization term as done in (Nguyen et al., 2010) which lead us to the following regularized optimization
problem:

1
arg min _α[T]_ _Kp[T]_ _[K][p][α][ −]_ [1] 1[T] _Kf_ _α +_ _[γ]_ Φ
_α_ 2np _nf_ 2

_[∥][α][∥][2]_

whose solution is trivially given by:


_−1_
_Kp[T]_ _[K][p]_ [+][ n][p][γI][n]tr _Kf[T]_ [1][ :=][ n][p] _H_ _[−][1]Kf[T]_ [1] (20)

_nf_




_αˆ =_ _[n][p]_

_nf_


The estimator of the likelihood ratio for an arbitrary location x is given by the following:

_pq(x)_ _−1_
_φn(x) = K(x, Xtr)ˆα =_ _[n][p]_ _K(x, Xtr)_ _Kp[T]_ _[K][p]_ [+][ n][p][γI][n]tr _Kf[T]_ [1] (21)
_pp(x)_ _nf_

_[≈]_ [ˆ]  

**Remark: In the following we shall exploit RBF kernels which are defined by the length scales σ.**

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
_c_ )
CUMSUM algorithm requires to compute the likelihood ratio _[p][f][ (][y][(][t][)][|][Y][ t][−][1]_

_pp(y(t)_ _Yc[t][−][1])_ [for each time][ t][. We denote][ p][p][ as the]
_|_
normal density and pf as the abnormal one (after the anomaly has occurred).

We shall proceed to express the conditional probability p(y(t) _Y1[t][−][1]) using our predictor. In particular it is always_
_|_
possible to express the optimal (unknown) one-step ahead predictor as:

_yˆt|t−1 = F_ _[∗](Yt[t]−K+1[) :=][ E][[][y][(][t][)][ |][ Y][ t]t−K+1[]]_ (22)

which is a deterministic function given the past of the time series (whose length is K). So that the data density
distribution can be written in innovation form (based on the optimal prediction error) as:

_y(t) = F_ _[∗](Yt[t]−K+1[) +][ e][(][t][)]_ (23)

where e(t) := y(t) − _F_ _[∗](Yt[t]−K+1[)][ is, by definition, the one step ahead prediction error (or][ innovation sequence][) of]_
_y(t) given its past. We therefore have: p(y(t)_ _Yt[t]_ _K+1[) =][ p][(][e][(][t][)][ |][ Y][ t]t_ _K+1[)][. Where][ e][(][t][)][ is the optimal prediction]_
_|_ _−_ _−_
error for each time t and is therefore indipendent on each time t.

**Remark: In practice we do not know F** _[∗]_ and we use our predictor learnt from normal data as a proxy. This implies
the prediction residuals are approximately independent on normal data (the predictor can explain data well), while the
prediction residuals are, in general, correlated on abnormal data.

To summarize: under normal conditions the joint distribution of Yc[t] [can be written as:]


_p(y(i)_ _Yc[i][−][1]) =_
_|_
_i=c_

Y

_t_

_p(y(i)_ _Yc[i][−][1]) =_
_|_
_i=c_

Y


_p(Yc[t][) =]_

_p(Yc[t][) =]_


_p(e(i))_ normal conditions (24)

_i=c_

Y

_t_

_p(e(i)_ _Ec[i][−][1])_ abnormal conditions (25)
_|_
_i=c_

Y


These two conditions in turn influence the log likelihood ratio test as follows: under H0 =⇒ [Q]i[t]=c _ppfp ((ee((ii))))_ [while]

_t_ _pf (e(i)_ _Ec[i][−][1])_
under Hc = _i=c_ _pp(e|(i))_ . The main issue here is the numerator under Hc: the distribution of residuals
_⇒_

changes at each time-stamp (it is a conditional distribution) and pf (e(i) _Ec[i][−][1]) is difficult to approximate (it requires_

Q _|_

the model of the fault). In the following we show that replacing pf (e(i) _Ec[i][−][1]) with pf_ (e(i)) allows us to compute a
_|_

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 pf** (e(i) _Ec[i][−][1]) as_
_|_
independent random variables pf (e(i)) for i = 1, ..., t allows us to compute a lower bound on the log likelihood log Ω[t]c
(i.e. the cumulative sum). We denote the cumulative sum of the log likelihood ratio using independent variables as
log Ω[¯] _[t]c_ [=][ P]i[t]=c [log][ p]p[f]p[ (](e[e]([(]i[i]))[))]

**Proposition A.2. Assume a change happens at time c so that Hc is true and the following log likelihood ratio holds**
_true: log Ω[t]c_ [=][ P]i[t]=c [log][ p][f][ (]p[e]p[(][i]([)]e[|]([E]i))c[i][−][1]) _. Then it holds log Ω[t]c_ Ω[t]c[.]

_[≥]_ [log ¯]

_Proof. By simple algebra we can write:_
_pf_ (e(i) | Ec[i][−][1]) = _[p][f]_ [(][e][(][i][)][ |][ E]c[i][−][1]) _pf_ (e(i)) _i_

_pp(e(i))_ _pf_ (e(i)) _pp(e(i))_ _∀_

Now recall the cumulative sum of the log-likelihood ratios taken under the current data generating mechanism pf (E1[t][)]
provides an estimate of the expected value of the log-likelihood ratio. Due to the correlated nature of data E1[t] [the]
samples are drawn from a multidimensional distribution of dimension t (a sample from this distribution is an entire
trajectory from c to t).

We now take the expectation of previous formula w.r.t. the ‘true’ distribution pf (E1[t][)][:]


_pf_ (e(i) _Ec[i][−][1])_
_|_

_pp(e(i))_

_pf_ (e(i) _Ec[i][−][1])_
_|_ + Epf (Ect[)][ log]

_pf_ (e(i))


Epf (Ect[)][Ω]c[t] [=][ E]pf (Ec[t][)][ log]

= Epf (Ect[)][ log]


_i=c_

_t_

_i=c_

Y


_pf_ (e(i))
_pp(e(i))_

_t_

_pp(e(i))_
_i=c_

Y


_i=c_


= MI _pf_ (Ec[t][);]



_pf_ (e(i)) + KL
_i=c_

Y 


_pf_ (e(i))
_i=c_

Y


_≥_ _KL_


_pf_ (e(i))
_i=c_

Y


_pp(e(i))_
_i=c_

Y


where we used the fact the mutual information is always non negative.

**How to approximate pre and post fault distributions: Both pp and pf are not known and their likelihood ratio**
need to be estimated from available data. From Section 5.1.1 we know how to approximate the likelihood ratio
given two set of data without estimating the densities. In our anomaly detection setup we define these two sets as:
_Ep := Et[t]−[−]n[n]p[f]−nf +1_ [and][ E][f][ :=][ E]t[t]−nf +1[. So that given current time][ t][ we look back at a window of length][ n][p] [+][ n][f] [.]
The underlying assumption is that under Hc normal data are present in Ep and abnormal ones in Ef . We estimate the
likelihood ratio _[p]p[f]p[ (](e[e]([(]t[t]))[))]_ [at each time][ t][ by assuming both][ E][p][ and][ E][f][ data are independent (see Proposition A.2) and]

cumulate their log as t increases.

**How do np and nf affect our detector? The choice of the windows length (np and nf** ) is fundamental and highly
influences the likelihood estimator. Using small windows makes the detector highly sensible to both point and sequential outliers, while larger windows are better suited to estimate only sequential outliers. We now assume np = nf and
study how small and large values affect the behaviour of our detector in simple working conditions.

In Figure 5 and Figure 6 we compute the cumulative sum of log likelihood ratios estimated from data on equally sized
windows. Intuitively any local minimum after a ‘large’ (depending on the threshold τ ) increase of the cumulative sum
is a candidate abnormal point.

In Figure 7 and Figure 8 we compare the cumulative sum of estimated likelihood ratios on data in which both sequential
and point outliers are present. In particular we highlight that large window sizes np and nf are usually not able to
capture point anomalies Figure 7 while using small window sizes allow to detect both (at the expenses of a more
sensitive detector) Figure 8.

21


-----

40 50 60 70 80 90 100

|Col1|Col2|
|---|---|
|5 4 3 2 1 0|cum Dat poin|
|||


cumsum
Data
point outlier

Time


3.5

3.0

2.5

2.0

1.5

Values 1.0

0.5

cumsum

0.0 Data

change point

0.5

40 50 60 70 80 90 100

Time


Figure 5: Change point: Cumulative sum (blue) obtained with our method in a synthetic example. We use
the cumulative sum of estimated likelihood ratios on data
in which a change point is present at t = 60. We use
_np = nf = 20 and kernel length scale=0.2_

Time


Figure 6: Point anomaly: Cumulative sum (blue) obtained with our method in a synthetic example. We use
the cumulative sum of estimated likelihood ratios on data
in which a point outlier is present at t = 60. We use
_np = nf = 2 and kernel length scale=2._

Values 2

|Col1|sum|Col3|
|---|---|---|
|Data anom 6 4 2 0|alous points||
||||


cumsum
Data
anomalous points

50 75 100 125 150 175 200


Time


4

2

0

Values

2

cumsum
Data

4

anomalous points

50 75 100 125 150 175 200


Figure 7: Large np and nf : Cumulative sum (blue)
obtained with our method in a synthetic example. We
use the cumulative sum of estimated likelihood ratios
on data which contain both change points (t = 60 and
_t = 200) and point outliers (t = 80 and t = 160). We_
use np = nf = 20 and kernel length scale=1.

A.6 DATASETS


Figure 8: Small np and nf : Cumulative sum (blue)
obtained with our method in a synthetic example. We
use the cumulative sum of estimated likelihood ratios
on data which contain both change points (t = 60 and
_t = 200) and point outliers (t = 80 and t = 160). We_
use np = nf = 3 and kernel length scale=5.


A.6.1 YAHOO DATASET

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


22


-----

A.6.2 NAB DATASET

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


A.6.3 CO2 DATASET

We test the prediction and interpretability capabilities of our model on the CO2 dataset from kaggle [2]. The main goal
here is to predict both trend and periodicity of CO2 emission rates on different years. Note this is not an Anomaly
detection task.
Table 3: Datasets summaries. We report some properties of the datasets used (see (Geiger et al., 2020) for mode
details).


**Yahoo** **NAB** **Kaggle**


A1 A2 A3 A4 Traffic Tweets CO2

# signals 67 100 100 100 7 10 9
# anomalies 178 200 939 835 14 33
point 68 33 935 833 0 0
sequential 110 167 4 2 14 33
# anomalous 1669 466 943 837 1560 15651

**Properties**

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

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 to one. As done in (Braei & Wagner, 2020) we note standardization is not equal to normalization,
where data are forced to belong to the interval (0, 1). Normalization is more sensitive to outliers, thus it would be
inappropriate to normalize our datasets, which contain outliers.

We do not apply any deseasonalizing or detrending pre-processing.

**Data splitting: We split each dataset into training and test sets preserving time ordering, so that the first data are used**
as train set and the following ones are used as test set. The data used to validate the model during optimization are last
10% of the training dataset. Depending on the experiment, we choose a different percentage in splitting train and test.
When comparing with (Braei & Wagner, 2020) we used 30% as training data, while when comparing to (Munir et al.,
2019) we use 40%. Such a choice is dictated by the particular (non uniform) experimental setup reported in (Braei &

2https://www.kaggle.com/txtrouble/carbon-emissions
3https://catalog.ldc.upenn.edu/LDC2008T19

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 metrics: 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 Y1[N] and the one-step ahead predictions _Y[ˆ]1[N]_ the RMSE is

defined as: _N1_ _Ni=1_ _y(i)_ .

_[∥][y][(][i][)][ −]_ [ˆ] _∥[2]_

q

As done in (Braei & Wagner, 2020) we compare different anomaly detection methods taking into account severalP
metrics. We use F1-Score which is defined as the armonic mean of Precision and Recall (see (Braei & Wagner, 2020;
Munir et al., 2019)) and another metric that is often used is receiver operating characteristic curve, ROC-Curve, and
its associated metric area under the curve (AUC). The AUC is defined as the area under the ROC-Curve. This metric
is particularly useful in our anomaly detection setting since it describes with an unique number true positive rate and
_false positive rate on different threshold values. We now follow (Braei & Wagner, 2020) to describe how AUC is_
computed. Let the true positive rate and false positive rate be defined, respectively, as: TPR = _[T P]P_ [and][ FPR][ =][ F P]N [,]

where TP stands for true positive, P for positive, FP 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 x and y axes are, respectively: FPR and TPR.
The resulting curve starts at the origin and ends in the point (1,1). The AUC is the area under this curve. In anomaly
detection, the AUC expresses the probability that the measured algorithm assigns a random anomalous point in the
time series a higher anomaly score than a random normal point.

**Hardware: We conduct our experiments on the following hardware setup:**

-  Processor: Intel(R) Core(TM) i9-10980XE CPU @ 3.00GHz

-  RAM: 128 Gb

-  GPU: Nvidia TITAN V 12Gb and Nvidia TITAN RTX 24Gb

**Hyper-parameters: All the experiments we carried out are uniform on the optimization hyper-parameters. In 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: l0 = 10, l1 = 100, l2 = 200

-  linear filters kernel lengths (N0, N1, N2): half predictor’s memory

In all experiments we either use a TCN composed of 3 hidden layers with 300 nodes per layer or a TCN with 8 layers
and 32 nodes per layer. Moreover we chose N3 = 5 (TCN kernels’ lengths) and relu activation functions (Bai et al.,
2018).

**Comparison with SOTA methods: We tested our model against other SOTA methods (Table 1) in a comparable**
experimental setup. In particular, we chose comparable window lengths and architecture sizes (same order of magnitude of the number of parameters) to make the comparison as fair as possible. For the hyper-parameters details
of any SOTA method we used we refer the relative cited reference. We point out that while the window length is a
critical hyper-parameter for the accuracy of many methods, our architecture is robust w.r.t. choice of window length:
thanks to our fading regularization, the user is required only to choose a window length larger than the optimal one
and then our automatic complexity selection is guaranteed to find the optimal model complexity given the available
data Section 4.1.

**Anomaly scores: When computing the F-score we use the predictions of the CUMSUM detector which we collect as**
a binary vector whose length is the same as the number of available data. Ones are associated to the presence of an
anomalous time instants while zeros are associated to normality.

When computing the AUC we need to consider a continuous anomaly score, therefore the zero-one encoded vector
from the CUMSUM is not usable. We compute the anomaly scores for each time instant as the estimated likelihood
ratios. Since we write the likelihood ratio as _[p]p[f]p_ [, it is large when data does not come from][ p][p][ (which we consider the]

reference distribution).

24


-----

Train Error

|Standard TCN STRIC No Fading STRIC|Col2|Standard TCN STRIC No Fading|Col4|Standard TCN STRIC No Fading STRIC|Col6|
|---|---|---|---|---|---|



20 40Memory60 80 100

Standard TCN
STRIC No Fading
STRIC

Train Error


Test Error


Generalization Gap

Standard TCN
STRIC No Fading
STRIC

20 40Memory60 80 100

Generalization Gap


1.2

1.0

0.8

0.6

0.4

0.2

0.0

2.5

2.0

1.5

1.0

0.5

0.0


Standard TCN
STRIC No Fading

20 40Memory60 80 100

STRIC

Test Error

|Standard TCN STRIC No Fading STRIC|Standard TCN STRIC No Fading STRIC|Col3|Standard TCN STRIC No Fading STRIC|Standard TCN STRIC No Fading STRIC|Col6|Col7|Standard TCN STRIC No Fading STRIC|Standard TCN STRIC No Fading STRIC|Col10|
|---|---|---|---|---|---|---|---|---|---|


20 40Memory60 80 100

Standard TCN Standard TCN
STRIC No Fading STRIC No Fading
STRIC STRIC

Train Error

|Standard TCN STRIC No Fading STRIC|Standard TCN STRIC No Fading STRIC|Col3|Standard TCN STRIC No Fading STRIC|Standard TCN STRIC No Fading STRIC|Col6|Standard TCN STRIC No Fading STRIC|Standard TCN STRIC No Fading STRIC|Col9|
|---|---|---|---|---|---|---|---|---|



20 40Memory60 80 100

Standard TCN
STRIC No Fading
STRIC

Train Error


20 40Memory60 80 100

Test Error

Standard TCN
STRIC No Fading
STRIC

20 40Memory60 80 100

Test Error


20 40Memory60 80 100

Generalization Gap

Standard TCN
STRIC No Fading
STRIC

20 40Memory60 80 100

Generalization Gap


0.8

0.6

0.4

0.2

0.0

1.2

1.0

0.8

0.6

0.4

0.2

0.0

|Standard TCN STRIC No Fading STRIC|Standard TCN STRIC No Fading STRIC|Col3|Standard TCN STRIC No Fading STRIC|Standard TCN STRIC No Fading STRIC|Col6|Standard TCN STRIC No Fading STRIC|Standard TCN STRIC No Fading STRIC|Col9|
|---|---|---|---|---|---|---|---|---|


20 40Memory60 80 100

Standard TCN Standard TCN
STRIC No Fading STRIC No Fading
STRIC STRIC

Train Error


20 40Memory60 80 100

Test Error


20 40Memory60 80 100

Generalization Gap


0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0.0

1.2

1.0

0.8

0.6

0.4

0.2

0.0

|Standard TCN STRIC No Fading STRIC|Col2|Col3|Col4|Standard TCN STRIC No Fading STRIC|Col6|
|---|---|---|---|---|---|


Standard TCN
STRIC No Fading
STRIC

Standard TCN
STRIC No Fading
STRIC


|Memory Train Error|Col2|Col3|Memory Test Error|Col5|Memory Generalization Gap|Col7|
|---|---|---|---|---|---|---|
|Standard TCN STRIC No Fading STRIC|Standard TCN STRIC No Fading STRIC||Standard TCN STRIC No Fading||Standard TCN STRIC No Fading STRIC||
||||Standard TCN STRIC No Fading||||


20 40Memory60 80 100

STRIC

Test Error


20 40Memory60 80 100

Standard TCN
STRIC No Fading
STRIC

Train Error

|Standard TCN STRIC No Fading STRIC|Col2|Col3|Col4|Standard TCN STRIC No Fading STRIC|Col6|
|---|---|---|---|---|---|



20 40Memory60 80 100


20 40Memory60 80 100

Generalization Gap

Standard TCN
STRIC No Fading
STRIC


1.2

1.0

0.8

0.6

0.4

0.2

0.0


20 40Memory60 80 100


20 40Memory60 80 100


Standard TCN
STRIC No Fading
STRIC

Standard TCN
STRIC No Fading
STRIC


Figure 9: Ablation studies on different datasets: Effects of interpretable blocks and fading regularization on 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 a 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) as soon as 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 metrics based on the predictor’s RMSE (training, test and generalization gap) as a
function of the memory of the predictor. We test our fading regularization on a variety of different datasets. In all
situations fading regularization helps improving test generalization and preserving the generalization gap (by keeping
it constant) as the model complexity increases. All plots show confidence intervals around mean values evaluated on
10 different random seeds.

In Table 4 we extend the results we show in the main paper by adding uncertainties (measured by standard deviations
on 10 different random seeds) to the values of train and test RMSE on different ablations of STRIC. Despite the high
variability across different datasets STRIC achieves the most consistent results (smaller standard deviations both on
training and testing).

Finally, in Table 5 we show the effects on different choices of the predictor’s memory npred and length of the anomaly
detectors windows ndet on the detection performance of STRIC. Note both F-score and AUC are highly sensible to the
choice of ndet: the best results are achieve for small windows. On the other hand when ndet is large the performance
drops. This is due to the type of anomalies present in the Yahoo benchmark: most of the them can be considered to be
point anomalies. In fact, as we showed in Appendix A.5, our detector is less sensible to point anomalies when a large
window ndet is chosen.

In Table 5 we also report the reconstruction error of the optimal predictor given it’s memory npred. Note small memory
in the predictor introduce modelling bias (higher training error) while a large memory does not (thanks to fading
regularization). As we observed in Appendix A.5 better predictive models provide the detection module with more
discriminative residuals: the downstream detection module achieves better F-scores and AUC.
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 hyperparameters). The effect of adding a linear interpretable model before a TCN improves generalization error. Fading
regularization has a beneficial effect in controlling the complexity of the TCN model and reducing the generalization
gap.


**TCN** **TCN + Linear** **TCN + Fading** **STRIC pred**


Train Test Train Test Train Test Train Test

Yahoo A1Yahoo A2Yahoo A3Yahoo A4CO2 DatasetNAB TrafficNAB Tweets **0.100.110.130.150.140.030.18 ± ± ± ± ± ± ± 0.06 0.02 0.01 0.01 0.02 0.01 0.05** 0.920.820.430.610.621.061.02 ± ± ± ± ± ± ± 0.06 0.02 0.01 0.01 0.02 0.02 0.05 0.130.160.190.150.040.200.10 ± ± ± ± ± ± ± 0.03 0.01 0.01 0.01 0.02 0.01 0.05 0.880.350.350.451.000.980.22 ± ± ± ± ± ± ± 0.03 0.02 0.01 0.01 0.02 0.02 0.05 0.440.200.150.170.180.620.47 ± ± ± ± ± ± ± 0.03 0.01 0.01 0.01 0.03 0.01 0.02 0.920.710.400.550.610.930.83 ± ± ± ± ± ± ± 0.03 0.01 0.01 0.01 0.03 0.01 0.02 0.430.140.190.230.330.630.70 ± ± ± ± ± ± ± 0.02 0.01 0.01 0.01 0.01 0.02 0.01 **0.620.300.220.240.410.740.77 ± ± ± ± ± ± ± 0.02 0.01 0.01 0.01 0.01 0.02 0.01**


Table 5: Sensitivity of STRIC to hyper-parameters: We compare STRIC on different anomaly detection benchmarks
datasets using different hyper-parameters: memory of the predictor npred and length of anomaly detector windows
_np = nf = ndet._

**Yahoo A1** **Yahoo A2** **Yahoo A3** **Yahoo A4**


F1 AUC F1 AUC F1 AUC F1 AUC

_npred = 10, ndet = 2_ 0.45 0.89 0.63 0.99 0.87 0.99 0.64 0.89
_npred = 100, ndet = 2_ **0.48** **0.9308** **0.98** **0.9999** **0.99** **0.9999** **0.68** **0.9348**
_npred = 10, ndet = 20_ 0.10 0.58 0.63 0.99 0.47 0.83 0.37 0.72
_npred = 100, ndet = 20_ 0.10 0.55 **0.98** **0.9999** 0.49 0.86 0.35 0.76

**Yahoo A1** **Yahoo A2** **Yahoo A3** **Yahoo A4**


Train Test Train Test Train Test Train Test


_npred = 10_ 0.44 0.62 0.16 0.31 0.22 0.23 0.25 0.26

**Models** _npred = 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 errors of a standard TCN model against our STRIC on the A3 Yahoo dataset. We train both
models using the same optimization hyper-parameters (as described in previous section). Note a plain TCN does not
necessarily capture the trend component in the test set.


Plain TCN

200 400 600 800 1000 1200 1400 1600


STRIC

200 400 600 800 1000 1200 1400 1600


train data
test data
train predictions
test predictions


200 400 600 800 1000 1200 1400 1600

Timestamps


200 400 600 800 1000 1200 1400 1600

Timestamps


Figure 10: We compare an off-the-shelf TCN against STRIC (time series predictor) on the Yahoo dataset A3 Benchmark. 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. In particular the generalization gap (measured using squared reconstruction error)
for the two models is: Standard TCN 0.3735 and STRIC 0.0135.


Plain TCN

500 600 700 800 900 1000 1100 1200

4 4 train data

test data

3 3 train predictions

test predictions

2 2

Value 1 1

0 0

1 1

2 2

3 3

Timestamps


STRIC

4 train data

test data

3 train predictions

test predictions

2

1

0

1

2

3

500 600 700 800 900 1000 1100 1200

Timestamps


Figure 11: Zoom on the second row of panels in Figure 10. We show the interface between train and test data both on
a plain TCN and on our STRIC predictor. A plain TCN overfits w.r.t. STRIC also when not trend is present.

A.8 STRIC VS 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 1 and Table 6. Table 6 is used to highlight the relative performance
of STRIC when the peformance are nearly saturated (e.g. Yahoo A2 and A3): we report the relative performance
indeces against TRAID for all the models we tested in Table 6. Both for F1 and AUC we report the following for each
comparing SOTA method: [AUC]1[method]AUC[−][AUC]STRIC[STRIC] 100 (similarly for F1).

_−_ _·_

To begin with, STRIC outperforms ‘traditional’ methods (LOF and One-class SVM) which are considered as baselines
models for comparing time series anomaly detectors.

27


-----

Table 6: Comparison with SOTA anomaly detectors: We compare STRIC with other anomaly detection methods
on 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 et al., 2020), LSTM (Braei & Wagner, 2020; Munir et al., 2019),
Wavenet (Braei & Wagner, 2020), Yahoo EGADS (Munir et al., 2019), 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., 2017) and DeepAnT (Munir et al., 2019) . STRIC outperforms most of the other
methods based on statistical models and based on DNNs. Same as Table 1, here we report the relative improvements
w.r.t. STRIC (the higher the better).


**Relative** **F1-score**
**improvement** **over**
**STRIC in %**


**Yahoo A1** **Yahoo A2** **Yahoo A3** **Yahoo A4** **NAB Tweets** **NAB Traffic**


ARIMA -20 -88 -42 **6** -33 -37
LSTM -7 -33 -60 -  21
Yahoo EGADS -1 -95 -78 -54
OmniAnomaly -1 -60 -45 -11 -6 -10
Twitter AD **0** -98 -85 -53
TanoGAN -11 -85 -73 -13 -36 -44
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-**
**provement** **over**
**STRIC in %**


**Yahoo A1** **Yahoo A2** **Yahoo A3** **Yahoo A4** **NAB Tweets** **NAB Traffic**


MA -47 -98 -98 **379**
ARIMA -45 -99 -99 124
LOF -28 -99 -99 -81 -32 -44
Wavenet -60 -99 -99 -84
LSTM -63 -99 -99 -84
GOAD -37 -99 -99 -51 -19 -12
DeepAnT -32 -99 -99 -53 -23 -13

**Models** STRIC (ours) **0** **0** **0** 0 **0** **0**

**Comparison with other Deep Learning based methods: STRIC outperforms most of the SOTA Deep Learning**
based methods reported in Table 1: TadGAN, TAnoGAN, DeepAnT and DeepAR (the last one is a SOTA time series predictor). Note the relative improvement of STRIC is higher on the Yahoo dataset where statistical models
outperforms deep learning based ones. We believe this is due to both fading regularization and the seasonal-trend
decomposition performed by STRIC.

Despite the general applicability of GOAD (Bergman & Hoshen, 2020) this method has not been designed to handle
time series, but images and tabular data. “Geometric” transformations which have been considered in GOAD and
actually have inspired it (rotations, reflections, translations) might not be straightforwardly applied to time series.
Nevertheless, while we have not been able to find in the literature any direct and principled extension of this work
to the time series domain, we have implemented and compared against (Bergman & Hoshen, 2020) by extending the
main design ideas of GOAD to time-series. So that we applied their method on lagged windows extracted from time
series (exploiting the same architectures proposed for tabular data case with some minor modifications). We report the
results we obtained by running the GOAD’s official code on all our benchmark datasets. Overall, STRIC performs (on
average) 70% better than GOAD on the Yahoo dataset and 15% better on the NAB dataset.


A.8.1 DETAILS ON THE NYT EXPERIMENT

Figure 4 shows the normalized anomaly score computed by STRIC on the NYT dataset, following the setup described
in Appendix A.6.4. Some additional insights can be gained by zooming in around some of the detected change-points.
In Figure 12 (left), we see that the anomaly score (blue line) rapidly increases immediately after the 9/11 attack and
reaches its peak some days later. Such delay is inherently tied to our choice of time scale, that privileges the detection
of prolonged anomalies as opposed to single-day anomalies (which are not meaningful due to the high variability of


28


-----

Figure 12: A closer look at some of the change-points detected by STRIC. Left: Normalized anomaly score (blue
line) and normalized frequency of the “Terrorism” descriptor (orange line) around the 9/11 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
change-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 on any human annotations. However,
they can help interpreting the change-points found by STRIC.

In Figure 12 (right), we can observe that the anomaly score (blue line) is higher in the months around the 2004 U.S.
election and immediately after the inauguration day. However, the highest values for the anomaly score occur around
the end of 2004, shortly after the Indian Ocean tsunami. Indeed, this is reflected by an abrupt increase of the frequency
of descriptors like “Earthquakes” (orange line) and “Tsunami”.

We note this experiment is qualitative and unfortunately we are not aware of any ground truth or metrics (e.g., in
Rayana & Akoglu (2015) a similar qualitative result has been reported on the NYT dataset). We therefore tested
STRIC against a simple baseline which uses PCA on BERT features and a threshold to detect anomalies. Despite being
a simple baseline, this method prooved to be highly applied in practice due to its simplicity (Bl´azquez-Garc´ıa et al.,
2020). The PCA + threshold baseline is able to pick up some events (2000 election, 9/11 attack, housing bubble) but
is otherwise more noisy than STRIC’s anomaly score. This is likely due to the lack of a modeling of seasonal/periodic
components. For instance, the anomaly score of the simple baseline contains many false alarms which are related to
normal weekly periodicity that is not easily modeled by the baseline. This does not affect STRIC’s predictions, since
normal weekly periodicity is directly modeled and identified as normal behaviour.

29


-----

