SARAD: Spatial Association-Aware Anomaly Detection and Diagnosis for Multivariate Time Series

Authors: zhihao  dai
ligang  he
shuang-hua  yang
matthew  leeke
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 18:09:20

Abstract: Anomaly detection in time series data is fundamental to the design, deployment, and evaluation of industrial control systems. Temporal modeling has been the natural focus of anomaly detection approaches for time series data. However, the focus on temporal modeling can obscure or dilute the spatial information that can be used to capture complex interactions in multivariate time series. In this paper, we propose SARAD, an approach that leverages spatial information beyond data autoencoding errors to improve the detection and diagnosis of anomalies. SARAD trains a Transformer to learn the spatial associations, the pairwise inter-feature relationships which ubiquitously characterize such feedback-controlled systems. As new associations form and old ones dissolve, SARAD applies subseries division to capture their changes over time. Anomalies exhibit association descending patterns, a key phenomenon we exclusively observe and attribute to the disruptive nature of anomalies detaching anomalous features from others. To exploit the phenomenon and yet dismiss non-anomalous descent, SARAD performs anomaly detection via autoencoding in the association space. We present experimental results to demonstrate that SARAD achieves state-of-the-art performance, providing robust anomaly detection and a nuanced understanding of anomalous events.

Main paper: 
 introduction: 
 Time series anomaly detection is critical for industrial automation (Rieth et al., 2018) , intrusion detection (Mathur and Tippenhauer, 2016) , and healthcare sensing (Goldberger et al., 2000) . Anomaly detection in these contexts is typically treated as an unsupervised learning problem, owing to the novelty of anomalies and the scarcity of labeled anomalies.Temporal modeling is the mainstream basis of current time series anomaly detectors. By learning the dependencies between discrete time steps, temporal modeling can pinpoint the time spans of anomalies. During anomalies, unseen and peculiar temporal dependence patterns degrade data autoencoding (Wang et al., 2023b) or autoregression (Zhao et al., 2020) performance, thereby enabling detection. Alternatively, irregular temporal representations can be driven to breach learned enclosing h by final L-th layer's MHSA are averaged across heads to derive A before (1b), during (1c), and after (1d) the anomaly. Darker cells have larger values. Anomalous features (#12 and #15) are highlighted with red bounding boxes. The reduction-only changes from before or after the anomaly to during the anomaly are shown in 1e and 1f, i.e., ReLU(A pre -A in ) and ReLU(A post -A in ). The anomaly leads to association reductions on anomalous features, prominently column-wise on A.hyperspheres (Shen et al., 2020) , resulting in high anomaly scores that enable detection. Despite its temporal precision in anomaly detection, temporal modeling either assumes feature independence or combines variables of diverse physical nature, the former simplifying the modeling and the latter mitigating the multicollinearity issue. Such assumptions lead to either omission or dilution of spatial information crucial to anomaly detection. Specifically, it overlooks the long-time-range spatial associations, the relationships between various features which ubiquitously characterize normal behaviors of multivariate time series. Where anomaly detection pinpoints the temporal locations of an anomaly, anomaly diagnosis identifies the spatial locations, i.e., the anomalous feature set, of an anomaly. Temporal methods also restrict diagnostic capabilities, as the lack of spatial information mismatches autoencoding-based or autoregression-based anomaly criterion, which de facto measures temporal novelty, with its objective of capturing spatial novelty. Furthermore, time series anomalies frequently dissolve spatial associations, motivating anomaly detection in the association space. Using an vanilla Transformer (Vaswani et al., 2017) , we investigate the changes in spatial associations throughout anomalies. Applied on transposed time windows (the spatial dimension comes before the temporal), an encoder-only Transformer is trained to minimize reconstruction errors on unlabeled N -variate time series and, by doing so, learns to model the multivariate series spatially via the Multi-Head Self-Attention (MHSA) illustrated in Figure 2 . MHSA at each stacked l-th layer computes an intermediate association mapping A l h ∈ R N ×N per h-th head, mapping back input X to produce attention scores. The last layer's mapping A L h thus effectively captures the contributions of k-th feature to the reconstruction of j-th feature at each location (j, k), not least for its architectural proximity to the reconstructed output. Recent research (Liu et al., 2024 ) also highlights the important role MHSA plays in capturing the inter-feature associative relationships when applied on the multi-variate dimension. As new associations emerge and old ones dissolve over time, Figure 1 shows the association changes on a real-world benchmark. We observe that anomalies exhibit reductions for anomalous features, a phenomenon we herein coin as Spatial Association Reduction (SAR). The rationale is that anomalies either originate from or result in dissolution of pre-existing associations, detaching anomalous features from their non-anomalous counterparts. Additionally, we make the observation that SAR is most prominent column-wise on A L h , since each j-th column characterizes the dropouts of j from associating with other, mostly non-anomalous, features. Due to lack of explicit spatial information, temporal modeling is inadequate for exploiting SAR. More examples are given in Appendix C.From a spatial modeling perspective, we propose SARAD to leverage spatial information and to exploit SAR for enabling robust time series anomaly detection and diagnosis. For quantifying anomalous spatial novelty in the data space, we train a Transformer on transposed time windows as an autoencoder. To capture the spatial association progression, the reduction-only changes of associations over time, the data reconstruction divides the input window by time into two halves to be processed in parallel. Consequently, the progression is the non-negative backward difference of the intermediate association mappings via MHSA. Subseries division circumvents memory storage of latest association mappings and enables time window shuffling during training, which reduces order bias, enhances generalization, and prevents catastrophic forgetting. For quantifying anomalous reduction novelty, we train a Multi-Layer Perceptron (MLP) as an antuencoder on progression in the association space. Whereas progression encompass all association reduction, autoencoding rules out those not caused by anomalies. The reconstruction errors via the data module measure data-only anomalous deviation from expected system behaviors and falter when such deviations are not prominent, e.g., at the start of an anomaly. The reconstruction errors via the progression module are sensitive to change in spatial associations, thus complementing the former. We develop a joint anomaly detection criterion that combines both. Experiments show SARAD delivers stateof-the-art detection and diagnosis performance with architectural elegance. Code is available at https://github.com/daidahao/SARAD/. We summarize our contributions as follows.• We reveal and extract spatial association descending patterns of time series anomalies with a bespoke Transformer and subseries division. The former learns the pairwise inter-feature associations via autoencoding in the data space and the latter enables shuffled autoencoding training and memory-less progression aggregation.• We propose progression autoencoding to quantify anomalous descent in the association space and a joint detection criterion in both data and association spaces, which complement each other.• Experimentally, SARAD performs state-of-the-art anomaly detection and diagnosis on multivariate time series and ablation studies support our design choices. 
