Graph Neural Network-Based Anomaly Detection in Multivariate Time Series

In plain words

A water plant or factory has many sensors, and they affect each other. Their readings are linked in complex, non-linear ways. This paper detects attacks and faults by checking whether those expected links still hold.

The method has two parts. First it learns which sensors should be linked, from how similar their learned fingerprints are. Then it predicts what each sensor should read next, using the readings of its linked sensors. When the real reading moves away from the prediction, that counts as a deviation.

The score for a moment is driven by the sensor that deviates most, and the sensor's links show which other sensors it normally moves with. That gives an operator a place to start looking. The authors compare their results with several standard anomaly detectors and report better numbers on two water-treatment datasets. Those numbers are the authors' own.

TL;DR

• Learns a directed graph of sensor dependencies from similarity of learned sensor embeddings, keeping the top k neighbours per sensor §3.3 §3.4.
• Forecasts each sensor's next value with graph attention over its neighbours and a per-sensor embedding gate §3.5.
• Flags a time tick when its smoothed score, the largest per-sensor deviation normalised by median and IQR, exceeds a threshold set on validation data §3.6 §4.3.
• Author-reported. Trained on normal data only. Reported F1 is 0.81 on SWaT and 0.57 on WADI. The SWaT and WADI data are not in the repository and must be requested from iTrust README.md.

Problem

Anomaly detection in high-dimensional sensor data is usually unsupervised, because labelled attacks are scarce §1. Classical linear and distance methods model sensor relationships too simply §1 §2. Deep autoencoders, GANs and LSTM models do not learn explicitly which sensors relate to each other, so they struggle to say which relationship broke §1 §2. Graph neural networks take a graph as input, but here the graph is unknown and must be learned §1.

Method

The model has four components: sensor embedding, graph structure learning, graph-attention forecasting and graph deviation scoring §3.2. Figure 1 gives an overview fig. 1. The paper's experiments use a sliding window of size w = 5 §4.4.

1. Input windows §3.1. At time t, the input of sensor i is its last w values. The target is the value at time t eq. 4. Training data are assumed to be normal §3.1.

2. Sensor embedding §3.3. Each sensor i has a learned vector v_i in R^d. These vectors start random and train with the rest of the model §3.3. Similar vectors should mean related behaviour §3.3.

3. Graph structure learning §3.4. Each sensor i has a candidate set C_i of possible neighbours eq. 1. Without prior knowledge, C_i is every other sensor §3.4. For each candidate j, the cosine similarity e_ji between v_i and v_j is computed eq. 2. The top k candidates by this similarity become neighbours, giving a directed edge j to i in the adjacency matrix A eq. 3 §3.4. The paper notes that edges need not be symmetric §3.4.

4. Graph-attention forecasting §3.5. Each sensor's window is projected with a shared matrix W eq. 5. The representation z_i sums the projected windows of i and its neighbours, weighted by attention, and applies ReLU eq. 5. The attention for neighbour j uses g_i, which is the concatenation of v_i and W x_i eq. 6 eq. 7. Scores are passed through LeakyReLU and softmax over N(i) together with i itself eq. 7 eq. 8. The attention therefore depends on both sensors' embeddings and their projected windows §3.5.

5. Output layer §3.5. Each z_i is multiplied element-wise by v_i §3.5. The gated vectors of all sensors go into stacked fully connected layers that output the N predicted next values eq. 9.

6. Training loss §3.5. The squared error between prediction and observation is summed over sensors and averaged over time eq. 10.

7. Deviation scoring §3.6. The absolute error of each sensor is Err_i(t) eq. 11. Each sensor's error is normalised by its median and IQR eq. 12. The score for a tick is the maximum over sensors, so one sensor is enough to raise it eq. 13 §3.6. A simple moving average of these scores gives the smoothed score §3.6. A tick is an anomaly when the smoothed score exceeds a threshold, which the paper sets to the maximum smoothed score on validation data, to avoid adding a hyperparameter §3.6 §4.3.

8. Explanation §3.6 §4.8. The deviating sensor with the highest score is the central sensor. Its strongest attention neighbours are shown, and its predicted and observed values are compared over the attack period fig. 3 §4.8. Sensor embeddings can also be projected with t-SNE to show groups of similar sensors fig. 2 §4.7.

# Learn directed neighbours from embedding similarity  [§3.4] [eq. 2] [eq. 3]
# Input:  embeddings (N, d), k
# Output: neighbours (N, k), sensor indices

def learn_neighbours(embeddings, k):
    # Cosine similarity between every pair of sensor embeddings  [eq. 2]
    unit = embeddings / norm(embeddings, axis=1, keepdims=True)  # (N, d)
    similarity = unit @ unit.T                                   # (N, N)
    # Keep the k most similar sensors for each sensor  [eq. 3]
    neighbours = top_k_indices(similarity, k)                    # (N, k)
    return neighbours

