{'fraction_non_alphanumeric': 0.03556877323420074, 'fraction_numerical': 0.010468401486988847, 'mean_word_length': 4.717735079068186, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On Anomaly Interpretation via Shapley Values

Anomaly localization is an essential problem as anomaly detection is. Because a rigorous localization requires a causal model of a target system, practically we often resort to a relaxed problem of anomaly interpretation, for which we are to obtain meaningful attribution of anomaly scores to input features. In this paper, we investigate the use of the Shapley value for anomaly interpretation. We focus on the semi-supervised anomaly detection and newly propose a characteristic function, on which the Shapley value is computed, specifically for anomaly scores. The idea of the proposed method is approximating the absence of some features by minimizing an anomaly score with regard to them. We examine the performance of the proposed method as well as other general approaches to computing the Shapley value in interpreting anomaly scores. We show the results of experiments on multiple datasets and anomaly detection methods, which indicate the usefulness of the Shapley-based anomaly interpretation toward anomaly localization.


Introduction
Anomaly detection is one of the most prevalent yet difficult problems in data mining and machine learning, and data-driven anomaly detection methodologies have been intensively studied for decades [5,20]. Despite such active researches on anomaly detection, an important subsequent problem of anomaly localization has not necessarily been discussed together. Here, anomaly localization means to identify anomalous features within an anomalous data point, whereas anomaly detection means to find anomalous data points from unlabeled ones. Indeed, the localization problem is as essential as detection because we cannot fully utilize detection results without any clues to understanding why specific instances have been considered to be anomalous.
In principle, a rigorous localization of anomalies can only be performed given some physical or causal models of the target systems, which are not available in many practices. Instead, we often want to obtain partial clues for localizing anomalies by solving somewhat ill-posed problems of anomaly interpretation. A simple yet popular method to interpret anomaly detection is to examine the decomposed or marginal values of anomaly scores used in detection. For example, reconstruction errors, which are often used as an anomaly score, are inherently decomposable into individual features. Also, we may compute the energies of marginal distributions if the model is simple enough. Such straightforward decompositions or marginals are useful in some practices but do not always provide arXiv:2004.04464v1 [cs.
LG] 9 Apr 2020 correct interpretation. For instance, non-anomalous features may also have large reconstruction errors when other features that were correlated in normalcy are anomalous. Moreover, an essential challenge is that we cannot easily compute such straightforward attributions of the anomaly scores based on many modern nonlinear models.
Informally, we can cast the anomaly interpretation as a problem to compute each feature's contribution to the anomalousness. In other words, given some anomaly score values that have been used for anomaly detection, we would like to attribute such values to each input feature for interpreting why a high score has occurred. This is in line with the problem of interpreting machine learning outputs in general, which has been actively studied by many researchers. One of the most famous approaches to this end is the use of the Shapley value, which has been originally discussed in game theory [22] and recently applied to statistical machine learning [16,24,8,18,26,19,2]. Attribution based on the Shapley value is popular mainly because it possesses nice properties such as the insensitivity to irrelevant features. Given the success in general machine learning interpretation, it seems promising to use the Shapley value also for anomaly interpretation, and in fact, a few results have already been reported in this context [3,10,27]. These results are affirmative for using the Shapley value for anomaly interpretation but the evaluations are mostly only qualitative. Moreover, they adopt only general techniques in defining and computing the Shapley value, which are not designed specifically for anomaly scores.
In this paper, we investigate the possibility of using the Shapley value for anomaly interpretation toward anomaly localization. We focus on the semisupervised anomaly detection, which has been actively studied and used in many applications. We newly propose a way to compute characteristic function, upon which the Shapley value is defined, explicitly designed for anomaly scores. This is advantageous to fully explore the usefulness of the Shapley value for anomaly interpretation. The idea of the proposed characteristic function lies in approximating the absence of features by minimizing an anomaly score in the proximity of the data point to be interpreted. We also present practical procedures to compute such a characteristic function.
In addition to the proposed method, we also examine other approaches to defining characteristic function of the Shapley value (and its variants), such as the integrated gradient [26] and kernel SHAP [18]. Moreover, we compare the Shapley-based interpretation to other anomaly interpretation methods [23,31] that can be used in semi-supervised settings. The performance of these methods is examined via experiments on synthetic and real anomalies using several realworld datasets and popular anomaly detection methods.

Background
In this section, we introduce the concept of semi-supervised anomaly detection and related notations. We also review the notion of the Shapley value and the practical challenges in its application.

Semi-Supervised Anomaly Detection
Semi-supervised anomaly detection 3 is a problem to find anomalies given data whose labeled part comprises only normal behaviors of a target system. A typical solution consists of two phases: training and test. First, in the training phase, a model (e.g., density estimation, subspace-based models, clusters, and other generative models) is learned on the data of normal behaviors. Then, in the test phase, a score e(x) : X → R (e.g., reconstruction errors, negative log-likelihoods, and distance from cluster centers, etc.) is computed using the learned model, where x ∈ X is a test data point. We call e(x) anomaly score because a high value of e(x) implies the possibility that x is anomalous. In practice, an alarm is issued if e(x) exceeds some threshold value e th , and possibly x is further inspected by domain experts. Semi-supervised anomaly detection methods are useful in many practices as they do not assume labeled anomalies and generally do not require storing the whole training data in the test phase.
Popular techniques. While we do not go into the details here, let us review some well-known approaches to semi-supervised anomaly detection. For more details, readers can consult excellent surveys such as Chandola et al. (2009) [5] and Pimentel et al. (2014) [20]. One of the most common approaches is the use of probabilistic generative models such as Gaussian mixture models (GMMs). In the training phase we learn a probabilistic density p(x) under some model. In the test phase, given a test data point x, we compute the negative log-likelihood (i.e., energy) as an anomaly score, that is, Another popular choice is the use of autoencoder-like models including denoising autoencoders [29] and variational autoencoders (VAEs) [13] as well as principal component analysis. After learning an encoder f : X → L and a decoder g : L → X , where L is a latent representation space, we can use the reconstruction error of x as an anomaly score, that is, Moreover, these two views are not mutually exclusive. For example, Zong et al. (2018) [32] proposed to use the energy of a GMM learned on the latent representations and the reconstruction errors of an autoencoder for anomaly detection.
Sometimes anomaly scores can be decomposed into each input feature component. For example, the reconstruction error is inherently decomposable as it is computed as the sum of the errors of each feature. As for the energy-based anomaly scores, we may compute the energy of marginal distributions and use them for comparing features. The marginal energy will work well when anomalies occur independently for each feature, but such an assumption may be too restrictive in some practices.

Shapley Value for Interpreting Machine Learning Outputs
The Shapley value [22] is a notion discussed in the study of coalitional games. Now let v : 2 n → R be a set function called characteristic function that represents a game's gain obtained from each coalition, and let N = {1, . . . , n} denote the set of all players. The Shapley value of a game (v, N ) is a way to distribute the total gain v(N ) to each player in accordance with each contribution. The Shapley value of a player i ∈ N , namely φ i , is the weighted average of the marginal contributions, that is, It is preferred as an attribution method because it satisfies several nice properties such as efficiency (the Shapley values sum to the total gain), symmetry (equal contributions result in equal Shapley values), and dummy (Shapley values are zero for non-contributing players). While the Shapley value has been originally studied in the context of the game theory, it has recently been utilized for explaining prediction results of statistical machine learning (e.g., [16,24,8,18,26,19,2,3,10,27]). In such a usage, each input feature corresponds to each player in a game, and the gain of the game means the output prediction by a machine learning method.

Challenges in applications.
Major challenges in utilizing the Shapley value for interpreting machine learning outputs lie in the following two points: (i) How to compute the summation for O(2 n ) terms (i.e., S⊆N {i} in Eq. (1))? (ii) How to define a characteristic function v?
The former is a general difficulty that occurs in computing the Shapley value and not limited to machine learning interpretation. A common solution is via Monte Carlo approximations (see, e.g., [4,18]). In this work, we adopted a Monte Carlo approximation based on the reformulation of Eq. (1) as a weighted least squares problem [6,18]. We defer its concrete procedures to Section 3.3. The latter, how to define v, rises as a challenge specific to interpreting machine learning outputs because v(S) should simulate the absence of the features not in S. This is, in principle, a hard problem that admits no unique solutions because retraining models for exponentially many subsets of N is often computationally prohibitive, and comparing the outputs of retrained models is not straightforward.
In most existing studies on the Shapley value for machine learning outputs, they define a characteristic function v in either or both of the following two approaches. The first approach [26,18,2] is to replace the values of the features not in S by some reference (i.e., baseline) values. Now suppose X = R d , and let h : R d → Z be a machine learning model from the input space to an output space Z. Then the reference-based definition of v is where [x S r S c ] denotes a vector created by replacing x i (i-th element of x) by r i for i ∈ S c ; e.g., Here, S c = N \S is the complement of S, and r ∈ R d is a vector comprising some reference values of the features. A challenge in this reference-based approach is that we need to choose a good reference vector r. General options often employed in practice include zero values (often for images) and statistics of training data (e.g., mean). Another approach [24,8,18] is to marginalize the features not in S. That is, we compute the conditional expectation of h given only where x S denotes the subvector of x corresponding to the indices in S. An obvious challenge in this approach is the computation of the conditional expectation, which is intractable in general. Typically it is approximated via nearest neighbors using training data. Meanwhile, it has been discussed [25] that the Shapley value based on such an expectation-based approach loses some nice properties, including the dummy property, because it depends not only on the target function h but also on the data distribution. Note that, if the independence of features is assumed, this expectation-based approach reduces to (the average of) the reference-based one with multiple reference vectors [18,25].

Shapely Value of Anomaly Scores
We aim to interpret anomaly detection results using their Shapley value. This idea is in line with the existing approaches to interpreting general machine learning outputs via the Shapley value [16,24,8,18,26,19,2], and in fact, some researchers have already reported the use of the Shapley value for anomaly scores interpretation [3,10,27]. However, their definition of the characteristic function is not specialized for anomaly scores, which may limit the performance of the Shapley value as a means of anomaly interpretation.
In the remainder of this section, we first discuss the potential issues in the existing methodologies. Afterward, we present our proposal to define a characteristic function suitable for anomaly scores and show the overall algorithm of the proposed method.

Issues in Existing Approaches
In literature, several researchers discussed the use of Shapley values for interpreting anomaly detection results. Below we briefly review three pieces of related studies [3,10,27].
Two of them are based on the reference-based approach in Eq. (2) to defining the characteristic function, v, for which we need to choose a good reference value r. Since the goodness of a reference depends on target values x and target features S, it should be determined adaptively to both x and S. Antwarg et al. (2019) [3] suggested to interpret the reconstruction of autoencoders using kernel SHAP [18], in which a reference value does not depend on x nor S in principle. Giurgiu and Schumann (2019) [10] proposed to choose a reference value adaptively using the influence weights between a target data point and training data. Such a reference value is chosen adaptively to x, but the same value is used for every S. Also note that their method requires to store enough portion of training data, which may be undesirable in some applications of semisupervised anomaly detection.
In contrast, Takeishi (2019) [27] explicitly adopted the expectation-based approach in Eq. (3). Specifically, they restricted the target model to be the probabilistic principal component analysis and exactly computed the conditional expectation of its reconstruction errors. This approach is free from choosing a good reference value, but the type of applicable anomaly detection models is too restrictive. We can consider an approximation of the conditional expectation via Monte Carlo methods using some approximated probability density, but approximating the conditional distribution under the presence of an anomaly, which is usually out-of-distribution, is not realistic.
Given the above issues in the existing methods, while we adopt the referencebased approach (Eq. (2)) for defining v, we newly propose a way to determine a reference value r adaptively to both x and S. The proposed method can determine r given only x and S, without referring to the training dataset.

Characteristic Function for Anomaly Scores
We propose to define a characteristic function v for anomaly scores as follows. Now recall that we have to design v(S; x) that simulates the absence of the features not in set S. In other words, in our case v(S; x) should represent how anomaly x S is, ignoring the rest, x S c . To this end, we define v(S; x) as the lowest value of the anomaly score e(·) achieved in the neighborhood of x with x S being fixed. Rephrasing this idea more formally, we propose v(S; x) = min where M x ⊂ X is a compact neighborhood of x ∈ X . By defining the characteristic function as in Eq. (4), we can approximately examine how e(x) is high solely due to the features in S. This approach is different from the general methods (e.g., [24,8,26,18,2]) in which the absence of some features is simulated via the replacement with predefined reference values or conditional expectation. Our method basically follows the reference-based approach but is free of the difficulty in preparing good reference values in advance independent of x S . Rather, the proposed method automatically determines the reference values depending on x S via solving the optimization problem in Eq. (4).
Relaxation. Practically, it is unrealistic to determine M x manually for each x because it should depend on both the geometry of X and the property of e(x). Hence, we propose a relaxation of Eq. (4) that admits a less complexity of manual tuning. Instead of taking the minimum in a compact neighborhood as in Eq. (4), we suggest to use a local minimum value of e(·) around x with x S being fixed. That is, we take a local minimum of with regard to y under a constraint y S = x S , where γ ≥ 0 is a hyperparameter. The second term is for ensuring that a local minimum point in the proximity of the original x is chosen. In summary, we suggest a relaxed definition of characteristic function, namelyv, aŝ where the min is considered only in the sense of local minimum. This relaxation is advantageous, especially when e(·) is non-convex, which is often the case with modern anomaly detection methods using nonlinear models. A local minimum point of Eq. (5) is not unique in general, but it would not matter practically; we can approximate to some extent the absence of the features in S c by using an anomaly score smaller than the original e(x) that is obtained changing only x S c .
Further relaxation. To computev in Eq. (6), we need to computex(S) by solving a (local) minimization problem for every S. This is usually prohibitive for a moderate number of features, so we propose a further relaxation. The idea is to use an "ersatz" ofx(S) in Eq. (6). We suggest to use a surrogate reference valuex(S) computed as:x In words, instead of computingx(S) directly for a set S, we use the average ofx of the empty set and the singletons, that is,x(∅) andx({i}) for i ∈ S. Accordingly, we define a further relaxed characteristic function: v(S; x) = e(x(S)).
For computingx(S) (and thusṽ(S; x)), for every S we can reuse the solutions of dim(X ) + 1 minimization problems: where again note that the min is considered in the sense of local minimum. As obtaining a local minimum point close to x requires only moderate efforts, we can quickly computex(S) for moderate-dimensional datasets from dim(X ) = 100 to dim(X ) = 500. Developing more efficient methods for higher-dimensional data remains an open challenge.

Overall Algorithm
The Shapley value-based anomaly interpretation does not premise specific types of anomaly detection methods. Base statistical models can be GMMs, autoencoders, their combinations or ensembles, or anything else. Anomaly scores can be negative log-likelihoods, reconstruction errors, or others. For the method proposed in Section 3.2, the only condition is that ∇ x e(x) can be (approximately) computed with admissible computational costs, which holds in most methods.
We show the overall procedures to compute the Shapley values of an anomaly score in Algorithm 1. The procedures in Lines 4-8 are based on the formulation of Shapley value computation as a weighted least squares problem [6,18]. We recommend to see [1] for more details. For the optimization in Line 2, we used a quasi-Newton method that admitted a fast convergence. The only hyperparameter specific to the proposed method is γ. We empirically found that the performance was not sensitive to γ in a range 0.001 ≤ γ ≤ 0.1 (see Section 5.4 for details). As for the number of samples m, we used the value recommended in the implementation by Lundberg et al. (2017) [18], that is, m = 2 dim(X ) + 2 11 .
Lastly, let us discuss how to use the computed {φ} for anomaly interpretation. Similarly to other general usages, a large value of φ i implies that the i-th feature is contributing to the anomaly score. Here an issue is the sign of φ. We observed that in many cases, the computed φ's were only zero or positive valued. Some cases resulted in negative φ, but we empirically found that regarding only positive-valued φ as anomalous (regardless of absolute values) led to better performance in our experiments. There is no definite strategy on this matter, so how to utilize φ finally should be tried and discussed in accordance with each application.

Related Work
As discussed earlier in Section 3.1, the use of the Shapley value for anomaly interpretation has already been studied [3,10,27], though their methods are not  [14] proposed to extract explanation rules from random forests that separate anomalies and normal training data, for which enough portion of training data has to be stored. Apart from the semi-supervised setting, many researchers have been working on the interpretation of unsupervised anomaly detection (e.g., [15,12,7,30,17]), in which one is given a unlabeled dataset comprising (hopefully) many normal data points and a few anomalous points. These methods are not directly applicable to the semi-supervised detection setting, but it is of great interest to exploit the ideas also in this context.

Experiments
We conducted experiments using several datasets and popular detection methods. The evaluation is both on synthetic anomalies and real anomalies.

General Settings
Datasets and preprocessing. We used public datasets listed in Table 1. The U2R dataset is taken from the NSL-KDD dataset [28], which is a modified version of the KDDCup'99 data. We extracted the U2R attack type part of NSL-KDD and eliminated categorical features and a constant-valued feature, which resulted in the 10-dimensional dataset. The other datasets are from the ODDS repository [21]. As for the Lympho dataset, we converted the categorical features into binary features by one-hot encoding, which resulted in the 59-dimensional dataset.
We split each dataset into a training set D train , a validation set D valid , and a test set pair D test = (D norm test , D anom test ). We used the whole anomaly part of each dataset as D anom test and randomly chose the same size of normal data for D norm test .
Then, D train and D valid were created from the remaining normal data points. Hence, D train , D valid , and D norm test comprise only normal data points, whereas D anom test is made of only anomalous data points. We normalized the real-valued datasets using the mean and the standard average of each training set.
Base anomaly detectors. The anomaly scores e(x) we examined in the experiments are as follows: -The energy (i.e., negative log-likelihood) of GMMs (hereafter GMM ).
-The reconstruction error (i.e., squared distance for the real-valued datasets or cross entropy for the binary-valued) of VAEs (VAE-r). -The negative ELBO of VAEs (VAE-e).
-The energy of deep autoencoding Gaussian mixture models [32] (DAGMM ). The details of the model architectures are presented in Appendix. We selected the models (e.g., the number of mixture components of GMMs) examining the detection performance on D test to pass "best condition" detection results to the interpretation part. Selection using test data is admitted here because our interest is not in improving the detection performance.
Interpretation methods. In addition to the Shapley value, we tried the following methods for interpreting detected anomalies: -Marginal scores for individual features (namely marg). We computed the energy of marginals and the reconstruction error of each feature for GMM and VAE-r, respectively. Note that marg cannot be computed for VAE-e and DAGMM without some approximation or sampling procedures.  [31]. We used the settings suggested in their original paper.
As for the Shapley value, we tried the following variants other than the proposal: -The integrated gradients method (IG) [26], which can be regarded as a realization of the Aumann-Shapley value. We used the mean of training data as a reference value. We refer to the strategy proposed in Section 3 as anomaly score Shapley (ASH). For ASH, we set γ = 0.01 for every experiment unless otherwise stated.

Synthetic Anomaly
Using the normal half of test data, D norm test , we investigated the quantitative performance of the interpretation methods listed above under controlled conditions. That is, we generated synthetic anomalies in D norm test by perturbing some features of normal data and tried to localize them via interpretation results.
Settings. We generated synthetic anomalies as follows. First, we chose a data point from D norm test and randomly selected d anom feature indices. Then, for realvalued datasets, we perturbed those feature values of the chosen data point by adding noise drawn independently and uniformly from [−2, −1] ∪ [1,2]. This magnitude of noise was set because every real-valued dataset was normalized to have the standard deviation being one. For the binary-valued dataset (i.e., Lympho), instead of adding the uniform noise, we flipped (i.e., 0 → 1, 1 → 0) the values of the d anom selected features. Then, we computed attributions by the interpretation methods. We repeated the above procedures for 100 times for each dataset and each interpretation method. Results. Every anomaly detector resulted in a similarly good detection performance, so we can focus only on the interpretation performance. Here we show the results for d anom = 1 and d anom = 3. For d anom = 1, we report the mean reciprocal rank of the ground truth (i.e., the perturbed feature) for the 100 trials for each interpretation method. A large mean reciprocal rank value indicates that the ground truth features tend to be ranked at the top, which implies the interpretation is correct as localization. In Table 2 (in the previous page), the mean reciprocal rank values are listed for all the datasets and the detectors. The Shapley-based methods, such as IG and ASH, perform relatively well, and the proposed approach, ASH, is better than the other Shapley-based ones in most cases. The simple strategy of marg works well in some cases but also fails in some other cases.
For d anom = 3, we report the averages of the area under the receiver operating characteristic curves (AUROC) computed by setting the perturbed features as ground truths. In Table 3, the AUROC values are listed for all the settings. The overall tendency is similar to the d anom = 1 case. While the performance differences are smaller than the previous case, the proposed strategy, ASH, is better than the other methods in most settings.  Discussion. The simple strategy of marg worked well probably because the perturbations were independent of the original data. However, it failed in some cases and can only be computed for limited types of anomaly scores. Within the Shapley-based methods (IG, KSH, wKSH, and ASH), the proposed strategy of ASH resulted in a somewhat stably good performance. Nonetheless, as it sometimes fails while the others are good (e.g., on U2R with GMM in Table 2), using these methods in ensembles may lead to a good performance in practice. As a general observation, DAGMM seems relatively difficult to interpret with more than tens of features, probably because of its non-additive nature. Such "interpretation hardness" will be interesting to elaborate on.

Real Anomaly
Using the anomalous half of test data, D anom test , we examined the performances in more realistic situations. As we usually do not know which features are the ground truths (i.e., anomalous) in real anomalies, we prepared surrogate ground truths using a method [11] like a two-sample test.
Settings. For each dataset, we identified features significantly different between D norm test and D anom test using the method [11], which is like a two-sample test. The features identified as such can work as surrogate ground truths because they are based on the full information of both normal and anomalous data. In contrast, the interpretation methods only refer to a single data point in D anom test in principle. Therefore, if the outputs of the interpretation methods well coincide with such identified features, we can anticipate that the given interpretation is meaningful to some extent in the sense that it agrees with more information-rich results. We computed the agreement of these two for the whole D anom test for each dataset.
Results. In Table 4, we report the mean reciprocal rank of the most different feature (determined as above) of each dataset for each interpretation method. The results are shown only for Thyroid, BreastW, U2R, and Lympho, because referring to a single feature for the higher-dimensional datasets is less meaningful. Compared to the simple methods like marg, the Shapley-based methods, such as IG, KSH, and ASH, agree relatively well with the different features identified upon the full information of the two classes. A specific result is also visually inspected in Figure 1, wherein the distributions of the attributions within D anom test of the U2R dataset are drawn. The features identified as being significantly different in U2R were #9, #7, #10, and #8. In Figure 1, we can see marg, IG, and ASH successfully offer large attributions to features #9,#10, and #8 (while #7 is somewhat missed). However, marg also attributes a little to the feature #1, which may be misleading.
Discussions. The results in Table 4 are less coherent than those of the synthetic anomaly experiments, which is natural as the surrogate ground truths are just surrogate. While detailed evaluation has to be done for each use case, the results imply the usefulness of the Shapley-based anomaly interpretation.

Sensitivity to γ
In Table 5, we report additional results of the synthetic anomaly experiment (in Section 5.2) using the GMM detector with γ being varied from 0.001 to 10. We can observe that the performance in mean reciprocal rank is mostly insensitive to γ in 0.001 ≤ γ ≤ 0.1.

Conclusion
We discussed the use of Shapley value for anomaly interpretation, which will be useful for anomaly localization when no causal model of target systems is provided. We proposed a characteristic function, on which the Shapley value is computed, specifically for anomaly scores. We examined the performance of the proposed method and other approaches using several datasets and anomaly detection methods, and the results imply the usefulness of the Shapley-based methods for anomaly interpretation. Remaining technical challenges include developing more efficient methods to compute the characteristic function.
three non-musk classes are used as normal data, and two musk classes are as anomalous data.
Arrhythmia is originally from the UCI machine learning repository [9], and we used the reformatted dataset provided as a part of the ODDS repository [21]. The reformatted dataset comprises 274 real-valued features. The original problem is a 16-class classification to distinguish between the presence and absence of cardiac arrhythmia. For anomaly detection purposes, the eight smallest classes are used as anomalous data.

A.2 Anomaly detection methods
We present more details on the anomaly detection methods, including the model architectures and the definition of the anomaly scores.
VAEs. As for VAE-r and VAE-e [13], the encoder is a multilayer perceptron with one hidden layer, whereas the decoder is with two hidden layers. Every activation function is the softplus function. We selected best values of the dimensionality of the latent variable, namely dim(z), and the number of the hidden units of the multilayer perceptrons, namely dim(MLP). The candidate of dim(z) was the rounded values of 0.2d, 0.4d, 0.6d, and 0.8d, where d is the dimensionality of each dataset. The candidate of dim(MLP) was the rounded values of 0.5d, d, and 2d. The loss function used for learning was the negative ELBO with the mean squared loss for the real-valued datasets or with the cross-entropy loss for the binary-valued dataset. We used the Adam optimizer with a learning rate 0.001 and stopped the optimization observing the validation set loss. The corresponding anomaly scores are: respectively for VAE-r and VAE-e. Here, θ and ψ denote the sets of parameters of the decoder and the encoder, respectively. Also, g and f are the outputs (corresponding to mean parameters) of the decoder and the encoder, respectively.
DAGMM. As for DAGMM [32], the architectures of the encoder and decoder are the same with the VAEs above. What is specific to DAGMMs is the estimation network, which outputs the estimation of cluster assignment of a GMM learned on the latent representations and the reconstruction error values. In our experiments, the estimation network is a multilayer perceptron with one hidden layer using the softplus function as activation. The candidates of the hyperparameters were the same as the above cases, both for the autoencoder part and the GMM part. We used the Adam optimizer with a learning rate 0.0001 and stopped the optimization observing the validation set loss.

B Additional Results for Synthetic Anomaly Experiment
In Tables 6 and 7, we report the hits@n values (n = 1 and n = 3) for the synthetic anomaly experiment with d anom = 1.

C Additional Results for Real Anomaly Experiment
In Figures 2-5, we show the boxplots of the attributions computed by every interpretation method on D anom test of the U2R dataset.