Anomaly Detection of Time Series with Smoothness-Inducing Sequential Variational Auto-Encoder

In plain words

Sensor data and server logs often contain faulty readings and real anomalies. A model trained to reproduce everything it sees also learns those bad points, so it stops flagging them. The usual way to score an anomaly is to ask how badly the model reconstructs a point, and that score breaks when the normal noise level changes over time.

This paper's model is a recurrent version of a variational autoencoder. For every time step it learns a range of likely values and how much noise to expect there, so the noise level can change over time. It adds one extra rule: the predicted distribution should change smoothly from one time step to the next. That rule stops the model from copying single spikes. Points the model finds unlikely are flagged. The authors compare two scoring methods, a probability-based one and a plain error-based one.

TL;DR

• SISVAE is a sequential VAE for multi-dimensional, unsupervised time-series anomaly detection §1 §4.
• Each time step gets its own predicted mean and variance; a KL-based smoothness term links consecutive predictions §4.2 eq. 18 eq. 19.
• Anomalies are scored by Monte Carlo reconstruction probability (preferred) or absolute reconstruction error §4.4 alg. 2.
• Author-reported: highest AUPRC on four of five real datasets against the baselines in Table IV tab. 4.
• Eq. 20 as printed has a sign that conflicts with the text's "penalty" wording, and the code is not public; check the sign before implementing (see Pitfalls).

Problem

Multi-dimensional time series are often contaminated with anomalies, and the goal is to recover the true signal and detect anomalous points without labels §3.1 eq. 1. Standard smoothness-based models assume constant noise, which fails for non-stationary series §2. Deterministic regularisers do not carry over to probabilistic VAE models, because the observations are random variables there §4.2. Constant-variance models are also restricted to simple linear forms with conjugate distributions §2.

The paper assumes unsupervised training with no anomaly labels and point-level detection §1 §5.1. Anomaly-free training data are not assumed.

Method

Background. A smoothness prior on a signal is motivated by a Gaussian observation model with a penalty on high-order differences §3.2 eq. 2 eq. 3 eq. 4 eq. 5. A VAE maximises the evidence lower bound (ELBO), with the reparameterisation trick giving the SGVB estimator §3.3 eq. 7 eq. 8:

\mathcal{L}(\theta,\phi;\mathbf{x}) = -D_{KL}\big(q_\phi(\mathbf{z}\mid\mathbf{x}) \,\|\, p_\theta(\mathbf{z})\big) + \mathbb{E}_{q_\phi(\mathbf{z}\mid\mathbf{x})}\big[\log p_\theta(\mathbf{x}\mid\mathbf{z})\big] \quad [\text{eq. 7}]

Step 1: sequential VAE density model. The joint factorises over time, with the latent prior and likelihood conditioned on past states §4.1.1 eq. 9:

p_\theta(\mathbf{x}_{\le T},\mathbf{z}_{\le T}) = \prod_{t=1}^{T} p_\theta(\mathbf{x}_t \mid \mathbf{z}_{\le t}, \mathbf{x}_{<t})\, p_\theta(\mathbf{z}_t \mid \mathbf{x}_{<t}, \mathbf{z}_{<t}) \quad [\text{eq. 9}]

• The decoder outputs a per-time-step diagonal Gaussian for x_t §4.1.1 eq. 10.
• The latent prior at step t is a Gaussian whose mean and variance come from h_{t-1} §4.1.1 eq. 11.
• The recurrent state is a GRU over the encoded latent §4.1.1 eq. 12 eq. 13.

[\boldsymbol{\mu}_{x,t}, \boldsymbol{\sigma}_{x,t}] = \varphi^{\mathrm{dec}}_\theta\big(\varphi^{\mathbf{z}}_\theta(\mathbf{z}_t), \mathbf{h}_{t-1}\big), \quad p_\theta(\hat{\mathbf{x}}\mid\mathbf{z}_t) = \mathcal{N}\big(\boldsymbol{\mu}_{x,t}, \mathrm{diag}(\boldsymbol{\sigma}^2_{x,t})\big) \quad [\text{eq. 10}]

