Spatio-Temporal Anomaly Detection with Graph Networks for Data Quality Monitoring of the Hadron Calorimeter

Abstract

The Compact Muon Solenoid (CMS) experiment is a general-purpose detector for high-energy collision at the Large Hadron Collider (LHC) at CERN. It employs an online data quality monitoring (DQM) system to promptly spot and diagnose particle data acquisition problems to avoid data quality loss. In this study, we present a semi-supervised spatio-temporal anomaly detection (AD) monitoring system for the physics particle reading channels of the Hadron Calorimeter (HCAL) of the CMS using three-dimensional digi-occupancy map data of the DQM. We propose the GraphSTAD system, which employs convolutional and graph neural networks to learn local spatial characteristics induced by particles traversing the detector and the global behavior owing to shared backend circuit connections and housing boxes of the channels, respectively. Recurrent neural networks capture the temporal evolution of the extracted spatial features. We validate the accuracy of the proposed AD system in capturing diverse channel fault types using the LHC collision data sets. The GraphSTAD system achieves production-level accuracy and is being integrated into the CMS core production system for real-time monitoring of the HCAL. We provide a quantitative performance comparison with alternative benchmark models to demonstrate the promising leverage of the presented system.

1. Introduction

Deep learning (DL) has become increasingly prevalent for anomaly detection (AD) applications for reliability, safety, and health monitoring in several domains with the proliferation of sensor data in recent years [1,2,3]. AD has been applied for a diverse set of tasks, including but not limited to machinery fault diagnosis and prognosis [4,5], electronic device fault diagnosis [6,7,8,9], medical diagnosis [10,11,12,13], cybersecurity [14,15,16], crowd monitoring [17,18,19,20,21,22,23], traffic monitoring [24,25], environment monitoring [26], the Internet of things [3,27], and energy and power management [28,29]. AD aims to determine anomalies depending on the setting and application domain [2]. An anomaly is generally an odd observation—abnormalities, deviants, outliers, discords, failures, intrusions, exceptions, aberrations, peculiarities, or contaminants—from a bulk of observations often indicating peculiar underlying incidents [1]. AD methods can be categorized as supervised or unsupervised approaches: (1) supervised approaches require annotated ground-truth anomaly observations, and (2) unsupervised methods do not require labeled anomaly data and are more generally pragmatic in many real-world application settings, as data annotation is expensive. Unsupervised AD models trained with only healthy observations are often categorized as semi-supervised approaches.

Deep learning has become effective for AD modeling because of its capability to capture complex structures, extract end-to-end automatic features, and scale for large data sets [1,2]. Several DL models have been proposed in the literature for diverse data types, such as structural [1], time series [7,8,9,12,13,16,27,29,30,31,32,33,34,35,36,37,38], image [10,26], graph network data [14,15,24,25,39], and spatio-temporal [10,14,15,17,18,19,20,21,22,24,25,39]. Spatio-temporal (ST) data are commonly collected in diverse domains, such as visual streaming data [17,18,19,20,21,22,23], transportation traffic flows [24,25], sensor networks [14,15,39], geoscience [26], medical diagnosis [10], and high-energy physics [40,41]. A unique quality of ST data that differentiates it from other classic data is the presence of dependencies among measurements induced by the spatial and temporal attributes, where data correlations are more complex to capture by conventional techniques [42]. ST anomaly is thus defined as a data point or cluster of data points that violate the nominal ST correlation structure of the healthy ST data [10,14,15,17,18,19,20,21,22,24,25,39]. The wide range of unsupervised DL AD methods discover anomalies in temporal context using density clustering on latent space [9], data reconstruction [9,13,30], and prediction [16,27,30,33,34]. Variants of recurrent neural networks (RNNs) [7,8,9,13,21,22,24,27,34,35], convolutional neural networks (CNNs) [9,18,19,20,21,27,30,33,34], generative adversarial networks (GANs) [12,29,35,36], graph neural networks (GNNs) [24,25,30,37,38], and transformers [37] have been explored and achieved competitive performance for multivariate temporal or ST AD.

The Large Hadron Collider (LHC) is the largest particle collider ever built globally. It is designed to conduct experiments in physics and increase our understanding of the universe, expecting that new findings will lead to practical applications. The LHC is a two-ring superconducting hadron accelerator and collider capable of accelerating and colliding beams of protons and heavy ions with the unprecedented luminosity of ${10}^{34}$ cm${}^{-2}$s${}^{-1}$ and ${10}^{27}$ cm${}^{-2}$s${}^{-1}$, respectively, at a velocity close to the speed of light—$3\times {10}^8$ ms${}^{-1}$ [43,44]. The Compact Muon Solenoid (CMS) experiment is a general-purpose detector for high-energy physics (HEP) at the LHC [40]. The CMS employs a data quality monitoring (DQM) system to guarantee high-quality physics data through online monitoring that provides live feedback during data acquisition and offline monitoring that certifies the data quality after offline processing [45]. The online DQM identifies emerging problems using a reference distribution and predefined tests to detect known failure modes using summary histograms, such as a digi-occupancy map of the CMS calorimeters [46,47]. A digi-occupancy map contains a histogram record of particle hits of the data-recording channels of the calorimeters. The calorimeters could have several flaws, such as issues with the front-end particle sensing scintillators, digitization and communication systems, backend hardware, and algorithms, which are usually reflected in the digi-occupancy map. The growing complexity of the detector and the physics experiments make data-driven AD systems essential tools for the CMS to identify and localize detector anomalies automatically. The CMS detector consists of a tracker to reconstruct particle paths accurately, two calorimeters—the electromagnetic (ECAL) and the hadronic (HCAL) calorimeters to detect electrons, photons, and hadrons, respectively—and several muon detectors. The synergy in AD has thus far achieved promising results on spatial 2D histogram maps of the DQM for the ECAL [48] and the muon detectors [49].

