# Black-Box Anomaly Attribution ††thanks: Currently under review. The original
version was received on May 8, 2021.

Tsuyoshi Idé Naoki Abe
IBM Research, Thomas J. Watson Research Center
1101 Kitchawan Rd., Yorktown Heights, NY 10598, USA
<EMAIL_ADDRESS>
###### Abstract

When the prediction of a black-box machine learning model deviates from the
true observation, what can be said about the reason behind that deviation?
This is a fundamental and ubiquitous question that the end user in a business
or industrial AI application often asks. The deviation may be due to a sub-
optimal black-box model, or it may be simply because the sample in question is
an outlier. In either case, one would ideally wish to obtain some form of
attribution score — a value indicative of the extent to which an input
variable is responsible for the anomaly. In the present paper we address this
task of “anomaly attribution,” particularly in the setting in which the model
is black-box and the training data are not available. Specifically, we propose
a novel likelihood-based attribution framework we call the “likelihood
compensation (LC),” in which the responsibility score is equated with the
correction on each input variable needed to attain the highest possible
likelihood. We begin by showing formally why mainstream model-agnostic
explanation methods, such as the local linear surrogate modeling and Shapley
values, are not designed to explain anomalies. In particular, we show that
they are “deviation-agnostic,” namely, that their explanations are blind to
the fact that there is a deviation in the model prediction for the sample of
interest. We do this by positioning these existing methods under the unified
umbrella of a function family we call the “integrated gradient family.” We
validate the effectiveness of the proposed LC approach using publicly
available data sets. We also conduct a case study with a real-world building
energy prediction task and confirm its usefulness in practice based on expert
feedback.

_K_ eywords Explainable AI $\cdot$ Anomaly attribution $\cdot$ Shapley Value
$\cdot$ Integrated gradient $\cdot$ LIME

##  1 Introduction

With the recent advances in machine learning algorithms and their wide-spread
deployment, automated anomaly detection has become a critical component in
many modern industrial systems. In its most ambitious form, it is coupled with
the idea of a “digital twin,” which has recently emerged in the manufacturing
industry [Tao et al., 2018, Fuller et al., 2020]. The “digital twin” captures
the behavior of an entire production system with a machine learning model. As
the model is a replica of the system under the normal operating conditions,
any significant deviation from its prediction implies the presence of some
anomalies. The expectation then is that critical incidents can be prevented by
leveraging the digital twin model in a cycle of prediction, evaluation, and
risk mitigation.

Despite the initial optimism, however, broad deployment of digital twins is
still far from reality. One of the biggest concerns, as perceived by the end-
user, is the lack of explainability, and hence actionability, of a typical
digital twin model. Consider, for instance, the use-case of building energy
management (see Sec. 6 for more details). The energy consumption of a
building, $y$, is predicted with a regression function $y=f(\bm{x})$, where
$\bm{x}$ is a vector of measurements such as the outside temperature and
humidity. The monitoring system typically consists of a few sub-components
including those for HVAC (heating, ventilating, and air conditioning), sensing
and data collection, and data management. Since they are often developed by
different vendors using proprietary technologies, the end-user may not have
first-hand knowledge on the prediction model and the training data used to
train it, even when they have the ownership of the overall system. In such a
scenario, the prediction system acts as black-box, often leaving the end-user
in the dark in cases of anomalies, i.e. prediction deviations.

In this paper, we address the task of anomaly attribution in the black-box
regression setting. We assume the model is “doubly black-box,” meaning that we
neither have access to the parametric form of the regression function
$f(\bm{x})$ nor its training data (see Table 1 for more details). In this
especially challenging setting, we strive to provide a reliable “attribution
score” for each of the input variables with regard to the anomaly in question.
See Fig. 1 for illustration. We call this task “anomaly attribution.” A large
deviation from the observation may be due to sub-optimal model training, or
simply because the observed sample is an outlier. In either case, the
attribution score indicative of the extent to which an input variable is
responsible for the anomalous output will be useful in aiding the end-user
make better decisions on what action to take, e.g. using their knowledge on
the system.

Figure 1: Problem setting and motivation. (a) Given a black-box regression
model and anomalous sample(s), our goal is to quantify input variables’
responsibility without using training data. (b) Existing attribution methods
attempt to explain either the local gradient or the increment from a reference
point $\bm{x}^{0}$, rather than the deviation.

Anomaly attribution has previously been studied as a sub-task of anomaly
detection. For instance, in subspace-based anomaly detection, computing each
variable’s responsibility has been part of the standard procedure in the
white-box setting [Chandola et al., 2009, Jiang et al., 2011, Dang et al.,
2013b, Dang et al., 2013a, Micenková et al., 2013]. On the other hand, in the
field of explainable artificial intelligence (XAI), growing attention has been
paid on “post-hoc” explanations of black-box prediction models. Among numerous
XAI techniques, as conveniently summarized in recent review articles [Burkart
and Huber, 2021, Molnar, 2022, Samek et al., 2019, Arrieta et al., 2020,
Speith, 2022], there are a few methods that address model-agnostic feature
attribution in the regression setting: 1) Local linear surrogate modeling,
which is best known under the name LIME (Local Interpretable Model-agnostic
Explanations) [Ribeiro et al., 2016], 2) Shapley value (SV), which was first
introduced to the machine learning community by [Štrumbelj and Kononenko,
2010], and 3) integrated gradient (IG) [Sundararajan et al., 2017]. We
collectively call these approaches the function-based approach as they work
with the prediction function $f(\bm{x})$ directly for attribution. One obvious
limitation of these methods is that they are designed to explain what
$f(\bm{x})$ looks like at a test sample $\bm{x}^{t}$, rather than to explain
the deviation $f(\bm{x}^{t})-y^{t}$.

In this paper, we propose a novel likelihood-based framework to anomaly
attribution in the doubly black-box regression setting. Likelihood is a
legitimate starting point for explaining deviation as it is a canonical metric
of non-anomalousness. We begin by pointing out that existing attribution
methods, LIME, SV, and IG, are inherently deviation-agnostic, and, in fact,
are not appropriate for anomaly explanation. Interestingly, one can show that
they are derived from a unified framework we call the integrated gradient
family, and the deviation-agnostic property is a general characteristic of the
IG family. Based on the solid understanding of the function-based approach, we
propose the new notion of likelihood compensation (LC), which seeks a local
perturbation that achieves the highest possible likelihood. As explained later
(see Fig. 2 for a preview), LC as “deviation measured horizontally” is a
useful attribution score as it is interpreted as an action that might be taken
to bring back the outlying sample to normalcy. To the best of our knowledge,
this is the first principled framework for model-agnostic anomaly attribution
in the regression setting.

The notion of LC was first introduced in our conference paper [Idé et al.,
2021]. This paper expands it by re-positioning the entire approach in a
unified framework. In doing so, we formally show, for the first time, a close
relationship among the function-based attribution algorithms, including the
equivalence between SV and EIG. We also conduct a systematic empirical study
using new datasets on the deviation agnostic property, score variability
issues, and the consistency among different attribution approaches. In the
last subsection, we present a real-world use-case, where the theory of LC made
a substantial difference.

##  2 Related Work

This section summarizes prior works in the context of anomaly attribution.

###  2.1 General background: Anomaly attribution in doubly black-box setting

As mentioned earlier, anomaly attribution has been studied as a sub-task of
anomaly detection in the machine learning community, typically in the white-
box unsupervised setting. In the modern deployment of AI systems, however, a
black-box situation often arises [Li et al., 2022]. Table 1 compares white-
and black-box settings, although the definition of a black-box model varies in
the literature. In the context of XAI, the white- and grey-box cases are
typically associated with deep neural networks (DNNs). Since the number of
network parameters is extremely large, the model is often viewed as black-box
even with full access to the internal network parameters. Saliency maps
[Simonyan et al., 2014, Selvaraju et al., 2017] and layer-wise relevance
propagation [Montavon et al., 2019] are well-known DNN-specific attribution
methods that fall into these categories.

The doubly black-box scenario occurs when the end-user has access only to the
model’s API (application programming interface) or does not have full
understanding of the algorithm used. The latter can occur even when the source
code is available. Since the nature of dependency on the internal parameters
is unknown, for any XAI approach to be applicable to this setting it must be
model-agnostic, making most of the DNN-specific methods inapplicable.

Table 1: Comparison of different XAI settings. Our focus is the doubly black-box case. | model API access | model internal access | training data access
---|---|---|---
white-box | yes | yes | yes
grey-box | yes | yes | no
doubly black-box | yes | no | no

###  2.2 Local linear modeling, Shapley value (SV), and integrated gradient
(IG)

As discussed before, we focus on model-agnostic post-hoc anomaly attribution
in the doubly black-box regression setting. Here, 1) local linear surrogate
modeling, 2) Shapley value (SV), and 3) integrated gradient (IG) are the main
existing approaches that are potentially applicable to our task. Let us
quickly review recent works of these approaches.

Local linear modeling has been extensively used for attribution for decades,
often under the name of sensitivity analysis [Abhishek and Kamath, 2022].
Recent applications to anomaly explanation include [Giurgiu and Schumann,
2019] and [Zhang et al., 2019]. The former used LIME for attribution inspired
by the kernel SHAP approach [Lundberg and Lee, 2017]. The latter also used
LIME for attribution with a new regularization approach. Many gradient-based
attribution methods can be viewed as local linear modeling, although many of
them are DNN-specific, assuming full access to the internal parameters.

SV is one of the most popular attribution methods in the AI community, and the
last few years have seen many attempts to apply SV to various industrial
domains. For instance, [Hwang and Lee, 2021] proposed to use SV for sensor
fault diagnosis, [Mariadass et al., 2022] used SV to explain unexpected
observations in crop yield analysis, and [Antwarg et al., 2021] used SV to
explain unusual warranty claims.

IG [Sundararajan et al., 2017] is another generic input attribution approach
potentially applicable to the black-box setting. The application of IG to
anomaly explanation can be found in, e.g., [Sipple, 2020, Sipple and Youssef,
2022].

Table 2 compares the key properties of LC, the proposed algorithm, with LIME,
SV, and IG along with two additional methods: The expected integrated gradient
(EIG), which generalizes IG by taking the expectation with respect to the
baseline input (see Sec. 4.1 for the detail) and the $Z$-score, which
quantifies the deviation of each input variable from its expected value,
independently of $y$. Since the true data distribution is unknown in general,
the expectation has to be computed as the empirical approximation on the
training data. The same applies to SV, which make them inapplicable to the
doubly black-box setting, as shown in the ‘training-data-free’ column. IG does
not need the training data, but it does need an extra piece of information of
the baseline input (the ‘baseline-free’ column). Except for those domains in
which a universally accepted data pre-processing method has been established,
it is generally hard to choose a valid baseline input. The dependence on the
baseline input is considered a major factor that limits practical utility of
IG in anomaly attribution.

One fundamental issue with the existing approaches is that they are deviation-
agnostic (the ‘$y$-sensitive’ column), which will be mathematically shown in
Sec. 4. This means that LIME, SV, and (E)IG, in fact, do not explain the
reasons for the sample of interest to be anomalous: As illustrated in Fig. 1
(b), they explain either the local gradient or the increment of $f(\bm{x})$,
rather than what may have caused the deviation $f(\bm{x})-y$ at a test point
$(\bm{x},y)=(\bm{x}^{t},y^{t})$. This situation remains unchanged even if we
apply these methods to the modified function $F(\bm{x},y)\triangleq
f(\bm{x})-y$, as discussed later.

Table 2: Comparison of different anomaly attribution methods in the regression setting.  | training-data-free | baseline-free |  $y$-sensitive | reference point
---|---|---|---|---
LIME | yes | yes | no | infinitesimal vicinity
SV | no | yes | no | globally distributional
IG | yes | no | no | arbitrary
EIG | no | yes | no | globally distributional
Z-score | no | yes | no | global mean of predictors
LC | yes | yes | yes | maximum likelihood point

###  2.3 Unified attribution framework

One of our contributions is establishing a unified framework for many of the
popular “function-based” attribution methods. Prior work along this line
includes [Deng et al., 2021], which attempts to characterize IG using Taylor
expansion and proposes to use the expectation to neutralize the need for a
specific baseline input. [Sundararajan and Najmi, 2020] is another important
work aiming at building a unified attribution framework. The authors pointed
out that there can be a few different definitions for SV, such as the baseline
SV and expected SV, and discussed the relationship with IG in a qualitative
manner. [Lundberg and Lee, 2017] reintroduced the SV-based attribution method
originally proposed by [Štrumbelj and Kononenko, 2010] and proposed a hybrid
method that lies between SV and LIME. Also, [Kumar et al., 2020] analyzed SV’s
risk of producing misleading attribution because of the gap between
conditional and marginal distributions, while [Zhou et al., 2022] conducted a
systematic empirical study to compare different attribution methods including
LIME and SV.

Inspired by these works, we go one step further in this paper: Using power
expansion, we mathematically show certain equivalence properties of the
function-based attribution methods, including the equivalence between SV and
EIG (Theorem 4), which naturally lead to the notion of the IG family. We then
point out their two fundamental limitations in anomaly attribution that have
been hitherto unnoticed. One is the deviation-agnostic property, and the other
is the explicit or implicit dependency on arbitrary baseline points (as
summarized in the ‘reference point’ column in Table 2). An in-depth analysis
of the latter brings us to the new notion of likelihood-based attribution
proposed in this paper, as discussed in Sec. 5.1.

##  3 Problem Setting

As mentioned earlier, we focus on the task of anomaly attribution in the
regression setting rather than classification or unsupervised settings. Figure
1 (a) summarizes the overall problem setting. Suppose we have a
(deterministic) regression model $y=f(\bm{x})$ in the doubly black-box setting
(see Table 1). Neither the training data set $\mathcal{D}_{\mathrm{train}}$
nor the (true) distribution of $\bm{x}$ is available. Throughout the paper,
the input variable $\bm{x}\in\mathbb{R}^{M}$ and the output variable
$y\in\mathbb{R}$ are assumed to be noisy real-valued, where $M$ is the
dimensionality of the input vector. We also assume that queries to get the
response $f(\bm{x})$ can be performed cheaply at any $\bm{x}$.

In practice, anomaly attribution is typically coupled with anomaly detection:
When we observe a test sample $(\bm{x},y)=(\bm{x}^{t},y^{t})$, we first
compute an anomaly score $a^{t}=a(\bm{x}^{t},y^{t})$ to quantify how anomalous
it is. Then, if $a^{t}\in\mathbb{R}$ is high enough, we go to the next step of
anomaly attribution. In this scenario, the task of anomaly attribution is
defined as follows.

######  Definition 1 (anomaly attribution).