# One forecast for all sensors  [§3.5] [eq. 4] to [eq. 9]
# Input:  windows (N, w), embeddings (N, d), neighbours (N, k)
# Output: predicted next values (N,)

def forecast(windows, embeddings, neighbours, W, a, f_theta):
    # Shared linear projection of each sensor's window  [eq. 5]
    projected = windows @ W.T                                    # (N, d)
    # Key of each sensor: its embedding joined with its projection  [eq. 6]
    keys = concat(embeddings, projected, axis=1)                 # (N, 2d)
    features = zeros((N, d))
    for i in range(N):
        # Neighbours of i, with i itself included  [eq. 8]
        group = append(neighbours[i], i)
        # Attention score of each member of the group  [eq. 7]
        scores = [leaky_relu(a @ concat(keys[i], keys[j])) for j in group]
        # Normalise the scores over the group  [eq. 8]
        weights = softmax(scores)
        # Weighted sum of projected windows, then ReLU  [eq. 5]
        features[i] = relu(
            sum(w_m * projected[j] for w_m, j in zip(weights, group))
        )
    # Gate each sensor's features by its own embedding  [eq. 9]
    gated = features * embeddings                                # (N, d)
    # Stacked fully connected layers over all gated vectors  [eq. 9]
    predicted = f_theta(flatten(gated))                          # (N,)
    return predicted

# Training objective  [§3.5] [eq. 10]
# Input:  predicted (B, N), observed (B, N)
# Output: scalar loss

def training_loss(predicted, observed):
    # Squared error for each sensor and window  [eq. 10]
    errors = (predicted - observed) ** 2                         # (B, N)
    # Sum over sensors, mean over windows  [eq. 10]
    return mean(sum(errors, axis=1))

# Overall anomaly score per time tick  [§3.6] [eq. 11] to [eq. 13]
# Input:  predicted (T, N), observed (T, N), median (N,), iqr (N,)
# Output: smoothed score (T,)

def anomaly_scores(predicted, observed, median, iqr):
    # Absolute error of each sensor at each time tick  [eq. 11]
    errors = abs(observed - predicted)                           # (T, N)
    # Robust per-sensor normalisation  [eq. 12]
    normalised = (errors - median) / iqr                         # (T, N)
    # The largest sensor deviation represents the tick  [eq. 13]
    scores = normalised.max(axis=1)                              # (T,)
    # Smooth with a moving average; the window is not stated  [§3.6]
    return moving_average(scores)

# Threshold from held-out normal data  [§3.6] [§4.3]
# Input:  smoothed scores on validation data
# Output: threshold

def threshold(validation_scores):
    # Largest smoothed score seen on validation data  [§4.3]
    return max(validation_scores)

Hyperparameters

Values from the paper unless marked "code". The code values are the ones in run.sh, which the repository's README uses as its example run.

Results

Author-reported. Setup: two water-treatment test beds with simulated attacks, SWaT (51 features) and WADI (127 features). Both are down-sampled to one value per 10 seconds §4.1. The paper uses two weeks of normal data for training and tests on the attack period §4.1.

Dataset statistics tab. 1

Detection accuracy on the test set tab. 2. Precision and recall are percentages, F1 is a fraction. Part of the baseline results are taken from the MAD-GAN paper, the table caption says, without naming which rows tab. 2.

The text reads precision as 0.99 on SWaT and 0.98 on WADI §4.5, which matches the table after rounding. It reports 54% higher F-measure than the next best baseline on WADI §4.5. The table's WADI F1 of 0.57 against 0.37 gives that ratio.

Ablation tab. 3 §4.6

Case study §4.8 fig. 3. A flow sensor on WADI was attacked with false readings inside its normal range. The model's highest-scoring sensor was a valve status sensor that is normally strongly related to the attacked flow sensor. The predicted value for the valve rose with the flow, but the observed value stayed flat. fig. 2 shows embedding clusters of sensors with similar roles on WADI §4.7.

Paper vs code

The official repository is https://github.com/d-ailin/GDN (MIT), linked from the paper's footnote §1. The commit is 9853899. Its LICENSE copyright line reads "d-ailin" LICENSE.