Previous studies only considered extreme anomalies, such as no reading, dead, and high-noise, hot-particle-sensing calorimeter channels. Detecting degrading channels is essential for quality deterioration monitoring and early intervention, but they are often challenging to capture; for instance, the improperly tuned bias voltage on the HCAL physics-particle-sensing channels caused nonuniformity in the hit map of the DQM, but the channels were neither dead nor hot [50]. The calorimeter channels may degrade with a subtle abnormality before reaching extreme channel fault status. Capturing such subtle anomalies, e.g., a slow system degradation, makes temporal AD models appealing for early anomaly prediction before ultimate system failure. Time-aware models extract temporal context to enhance AD performance. A few efforts have thus far been focused on temporal models despite the acknowledged potential in the future automation technology challenges at the LHC [7,48]. Our study focuses on DQM automation through time-aware AD modeling using digi-occupancy histogram maps of the HCAL. The digi-occupancy data of the HCAL are 3D due to its depthwise calorimeter segmentation. It poses multidimensional challenges, and it is relatively unexplored in ML endeavors. The particle hit map data of the HCAL are highly dependent on the collision luminosity—a measure of how many collisions are happening in a particle accelerator—and the number of particles traversing the calorimeter. The effort on data normalization that enhances the learning generalization of machine learning models is still limited.

In this study, we address the above gaps while investigating the performance of temporal DL models in enhancing AD for the HCAL DQM system. We propose to detect anomalies of the HCAL particle-sensing channels through a semi-supervised AD system—GraphSTAD—from spatial digi-occupancy maps of the DQM. Anomalies can be unpredictable and come in different patterns of severity, shape, and size, often limiting the availability of labeled anomaly data covering all possible faults. We employ a semi-supervised approach for the AD system; the concept for the AD is that an autoencoder (AE) trained to reconstruct healthy digi-occupancy maps would adequately reconstruct the healthy maps, whereas it would yield a high reconstruction error for maps with anomalies. Since abnormal events can have a spatial appearance and temporal context, we combine both the spatial and temporal features—spatio-temporal—for the AD [14,17,18,19,20,21,22,23,24,26]. The spatial nature of the digi-occupancy map of the HCAL may exhibit irregularity; although adjacent channels with the Euclidean distance are exposed to collision article hits around their region, the channels may belong to different backend circuits, resulting in a non-Euclidean spatial behavior on the digi measurements. The GraphSTAD system captures the behavior of channels from regional collision particle hits, and electrical and environmental characteristics due to the shared backend circuit of the channels to effectively detect the degradation of faulty channels. The AD system attains these utilities using a deep AE model that learns the local spatial behavior, the physical-connectivity-induced shared behavior, and the temporal behavior through convolutional neural networks (CNNs), graph neural networks (GNNs), and recurrent neural networks (RNNs), respectively.

We evaluate our proposed AD approach in detecting spatial faults and temporal discords on digi-occupancy maps of the HCAL. We simulate different realistic types of anomalies—dead channels without registered hits and hot channels dominated by electronic noise—resulting in a much higher hit count than expected, and degraded channels with deteriorated particle detection efficiency, resulting in lower hit counts than expected, to analyze the effectiveness of the AD model. The results demonstrate promising performance in detecting and localizing the anomalies. We further validate the efficacy in detecting real anomalies and discuss comparisons to benchmark models and the existing DQM system.

We briefly describe the DQM and HCAL systems in Section 2, and our data sets in Section 3. Section 4 explains the methodology of the proposed GraphSTAD model, and Section 5 presents the performance evaluation and result discussion. Finally, we summarize the contribution of our study in Section 6.

2. Background

This section describes the DQM and HCAL systems of the CMS experiment.

2.1. Data Quality Monitoring of the CMS Experiment

The detector and collision data’s offline processing complexity requires continuous data quality monitoring. Shifters and physicists at the CMS monitor the collision quality and select data usable for analysis; they look for unexpected issues that could affect the data quality, e.g., noise spikes, dead areas of the detector, and calibration problems [51]. The DQM provides feedback on detector performance and data reconstruction; it generates a list of certified data for physics analyses—the “Golden JSON” [45]. The DQM employs online and offline monitoring mechanisms: (1) the online monitoring is a real-time DQM during data acquisition, and (2) the offline monitoring—after 48 h since the collisions were recorded—provides the final fine-grained data quality analysis for data certification. The online DQM populates a set of histogram-based maps on a selection of events and provides summary plots with alarms that DQM experts inspect to spot problems. The digi-occupancy maps—one of the maps generated by the online DQM—incorporate particle hit histogram records of the particle readout channel sensor of the calorimeters. A digi—also called hit—is a reconstructed and calibrated collision physics signal of the calorimeter. Various faults in the calorimeter affecting the front-end hardware and software components appear in the digi-occupancy map. Previous efforts by [45,48,49,52] demonstrate the promising AD efficacy of using digi-occupancy maps for calorimeter channel monitoring using machine learning. However, end-to-end DL with temporal models is relatively unexplored [48,49].