Given a black-box regression model $y=f(\bm{x})$, compute the score for each
input variable indicative of the extent to which an input variable is
responsible for the sample being anomalous.

We can readily generalize the problem to that of collective anomaly detection
and attribution. Specifically, given a test data set
$\mathcal{D}_{\mathrm{test}}=\\{(\bm{x}^{t},y^{t})\mid
t=1,\ldots,N_{\mathrm{test}}\\}$, where $t$ is the index for the $t$-th test
sample and $N_{\mathrm{test}}$ is the number of test samples, we can consider
an anomaly score as well as attribution score for the whole test set
$\mathcal{D}_{\mathrm{test}}$. We will see later an example where a daily
attribution score is computed from 24 hourly observed measurements.

The standard approach to computing the anomaly score is to use the negative
log-likelihood of the test sample(s) (See, e.g., [Lee and Xiang, 2000,
Yamanishi et al., 2000, Staniford et al., 2002, Noto et al., 2010]). Assume
that, from the deterministic regression model, we can somehow obtain
$p(y|\bm{x})$, a probability density over $y$ given the input signal $\bm{x}$.
Under the i.i.d. assumption, the anomaly score can be written as

$\displaystyle a(\bm{x}^{t},y^{t})=-\ln p(y^{t}\mid\bm{x}^{t}),\quad\mbox{or,
}\quad
a(\mathcal{D}_{\mathrm{test}})=-\frac{1}{N_{\mathrm{test}}}\sum_{t\in\mathcal{D}_{\mathrm{test}}}\ln
p(y^{t}\mid\bm{x}^{t}),$ (1)

corresponding to the single sample case and collective case, respectively.
Obviously, one challenge here is how to estimate $p(y\mid\bm{x})$ from the
deterministic regression function. We provide one such approach in Sec. 5.2.

Given an anomalous sample $(\bm{x}^{t},y^{t})$ and the distribution
$p(y\mid\bm{x})$, computing the anomaly score is straightforward. However,
computing anomaly attribution score is more challenging. This is in some sense
an inverse problem: The function $f(\bm{x})$ readily gives an estimate of $y$
from $\bm{x}$, but, in general, there is no obvious way to do the reverse in
the multivariate case. When an estimate $f(\bm{x}^{t})$ looks ‘bad’ in light
of an observed $y^{t}$, what can we say about the contribution, or
responsibility, of the respective input variables? Section 5 provides our
proposed answer to this question.

###  3.1 Notation

We use boldface to denote vectors. The $i$-th dimension of a vector
$\bm{\delta}$ is denoted as $\delta_{i}$. The $\ell_{1}$ and $\ell_{2}$ norms
of a vector are denoted by $\|\cdot\|_{1}$ and $\|\cdot\|_{2}$, respectively,
and are defined as $\|\bm{\delta}\|_{1}\triangleq\sum_{i}|\delta_{i}|$ and
$\|\bm{\delta}\|_{2}\triangleq\sqrt{\sum_{i}\delta_{i}^{2}}$. The sign
function $\mathrm{sign}(\delta_{i})$ is defined as being $1$ for
$\delta_{i}>0$, and $-1$ for $\delta_{i}<0$. For $\delta_{i}=0$, the function
takes an indeterminate value in $[-1,1]$. For a vector input, the definition
applies element-wise, yielding a vector of the same size as the input vector.

We distinguish between a random variable and its realizations via the absence
or presence of a superscript. For notational simplicity, we use $p(\cdot)$ as
a proxy to represent different probability distributions, whenever there is no
confusion. For instance, $p(\bm{x})$ is used to represent the probability
density of a random variable $\bm{x}$ while $p(y|\bm{x})$ is a different
distribution of another random variable $y$ conditioned on $\bm{x}$. The
Gaussian distribution of a scalar variable $y$ is defined as

$\displaystyle\mathcal{N}(y\mid
m,\sigma^{2})\triangleq\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left\\{-\frac{(y-m)^{2}}{2\sigma^{2}}\right\\}$
(2)

where $m$ is the mean and $\sigma^{2}$ is the variance. The multivariate
Gaussian distribution is defined in a similar way.

##  4 Limitations of Function-Based Anomaly Attribution

To motivate the likelihood-based attribution approach presented in the next
section, this section shows fundamental limitations of the existing function-
based attribution approaches: IG, SV, and LIME are inherently deviation-
agnostic and are not appropriate for anomaly attribution. For simplicity, we
assume for now that the derivative of the black-box regression function $f$ is
computable somehow to an arbitrary order. We discuss numerical gradient
estimation approaches in Sec. 5.4.

###  4.1 Deviation-agnostic property of integrated gradient

#### Definitions

The notion of integrated gradient (IG) was first introduced to the AI
community by Sundararajan et al. [Sundararajan et al., 2017] as a method for
axiomatic derivation of an input attribution method. For a test sample at
$\bm{x}^{t}$, IG of the black-box regression function $f$ for the $i$-th
variable is defined by

$\displaystyle\mathrm{IG}_{i}(\bm{x}^{t}\mid\bm{x}^{0})$
$\displaystyle\triangleq(x_{i}^{t}-x_{i}^{0})\int_{0}^{1}\mathrm{d}\alpha\
\left.\frac{\partial f}{\partial
x_{i}}\right|_{\bm{x}^{0}+(\bm{x}^{t}-\bm{x}^{0})\alpha},$ (3)

where $\bm{x}^{0}$ is called the baseline input, a parameter representing a
“default” value of the input. Despite the intimidating look, the integration
term simply computes an averaged gradient w.r.t. $x_{i}$ on the line from
$\bm{x}^{0}$ to $\bm{x}^{t}$. Hence,
$\mathrm{IG}_{i}(\bm{x}^{t}\mid\bm{x}^{0})$ is $x_{i}$’s contribution to the
increment of $f$ when moving from $\bm{x}^{0}$ to $\bm{x}^{t}$. Despite the
term “gradient,” IG is not a gradient but represents the increment (See Fig. 1
(b) and Table 2). This fact becomes clearer if we expand $\partial f/\partial
x_{i}$ into the Taylor series w.r.t. $\alpha$ and perform integration:

$\displaystyle\mathrm{IG}_{i}(\bm{x}^{t}\mid\bm{x}^{0})$
$\displaystyle=\frac{\partial f(\bm{x}^{0})}{\partial
x_{i}}\Delta_{i}+\frac{1}{2!}\sum_{j=1}^{M}\frac{\partial^{2}f(\bm{x}^{0})}{\partial
x_{i}\partial
x_{j}}\Delta_{i}\Delta_{j}+\frac{1}{3!}\sum_{j,k=1}^{M}\frac{\partial^{3}f(\bm{x}^{0})}{\partial
x_{i}\partial x_{j}\partial x_{k}}\Delta_{i}\Delta_{j}\Delta_{k}+\ldots,$ (4)

where we have defined $\Delta_{i}\triangleq x^{t}-x^{0}_{i}$, etc. The right
hand side is simply the collection of differential increments over different
orders in the power expansion. Intuitively, the $i$-th attribution score gets
large if the $i$-th gradient is large and if the test sample is far from the
baseline point along the $i$-th axis.

As pointed out by [Sipple, 2020], one of the major issues of IG is the need
for the baseline input. Since real-world data may often follow a multi-peaked
distribution (see Fig. 14 for an actual example), there may not exist a
clearly defined default value. One natural approach for addressing this issue
is to integrate out $\bm{x}^{0}$ using a distribution $P(\bm{x})$:

$\displaystyle\mathrm{EIG}_{i}(\bm{x}^{t})$
$\displaystyle\triangleq\\!\\!\int\\!\mathrm{d}\bm{x}^{0}P(\bm{x}^{0})\mathrm{IG}_{i}(\bm{x}^{t}\mid\bm{x}^{0})=\\!\\!\int\\!\mathrm{d}\bm{x}^{0}P(\bm{x}^{0})(x_{i}^{t}-x_{i}^{0})\int_{0}^{1}\\!\\!\mathrm{d}\alpha\left.\frac{\partial
f}{\partial x_{i}}\right|_{\bm{x}^{0}+(\bm{x}^{t}-\bm{x}^{0})\alpha}\\!\\!,$
(5)

which we call the expected integrated gradient (EIG). Obviously, EIG is
reduced to IG as a special case when $P(\bm{x})$ is chosen to be Dirac’s delta
function. When computing EIG, $P(\bm{x})$ should ideally be the true
distribution or its empirical approximation using the training dataset.
Unfortunately, neither of them is available in our setting (see Table 2).

#### IG and EIG for deviation

In the context of anomaly attribution, we are interested in explaining the
deviation $f(\bm{x})-y$ rather than $f(\bm{x})$ itself. Let us define a new
function $F(\bm{x},y)\triangleq f(\bm{x})-y$ and consider EIG for this
function. As an input attribution method, the deviation version of EIG,
denoted as a two-place function $\mathrm{EIG}_{i}(\bm{x}^{t},y^{t})$, is
defined by

$\displaystyle\mathrm{EIG}_{i}(\bm{x}^{t},y^{t})$
$\displaystyle\triangleq\int\mathrm{d}y^{0}\int\mathrm{d}\bm{x}^{0}\
P(\bm{x}^{0},y^{0})\mathrm{IG}_{i}(\bm{x}^{t},y^{t}\mid\bm{x}^{0},y^{0})$ (6)
$\displaystyle\mathrm{IG}_{i}(\bm{x}^{t},y^{t}\mid\bm{x}^{0},y^{0})$
$\displaystyle\triangleq(x_{i}^{t}-x_{i}^{0})\int_{0}^{1}\mathrm{d}\alpha\left.\frac{\partial
F}{\partial x_{i}}\right|_{\bm{x}^{0}+(\bm{x}^{t}-\bm{x}^{0})\alpha,\
y^{0}+(y^{t}-y^{0})\alpha}$ (7)

for $i=1,\ldots,M$, where $P(\bm{x},y)$ is the joint distribution between
$\bm{x}$ and $y$. The following property holds:

######  Theorem 1.

IG and EIG are deviation-agnostic.

###### Proof.

Since $\frac{\partial F}{\partial x_{i}}=\frac{\partial f}{\partial x_{i}}$,
the statement about IG obviously holds. For EIG, the integration w.r.t.
$y^{0}$ produces $\int\mathrm{d}y^{0}P(\bm{x}^{0},y^{0})=P(\bm{x}^{0})$,
yielding $\mathrm{EIG}_{i}(\bm{x}^{t},y^{t})=\mathrm{EIG}_{i}(\bm{x}^{t})$. ∎

This means that IG and EIG explain what the regression surface looks like at
$\bm{x}^{t}$ regardless of the nature of the deviation. Hence, they may not be
the best approach if our interest is in explaining the deviation.

#### Lower-order approximations

The power expansion approach introduced in Eq. (4) can be used also for EIG.
In this case, however, the expansion should be around $\bm{x}^{t}$:

$\displaystyle\left.\frac{\partial f}{\partial
x_{i}}\right|_{\bm{x}^{0}+\alpha\bm{\Delta}}\\!\\!\\!\\!\\!\\!\\!\\!=\frac{\partial
f(\bm{x}^{t})}{\partial
x_{i}}+(\alpha-1)\sum_{j=1}^{M}\frac{\Delta_{j}}{1!}\frac{\partial^{2}f(\bm{x}^{t})}{\partial
x_{i}\partial
x_{j}}+(\alpha-1)^{2}\sum_{j,k=1}^{M}\frac{\Delta_{j}\Delta_{k}}{2!}\frac{\partial^{3}f(\bm{x}^{t})}{\partial
x_{i}\partial x_{j}\partial x_{k}}+\ldots,$ (8)

where $\bm{\Delta}\triangleq\bm{x}^{t}-\bm{x}^{0}$ and $\Delta_{i}\triangleq
x^{t}-x^{0}_{i}$. We have used the fact
$(\bm{x}^{0}+\alpha\bm{\Delta})-\bm{x}^{t}=(\alpha-1)\bm{\Delta}$. This
expansion allows performing the integration w.r.t. $\alpha$ analytically:

$\displaystyle\mathrm{EIG}_{i}(\bm{x}^{t})=\langle\Delta_{i}\rangle\frac{\partial
f(\bm{x}^{t})}{\partial
x_{i}}-\sum_{j=1}^{M}\frac{\langle\Delta_{i}\Delta_{j}\rangle}{2!}\frac{\partial^{2}f(\bm{x}^{t})}{\partial
x_{i}\partial
x_{j}}+\sum_{j,k=1}^{M}\frac{\langle\Delta_{i}\Delta_{j}\Delta_{k}\rangle}{3!}\frac{\partial^{3}f(\bm{x}^{t})}{\partial
x_{i}\partial x_{j}\partial x_{k}}-\ldots,$ (9)

where $\langle\cdots\rangle\triangleq\int\mathrm{d}\bm{x}\ \cdots P(\bm{x})$.
The first and the second terms provide the first- and second-order
approximations of EIG, respectively.

#### Sum rules

Equations (4) and (9) readily lead to the following important property:

$\displaystyle\sum_{i=1}^{M}\mathrm{IG}_{i}(\bm{x}^{t}\mid\bm{x}^{0})=f(\bm{x}^{t})-f(\bm{x}^{0}),\quad\quad\sum_{i=1}^{M}\mathrm{EIG}_{i}(\bm{x}^{t})=f(\bm{x}^{t})-\left\langle
f\right\rangle,$ (10)

###### Proof.

Equation (4) is the same as the power expansion of $f(\bm{x})$ around
$\bm{x}^{0}$ evaluated at $\bm{x}=\bm{x}^{t}$ except for $f(\bm{x}^{0})$, the
first term of the Taylor series. Hence, the first equation holds. The same
argument applies to Eq. (9) to prove the second equation. ∎

These sum rules allow the interpretation that $\mathrm{IG}_{i}$ and
$\mathrm{EIG}_{i}$ are the share of the $i$-th variable in the total change
$f(\bm{x}^{t})-f(\bm{x}^{0})$ and $f(\bm{x}^{t})-\langle f\rangle$,
respectively. In EIG, this sum rule implies that EIG is to contrast the local
output $f(\bm{x}^{t})$ at the test point $\bm{x}^{t}$ with the global mean. If
a certain local distribution at $\bm{x}=\bm{x}^{t}$ is used for $P(\bm{x})$,
we have $\langle f\rangle\approx f(\bm{x}^{t})$, resulting in a meaningless
attribution score. In SV, the corresponding property (16) is called the
efficiency [Roth, 1988].

###  4.2 Deviation-agnostic property of Shapley value

#### Definition

The Shapley value (SV) is one of the most popular attribution metrics in the
AI community. There are a few different versions in SV in the literature,
depending on how the absence of variables is defined. Here we adopt the
definition of the conditional expectation SV [Sundararajan and Najmi, 2020]:

$\displaystyle\mathrm{SV}_{i}(\bm{x}^{t})$
$\displaystyle=\frac{1}{M}\sum_{k=0}^{M-1}\binom{M-1}{k}^{-1}\sum_{\mathcal{S}_{i}:|\mathcal{S}_{i}|=k}\left[\langle
f\mid x_{i}^{t},\bm{x}_{\mathcal{S}_{i}}^{t}\rangle-\langle
f\mid\bm{x}_{\mathcal{S}_{i}}^{t}\rangle\right].$ (11)

Due to the combinatorial nature, this definition appears rather complicated.
Here $\mathcal{S}_{i}$ denotes any subset of the variable indices
$i\in\\{1,\ldots,M\\}$ that does not include $i$ and $|\mathcal{S}_{i}|$ is
its size. The second summation runs over all possible choices of
$\mathcal{S}_{i}$ under the constraint $|\mathcal{S}_{i}|=k$ from the first
summation. We also define the complement $\bar{\mathcal{S}}_{i}$, which is the
subset of $\\{1,\ldots,M\\}$ excluding $i$ and $\mathcal{S}_{i}$. For example,
if $M=12,i=3$ and $\mathcal{S}_{i}=\\{1,2\\}$, the complement
$\bar{\mathcal{S}}_{i}$ will be $\\{4,5,\ldots,12\\}$. Corresponding to this
division, we rearrange the $M$ variables as
$\bm{x}=(x_{i},\bm{x}_{\mathcal{S}_{i}},\bm{x}_{\bar{\mathcal{S}}_{i}})$.

In Eq. (11), the prefactor $\frac{1}{M}$ is there to average over the possible
choices of $|\mathcal{S}_{i}|$. Similarly, the binomial coefficient
$\binom{M-1}{|\mathcal{S}_{i}|}^{-1}$ is to average over all the choices of
$\mathcal{S}_{i}$, which is given by the number of combinations of choosing
$|\mathcal{S}_{i}|$ variables from the $M-1$ variables excluding $i$. This
means that SV is essentially the average of $\left[\langle f\mid
x_{i}^{t},\bm{x}_{\mathcal{S}_{i}}^{t}\rangle-\langle
f\mid\bm{x}_{\mathcal{S}_{i}}^{t}\rangle\right]$, where

$\displaystyle\langle f\mid x_{i}^{t},\bm{x}_{\mathcal{S}_{i}}^{t}\rangle$
$\displaystyle\triangleq\int\\!\\!\mathrm{d}\bm{x}\;P(\bm{x})f(x_{i}=x_{i}^{t},\bm{x}_{\mathcal{S}_{i}}=\bm{x}^{t}_{\mathcal{S}_{i}},\bm{x}_{\bar{\mathcal{S}}_{i}}),$
(12) $\displaystyle\langle f\mid\bm{x}^{t}_{\mathcal{S}_{i}}\rangle$
$\displaystyle\triangleq\int\\!\\!\mathrm{d}\bm{x}\;P(\bm{x})f(x_{i},\bm{x}_{\mathcal{S}_{i}}=\bm{x}_{\mathcal{S}_{i}}^{t},\bm{x}_{\bar{\mathcal{S}}_{i}}).$
(13)

Here $P(\bm{x})$ is the true distribution of $\bm{x}$, which is not available
in our setting (see Table 2). In Eq. (12), the integration is reduced to the
expectation over the marginal distribution of
$\bm{x}_{\bar{\mathcal{S}}_{i}}$. Note that these quantities have to capture
some of the global properties of the data generating mechanism. If $P(\bm{x})$
were a localized distribution at $\bm{x}^{t}$, the difference would simply be
zero.

#### SV for deviation

Similarly to $\mathrm{EIG}_{i}(\bm{x}^{t},y^{t})$ in Eq. (6), we define
$\mathrm{SV}_{i}(\bm{x}^{t},y^{t})$ for the function
$F(\bm{x},y)=f(\bm{x})-y$. As an input attribution method, we need to replace
Eqs. (12) and (13) with

$\displaystyle\langle F\mid
x_{j}^{t},\bm{x}_{\mathcal{S}_{j}}^{t},y^{t}\rangle$
$\displaystyle\triangleq\int\mathrm{d}y\int\mathrm{d}\bm{x}\;P(\bm{x},y)F(x_{i}=x_{i}^{t},\bm{x}_{\mathcal{S}_{i}}=\bm{x}^{t}_{\mathcal{S}_{i}},\bm{x}_{\bar{\mathcal{S}}_{i}},y=y^{t}),$
(14) $\displaystyle\langle F\mid\bm{x}^{t}_{\mathcal{S}_{j}},y^{t}\rangle$
$\displaystyle\triangleq\int\mathrm{d}y\int\mathrm{d}\bm{x}\;P(\bm{x},y)F(x_{i},\bm{x}_{\mathcal{S}_{i}}=\bm{x}^{t}_{\mathcal{S}_{i}},\bm{x}_{\bar{\mathcal{S}}_{i}},y=y^{t}),$
(15)

respectively, to get $\mathrm{SV}_{i}(\bm{x}^{t},y^{t})$. Again, the following
property holds:

######  Theorem 2.

SV is deviation-agnostic.

###### Proof.

Since $F$ is linear in $y$, we can easily see that $\langle F\mid
x_{j}^{t},\bm{x}_{\mathcal{S}_{j}}^{t},y^{t}\rangle=\langle f\mid
x_{j}^{t},\bm{x}_{\mathcal{S}_{j}}^{t}\rangle-y^{t}$ and $\langle
F\mid\bm{x}^{t}_{\mathcal{S}_{j}},y^{t}\rangle=\langle
f\mid\bm{x}^{t}_{\mathcal{S}_{j}}\rangle-y^{t}$ hold, which implies
$\mathrm{SV}_{i}(\bm{x}^{t},y^{t})=\mathrm{SV}_{i}(\bm{x}^{t})$. ∎

#### Sum rule

Finally, SV meets the condition called the efficiency [Roth, 1988]:

$\displaystyle\sum_{i=1}^{M}\mathrm{SV}_{i}(\bm{x}^{t})=f(\bm{x}^{t})-\langle
f\rangle.$ (16)

###### Proof.

See Appendix A. ∎

The efficiency condition is exactly the same as EIG’s sum rule in Eq. (10).
Similarly to the case of EIG, the efficiency condition implies that SV is
essentially the share of the differential increment between the local value
$f(\bm{x}^{t})$ and the global mean. Hence, $P(\bm{x})$ cannot be a local
approximation as the one used in LIME.

###  4.3 Deviation-agnostic property of LIME

In general, the local linear surrogate modeling approach fits a linear
regression model locally to explain a black-box function in the vicinity of a
given test sample $(\bm{x}^{t},y^{t})$. Algorithm 1 summarizes the local
anomaly attribution procedure based on this approach to explain the deviation
$f(\bm{x})-y$.

Algorithm 1 Local linear surrogate modeling for anomaly attribution

1:Regression model $f(\bm{x})$, test point $(\bm{x}^{t},y^{t})$,
regularization parameter $\nu$.

2:Randomly populate $N_{s}$ points
$\\{\bm{x}^{t[1]},\ldots,\bm{x}^{t[N_{s}]}\\}$ in the vicinity of $\bm{x}^{t}$
($N_{s}\sim 1000)$.

3:Compute the deviation $z^{t[n]}\triangleq f(\bm{x}^{t[n]})-y^{t}$ for all
$n$.

4:Fit a linear model $z=\beta_{0}+{\bm{\beta}}^{\top}\bm{x}$ using $\nu$ on
$\\{(\bm{x}^{t[n]},z^{t[n]})\mid n=1,\ldots,N_{s}\\}$.

5:return ${\bm{\beta}}$, which is the local attribution score at
$(\bm{x}^{t},y^{t})$.

In LIME, an $\ell_{1}$-regularized model is used to get a sparse and thus
easy-to-interpret score. Rather surprisingly, despite the modification to fit
$f(\bm{x})-y$ rather than $f(\bm{x})$, the following property holds:

######  Theorem 3.

LIME is deviation-agnostic.

###### Proof.

With $\nu$ being the $\ell_{1}$ regularization strength, the lasso loss
function for LIME is written as

$\displaystyle\Psi(\bm{\beta},\beta_{0})$
$\displaystyle=\frac{1}{N_{s}}\sum_{n=1}^{N_{s}}(z^{t[n]}-\beta_{0}-\bm{\beta}^{\top}\bm{x}^{t[n]})^{2}+\nu\|\bm{\beta}\|_{1},$
$\displaystyle=\frac{1}{N_{s}}\sum_{n=1}^{N_{s}}(f(\bm{x}^{t[n]})-(y^{t}+\beta_{0})-\bm{\beta}^{\top}\bm{x}^{t[n]})^{2}+\nu\|\bm{\beta}\|_{1},$

which is equivalent to the lasso objective for LIME with the intercept
$y^{t}+\beta_{0}$. Since the lasso objective is convex, the solution
$\bm{\beta}$ is unique. With an adjusted intercept, the attribution score
$\bm{\beta}$ remains the same. ∎

In the local linear surrogate modeling approach, the final attribution score
can vary depending on the nature of the regularization term. For theoretical
analysis below, we use a generic algorithm by setting $\nu\to 0_{+}$ in
Algorithm 1, and call the resulting attribution score $\mathrm{LIME}^{0}_{i}$
for $i=1,\ldots,M$. As is well-known, $\mathrm{LIME}^{0}_{i}$ is a local
estimator of $\partial f/\partial x_{i}$ at $\bm{x}=\bm{x}^{t}$.

###  4.4 Unifying LIME and SV into IG

We showed that (E)IG, SV, and LIME all share the same deviation-agnostic
property. We also showed that EIG and SV satisfy the same sum rule. These
findings suggest that they may share a common mathematical structure. This is
indeed the case, as discussed below.

#### Relationship between SV and EIG

The combinatorial definition SV is a major obstacle in getting deeper insights
into what it really represents. With that in mind, we look at the definition
(11) from a somewhat different angle. We again use the Taylor expansion around
$\bm{x}^{t}$ for $f$ in the integrand of Eqs (12) and (13), which leads to

$\displaystyle\langle f\mid
x_{i}^{t},\bm{x}_{\mathcal{S}_{i}}^{t}\rangle-\langle
f\mid\bm{x}_{\mathcal{S}_{i}}^{t}\rangle=\langle\Delta_{i}\rangle\frac{\partial
f(\bm{x}^{t})}{\partial
x_{i}}-\frac{1}{2}\langle\Delta_{i}^{2}\rangle\frac{\partial^{2}f(\bm{x}^{t})}{\partial
x_{i}^{2}}-\sum_{k\in\bar{\mathcal{S}}_{i}}\langle\Delta_{i}\Delta_{k}\rangle\frac{\partial^{2}f(\bm{x}^{t})}{\partial
x_{i}\partial x_{k}}-\ldots.$ (17)

The first and second terms on the r.h.s. do not depend on the choice of
$\mathcal{S}_{i}$, given $i$. In the third term, a $k\in\\{1,\ldots,M\\}$
($k\neq i$) will not be included in $\bar{\mathcal{S}}_{i}$ if it is chosen in
${\mathcal{S}}_{i}$. Thus, for a given value of $|\mathcal{S}_{i}|$, the total
number of appearances of the $k$ in $\sum_{\mathcal{S}_{i}}$ is
$\binom{M-2}{|\mathcal{S}_{i}|}$ because it is the same as the number of
combinations of choosing $|\mathcal{S}_{i}|$ variables out of the $M-2$
variables excluding the $k$ in addition to the $i$. Using the following
identity on the quotient of the binomial coefficients (see, e.g., Chap.4 of
[Gross, 2016])

$\displaystyle\sum_{a=0}^{M-2}\binom{M-1}{a}^{-1}\binom{M-2}{a}=\frac{M}{2},$
(18)

we have the second-order approximation of SV as

$\displaystyle\mathrm{SV}_{i}(\bm{x}^{t})$
$\displaystyle\approx\langle\Delta_{i}\rangle\frac{\partial
f(\bm{x}^{t})}{\partial
x_{i}}-\frac{1}{2}\sum_{k=1}^{M}\langle\Delta_{i}\Delta_{k}\rangle\frac{\partial^{2}f(\bm{x}^{t})}{\partial
x_{i}\partial x_{k}},$ (19)

which is exactly the same as the first two terms in the EIG expansion in Eq.
(9). We have just completed the proof of the following theorem, which reveals
what is behind the seemingly complicated combinatorial definition of SV in Eq.
(11).

######  Theorem 4 (Equivalence of SV to EIG).

The Shapley value is equivalent to the expected integrated gradient up to the
second order.

To the best of our knowledge, this is the first result directly establishing
the correspondence between SV and IG. In Sec. 6, we empirically show that
indeed SV and EIG systematically give similar attribution scores.

#### Relationship between LIME and EIG

LIME as a local linear surrogate modeling approach differs from EIG and SV in
two regards. First, LIME does not need the true distribution $P(\bm{x})$.
Instead, it uses a local distribution to populate local samples. Second, LIME
is defined as the gradient, not a differential increment. These observations
lead us to an interesting question: Is the derivative of EIG in the local
limit the same as the LIME attribution score?

To answer this question affirmatively, consider a local distribution around
$\bm{x}^{t}$ in the following form:

$\displaystyle
P_{\eta}(\bm{x}^{0}\mid\bm{x}^{t})=\mathcal{N}(\bm{x}^{0}\mid\bm{x}^{t},\eta\mathsf{I}_{M})\quad\mbox{with}\quad\eta\to
0,$ (20)

where $\mathsf{I}_{M}$ is the $M$-dimensional identity matrix. With this
distribution, we have

$\displaystyle\langle
x^{0}_{i}-x_{i}^{t}\rangle=0,\quad\langle(x^{0}_{i}-x_{i}^{t})(x^{0}_{k}-x_{k}^{t})\rangle=\delta_{i,k}\eta,$

where $\delta_{i,k}$ is Kronecker’s delta function that takes 1 only when
$i=k$ and 0 otherwise. Notice that the second order term is proportional to
$\eta$ and vanishes as $\eta\to 0$. The other higher-order terms that appear
in the power expansion are either zero or vanish as $\eta\to 0$. An immediate
consequence from the expression (8) is

$\displaystyle\lim_{\eta\to 0}\mathrm{EIG}_{i}(\bm{x}^{t})=0,$ (21)

which confirms the previous discussion for Eq. (10). On the other hand, the
derivative of EIG becomes