\mathbf{z}_t \sim \mathcal{N}\big(\boldsymbol{\mu}_{0,t}, \mathrm{diag}(\boldsymbol{\sigma}^2_{0,t})\big), \quad [\boldsymbol{\mu}_{0,t}, \boldsymbol{\sigma}_{0,t}] = \varphi^{\mathrm{prior}}_\theta(\mathbf{h}_{t-1}) \quad [\text{eq. 11}]

\mathbf{h}_t = f_\theta\big(\varphi^{\mathbf{z}}_\theta(\mathbf{z}_t), \mathbf{h}_{t-1}\big) \quad [\text{eq. 12}]

The GRU is written as follows in the paper §4.1.1 eq. 13. It is reproduced as printed, and the typos it contains are listed under Pitfalls:

\begin{aligned} \mathbf{y}_t &= \varphi^{\mathbf{z}}_\theta(\mathbf{z}_t) \\ \mathbf{r}_t &= \sigma(\mathbf{W}_r \mathbf{y}_t + \mathbf{U}_r \mathbf{h}_{t-1} + \mathbf{b}_r) \\ \mathbf{s}_t &= \sigma(\mathbf{W}_s \mathbf{y}_t + \mathbf{U}_s \mathbf{h}_{t-1} + \mathbf{b}_r) \\ \tilde{\mathbf{h}}_t &= \tanh(\mathbf{W}_h + \mathbf{r}_t (\mathbf{U}_h \mathbf{h}_{t-1}) + \mathbf{b}_h) \\ \mathbf{h}_t &= (1 - \mathbf{s}_t)\tilde{\mathbf{h}}_t + \mathbf{s}_t \mathbf{h}_{t-1} \end{aligned} \quad [\text{eq. 13}]

Step 2: inference model. The posterior factorises over time, and each latent depends on the current observation and h_{t-1} §4.1.2 eq. 15 eq. 16:

q_\phi(\mathbf{z}_{\le T}\mid\mathbf{x}_{\le T}) = \prod_{t=1}^{T} q_\phi(\mathbf{z}_t \mid \mathbf{x}_{\le t}, \mathbf{z}_{<t}) \quad [\text{eq. 15}]

[\boldsymbol{\mu}_{z,t}, \boldsymbol{\sigma}_{z,t}] = \varphi^{\mathrm{enc}}_\phi\big(\varphi^{\mathbf{x}}_\theta(\mathbf{x}_t), \mathbf{h}_{t-1}\big) \quad [\text{eq. 16}]

Training uses AEVB §4.1.2 eq. 14.

Step 3: smoothness-regularised objective. The sequential ELBO is written as a sum over time of an inference term and a reconstruction term §4.2 eq. 17. The paper adds a smoothness term, the KL divergence between the reconstruction distributions at consecutive steps §4.2 eq. 18 eq. 19. The per-step Gaussian KL for eq. 19 is:

D_{KL}\big(\mathcal{N}(\mu_{t-1},\sigma_{t-1}) \,\|\, \mathcal{N}(\mu_t,\sigma_t)\big) = \log\frac{\sigma_t}{\sigma_{t-1}} + \frac{\sigma_{t-1}^2 + (\mu_{t-1}-\mu_t)^2}{2\sigma_t^2} - \frac{1}{2} \quad [\text{eq. 19}]

The full learning objective, as printed, is §4.2 eq. 20:

\mathcal{L}_{\mathrm{SISVAE}} = \sum_{t=1}^{T} \Big\{ -D_{KL}\big(q_\phi(\mathbf{z}_t\mid\mathbf{x}_{\le t},\mathbf{z}_{<t}) \,\|\, p_\theta(\mathbf{z}_t\mid\mathbf{x}_{<t},\mathbf{z}_{<t})\big) + \log p_\theta(\mathbf{x}_t\mid\mathbf{z}_{\le t},\mathbf{x}_{<t}) + \lambda\, D_{KL}\big(p_\theta(\hat{\mathbf{x}}_{t-1}\mid\cdot) \,\|\, p_\theta(\hat{\mathbf{x}}_t\mid\cdot)\big) \Big\} \quad [\text{eq. 20}]