The purpose of leveraging the DQM through machine learning is to address particular challenges: (1) the latency of human intervention and thresholds require sufficient statistics; (2) the volume of data a human can process in a finite time is limited; (3) rule-based approaches do not scale and assume limited potential failure scenarios; (4) dynamic running conditions change reference samples; (5) the effort to train human shifters who monitor DQM dashboards and maintain instructions is expensive. Developing machine learning models for the DQM comes with some impediments despite the potential promises; data normalization to handle variation in experimental settings, the granularity of the failures to spot, and limited availability of the ground-truth labels are among the challenges [49].

We extend the efforts in AD with ST modeling of the digi-occupancy maps of the DQM for the HCAL. Several promising ST AD models have been proposed in the literature in diverse domains [10,14,15,17,18,19,20,21,22,23,24,25,26,39]. The previous AD studies on video data sets [18,19,20,21] focus on CNNs for regular spatial feature extraction, and GNNs are gaining popularity for sensor and traffic flow data [24,25] that exhibit irregular spatial attributes with a non-Euclidean distance among nodes. GNNs have recently achieved promising results at the LHC [41,53] and outperformed CNN in learning irregular calorimeter geometry [54] and in pileup mitigation [55]. The spatial characteristics of the HCAL channels exhibit a regular spatial positioning of particle hits in the calorimeter and an irregularity in measurement due to adjacent channels may share different backend circuits. Our proposed study presents an AD model for the DQM by integrating both CNNs and GNNs [56,57] to capture Euclidean and non-Euclidean spatial characteristics, respectively, and an RNN for temporal learning.

2.2. Readout Boxes of the Hadron Calorimeter

The HCAL is a specialized calorimeter to capture hadronic particles. The calorimeter is composed of multiple subsystems such as HCAL Endcap (HE), HCAL Barrel (HB), HCAL Forward (HF), and HCAL Outer (HO) (see Figure 1).

The HCAL subsystems are made of readout boxes (RBXes) to house the data acquisition electronics. The RBXes provide high voltage, low voltage, backplane communications, and cooling to the data acquisition electronics. The HE—the use-case of our study—consists of 36 RBXes arranged on the plus and minus hemispheres of the CMS. Its front-end particle detection system is built on brass and plastic scintillators, and the produced photon is transmitted via the wavelength-shifting fibers to silicon photomultipliers (SiPMs) (see Figure 2). Each RBX comprises 4 readout modules (RMs) for signal digitization [59]; each RM has 48 SiPMs and 4 readout cards, each including 12 charge-integrating and -encoding channels (QIE11 ASICs) connected to corresponding SiPMs and a field-programmable gate array (Microsemi Igloo2 FPGA). A QIE integrates the charge from a SiPM at 40 MHz, and the FPGA serializes and encodes the data from 12 QIE channels (see Figure 2). The encoded data are optically transmitted to the backend system via the CERN versatile twin transmitter (VTTx) at 4.8 Gbps. The HE system has 17 detector scintillator layers that are read out in seven groups—hereafter referred to as $depths$; the light from the scintillators in any given group is optically added together by sending it to a single SiPM. Additional channels enable a more refined depth segmentation, ideal for precisely calibrating the depth-dependent radiation damage on the HCAL [45].

3. Data Set Description

We employed digi-occupancy map data of the online DQM system to train and validate the proposed AD system. The collision data of the LHC are aggregated into runs, each containing thousands of lumisections. A lumisection (LS) corresponds to approximately 23 s of data acquisition and comprises hundreds or thousands of collision events containing particle hit records. The digi-occupancy maps generated by the online DQM contain particle hit histogram records of the particle readout channel sensor of the calorimeters. Several faults in the calorimeter affecting the front-end particle-sensing scintillators, the digitalization and communication systems, the backend hardware, and the algorithms usually appear in the digi-occupancy map. The value of the digi-occupancy varies with the received luminosity recorded by the CMS—hereafter referred to as the luminosity—and the number of events—particles traversing the calorimeter. The maps from a sequence of LSs constitute the attribution of ST data with correlated spatial and temporal relations [42].

The digi-occupancy map root-file data sets were collected in 2018 during the LHC RUN-2 collision by the CMS experiment. The data set, from the CMS ZeroBias Primary Dataset, contains approximately 20,000 LSs from 20 different healthy runs prescrutinized by the CMS certifiers and declared in the “Golden JSON” of the DQM as of good quality [60]. The digi-occupancy map data of the HCAL have 3D spatial dimensions with $\eta$ $\varphi$, and $depth$ axes and contain digi histogram records of the physics readout channel sensor of the calorimeter referenced by $i\eta =[-32,32]$, $i\varphi =[1,72]$, and $depth=[1,7]$ axes (see Figure 3). The maps—one per LS—were populated with the per-LS received luminosity up to 0.4 ${\mathrm{pb}}^{-1}$ and the number of events up to 2250. Our working data set contains about 20,000 map samples, each with a dimension of $[i\eta =64\times i\varphi =72\times depth=7])$.

4. Methodology

This section presents the proposed GraphSTAD approach for HCAL monitoring using digi-occupancy maps.