related work: 
 Influenced by the dominance of temporal modeling in time series forecasting (Wang et al., 2023a; Zhang et al., 2023; Wu et al., 2021) , temporal modeling is also prevalent in time series anomaly detection. Recurrent neural networks such as LSTM (Hochreiter and Schmidhuber, 1997) have innate capabilities for handling sequential data. These approaches use hidden states for past input memorization, enabling detection (Li et al., 2019; Malhotra et al., 2015) and diagnosis (Qian et al., 2021) . Transformer (Vaswani et al., 2017) network is widely adopted (Fan et al., 2023; Xu et al., 2022) approach that is commonly applied to model temporal associations between different time points using its attention mechanism. Linear regression (Zeng et al., 2023) and MLP (Wang et al., 2024; Audibert et al., 2020) directly model temporal dependencies. TranAD (Tuli et al., 2022) replaces the MLP in Audibert et al. (2020) with a Transformer, making the detection criterion more robust through its adversarial training paradigm. Temporal modeling, however, is restricted by the exceptionally small receptive field in time and adversely impacted by the timestamp misalignment across features. In the context of anomaly detection, temporal modeling helps capture anomalous temporal associations (Xu et al., 2022; Yang et al., 2023) , but offers limited detection capabilities in absence of spatial information. In a diagnostic context, temporal detectors mismatch anomaly criterion of temporal novelty with spatial interpretation.Spatial associations characterize the multivariate time series commonly found in such supervisory systems for industrial control. The relationships range from strongly correlated, e.g., due to spatial proximity, to fully independent, e.g., due to mechanical disconnection. For forecasting, iTransformer (Liu et al., 2024) applies Transformer on the transposed time series to enable direct spatial modeling. Crossformer (Zhang and Yan, 2023) screens the time series through custom Two-Stage Attention layers for more efficient spatial modeling. In terms of detection, GDN (Deng and Hooi, 2021) learns a directed graph of features for the prediction of last time points, whose errors serve as anomaly scores. GDN is partially limited by a mismatch between its singe-timestamp prediction target and the prevalent range-wise anomalies as well as unstable Top-K node selection during training. InterFusion (Li et al., 2021) learns compressed spatial and temporal dependencies, using a hierarchical Variational Auto-Encoder (Kingma and Welling, 2014) to reconstruct the series. Neither inspects temporal changes in associations throughout anomalies.On another front, Isolation Forest (IF) models build a binary decision tree ensemble by partitioning either the data space (Liu et al., 2008) or the deep embedding space (Xu et al., 2023) formed by randomized neural networks. They are constrained by the lack of temporal and spatial (in the former case) or spatial (in the latter case) information, and their anomaly scores are not reflective of the degrees of anomalies.We emphasize anomalous association descending patterns towards better time series detection and diagnosis. Different from previous work, we explicitly utilize the reduction in spatial associations over time during an anomaly, an insight we derived from the cyber-physical defense space. Dynamic watermarking (Satchidanandan and Kumar, 2017) and similar defense techniques (Dai et al., 2023) overlay actuation with private signals to reveal attacks resulting in correlational breakdowns. While their approaches are intrusive and actively alter system behaviors, our detector remains non-intrusive, passively monitors the spatial associations, and is applicable to any supervisory system.We refer to spatiality in this work as the multi-dimensional vector nature inherent to multivariate time series data. The terminology is also used in literature on time series related tasks (Gangopadhyay et al., 2021; Zheng et al., 2023) . We note that spatiality may carry different meanings in other AI contexts, such as geographic positions or characteristics on Earth. We differentiate those meanings from our definition of spatiality, which traces its root to the spatial distribution of sensors and actuators in control systems where time series are routinely collected. 