$\displaystyle\frac{\partial\mathrm{EIG}_{i}(\bm{x}^{t})}{\partial x_{i}}$
$\displaystyle=\int\mathrm{d}\bm{x}^{0}\
P_{\eta}(\bm{x}^{0}\mid\bm{x}^{t})\left[\frac{\partial f}{\partial
x_{i}}+(x^{t}_{i}-x^{0}_{j})\frac{\partial^{2}f}{\partial
x_{i}^{2}}\right]_{\bm{x}^{0}+(\bm{x}^{t}-\bm{x}^{0})\alpha}\to\frac{\partial
f(\bm{x}^{t})}{\partial x_{i}}$ (22)

as $\eta\to 0$. As this limit is equivalent to $\bm{x}^{0}\to\bm{x}^{t}$ in IG
in Eq. (3), the l.h.s. approaches the derivative of IG. Since the local linear
surrogate modeling estimates local gradient, we have just proved the following
theorem:

######  Theorem 5 (LIME and IG).

The derivative of IG and EIG is equivalent to LIME:

$\displaystyle\mathrm{LIME}^{0}_{i}(\bm{x}^{t})=\lim_{\eta\to
0}\frac{\partial\mathrm{EIG}_{i}(\bm{x}^{t})}{\partial
x_{i}}=\lim_{\bm{x}^{0}\to\bm{x}^{t}}\frac{\partial\mathrm{IG}_{i}(\bm{x}^{t}\mid\bm{x}^{0})}{\partial
x_{i}},$ (23)

where $P(\bm{x}^{0})=\mathcal{N}(\bm{x}\mid\bm{x}^{t},\eta\mathsf{I}_{M})$ is
used in the definition of EIG in Eq. (5).

Since EIG, SV, and LIME can be derived from or associated with IG as shown
above, it is legitimate to say that they are in the integrated gradient
family.

###  4.5 Summary of limitations

We have shown that IG, EIG, SV, and LIME are deviation-agnostic in Theorems 1,
2, and 3. Since our goal is to provide an actionable explanation on a detected
anomaly, this can be a serious issue. We have also shown that SV and LIME are
derived from or associated with IG in Theorems 4 and 5. As suggested by IG’s
Taylor series representation in Eq. (4), the attribution score of the IG
family is essentially the differential increment of $f(\bm{x})$ when going
from the baseline input $\bm{x}^{0}$ to the test point $\bm{x}^{t}$ (see Fig.
1 (b)). As already discussed, the biggest issue here is the arbitrariness of
the baseline input. EIG and SV attempt to neutralize it by integrating out
$\bm{x}^{0}$ in exchange for the demanding requirement on the availability of
the global distribution $P(\bm{x})$.

These two issues – the deviation-agnostic property and the explicit or
implicit dependency on the arbitrary baseline input – are inherent to the IG
family. A key idea in the proposed framework is to leverage the point that
gives the highest possible likelihood in the vicinity of a test sample as the
reference point for attribution (c.f. Table 2). In the next section, we will
show that this idea indeed succeeds in eliminating these issues.

##  5 Likelihood Compensation

This section presents a novel likelihood-based framework for anomaly
attribution.

###  5.1 Seeking reference point through likelihood

#### Definition of LC

In a typical anomaly detection scenario, samples in the training dataset are
assumed to have been collected under normal conditions, and hence, the learned
function $y=f(\bm{x})$ represents normalcy as well. As discussed in Sec. 3,
the canonical measure of anomalousness is negative log likelihood $-\ln
p(y\mid\bm{x})$. A low likelihood value signifies anomaly, and vice versa.
From a geometric perspective, on the other hand, being an anomaly implies
deviating from a certain normal value. We are interested in integrating these
two perspectives.

Suppose we just observed a test sample $(\bm{x}^{t},y^{t})$ being anomalous
because of a low likelihood value. Given the regression function
$y=f(\bm{x})$, there are two possible geometric interpretations on the
anomalousness (see Fig. 2). One is to start with the input
$\bm{x}=\bm{x}^{t}$, and observe the deviation $f(\bm{x}^{t})-y^{t}$. In some
sense, $(\bm{x},y)=(\bm{x}^{t},f(\bm{x}^{t}))$ is a reference point against
which the observed sample $(\bm{x}^{t},y^{t})$ is judged. The other is to
start with the output $y=y^{t}$, and move horizontally, looking for a
perturbation $\bm{\delta}$ such that $\bm{x}=\bm{x}^{t}+\bm{\delta}$ gives the
maximum possible fit to the normal model. In this case, the reference point is
$(\bm{x}^{t}+\bm{\delta},y^{t})$ and $\bm{\delta}$ is a “horizontal
deviation.” Since $\bm{\delta}$ is supposed to be zero if the sample is
perfectly normal, each component $\delta_{1},\ldots,\delta_{M}$ can be viewed
as a value indicative of the responsibility of each input variable.

Based on the intuition above, we propose a new likelihood-based attribution
scoring framework with the following optimization problem:

$\displaystyle\bm{\delta}^{*}=\arg\max_{\bm{\delta}}\left\\{\ln
p(y^{t}\mid\bm{x}^{t}+\bm{\delta})\right\\}\quad\mbox{subject to
}\quad\bm{x}^{t}+\bm{\delta}\in\mathrm{vic}(\bm{x}^{t}),$ (24)

where $\mathrm{vic}(\bm{x}^{t})$ reads “in the vicinity of $\bm{x}^{t}$.” We
call the $\bm{\delta}^{*}$ the likelihood compensation (LC), as it compensates
for the loss in likelihood incurred by the anomalous prediction. Notice that
$\bm{\delta}^{*}$ is defined through $p(y\mid\bm{x})$, and hence, the
randomness of $y$ is automatically taken into account. In other words, $y^{t}$
does not have to be absolutely correct.

In the collective anomaly detection/attribution case corresponding to
$a(\mathcal{D}_{\mathrm{test}})$ in Eq. (1), the LC score is defined as

$\displaystyle\bm{\delta}^{*}=\arg\max_{\bm{\delta}}\left\\{\frac{1}{N_{\mathrm{test}}}\sum_{t=1}^{N_{\mathrm{test}}}\ln
p(y^{t}\mid\bm{x}^{t}+\bm{\delta})\right\\}\quad\mbox{subject to
}\quad\bm{x}^{t}+\bm{\delta}\in\mathrm{vic}(\bm{x}^{t}),$ (25)

which obviously includes Eq. (24) as a special case. In the collective case,
the resulting attribution is an averaged explanation for
$\mathcal{D}_{\mathrm{test}}$. For example, it explains what was wrong with
yesterday overall, rather than explaining about a specific moment in the day.
Section 6.6 provides such a real-world scenario.

The LC score has unique features compared to that of the IG family. First, it
is deviation-sensitive. This is obvious because the log-likelihood itself is
the measure of anomalousness, and the deviation should be reflected in the
anomaly score. Second, it is more principled as a solution to anomaly
attribution because both anomaly detection and attribution are formalized on a
common ground of likelihood. Third, it provides richer information on the
input variables. $\bm{\delta}$ is a correction to the input to get the
reference point having the highest possible likelihood in the vicinity of
$\bm{x}^{t}$, admitting an interpretation like “the sample would have been
normal if the input value had been $\bm{x}^{t}+\bm{\delta}$.” This is one way
of analyzing anomalies with a counterfactual hypothesis, a unique property
lacking in the IG family.

#### Optimization problem for LC

To compute the LC score, we need to specify the functional form of
$p(y\mid\bm{x})$. We employ a Gaussian-based observation model

$\displaystyle p(y^{t}\mid\bm{x}^{t}+\bm{\delta})=\mathcal{N}(y^{t}\mid
f(\bm{x}^{t}+\bm{\delta}),\sigma^{2}(\bm{x}^{t})).$ (26)

We discuss how to estimate the variance $\sigma^{2}(\bm{x}^{t})$ in Sec. 5.2.

In addition, we need to define the vicinity. The vicinity constraint can be
incorporated as regularization on $\bm{\delta}$. This can be problem-specific.
For example, if there is an infeasible region in the domain of $\bm{x}$, the
regularization should penalize such a choice of $\bm{\delta}$ that makes
$\bm{x}^{t}+\bm{\delta}$ infeasible. If there are no known constraints in the
domain of $\bm{x}$, it should be designed to properly address a well-known
issue of $\ell_{1}$-regularization in the original LIME: In the presence of
multiple correlated explanatory variables, lasso tends to pick one at random
[Roy et al., 2017], which can be problematic in attribution. Here, we propose
to use the elastic net regularization [Hastie et al., 2009]. Now the
optimization problem on a perturbation $\bm{\delta}$ is written as:

$\displaystyle\bm{\delta}^{*}=\arg\min_{\bm{\delta}}\left\\{\frac{1}{N_{\mathrm{test}}}\sum_{t=1}^{N_{\mathrm{test}}}\frac{\left[y^{t}-f(\bm{x}^{t}+\bm{\delta})\right]^{2}}{2\sigma^{2}(\bm{x}^{t})}+\frac{1}{2}\lambda\|\bm{\delta}\|_{2}^{2}+\nu\|\bm{\delta}\|_{1}\right\\}.$
(27)

This is the main problem formulation studied in this paper. Note that this
includes $N_{\mathrm{test}}=1$ as a special case.

#### Is Gaussian general enough?

Here we briefly discuss how the Gaussian-based formulation (26) does not
result in much loss of generality. Clearly, the distribution of
$f(\bm{x}^{t})$ is not always Gaussian in general. Notice, however, Eq. (26)
says that the deviation or the error $f(\bm{x}^{t})-y^{t}$ should follow
Gaussian. This is exactly the same situation when Carl Friedrich Gauss
invented Gaussian-based fitting [Brereton, 2014]: Planetary motions do not
follow Gaussian but the error does. See Fig. 14 for a real example in our
context.

#### What if $y^{t}$ is incorrect?

We have provided an intuition of LC as the deviation measured horizontally.
This may lead to a question of what if $y^{t}$ is incorrect. As commented
below Eq. (24), our framework views $y$ as a random variable and it does not
have to be absolutely correct. In fact, the optimization problem (27) shows
that the resulting attribution score $\bm{\delta}^{*}$ does depend on the
variance of $y$, given $\bm{x}$. Of course, it is possible that the error in
$y^{t}$ happens to go far beyond the reasonable range assumed in
$p(y\mid\bm{x})$. In that case, however, the attribution problem itself would
be ill-posed to any attribution methods.

#### Relationship with adversarial training

The optimization problem of LC (24) can be rewritten as

$\displaystyle\min_{\bm{\delta}:\
\bm{x}^{t}+\bm{\delta}\in\mathrm{vic}(\bm{x}^{t})}\left\langle\mbox{Loss}(\bm{x}^{t}+\bm{\delta}\mid
y^{t},\bm{\theta})\right\rangle,$ (28)

where $\mathrm{Loss}$ denotes the loss function, which is the negative log
likelihood in our case, and $\bm{\theta}$ is the model parameters that are
actually not accessible in our doubly black-box setting. Also,
$\langle\cdot\rangle$ denotes empirical average over
$\\{(\bm{x}^{t},y^{t})\\}$. This form is reminiscent of the min-max problem in
adversarial training [Madry et al., 2018, Qin et al., 2019]:

$\displaystyle\min_{\bm{\theta}}\max_{\bm{\delta}:\
\bm{x}^{t}+\bm{\delta}\in\mathrm{vic}(\bm{x}^{t})}\left\langle\mbox{Loss}(\bm{x}^{t}+\bm{\delta}\mid
y^{t},\bm{\theta})\right\rangle,$ (29)

where $y^{t}$ is typically the class label of the $t$-th training sample,
unlike ours. The similarity is obvious, but they are working towards the
opposite directions. LC’s starting point is that the sample is anomalous, and
$\bm{\delta}$ is to bring it back to a normal point. In contrast, adversarial
training assumes the samples are normal, and $\bm{\delta}$ is to make the
normal sample as adversarial as possible by changing the output significantly.
Also a requirement for an example to be adversarial is that the change should
be imperceptible, which is not important in our case.

Figure 2: Illustration of likelihood compensation (LC). (a) For a given test
sample $(y^{t},\bm{x}^{t})$, LC seeks a perturbation $\bm{\delta}$ that
achieves the best possible fit with the black-box regression model
$f(\bm{x})$. (b) The iterative updates Eqs. (37)-(38) converge when the
deviation or the (smoothed) gradient vanishes. See Sec. 5.3 for more details.

###  5.2 Deriving probabilistic prediction model

So far we have assumed the predictive distribution $p(y\mid\bm{x})$ is given
as Eq. (26), i.e.,

$\displaystyle p(y\mid\bm{x})=\mathcal{N}(y\mid
f(\bm{x}),\sigma^{2}(\bm{x})).$ (30)

In this Gaussian observation model, the only parameter to be estimated is
$\sigma^{2}(\bm{x})$. If there are too few test samples, we have no choice but
to set $\sigma^{2}(\bm{x})$ to a constant using prior knowledge. Otherwise, we
can obtain an estimate of $\sigma^{2}(\bm{x})$ using a subset of
$\mathcal{D}_{\mathrm{test}}$ in a cross-validation (CV)-like fashion as
follows.

Let $\mathcal{D}^{t}_{\mathrm{ho}}=\\{(\bm{x}^{(n)},y^{(n)})\mid
n=1,\ldots,N_{\mathrm{ho}}\\}\subset\mathcal{D}_{\mathrm{test}}$ be a held-out
(‘ho’) data set that does not include a given test sample
$(\bm{x}^{t},y^{t})$. Here $N_{\mathrm{ho}}$ is the number of samples in it.
For the observation model Eq. (26) and the test sample $\bm{x}^{t}$, we
consider a locally weighted version of maximum likelihood:

$\displaystyle\max_{\sigma^{2}}\sum_{n=1}^{N_{\mathrm{ho}}}w_{n}(\bm{x}^{t})\ln
p(y^{(n)}\mid\bm{x}^{(n)})=\max_{\sigma^{2}}\sum_{n=1}^{N_{\mathrm{ho}}}w_{n}(\bm{x}^{t})\left\\{\ln\frac{1}{\sqrt{2\pi\sigma^{2}}}-\frac{(y^{(n)}-f(\bm{x}^{(n)}))^{2}}{2\sigma^{2}}\right\\},$
(31)

where $w_{n}(\bm{x}^{t})$ is the similarity between $\bm{x}^{t}$ and
$\bm{x}^{(n)}$. A reasonable choice of this is

$\displaystyle
w_{n}(\bm{x}^{t})=w_{0}+\exp\left\\{-\frac{1}{2\eta_{0}^{2}}\|\bm{x}^{(n)}-\bm{x}^{t}\|^{2}\right\\},$
(32)

where $w_{0}$ and $\eta^{2}$ are constants. The maximizer of Eq. (31) can be
easily found by taking the derivative w.r.t. $\sigma^{-2}$. The solution is
given by