The paper sets the smoothness term for t = 1 to use the same distribution on both sides §4.2. The smoothness term is not given in Algorithm 1 alg. 1.

Step 4: training on chunks. Long series are cut into windows of length W with step s, one window per training item §4.3.1 eq. 21 eq. 22. The paper gives W = 120 and does not give s §5.1. Algorithm 1 processes each chunk time step by time step, with a zero initial state alg. 1.

# One SISVAE training loss on a mini-batch of chunks  [alg. 1] [eq. 20]
# Input:  chunks (P, M, W) cut by sliding windows  [§4.3.1]
#         lam: smoothness weight (0.5 in the paper)  [§5.1]
# Output: scalar loss to minimise

def training_loss(chunks, lam):
    # h_0 is zero  [alg. 1]
    h_prev = zeros(P, h_dim)
    # Accumulates the per-time-step terms of eq. 20  [eq. 20]
    total = 0
    prev_mu_x = None
    prev_sigma_x = None
    for t in range(1, W + 1):
        # Observation at step t, shape (P, M)  [alg. 1]
        x_t = chunks[:, :, t]
        # Inference: mean and std of z_t from x_t and h_{t-1}  [eq. 16]
        mu_z, sigma_z = encoder(features_x(x_t), h_prev)
        # Reparameterised sample of z_t  [alg. 1]
        z_t = mu_z + sigma_z * normal_noise(P, latent_dim)
        # Prior of z_t from h_{t-1}  [eq. 11]
        mu_0, sigma_0 = prior_net(h_prev)
        # Decoder: mean and std of the reconstruction of x_t  [eq. 10]
        mu_x, sigma_x = decoder(features_z(z_t), h_prev)
        # Inference term: KL(q(z_t) || p(z_t)), first term of eq. 20
        kl_latent = gaussian_kl(mu_z, sigma_z, mu_0, sigma_0)
        # Reconstruction term: log p(x_t | z_t), second term of eq. 20
        recon = gaussian_log_lik(x_t, mu_x, sigma_x)
        # Smoothness term: KL between consecutive reconstructions  [eq. 19]
        # At t = 1 the paper uses the same distribution twice  [§4.2]
        if t == 1:
            smooth = 0
        else:
            smooth = gaussian_kl(prev_mu_x, prev_sigma_x, mu_x, sigma_x)
        # Eq. 20 as printed: the smoothness term has a plus sign.
        # See Pitfalls: with the bound maximised, this sign rewards
        # non-smooth outputs.
        step_term = -kl_latent + recon + lam * smooth
        total = total + step_term
        # Recurrence: h_t from z_t and h_{t-1}  [eq. 12] [eq. 13]
        h_prev = gru(features_z(z_t), h_prev)
        prev_mu_x = mu_x
        prev_sigma_x = sigma_x
    # Eq. 14 maximises the bound, so return its negative.
    # Batch reduction (mean or sum) is not stated.
    return -mean(total)

Step 5: anomaly scoring. Two scores are described §4.4:

• Absolute reconstruction error: e_{m,t} = |x_{m,t} − μ_{m,t}| eq. 23. The paper notes this assumes constant variance §4.4 fig. 3.
• Reconstruction probability: the density p(x) is the expectation of p(x|z) under the encoder eq. 24. Because z is sequentially dependent, the paper uses Monte Carlo with L draws §4.4 alg. 2. The paper names this Sequential Monte Carlo, but Algorithm 2 has no weights or resampling step (see Pitfalls).

# Monte Carlo anomaly score as written in the paper  [alg. 2] [§4.4]
# Input:  trained encoder and decoder, chunk (M, W), L draws
# Output: score (M, W); anomaly if score > alpha  [alg. 2]