method: 
 The problems of anomaly detection and diagnosis are specified as follows.Anomaly Detection Given a N -feature time series T = {x 1 , • • • , x N } where x n ∈ R T is of the same length T , the objective is to predict the anomaly label y t ∈ {0, 1} at each timestamp t.Given the same time series T , the objective is to predict the diagnosis label g t ⊆ [N ], the set of anomalous features at each timestamp t.SARAD comprises two sequential modules; a Transformer for time series data reconstruction and a MLP network for spatial progression reconstruction. Table 1 decomposes the system framework of SARAD. The Transformer temporally divides by 2 and reconstructs the input time series to learn pairwise inter-feature associations and to enable order-free memory-efficient progression aggregation.The MLP reconstructs the aggregated progression to quantify anomalous association reduction while dismissing non-anomalous reduction. Towards robust anomaly detection, reconstruction errors from the two modules jointly serve as a criterion, sensitive to data-only anomalous deviation and anomalous association reduction.In light of restricted capabilities of temporal detectors, here we adapt Transformer to spatially reconstruct the series data. The data module contains two components, Subseries Split & Merge and Subseries Reconstruction, shown in the first and second columns in Table 1 . The former wraps around the second by temporally splitting a multivariate input series in half at its beginning and temporally merging at its end. Subseries division enables capturing of spatial progression within a single time window. Without the former, the model must store in memory the last association mappings at each step and keep to the time ordering during training, which is prone to overfitting and min || Ŝ -S|| 2 2 Detection Criterion || X -X|| 2 2 || Ŝ -S|| 2 2 Diagnosis Criterion || X(j,•) -X (j,•) || 2 2 || Ŝ(•,j) -S (•,j) || 2 2catastrophic forgetting. The latter is an encoder-only Transformer composed of an embedding layer, a L-layer attention-based encoder, and finally a linear projection layer. Encoder-only Transformers are commonly found in Transformer-based detectors (Kang and Kang, 2024; Kim et al., 2023) due to its simplicity and the length uniformity of the target output, i.e., the reconstructed series.Ours exclusively models spatial associations, unconventional to the aforenamed detectors and most temporal forecasters (Nie et al., 2023; Zhou et al., 2022; Liu et al., 2022 ) and yet more aligned with recent spatial-aware forecasters (Liu et al., 2024; Zhang and Yan, 2023) . At each encoding layer, MHSA computes pairwise association mappings, a representation of inter-feature dependencies which ubiquitously characterize the multivariate time series and are crucial to anomaly detection.We suppose the input series is a time window X ∈ R 2W ×N of length 2W , where 2 is for the convenience of a temporal split. Before reconstruction, X is split into two half multivariate subseries of equal temporal length:X = {X 1 ∈ R W ×N , X 2 ∈ R W ×N }.After subseries reconstruction, the two reconstructed subseries are concatenated to form the full reconstructedX = { X1 ∈ R W ×N , X2 ∈ R W ×N }.Embedding To lead subseries reconstruction, each X i is embedded asX 0 i = E i + M , wherein E i = Linear(X T i ) ∈ R N ×D and M = {m i ∈ R D |i ∈ [N ]} is a learnable feature-level embedding.Spatial-Aware Encoding A stack of L Transformer encoding layers is used to encode the series in the D-length attention space. Each layer is stacked with MHSA and MLP with residual connections:EQUATIONwhere X l-1 i , Z l i , X l i ∈ R N ×D are l-1 layer's output, l-th layer's hidden state and output respectively and LN (•) is Layer Normalization (Ba et al., 2016) . It is spatial-aware because MHSA explicitly learns a pairwise inter-feature association mapping to exchange information in between feature representations. Notably within the MHSA with H heads as shown in Figure 2 , each (j, k)-th element on the association mapping A h ∈ R N ×N of its h-th head computes how much j-th feature's attention scores should originate from the k-th feature's key. From a broader perspective of data reconstruction, it is the quantification of the residual impact of k-th feature's originals on the j-th feature's reconstructions. We refactor MHSA implementation to enable parallel encoding of subseries.Linear Projection Output from the last encoding layer X L i is linearly projected and transposed to derive the reconstructed subseries Xi ∈ R W ×N .To exploit SAR caused by anomalies, the module first extracts and aggregates the spatial progression, the non-negative backward difference of the association mappings via MHSA. In line with general anomaly detectors (Aggarwal, 2013) , the module conducts autoencoding in the association space to quantify anomalous SAR and to dismiss non-anomalous SAR. Anomalous SAR occurs when, say, a compromised sensor's readings are no longer correlating with its spatially adjacent or mechanically related counterparts.We define association progression S l h ∈ R N ×N at the h-th attention head in the l-th layer to be the non-negative backward difference in association mappings {A l 1,h , A l 2,h }:EQUATIONwhere ReLU(•) passes through only non-negative values and outputs zeros otherwise.Progression Aggregation To center the detection on association dropouts, we aggregate the column sums of progression S l h from all attention heads in the final L-th layer to form S ∈ R H×N :EQUATIONWe recall from the data module, each k-th column in A l i,h quantifies the impact of k-th feature on all features' reconstruction. Taking the sum per each k-th column, we measure with S the dropout rates of k from participating all features' reconstruction. As we have observed in Section 1, SAR at the column level is indicative of time series anomalies, more so than at the row level. The last layer's progression is focused not least for its proximity to the final reconstructed output, whereafter no more information is exchanged between features. Liu et al. (2024) manifests that final layer's mappings resemble closely with the inter-feature correlations of the target, in our case, the reconstructed.Autoencoding With S flattened as an one-dimensional vector, a 2-layer MLP is trained to output the reconstructed and reshaped Ŝ ∈ R H×N , synchronously with the data module training.Training Objective We train an end-to-end model with a joint minimization objective:L R = || X -X|| 2 2 , L S = || Ŝ -S|| 2 2 , L = L R + λ L S L S (4)where L R is the data reconstruction loss, L S the progression reconstruction loss, and λ L S a weight hyper-parameter. Gradients are stopped from flowing into S to prevent updates to the data module and collapses in association representation. We are training two anomaly detectors simultaneously, one working in the original data space, the other in the spatial progression space.Anomaly Detection Criterion For an input series X, the anomaly score s is a scalar defined to be:EQUATION)where r is the data reconstruction error, p the progression reconstruction error, µ r , µ p the means of r, p on the validation set, and σ r , σ p the standard deviation of r, p. The criterion takes into account the normalized errors in the data space and the progression space, each of which quantifies the anomalous magnitude in respective spaces and complment the other. Anomaly Diagnosis Criterion For an input series X and its j-th feature, its anomaly score s j is a scalar defined to be:EQUATIONwhere r j is feature j's data reconstruction error. The criterion is sensitive to spatial novelty.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Answer: [NA]• It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments: 
 SARAD is compared against state-of-the-art detectors on real-world benchmarks for detection and diagnosis, the latter only when diagnostic labels are available.Datasets We evaluate on four real-world datasets collected under industrial control and service monitoring settings. These dataset are: 1) Server Machine Dataset (SMD) (Su et al., 2019b,a) , 2) Pooled Server Metrics (PSM) dataset (Abdulaal et al., 2021a,b) 3) Secure Water Treatment (SWaT) dataset (Mathur and Tippenhauer, 2016; iTrust, 2023) , and 4) Hardware-In-the-Loop-based Augmented ICS (HAI) dataset (Shin et al., 2021b,a) . All training sets contain only unlabeled data and the test sets contain data with anomaly labels. Anomalies range from service outages to external cyber-physical attacks. We summarize the statistics of the datasets in Table 2 . Descriptions of each dataset are detailed in Appendix E.Detection Metrics Real-world benchmarks are rife with range-wise anomalies spanning consecutive time points (Wagner et al., 2023) . We use the range-based metrics proposed in (Paparrizos et al., 2022) . Compared against their point-based counterparts, they provide robustness to labeling delay and scoring noises as well as performant detector separability and series consistency. We compute the threshold-independent AUC-ROC and AUC-PR scores to be rid of thresholding impact and fully parameter-free Volume Under the Surface (VUS) AUC-ROC and AUC-PR scores. Full details are discussed in Appendix I.Diagnosis Metrics Consistent with previous works (Tuli et al., 2022; Zhao et al., 2020) , we use common metrics such as Hit Rate (HR) (Su et al., 2019b) and Normalized Discounted Cumulative Gain (NDCG) (Järvelin and Kekäläinen, 2002) where diagnosis labels are available. At the range level, we measure the Interpretation Score (IPS) initially proposed in Li et al. (2021) and here expanded to fit the P % parameterization. Full details are discussed in Appendix J.Baselines We compare SARAD against state-of-the-art anomaly detection baselines, including Isolation Forest-based IF (Liu et al., 2008) , Deep IF (DIF) (Xu et al., 2023) ; MLP-based USAD (Audibert et al., 2020) ; graph-based GDN (Deng and Hooi, 2021) ; LSTM-based MAD-GAN (Li et al., 2019) ; CNN-based DiffAD (Xiao et al., 2023); and Transformer-based TranAD (Tuli et al., 2022) , ATF-UAD (Fan et al., 2023) , AT (Xu et al., 2022) , DCdetector (Yang et al., 2023) . Noticeably, GDN employs explicit spatial modeling in its graph construction although spatial associations are not directly involved in anoamly scoring. MAD-GAN emphasizes on anomaly detection within cyber-physical systems. All baselines are trained using official implementations where available and recommended hyperparameters from respectively papers are used. Anomaly Detection Table 3 shows the anomaly detection performance in metrics defined in Section 4.1. It demonstrates that, despite its architectural elegance, SARAD either outperforms all baselines by significant margins on the threshold-independent VUS-ROC scores (SMD: +15.51% MAD-GAN, HAI: +9.92% DiffAD) or performs on par with current best detectors (PSM: -1.36% GDN, SWaT: -0.36% DIF). IF scrutinizes the distributional shifts of anomalies with random data partitions and delivers consistent performance across datasets. DIF extends IF into randomized deep representation spaces and archives decent improvements due to more flexible partitions and temporally local information extraction via dilated convolutions. Temporal modeling methods such as DiffAD, ATF-UAD, and AT rely solely or heavily on reconstruction errors and when the errors do not correspond the the underlying anomalies their performance plummet. Adversarial training in USAD and MAD-GAN amplifies reconstruction errors of anomalies to mitigate but not eliminate such issues and thus suffer less performance drops. In contrast, our SARAD additionally accounts for the SAR frequent with anomalies and independent of data distributional shifts, thus outperforming all. SARAD also overpasses GDN, which despite its explicit spatial modeling adopts prediction errors as its sole detection criterion, limiting its performance. SARAD's top performance on SWaT and HAI underlines its ability to unravel complex spatial associations even in complex large-scale systems. Standard deviations of Table 3 are reported in Appendix K.Anomaly Diagnosis Table 4 shows the anomaly diagnosis performance in metrics defined in Section 4.1. DiffAD uses a subset of SMD features and thus is discarded from comparisons for fairness. SARAD outperforms baselines on the point-based HR@150% (SMD: +26.67% TranAD, SWaT: +5.10% USAD, HAI: +3.81% GDN) and NDCG@150% (SMD: +27.97% TranAD, SWaT: +4.76% GDN, HAI: +2.93% GDN). SARAD also outperforms on the range-based IPS@150% on most datasets (SMD: 43.19% TranAD, SWaT: 33.43% TranAD). Unlike SMD which performs forensic diagnosis to label anomalous features, SWaT and HAI label only the origins of cyberattacks as diagnosis labels. Consequentially, attack origins sometimes might not behave anomalously, e.g., attacks had failed, or the full set of anomalous features were not identified, thus diminishing the performance numbers on SWaT and HAI. SARAD generally outperforms detectors underpinned by temporal modeling due to its sensitivity to spatial associative changes. SARAD also outperforms spatial detectors such as GDN whose prediction errors limit its temporal scope to a single time point. Standard deviations of Table 4 are reported in Appendix L.Visualization Figure 3 visualizes a real-world anomaly example via SARAD. Our detector captures the significant SAR caused by the anomalous features. The loose reconstruction of the progression raises the progression-based score p and, in turn, the joint detection score s. Taking a broader view of the series in Fig. 3h , SAR significantly raises the scores at the start of the anomaly, even when the data-based errors r are small. SARAD exploits SAR to achieve more robust anomaly detection.Table 4 : Anomaly diagnosis performance. Point-based HR@P %, NDCG@P %, and range-based IPS@P % are reported. All values are average percentages from five random seeds.Method HR@P % ND@P % IPS@P % HR@P % ND@P % IPS@P % HR@P % ND@P % IPS@P % P(e) (f) (g) (h) rj pj sj (i)Figure 3 : Visualization of applying SARAD for detection on SMD. 3a shows the raw time series right before and during an anomaly p i (colored in red). An input time window for SARAD is bounded in the black box. 3b and 3c show the average association mapping ĀL via final L-th layer's MHSA. 3d shows the aggregated progression S according to Eq. 3. 3e is its reconstruction. 3f, 3g, 3h show the scores p, r, and joint s according to Eq. 5 per feature. 3i shows the anomaly scores for 3a's segment. Anomalous features (#1, #9, #10, #12, #13, #14 ,and #15) are highlighted with red bounding boxes.Complexity and Time Overheads SARAD incurs 32 mins for training and and 0.39 ms for inference per sample on HAI, the largest dataset, falling far below the data collection time and sampling frequency. Those numbers are comparable with baselines and detailed in Appendix N.Spatial Progression Reconstruction To evaluate the effectiveness of the progression module, we perform ablation studies on its submodules in Table 5 . Standard deviations are reported in Appendix K. Removing the ReLU, i.e., to capture both association increases and reductions, in progression loses the focus on asscoation reduction and impairs the detection performance. Replacing the column sum operation in aggregation with the row sum which characterizes the disconnection of anomalous features from others and is shown to be less effectiveness than the column sum representing the drop out rates. Fully concatenating without sum operation dilutes the reduction patterns and significantly hurts the detection, at a cost of complexity. For the detection submodule, using the progression directly instead of the reconstruction errors registers reductions as anomalies directly and underperforms except on SMD due to its inability to rule out normal reduction patterns. 6 compares the detection performance using Eq. 5 (Joint), using only data-based r (DR), and using only progression-based p (SPR). While the data reconstruction is a robust criterion of anomalousness, SARAD embeds the spatial information into the joint criterion and outperforms either single criterion overall. Standard deviations are reported in Appendix K.Additional ablation studies on the choice of diagnosis criterion are detailed in Appendix D. The broader impact of the work presented rests in the increasingly pervasive nature of industrial control systems, intrusion detection systems, and remote monitoring solutions in healthcare contexts, all of which commonly utilize some form of anomaly detection. Further to the time series analysis that is commonplace, the work presented in this paper demonstrates how Transformer can be used to learn the spatial associations that ubiquitously characterize these feedback-controlled systems, supplementing time series analysis to provide state-of-the-art performance. As such, the work presented has broad applicability, whilst explicitly targeting automated industrial control systems.While the model size scales linearly with the number of features, the time complexity of SARAD is quadratic with respect to the features. SARAD could incur significant training and inference overheads when the supervisory system is extensively large. We caution that the overheads of the largest dataset in our experiments fall well below data collection overhead and sampling frequency (see Appendix N). To scale, we will explore hierarchical time series anomaly detection via clustering. Another limitation of this work is the scarcity of forensically labeled datasets like SMD for anomaly diagnosis, not least due to the intensive labor and domain knowledge implied. To address that gap, we will explore publicly available audit and operational time series data for sources. (a) Raw time series before, during, and after an anomaly p i . # 1 # 9 # 1 0 # 1 2 # 1 3 # 1 4 # 1 5SWaT HAI Method Criter. HR@P % ND@P % IPS@P % HR@P % ND@P % IPS@P % HR@P % ND@P % IPS@P % Concerning the rationality of data-only diagnosis criterion in Eq. 6, we consider an alternate joint diagnosis criterion in line with the detection criterion in Eq. 5.r j = || X(j,•) -X (j,•) || 2 2 , p j = || Ŝ(•,j) -S (•,j) || 2 2 , s j = (r j -µ rj )/σ rj + (p j -µ pj )/σ pj (7) where r j is feature j's data reconstruction error, p j its progression reconstruction error, µ rj , µ pj the means of r j , p j on the validation set, and σ rj , σ pj the standard deviation of r j , p j . Table 7 compares the diagnosis performance using only r j (SARAD), using only p j (SPR), and using Eq. 7 (Joint). Unlike in anomaly detection, the great discrepancy between r j and p j more than often degrades the performance of the joint criterion. SARAD uses r j only which produces suboptimal and yet reliable performance in the longer anomalous horizons. Standard deviations are reported in Appendix L.We include four real-world datasets collected under industrial control and IT service monitoring settings for evaluations. Anomalies range from IT service outages to external cyber-physical attacks against control systems.1. Server Machine Dataset (SMD) (Su et al., 2019b,a) is a server metric dataset from a large-scale IT company. Engineers annotated anomalous events in the second half of the data with indicator-level attributions. 2. Pooled Server Metrics (PSM) dataset (Abdulaal et al., 2021a,b) captures key performance indicators of servers on an online shopping platform. Website engineers annotated anomalous events for data in the last eight weeks. 3. Secure Water Treatment (SWaT) dataset (Mathur and Tippenhauer, 2016; iTrust, 2023) contains sensor readings and actuator status on a minuscule real-world water treatment system during a six-day normal operational period. A knowledgeable attacker performed 36 cyber-physical attacks during a five-day attack period and they are labelled as anomalous accordingly. 4. HIL-based Augmented ICS (HAI) dataset (Shin et al., 2021b,a) records measurements and control actions within a Hardware-In-the-Loop (HIL) dual power (steam-turbine and hydropower) generation testbed during its two-week operation. Both single-point primitive and multi-point combined attacks are performed on the testbed to emulate a threat actor with cyber-physical capacities. We use the 21.03 version of HAI.All training sets contain only unlabeled data and the test sets contain data with anomaly labels. Statistics of the datasets are given in Table 2 .We further characterize the detection datasets by the lengths of the anomalous events. Figure 8 shows the empirical cumulative distribution function of the anomalous lengths. SWaT has the longest median length of 447 among the four datasets considered, followed by HAI (162), SMD (11), and lastly PSM (5). The very short lengths on SMD and PSM benefit temporal detectors which tend to embed a single or few time points (see Appendix M), whereas SARAD adopts a half time window embedding strategy. The catch is that SARAD can learn spatial relationships with temporal aggregated information per feature, while temporal detectors could not, bringing about benefits of performing anomaly detection in the spatial association space. A more scalable approach for temporal aggregation is to be explored in the future, though variable window sizes or subseries splits might be implied.Figure 8 shows the empirical cumulative distribution function of the diagnosis label lengths. Whereas SMD forensically labels features which deviate from their normal behavioral patterns as anomalous, SWaT and HAI only label the points of attacks as anomalous as the their creators have advanced knowledge of such attacks. The latter labeling strategy results in incomplete sets of anomalous features and diminishes the diagnosis performance of all models, as evidence in Table 4 (a) SMD.1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75 3.00 We implement SARAD in Python using pyTorch library (Paszke et al., 2019) and Hydra framework (Yadan, 2019) . All experiments are run on a single NVIDIA A10 (24GB) GPU. Adam optimizer (Kingma and Ba, 2015) is used and learning rate is halved every epoch for 3 epochs to prevent over-fitting. The time window size is 2W = 100. The data reconstruction module has H = 8 attention heads per layer with attention length D = 512 and hidden length D F F = 2048.For hyperparameter tuning, training set is temporally partitioned into 80% for training and 20% for validation. On each dataset we first perform TPE sampling (Bergstra et al., 2011) for number of encoding layers L ∈ {3, 5} and learning rate ∈ [10 -4 , 10 -2 ] to derive the best data reconstruction loss L R on the validation set. The progression module by default has hidden length of D P = 64. We then perform TPE sampling to search weight λ L S ∈ [10 -2 , 10 2 ] for the progression reconstruction loss L S on the validation set.During the review period, code is anonymized and openly available at https://github.com/ daidahao/SARAD/ with specific instructions and scripts to reproduce experimental results. All data used in our experiments can be openly accessed from public repositories or requested via original authors' websites. Full links are provided in Appendix E. We examine the hyperparameter sensitivity of SARAD's detection performance. Concretely, we consider the effects of the sliding window size (default is 100), the number of training epochs (3), the attention length of Transformer encoding layers D (512), and the hidden length of the progression reconstruction module D P (64). Figure 10 and 11 present the results. On the window size, a sliding window too small confines temporal modeling to small temporal receptive fields and contains anomaly detection in the association space. However, a sliding window too large incurs higher computational costs, although unlike temporal modeling the costs here are linear. On datasets with shorter anomalous lengths such as SMD and PSM, the anomalous patterns are diluted even further, resulting in performance degradation.On the number of training epochs, fewer epochs lead to model underfitting, and yet overfitting is largely prevented with more epochs due to the aggressive learning rate halving per epoch. On the attention length D, a larger Transformer is prone to overfitting with visible performance drops on SMD and PSM as D passes the default 512, both of whose monitored systems are smaller in scale. On the hidden length D P , a more complex progression anomaly detector does not adversely impact the performance, suggesting that the association space is less prone to detection overfitting than the data space. Conventional metrics such as precision, recall, F1 and the threshold-independent Area Under the Curve (AUC) scores are commonly point-based, i.e., predicted labels are scored individually by time points (Zhou et al., 2023; Zhang et al., 2022; Xu et al., 2018) . Real-world benchmarks are rife with range-wise anomalies spanning consecutive time points (Wagner et al., 2023) . Point-based metrics are generally ill-suited for evaluating detection performance due to the continuous-discrete conversion and the series-label misalignment (labeling anomalies precisely is hard) (Garg et al., 2022; Tatbul et al., 2018) . Here, we use the range-based metrics proposed in (Paparrizos et al., 2022) . Compared against their point-based counterparts, they provide robustness to labeling delay and scoring noises as well as performant detector separability and series consistency.Given a set of anomalous ranges P = {p i = (s i , e i )} wherein each anomaly p i starts at timestamp s i and ends at e i , we enclose each range with uniform l/2-length preceding and succeeding buffers.Given the anomaly label y t ∈ {0, 1} at each timestamp t, we derive a new soft label ỹt ∈ [0, 1] as per the minimum temporal distance of t to any anomaly p i ∈ P:EQUATIONwhere l is the buffer length, normally set to the median segment length in P. Within the buffers, ỹt monotonically increases from √ 2/2 to 1 as the distance decreases. With the new soft label series Ỹ = {ỹ t }, we define the True Positives (TP), False Positives (FP), True Negatives (TN), and False Negatives (FN) accordingly:T P = ỸT • Ŷ, F P = (1 -Ỹ) T • Ŷ, T N = (1 -Ỹ) T • (1 -Ŷ), F N = ỸT • (1 -Ŷ) (9)We then compute the threshold-independent AUC for the Receiver Operating Characteristic (AUC-ROC), i.e., TP rate vs FP rate, and the Precision-Recall (AUC-PR) curves respectively to be rid of thresholding impact. Fully parameter-free Volume Under the Surface (VUS) scores for AUC-ROC (VUS-ROC) and AUC-PR (VUS-PR) are also computed under different buffer lengths l ∈ [0, 2l].In line with previous works (Tuli et al., 2022; Zhao et al., 2020) , we use common metrics such as Hit Rate (HR) (Su et al., 2019b) and Normalized Discounted Cumulative Gain (NDCG) (Järvelin and Kekäläinen, 2002) to measure performance where diagnosis labels are available. Given a set of anomalous features G i ⊆ [N ] at an anomalous timestamp t ∈ p i ∈ P as a diagnosis ground-truth and the set of top k-ranked features Γ t @P % according to Eq. 5 where k = ⌈|G i | × P %⌉, say k = 5 when |G i | = 3 and P = 150, the HR at P % (P ≥ 100) features is the overlap ratio between the two:EQUATIONIn information retrieval, DCG measures the cumulative utility of retrieved documents by their ranking order up to a certain position. NDCG normalizes the DCG by the maximum possible DCG. They are parameterized by P % to determine the location in our evaluation and calculated as follows.DCG t @P % = k j=1 r j log 2 (j + 1),IDCG t = |Gi| j=1 1 log 2 (j + 1), NDCG t @P % = DCG t @P % IDCG t (11)where r j ∈ {0, 1} is the relevance value of the j-th element and, in this case, the membership of Γ t @P %'s j-th feature in G i . NDCG has a value strictly between 0 and 1.At the range level, we measure the Interpretation Score (IPS) initially proposed in Li et al. (2021) and here expanded to fit the P % parameterization. For each anomalous range p i , the IPS score is:EQUATIONwhere Ω i @P % is the top k-ranked features according to max t∈pi s j,t , the j-th feature's maximum anomaly score during p i and k = ⌈|G i | × P %⌉. It is the HR equivalence at the range level as per the highest anomalous scores per feature. A PR V ROC V PR A ROC A PR V ROC V PR A ROC A PR V ROC V PR A ROC A PR V ROC V PROurs -0.98 0.75 1.01 0.76 1.07 0.69 1.09 0.69 0.70 1.34 0.60 1.05 0.47 1.08 0.56 1.11 Prog. (Eq.2) no ReLU 1.00 0.45 1.00 0.45 0.76 0.61 0.73 0.61 0.52 5.17 0.51 4.87 0.19 0.33 0.17 0.33 Aggr. (Eq.3) row sum 0.87 0.61 0.88 0.62 1.30 1.30 1.29 1.28 1.39 3.62 1.40 3.61 0.22 0.67 0.21 0.60 Aggr. (Eq.3) no sum 6.34 1.21 6.37 1.20 0.71 0.21 0.69 0.21 0.48 0.98 0.48 1.00 0.51 0.42 0.56 0.39 Detection S directly 0.70 0.51 0.69 0.50 0.85 0.42 0.85 0.42 1.17 1.20 1.12 1.12 1.74 4.39 1.81 4.18Table 8 reports the standard deviations of anomaly detetion performance as reported in Table 3 in Section 4.2. 4 . Standard deviations of point-based HR@P %, NDCG@P %, and range-based IPS@P % are reported. All values are percentages.Method HR@P % ND@P % IPS@P % HR@P % ND@P % IPS@P % HR@P % ND@P % IPS@P % 7 . Standard deviations of point-based HR@P %, NDCG@P %, and range-based IPS@P % are reported. All values are percentages.PMethod HR@P % ND@P % IPS@P % HR@P % ND@P % IPS@P % HR@P % ND@P % IPS@P % 0.00 0.00 0.00 0.44 0.43 0.00 0.00 0.00 1.95 1.01 0.00 0.00 0.00 5.39 1.35 3.33 3.13 3.18 5.81 1.59 DCdetector 0.00 0.00 0.00 0.78 0.13 0.00 0.00 0.00 0.32 0.15 0.00 0.00 0.00 0.04 0.21 N/A N/A N/A N/A N/A USAD 2.10 1.91 1.17 2.68 1.20 0.44 0.43 0.50 0.61 0.27 3.30 1.41 0.70 2.64 1.25 13.69 6.17 4.12 4.73 PIn addition to the model evaluations using range-based metrics such as VUS-ROC and VUS-PR in Sections 4.2 and 4.3, we conduct point-based evaluations herein using point-based metrics exclusively. Table 13 reports report the Precision, Recall, and F1 scores under the threshold where the method of interest achieves the best F1 score. To mitigate the impact of thresholding protocol, Table 13 also reports the threshold-independent AUC-ROC and AUC-PR scores. Table 14 reports the standard deviations.SARAD achieves state-of-the-art performance on most datasets except PSM. The enlarged temporal receptive field of SARAD (half input window) contributes to its underperformance by point-based metrics, especially on datasets where anomalous ranges are short (see Table 2 ), when compared against most others' receptive filed of single or few time points (typical 1D Conv kernel size is 3). We again underline that most real-world anomalies are continuous and point-based metrics are mismatched for such anomalies (see Appendix I for discussion). = 708K T = 132K T = 497K T = 922K N = 38 N = 25 N = 51 N = 79 T = 1 min T = 1 min T = 1 s T =We trained all baselines using official implementations where available and recommended hyperparameters from respectively papers are used. Some baselines, such as DiffAD and AT, have dataset-specific hyperparameters and here we adopted them as well. The open-accessed URLs of the baselines used are listed as followed.• IF (ICDM'08) (Liu et al., 2008) : https://github.com/xuhongzuo/deep-iforest.• DIF (TKDE'23) (Xu et al., 2023) : https://github.com/xuhongzuo/deep-iforest.• TranAD (VLDB'22) (Tuli et al., 2022) : https://github.com/imperial-qore/ TranAD. • ATF-UAD (NN'23) (Fan et al., 2023) : https://github.com/wzhSteve/ATF-UAD.• AT (ICLR'22) (Xu et al., 2022) : https://github.com/thuml/Anomaly-Transformer.• DCdetector (KDD'23) (Yang et al., 2023) : https://github.com/DAMO-DI-ML/ KDD2023-DCdetector. • USAD (KDD'20) (Audibert et al., 2020) : https://github.com/manigalati/usad.• GDN (AAAI'21) (Deng and Hooi, 2021) : https://github.com/d-ailin/GDN.• MAD-GAN (ICANN'19) (Li et al., 2019) : https://github.com/LiDan456/MAD-GANs.Official implementation was migrated to pyTorch for uniform environmental set-up. See our codebase for details. • DiffAD (KDD'23) (Xiao et al., 2023) : https://github.com/ChunjingXiao/DiffAD. Figure 12 visualizes the two components of the joint detection criterion s in Eq. 5, i.e., the data reconstruction error r and the progression reconstruction error p. Balanced resampling is applied here. Recall that s is the sum of normalized r and p. Most anomalous samples (input series) either has high r or high p, and oftentimes both. The former measures the magnitude of anomalousness in the data space, the latter in the spatial association space. The basis underpins the formalization of the joint detection criterion. Figure 13 the distributions of the joint anomaly scores. On SMD and PSM, the scores of the normal and anomalous samples (input series) overlap more heavily than on SWaT and HAI. The observations here highlight the difficulty of anomaly detection in the service monitoring space, more so than in industrial control where measurements and actuation result from well-defined control logics. The hardness is evident is Table 3 of detection performance and 4 of diagnosis performance. TPR↑ (@FPR) FPR↓ (@TPR) TPR↑ (@FPR) FPR↓ (@TPR) TPR↑ (@FPR) FPR↓ (@TPR) TPR↑ (@FPR) FPR↓ (@TPR) TPR↑ (@FPR) FPR↓ (@TPR) TPR↑ (@FPR) FPR↓ (@TPR) TPR↑ (@FPR) FPR↓ (@TPR) TPR↑ (@FPR) FPR↓ (@TPR) Effectiveness of each anomaly detector is influenced by the selection of thresholds. Previously, we reported the range-based and point-based threshold-independent performance metrics. Here, we study the influence of threshold selection on the detection performance. Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: We disclose information on compute resources for experiments in Appendix F.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Answer: [Yes] Justification: The research conform with the NeurIPS Code of Ethics.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [Yes] Justification: We discuss the broader impacts in Appendix A.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: The paper poses no risks of misuse.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: We cite and properly credit all datasets, baselines, and libraries used in this paper (see Appendices E and O for details).Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. 
 
 Appendix:  