$\displaystyle\sigma^{2}(\bm{x}^{t})=\sum_{n=1}^{N_{\mathrm{ho}}}\frac{w_{n}(\bm{x}^{t})}{\sum_{m}w_{m}(\bm{x}^{t})}\left[y^{(n)}-f(\bm{x}^{(n)})\right]^{2}.$
(33)

This has to be computed for each $\bm{x}^{t}\in\mathcal{D}_{\mathrm{test}}$.
When LC scores are compared over different $t$’s, too much variability in
$\sigma^{2}(\bm{x}^{t})$ tends to obfuscate meaningful signals. For
standardized data, $\eta_{0}=1$ and $w_{0}\gtrsim 5$ would be a reasonable
choice.

###  5.3 Deriving updating equation

Although seemingly simple, solving the optimization problem (27) is generally
challenging. Due to the black-box nature of $f$, we do not have access to the
parametric form of $f$, let alone the gradient. In addition, as is the case in
deep neural networks, $f$ can be non-smooth (see the red curves in Fig. 2
(b)), which makes numerical estimation of the gradient tricky.

To derive an optimization algorithm, we first note that there are two origins
of non-smoothness in the objective function in (27). One is inherent to $f$
while the other is due to the added $\ell_{1}$ penalty. To separate them, let
us denote the objective function in Eq. (27) as
$J(\bm{\delta})+\nu\|\bm{\delta}\|_{1}$, where

$\displaystyle
J(\bm{\delta})\triangleq\frac{1}{N_{\mathrm{test}}}\sum_{t=1}^{N_{\mathrm{test}}}\frac{\left[y^{t}-f(\bm{x}^{t}+\bm{\delta})\right]^{2}}{2\sigma^{2}_{t}}+\frac{1}{2}\lambda\|\bm{\delta}\|_{2}^{2}.$
(34)

Since we are interested only in a local solution in the vicinity of
$\bm{\delta}=\bm{0}$, it is natural to adopt an iterative update algorithm
starting from $\bm{\delta}\approx\bm{0}$. Suppose that we have an estimate
$\bm{\delta}=\bm{\delta}^{\mathrm{old}}$ that we wish to update. If we have a
reasonable approximation of the gradient in its vicinity, denoted by
$\llangle\nabla J(\bm{\delta}^{\mathrm{old}})\rrangle$, the next estimate can
be found by

$\displaystyle\bm{\delta}^{\mathrm{new}}=\arg\min_{\bm{\delta}}$
$\displaystyle\left\\{J(\bm{\delta}^{\mathrm{old}})+(\bm{\delta}-\bm{\delta}^{\mathrm{old}})^{\top}\llangle\nabla
J(\bm{\delta}^{\mathrm{old}})\rrangle+\frac{1}{2\kappa}\|\bm{\delta}-\bm{\delta}^{\mathrm{old}}\|_{2}^{2}+\nu\|\bm{\delta}\|_{1}\right\\}$
(35)

in the spirit of the proximal gradient [Parikh et al., 2014], where $\kappa$
is a hyperparameter representing the learning rate. Notice that the first
three terms in the curly bracket correspond to a second-order approximation of
$J(\bm{\delta})$ in the vicinity of $\bm{\delta}^{\mathrm{old}}$. We find the
best estimate under this approximation.

Fortunately, the r.h.s. has an analytic solution. By differentiating the
objective in Eq. (35) w.r.t. $\bm{\delta}$ and equating the result to
$\bm{0}$, we have the condition of optimality as

$\displaystyle\bm{\delta}=\bm{\phi}-\kappa\nu\
\mathrm{sign}(\bm{\delta}),\quad\quad\mbox{where}\quad\bm{\phi}\triangleq\bm{\delta}^{\mathrm{old}}-\kappa\llangle\nabla
J(\bm{\delta}^{\mathrm{old}})\rrangle.$ (36)

For the $i$-th dimension, if $\phi_{i}>\kappa\nu$ holds, we have
${\phi}_{i}\pm\kappa>0$ and thus $\phi_{i}-\kappa\nu\
\mathrm{sign}(\delta_{i})$ must be positive. By setting
$\mathrm{sign}(\delta_{i})=1$, we conclude $\delta_{i}=\phi_{i}-\kappa\nu$ in
this case. Similar arguments easily verify $\delta_{i}=\phi_{i}+\kappa\nu$ for
$\phi_{i}<-\kappa\nu$. An interesting situation arises when
$|\phi_{i}|\leq\kappa\nu$. Remember that the sign function takes an
indeterminate value within $[-1,1]$ at zero. If $\delta_{i}>0$ is assumed,
$\mathrm{sign}(\delta_{i})=+1$ and the r.h.s. of Eq. (36) must be
$\phi_{i}-\kappa\nu$, which is negative and contradicts the positivity
assumption. Thus, the only possible choice is $\delta_{i}=0$. To summarize,
the solution of Eq. (35) is given by

$\displaystyle\delta_{i}=\begin{cases}\phi_{i}-\kappa\nu,&\phi_{i}>\kappa\nu\\\
0,&|\phi_{i}|\leq\kappa\nu\\\
\phi_{i}+\kappa\nu,&\phi_{i}<-\kappa\nu\end{cases}.$ (37)

Performing differentiation, we see that $\bm{\phi}$ is given by

$\displaystyle\bm{\phi}$
$\displaystyle=(1-\kappa\lambda)\bm{\delta}^{\mathrm{old}}+\kappa\frac{1}{N_{\mathrm{test}}}\sum_{t=1}^{N_{\mathrm{test}}}\left\\{\frac{y^{t}-f(\bm{x}^{t}+\bm{\delta})}{\sigma^{2}_{t}}\right\\}\left\llangle\frac{\partial
f(\bm{x}^{t}+\bm{\delta})}{\partial\bm{\delta}}\right\rrangle.$ (38)

Note that $f(\bm{x}^{t}+\bm{\delta})$ is readily available at any
$\bm{\delta}$ without approximation.

Here we provide some intuition behind the updating equation (38). Convergence
is achieved when either the deviation $y^{t}-f$ or the gradient
$\llangle\partial f/\partial\bm{\delta}\rrangle$ vanishes at
$\bm{x}^{t}+\bm{\delta}$. These situations are illustrated in Fig. 2 (b). As
shown in the figure, $\delta_{i}$ corresponds to the “horizontal deviation”
along the $x_{i}$ axis between the test sample and the regression function. If
there is no horizontal intersection on the regression surface it seeks the
zero gradient point based on a smooth surrogate of the gradient.

###  5.4 Estimating smooth gradient

The final step is to estimate the smooth surrogate of the gradient
$\llangle\partial f/\partial\bm{\delta}\rrangle$ in the vicinity of
$\bm{x}_{\delta}\triangleq\bm{x}^{t}+\bm{\delta}$. To handle the potential
non-differentiability of $f$, we define the gradient as the local mean of the
slope function $[f(\bm{x}_{\delta}+h\bm{e}_{i})-f(\bm{x}_{\delta})]/h$, where
$h$ is a small perturbation and $\bm{e}_{i}$ is a unit vector which assumes
value 1 in the $i$-th entry and 0 otherwise. Let $p(h\mid\bm{x}_{\delta})$ be
an assumed local distribution for $h$ around $\bm{x}_{\delta}$. Now we define
the local gradient as

$\displaystyle\left\llangle\frac{\partial
f(\bm{x}_{\delta})}{\partial\delta_{i}}\right\rrangle$
$\displaystyle\triangleq\int\mathrm{d}h\
p(h\mid\bm{x}_{\delta})\frac{f(\bm{x}_{\delta}+h\bm{e}_{i})-f(\bm{x}_{\delta})}{h}\approx\frac{1}{N_{s}}\sum_{n=1}^{N_{s}}\frac{f(\bm{x}_{\delta}+h^{[n]}\bm{e}_{i})-f(\bm{x}_{\delta})}{h^{[n]}},$
(39)

where $h^{[n]}$ is the $n$-the sample from $p(h\mid\bm{x}_{\delta})$ and
$N_{s}$ is the number of perturbations generated. The second approximate
equality is due to Monte Carlo approximation, which is guaranteed to converge
as $N_{s}\to\infty$. One reasonable choice for the local distribution is
$p(h\mid\bm{x}_{\delta})=\mathcal{N}(h\mid\bm{x}_{\delta},\eta^{2})$ with
$\eta^{2}$ being the standard deviation of the perturbation. In this case,
perturbations that happen to be zero numerically have to be excluded from the
computation.

#### Alternative approaches to smooth gradient estimation

It is worth noting that local gradient estimation has been studied in
evolutionary computation for years [Salomon, 1998, Salomon and Arnold, 2009].
The key idea is to leverage the notion of Gaussian smoothing of a potentially
non-continuous function:

$\displaystyle f_{\eta}(\bm{x}_{\delta})$
$\displaystyle\approx\int\mathrm{d}\bm{h}\
\mathcal{N}(\bm{h}\mid\bm{0},\eta^{2}\mathsf{I}_{M})f(\bm{x}_{\delta}+\bm{h}),$
(40)

where $\mathsf{I}_{M}$ is the $M$-dimensional identity matrix. As $f_{\eta}$
can be viewed as a locally smoothed version of $f$, the gradient of $f_{\eta}$
is a reasonable estimate of $\llangle\partial f/\partial\bm{\delta}\rrangle$.
By using integration by parts and Monte Carlo estimation, the Gaussian
smoothing approach gives

$\displaystyle\frac{\partial f_{\eta}(\bm{x}_{\delta})}{\partial
x_{i}}\approx\frac{1}{N_{s}}\sum_{n=1}^{N_{s}}\frac{h_{i}^{[n]}}{\eta^{2}}\left[f(\bm{x}_{\delta}+h_{i}^{[n]}\bm{e}_{i})-\bar{f}(\bm{x}_{\delta})\right],\quad\bar{f}(\bm{x}_{\delta})\triangleq\frac{1}{N_{s}}\sum_{n}f(\bm{x}_{\delta}+\bm{e}_{i}h^{[n]}),$
(41)

which is another reasonable estimate of the local gradient. From the
perspective of numerical computation, however, there is no compelling reason
to use this expression instead of (39) because (41) tends to have a much
larger variance than (39).

Another reasonable approach is to locally fit a linear function and use the
coefficients as a surrogate of the gradient, as proposed by [Idé et al.,
2021]. This provides an estimate as accurately as the direct slope estimation
approach (39) does. However, one issue is that the resulting estimation
formula is not linear in $f$ and hence is not eligible for fast vectorized
computation. This can be problematic when repeated gradient estimation is
required.

Algorithm 2 Likelihood Compensation

1:Black-box regression model $f(\bm{x})$, test data
$\mathcal{D}_{\mathrm{test}}$, and parameters $\lambda,\nu,\kappa$.

2:for all $\bm{x}^{t}\in\mathcal{D}_{\mathrm{test}}$ do

3: Compute $\sigma^{2}_{t}$ with Eq. (33).

4:end for

5:Randomly initialize $\bm{\delta}\approx\bm{0}$.

6:repeat

7: Set $\bm{g}=\bm{0}$.

8: for all $\bm{x}^{t}\in\mathcal{D}_{\mathrm{test}}$ do

9: Compute $\left\llangle\frac{\partial
f(\bm{x}^{t}+\bm{\delta})}{\partial\bm{\delta}}\right\rrangle$ with Eq. (39)

10: Update $\bm{g}\leftarrow\bm{g}+\left\llangle\frac{\partial
f(\bm{x}^{t}+\bm{\delta})}{\partial\bm{\delta}}\right\rrangle\frac{y^{t}-f(\bm{x}^{t}+\bm{\delta})}{N_{\mathrm{test}}\sigma^{2}_{t}}$.

11: end for

12: $\bm{\phi}=(1-\kappa\lambda)\bm{\delta}+\kappa\bm{g}$.

13: Find $\bm{\delta}$ with Eq. (37).

14:until convergence.

15:return $\bm{\delta}$

###  5.5 Algorithm Summary

Algorithm 2 summarizes the iterative procedure for finding $\bm{\delta}$. The
most important parameter is the $\ell_{1}$ regularization strength $\nu$,
which has to be hand-tuned depending on the business requirements of the
application of interest. On the other hand, the $\ell_{2}$ strength $\lambda$
controls the overall scale of $\bm{\delta}$. It can be fixed to some value
between 0 and 1. In our experiments, it was adjusted so its scale is on the
same order as LIME’s output for consistency. It is generally recommended to
rescale the input variables to have zero mean and unit variance before
starting the iteration (assuming $N_{\mathrm{test}}\gg 1$), and retrieve the
scale factors after convergence. For the learning rate $\kappa$, in our
experiments, we fixed $\kappa=0.1$ and shrank it (geometrically) by a factor
of 0.98 in every iteration.

In addition to the parameters listed in Algorithm 2, gradient estimation by
Eq. (41) requires two minor parameters, $N^{\mathrm{s}},\eta$. In our
experiments, we used $N^{\mathrm{s}}=10$, which was confirmed to provide
sufficient convergence.

##  6 Experiments

This section presents empirical evaluation of the proposed anomaly attribution
framework. For comprehensive coverage, we use five datasets with different
statistical complexities, including one from a real business use case. The
goals of this evaluation are to 1) provide a clear picture of what deviation-
sensitivity of an attribution method buys us using a simple synthetic model;
2) demonstrate LC’s capability of providing directly interpretable attribution
scores; 3) point out inherent issues with the IG family in anomaly
attribution; 4) quantitatively analyze the consistency and inconsistency among
different attribution methods; and 5) demonstrate how LC was able to make a
difference in a real business scenario. For the reader’s convenience, we
summarize the datasets we used in Table 3.

Table 3: Summary of the datasets used. ‘NA’ denotes ‘not available.’ | $N_{\mathrm{train}}$ | $N_{\mathrm{test}}$ | $M$ | $f(\bm{x})$ | characteristics
---|---|---|---|---|---
2D sinusoidal | $\infty$ | 1 | 2 | analytic | smooth, periodic, noise-free
Boston Housing | 506 | 1 | 13 | RF | multimodal, clustered, partially discontinuous
California Housing | $20\,640$ | 1 | 9 | GBT | large, uni- or bimodal, partially discontinuous
Diabetes | 442 | 1 | 10 | DNN | unimodal, semi-discrete
Building Energy | NA | 24 | 12 | commercial | noisy, periodic, online, real building HVAC

###  6.1 Baselines