def anomaly_score(chunk, num_draws):
    # Score starts at zero for every variable and time step  [alg. 2]
    score = zeros(M, W)
    # Repeat the full pass L times, each with fresh noise  [alg. 2]
    for draw in range(num_draws):
        # Algorithm 2 does not state h_0; zero is assumed as in Alg. 1
        h_prev = zeros(h_dim)
        for t in range(1, W + 1):
            x_t = chunk[:, t]
            # Encode x_t with h_{t-1} and draw z_t  [alg. 2]
            mu_z, sigma_z = encoder(features_x(x_t), h_prev)
            z_t = mu_z + sigma_z * normal_noise(latent_dim)
            # Decode the reconstruction distribution of x_t  [alg. 2]
            mu_x, sigma_x = decoder(features_z(z_t), h_prev)
            # Add -log p(x_t) / L to this step's score  [alg. 2]
            nll = -gaussian_log_lik(x_t, mu_x, sigma_x)
            score[:, t] = score[:, t] + nll / num_draws
            # Advance the recurrent state  [alg. 2]
            h_prev = gru(features_z(z_t), h_prev)
    # Threshold: anomaly where score > alpha  [alg. 2]
    return score

Real datasets are sliced into non-overlapping chunks for scoring, and each chunk is reconstructed separately §4.4.

Preprocessing per dataset. Each sequence is standardised to a standard normal distribution before training §5.3. SISVAE-p takes the multivariate matrix through sliding windows §5.1. Donut, which is univariate by design, was run on flattened sequences §5.1.

Hyperparameters

Results

Author-reported. Setup: five real datasets (Yahoo S5 A1–A4 and μPMU), unsupervised training with no anomaly labels, point-level metrics AUROC, AUPRC and best F1 over thresholds §5.1 §5.3. No variance across runs is reported in the table.

Datasets tab. 3:

Real-world performance tab. 4. Values are as printed; the columns are reordered from the paper's layout, with this paper's models first. Bold marks the best value per dataset and metric in the paper's table. The table caption says SISVAE-e is left out of the paper's text comparisons, but it is in the table.

Headline (author-reported): SISVAE-p has the best AUPRC on A1, A2, A3 and A4, and the best F1 on A1–A4 tab. 4. On μPMU, Donut is best for AUPRC (0.312 against 0.310) and F1 (0.389 against 0.367) tab. 4. Best AUROC is held by SISVAE-p on A3 and μPMU, by Donut on A2 and A4, and by VAE-s on A1 tab. 4.

Other results, all in figures only (described, not tabulated):

• Synthetic data: as the anomaly proportion grows, SISVAE-0 and STORN drop significantly, while the regularised models hold their AUPRC §5.2 fig. 4a. λ = 0.5 is best on the 2 % set §5.2 fig. 4b.
• Convergence: the KL-regularised model fluctuates less than the unregularised and mean-regularised models after about 100 iterations §5.2 fig. 4c fig. 4d.
• Ranking across datasets and metrics, with standard deviations of ranks: fig. 5. The paper says the model outperforms the others notably.
• Precision@K for K = 10, 50, 200: SISVAE-p stays stable as K increases; Donut fails on A2–A4 at K = 10 §5.3 fig. 8.
• Mackey-Glass (univariate): SISVAE-p beats Donut for AUPRC and ROC after Donut is tuned to z_dim = 5 §5.5 fig. 13.

Paper vs code

No official code. The paper gives no code link, and the arXiv metadata has no comment. GitHub searches: SISVAE returned four repositories (Gabriel210121/SISVAE, Gabriel210121/Sisvaen, manuel9716/sisvae-2.0-dev, DavidTan-cloud/wqf7023-industrial-anomaly-detection-radar), all with owners who are not among the paper's authors, and with no license. None was used. smoothness inducing sequential variational and time series anomaly detection sequential VAE smoothness Longyuan returned no repositories.

Note: the paper says it used "the code open sourced by the authors" for Donut §5.1. That refers to Donut's authors' code, not this paper's code.

Pitfalls