Decision: The authors proposed algorithm, SARAD, employs a transformer-based data reconstruction approach to capture inter-feature associations. By analyzing changes in these associations over time, the algorithm identifies anomalies. The strength of the paper lies in novel formulation, new anomaly score, and a Detailed experimental analysis. The weakness lies in lack of motivation, lack of clarity about contribution, and missing benchmarks. Most of them have been addressed in rebuttal. All reviewers are consensus about the acceptance of the paper after rebuttal.  Congratulations to the authors! I am requesting the author to carefully incorporate the suggestion provided by the review in the camera-ready version.

Reviews: review: paper_summary: paper_summary: Paper proposes a deep learning based anomaly detection method for multi-variate time series data. The proposed method has two components in a single neural network. First component is a traditional auto encoder approach in which the input time series is broken into two parts (along time) and the model is trained to use the first part to generate the second part, and identify anomalies as a difference between the predicted and true output. The second component trains a transformer to generate the associations across the different features (akin to cross-correlation) and then train an MLP to predict the associations and identify anomalies if the predict associations are different than the observed associations. Thus, the method captures a broader range of temporal anomalies. 

Experimental results on a variety of benchmark data sets and comparisons with state of art methods demonstrate that the proposed method (SARAD) is able to identify anomalies that might not be apparent to approaches that do not consider the associations.
strength_weakness: strength_weakness: A terminology issue - spatial often refers to data with positional coordinates (geographic or xy on a domain) which indicates spatial relationships between observations. On top of that, each observation would be represented as a multi-dimensional feature vector. In this paper, the authors use the term spatial to refer to the multi-dimensional feature vector representation of data. I found that very confusing. I think it might be better to fix the terminology to improve readability of the paper.
questions: questions: Are there any types of anomalies that would still not be captured by this method?