• Adam betas. Paper: (0.9, 0.99) §4.4. Code passes no betas, so PyTorch's default (0.9, 0.999) applies train.py.
• Epochs and patience. Paper: up to 50 epochs, patience 10 §4.4. Code: run.sh uses 30 epochs run.sh and early stopping waits 15 epochs train.py. The argparse default is 100 main.py.
• Per-dataset settings. Paper: d, k and hidden width differ between WADI and SWaT §4.4. The repository has one configuration, run.sh, with d = 64, k = 5 and hidden width 128 for any dataset run.sh run.sh run.sh. The per-dataset values in the paper are not reproduced by the shipped script.
• Hidden layers. Paper: "hidden layers of 128 (64) neurons" §4.4. Code: run.sh sets out_layer_num=1 run.sh. With one layer the output is a single linear map per sensor and the hidden width is not used GDN.py.
• Candidate set and self-loops. Paper: the top k come from the candidate set, which excludes the sensor itself eq. 1 eq. 3 §3.4. Code: the top k are taken over all sensors, including the sensor itself GDN.py. The self-loop is removed and added back in the layer graph_layer.py. A sensor therefore usually gets k − 1 learned neighbours plus itself.
• Prior graph unused. The code builds a fully connected candidate graph net_struct.py main.py. The batched version is computed GDN.py but never used in message passing. The prior-candidate option the paper describes is therefore not active in the released code.
• Output layer. Paper: a stacked fully connected network over the concatenation of all N gated vectors eq. 9. Code: the output layer runs on each sensor's vector separately, with weights shared across sensors GDN.py GDN.py. Sensors mix only through the graph layer.
• Layers not in the paper. Code adds batch normalisation and ReLU after aggregation GDN.py. It adds a bias to the graph layer graph_layer.py graph_layer.py. It adds batch normalisation, ReLU and dropout of 0.2 before the output layer GDN.py GDN.py GDN.py. Equation 5 has only a ReLU and no bias eq. 5.
• Embedding initialisation. Paper: random §3.3. Code: Kaiming uniform GDN.py.
• Loss scale. Paper: sum of squared errors over sensors eq. 10. Code: mean over every element, which also divides by N train.py train.py.
• Median and IQR source. Paper: computed "across time ticks" §3.6. Code: computed from the test-set errors of each sensor evaluate.py. Validation scores are normalised with their own statistics evaluate.py.
• Smoothing. Paper: simple moving average, window not stated §3.6. Code: causal mean over four values, the current one and three before it. The first three scores are set to zero evaluate.py.
• Threshold mode. Paper: validation maximum §3.6 §4.3. Code: the validation mode does this evaluate.py main.py. The default report='best' instead picks the threshold that maximises F1 on the test labels run.sh main.py evaluate.py.
• Preprocessing labels. Paper: each 10-second bin takes its most common label §4.1. Code: a bin is an attack if any sample in it is an attack process_swat.py process_wadi.py.
• Normalisation. The paper does not describe any scaling. The preprocessing scripts apply min–max scaling fitted on the training data, after filling missing values with the column mean and then with zero process_swat.py process_swat.py process_wadi.py.

Pitfalls

• Oracle threshold. The default report='best' chooses the threshold on the test labels, giving the best F1 on the test set evaluate.py. A deployable number needs report='val'. The paper says the threshold is the validation maximum §4.3, but it does not say which mode produced Table 2. Treat the table as possibly optimistic.
• Test-set normalisation. Median and IQR come from the test errors evaluate.py. For a live stream these statistics would have to come from normal data instead. The paper's wording does not settle this.
• Tiny IQR. The code adds 1e-2 to the IQR evaluate.py. The paper's formula has no such term eq. 12. A sensor with a very small IQR gives very large normalised scores, so check the scale per sensor.
• Startup gap. The first three smoothed scores are left at zero, so the first three ticks of each test series get no real score evaluate.py.
• Neighbour count. Because the self index is usually among the top k, k = 5 in run.sh gives four real neighbours GDN.py. Compare k with the paper's wording only after this correction.
• Validation split. Paper: not stated, and the paper marks it as an open question §4.3. Code: 20% of the training windows as one contiguous block at a random start, using the seed main.py main.py. Adjacent windows share time steps, so the split is not strictly disjoint in time.
• Model selection. The checkpoint kept is the one with the lowest validation MSE, not the best F1 train.py.
• Data access. SWaT and WADI must be requested from iTrust README.md. The repository ships only a demo subset, the msl directory msl. Table 2 cannot be reproduced from the repository alone.
• Baseline count. The paper says it compares with five baselines §4.2 but lists seven §4.2. The baseline code is not in the repository, and the caption says some baseline numbers come from another paper tab. 2.
• No variance. Results come from one configuration, with a single seed in the released script run.sh. No standard deviations are reported tab. 2.
• Cost. Each forward pass builds an N × N similarity matrix GDN.py and a batched edge list whose size grows with batch size times N times k GDN.py. The pinned library versions are from 2020 README.md and were not tested here.
• Equation numbers. The equation numbers above were counted in order from the TeX source. The rendered PDF was not checked.

Lineage

• Builds on: graph convolution networks §2 and graph attention networks §2, which use a fixed graph as input. This paper learns the graph and adds sensor embeddings to attention §3.4 §3.5.
• Compared against: PCA, KNN, feature bagging, autoencoders, DAGMM, LSTM-VAE and MAD-GAN §4.2 tab. 2.
• Superseded by: not stated in this card's sources.