## 4 Proposed Method
### Model Training and Inference  
Before training the model, we preprocess the raw data by filling missing values with zeros and normalizing data using MinMaxScaler. Then, we transform all observations into window samples, as the model inputs using sliding window. Depending on the monitoring strategies of target systems, MTS data can be either regular or irregular. In the irregular MTS, the data may have various sampling durations or time intervals; thus, it may degrade the performance of anomaly detection methods [35]. However, in the scope of our work, we assume that the data has the same time intervals and lengths among metrics. In practice, the generated data may not strictly follow Gaussian distributions, which leads to the difficulty of optimizing the reconstruction error \(\mathbb{E}_{q_{\theta}(z|x)}[\log p_{\theta}(x|z)]\) in the ELBO. As a result, the model tends to optimize the Kullback-Leibler divergence. That may cause posterior collapse, in which the posterior of the latent variable is equal to its prior. To solve the problem, we adopted \(\beta-\)VAE [36] to put more weight on the first term of the ELBO while optimizing. We set \(\beta=0.2\) as in [37] and observe that this value works well in our experiment. The modified ELBO is defined as:  
\[\mathcal{L}(x,\theta,\phi)=\mathbb{E}_{q_{\theta}(z|x)}[\log p_{\theta}(x|z)] -\beta KL[q_{\phi}(z|x)||p_{\theta}(z)] \tag{5}\]  
If we sample directly from the stochastic variable \(z\sim\mathcal{N}(\mu_{z},\sigma_{z}^{2})\), computing the ELBO when training VAE will be infeasible because the process of sampling is not differentiable. Thus, the stochastic gradient descent cannot be applied to the model. To solve this problem, SGVB proposed the re-parameterization trick, which creates a random variable \(\epsilon\sim\mathcal{N}(0,\mathbf{I})\) for sampling independently with \(z\), such that \(z\) can be rewritten as \(z(\epsilon)=\mu_{z}+\sigma_{z}\epsilon\). As recommended in SGVB, sampling one sample is sufficient for computing the ELBO in the training model; thus, we also set \(L=1\). In the inference phase, we calculate the reconstruction probability as the anomaly score, which can be referred to in [1,10]. To classify the observation at timestep \(t\) into normal or abnormal, the window sample \(w_{t}=\{x_{t-L+1},\ldots,x_{t}\}\) is used to infer from the trained model. Although the model output has all the scores corresponding to each observation in the window sample \(w_{t}\), we only use the score of the last observation \(x_{t}\) for anomaly detection. During detection, a window sample may contain abnormal points, and these points may cause a bias to the embedded representation \(z\), then affect the reconstruction probability. To reduce the effect of anomalies in the input, we adopted MCMC imputation [10,38] to replace anomalous points with normal points. Based on the imputation, a more reasonable embedded representation could be achieved, thus helping the model to detect anomalies more accurately. The procedure for MCMC imputation is described as follows: given an input \(x\); the embedded representation \(z\) is sampled from \(q_{\phi}(z|x)\); then, a reconstruction is obtained from \(pnet\). The last observation in \(x\) is replaced by the last observation in the reconstruction to achieve a new input \(\hat{x}\). We iterate this procedure \(M\) times to reduce the bias and achieve a better representation. Note that the output of current iteration will be the input of the next iteration. After finishing the  
Figure 6: Long-scale module, window size = 30, strides = 2, kernel size = 15 (first layer), no zero padding. MCMC imputation, the reconstruction probability can be calculated using a Monte Carlo integration, which is similar to Equation (4). The anomaly score is:  
\[score_{anomaly}=-\mathbb{E}_{q_{\phi}(z|x)}[\log p_{\theta}(x|z)]=-\frac{1}{L} \sum_{l=1}^{L}[\log p_{\theta}(x|z^{(l)})] \tag{6}\]  
where \(z^{(l)}\) is sampled from \(q_{\phi}(z|\tilde{x})\).