Similar to above, what kind of false positives could be captured by adding the association part?
limitations: limitations: The paper has addressed the limitations - I have suggested a couple above. Paper does not have a direct societal impact so the authors have a reasonable response to that.
paper_summary: paper_summary: In this paper, the authors focus on incorporating spatial information for time series anomaly detection and diagnosis. The proposed algorithm, SARAD, employs a transformer-based data reconstruction approach to capture inter-feature associations. By analyzing changes in these associations over time, the algorithm identifies anomalies under the assumption that anomalous features cause a reduction in the perceived association.
strength_weakness: strength_weakness: The related work section lacks a detailed discussion of all the baseline methods used for comparison in result analysis. For a more comprehensive understanding and to highlight the contributions of the paperwork, it would be better to incorporate those methods and discuss their limitations.

The third paragraph in the Introduction section is not easy to follow.  Providing more background information and reorganizing would be recommended.
questions: questions: I'm interested in understanding how the SARAD method distinguishes reductions in spatial associations that indicate anomalies from those that could be considered normal variations. Please clarify the types of anomalies the SAR phenomenon is effective for.

The authors could have enhanced their study by comparing SARAD with existing methods (e.g., InterFusion), which integrate both spatial and temporal information for anomaly detection
limitations: limitations: yes
paper_summary: paper_summary: In this paper, the authors aim to address the problem of anomaly detection and diagnosis for time series data. The authors consider that the existing methods may obscure or dilute the spatial information and the interaction between different variates, they propose the SARAD model to capture these interactions with the transformer. The authors evaluate the proposed method on several datasets and achieve good performance.
strength_weakness: strength_weakness: However, there are some problems to solve  
1.	The authors claim that existing methods may assume feature independence or combine variables of a diverse physical nature. Indeed, the recent methods for time series analysis usually employ independent channel assumptions. However, several methods like the conventional RNN and the graph-based methods do not assume that the features are independent. Therefore, the motivation of this paper might not be convincing.  
2.	Moreover, Figure 1 might be not clear to explain the motivation of the paper. Specifically, Figure 1(a) is not informative and Figure 1(b)(c)(d) are almost the same.  
3.	The contribution of the proposed method is limited. There are other methods that use transformers as the backbone networks like [1]. It is suggested that the authors should provide a detailed discussion and compare with it.  