We compare LC with five possible alternatives111The Python implementation will
be made available upon acceptance of the paper.: LIME [Ribeiro et al., 2016],
SV [Štrumbelj and Kononenko, 2014], IG [Sipple, 2020], EIG as defined in Eq.
(5), and the $Z$-score, as summarized in Table 2. For anomaly attribution,
LIME, SV, IG, and EIG are applied to the deviation $f(\bm{x})-y$ rather than
$f(\bm{x})$. The $Z$-score is one of the standard univariate outlier detection
metrics in the unsupervised setting, and defined as
$Z_{i}\triangleq(x_{i}^{t}-m_{i})/\sigma_{i}$ for the $i$-th variable, where
$m_{i},\sigma_{i}$ are the mean and the standard deviation of $x_{i}$,
respectively. In SV, we used the same sampling scheme as that proposed in
[Štrumbelj and Kononenko, 2014] to handle combinatorial complexity. In IG and
EIG, we used the trapezoidal rule with 100 equally-spaced intervals to perform
the integration w.r.t. $\alpha$. In IG, EIG, and LC, we used the same gradient
estimation algorithm in Eq. (39).

We listed SV, EIG, and the $Z$-score here for comparison purposes despite the
fact that they are not actually applicable in our doubly black-box setting
(see Table 2). We excluded other contrastive and counterfactual methods such
as [Wachter et al., 2017] as they require white-box access to the model and
are predominantly used in classification settings. In the real-world case
study in Sec. 6.6, we validated our approach with feedback from domain experts
as opposed to crowd-sourced user studies with lay users. In industrial
applications, the end-user’s interests can be highly specific to the
particular business needs and the system’s inner workings tend to be difficult
for non-experts to understand and simulate.

Figure 3: 2D Sinusoidal Curve with the $x_{2}=0$ slice on the top. The points
A and B are at $y^{t}=1$ and $-1$, respectively, while they are at the same
$\bm{x}^{t}=(0.5,0)$.

###  6.2 Deviation-sensitivity

#### 2D Sinusoidal Curve data

To illustrate LC’s deviation-sensitive property, we computed attribution
scores for a regression curve defined by a two-dimensional (2D) sinusoidal
function

$\displaystyle f(\bm{x})$ $\displaystyle=2\cos(\pi x_{1})\sin(\pi x_{2}).$
(42)

We defined two test points: A is at $(\bm{x}^{t},y^{t})=((0.5,0),+1)$ and B is
at $((0.5,0),-1)$, as shown in Fig. 3. Computation of SV, EIG, and the
$Z$-score needs the true distribution $P(\bm{x})$. We used the empirical
approximation for $P(\bm{x})$ by randomly generating samples with the uniform
distribution in $[-4,4]^{2}$, which was set to be wide enough to simulate
sparse sample distributions, as opposed to Gaussian-like unimodal
distributions. We created 10 held-out datasets, each of which consists of
$N_{\mathrm{ho}}=100$ samples. Those data sets were used also to compute the
mean and standard deviation for the $Z$-score. For IG, we gave two baseline
inputs: $(0,0)$ and $(0,1)$, resulting in two IG scores for each dataset. The
regularization parameters are set to be negligible values as regularization is
unimportant in this low-dimensional setting.

Figure 4: Comparison of attribution scores on the 2D Sinusoidal Curve at two
test points (A and B in Fig. 3). The scores were evaluated 10 times over
randomly generated datasets. Only the LC scores differ between A and B.

Figure 4 compares the attribution score, where the mean and the standard
deviation over the ten trials is shown. Randomness in LC, LIME and IGs comes
only from gradient estimation and is negligible, while EIG, SV and the
$Z$-score are directly impacted by the variability of the samples. The
attribution scores are normalized by dividing by a scaling factor to make them
mutually comparable. The most conspicuous observation from Fig. 4 is that only
LC can distinguish the direction of the deviation between A and B; All the
other methods produce exactly the same score between A and B due to the
deviation-agnostic property. As seen from Fig. 3, $\delta_{1}$ (the LC score
for $x_{1}$) was negative for point A. This can be easily understood from the
definition of LC as ‘horizontal deviation.’ For the given $y^{t}=1$ value, the
location of the point is unusual; it should have been a little more to the
left (so it gets the maximum reward of likelihood). Similarly, for point B,
$x_{1}=0.5$ was unusual for $y^{t}=-1$ and it should have been a little more
to the right to match $y^{t}=-1$.

We see that SV, EIG, and $Z$-scores have significant variability in contrast
to LC, LIME and IGs. As discussed in Sec. 4.1, EIG is reduced to IG when the
probability mass of $P(\bm{x})$ is concentrated at the baseline input of IG.
In this case, however, the distribution has been chosen to be a broad uniform
distribution, which is at the opposite end of the spectrum from Dirac’s delta
function. In such a case, EIG generally has a large variability. For the cases
where $P(\bm{x})$ is close to unimodal, see Sec. 6.4.

We also see that the attribution scores of IG significantly vary depending on
the choice of the baseline input, which is either $(0,0)$ or $(0,1)$. EIG
eliminates the need for arbitrary input by expectation. However, empirical
expectation resulted in large error bars. These characteristics of the
baseline methods, along with their deviation-agnostic property, make LC the
preferred method for the task of anomaly attribution.

###  6.3 Direct interpretability of LC

#### Boston Housing data

As an example of a real-world attribution task, we next use Boston Housing
data [Belsley et al., 2005], a well-known benchmark data for regression. The
task is to predict $y$, the median home price (‘MEDV’) of the districts in
Boston, with $\bm{x}$, the input vector of size $M=13$, such as the percentage
of the lower status of the population (‘LSTAT’) and the average number of
rooms (‘RM’)222We excluded a variable named ‘B’ from attribution for ethical
concerns [scikit-learn 1.1.3, 2022].. The total number of samples is 506. As
one might expect, the data is very noisy. Figure 5 shows scatter plots between
$y$ and six selected input variables, where multimodal, clustered structures
are observed. We standardized the entire dataset for each variable to be zero
mean and unit variance, then we held out $20$% of the data as
$\mathcal{D}_{\mathrm{test}}$, and trained the random forest (RF) [Hastie et
al., 2009] on the rest. Viewing it as a black-box regression model
$f(\bm{x})$, we computed the anomaly score with Eq. (1), based on
$p(y\mid\bm{x})$ estimated with Eqs. (30)-(33), where $w_{n}$ is set to a
constant on the held-out samples. Note that the training dataset is supposed
to be unavailable in the doubly black-box setting. We explicitly use the
training data here for comparison purposes.

Figure 5: Boston Housing: Pairwise scatter plot between $y$ (MEDV) and
selected input variables. The square and triangle show the detected first and
second outliers in the test data.

Figure 6 (a) shows the computed anomaly scores. The first and second outliers
are highlighted with the square and triangle symbols, respectively, which have
also been shown in Fig. 5 with the same symbols. Unlike the general
expectation from the word ‘outlier,’ those anomalous samples are not
necessarily in the low-density region in the scatter plot. One reason for this
is that the pairwise scatter plot visualize $p(y^{t}\mid x^{t}_{i})$ for each
$i$ rather than $p(y^{t}\mid\bm{x}^{t})$, on which the anomaly score was
calculated. It is well-known that the marginal distribution makes
interpretation tricky when deviations are important [Molnar, 2022, Kumar et
al., 2020]. We strive to get a more sophisticated explanation specific to a
test sample rather than aggregated signals in the marginals.

Figure 6: Anomaly score for (a) Boston Housing, (b) California Housing, and
(c) Diabetes datasets. Top outliers are highlighted.

Figure 7 shows attribution scores computed for the first and second outliers.
We used $\lambda=0.5,\nu=0.1,\kappa=0.1$ for LC, which were adjusted so the
entropy of the absolute score distribution is roughly the same as that of IG.
The same regularization parameter was used for LIME. As discussed in Sec. 4,
IG, EIG, and SV are pseudo-local methods in the sense that they require either
a global distribution or a baseline input outside the vicinity of the test
input, while LC and LIME are purely local attribution framework. To make LC
and LIME comparable to the others, we used a relatively large $\eta=1$ for
gradient estimation. For IG, the baseline input was equated to the mean of
$\bm{x}$. In Fig. 7 (1), all the methods except for the $Z$-score gave similar
score distributions, where LSTAT and RM have dominating weights. Close
inspection shows, however, that there are interesting differences in their
signs. In the gradient-based approach LIME, LSTAT and RM give negative and
positive signs, respectively, following the global trend in the scatter plots
in Fig. 5. As discussed in the previous subsection, LIME does not take account
of $y$ of the test point. If a variable has a globally monotonic distribution,
LIME tends to provide a similar score, regardless of the $y$ value. In
contrast, LC’s attribution score is directly interpretable. In this case, LC
was negative for RM because a negative shift would give a better fit: “The
number of rooms is a bit too large for the (relatively low) price.” In
addition, LC’s attribution score represents an actual shift. We can readily
confirm that the red square would fall into the densest region of the RM-MEDV
scatter plot with a shift by about 0.5 to the left in Fig. 5.

Figure 7: Boston Housing: Comparison of the attribution scores for (1) the
first and (2) the second outliers highlighted with the square (red) and
triangle (blue) symbols, respectively, in Fig. 5.

#### California Housing data

To confirm general applicability of LC, we next used the California Housing
dataset [Pace and Barry, 1997]. This is a relatively large data set of
$20\,640$ samples with 9 predictor variables, of which we log-transformed four
variables (‘AveRooms’, ‘AveBedrms’, ‘Population’, ‘AveOccup’). The task is to
predict the median house value of small geographical segments using predictor
variables such as the median household income of each segment. We randomly
held out 20% of the samples after standardization and trained gradient boosted
trees (GBT) [Friedman, 2002] on the rest. Following the same procedure as that
for the Boston Housing dataset, we computed anomaly scores for the held-out
samples as shown in Fig. 6 (b). Then we took the one with the highest anomaly
score (highlighted with the red rectangle) as the test sample
$(\bm{x}^{t},y^{t})$.

Figure 8 shows pairwise scatter plots of selected variables against the target
variable (‘MedHouseVal’). An interesting bi-modal structure is observed in the
two variables. As clearly seen from the figure, the test sample is off the
main cluster in many variables. From the figure, we see, for example, that the
median income of the segment (‘MedInc’) is a bit too small and the average
number of household members (‘AveOccup’) is a bit too large. One interesting
question is whether these observations are consistent with LC’s attribution.

To answer this question, we compare computed attribution scores in Fig. 9. We
used $\lambda=0.4,\nu=0.2,\kappa=0.1$ for LC, which were adjusted using the
same approach as the Boston case. The same regularization parameter was used
for LIME. IG’s baseline input was set to be the origin (the population mean of
$\bm{x}$). For EIG, SV, and the $Z$-score, we used empirical approximation
using 100 bootstrapped samples from the training data to simulate a ‘semi-
doubly black-box’ situation, where only a limited number of test samples are
available (see Sec. 6.4 for more detail). The bootstrap approach allows
estimating the distribution of attribution scores (see the next subsection for
the detail). The $Z$-score confirms that MedInc and AveOccup are significantly
smaller and larger than the mean, respectively. However, the attribution
scores for the latter are almost zero in all the five methods. This is a clear
example where the deviation in $\bm{x}$ does not necessarily mean being an
outlier in regression. Specifically, since AveOccup is almost unrelated with
the target variable, as suggested by the distribution in Fig. 8, any shifts in
that variable should not decisively explain the anomalousness.

For MedInc, on the other hand, LC pinpoints that it is the biggest contributor
and a large positive attribution score provides a counterfactual explanation:
If MedInc were a bit higher, the sample would have looked less unusual. In
other words, the median income was a bit too low for the median house value.

Figure 8: California Housing: Pairwise scatter plot between $y$
(‘MedHouseVal’) and selected input variables. The square symbol indicates the
top outlier in Fig. 6 (b).  Figure 9: California Housing: Comparison of
attribution scores for the outlier highlighted with the red square in Fig. 8.

###  6.4 Variability of EIG and SV

#### Boston Housing

Let us go back to Figs. 7 (1) and (2) of Boston Housing. Between these two
outliers, LC, LIME, and IG exhibit consistent scores. Interestingly, this is
not the case in SV and EIG: IG’s large weight on RM vanishes in EIG and SV,
resulting in significantly different attributions. Since the baseline input of
IG was set to be the mean of $\bm{x}$, one would expect largely consistent
attributions between IG and EIG (and hence, SV, according to Theorem 4).
Recall Fig. 1 (b) and the sum rule in Eq. (10). In this particular case, IG
has a large positive weight on RM because the mean is on the left to the test
points, and, as going from the mean to the test points, $f(\bm{x})$ value goes
up along the globally positive gradient. In (2), the test point is closer to
the mean than (1), as indicated by their $Z$-scores; Its RM value is located
in the densest region of RM. In such a case, the contribution of the variable
tends to be unpredictable because the paths from a test point can be in almost
any direction, even from the right, and averaging out those paths can result
in either a large or a near-zero value. In general, the attribution of EIG can
be counter-intuitive when the test input is close to the population mean of
$\bm{x}$ and the local gradient is steep.

This “vanishing weight” issue originates mainly from the fact that IG and EIG
(and thus, SV) explain the increment rather than the deviation. We see the
test points as outliers Fig. 7 (2) because of the large deviations in
$\bm{y}$, i.e. the vertical shifts. However, the deviation-agnostic methods do
not see the test points in that way. This relatively simple example
empirically demonstrates the subtle but intricate difficulty in using
deviation-agnostic algorithms for anomaly attribution.

Figure 10: Boston Housing: Comparison of attribution scores with simulated
score variability for the outlier highlighted with the blue triangle in Fig.
5.

As mentioned earlier, EG, SV, and the $Z$-score are basically inapplicable to
anomaly attribution in the doubly black-box setting. If there is some amount
of samples, however, one can compute these quantities somehow by taking the
empirical average, as we have done above. One question here is whether those
methods can produce stable attribution scores under a limited number of
samples. To simulate such a situation, we created 10 sets of bootstrapped
datasets of $N_{\mathrm{b}}=100$, generated from the training set of the
Boston Housing data. Between the two outliers in Fig. 7, we focused on the
second one as it has much more diversity in the score. The expectation in EIG
and SV was performed using the empirical approximation by the bootstrapped
samples. The distribution of the attribution scores are shown in Fig. 10 for
EIG, SV, and the $Z$-score, where the points denote computed score values of
the 10 bootstrap rounds. The curves were estimated Gaussian-based kernel
density estimation333We used the KernelDensity implementation of scikit-learn
1.2.0.

$\displaystyle
p(s_{i})=\frac{1}{N_{\mathrm{b}}}\sum_{l=1}^{N_{\mathrm{B}}}\mathcal{N}(s_{i}\mid
s_{i}^{(l)},\eta_{\mathrm{b}}^{2}),$ (43)