• Sign of the smoothness term. Eq. 20 adds λ·KL to the terms of the ELBO eq. 17 eq. 20. Eq. 14 maximises the bound, so the printed form rewards divergence between consecutive reconstructions. The text describes the term as a penalty on non-smooth outputs §4.2 §5.2. No code is available to settle this. The card's pseudocode keeps the printed sign and flags it; check against the text's intent before implementing.
• Wrong equation cited for λ. The λ sweep points to eq. 17 ("Eq. objective") §5.2, which has no smoothness term. The objective with λ is eq. 20.
• Monte Carlo score as written. Eq. 25 puts the log outside the average without a log inside eq. 25. Algorithm 2 averages the negative log-likelihoods over L draws alg. 2. These are different scores: eq. 25 is the log of an average, while Algorithm 2 averages the logs. The paper does not say which one it used, but Algorithm 2 is the implementable one.
• "Sequential" Monte Carlo is plain sampling. The paper names SMC §4.4 but Algorithm 2 has no particle weights and no resampling; each draw comes from the encoder alg. 2.
• GRU equations contain typos. The s_t gate uses b_r where b_s is likely meant eq. 13. The candidate state's W_h term has no input argument, unlike the other gates eq. 13. The paper does not state which form was implemented.
• Inference density notation. Eq. 16 writes p(z_t | x_t) where q_φ is meant eq. 16.
• Chunking. Training uses sliding windows with step s, which is not stated. Scoring uses non-overlapping chunks §4.4. How the tail of a series is handled is not stated.
• Smoothness at t = 1. The paper uses the same distribution on both sides, which makes that term zero §4.2.
• Threshold selection on labelled test data. The reported F1 is the best value over all thresholds, chosen with the test labels §5.3. No validation split is described. This biases F1 upward. AUROC and AUPRC do not depend on a chosen threshold, but they still use the labels.
• Donut was run without labels and in a flattened protocol. Donut's own method uses labels to learn normal patterns; the paper removed them for fairness §5.1. Donut's numbers here are therefore not its published setup. The paper chose the flattened protocol because it did better on this data §5.1.
• Single-run reporting. No seed count, variance or confidence interval is given in Table IV. Figure 5 shows rank variability, not run-to-run variance fig. 5.
• Non-stationary variance. The paper says the variance is modelled per time step §1 §4.4. How σ is kept positive is not stated.
• Unstated optimiser and batch details. Mini-batch size, Adam betas, learning-rate schedule and gradient clipping are not stated. The paper describes its training as "200 epochs" without giving the batch size.
• Mackey-Glass not fully specified. τ and n in eq. 26 have no values §5.5.
• Sensitivity to anomaly rate and λ. The paper shows that the regularised model is more stable at higher anomaly proportions §5.2 fig. 4a. The λ = 0.5 value was tuned on synthetic data with labels §5.2 §5.1. Whether real-data runs were tuned is not stated.
• Dataset access. The real datasets are Yahoo S5 Webscope and an open μPMU dataset. The paper gives only footnote URLs and does not state their terms of use §5.3.
• Old framework. The code targets PyTorch 0.4 §5.1. Porting requires checking the RNN and sampling primitives used.

Lineage

• Builds on: the VAE and AEVB [arXiv:1312.6114] §2 eq. 7 eq. 14; the SGVB estimator §3.3 eq. 8; sequential VAEs such as VRNN and STORN §2. STORN's arXiv id as listed in the paper's bibliography: 1602.07109.
• Smoothness regularisation idea: classical smoothing and de-trending for deterministic models, which the paper extends to a probabilistic setting §2 §4.2.
• Compared against: Donut §5.1 (VAE with sliding windows, the modified ELBO), STORN §5.1 (sequential VAE with independent latent sampling), VRNN (cited in §2 and not compared in the experiments §5.1), and the statistical baselines HA, ARMA and LDS §5.1 tab. 2 tab. 4. SISVAE-0 and VAE-s are the paper's own ablations §5.1.
• Superseded by: not stated in this wiki.