[1] Xu, Jiehui, et al. "Anomaly transformer: Time series anomaly detection with association discrepancy." arXiv preprint arXiv:2110.02642 (2021).
questions: questions: Please refer to weaknesses
limitations: limitations: Please refer to weaknesses
paper_summary: paper_summary: This manuscript proposes a spatial association reduction method for anomaly detection and diagnosis (SARAD) where anomalous features have low association values in reconstructing the multivariate time-series. The association values are derived from an attention between the features and exhibit changes with time which enables the detection of anomalous features and the time period of anomalies. SARAD is built on two models: a transformer for data reconstruction and an MLP for the reconstruction of spatial progression. The last anomaly score is a combination of both reconstruction scores.
strength_weakness: strength_weakness: - In my opinion a main issue regarding the anomaly diagnosis should be addressed. For example, in Fig. 1 there are three features (9, 13 and 16) correlated with 12 and 15. An additional example is from Fig. 5 where features 13 and 16 are highly correlated with 12 and 15. What distinguishes the detected features from their correlated features?
- There was no discussion about the identity shortcut problem. Such a problem is typical for reconstruction-based methods where the model learns to reconstruct anomalies and thus it becomes harder to separate them from the normal ones. Maybe I missed it in the manuscript, but can you elaborate on this issue? How do you address this issue in your method?
- The Point-based evaluation on PSM dataset indicates some weaknesses regarding the sensitivity to the window size and short-term anomalies (see Table 13). The clustering for SMD and PSM is also not of the same quality compared to other datasets (see Figs. 12 and 13).
- The model size/parameters of the proposed method exceed other baselines by big margines (except for DiffAD).
questions: questions: - How is the threshold defined from the anomaly scores to detect anomalies? 
- Can you please provide a brief explanation why progression-based $p$ (SPR) is less effective for SWaT?