where $s_{i}$ denotes the attribution score of the $i$-th variable,
$s_{i}^{(l)}$ denotes the computed $s_{i}$ value on the $l$-th bootstrapped
replica (dataset), and $N_{\mathrm{B}}$ is the number of bootstrapped
replicas, which is 10 in this case. Also, the bandwidth $\eta_{\mathrm{b}}$
was set to 4% of the range for each variable. As LC, LIME and IG do not need
$P(\bm{x})$, their scores are the same as those presented in Fig. 7. The same
approach was used to draw the curves in Fig. 9.

From Fig. 10, we first see that there is a systematic similarity between SV’s
and EIG’s scores. This is an empirical confirmation of Theorem 4. As SV is
computationally demanding in high-dimensional data due to its combinatorial
nature, one can use EIG to approximate SV in practice. Given the vanishing
weight issue of EIG, even IG could be used as long as prior knowledge
naturally defines the baseline input, although these methods are deviation-
agnostic and not the best choice for anomaly attribution. Second, in SV and
EIG, the score distributes around zero in many variables. This underscores the
need for a sparsity-enforcing mechanism in attribution. In contrast, LC and
LIME are not distracted by pseudo-signals thanks to the sparsity in their
scores. Third, the variability of the score can be extremely large in some
variables. In fact, in RM, the standard deviation of the distribution is even
larger than the absolute score itself. This means that RM’s attribution score
is not at all trustworthy. The large variance issue in RM is another
manifestation of the vanishing weight issue discussed in the previous
subsection. The above observations remained unchanged when we increase
$N_{\mathrm{b}}$ e.g. to 200.

###  6.5 Consistency among attribution methods

Figures 9 and 10 indicate that LC and the alternative attribution methods are
largely consistent at least in the top attributions. To quantitatively
evaluate the consistency among different attribution methods, we computed the
following four metrics. The first and second metric is Kendall’s $\tau$ and
Spearman’s $\rho$ computed on two absolute attribution score vectors. They are
typically called the rank correlation coefficients and take a value of 1 if
the order of the two absolute scores are the same regardless of their values.
The third metric is what we call the sign match ratio, which takes on 1 when
all the signs are consistent between corresponding vector elements. When
comparing an attribution score vector $\bm{u}$ against a reference score
vector $\bm{r}$, the sign match ratio is defined as

$\displaystyle(\mbox{sign match ratio})$ $\displaystyle\triangleq
1-\frac{1}{M}\sum_{i=1}^{M}\mathbb{I}\left(\operatorname{sign}(r_{i})\operatorname{sign}(u_{i})=-1\right),$
(44)

where $\mathbb{I}(\cdot)$ is the indicator function that takes on 1 when the
argument is true, 0 otherwise. We define $\operatorname{sign}(0)=0$ in this
case. Note that this favors sparse attribution scores: If $\bm{r}=\bm{0}$,
then the score is always 1 regardless of $\bm{u}$. Finally, the fourth metric
is what we call hit25, which gives 1 when the top 25% of the absolute entries
perfectly match between $\bm{r}$ and $\bm{u}$, and 0 if none of the top 25%
members of $\bm{r}$ is included in that of $\bm{u}$. As hit25 depends on
neither the sign nor the rank, it quantifies simply the match of top
contributors.

Figure 11: Diabetes: Pairwise scatter plot between $y$ (‘progression’) and
selected input variables. The square highlights the detected top outlier in
Fig. 6 (c).  Figure 12: Diabetes: Comparison of attribution scores with
simulated score variability for the top outlier highlighted with the square in
Fig. 11.

#### Diabetes data

To study the consistency among the attribution methods in a more comprehensive
fashion, we added another benchmark dataset called the Diabetes dataset [Efron
et al., 2004], which contains 442 samples of one real-valued target variable
(‘progression’) and $M=10$ predictors including the body-mass index (‘bmi’),
the average blood pressure (‘bp’), and eight other biomarkers. For this
dataset, we held-out 20% of samples after the mini-max scaling with in [0,1]
and trained a deep neural network (DNN) on the rest. The resulting model
identified has two hidden layers of 32 and 8 neurons with the ReLU (rectified
linear unit) activation. The average test $R^{2}$ score was $0.54$. We thought
of the trained DNN as a black-box regression function $y=f(\bm{x})$. Following
the same procedure as the other datasets, we computed the anomaly score
presented in Fig. 6 (c), where the top outlier is highlighted with the red
rectangle. Figure 11 shows where the outlier point is located in the pairwise
scatter plots against the target variable.

For this outlier as the test point $(\bm{x}^{t},y^{t})$, we computed
attribution scores as shown in Fig. 12. We used the same parameters as those
for the California Housing experiment. The baseline input for IG is also at
the origin although the origin is no longer the population mean due to the
min-max scaling. In this particular test point, LIME, IG, EIG, and SV look
quite consistent. LIME simply captures the overall positive trend. IG also
captures the positive increment when going from zero (the baseline input in
this case) to the test point except for s2 (and s1, which is not shown). As
the population mean is close to zero, EIG and SV follow a similar trend.

#### Comprehensive consistency analysis

Table 4: Result of consistency analysis. The mean and the standard deviation of the metrics are shown in each cell, where 1 represents the highest consistency.  |  | LC | LIME | IG | EIG | SV |  $Z$-score
---|---|---|---|---|---|---|---
Boston | $\tau$ | $1.00\pm 0.00$ | $0.49\pm 0.31$ | $0.61\pm 0.09$ | $0.40\pm 0.07$ | $0.53\pm 0.16$ | $0.15\pm 0.34$
$\rho$ | $1.00\pm 0.00$ | $0.59\pm 0.32$ | $0.72\pm 0.07$ | $0.49\pm 0.09$ | $0.61\pm 0.19$ | $0.18\pm 0.39$
sign | $1.00\pm 0.00$ | $0.68\pm 0.10$ | $0.83\pm 0.08$ | $0.80\pm 0.09$ | $0.86\pm 0.06$ | $0.74\pm 0.12$
hit25 | $1.00\pm 0.00$ | $0.80\pm 0.18$ | $0.80\pm 0.30$ | $0.60\pm 0.15$ | $0.67\pm 0.24$ | $0.27\pm 0.28$
California | $\tau$ | $1.00\pm 0.00$ | $0.87\pm 0.04$ | $0.75\pm 0.10$ | $0.50\pm 0.12$ | $0.52\pm 0.13$ | $0.08\pm 0.15$
$\rho$ | $1.00\pm 0.00$ | $0.93\pm 0.03$ | $0.89\pm 0.04$ | $0.65\pm 0.15$ | $0.67\pm 0.15$ | $0.16\pm 0.21$
sign | $1.00\pm 0.00$ | $0.48\pm 0.04$ | $0.48\pm 0.04$ | $0.62\pm 0.13$ | $0.66\pm 0.15$ | $0.66\pm 0.05$
hit25 | $1.00\pm 0.00$ | $1.00\pm 0.00$ | $0.90\pm 0.22$ | $0.80\pm 0.27$ | $0.80\pm 0.27$ | $0.30\pm 0.44$
Diabetes | $\tau$ | $1.00\pm 0.00$ | $0.67\pm 0.08$ | $0.53\pm 0.12$ | $0.52\pm 0.21$ | $0.54\pm 0.19$ | $-0.04\pm 0.19$
$\rho$ | $1.00\pm 0.00$ | $0.78\pm 0.07$ | $0.67\pm 0.13$ | $0.66\pm 0.17$ | $0.66\pm 0.19$ | $-0.09\pm 0.27$
sign | $1.00\pm 0.00$ | $0.85\pm 0.22$ | $0.68\pm 0.11$ | $0.65\pm 0.14$ | $0.65\pm 0.14$ | $0.65\pm 0.16$
hit25 | $1.00\pm 0.00$ | $0.90\pm 0.22$ | $0.80\pm 0.27$ | $0.80\pm 0.27$ | $0.80\pm 0.27$ | $0.20\pm 0.27$

On the other hand, LC looks like the LIME attribution scores with the opposite
sign in Fig. 12, which was also the case in Fig. 10. One interesting question
here is whether or not there is a systematic correspondence of LC’s
attribution to the other methods, especially when we ignore the signs. Among
the four metrics, $\tau,\rho$, and hit25 are about the absolute scores. In the
three data sets, we picked top five outliers and computed the attribution
scores for them. For each, we computed the four metrics with the reference
$\bm{r}$ being the LC’s score. Those metrics are designed to capture the
consistency between the top contributors, disregarding minor contributors. LC
is the only method that has both a built-in sparsity-enforcing mechanism and
the deviation-sensitive property, making it an appropriate choice for the
reference.

The result is summarized in Table 4. As expected, hit25 has generally high
scores, apart from the $Z$-score. This suggests that those attribution methods
are a useful tool to select important features. Even in the other metrics
including the sign match ratio, they produce reasonably consistent
attributions in some cases. However, some 20-30% of cases are still not
necessarily consistent, which is a natural consequence that LC is deviation-
sensitive but the others are not.

###  6.6 Real-world business use-case: building energy management

Finally, we provide an example about how LC can make a difference in a real
business use-case. Collaborating with IBM IoT Business Unit, we obtained
energy consumption data for an office building in India. The total wattage $y$
is predicted by a black-box commercial prediction tool as a function of
weather-related (temperature, humidity, etc.) and time-related variables (time
of day, day of week, month, etc.). There are two intended usages of the
predictive model. One is near future prediction with short time windows for
optimizing HVAC (heating, ventilating, and air conditioning) system control.
The other is retrospective analysis over the last few months for the purpose
of planning long-term improvement of the building facility and its management
policies. In the retrospective analysis, it is critical to get clear
explanation on unusual events.

At the beginning of the project, we interviewed 10 professionals on what kind
of model explainability would be most useful for them. Their top priority
capabilities were uncertainty quantification in forecasting and anomaly
diagnosis in retrospective analysis. Our choices in the current research
reflect these business requirements.

Figure 13: Building Energy: Outlier score computed with Eq. (1) for the test
data.

#### Anomaly detection

We obtained a one month worth of test data with $M=12$ input variables
recorded hourly. We first validated the assumed Gaussian observation model by
performing kernel density estimation for the probability density function
(p.d.f.) of $y$. As shown in Fig. 14, the p.d.f. of $y$ itself is double-
peaked, corresponding to different consumption patterns between night and day.
On the other hand, the p.d.f. of $y^{t}-f(\bm{x}^{t})$ in the right panel is
single-peaked, which confirms the validity of the Gaussian-based model in Eq.
(26).

Figure 14: Building Energy: Estimated probability densities. Although raw
wattage value $y$ (left) is far from Gaussian, the deviation $y-f(\bm{x})$ is
well approximated by the Gaussian.

Next, we computed the anomaly score by Eq. (1) for each sample
$(y^{t},\bm{x}^{t})$ under the Gaussian observation model. The variance
$\sigma^{2}_{t}$ was computed using Eq. (33) for each $t$ by leaving
$(y^{t},\bm{x}^{t})$ out from the dataset. The computed anomaly score is shown
in Fig. 13. We see that there are two periods showing conspicuous anomalies,
namely, October 2 and 18-20. An important business question was who or what
may be responsible for these anomalies.

Figure 15: Building Energy: Comparison of the explainability scores computed
for the test data.

#### Anomaly attribution

To obtain insights regarding the detected anomalies, we computed the LC score
as shown in Fig. 15, where we computed $\bm{\delta}$ each day with
$N_{\mathrm{test}}=24$ in Eq. (27), and visualized $\|\bm{\delta}\|_{2}^{2}$.
We also compared two alternative methods that were applicable: For the
$Z$-score, we visualized the daily mean of the absolute values. For LIME, we
computed regression coefficients for every sample, and visualized the
$\ell_{2}$ norm of their daily mean. We used $(\nu,\lambda)=(0.1,0.5)$, which
was determined by the level of sparsity and scale preferred by the domain
experts. IG, EIG, and SV need either a baseline input or training data, and
were inapplicable to this real-world setting.

As shown in the plot, the LC score clearly highlights a few variables whenever
the outlier score is exceptionally high in Fig. 13, while the $Z$-score and
LIME do not provide much information beyond the trivial weekly patterns. The
pattern of LIME was very stable over $0<\nu\leq 1$, showing empirical evidence
of the deviation-agnostic property. On the other hand, the $Z$-score
sensitively captures the variability in the weather-related variables, but it
fails to explain the deviations in Fig. 13. This is understandable because the
$Z$-score does not reflect the relationship between $y$ and $\bm{x}$. The
artifact seen in the “daytype” variables is due to the one-hot encoding of the
day of week.

With LC, the strongest signal is observed around October 19 (Thursday) in Fig.
13. The variables highlighted are ‘timeofday’, ‘daytype_Sa’, and ‘daytype_Su’,
implying that those days had an unusual daily wattage pattern for a weekday
and looked more like weekend days. Interestingly, it turned out that the 19th
was a national holiday in India and many workers were off on and around that
date. The other anomalous period on October 2 was also a national holiday.
Thus we conclude that the anomaly is most likely not due to any faulty
building facility, but due to the model limitation caused by the lack of full
calendar information. Though simple, such pointed insights made possible by
our method were highly appreciated by the professionals.444LC has been
productized as IBM’s software offering.

##  7 Conclusions

We have proposed a new framework for model-agnostic anomaly attribution in the
doubly black-box regression setting. We mathematically proved that integrated
gradient IG), local linear surrogate modeling (LIME), and Shapley values (SV)
are inherently deviation-agnostic and thus, cannot be a viable solution for
anomaly attribution. We have clarified a mathematical structure leading to the
deviation-agnostic property using a power expansion technique. Unlike these
methods, the proposed likelihood compensation approach is built upon the
maximum likelihood principle, and is capable of capturing specific
characteristics of anomalies observed. We conducted a comprehensive empirical
study using benchmark datasets to verify our mathematical characterization
such as an equivalence between SV and IG. We also validated the proposed
method on a real-world use-case of building energy management where—based on
expert feedback received—the proposed LC method offers significant practical
advantages over existing methods.