There is a lack of adequate labeled anomaly data covering all possible channel fault scenarios for the HCAL, and the anomalies may follow unpredictable patterns with different severity, shape, and size. We thus employed a semi-supervised approach for the AD system—GraphSTAD system; we trained a deep AE model to reconstruct healthy digi-occupancy maps with low contamination of anomalies. We present an ST reconstruction AE to detect abnormality in the HCAL channels using reconstruction deviation scores on ST digi-occupancy maps from consecutive lumisections (see Figure 4). The AE combines CNNs, GNNs, and RNNs to capture ST characteristics of digi-occupancy maps. The spatial feature extraction of the CNNs is leveraged with GNNs to learn circuit and housing-connectivity-induced spatial behavior irregularities among the HCAL sensor channels. There are approximately 7000 channels—pixels—on the digi-occupancy map of the HCAL endcap subsystem, housed in 36 RBXes. The channels in a given RBX are susceptible to system faults in the RBX due to the shared backbone circuit and environmental factors like temperature and humidity. Our proposed GraphSTAD employs GNNs in its spatial feature extraction network pipeline to capture the characteristics of the HCAL channels owing to their shared physical connectivity to a given RBX. GNNs have recently achieved promising results in several applications at the LHC [41,53] and outperformed CNNs in learning irregular calorimeter geometry [54] and in pileup mitigation [55]. The GraphSTAD system exploits both CNNs and GNNs [56,57] to capture Euclidean and non-Euclidean spatial characteristics of the HCAL channels, respectively.

4.1. Data Preprocessing

This section explains the data preprocessing stages of our proposed AD approach: (1) digi-occupancy renormalization, and (2) graph adjacency matrix generation.

4.1.1. Digi-Occupancy Map Renormalization

The digi-occupancy ($\gamma$) map data of the HCAL vary with the received luminosity ($\beta$) and the number of events ($\xi$) (see Figure 5). We devised a renormalization of $\gamma$ through a regression model $\mathcal{R}$ to have a consistent quantity interpretation of $\gamma$ and build an AD model that robustly generalizes previously unseen run settings—$\beta$ and $\xi$ variations. The model $\mathcal{R}$ estimates the renormalizing ${\overline{\gamma }}_s$ at the sth LS using $\beta$ and $\xi$ as:

$${\overline{\gamma }}_s=\mathcal{R}(\xi ,\beta )$$

(1)

The model $\mathcal{R}$ is trained to minimize the MSE cost function, $\mathbb{E}\left[{({\gamma }_s-{\overline{\gamma }}_s)}^2\right]$, where ${\gamma }_s$ is calculated as:

$${\gamma }_s=\sum _{\forall i}^{}\gamma (s,i)$$

(2)

where $\gamma (s,i)$ is the digi-occupancy of the ith channel in the map at the sth LS. Finally, the per-channel $\gamma (s,i)$ is renormalized by its corresponding ${\overline{\gamma }}_s$ as:

$$\widehat{\gamma }(s,i)=\frac{K\gamma (s,i)}{{\overline{\gamma }}_s}$$

(3)

where $\widehat{\gamma }$ is the renormalized $\gamma$, and K is a scaling factor to compensate for the difference in the number of channels on the depth axes.

We employ fully connected (FC) neural networks to build the regression model to effectively capture the nonlinear relationships illustrated in Figure 5:

$$input(\xi ,\beta )\to \mathrm{ReLU}\left(\mathrm{FC}\right(64\left)\right)\to \mathrm{ReLU}\left(\mathrm{FC}\right(64\left)\right)\to \mathrm{ReLU}\left(\mathrm{FC}\right(7\left)\right)\to output\left({\overline{\gamma }}_s\right)$$

(4)

Figure 6 depicts the data distribution of ${\gamma }_s$ before and after renormalization with $\mathcal{R}$. The renormalization has successfully handled the discrepancies on the ${\gamma }_s$ from several runs; it overlaps and centers distributions of ${\widehat{\gamma }}_s$ and minimizes the outliers.

4.1.2. Adjacency Matrix Generation for Graph Network

We employed an undirected graph network $\mathcal{G}(\mathcal{V},\Theta )$ to represent the calorimeter channels in a graph network based on their connection to a shared RBX system. The graph $\mathcal{G}$ contained nodes $\upsilon \in \mathcal{V}$, with edges $({\upsilon }_i,{\upsilon }_j)\in \Theta$ in a binary adjacency matrix $\mathcal{A}\in {\mathbb{R}}^{M\times M}$, where M is the number of channel nodes. An edge indicated the channels sharing the same RBX as:

$$A({\upsilon }_i,{\upsilon }_j)=\left\{\begin{array}{cc}1,\hfill & \mathrm{if}\Omega \left({\upsilon }_i\right)=\Omega \left({\upsilon }_j\right)\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(5)

where $\Omega \left(\upsilon \right)$ returns the RBX ID of the channel $\upsilon$.

There are about 7K channels in a graph representation of the digi-occupancy map of the HE calorimeter—each RBX network contains roughly 190 nodes. We retrieved the channel to RBX mapping from the calorimeter segmentation map of the HE.

4.2. Anomaly Detection Modeling with Autoencoder Model

We denote the AE model of the GraphSTAD system as $\mathcal{F}$. The ST data, $\mathrm{X}\in {\mathbb{R}}^{T\times N_{i\eta }\times N_{i\varphi }\times N_d\times N_f}$, are represented as a sequence in a time window $t_x\in [t-T,t]$, where $N_{i\eta }\times N_{i\varphi }\times N_d$ is the spatial dimension corresponding to the $i\eta$, $i\varphi$, and $depth$ axes, respectively, and $N_f=1$ is the number of input variables—only a digi-occupancy quantity in the spatial data. ${\mathcal{F}}_{\theta ,\omega }:\mathrm{X}\to \overline{\mathrm{X}}$ is parameterized by $\theta$ and $\omega$ and attempts to reconstruct the input ST data $\mathrm{X}$ and outputs $\overline{\mathrm{X}}$. The encoder network of the model $E_{\theta }:\mathrm{X}\to \mathbf{z}$ provides the low-dimension latent space $\mathbf{z}=E_{\theta }\left(\mathrm{X}\right)$, and the decoder $D_{\omega }:\mathbf{z}\to \overline{\mathrm{X}}$ reconstructs the ST data from $\mathbf{z}$–$\overline{\mathrm{X}}=D_{\omega }\left(\mathbf{z}\right)$ as:

$$\overline{\mathrm{X}}={\mathcal{F}}_{\theta ,\omega }\left(\mathrm{X}\right)=D_{\omega }\left(E_{\theta }\left(\mathrm{X}\right)\right)$$

(6)

The channel anomalies can be transients—short-lived and impacting only a single digi-occupancy map—or persist over time—affecting a sequence of maps. The spatial reconstruction error e to detect a transient anomaly is calculated as:

$$e_i=|x_i-{\overline{x}}_i|$$

(7)

where $x_i\in \mathrm{X}$ and ${\overline{x}}_i\in \overline{\mathrm{X}}$ are the input and reconstructed digi-occupancy of the ith channel. $e_i$ detects a channel abnormality occurrence on isolated maps. We opted for an aggregated error in a time window T using the mean absolute error (MAE) to capture a time-persistent anomaly as:

$$e_{i,MAE}=\frac{1}{T}\sum _{t^{\prime }=t-T}^t{e}_i\left(t^{\prime }\right)$$

(8)

We standardized $e_i$ to regularize the reconstruction accuracy variations among the channels, allowing a single AD decision threshold $\alpha$ to all the channels in the spatial map:

$$s_i=\frac{e_i}{{\sigma }_i}$$

(9)

where ${\sigma }_i$ is the standard deviation of $e_i$, or $e_{i,MAE}$ if the time window is considered, on the training data set. The anomaly flags $a_i$ are generated after applying $\alpha$ to the anomaly scores—$a_i=s_i>\alpha$. $\alpha$ is a tunable constant that controls the detection sensitivity.

4.3. Autoencoder Model Architecture

Convolutional neural networks have achieved state-of-the-art performance in several applications, mainly with image data [18,19,20,21,24]. The shared nature of the kernel filters of CNNs substantially reduces the number of trainable parameters in the model compared to fully connected neural networks. Directly supplying the learned spatial features to temporal neural networks such as RNNs could become inherently challenging due to the considerable computational demand for high-dimensional data. We employed CNNs and GNNs with a pooling mechanism to extract relevant features from high-dimensional spatial data followed by RNNs to capture temporal characteristics of the extracted features (see Figure 7). We integrated the variational layer [61] at the end of the encoder for overfitting regularization by enforcing continuous and normally distributed latent representations [9,31,62,63].

The CNN of the encoder has $L_c$ networks, each containing $Conv3D(·,kernel\_size=[3\times 3\times 3\left]\right)$ for regular spatial learning followed by batch normalization (BN) for network weight regularization and faster convergence, ReLU for nonlinear activation, and MaxPooling3D for spatial dimension reduction. The model can be summarized as:

$$y_t^c,{\psi }_t^c=Pool(ReLU(BN(Conv3D(x_t^l,N_c^l)))){|}_{l=1,\cdots ,L_c}$$

(10)

where $x_t^l$ is the input spatial $\gamma$ map data at time step t, and $N_c^l$ is the feature size of the lth network. $y_t^c$ is the extracted feature set of the CNN at t. $Pool(·)$ denotes $MaxPooling3D(·,stride=[2\times 2\times 2\left]\right)$. ${\psi }_t^c$ holds the pooling spatial location indices of the MaxPooling3D layers to be used later for upsampling in the decoder during map reconstruction. The final extracted feature set ${\mathrm{Y}}_c\in {\mathbb{R}}^{T\times N_c}$ of the CNN is an aggregation of all $y_t^c$ in the time window T, concatenated on the time dimension:

$${\mathrm{Y}}_c=[y_1^c,y_2^c,\dots ,y_T^c]$$

(11)

We used $L_c=4$ to map the input spatial dimension $[64\times 72\times 7]$ into $[4\times 4\times 1]$, which yielded a reduction factor of $2^{L_c}$ and expanded the feature space of the input from $N_f=1$ to $N_c=128$. $N_c^{\prime }:[4\times 4\times 1\times 128]=2048$ spatial features were generated after reshaping.

The GNN of the encoder has $L_g$ networks of a graph convolutional network (GCN) with a ReLU activation, and a final global attention pooling [64]. The networks are summarized as:

$$\begin{array}{c}{y}_t^g=Pool(ReLU(GCN(x_t^l,N_g^l){|}_{l=1,\cdots ,L_g}))\\ {\mathrm{Y}}_g=[y_1^g,y_2^g,\dots ,y_T^g]\end{array}$$

(12)

where the GCN layers have a feature size of $N_g^l$, and $Pool(·)$ signifies the $GlobalAttentionPooling(·)$ at the end of the GNN. GlobalAttentionPooling aggregates the graph node features with an attention mechanism to obtain the final feature set of the GNN ${\mathrm{Y}}_g\in {\mathbb{R}}^{T\times N_g}$. Similar to the CNN, we set $L_g=4$ and $N_g=128$ to generate the ${\mathrm{Y}}_g$.

The encoded ST feature set $\zeta \in {\mathbb{R}}^{1\times N_z}$ is obtained by learning the temporal context on the extracted spatial features $\mathrm{Y}=[{\mathrm{Y}}_c,{\mathrm{Y}}_g]$ with two layers of long short-term memory (LSTM) as:

$$\zeta =LSTM(\mathrm{Y},N_r^l){|}_{l=1,2}$$

(13)

where $N_r^l$ is the feature size of the lth LSTM layer. The last layer ($N_r^2=N_z=32$) generates the low-dimensional latent representation of the encoder. The VAE layer of the encoder generates the normally distributed representation latent features $\mathbf{z}$ as:

$$\begin{array}{c}\mathbf{z}={\mu }_z+{\sigma }_z\odot ϵ\end{array}$$

(14)

where ⊙ signifies an element-wise product with a standard normal distribution sampling $ϵ\sim \mathcal{N}(0,1)$[62]. ${\mu }_z$ and ${\sigma }_z$ of the VAE are implemented with FC layers taking $\zeta$ as input.

The decoder network of the AE is comprised of an RNN and a CNN to reconstruct the target ST data from the latent features. The decoding starts with a temporal feature reconstruction using an LSTM network as:

$$\overline{\zeta }=LSTM(\mathbf{z},N_r^l){|}_{l=1,2}$$

(15)

where $\overline{\zeta }$ is the reconstructed temporal feature set from the latent space $\mathbf{z}$. A spatial reconstruction follows for each time step t through a multilayer deconvolutional neural network. Each network starts with $MaxUnpooling3D(·,stride=[2\times 2\times 2],{\psi }_l^c)$ to upsample the spatial data using localization indices ${\psi }_l^c$ from the lth MaxPooling3D of the encoder followed by a deconvolutional layer ($Deconv3D(·,kernel\_size=[3\times 3\times 3\left]\right)$) [65], a BN, and a ReLU. Eventually, $Deconv3D(·,kernel\_size=[1\times 1\times 1\left]\right)$ is incorporated for the final output stabilization. The decoder network is summarized as:

$$\begin{array}{c}{\overline{x}}_t=ReLU{(BN(Deconv3D(Unpool(\overline{{\zeta }_t},{\psi }_t^c),N_c^l))|}_{l=1,\cdots ,L_c}\\ {\overline{x}}_t=ReLU(Deconv3D({\overline{x}}_t,N_f))\end{array}$$

(16)

where ${\overline{x}}_t$ is the reconstructed spatial data, and $Unpool(·)$ denotes $MaxUnpool3D(·)$. The final reconstructed ST data $\overline{\mathrm{X}}\in {\mathbb{R}}^{T\times N_{i\eta }\times N_{i\varphi }\times N_d\times N_f}$ are obtained as:

$$\overline{\mathrm{X}}=[{\overline{x}}_1,{\overline{x}}_2,...,{\overline{x}}_T]$$

(17)

4.4. Model Training

We trained the AE on healthy digi-occupancy maps of LHC collision runs. The modeling task became a multivariate learning problem since the target data contained readings from multiple calorimeter channels in the spatial digi-occupancy map. An appropriate scaling of the spatial data was thus necessary for effective model training; we further normalized the spatial data per channel into a range of $[0,1]$. We also observed that the $\gamma$ distribution of the channels at the first depth of the spatial map was different from the channels at the higher depths (see Figure 3); a distribution imbalance on target channel data affects model training efficacy when well-known statistical algorithms, e.g., MSE, are employed as loss functions. The MSE loss minimizes the cost of the entire space, and it may converge to a nonoptimal local minimum in the presence of an imbalanced data distribution; this phenomenon is known as the class imbalance challenge in machine learning classification problems. A widely used remedy is to employ a weighting mechanism—assigning weights to the different targets. We applied a weighted MSE loss function to scale the loss from the different distributions, the $depth\in 1$ and $depth\in 2,\cdots ,7$:

$${\mathcal{L}}^{\prime }=\sum _j\frac{{\varsigma }_j}{M_j}\sum _{i\in C_j}{(x_i-{\overline{x}}_i)}^2$$

(18)

where $x_i$ is the $\widehat{\gamma }$ of the ith channel in the jth group set $C_j$, $M_j$ is the number of channels in $C_j$, and ${\varsigma }_j$ is the weight factor of the MSE loss of the jth group. We holistically set ${\varsigma }_1=0.4$ and ${\varsigma }_2=1$ after experimenting with different $\varsigma$ values.

The VAE regularized the training MSE loss using the KL divergence loss $D_{KL}$ to achieve the normally distributed latent space as:

$$\mathcal{L}=\underset{W\in \mathbb{R}}{\mathrm{argmin}}\left\{{\mathcal{L}}^{\prime }-\beta D_{KL}\left[\mathcal{N}({\mu }_z,{\sigma }_z)\parallel \mathcal{N}(0,I)\right]+\rho {\parallel W\parallel }_2^2\right\}$$

(19)

where $\mathcal{N}$ is a normal distribution with zero mean and unit variance, and $\parallel .\parallel$ is the Frobenius norm of the $L_2$ regularization for the trainable model parameters W. $\beta =0.003$ and $\rho ={10}^{-7}$ are tunable regularization hyperparameters. We finally used the Adam optimizer with superconvergence via one-cycle learning rate scheduling [66] for training.

5. Experimental Results and Discussion

AD studies for the DQM inject simulated anomalies into good data to validate the effectiveness of the developed models since a small fraction of the data is affected by real anomalies [48]. Likewise, we trained the GraphSTAD autoencoder model using four GPUs on 10,000 digi-occupancy maps—from LS sequence number $[1,500]$—and evaluated on LSs $[500,1500]$ injected with synthetic anomalies simulating real dead, hot, and degraded calorimeter channels. We employed early stopping using 20% of the training dataset to estimate the validation loss during each training epoch (see Figure 8). The model training achieved good fitting and generalization, as demonstrated by the low loss and closeness between the training and validation losses.

Figure 9 demonstrates the capability of the proposed ST AE in reconstructing normal digi-occupancy maps from a sequence of lumisections. The AE accomplished a promising reconstruction ability on the ST digi-occupancy data. A high reconstruction accuracy on the healthy data is essential to reduce false-positive flags when a semi-supervised AE is employed for AD application. We further discuss the reconstruction error distribution comparison on the healthy and abnormal channels in Section 5.1.2.

We discuss below the performance of our proposed model, comparisons with benchmark models, detection results on real faulty channels, and model complexity cost.

5.1. Anomaly Detection Performance

We created synthetic anomalies to simulate dead, hot, and degraded channels and then injected them into healthy digi-occupancy maps. We subsequently evaluated the ability of the AD to detect the injected anomalies. The anomaly generation algorithm involved three steps: (1) a selection of a random set of LSs $\tau \in [500,1500]$ from the test set, (2) a random selection of spatial locations $\phi$ for each $\tau$, where $\phi \in \{i\eta \times i\varphi \times depth\}$ on the HE axes (see Figure 3), and (3) injection of anomalies. The anomalies were simulated using degrading factor $R_D$ with ${\gamma }_a=R_D{\gamma }_h$, where ${\gamma }_a$ and ${\gamma }_a$ are the healthy and anomaly channel $\gamma$ values: dead ($R_D=0$, and ${\gamma }_a=0$), hot ($R_D>1$, and ${\gamma }_a>>{\gamma }_h$), and degraded ($0<R_D\in <1$, and $0\le {\gamma }_a<{\gamma }_h$). We kept the same $\tau$ and $\phi$ as for the generated anomalies for consistency when evaluating the AD performance of the different anomaly types.

5.1.1. Detection of Dead and Hot Channels

We evaluated the AD accuracy on dead—${\gamma }_a=0,R_D=0$—and hot—${\gamma }_a=R_D{\gamma }_h,R_D=200\%$—channels on the 10,000 maps—5000 maps for each anomaly type.

Table 1. AD on dead and hot channel anomalies on isolated digi-occupancy maps.

Anomaly Type	Captured Anomalies	P	R	F1	FPR
Dead Channel	99%	0.999	0.99	0.995	6.722 × 10−6
	95%	1.000	0.95	0.974	3.102 × 10−6
	90%	1.000	0.90	0.947	2.068 × 10−6
Hot Channel	99%	0.999	0.99	0.994	9.113 × 10−6
	95%	1.000	0.95	0.974	1.939 × 10−6
	90%	1.000	0.90	0.947	1.196 × 10−6

P—precision, R—recall, F1—F1-score, FPR—false positive rate.

and

Table 2. AD on time-persistent dead and hot channel anomalies.

Anomaly Type	Captured Anomalies	P	R	F1	FPR
Dead Channel	99%	0.999	0.99	0.995	7.691 × 10−6
	95%	1.000	0.95	0.974	2.715 × 10−6
	90%	1.000	0.90	0.947	1.616 × 10−6
Hot Channel	99%	0.999	0.99	0.995	5.461 × 10−6
	95%	1.000	0.95	0.974	1.357 × 10−6
	90%	1.000	0.90	0.947	7.756 × 10−7

present the AD performance on transient anomalies—short-lived in isolated maps—and time-persisting anomalies—encroaching consecutive maps in a time window—respectively. Our model achieved a good accuracy with precise localization of the faulty channels—a 0.99 precision when capturing 99% of the faulty channels. Time-persistent anomalies were easier to detect; the FPR generally improved by 13–23% and 28–40% for the dead and hot anomalies, respectively, compared to the short-lived anomalies on isolated LSs. We observed that most false positives (FPs) occurred on channels with a low expected ${\gamma }_h$, where the model achieved a relatively lower reconstruction accuracy. The performance was not entirely unexpected since we trained the AE to minimize a global MSE loss function (19). The reconstruction errors became relatively high for channels with low $\gamma$ ranges that limited the effectiveness in distinguishing the anomalies when capturing 99% of the time-persistent dead channels using (8).

We monitored roughly 31.28 million HE sensor channels, of which 335,000 (1.07%) were simulated abnormal channels, from the 5000 maps on the isolated map evaluation in Table 1. The monitored channels grew to 156 million with 1.68 million (1.07%) anomalies for the evaluation of time-persistent anomalies in Table 2 using five time-window maps resulting in 25,000 maps.

5.1.2. Detection of Degrading Channels

Table 3. AD on time-persistent degraded channels.

Anomaly Type	${\mathit{R}}_{\mathit{D}}$	FPR (90%)	FPR (95%)	FPR (99%)
Degraded Channel	80%	1.636 × 10−3	3.614 × 10−3	2.988 × 10−2
	60%	1.329 × 10−4	3.834 × 10−4	1.550 × 10−3
	40%	8.405 × 10−6	2.764 × 10−5	2.242 × 10−4
	20%	2.263 × 10−6	5.173 × 10−6	2.505 × 10−5
	0%	9.699 × 10−7	1.778 × 10−6	6.142 × 10−6

presents the AD accuracy of time-persistent degraded channels simulated with $R_D=[80\%,60\%,40\%,20\%,0\%]$; $R_D=0\%$ corresponds to a dead channel. We injected the generated channel faults into 1000 maps for each decay factor. We monitored around 156 million channels, of which 1.74 million (1.11%) were abnormal channels, from the total of 25,000 digi-occupancy maps—5000 maps per time window. The AD system demonstrated a promising potential in detecting degraded channel anomalies. The FPR to capture 99% of the anomaly was 2.988%, 0.155%, 0.022%, 0.002%, and 0.001% when channels operated at 80%, 60%, 40%, 20%, and 0% of their expected capacity, respectively.

The relatively lower precision at $R_D=80\%$ indicated that there were still a few anomalies challenging to catch despite the very low FPR considering the accurate classification of numerous true-negative healthy channels (see Figure 10); the channels operating at $R_D=80\%$ were mostly inliers overlapping with the healthy operating ranges, and detecting them was difficult when the expected ${\gamma }_h$ of the channel was low. The significant improvement of the FPR by 88% and 95% when the number of the captured anomalies was reduced to 95% and 90%, respectively, demonstrated a small percentage of the channels caused the performance drop at $R_D=80\%$. Figure 11 illustrates the overlap regions on the distribution of the reconstruction errors of the healthy and faulty channels at the various $R_D$ values.

5.2. Performance Comparison with Benchmark Models

We quantitatively compared alternative benchmark models to validate the capability of GraphSTAD (see Figure 12). The benchmark AE models employed a similar architecture as the GraphSTAD AE but with different layers. The results demonstrated that the integration of the GNN had a significant performance improvement from 1.6 to 3.9 times in the FPR. The temporal models—with RNN—achieved a three- to fivefold boost over the nontemporal spatial AD model when capturing severely degraded channels. The GraphSTAD system had a substantial 25-time amelioration over the nontemporal model for subtle and inlier anomalies, e.g., channels deteriorated by 20% at $R_D=80\%$. Incorporating temporal modeling and a GNN enhanced degrading channel detection performance.

5.3. Detection of Real Anomalies in the HCAL

Our GraphSTAD system caught five real faulty HE channels in collision data RunId = 324841 using the digi-occupancy maps. The faulty channels were located at $[i\eta ,i\varphi ,depth]:[17,71,3]$, $[18,71,3]$, $[18,71,4]$, $[18,71,5]$, and $[28,71,4]$ and impacted 52 consecutive LSs (see Figure 13). Figure 13 and Figure 14 illustrate the detected faults fell into the dead channel category except the last one LS = 57, where the channels operated in a degraded state—the $\gamma$ was lower than expected. Detecting degraded channels is challenging since the $\gamma$ reading is nonextreme as in dead and hot channels, and the $\gamma$ drop overlaps with other false down-spikes (see LS > 57 in Figure 13). The down-spikes in the digi-occupancy for LS > 57 are due to a nonlinearity in the LHC—changes in collision run settings (see Figure 13b); our normalizing regression model successfully handled the fluctuation during prepossessing before causing false-positive alerts (see Figure 13a). Figure 15 and Figure 16 portray the spatial anomaly scores during the death and degraded status of the faulty channels; the high anomaly scores localized at the faulty channels demonstrated the GraphSTAD AD performance at a channel-level granularity. The existing production DQM system of the CMS uses rule-based and statistical methods and has also reported these abnormal channels in a run-level analysis; the results are only available at the end of the run after analyzing all the LSs for the run [46]. Our approach is adaptive to variability in the digi-occupancy maps and provides an anomaly localization that detects faulty channels, including nonextreme degraded channels, per lumisection granularity.

5.4. Cost of Model Complexity

We developed the models with PyTorch and trained them on four GPUs of NVIDIA Tesla V100 SXM3 32GB and an Intel(R) Xeon(R) Platinum 8168 CPU 2.70 GHz. We utilized a time window $T=5$ and batch size $B=8$ for training, and the dimension of a batch was $[B\times T\times N_{i\eta }\times N_{i\varphi }\times N_d\times N_f]$. The training time of the GraphSTAD model was approximately 45 s per epoch. The training iteration epoch 200 achieved good accuracy with a one-cycle learning rate schedule [66]. The nontemporal model—CNN + FC + VAE—was the fastest, and its superiority emanated from its nonrecurrent networks that only analyzed a single map instead of a sequential processing of five maps in a time window. The median inference time of the GraphSTAD system on a single GPU was roughly 0.05 s with a standard deviation of 0.006 s. The integration of the GNN made the inference relatively slower compared to the benchmark models (see Figure 17). The processing cost was within an acceptable range for the CMS production requirement since the input digi-occupancy map was generated at each lumisection with a time interval of 23 s.

6. Conclusions

In this study, we presented a semi-supervised anomaly detection system for the data quality monitoring system of the Hadron Calorimeter using spatio-temporal digi-occupancy maps. We extended the synergy of temporal deep learning developments for the CMS experiment. Our approach addressed modeling challenges, including digi-occupancy map renormalization, learning non-Euclidean spatial behavior, and degrading channel detection. We proposed the GraphSTAD system that combined convolutional, graph, and temporal learning networks to capture spatio-temporal behavior and achieve a robust localization of anomalies at a channel granularity on high-dimensional spatial data. The AD performance evaluation demonstrated the efficacy of the proposed system for channel monitoring. Our proposed AD system will facilitate monitoring and diagnostics of faults in the front-end hardware and software systems of the calorimeter. It will enhance the accuracy and automation of the existing DQM system, providing instant anomaly alerts on a broader range of channel faults in realtime and offline; the improved monitoring of the calorimeter will result in the collection of high-quality physics data. The methods and approaches discussed in this study are domain-agnostic and can be adopted in other spatio-temporal fields, particularly when the data exhibit regular and irregular spatial characteristics.