Minor:
- I find it a bit confusing to call it spatial association for time-series processing. Maybe it is better to call it channel/feature/variable.
- It is clearer if you visualize or indicate the ground truth for anomaly diagnosis alongside the detected features i.e., in Figs. 1, 3, and 4-7.
- Lines 84-85: Can you please make it clearer which spaces you are refereeing to? I think you mean data space and spatial progression space.
limitations: limitations: I think the issues about anomaly diagnosis and the identity shortcut should be discussed (see weaknesses).

Comments: title: Response to W1 & W4
Thank you
comment: Thanks for your responses. The clarifications are much appreciated. I stand by my current review rating.
**W1:**
You are correct in that there are non-anomalous features (9, 13, 16) correlating with anomalous features (12, 15) in Fig. 1. Similar correlated non-anomalous features can be found in Fig. 5. What distinguishes them from the anomalous features is the frequency of the patterns they showed during the anomaly period in the whole training set. Such patterns are frequent in the training set for features 9, 13, and 16, but not for features 12 and 15 (**see the following MSE experiment** for measuring frequencies). As the Transformer model was trained on the training set and learned the distributions of normal feature patterns, it does not register such frequent patterns of those correlated features as showing association reductions.

**MSE Experiment:**
We herein illustrate the frequencies of correlated non-anomalous features (9, 13, 16) in Fig. 1 relative to the anomalous (12, 15) using **Mean Squared Errors (MSEs)**. We first extract a subseries (of length 255) for each feature of interest during the anomaly period in Fig. 1 and refer to them as templates. We then compute the MSE between each template and each equal-length subseries from the training set (of size 566K) of the same feature. All features have already been normalized in data pre-processing. The table below shows the densities of MSE values for each feature on ranges starting from 0. It demonstrates a stark contrast between the anomalous features and the correlated non-anomalous features, with the majority of MSE values falling below 1 for the correlated and yet 0.00% for the anomalous features. The contrast remains stark as we shrink the ranges when plotting a histogram (which cannot be shown here due to OpenReview’s restrictions), with a dense concentration of MSEs near 0 for only the correlated features. 