## References

  * [Abhishek and Kamath, 2022] Abhishek, K. and Kamath, D. (2022).  Attribution-based XAI methods in computer vision: A review.  arXiv preprint arXiv:2211.14736.
  * [Antwarg et al., 2021] Antwarg, L., Miller, R. M., Shapira, B., and Rokach, L. (2021).  Explaining anomalies detected by autoencoders using shapley additive explanations.  Expert Systems with Applications, 186:115736.
  * [Arrieta et al., 2020] Arrieta, A. B., Díaz-Rodríguez, N., Del Ser, J., Bennetot, A., Tabik, S., Barbado, A., García, S., Gil-López, S., Molina, D., Benjamins, R., et al. (2020).  Explainable artificial intelligence (XAI): Concepts, taxonomies, opportunities and challenges toward responsible ai.  Information Fusion, 58:82–115.
  * [Belsley et al., 2005] Belsley, D. A., Kuh, E., and Welsch, R. E. (2005).  Regression diagnostics: Identifying influential data and sources of collinearity, volume 571.  John Wiley & Sons.
  * [Brereton, 2014] Brereton, R. G. (2014).  The normal distribution.  Journal of Chemometrics, 28(11):789–792.
  * [Burkart and Huber, 2021] Burkart, N. and Huber, M. F. (2021).  A survey on the explainability of supervised machine learning.  Journal of Artificial Intelligence Research, 70:245–317.
  * [Chandola et al., 2009] Chandola, V., Banerjee, A., and Kumar, V. (2009).  Anomaly detection: A survey.  ACM Computing Survey, 41(3):1–58.
  * [Dang et al., 2013a] Dang, X. H., Micenková, B., Assent, I., and Ng, R. T. (2013a).  Local outlier detection with interpretation.  In Joint European Conference on Machine Learning and Knowledge Discovery in Databases, pages 304–320. Springer.
  * [Dang et al., 2013b] Dang, X. H., Micenková, B., Assent, I., and Ng, R. T. (2013b).  Outlier detection with space transformation and spectral analysis.  In Proceedings of the 2013 SIAM International Conference on Data Mining, pages 225–233. SIAM.
  * [Deng et al., 2021] Deng, H., Zou, N., Du, M., Chen, W., Feng, G., and Hu, X. (2021).  A unified taylor framework for revisiting attribution methods.  In Proceedings of the AAAI Conference on Artificial Intelligence, pages 11462–11469.
  * [Efron et al., 2004] Efron, B., Hastie, T., Johnstone, I., and Tibshirani, R. (2004).  Least angle regression.  The Annals of statistics, 32(2):407–499.
  * [Friedman, 2002] Friedman, J. H. (2002).  Stochastic gradient boosting.  Computational statistics & data analysis, 38(4):367–378.
  * [Fuller et al., 2020] Fuller, A., Fan, Z., Day, C., and Barlow, C. (2020).  Digital twin: Enabling technologies, challenges and open research.  IEEE Access, 8:108952–108971.
  * [Giurgiu and Schumann, 2019] Giurgiu, I. and Schumann, A. (2019).  Additive explanations for anomalies detected from multivariate temporal data.  In Proceedings of the 28th ACM International Conference on Information and Knowledge Management, pages 2245–2248. ACM.
  * [Gross, 2016] Gross, J. L. (2016).  Combinatorial methods with computer applications.  CRC Press.
  * [Hastie et al., 2009] Hastie, T., Tibshirani, R., and Friedman, J. (2009).  The Elements of Statistical Learning: Data Mining, Inference, and Prediction.  Springer, 2 edition.
  * [Hwang and Lee, 2021] Hwang, C. and Lee, T. (2021).  E-SFD: Explainable sensor fault detection in the ics anomaly detection system.  IEEE Access, 9:140470–140486.
  * [Idé et al., 2021] Idé, T., Dhurandhar, A., Navrátil, J., Singh, M., and Abe, N. (2021).  Anomaly attribution with likelihood compensation.  In Proceedings of the AAAI Conference on Artificial Intelligence (AAAI-21), volume 35, pages 4131–4138.
  * [Jiang et al., 2011] Jiang, R., Fei, H., and Huan, J. (2011).  Anomaly localization for network data streams with graph joint sparse PCA.  In Proceedings of the 17th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 886–894.
  * [Kumar et al., 2020] Kumar, I. E., Venkatasubramanian, S., Scheidegger, C., and Friedler, S. (2020).  Problems with Shapley-value-based explanations as feature importance measures.  In International Conference on Machine Learning, pages 5491–5500. PMLR.
  * [Lee and Xiang, 2000] Lee, W. and Xiang, D. (2000).  Information-theoretic measures for anomaly detection.  In Proceedings 2001 IEEE Symposium on Security and Privacy. S&P 2001, pages 130–143.
  * [Li et al., 2022] Li, Z., Zhu, Y., and van Leeuwen, M. (2022).  A survey on explainable anomaly detection.  arXiv preprint arXiv:2210.06959.
  * [Lundberg and Lee, 2017] Lundberg, S. M. and Lee, S.-I. (2017).  A unified approach to interpreting model predictions.  In Advances in Neural Information Processing Systems, pages 4765–4774.
  * [Madry et al., 2018] Madry, A., Makelov, A., Schmidt, L., Tsipras, D., and Vladu, A. (2018).  Towards deep learning models resistant to adversarial attacks.  In International Conference on Learning Representations.
  * [Mariadass et al., 2022] Mariadass, D. A., Moung, E. G., Sufian, M. M., and Farzamnia, A. (2022).  Extreme gradient boosting (xgboost) regressor and shapley additive explanation for crop yield prediction in agriculture.  In 2022 12th International Conference on Computer and Knowledge Engineering (ICCKE), pages 219–224. IEEE.
  * [Micenková et al., 2013] Micenková, B., Ng, R. T., Dang, X.-H., and Assent, I. (2013).  Explaining outliers by subspace separability.  In 2013 IEEE 13th international conference on data mining, pages 518–527.
  * [Molnar, 2022] Molnar, C. (2022).  Interpretable machine learning – A Guide for Making Black Box Models Explainable.  Independently published.
  * [Montavon et al., 2019] Montavon, G., Binder, A., Lapuschkin, S., Samek, W., and Müller, K.-R. (2019).  Layer-Wise Relevance Propagation: An Overview, pages 193–209.  Springer.
  * [Noto et al., 2010] Noto, K., Brodley, C., and Slonim, D. (2010).  Anomaly detection using an ensemble of feature models.  In 2010 IEEE International Conference on Data Mining, pages 953–958. IEEE.
  * [Pace and Barry, 1997] Pace, R. K. and Barry, R. (1997).  Sparse spatial autoregressions.  Statistics & Probability Letters, 33(3):291–297.
  * [Parikh et al., 2014] Parikh, N., Boyd, S., et al. (2014).  Proximal algorithms.  Foundations and Trends in Optimization, 1(3):127–239.
  * [Qin et al., 2019] Qin, C., Martens, J., Gowal, S., Krishnan, D., Dvijotham, K., Fawzi, A., De, S., Stanforth, R., and Kohli, P. (2019).  Adversarial robustness through local linearization.  Advances in Neural Information Processing Systems, 32.
  * [Ribeiro et al., 2016] Ribeiro, M. T., Singh, S., and Guestrin, C. (2016).  Why should I trust you?: Explaining the predictions of any classifier.  In Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining, pages 1135–1144. ACM.
  * [Roth, 1988] Roth, A. E. (1988).  The Shapley value: Essays in honor of Lloyd S. Shapley.  Cambridge University Press.
  * [Roy et al., 2017] Roy, V., Chakraborty, S., et al. (2017).  Selection of tuning parameters, solution paths and standard errors for Bayesian lassos.  Bayesian Analysis, 12(3):753–778.
  * [Salomon, 1998] Salomon, R. (1998).  Evolutionary algorithms and gradient search: similarities and differences.  IEEE Transactions on Evolutionary Computation, 2(2):45–55.
  * [Salomon and Arnold, 2009] Salomon, R. and Arnold, D. V. (2009).  The evolutionary-gradient-search procedure in theory and practice.  In Nature-Inspired Algorithms for Optimisation, pages 77–101. Springer.
  * [Samek et al., 2019] Samek, W., Montavon, G., Vedaldi, A., Hansen, L. K., and Müller, K.-R. (2019).  Explainable AI: interpreting, explaining and visualizing deep learning, volume 11700.  Springer Nature.
  * [scikit-learn 1.1.3, 2022] scikit-learn 1.1.3 (2022).  sklearn.datasets.load_boston.  https://scikit-learn.org/1.1/modules/generated/sklearn.datasets.load_boston.html.
  * [Selvaraju et al., 2017] Selvaraju, R. R., Cogswell, M., Das, A., Vedantam, R., Parikh, D., and Batra, D. (2017).  Grad-CAM: Visual explanations from deep networks via gradient-based localization.  In Proceedings of the IEEE International Conference on Computer Vision (ICCV).
  * [Simonyan et al., 2014] Simonyan, K., Vedaldi, A., and Zisserman, A. (2014).  Deep inside convolutional networks: Visualising image classification models and saliency maps.  In 2nd International Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April 14-16, 2014, Workshop Track Proceedings.
  * [Sipple, 2020] Sipple, J. (2020).  Interpretable, multidimensional, multimodal anomaly detection with negative sampling for detection of device failure.  In Proceedings of the 37th International Conference on Machine Learning (ICML 20).
  * [Sipple and Youssef, 2022] Sipple, J. and Youssef, A. (2022).  A general-purpose method for applying explainable ai for anomaly detection.  In International Symposium on Methodologies for Intelligent Systems, pages 162–174. Springer.
  * [Speith, 2022] Speith, T. (2022).  A review of taxonomies of explainable artificial intelligence (xai) methods.  In 2022 ACM Conference on Fairness, Accountability, and Transparency, pages 2239–2250.
  * [Staniford et al., 2002] Staniford, S., Hoagland, J. A., and McAlerney, J. M. (2002).  Practical automated detection of stealthy portscans.  Journal of Computer Security, 10(1-2):105–136.
  * [Štrumbelj and Kononenko, 2010] Štrumbelj, E. and Kononenko, I. (2010).  An efficient explanation of individual classifications using game theory.  Journal of Machine Learning Research, 11(Jan):1–18.
  * [Štrumbelj and Kononenko, 2014] Štrumbelj, E. and Kononenko, I. (2014).  Explaining prediction models and individual predictions with feature contributions.  Knowledge and information systems, 41(3):647–665.
  * [Sundararajan and Najmi, 2020] Sundararajan, M. and Najmi, A. (2020).  The many shapley values for model explanation.  In International conference on machine learning, pages 9269–9278. PMLR.
  * [Sundararajan et al., 2017] Sundararajan, M., Taly, A., and Yan, Q. (2017).  Axiomatic attribution for deep networks.  In Proceedings of the 34th International Conference on Machine Learning-Volume 70, pages 3319–3328.
  * [Tao et al., 2018] Tao, F., Zhang, H., Liu, A., and Nee, A. Y. (2018).  Digital twin in industry: State-of-the-art.  IEEE Transactions on Industrial Informatics, 15(4):2405–2415.
  * [Wachter et al., 2017] Wachter, S., Mittelstadt, B., and Russell, C. (2017).  Counterfactual explanations without opening the black box: Automated decisions and the GDPR.  Harvard Journal of Law & Technology, 31:841.
  * [Yamanishi et al., 2000] Yamanishi, K., Takeuchi, J., Williams, G., and Milne, P. (2000).  On-line unsupervised outlier detection using finite mixtures with discounting learning algorithms.  In Proc. the Sixth ACM SIGKDD Intl. Conf. on Knowledge Discovery and Data Mining, pages 320–324.
  * [Zhang et al., 2019] Zhang, X., Marwah, M., Lee, I.-t., Arlitt, M., Goldwasser, D., et al. (2019).  ACE—An anomaly contribution explainer for cyber-security applications.  In Proceedings of the 2019 IEEE International Conference on Big Data (Big Data), pages 1991–2000. IEEE.
  * [Zhou et al., 2022] Zhou, Y., Booth, S., Ribeiro, M. T., and Shah, J. (2022).  Do feature attribution methods correctly attribute features?  In Proceedings of the AAAI Conference on Artificial Intelligence, pages 9623–9633.

##  Appendix A Efficiency of Shapley value

In this section, we prove the “efficiency” of SV in Eq. (16):

$\displaystyle\sum_{i=1}^{M}\mathrm{SV}_{i}(\bm{x}^{t})=f(\bm{x}^{t})-\langle
f\rangle.$

###### Proof.

In

$\displaystyle\sum_{i=1}^{M}\mathrm{SV}_{i}(\bm{x}^{t})=\frac{1}{M}\sum_{i=1}^{M}\sum_{k=0}^{M-1}\binom{M-1}{k}^{-1}\sum_{\mathcal{S}_{i}:|\mathcal{S}_{i}|=k}\left[\langle
f\mid x_{i}^{t},\bm{x}_{\mathcal{S}_{i}}^{t}\rangle-\langle
f\mid\bm{x}_{\mathcal{S}_{i}}^{t}\rangle\right].$ (A.1)

let us define

$\displaystyle I^{M}_{k}$
$\displaystyle\triangleq\frac{1}{M}\sum_{i=1}^{M}\binom{M-1}{k}^{-1}\sum_{\mathcal{S}_{i}:|\mathcal{S}_{i}|=k}\langle
f\mid x_{i}^{t},\bm{x}_{\mathcal{S}_{i}}^{t}\rangle,$ (A.2) $\displaystyle
J^{M}_{k}$
$\displaystyle\triangleq\frac{1}{M}\sum_{i=1}^{M}\binom{M-1}{k+1}^{-1}\sum_{{\mathcal{S}}_{i}:|{\mathcal{S}}_{i}|=k+1}\langle
f\mid\bm{x}^{t}_{{\mathcal{S}}_{i}}\rangle$ (A.3)

for $k=0,\ldots,M-2$. Equation (16) holds if $I^{M}_{k}-J^{M}_{k}=0$. In
$I^{M}_{k}$, for a given $k$, the number of distinct sets
$\\{i\\}\cup\mathcal{S}_{i}$ is $\binom{M}{k+1}$, and the summations in
$I^{M}_{k}$ runs over $M\times\binom{M-1}{k}$ terms in total. Hence, each
unique set appears

$\displaystyle M\times\binom{M-1}{k}\times\frac{1}{\binom{M}{k+1}}=k+1$ (A.4)

times. Following the same argument, we see that each unique set in $J^{M}_{k}$
appears $M-k-1$.

Let $\mathcal{S}(k+1)$ be a set of $k+1$ variable indices to represent either
$\\{i\\}\cup\mathcal{S}_{i}$ with $|\mathcal{S}_{i}|=k$ or $\mathcal{S}_{i}$
with $|\mathcal{S}_{i}|=k+1$. For each member in $\mathcal{S}(k+1)$,
$I^{M}_{k}$ gives a prefactor

$\displaystyle(k+1)\times\frac{1}{M}\binom{M-1}{k}^{-1}=\binom{M}{k+1}^{-1},$
(A.5)

and $J^{M}_{k}$ gives a prefactor

$\displaystyle(M-k-1)\times\frac{1}{M}\binom{M-1}{k+1}^{-1}=\binom{M}{k+1}^{-1},$
(A.6)

which are the same. Hence, we conclude $I^{M}_{k}-J^{M}_{k}=0$. ∎