|||||**MSE Range**||||
|--------------|-----------|----------|----------|-------------|----------|----------|------------|
||**Feature**|**(0, 1]**|**(1, 2]**|**(2, 3]**|**(3, 4]**|**(4, 5]**|**(5, $+\infty$)**|
|**Correlated**|9|**76.62%**|8.74%|4.60%|3.53%|2.46%|4.06%|
||13|**74.87%**|14.03%|3.58%|1.98%|1.01%|4.53%|
||16|**55.83%**|21.00%|7.60%|5.99%|2.71%|6.87%|
|**Anomalous**|12|0.00%|1.43%|**90.45%**|3.66%|0.65%|3.81%|
||15|0.00%|**62.69%**|14.00%|11.88%|5.06%|6.36%|

The MSE distributions observed here (and similarly for the Fig. 5 case) underline the usualness of correlated features’ patterns. Since other subseries in the training set share similar patterns for the correlated features, the Transformer model does not register those features as showing association reductions.

**W2:** 
The large model size generally does not give our model an advantage in performance. We refer to our Hyperparameter Sensitivity experiments in Appendix H. Figures 10(c) and 11(c) both show that on most datasets (except SWaT) the performance is insensitive to exceedingly small attention length $D$ (also denoted as $d_{model}$ in some Transformer literature). The model size is largely decided by $D$ since the majority of parameters reside in linear transformations inside the Transformer's self-attention modules, whose weights are of shape $D \times D$. We chose 512 as the default value for $D$, following the original Transformer paper and subsequent works. Compared to the default value of 512, a small $D$ of 64 does not limit the model's performance on most datasets while scaling down the model by roughly 1/64, as shown in Figures 10(c) and 11(c). The following table reports the numbers of model parameters when $D=64$ and relative scales to the defaults when $D=512$ (as in Table 15). The model sizes are now comparable to most baselines, while the performances are of similar levels (except SWaT) as when $D=512$ and thus remain state-of-the-art.


|                   | SMD    |       | PSM    |       | SWaT   |       | HAI    |       |
| ----------------- | ------ | ----- | ------ | ----- | ------ | ----- | ------ | ----- |
| Method            | Param. | Scale | Param. | Scale | Param. | Scale | Param. | Scale |
| **Ours** ($D=64$) | 198K   | 1/48  | 284K   | 1/56* | 212K   | 1/45  | 343K   | 1/46  |


One possible explanation for under-performance on SWaT is that SWaT has many long-range anomalies, which could require a large $D$ value to adequately encode rich spatial information during anomalies. Figures 10(c) and 11(c) show that $D=256$ (Param. count=2.46M, Scale=1/4) results in a much closer performance to the default $D=512$.

------
*Correction: It came to our attention when computing the model scales that Table 15 in the original version incorrectly reported our parameter count on PSM. The count should be 15.85M, not 1.59M. Apologies for the typo. All other parameter counts for our model are correct. We will correct the count in the final version.
We thank you for replying all questions raised by reviewer.
Thanks for the detailed response. The response has addressed my concerns. I have raised the score.
Thank you for the response. The authors have addressed all the issues menitoned in the review.
Thank you for your response. There are still two unanswered questions. I will try to rephrase my questions to make them clearer:

**W1.** For example, in Fig. 5, the spatial association reduction happens for two features 12 and 15. The question is: Why didn’t this reduction happen for features 13 and 16 as well even though they are highly correlated with features 12 and 15?

**W4.** My comment was about the model capacity rather than the training/inference time. As seen in Table 15, page 29, the proposed model exceeds other baselines by big margines (except for DiffAD). I think this gives the proposed model an advantage and doesn’t fully guarantee a fair comparison.