# Explainable Learning with Gaussian Processes

Kurt Butler<EMAIL_ADDRESS>
Guanchao Feng<EMAIL_ADDRESS>
Petar M. Djurić<EMAIL_ADDRESS>
Department of Electrical and Computer Engineering
Stony Brook University
Stony Brook, NY 11794, USA

###### Abstract

The field of explainable artificial intelligence (XAI) attempts to develop
methods that provide insight into how complicated machine learning methods
make predictions. Many methods of explanation have focused on the concept of
feature attribution, a decomposition of the model’s prediction into individual
contributions corresponding to each input feature. In this work, we explore
the problem of feature attribution in the context of Gaussian process
regression (GPR). We take a principled approach to defining attributions under
model uncertainty, extending the existing literature. We show that although
GPR is a highly flexible and non-parametric approach, we can derive
interpretable, closed-form expressions for the feature attributions. When
using integrated gradients as an attribution method, we show that the
attributions of a GPR model also follow a Gaussian process distribution, which
quantifies the uncertainty in attribution arising from uncertainty in the
model. We demonstrate, both through theory and experimentation, the
versatility and robustness of this approach. We also show that, when
applicable, the exact expressions for GPR attributions are both more accurate
and less computationally expensive than the approximations currently used in
practice. The source code for this project is freely available under MIT
license at https://github.com/KurtButler/2024_attributions_paper.

Keywords: explainable ai, feature attribution, gaussian processes, integrated
gradients, sparse approximation

###### Contents

  1. 1 Introduction
    1. 1.1 Related Work
    2. 1.2 Contributions
  2. 2 The Feature Attribution Problem
    1. 2.1 Motivation: Bayesian Linear Regression
    2. 2.2 Integrated Gradients
      1. 2.2.1 Geometric Interpretation
      2. 2.2.2 Axiomatic Interpretation
    3. 2.3 Linear Operators
  3. 3 Gaussian Process Regression
    1. 3.1 Covariance Functions
    2. 3.2 Derivatives and Integrals of GPs
    3. 3.3 Spectral Approximations to GPR
  4. 4 Attributions in GPR Models
    1. 4.1 SE and ARD-SE Kernels
    2. 4.2 Random Feature GPs
    3. 4.3 Comments on Other GP Models
      1. 4.3.1 Additive GPs
      2. 4.3.2 Nonstationary GPs
  5. 5 Experiments
    1. 5.1 Simulated Data: Heteroscedasticity of the Attribution GPs
    2. 5.2 Breast Cancer Prognosis: A General Example
    3. 5.3 Taipei Housing Data: Demonstrating Affine Scale Invariance
    4. 5.4 Wine: Approximating the Integral in IG Attributions
    5. 5.5 Wine: Comparing the Attributions of RFGPs and Exact GPR
  6. 6 Conclusion
  7. A Proofs
    1. A.1 Operations that Preserve Gaussianity
    2. A.2 Derivatives of GPs
    3. A.3 Integrals of GPs
    4. A.4 Attributions of a GP
  8. B Derivations for Specific Models
  9. C Additional Experiments
    1. C.1 Wine: Quadrature Rules

##  1 Introduction

With the widespread adoption of machine learning (ML) methods across academia
and industry, there has been an ever increasing desire to explain how complex
artificial intelligence (AI) models make predictions Samek et al. (2021);
Angelov et al. (2021). The ability to explain AI models, in this sense, is
important to promote user trust in AI systems, and to mitigate the misuse of
dubious predictions Hess and Ditzler (2023). However, the process of
interpreting a model to gain physical insight is often difficult unless the
model was intentionally designed to be interpretable Peng et al. (2022).
Interpretable models are often transparent in how predictions are made, as is
the case for linear models Arrieta et al. (2020) and shallow decision trees
Silva et al. (2020), or they contain components in their architecture that are
designed for interpretability, such as physics-informed neural networks Huang
and Wang (2022). The sentiment that many ML methods are ‘black boxes’ reflects
the inherent difficulty of interpreting complicated, data-driven models like
deep neural networks and kernel machines Zhang et al. (2018); Rudin (2019).
Despite this, the ML paradigm has remained immensely useful, and the
exploration of post-hoc explanations for black box models has been a popular
subject in recent years, often referred to as explainable artificial
intelligence (XAI).

A wide variety of distinct notions of explainability exist in the literature
Jin et al. (2022). Many XAI approaches can be described as being either model-
specific or model agnostic. Model-specific approaches attempt to formulate
explanations by exploiting the known architecture of a ML model. In the
context of neural networks, this includes methods that interpret activation
weights Springenberg et al. (2014); Bach et al. (2015), graphs Ying et al.
(2019), or attention mechanisms Voita et al. (2019) to explain how predictions
were made. Model-specific approaches also include the use of models that are
inherently interpretable, including decision trees Silva et al. (2020) and
kernel machines with Automatic Relevance Determination (ARD) Neal (1996).

In contrast, model-agnostic approaches seek notions of explainability that can
be applied to a wide range of models, with only generic assumptions about the
type of data or properties of the model. Explanations of this type might
measure the sensitivity of the prediction to perturbations at the input
Baehrens et al. (2010); Simonyan et al. (2014) or measuring changes in the
model prediction after masking or perturbing input features Robnik-Šikonja and
Kononenko (2008); Štrumbelj et al. (2009). In prediction tasks where images
are inputs, a popular tool of this type is that of saliency maps, which are
images that visually highlight the relevant features in a prediction Adebayo
et al. (2018). Another class of methods attempt to explain a prediction by
showing relevant similar examples that were observed in the training set,
sometimes called example-based explanations van der Waa et al. (2021). Reverse
image searching approaches attempt to find similar training examples to a
given prediction and have been of particular interest in the field of image
diffusion models Carlini et al. (2023). Finally, many XAI methods can be
understood through the invocation of an approximate, interpretable model, as
seen in Local Interpretable Model-agnostic Explanations (LIME) Ribeiro et al.
(2016) and SHapley Additive exPlanations (SHAP) Lundberg and Lee (2017).

Many XAI methods attempt to explain a prediction by associating a numerical
score with each input feature. A natural approach to assigning this score is
to attempt to measure the contribution of each feature to the output, that is,
how much of the resulting prediction do we _attribute_ to each input feature.
This approach to explainability, called _feature attribution_ , was
popularized by Štrumbelj and Kononenko (2014), who introduced Shapley values
to XAI. Shapley values were proposed by Shapley (1953) as a method for
determining what fraction of the total profit each player is owed in a
coalitional game. In the context of XAI, this was reinterpreted to quantify
how much each feature contributed to a numerical prediction. Since then,
several other feature attribution approaches have been proposed, which we
review in the next section.

In this paper, we will focus on Gaussian process regression (GPR) as a general
non-parametric machine learning paradigm and assess the utility of post-hoc
explanations for GPR models. This setting is intriguing because (a)
interesting theoretical results are obtained from computing the attributions
of a GPR model and (b) the resulting form of the GPR result is heavily
dependent on the observed data set. Since the GPR paradigm is data-driven, any
post-hoc explanation in this setting must be compatible with the lack of
explicit structural knowledge in the GPR approach. Acknowledging critiques of
post-hoc explanations and the XAI paradigm Rudin (2019); Ghassemi et al.
(2021); de Bruijn et al. (2022), we will argue that the attributions can
remain useful to characterize aspects of the learned model, as long as the
practitioners’ expectations are tempered to be in line with what information
can actually be recovered from the attributions.

###  1.1 Related Work

Since our primary notion of explainability in this work is feature
attribution, we briefly highlight existing related work in this area. The most
well-known approach to feature attribution is that of Shapley values,
originally developed from game theoretic principles for econometric
applications by Shapley (1953). Shapley values benefit from an axiomatic
definition that uniquely defines them, but this original formulation is only
valid for discrete features Aumann and Shapley (2015). In continuous domains,
the original axioms are insufficient to uniquely define the method
Sundararajan and Najmi (2020). Several extensions of the Shapley value to
continuous features exist, including Aumann-Shapley Aumann and Shapley (2015)
and Shapley-Shubik Friedman and Moulin (1999). Other approaches to bring
Shapley values to continuous domains have attempted the inclusion/deletion of
features Grömping (2007), or through the invocation of a proxy model that
locally explains the original model, including LIME Ribeiro et al. (2016) and
SHAP Lundberg and Lee (2017). The Integrated Gradients (IG) framework of
Sundararajan et al. (2017) generalizes the framework of Aumann-Shapley.
Sturmfels et al. (2020) and Lundberg and Lee (2017) considered variations on
the IG attribution theory where baselines are given by a distribution of
possible baselines.

While the majority of research in feature attribution has concentrated on
neural models, there are authors who explore Gaussian process (GP) models.
Chau et al. (2023) proposed a GP-based version of SHAP analysis to provide
uncertainty quantification of Shapley values. Seitz (2022) considered the
application of IG attributions to GP models for regression and classification.
Previous research has also focused on designing GPR models to be
interpretable, typically through kernel design. This can be seen in the work
on ARD Neal (1996) and additive kernels Duvenaud et al. (2011). Yoshikawa and
Iwata (2021) considered a GP-based hierarchical model in which predictions are
made using conditionally linear models of interpretable features, with the
model weights derived from nonlinear functions of other non-interpretable
features.

###  1.2 Contributions

This paper examines several aspects of attribution theory in the context of GP
models. We show that the attributions represent a decomposition of the
predictive model into a sum of component functions, and this decomposition
preserves the differentiability of the original model. When a GP probability
distribution is used to model the function space, we show that the attribution
functions are also distributed according to GPs.

We derive exact expressions for the attributions of various common GPR models.
The resulting expressions are straightforward to implement and check in code.
We also comment on extensions of these results to more sophisticated kernels
and GPR models.

Using the attribution GP, we gain the ability to measure the confidence of
each attribution. We demonstrate through experimental data cases where these
considerations are relevant to ML practitioners in medicine and industry.

We also compare the performance of exact attributions to various
approximations that have appeared in the literature. These approximations are
often asymptotic, and to the knowledge of the authors no experimental work has
been done to ascertain how the convergence of these approximate solutions
relates to the measured attributions.

##  2 The Feature Attribution Problem

Consider the setting of multivariable regression. Given observations $y_{n}$,
for $n=1,...,N$, of an unknown function with corresponding input feature
vectors $\mathbf{x}_{n}$ in $\mathbb{R}^{D}$, the goal is to estimate a
function $F:\mathbb{R}^{D}\to\mathbb{R}$ that approximates the data, that is,
$y_{n}\approx F(\mathbf{x}_{n})$. Once we have learned a particular model $F$
for the data, we can make new predictions, and we may be interested to know
how each of the $D$ input features contribute to a given prediction. For
example, if the function $F$ is supposed to represent a clinical prediction,
such as the time needed to recovery from an injury, then we may be interested
to know what features in the prediction model are the most important for
determining the patient’s outcome. The goal is that by understanding the
relative importance of features used for prediction, we may be able to explain
how the model leverages this information to make predictions.

As stated, the explainability problem is still vague, and there have been
various interpretations of this idea that make the problem more concrete. In
the baseline attribution problem Lundstrom et al. (2022), one ‘explains’ a
prediction $F(\mathbf{x})$ at a location $\mathbf{x}$ by measuring how it
differs from the prediction at a baseline location $\tilde{\mathbf{x}}$ (a
reference point) and allocating a fraction of the observed difference to each
input feature $x_{i}$. Generally, the attributions to each feature should add
up to the difference between the prediction $F(\mathbf{x})$ and the reference
prediction $F(\tilde{\mathbf{x}})$, which is called the _completeness_
property. Many feature attribution approaches use a similar approach, most
notably approaches based on Shapley values Shapley (1953). We will focus on
the baseline attribution problem in this paper.

Notation: Let $F:\mathbb{R}^{D}\to\mathbb{R}$ be a function representing a
regression model, $\mathbf{x}\in\mathbb{R}^{D}$ be a feature vector from which
we are interested in making a prediction, and
$\tilde{\mathbf{x}}\in\mathbb{R}^{D}$ be a reference point. The attribution of
the prediction to the feature $x_{i}$ is denoted by $attr_{i}(\mathbf{x}|F)$.
When the function $F$ is understood, we may suppress it from the notation and
write simply $attr_{i}(\mathbf{x})$. Using this notation, the completeness
property can be written as

$F(\mathbf{x})-F(\tilde{\mathbf{x}})=\sum_{i=1}^{D}attr_{i}(\mathbf{x}|F).$
(1)

We sometimes have the need to compute the attributions of a vector-valued
function $\mathbf{F}(\mathbf{x})=[F_{1}(\mathbf{x})\cdots
F_{Q}(\mathbf{x})]^{\top}$. In this case, we define the attribution vectors
$\mathbf{attr}_{i}(\mathbf{x}|\mathbf{F})$ by computing the attribution for
each component of $\mathbf{F}$:

$\mathbf{attr}_{i}(\mathbf{x}|\mathbf{F})\stackrel{{\scriptstyle\Delta}}{{=}}\begin{bmatrix}attr_{i}(\mathbf{x}|F_{1})\\\
\vdots\\\ attr_{i}(\mathbf{x}|F_{Q})\end{bmatrix},\qquad i=1,...,D.$ (2)

Finally, we use the notation $C^{k}(\mathcal{X})$ to denote the set of all
$k$-times continuously differentiable functions $\mathcal{X}\to\mathbb{R}$,
where $k=1,2,...$ or $k=\infty$. We may write $C^{k}$ when the functions’
domain is clear from context.

###  2.1 Motivation: Bayesian Linear Regression

To motivate our work in this paper, it will be instructive to briefly consider
the problem of Bayesian linear regression (Murphy, 2022, sec.11.7). Our
reasons are twofold: Firstly, linear models provide a canonical notion of
attribution which is generalized by the IG method, and second, it explicates
how uncertainty in the model gives rise to uncertainty in the attributions.

Consider the typical Bayesian linear regression model, in which predictions
$y$ at a point $\mathbf{x}$ are given a model of the form

$y|\mathbf{x},\mathbf{w}\sim\mathcal{N}(\mathbf{w}^{\top}\mathbf{x},\sigma^{2}),$

where the noise (or residual) is tacitly assumed to be zero-mean, and
$\sigma^{2}$ is its variance. Predictions from this model are made using the
expected value, so our predictive function is simply
$F(\mathbf{x})=\mathbb{E}(y|\mathbf{x},\mathbf{w})=\mathbf{w}^{\top}\mathbf{x}$.
After choosing an appropriate baseline point $\tilde{\mathbf{x}}$, the
difference of a given prediction to the baseline can always be expressed as a
sum over the feature index $i$,

$F(\mathbf{x})-F(\tilde{\mathbf{x}})=\mathbf{w}^{\top}(\mathbf{x}-\tilde{\mathbf{x}})=\sum_{i=1}^{D}w_{i}(x_{i}-\tilde{x}_{i}).$

In this way, the quantity $attr_{i}(\mathbf{x})=w_{i}(x_{i}-\tilde{x}_{i})$ is
the canonical choice for the attribution to the $x_{i}$ feature. Attributions
of this form satisfy the completeness property, and are straightforward to
interpret. The challenge of attribution methods for nonlinear models is to
generalize this notion in a manner that is consistent with the linear case.

Now let us consider the situation in which the model parameters $\mathbf{w}$
are random variables, and not fixed unknown parameters. In this way,
attributions in a Bayesian setting should account for the uncertainty in the
model. The proper way to do this is to take the attributions above, and to use
the posterior distribution over $\mathbf{w}$ to quantify our uncertainty in
$attr_{i}(\mathbf{x})$. We now do this for a simple example.

Suppose that we have a data set
$\mathcal{D}=\\{(\mathbf{x}_{n},y_{n})|n=1,...,N\\}$, which we compactly
represent in matrix form:

$\mathbf{X}\stackrel{{\scriptstyle\Delta}}{{=}}\begin{bmatrix}\mathbf{x}_{1}^{\top}\\\
\vdots\\\ \mathbf{x}_{N}^{\top}\end{bmatrix}\in\mathbb{R}^{N\times
D},\qquad\mathbf{y}\stackrel{{\scriptstyle\Delta}}{{=}}\begin{bmatrix}y_{1}\\\
\vdots\\\ y_{N}\end{bmatrix}\in\mathbb{R}^{N}.$

If we assume a Gaussian prior for the model weights,

$\displaystyle\mathbf{w}\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma}),$
(3)

for some fixed hyperparameters $\boldsymbol{\mu},\boldsymbol{\Sigma}$, then
the posterior distribution for $\mathbf{w}$ is also Gaussian and given by,

$\displaystyle\mathbf{w}|\mathbf{X},\mathbf{y}$
$\displaystyle\sim\mathcal{N}(\boldsymbol{\mu}^{\prime},\boldsymbol{\Sigma}^{\prime}),$
(4) $\displaystyle\boldsymbol{\Sigma}^{\prime}$
$\displaystyle=\left(\boldsymbol{\Sigma}^{-1}+\frac{1}{\sigma^{2}}\mathbf{X}^{\top}\mathbf{X}\right)^{-1},$
(5) $\displaystyle\boldsymbol{\mu}^{\prime}$
$\displaystyle=\boldsymbol{\Sigma}^{\prime}\left(\boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}+\frac{1}{\sigma^{2}}\mathbf{X}^{\top}\mathbf{y}\right).$
(6)

With the posterior distribution in mind, and observing that
$attr_{i}(\mathbf{x})=w_{i}(x_{i}-\tilde{x}_{i})$ is just a scaled copy of
$w_{i}$, it becomes easy to see thet the marginal distribution for each
attribution is given by

$\displaystyle attr_{i}(\mathbf{x})|\mathbf{X},\mathbf{y}$
$\displaystyle\sim\mathcal{N}(\mu_{i}^{\prime}(x_{i}-\tilde{x}_{i}),\Sigma_{ii}^{\prime}(x_{i}-\tilde{x}_{i})^{2}).$
(7)

From this expression, we can immediately observe two important facts: the
variance or uncertainty in the attribution comes from uncertainty in the
model, and this uncertainty scales based on the distance from the baseline in
the $x_{i}$ coordinate. In Section 4, we obtain similar results for the GPR
model, which is a Bayesian non-parametric regression method. We note that the
GPR case parallels the basic results shown here.

###  2.2 Integrated Gradients

Consider a function $F:\mathbb{R}^{D}\to\mathbb{R}$ which we use as a
predictive model and for which we want to compute the attributions. The
Integrated Gradient (IG) method defines the _attribution_ of the prediction to
the $i$th input feature, given a specific input $\mathbf{x}$, to be

$attr_{i}(\mathbf{x}|F)\stackrel{{\scriptstyle\Delta}}{{=}}(x_{i}-\tilde{x}_{i})\int_{0}^{1}\frac{\partial
F(\tilde{\mathbf{x}}+t(\mathbf{x}-\tilde{\mathbf{x}}))}{\partial x_{i}}dt.$
(8)

The expression in (8) can be motivated from two perspectives, one using
geometric ideas and the other using axioms. Both perspectives are useful for
gaining insight into how we interpret IG attributions.

####  2.2.1 Geometric Interpretation

The fundamental theorem of line integrals states that the line integral of a
gradient equals the difference of the function evaluated at the endpoints
(Lee, 2012, p. 291). For a continuously differentiable function
$F:\mathbb{R}^{D}\to\mathbb{R}$,

$\int_{\Gamma}\nabla F\cdot
d\mathbf{x}=F(\boldsymbol{\gamma}(1))-F(\boldsymbol{\gamma}(0)),$

where

$\nabla=\begin{bmatrix}\frac{\partial}{\partial
x_{1}}&\cdots&\frac{\partial}{\partial x_{D}}\end{bmatrix}^{\top}$

is the gradient operator, and $\Gamma$ is a path in $\mathbb{R}^{D}$
parameterized by $\boldsymbol{\gamma}(t)$ for $0\leq t\leq 1$. The left-hand
side of the equation can be expanded to be expressed in terms of the partial
derivatives, which reveals that the integral may be expressed as a sum of path
integrals. In particular, if we define the path
$\boldsymbol{\gamma}(t)=\tilde{\mathbf{x}}+t(\mathbf{x}-\tilde{\mathbf{x}})$,
then we may write

$\displaystyle\int_{\Gamma}\nabla F\cdot
d\mathbf{x}=\int_{0}^{1}\sum_{i=1}^{D}\frac{\partial
F(\boldsymbol{\gamma}(t))}{\partial\gamma_{i}}\frac{d\gamma_{i}}{dt}dt$
$\displaystyle=\int_{0}^{1}\sum_{i=1}^{D}\frac{\partial
F(\boldsymbol{\gamma}(t))}{\partial\gamma_{i}}(x_{i}-\tilde{x}_{i})dt$ (9)
$\displaystyle=\sum_{i=1}^{D}(x_{i}-\tilde{x}_{i})\int_{0}^{1}\frac{\partial
F(\boldsymbol{\gamma}(t))}{\partial x_{i}}dt.$ (10)

The terms in the sum in (10) are precisely the IG attributions in (8). This
perspective allows us to see that the completeness property, (1), is
automatically satisfied for any attributions that are defined by a path,
although we must pick a specific path to yield the IG attributions. The
integral $\int_{0}^{1}\partial F/\partial x_{i}dt$ measures the average value
of $\partial F/\partial x_{i}$ between the baseline $\tilde{\mathbf{x}}$ and
the prediction at $\mathbf{x}$, and acts as a nonlinear analogue to the
weights of a linear model.

If the function $F$ has constant partial derivatives, $\partial F/\partial
x_{i}=w_{i}$, for each $i$, then $F$ is a linear function. In this case,
$attr_{i}(\mathbf{x})=(x_{i}-\tilde{x}_{i})\int_{0}^{1}\partial F/\partial
x_{i}dt=w_{i}(x_{i}-\tilde{x}_{i})$, which agrees with what we defined for the
linear model. For nonlinear functions $F$, the IG attributions generalize the
linear case by considering the average value of the gradient between the
baseline and the new prediction.

Since the fundamental theorem of line integrals is not exclusive to
differentiable models, but also applies to any continuous, piecewise-
differentiable model, one can apply IG attributions to a much wider family of
models and maintain the same geometric motivation. This includes ReLU networks
Glorot et al. (2011); Krizhevsky et al. (2012), splines Wahba (1990), and GPs
with deep kernel learning Wilson et al. (2016).

####  2.2.2 Axiomatic Interpretation

In the geometric approach to defining the IG attributions, we arbitrarily
chose $\boldsymbol{\gamma}(t)$ to be the straight-line path
$\tilde{\mathbf{x}}+t(\mathbf{x}-\tilde{\mathbf{x}})$. In principle, any path
$\boldsymbol{\gamma}(t)$ that connects $\mathbf{x}$ and $\tilde{\mathbf{x}}$
could define a set of attributions, but a different choice of path can yield
different values for the attributions. To this end, Sundararajan et al. (2017)
and Lundstrom et al. (2022) investigated the uniqueness of IG attributions.
Building upon earlier work by Friedman (2004), they proved that IG
attributions, with
$\boldsymbol{\gamma}(t)=\tilde{\mathbf{x}}+t(\mathbf{x}-\tilde{\mathbf{x}})$,
are uniquely determined by the seven axioms listed below.

We note that several of these axioms are related to those that define Shapley
values Mitchell et al. (2022). However, Shapley values were originally defined
for discrete covariates, and additional axioms are required to uniquely extend
Shapley values to the continuous setting. If a different set of axioms are
selected, then the resulting attribution operator may be different and
represent a different attribution theory. Sundararajan and Najmi (2020)
further discuss how different axioms induce different theories of feature
attribution.

  1. 1.

Sensitivity(a). If $\mathbf{x}$ and $\tilde{\mathbf{x}}$ differ in one
coordinate ($x_{i}\neq\tilde{x}_{i}$ and $x_{j}=\tilde{x}_{j}$ for $j\neq i$),
and also $F(\mathbf{x})\neq F(\tilde{\mathbf{x}})$, then
$attr_{i}(\mathbf{x}|F)\neq 0$.

  2. 2.

Sensitivity(b), also known as Dummy Friedman (2004). If $F$ is a constant
function with respect to $x_{i}$, then $attr_{i}(\mathbf{x}|F)=0$. In
particular, this means that if $\partial F/\partial x_{i}=0$ everywhere, then
the attribution is zero.

  3. 3.

Implementation invariance. The attributions should depend only on the
mathematical properties of the predictive model, and not on a particular
architecture or implementation. If $F,G$ are two models such that for all
$\mathbf{x}$, $F(\mathbf{x})=G(\mathbf{x})$, then
$attr_{i}(\mathbf{x}|F)=attr_{i}(\mathbf{x}|G).$

  4. 4.

Linearity. If $F$ can be expressed as a weighted sum,
$F(\mathbf{x})=\sum_{m=1}^{M}w_{m}F_{m}(\mathbf{x})$, then the attributions of
$F$ are also a weighted sum:

$attr_{i}(\mathbf{x}|F)=attr_{i}\left(\mathbf{x}\bigg{|}\sum_{m=1}^{M}w_{m}F_{m}\right)=\sum_{m=1}^{M}w_{m}attr_{i}(\mathbf{x}|F_{m}).$

  5. 5.

Completeness. The attributions sum to the difference between the prediction
and baseline,

$\sum_{i=1}^{D}attr_{i}(\mathbf{x}|F)=F(\mathbf{x})-F(\tilde{\mathbf{x}}).$

  6. 6.

Symmetry-preserving. For a fixed pair $(i,j)$, define
$\sigma_{ij}(\mathbf{x})$ to be $\mathbf{x}$ but with the entries $x_{i}$ and
$x_{j}$ swapped. We say that $x_{i}$ and $x_{j}$ are _symmetric w.r.t. a
function $F$_ if $F(\sigma_{ij}(\mathbf{x}))=F(\mathbf{x})$. An attribution
method is _symmetry-preserving_ if

$attr_{i}(\mathbf{x}|F)=attr_{j}(\mathbf{x}|F),$

whenever $x_{i}=x_{j}$, $\tilde{x}_{i}=\tilde{x}_{j}$, and $x_{i}$ and $x_{j}$
are symmetric w.r.t. $F$.

  7. 7.

Non-decreasing positivity. A path $\boldsymbol{\gamma}(t)$ in $\mathbb{R}^{D}$
is called _monotone_ for each $i=1,...,D$, we have that either
$d\gamma_{i}(t)/dt>0$ or $d\gamma_{i}(t)/dt<0$ for all $t$. _Non-decreasing
positivity_ asserts that if $F$ is non-decreasing along _every_ monontone path
from $\tilde{\mathbf{x}}$ to $\mathbf{x}$, then $attr_{i}(\mathbf{x}|F)\geq 0$
for all $i=1,...,D$.

Given these seven axioms, uniqueness can be claimed when $F\in
C^{1}(\mathbb{R}^{D})$, that is, when $F$ is continuously differentiable
Lundstrom et al. (2022). Whether or not a machine learning model satisfies
requirements like $F\in C^{1}$ depends on the model architecture. The
differentiability requirement is often automatically satisfied when using GPR,
which we discuss in Section 3.2.

###  2.3 Linear Operators

The basic motivation for feature attribution so far has focused on the case in
which we have a fixed input $\mathbf{x}$ from which we want to make a
prediction. However, as we vary the point $\mathbf{x}$ continuously throughout
the input space, the attributions $attr_{i}(\mathbf{x})$ also change
continuously. For this reason, it may be fruitful to not only consider
$attr_{i}(\mathbf{x}|F)$ as a fixed quantity, but rather as a function of the
input $\mathbf{x}$.

We define the attribution operator $attr_{i}$ to be the linear operator that
takes a base function $F$ and produces the attribution function
$attr_{i}(\mathbf{x}|F)$. In this way, completeness can be viewed as a
decomposition of the predictive function $F$ into components that correspond
to the attributions to each feature. The attribution operators also preserve
smoothness, since $attr_{i}(\mathbf{x}|F)\in C^{\infty}(\mathbb{R}^{D})$
whenever $F\in C^{\infty}(\mathbb{R}^{D})$ holds.111To see this, notice that
IG attributions in (8) can be defined in 4 steps: taking a derivative,
precomposing with a path
$\boldsymbol{\gamma}(t)=\tilde{\mathbf{x}}+t(\mathbf{x}-\tilde{\mathbf{x}})$,
integrating against $t$, and then scaling by $x_{i}-\tilde{x}_{i}$. All of
these operations preserve the $C^{\infty}$ property. Changing the baseline
$\tilde{\mathbf{x}}$ will yield a different set of attribution operators, but
they share the same properties regardless of the particular choice of
$\tilde{\mathbf{x}}$. For this reason, the baseline $\tilde{\mathbf{x}}$ may
be considered as a hyperparameter.

In Section 4, linearity of the attribution operator is the main property
needed to derive results for GPR models. While we only analyze IG attributions
in detail, the operator theoretic perspective clarifies how our results can be
extended to other attribution theories.

##  3 Gaussian Process Regression

In this section, we briefly review the basics of GPR and establish our
notation. The basic GPR model assumes that the data to be predicted
$y_{n},n=1,...,N$, can be represented by an unknown function $F$ of known
input vectors $\mathbf{x}_{n}$ for each point, with a constant additive noise:

$y_{n}=F(\mathbf{x}_{n})+\varepsilon_{n},\qquad\varepsilon_{n}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathcal{N}(0,\sigma^{2}),\qquad n=1,...,N.$ (11)

The goal of GPR is to infer the function $F$ in a non-parametric and Bayesian
manner, given the set of training data
$\\{(\mathbf{x}_{n},y_{n});n=1,...,N\\}$. To learn the function, one must
first put a prior distribution over $F$. Generally, the prior is used to
encode qualitative properties of the function such as continuity,
differentiability or periodicity Rasmussen and Williams (2006). In GPR, the
prior for $F$ is taken to be a _Gaussian process_ with mean function $m$ and
covariance function $k$, denoted as

$F\sim\mathcal{GP}(m,k).$

The function $F({\mathbf{x}})$ is a GP if for any finite set of inputs
$\mathbf{x}_{n}$, the function’s values $F(\mathbf{x}_{n})$ have a jointly
Gaussian distribution specified by $m$ and $k$, that is,

$\begin{bmatrix}F(\mathbf{x}_{1})\\\ \vdots\\\
F(\mathbf{x}_{N})\end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}m(\mathbf{x}_{1})\\\
\vdots\\\
m(\mathbf{x}_{N})\end{bmatrix},\begin{bmatrix}k(\mathbf{x}_{1},\mathbf{x}_{1})&\cdots&k(\mathbf{x}_{1},\mathbf{x}_{N})\\\
\vdots&&\vdots\\\
k(\mathbf{x}_{N},\mathbf{x}_{1})&\cdots&k(\mathbf{x}_{N},\mathbf{x}_{N})\end{bmatrix}\right).$
(12)

It is conventional to introduce some notation here. We write

$\mathbf{m}\stackrel{{\scriptstyle\Delta}}{{=}}\begin{bmatrix}m(\mathbf{x}_{1})\\\
\vdots\\\
m(\mathbf{x}_{N})\end{bmatrix},\qquad\mathbf{K}\stackrel{{\scriptstyle\Delta}}{{=}}\begin{bmatrix}k(\mathbf{x}_{1},\mathbf{x}_{1})&\cdots&k(\mathbf{x}_{1},\mathbf{x}_{N})\\\
\vdots&&\vdots\\\
k(\mathbf{x}_{N},\mathbf{x}_{1})&\cdots&k(\mathbf{x}_{N},\mathbf{x}_{N})\end{bmatrix},\qquad\mathbf{k}(\mathbf{x})\stackrel{{\scriptstyle\Delta}}{{=}}\begin{bmatrix}k(\mathbf{x}_{1},\mathbf{x})\\\
\vdots\\\ k(\mathbf{x}_{N},\mathbf{x})\end{bmatrix},$

to make the equations more concise.

To learn functions using GPR, we consider what happens when we augment (12)
with a new data point $\mathbf{x}_{*}$, separate from the data used during
training, and we attempt to predict $F(\mathbf{x}_{*})$, the value of the
function at this location. From (11) and (12), we may write that

$\begin{bmatrix}\mathbf{y}\\\
F(\mathbf{x}_{*})\end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}\mathbf{m}\\\
m(\mathbf{x}_{*})\end{bmatrix},\begin{bmatrix}\mathbf{K}+\sigma^{2}\mathbf{I}_{N}&\mathbf{k}(\mathbf{x}_{*})\\\
\mathbf{k}(\mathbf{x}_{*})^{\top}&k(\mathbf{x}_{*},\mathbf{x}_{*})\end{bmatrix}\right),$
(13)

where $\mathbf{y}=[y_{1}\cdots y_{N}]^{\top}$ and $\mathbf{I}_{N}$ is the
$N\times N$ identity matrix. Since the values of $\mathbf{y}$ were already
observed, we can condition on this data to obtain a posterior probability
distribution for $F(\mathbf{x}_{*})$, which is again Gaussian,

$\displaystyle F(\mathbf{x}_{*})|\mathbf{y}$
$\displaystyle\sim\mathcal{N}\left(\mu_{*},\sigma^{2}_{*}\right),$ (14)
$\displaystyle\mu_{*}$
$\displaystyle=m(\mathbf{x}_{*})+\mathbf{k}(\mathbf{x}_{*})^{\top}(\mathbf{K}+\sigma^{2}\mathbf{I}_{N})^{-1}(\mathbf{y}-\mathbf{m}),$
(15) $\displaystyle\sigma^{2}_{*}$
$\displaystyle=k(\mathbf{x}_{*},\mathbf{x}_{*})-\mathbf{k}(\mathbf{x}_{*})^{\top}(\mathbf{K}+\sigma^{2}\mathbf{I})^{-1}\mathbf{k}(\mathbf{x}_{*}).$
(16)

We omit the derivation of these expressions, as they can be found in several
standard references (Murphy, 2022, p. 84).

Viewing (14) as a function of $\mathbf{x}_{*}$, we gain the ability to reason
about the function $F$ in a Bayesian manner. An important fact about these
expressions is that the posterior expected value of the function,
$\mathbb{E}(F(\mathbf{x}_{*})|\mathbf{y})$, can be expressed as a weighted sum
of kernel functions,

$\mathbb{E}(F(\mathbf{x}_{*})|\mathbf{y})=\mathbf{k}(\mathbf{x}_{*})^{\top}(\mathbf{K}+\sigma^{2}\mathbf{I}_{N})^{-1}(\mathbf{y}-\mathbf{m})=\sum_{n=1}^{N}k(\mathbf{x}_{*},\mathbf{x}_{n})\alpha_{n},$
(17)

where

$\begin{bmatrix}\alpha_{1}\\\ \vdots\\\
\alpha_{N}\end{bmatrix}\stackrel{{\scriptstyle\Delta}}{{=}}(\mathbf{K}+\sigma^{2}\mathbf{I}_{N})^{-1}(\mathbf{y}-\mathbf{m}).$
(18)

The weighted sum in (18) is in particular useful to analyze the GPR model,
allowing us to perform calculations that would otherwise be intractable.

###  3.1 Covariance Functions

In many applications, the default kernel used in GPR is the squared
exponential (SE) kernel Rasmussen and Williams (2006), given by

$k_{\text{SE}}(\mathbf{x},\mathbf{x}^{\prime})=\sigma^{2}_{0}\exp\left(-\sum_{i=1}^{D}\frac{(x_{i}-x_{i}^{\prime})^{2}}{2\ell^{2}}\right),$
(19)

where $\ell$ and $\sigma^{2}_{0}$ are hyperparameters. A minor modification
that improves the fitting power of the SE kernel is to let each input
dimension $x_{i}$ have its own length-scale parameter $\ell_{i}$, in which
case we have the automatic relevance detection squared exponential (ARD-SE)
kernel Neal (1996):

$k_{\text{ARD-
SE}}(\mathbf{x},\mathbf{x}^{\prime})=\sigma^{2}_{0}\exp\left(-\sum_{i=1}^{D}\frac{(x_{i}-x_{i}^{\prime})^{2}}{2\ell_{i}^{2}}\right).$
(20)

###  3.2 Derivatives and Integrals of GPs

For computing the IG attributions, we will need to integrate and differentiate
GPs. We already remarked that functions sampled from a GP will be
differentiable when the GP mean and kernel functions are also sufficiently
differentiable Adler and Taylor (2007). In this way, it also makes sense to
ask what is the probability distribution of $\partial F/\partial x_{i}$ when
$F\sim\mathcal{GP}$. The answer is well known; $\partial F/\partial x_{i}$ is
also distributed according to a GP:

######  Lemma 1 (Derivatives of GPs)

If $F\sim\mathcal{GP}(m,k)$, where $m\in C^{1}(\mathbb{R}^{D})$ and $k\in
C^{2}(\mathbb{R}^{D}\times\mathbb{R}^{D})$, then

$\frac{\partial F}{\partial x_{i}}\sim\mathcal{GP}\left(\frac{\partial
m}{\partial x_{i}},\frac{\partial^{2}k}{\partial x_{i}\partial
x_{i}^{\prime}}\right).$

This result is well-established in the literature Solak et al. (2002); Le Gall
(2016), but for clarity we derive it in Appendix A.

In a similar spirit, the integral of a GP is also Gaussian, but we should be
careful to clarify what type of integration we are considering. In one-
dimension, the indefinite integral of a GP requires taking a limit that may
not be well defined, depending on the mean and covariance functions of the GP.
However, we can define the definite integral fairly easily, which turns out to
be a non-stationary GP.

######  Lemma 2 (Integral of a GP)

Suppose $f\sim\mathcal{GP}(m,k)$, where $m\in C^{0}(\mathbb{R})$ and $k\in
C^{0}(\mathbb{R}\times\mathbb{R})$, and let $a\in\mathbb{R}$ be fixed. The
definite integral $F(x)=\int_{a}^{x}f(t)dt$ is a GP in $x$,

$F(x)=\int_{a}^{x}f(t)dt\sim\mathcal{GP}(M,K),$

where

$\displaystyle M(x)$ $\displaystyle=\int_{a}^{x}m(t)dt,$ $\displaystyle
K(x,x^{\prime})$ $\displaystyle=\int_{a}^{x}\int_{a}^{x^{\prime}}k(s,t)dsdt.$

The integral of a GP could be established more generally, but the one-
dimensional case is all we require in the current work. Non-stationarity
arises from the fact that when $x=a$, the integral is exactly 0, but otherwise
the integral can accumulate variance as the distance between $x$ and $a$
increases. As before, we prove this lemma in Appendix A.

###  3.3 Spectral Approximations to GPR

Owing to the poor scalability of the GPR model to large numbers of training
samples, several approximations are often employed in practice to apply GPR
models to large data sets Hensman et al. (2013). In this section, we consider
the class of _spectral approximations_ to the GPR model.

According to Bochner’s theorem, every stationary covariance function,
expressed as
$r(\mathbf{x}-\mathbf{x}^{\prime})=k(\mathbf{x},\mathbf{x}^{\prime})$, can be
represented as the Fourier transform of a finite positive measure Stein
(1999). When this positive measure has a spectral density $S$, the Wiener-
Khintchin theorem (Rasmussen and Williams, 2006, Chapter 4) states that $r$
and $S$ are Fourier duals. Since the measure represented by $S$ is finite and
positive, it can be normalized to yield a probability density function
$p(\mathbf{v})$. The random feature GP (RFGP) model proposed by Lázaro-
Gredilla et al. (2010) samples from this probability distribution to generate
_random features_ used to approximate the full GP kernel. Under this approach,
we approximate

$k(\mathbf{x},\mathbf{x}^{\prime})\approx\boldsymbol{\phi}(\mathbf{x})^{\top}\boldsymbol{\phi}(\mathbf{x}^{\prime}),$

using $M$ samples from the power spectral density of the covariance function,
where

$\boldsymbol{\phi}(\mathbf{x})=\begin{bmatrix}\sin(\mathbf{x}^{\top}\mathbf{v}_{1})\\\
\cos(\mathbf{x}^{\top}\mathbf{v}_{1})\\\ \vdots\\\
\sin(\mathbf{x}^{\top}\mathbf{v}_{M})\\\
\cos(\mathbf{x}^{\top}\mathbf{v}_{M})\end{bmatrix},$

and

$\mathbf{v}_{m}\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}p(\mathbf{v}),\qquad
m=1,...,M.$

The Wiener-Khintchin theorem is particularly useful because it allows us to
explicitly compute the density $p$ in a number of cases. For example, the
Fourier dual of a Gaussian function
$\exp(-\mathbf{x}^{\top}\mathbf{K}^{-1}\mathbf{x})$, where $\mathbf{K}$ is a
positive definite matrix, is again a Gaussian. Thus, we can produce random
features from SE and ARD-SE kernels by sampling from Gaussian distributions.
Having this sparse approximation of the kernel function, we can derive the
predictive model for RFGPs. From Lázaro-Gredilla et al. (2010), the predictive
model is given by

$\mathbb{E}(F(\mathbf{x}_{*}))=\boldsymbol{\phi}(\mathbf{x}_{*})^{\top}\mathbf{A}^{-1}\bf{\Phi}\mathbf{y},$
(21)

and

$\text{Cov}(F(\mathbf{x}_{*}),F(\mathbf{x}_{*}^{\prime}))=\sigma^{2}_{n}\boldsymbol{\phi}(\mathbf{x}_{*})^{\top}\mathbf{A}^{-1}\boldsymbol{\phi}(\mathbf{x}_{*}^{\prime}),$
(22)

where we have defined the matrices

$\displaystyle\mathbf{A}=\boldsymbol{\Phi}\boldsymbol{\Phi}^{\top}+\frac{M\sigma^{2}_{n}}{\sigma^{2}_{0}}\mathbf{I}_{2M},$
$\displaystyle\boldsymbol{\Phi}=\left[\boldsymbol{\phi}(\mathbf{x}_{1}),\boldsymbol{\phi}(\mathbf{x}_{2}),\cdots,\boldsymbol{\phi}(\mathbf{x}_{N})\right]\in\mathbb{R}^{2M\times
N},$

for convenience. According to the Wiener-Khintchin theorem, frequency vectors
for an ARD-SE kernel, (20), can be sampled as

$\mathbf{v}_{m}\sim\mathcal{N}\left(\mathbf{0},\text{diag}(\boldsymbol{\ell})^{-2}\right),\qquad
m=1,...,M,$

where $\boldsymbol{\ell}=[\ell_{1}\cdots\ell_{D}]^{\top}$. This approach
leaves us with $\dim(\mathbf{x})+2$ hyperparameters in this model,

$\boldsymbol{\theta}=\begin{bmatrix}\sigma^{2}_{n}\\\ \sigma^{2}_{0}\\\
\boldsymbol{\ell}\end{bmatrix}.$

We close this section by noting that other variants of the RFGP exist, taking
different approaches to sample random features. In general, RFGPs are not
usually optimal, because the particular sample of frequency vectors
$\mathbf{v}_{m}$ strongly determines the quality of the approximation.
Additionally, other approaches may yield feature-based GP models with the same
form as (21). In particular, Hilbert space approximations to GPR have a
similar form to RFGPs, but these models see improvements in performance
because their features are optimally selected as opposed to being randomly
sampled Riutort-Mayol et al. (2023).

##  4 Attributions in GPR Models

In this section, we introduce our main theorem: GP models produce GP
attributions. We then apply this theorem to several common GPR models to
derive explicit expressions for their attributions.

######  Theorem 1 (Integrated gradients preserve Gaussianity)

If $F\sim\mathcal{GP}(m,k)$, where we assume $m\in C^{1}(\mathbb{R}^{D})$ and
$k\in C^{2}(\mathbb{R}^{D}\times\mathbb{R}^{D})$, then the IG attributions
follow a GP:

$\displaystyle
attr_{i}(\mathbf{x}|F)\sim\mathcal{GP}(\mu_{i},\kappa_{i}),\qquad i=1,...,D,$
(23)

where

$\displaystyle\mu_{i}(\mathbf{x})$
$\displaystyle=(x_{i}-\tilde{x}_{i})\int_{0}^{1}\frac{\partial
m(\tilde{\mathbf{x}}+t(\mathbf{x}-\tilde{\mathbf{x}}))}{\partial x_{i}}dt,$
(24) $\displaystyle\kappa_{i}(\mathbf{x},\mathbf{x}^{\prime})$
$\displaystyle=(x_{i}-\tilde{x}_{i})(x_{i}^{\prime}-\tilde{x}_{i})\int_{0}^{1}\int_{0}^{1}\frac{\partial^{2}k(\tilde{\mathbf{x}}+s(\mathbf{x}-\tilde{\mathbf{x}}),\tilde{\mathbf{x}}+t(\mathbf{x}^{\prime}-\tilde{\mathbf{x}}))}{\partial
x_{i}\partial x_{i}^{\prime}}dtds.$ (25)

Furthermore, the attributions satisfy a GP version of the completeness axiom:

$\sum_{i=1}^{D}attr_{i}(\mathbf{x}|F)=F(\mathbf{x})-F(\tilde{\mathbf{x}})\sim\mathcal{GP}(\tilde{m},\tilde{k}),$
(26)

where

$\displaystyle\tilde{m}(\mathbf{x})$
$\displaystyle=m(\mathbf{x})-m(\tilde{\mathbf{x}}),$ (27)
$\displaystyle\tilde{k}(\mathbf{x},\mathbf{x}^{\prime})$
$\displaystyle=k(\mathbf{x},\mathbf{x}^{\prime})+k(\tilde{\mathbf{x}},\tilde{\mathbf{x}})-k(\mathbf{x},\tilde{\mathbf{x}})-k(\mathbf{x}^{\prime},\tilde{\mathbf{x}}).$
(28)

The basic motivation for this theorem is that since $attr_{i}$ is a linear
operator, it preserves Gaussianity Solak et al. (2002). The proof is given in
Appendix A.4. Here we decompose the attribution operator into a sequence of
operations through which we can compute the form of the GP. Other attribution
theories, formulated in terms of linear operators as mentioned in 2.3, would
admit proofs in a similar manner. It should be noted that the attribution
functions are always _heteroscedastic_ , meaning that their variance is
location-dependent. This can be seen by observing that
$\kappa_{i}(\mathbf{x},\mathbf{x})$ generally grows as $x_{i}-\tilde{x}_{i}$
increases, but also that
$\text{Var}(attr_{i}(\mathbf{x}))=\kappa_{i}(\mathbf{x},\mathbf{x})=0$
whenever $x_{i}=\tilde{x}_{i}$.

Knowing the expressions for the mean and covariance of the GPR posterior
predictive distribution, we may derive expressions for the attributions of a
GPR model. For a bit of notation, recall from $\eqref{eq:vectorattr}$ that the
attribution vector $\mathbf{attr}_{i}(\mathbf{x}|\mathbf{k}(\mathbf{x}))$ is a
vector-valued function whose entries are the attribution operator applied to
each entry of $\mathbf{k}(\mathbf{x})$, that is,

$\mathbf{attr}_{i}(\mathbf{x}|\mathbf{k})=\begin{bmatrix}attr_{i}(\mathbf{x}|k(\mathbf{x},\mathbf{x}_{1}))\\\
\vdots\\\ attr_{i}(\mathbf{x}|k(\mathbf{x},\mathbf{x}_{N}))\end{bmatrix}.$

Using this notation, we present the attributions of the GPR predictive
posterior model in (14) via the following corollary.

######  Corollary 1 (Attributions of the GPR Posterior-Predictive)

If $F$ is distributed according to the GP posterior-predictive distribution
(14), (15), and (16), then the IG attributions $attr_{i}(\mathbf{x})$ are GPs
such that

$\mathbb{E}(attr_{i}(\mathbf{x}))=attr_{i}(\mathbf{x}|m)+\mathbf{attr}_{i}(\mathbf{x}|\mathbf{k})^{\top}(\mathbf{K}+\sigma^{2}\mathbf{I}_{N})^{-1}(\mathbf{y}-\mathbf{m}),$
(29)

and

$\displaystyle\text{Cov}(attr_{i}(\mathbf{x}),attr_{i}(\mathbf{x}^{\prime}))=(x_{i}-\tilde{x}_{i})$
$\displaystyle(x_{i}^{\prime}-\tilde{x}_{i})\int_{0}^{1}\int_{0}^{1}\frac{\partial^{2}k(\tilde{\mathbf{x}}+s(\mathbf{x}-\tilde{\mathbf{x}}),\tilde{\mathbf{x}}+t(\mathbf{x}^{\prime}-\tilde{\mathbf{x}}))}{\partial
x_{i}\partial x_{i}^{\prime}}dsdt$ (30)
$\displaystyle-\mathbf{attr}_{i}(\mathbf{x}|\mathbf{k})^{\top}(\mathbf{K}+\sigma^{2}\mathbf{I})^{-1}\mathbf{attr}_{i}(\mathbf{x}^{\prime}|\mathbf{k}),$

The tractability of the expressions in the corollary depends precisely on the
choice of mean and covariance functions used in the GP prior. In the following
subsections, we compute these results for some noteworthy GPR models.

###  4.1 SE and ARD-SE Kernels

The attributions for the GPR model can be derived in closed form when we
assume an ARD-SE kernel (20), although the derivation can be somewhat tedious.
Assuming a zero mean prior over $F$ with ARD-SE kernel, the attributions of
the GPR model are given by

$\displaystyle\mathbb{E}(attr_{i}(\mathbf{x}))$
$\displaystyle=\sum_{n=1}^{N}\alpha_{n}A_{n,i}(\mathbf{x}),$ (31)
$\displaystyle\text{Var}(attr_{i}(\mathbf{x}))$
$\displaystyle=B_{i}(\mathbf{x})-\sum_{n=1}^{N}\sum_{m=1}^{N}\Lambda_{n,m}A_{n,i}(\mathbf{x})A_{m,i}(\mathbf{x}),$
(32)

where we have defined

$\displaystyle A_{n,i}(\mathbf{x})$
$\displaystyle=\frac{e^{\frac{-a-b-c}{2}}\left[\sqrt{2\pi}e^{\frac{(2a+b)^{2}}{8a}}\left(\text{erf}\left(\frac{b}{2\sqrt{2a}}\right)-\text{erf}\left(\frac{2a+b}{2\sqrt{2a}}\right)\right)(bd-2af)+4\sqrt{a}d\left(e^{\frac{a+b}{2}}-1\right)\right]}{4a^{3/2}},$
$\displaystyle B_{i}(\mathbf{x})$
$\displaystyle=(x_{i}-\tilde{x}_{i})^{2}\left(\frac{\sqrt{2\pi}\text{erf}\left(\frac{\sqrt{a}}{\sqrt{2}}\right)(aw+v)}{a^{3/2}}-\frac{2e^{-\frac{a}{2}}\left(e^{a/2}-1\right)(aw+2v)}{a^{2}}\right),$

where $\text{erf}(\cdot)$ is the error function and

$\displaystyle a$
$\displaystyle=(\mathbf{x}-\tilde{\mathbf{x}})^{\top}\mathbf{L}^{-2}(\mathbf{x}-\tilde{\mathbf{x}}),\quad$
$\displaystyle\mathbf{L}=\text{diag}(\ell_{1},\cdots,\ell_{D}),$
$\displaystyle b$
$\displaystyle=2(\mathbf{x}-\tilde{\mathbf{x}})^{\top}\mathbf{L}^{-2}(\tilde{\mathbf{x}}-\mathbf{x}_{n}),\quad$
$\displaystyle
v=-\sigma^{2}_{0}\frac{(x_{i}-\tilde{x}_{i})^{2}}{\ell_{i}^{4}},$
$\displaystyle c$
$\displaystyle=(\tilde{\mathbf{x}}-\mathbf{x}_{n})^{\top}\mathbf{L}^{-2}(\tilde{\mathbf{x}}-\mathbf{x}_{n}),$
$\displaystyle w=\frac{\sigma^{2}_{0}}{\ell_{i}^{2}},$ $\displaystyle d$
$\displaystyle=-\sigma^{2}_{0}\frac{(x_{i}-\tilde{x}_{i})^{2}}{\ell_{i}^{2}},$
$\displaystyle[\Lambda_{n,m}]=(\mathbf{K}+\sigma^{2}_{n}\mathbf{I})^{-1},$
$\displaystyle f$
$\displaystyle=-\sigma^{2}_{0}\frac{(x_{i}-\tilde{x}_{i})(\tilde{x}_{i}-x_{n,i})}{\ell_{i}^{2}},$
$\displaystyle[\alpha_{n}]=(\mathbf{K}+\sigma^{2}_{n}\mathbf{I})^{-1}\mathbf{y},$

are constants which implicitly depend upon $\mathbf{x}$, $\tilde{\mathbf{x}}$,
$\mathbf{x}_{n}$ and $i$. The derivation of these expressions is provided in
Appendix B.

###  4.2 Random Feature GPs

Based on (21), it is easy to see that IG attributions for the RFGP model can
also be derived. The mean attribution and variance can be given as

$\displaystyle\mathbb{E}(attr_{i}(\mathbf{x}))$
$\displaystyle=(x_{i}-\tilde{x}_{i})\boldsymbol{\zeta}(\mathbf{x})\mathbf{A}^{-1}\boldsymbol{\Phi}\mathbf{y},$
(33) $\displaystyle\text{Var}(attr_{i}(\mathbf{x}))$
$\displaystyle=(x_{i}-\tilde{x}_{i})^{2}\sigma^{2}_{n}\boldsymbol{\zeta}(\mathbf{x})^{\top}\mathbf{A}^{-1}\boldsymbol{\zeta}(\mathbf{x}),$

where for clarity, we introduced the notation

$\boldsymbol{\zeta}(\mathbf{x})\stackrel{{\scriptstyle\Delta}}{{=}}\begin{bmatrix}\frac{v_{1,i}}{\mathbf{v}_{1}^{\top}(\mathbf{x}-\tilde{\mathbf{x}})}(\cos(\mathbf{v}_{1}^{\top}\mathbf{x})-\cos(\mathbf{v}_{1}^{\top}\tilde{\mathbf{x}}))\\\
\frac{v_{1,i}}{\mathbf{v}_{1}^{\top}(\mathbf{x}-\tilde{\mathbf{x}})}(\sin(\mathbf{v}_{1}^{\top}\mathbf{x})-\sin(\mathbf{v}_{1}^{\top}\tilde{\mathbf{x}}))\\\
\vdots\\\
\frac{v_{M,i}}{\mathbf{v}_{M}^{\top}(\mathbf{x}-\tilde{\mathbf{x}})}\left(\cos(\mathbf{v}_{M}^{\top}\mathbf{x})-\cos(\mathbf{v}_{M}^{\top}\tilde{\mathbf{x}})\right)\\\
\frac{v_{M,i}}{\mathbf{v}_{M}^{\top}(\mathbf{x}-\tilde{\mathbf{x}})}\left(\sin(\mathbf{v}_{M}^{\top}\mathbf{x})-\sin(\mathbf{v}_{M}^{\top}\tilde{\mathbf{x}})\right)\end{bmatrix}\in\mathbb{R}^{2M\times
1},$

and

$\mathbf{A}=\boldsymbol{\Phi\Phi}^{\top}+\frac{M\sigma^{2}_{n}}{\sigma^{2}_{0}}\mathbf{I}_{2M}.$

The derivation of these expressions is provided in Appendix B.

###  4.3 Comments on Other GP Models

While we only provided closed-form results for the SE and ARD-SE kernels,
Theorem 1 applies generally to all GPR models. We now briefly discuss how
these results extend to other GP models.

####  4.3.1 Additive GPs

In GP kernel design, it is common to add and multiply kernels to produce new
composite kernels (Gibbs, 1998, p.19). If $\mathcal{GP}(m,k)$ and
$\mathcal{GP}(m^{\prime},k^{\prime})$ are two GPs, then the direct sum GP is
defined to be $\mathcal{GP}(m+m^{\prime},k+k^{\prime})$. Using Corollary 1,
analytical expressions for the direct sum GP can be obtained whenever the
component kernels $k$ and $k^{\prime}$ are tractable. Duvenaud et al. (2011)
defines the _full additive kernel_ by summing over all possible combinations
of kernels that one may obtain through including or exclusion of features.
Given a choice of base kernel $k_{i}$ for each input feature $x_{i}$, the full
additive kernel is defined to be

$k_{\textrm{add}}(\mathbf{x},\mathbf{x}^{\prime})=\sum_{d=1}^{D}\sum_{1\leq
i_{1}<\cdots<i_{d}\leq
D}\sigma^{2}_{d}\prod_{i=i_{1},...,i_{d}}k_{i}(x_{i},x_{i}^{\prime}),$

where $\sigma^{2}_{d},d=1,...,D,$ are hyperparameters. Because of the
linearity property of the attribution operator, attributions can be written in
closed form for the full additive kernel whenever the integrals of the base
kernel are tractable.

####  4.3.2 Nonstationary GPs

_Nonstationary GPs_ are GPs whose covariance functions are not translation-
invariant, and thus may exhibit different amplitudes, length scales and other
properties depending on one’s location in the input space. Several approaches
to using nonstationary GPs for inference exist: scaling the kernel by
parametric functions to introduce heteroscedasticity (Gibbs, 1998, pp.16-18),
deep kernel learning Wilson et al. (2016) and deep GPs Damianou and Lawrence
(2013). Theorem 1 applies to directly nonstationary GPs with no modification.

##  5 Experiments

We now validate our theory using several real world and simulated data sets.
We provide several representative examples here, and we show additional
results in the appendix.

###  5.1 Simulated Data: Heteroscedasticity of the Attribution GPs

In this example, we examine how the variance of an attribution grows as we
move farther away from the baseline. To visualize this in a controlled manner,
we introduce a synthetic data set to validate our theory. Consider the
following generative model, where inputs $X_{1}$ and $X_{2}$ combine to
produce an output $Y$:

$\displaystyle X_{1}$ $\displaystyle\sim\mathcal{U}(0,10),$ (34)
$\displaystyle X_{2}$ $\displaystyle\sim\mathcal{U}(0,10),$ $\displaystyle
N_{Y}$ $\displaystyle\sim\mathcal{N}(0,1),$ $\displaystyle Y$
$\displaystyle=\sin(X_{1})\sin(2X_{2})+\frac{1}{2}N_{Y}.$

We sample 500 input-output pairs to train a GPR model with ARD-SE kernel. The
mean predictions of this model are visualized as a function of $X_{1}$ and
$X_{2}$ in Figure 1. To study attributions, we will consider the line segment
parameterized by

$\boldsymbol{\beta}(t)=t\begin{bmatrix}1\\\ 1\end{bmatrix},\qquad 0\leq t\leq
10.$ (35)

As we increase $t$, $\boldsymbol{\beta}(t)$ moves along the diagonal of the
square $[0,10]\times[0,10]$. If we set $\tilde{\mathbf{x}}=\mathbf{0}$, then
the distance between $\boldsymbol{\beta}(t)$ and $\tilde{\mathbf{x}}$
increases with $t$, and we expect that the attributions will increase in
variance as we move farther away from the baseline.

In Figure 1, we plot $attr_{i}(\boldsymbol{\beta}(t)|F)$ against $t$, for
$i=1,2$. We observe that for both attribution functions, the variance of the
attribution increases with $t$. As $t\to 0$, both the mean and variance of the
attributions tend towards zero. The rate at which the variance of the
attribution grows depends on the uncertainty of the GPR predictive model.

Figure 1: A demonstration that the uncertainty in an attribution grows with
the distance from the baseline, using data simulated with the model in (34).
The baseline is $\tilde{\mathbf{x}}=\mathbf{0}$ and we compute attributions
for predictions along the path $\boldsymbol{\beta}(t)=t[1,1]^{\top}$.

###  5.2 Breast Cancer Prognosis: A General Example

In this case study, we consider the breast cancer prognostic data set of
Wolberg et al. (1995), available through the UCI Machine Learning Repository.
In this data set, there 30 features (10 measurements from 3 different
recording sites) obtained from digitized images of nucleation sites in breast
tissue. The goal of analyzing this data set is to predict either the
recurrence time, or the disease-free time, of a patient after being screened.
To this end, we trained a GPR model with an ARD-SE kernel to predict the
disease-free time of patients in this data set. The goal of XAI in this
setting is to understand, given our GPR predictive model, what features were
contributing to the prediction, or in other words, which measurements at which
recording sites best informed our prediction. We define the baseline
$\tilde{\mathbf{x}}$ to be average feature vector across all patients. In
practice, feature attributions are sensitive to the choice of baseline point,
and careful consideration is required to select a meaningful reference. As an
example, we adopt the relatively simple approach of defining the baseline to
be the average patient, and thus the IG method will attribute the change of a
patient’s disease-free time according to how an individual patient’s
measurements differ from the average patient.

In Figure 2(a), we visualize all attributions for a single patient. In Figure
2(b), we visualize how two of these attributions vary across 50 different
patients. Attributions in this case measure how variations in each feature
contribute to the difference between the prediction for a given patient versus
a reference patient with the baseline features. Thus, the value of a feature
could contribute positively or negative towards this difference from the
baseline. We observe that the GPR model with ARD-SE kernel prefers to
attribute its predictions to only a few features, because the ARD-SE kernel
automatically performs feature selection. Across patients, we observe that a
feature might contribute strongly or weakly to a prediction. As the magnitude
of an attribution increases in a particular feature, so too does the
uncertainty in the attribution, which is something we observed in the Bayesian
linear regression model of Section 2.1

((a)) For 6 patients, we computed the feature attributions and their
uncertainties. Feature names are only shown for consistently significant
contributors to the prediction. The shaded region expresses the mean
attribution and the error bars show one standard deviation in uncertainty.

((b)) Two features are selected to visualize how confidence and mean
attribution vary across 50 patients. As in Figure 1, the midline represents
the mean and the height of the cross is one standard deviation of
$attr_{i}(\mathbf{x}|F)$ for each patient.

Figure 2: Feature attribution for patients from the breast cancer data set.
Predictions were made using a GPR model with ARD-SE kernel.

###  5.3 Taipei Housing Data: Demonstrating Affine Scale Invariance

In this example, we consider real estate data in New Taipei City obtained from
the UCI Machine Learning Repository mis (2018) in order to demonstrate a
property of IG attributions called affine scale invariance (ASI). In this data
set, there are 6 features used to predict the value of a house. After training
a GPR model with ARD-SE kernel to predict the house value, we wish to
attribute these predictions to each of the input features to understand how
this information contributed to the predicted value. We take the baseline
$\tilde{\mathbf{x}}$ to be the average across all houses, and consider a
prediction for a new house excluded from the training set.

As in many ML applications, we normalized the input data before training the
GPR model. However, we would like to know if this procedure has distorted our
predicted attributions. Fortunately, the ASI property says that the IG
attributions are invariant to scaling and translations in the input feature
space Lundstrom et al. (2022). To formalize this statement, let us consider a
mapping $T(\mathbf{x})=\mathbf{\lambda}\odot\mathbf{x}+\mathbf{v}$, where
$\odot$ represents a Hadamard (entrywise) product of vectors and all
$\lambda_{i}\neq 0$. The mapping $T$ scales each coordinate $x_{i}$ by
$\lambda_{i}$ and then translates by $v_{i}$. ASI asserts that

$attr_{i}(T(\mathbf{x}),T(\tilde{\mathbf{x}})|F\circ
T^{-1})=attr_{i}(\mathbf{x},\tilde{\mathbf{x}}|F),\qquad i=1,...,D,$

for any such transformation $T$.

An important note about this result is that while the attributions do not
theoretically depend upon normalization, the training of a ML model still may
benefit from normalization, and this is especially true when using GPR. Thus,
the ASI property suggests that normalization doesn’t distort attributions by
itself, but rather any change in the resulting attributions will be due to the
change in how the ML model tunes it’s parameters.

In Figure 3, we consider the attributions for a particular house both with and
without normalization. We observe that although the input feature vectors
appear noticeably different, the resulting attributions are identical. The
normalized model also retains a degree of interpretability; we see that the
given house has an above average distance to the MRT (subway station), and as
a result this contributed negatively to the house value.

Figure 3: Feature attributions for a prediction using the Taipei Housing Data.

###  5.4 Wine: Approximating the Integral in IG Attributions

In this experiment, we used the wine quality data set of Cortez et al.
(2009a), available through the UCI Machine Learning Repository Cortez et al.
(2009b). The goal of feature attribution with this data set is to attribute
the predicted wine quality to the 11 physiochemical features that were used to
make the prediction. Wines are scored on a scale of 1 to 10, with $5$ being an
average quality wine. Prior to inference, the data set was randomly shuffled,
the input features were normalized, and the baseline $\tilde{\mathbf{x}}$ was
chosen to be average of all wines which had average quality. We consider
prediction for a particular instance $\mathbf{x}$, whose quality was above-
average, which we omitted from the training set.

In this experiment, we examine the effect of approximating the integral in the
IG attribution formula (8). Sundararajan et al. (2017) do not compute this
integral directly, but rather they approximate it using a summation:

$(x_{i}-\tilde{x}_{i})\int_{0}^{1}\frac{\partial
F(\tilde{\mathbf{x}}+t(\mathbf{x}-\tilde{\mathbf{x}}))}{\partial
x_{i}}dt\approx\frac{(x_{i}-\tilde{x}_{i})}{L}\sum_{l=1}^{L}\frac{\partial
F(\tilde{\mathbf{x}}+(l/L)(\mathbf{x}-\tilde{\mathbf{x}}))}{\partial x_{i}}.$
(36)

This technique is known as the _right-hand rule_ for numerical integration.
There exist other approaches for numerical approximation of the integral,
called quadrature rules Heath (2018), which we compare to the right-hand rule
in Appendix C.1. Similar to the exact result, the approximation yields a
Gaussian distribution for the attribution Seitz (2022). This approximation
closely approximates the distribution of the exact result in Theorem 1. The
quality of the approximation depends on the number of partitions $L$ used to
compute the sum.

In Figure 4, we visualize the convergence of the approximation to the exact
GPR result as $L$ is increased. We consider the wine quality example as well
as simulated data generated by (34), to assess how this convergence might
depend on the data set and the learned model. When using the right-hand rule,
$L$ is also the number of function evaluations needed to compute the
approximation. In general, the exact result only needs to be evaluated once,
so using a large $L$ to compute the attribution might become expensive if
computed naively. How large $L$ should be to obtain a good approximation will
depend on the data set and the complexity of the learned model. We also note
that the approximation quality depends on the proximity to
$\tilde{\mathbf{x}}$, as points farther from the baseline will require more
partitions to achieve the same resolution as points close to
$\tilde{\mathbf{x}}$.

((a)) Results from the Wine quality data set.

((b)) Results from the simulated system, (34).

Figure 4: Convergence of attributions with the approximated integral, (36), to
the exact attributions, as the number of partitions $L$ increases.

###  5.5 Wine: Comparing the Attributions of RFGPs and Exact GPR

Much like the previous section, we used the wine quality data set to contrast
an approximation with the exact attributions in Theorem 1. However, in this
case we consider the RFGP approximation to the GPR kernel function.

In Figure 5, we consider a prediction from the wine quality data set, and we
visualize the distribution of the attribution to a feature using the exact GPR
result, and then how the RFGP approximation compares. We observe that as $M$
increases, the individual samples from the RFGP model become more similar, but
often there is still a bias in the RFGP attributions, although this bias
appears to decrease as we ad frequencies to the model. In addition to
individual RFGPs, we consider a marginalized RFGP model, under which we
averaged the distributions of 500 separate RFGP models.

Figure 5:  Attributions to alcohol in the Wine quality data set. As the number
of frequencies used in the RFGP approximation, $M$, increases, the densities
of the RFGP attributions more closely align with the exact GPR attributions.
We show both individual realizations of the RFGP model (the blue and red
curves), and a marginalized model which integrates over 500 independent RFGP
models (the red shaded region).

##  6 Conclusion

In this work, we develop a bridge between the theory of feature attribution
and GPR. By representing attributions as distributions, we can quantify a
model’s confidence in its feature attributions. When the selected GPR model
admits analytically tractable attributions, using closed-form expressions for
the attributions is more computationally efficient, and it mitigates the need
for computationally expensive or otherwise inaccurate approximations. While
non-parametric methods like GPR, at face value, offer a non-interpretable
approach to ML, our study of feature attribution may provide more trust in GPR
models. Furthermore, these results will extend the theoretical and practical
scope of GPR models to more XAI domains.

Acknowledgments and Disclosure of Funding

This work was supported by the National Science Foundation under Awards
2021002 and 2212506.

##  Appendix A Proofs

In this appendix, we provide mathematical proofs for the lemmas and main
theorem of the paper. We divide Appendix A. into four sections: In Section
A.1, we introduce some additional technical lemmas. In Section A.2, we derive
the expression for the derivatives of a GP, and in Section A.3, we do the same
for definite integrals of a GP. Finally, in Section A.4 we prove that the IG
attributions of a GPR model are GPs.

We use the notation $C^{\ell}(\mathcal{X})$ to denote the set of all
$\ell$-times continuously differentiable functions $\mathcal{X}\to\mathbb{R}$,
where $\ell=1,2,...$ or $\ell=\infty$. We may write $C^{\ell}$ when the
functions’ domain is clear from context.

###  A.1 Operations that Preserve Gaussianity

We say that a random process $(X(t))_{t\in T}$, with index set $T$, is a GP if
for every finite set of indices $t_{1},...,t_{D}$, the vector
$(X(t_{1}),...,X(t_{D}))$ has a multivariate (possibly degenerate) Gaussian
distribution. We use the notation $X(t)\sim\mathcal{GP}(m,k)$, where
$m(t)=\mathbb{E}(X(t))$ and $k(t,t^{\prime})=\text{Cov}(X(t),X(t^{\prime}))$.
From this definition, we can immediately show that scaling and precompositions
preserve Gaussianity:

######  Lemma 3 (Scaling preserves Gaussianity)

Let $X(t)\sim\mathcal{GP}(m,k)$ be a GP with index set $T$ and let
$g:T\to\mathbb{R}$ be a function. Then $g(t)X(t)$ is a GP, and for any finite
set of indices $t_{1},...,t_{D}$, we have that

$\begin{bmatrix}g(t_{1})X(t_{1})\\\ \vdots\\\
g(t_{D})X(t_{D})\end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}g(t_{1})m(t_{1})\\\
\vdots\\\
g(t_{D})m(t_{D})\end{bmatrix},\begin{bmatrix}g(t_{1})^{2}k(t_{1},t_{1})&\cdots&g(t_{1})g(t_{D})k(t_{1},t_{D})\\\
\vdots&&\vdots\\\
g(t_{D})g(t_{1})k(t_{D},t_{1})&\cdots&g(t_{D})^{2}k(t_{D},t_{D})\end{bmatrix}\right).$
(37)

Proof  For any set of indices $t_{1},...,t_{D}$, the pointwise scaling by
$g(t)$ can be represented in matrix form,

$\begin{bmatrix}g(t_{1})X(t_{1})\\\ \vdots\\\
g(t_{D})X(t_{D})\end{bmatrix}=\begin{bmatrix}g(t_{1})&&\\\ &\ddots&\\\
&&g(t_{D})\end{bmatrix}\begin{bmatrix}X(t_{1})\\\ \vdots\\\
X(t_{D})\end{bmatrix}.$

Multiplying a Gaussian vector by a matrix preserves Gaussianity, so the GP
property is upheld for the scaled process.

######  Lemma 4 (Precompositions preserve Gaussianity)

Let $X(t)\sim\mathcal{GP}(m,k)$ be a GP with index set $T$, and let
$\gamma:S\to T$ be a function between sets. Then $X(\gamma(s))$ is a GP
indexed by $S$, and for any finite set of indices $s_{1},...,s_{D}$,

$\begin{bmatrix}X(\gamma(s_{1}))\\\ \vdots\\\
X(\gamma(s_{D}))\end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}m(\gamma(s_{1}))\\\
\vdots\\\
m(\gamma(s_{D}))\end{bmatrix},\begin{bmatrix}k(\gamma(s_{1}),\gamma(s_{1}))&\cdots&k(\gamma(s_{1}),\gamma(s_{D}))\\\
\vdots&&\vdots\\\
k(\gamma(s_{D}),\gamma(s_{1}))&\cdots&k(\gamma(s_{D}),\gamma(s_{D}))\end{bmatrix}\right).$
(38)

Proof  Let $t_{i}=\gamma(s_{i})$ for $i=1,...,D$. Then $t_{1},...,t_{D}$ is a
set of indices in $T$ and hence

$\begin{bmatrix}X(\gamma(s_{1}))\\\ \vdots\\\
X(\gamma(s_{D}))\end{bmatrix}=\begin{bmatrix}X(t_{1})\\\ \vdots\\\
X(t_{D})\end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}m(t_{1})\\\ \vdots\\\
m(t_{D})\end{bmatrix},\begin{bmatrix}k(t_{1},t_{1})&\cdots&k(t_{1},t_{D})\\\
\vdots&&\vdots\\\ k(t_{D},t_{1})&\cdots&k(t_{D},t_{D})\end{bmatrix}\right).$

The only thing to notice here is that if $\gamma$ is not injective (one-to-
one), then the new GP is degenerate, regardless if the original GP was non-
degenerate.

###  A.2 Derivatives of GPs

As noted in the main text, differentiation and integration also preserve
Gaussianity. To show this, we use the following proposition from (Le Gall,
2016, p.3):

###### Proposition 1

Let $(X_{n})_{n\geq 1}$ be a sequence of real random variables such that, for
every $n\geq 1$, $X_{n}$ follows a $\mathcal{N}(m_{n},\sigma^{2}_{n})$
distribution. Suppose that $X_{n}$ converges in $L^{2}$ to $X$. Then:

  1. 1.

The random variable $X$ follows the $\mathcal{N}(m,\sigma^{2})$ distribution,
where $m=\lim_{n\to\infty}m_{n}$ and
$\sigma^{2}=\lim_{n\to\infty}\sigma^{2}_{n}$.

  2. 2.

This convergence also holds in all $L^{p}$ spaces, $1\leq p<\infty$.

We note that the $L^{2}$ convergence is automatically satisfied whenever the
sequences $m_{n}$ and $\sigma^{2}_{n}$ converge to something reasonable (in
essence, when $0<\sigma^{2}<\infty$). Using this proposition, we can now prove
results for derivatives and integrals of GPs.

Proof  of Lemma 1 To start, we recall the typical definition of a partial
derivative:

$\frac{\partial F}{\partial x_{i}}=\lim_{h\to
0}\frac{F(\mathbf{x}+h\mathbf{e}_{i})-F(\mathbf{x})}{h},$

where $\mathbf{e}_{i}$ is the unit vector in the $i$-th coordinate direction.
To take the limit, define a sequence $h_{n}=1/n$, for $n\in\mathbb{Z}_{>0}$,
and consider the sequence of functions
$F_{n}(\mathbf{x})=(F(\mathbf{x}+h_{n})-F(\mathbf{x}))/h_{n}$. We see that
each $F_{n}$ is a GP, and thus we have a sequence of GPs,
$F_{n}\sim\mathcal{GP}(m_{i,n},k_{i,n})$. By Proposition 1, a sequence of
Gaussian variables converges to a Gaussian if the mean and variance converge.
In this case,

$\lim_{n\to\infty}m_{i,n}=\lim_{n\to\infty}\mathbb{E}\left(\frac{F(\mathbf{x}+\frac{1}{n}\mathbf{e}_{i})-F(\mathbf{x})}{1/n}\right)=\lim_{n\to\infty}\frac{m(\mathbf{x}+\frac{1}{n}\mathbf{e}_{i})-m(\mathbf{x})}{1/n}=\frac{\partial
m}{\partial x_{i}},$

which exists since $m\in C^{1}$. For the covariance, we have that

$k_{i,n}(\mathbf{x},\mathbf{x}^{\prime})=\text{Cov}\left(\frac{F(\mathbf{x}+\frac{1}{n}\mathbf{e}_{i})-F(\mathbf{x})}{1/n},\frac{F(\mathbf{x}^{\prime}+\frac{1}{n}\mathbf{e}_{i})-F(\mathbf{x}^{\prime})}{1/n}\right),$

which can then be expressed as

$k_{i,n}(\mathbf{x},\mathbf{x}^{\prime})=\frac{k(\mathbf{x}+\frac{1}{n}\mathbf{e}_{i},\mathbf{x}^{\prime}+\frac{1}{n}\mathbf{e}_{i})-k(\mathbf{x},\mathbf{x}^{\prime}+\frac{1}{n}\mathbf{e}_{i})-k(\mathbf{x}+\frac{1}{n}\mathbf{e}_{i},\mathbf{x}^{\prime})+k(\mathbf{x},\mathbf{x}^{\prime})}{n^{2}}.$

Taking the limit then yields

$\lim_{n\to\infty}k_{i,n}(\mathbf{x},\mathbf{x}^{\prime})=\frac{\partial^{2}k}{\partial
x_{i}\partial x_{i}^{\prime}},$

which exists because $k\in C^{2}$. Thus, by Proposition 1 we have that
$\partial F/\partial x_{i}$ is Gaussian whose mean and covariance are
specified by the functions $\partial m/\partial x_{i}$ and
$\partial^{2}k/\partial x_{i}\partial x_{i}^{\prime}$. To conclude the proof,
notice that for any finite set of points $\mathbf{x}_{1},...,\mathbf{x}_{L}$,
the vector $[F(\mathbf{x}_{1})\cdots F(\mathbf{x}_{L})]$ will again be
multivariate Gaussian with mean and covariance controlled by these functions.
Since this property holds arbitrarily, the Kolmogorov extension theorem
asserts that a GP exists with the desired mean and covariance functions.

We might say that a GP is $\ell$-times continuously differentiable if the mean
and covariance functions are $C^{\ell}$ and $C^{2\ell}$, respectively, since
this ensures that functions sampled from the GP are $C^{\ell}$ almost-surely.
Using this notion of differentiability, we can discuss how scaling and
precomposition affect differentiability. The equations (37) and (38) show that
the output mean and covariance functions are $C^{\ell}$ whenever $g$, $\gamma$
and the original $m$ and $k$ functions are $C^{\ell}$. Hence, scaling and
precompositions by smooth functions preserve both smoothness and Gaussianity.
Lemma 1 shows that differentiating a $C^{\ell}$ GP yields a $C^{\ell-1}$ GP.

###  A.3 Integrals of GPs

The integration of a GP is established using similar techniques to
differentiation. For computing attributions, we only need to know how to
compute a one-dimensional integral, over a line segment.

Proof of Lemma 2 Since the integral is being taken over a line segment, it
suffices to consider the Riemannian definition of the integral as the limit
over increasingly fine partitions of the line. Hence we define

$F_{n}(x)=\sum_{i=0}^{n}\frac{1}{n}f\left(a+\frac{i(x-a)}{n}\right),$

and $F(x)$ to be the limit of $F_{n}(x)$ as $n\to\infty$. Similar in the proof
for derivatives of a GP, we need to compute the mean and covariance functions
$F_{n}(x)$, show that they converge, and apply Proposition 1. For the mean,
the summation becomes

$\mathbb{E}(F_{n}(x))=\sum_{i=0}^{n}\frac{\mathbb{E}(f(a+i(x-a)/n))}{n}=\sum_{i=0}^{n}\frac{m(a+i(x-a)/n)}{n},$

and since $m$ is continuous on the interval, the Riemann sum converges to the
integral. Thus,

$\mathbb{E}(F_{n}(x))\to\int_{a}^{x}m(t)dt=M(x),$

as $n\to\infty$, and we denote the result as a function $M$ of $x$. For the
covariance, we have that

$\displaystyle\text{Cov}\left(F_{n}(x),F_{n}(x^{\prime})\right)$
$\displaystyle=\text{Cov}\left(\sum_{i=0}^{n}\frac{f(a+i(x-a)/n)}{n},\sum_{j=0}^{n}\frac{f(a+j(x^{\prime}-a)/n)}{n}\right)$
$\displaystyle=\sum_{i=0}^{n}\sum_{j=0}^{n}\frac{1}{n^{2}}\text{Cov}\left(f\left(a+\frac{i(x-a)}{n}\right),f\left(a+\frac{j(x^{\prime}-a)}{n}\right)\right)$
$\displaystyle=\sum_{i=0}^{n}\sum_{j=0}^{n}\frac{1}{n^{2}}k\left(a+\frac{i(x-a)}{n},a+\frac{j(x^{\prime}-a)}{n}\right),$

which converges, again because $k$ is continuous, to the double integral,
which we define to be $K$, i.e.,

$\lim_{n\to\infty}\text{Cov}\left(F_{n}(x),F_{n}(x^{\prime})\right)=\int_{a}^{x}\int_{a}^{x^{\prime}}k(s,t)dsdt=K(x,x^{\prime}),$

as $n\to\infty$. Since the mean and covariance functions converge, we can
conclude the convergence of $F_{n}(x)\to F(x)$ for a given point $x$. As in
the proof of Lemma 1, for any finite set of points $x_{1},...,x_{L}$, we want
to show that the joint distribution of $F(x_{1}),...,F(x_{L})$ coincides with
the Gaussian distribution specified by the functions $M$ and $K$. Notice that
for fixed $n$,

$\mathbf{F}_{n}(x_{1},...,x_{L})\stackrel{{\scriptstyle\Delta}}{{=}}\sum_{i=1}^{n}\frac{1}{n}\begin{bmatrix}f\left(a+\frac{i(x_{1}-a)}{n}\right)\\\
\vdots\\\
f\left(a+\frac{i(x_{L}-a)}{n}\right)\end{bmatrix}=\begin{bmatrix}\sum_{i=1}^{n}\frac{1}{n}f\left(a+\frac{i(x_{1}-a)}{n}\right)\\\
\vdots\\\
\sum_{i=1}^{n}\frac{1}{n}f\left(a+\frac{i(x_{L}-a)}{n}\right)\end{bmatrix}=\begin{bmatrix}F_{n}(x_{1})\\\
\vdots\\\ F_{n}(x_{L})\end{bmatrix},$

and thus

$\lim_{n\to\infty}\mathbf{F}_{n}(x_{1},...,x_{L})=\lim_{n\to\infty}\begin{bmatrix}F_{n}(x_{1})\\\
\vdots\\\ F_{n}(x_{L})\end{bmatrix}=\begin{bmatrix}F(x_{1})\\\ \vdots\\\
F(x_{L})\end{bmatrix}.$

From the above calculations, we then know that
$\lim_{n\to\infty}\mathbf{F}_{n}(x_{1},...,x_{L})$ obeys a multivariate
Gaussian law with mean and covariance determined by $M$ and $K$, respectively.
By the Kolmogorov extension theorem, there exists a GP determined by these
functions.

###  A.4 Attributions of a GP

Given the previous lemmas, we are now equipped to compute the IG attributions
of a GPR model and prove the theorem in Section 4. Let $F(\mathbf{x})$ be a GP
over $\mathbb{R}^{D}$. We fix a baseline point $\tilde{\mathbf{x}}$ and recall
the definition of the $attr_{i}$ operator:

$attr_{i}(\mathbf{x}|F)=(x_{i}-\tilde{x}_{i})\int_{0}^{1}\frac{\partial
F(\tilde{\mathbf{x}}+t(\mathbf{x}-\tilde{\mathbf{x}}))}{\partial x_{i}}dt.$

The attribution is defined by four sequential operations,

  1. 1.

Taking the partial derivative $\partial/\partial x_{i}$ of $F$,

  2. 2.

Precomposition of $\partial F/\partial x_{i}$ with the function
$\boldsymbol{\gamma}(t)=\tilde{\mathbf{x}}+t(\mathbf{x}-\tilde{\mathbf{x}})$,

  3. 3.

Integration from $t=0$ to $t=1$, and

  4. 4.

Scaling the result by $(x_{i}-\tilde{x})$.

Together, we observe that these operations all preserve Gaussianity, and thus
to find the attribution GP it suffices to compute its mean and covariance,
whenever they exist. As a simple condition to ensure existence, we require
that $m\in C^{1}(\mathbb{R}^{D})$ and $k\in
C^{2}(\mathbb{R}^{D},\mathbb{R}^{D})$, which guarantees that we can
differentiate the GP and later integrate it.

Proof of Theorem 1 The existence of the attribution GP follows directly from
the lemmas. We show the computation for the reader’s convenience, implicitly
invoking Lemmas 1, 2, 3, and 4:

$\displaystyle\mathbb{E}(F(\mathbf{x}))$ $\displaystyle=m(\mathbf{x}),$
$\displaystyle\mathbb{E}\left(\frac{\partial F(\mathbf{x})}{\partial
x_{i}}\right)$ $\displaystyle=\frac{\partial m(\mathbf{x})}{\partial x_{i}},$
$\displaystyle\mathbb{E}\left(\frac{\partial
F(\boldsymbol{\gamma}(t))}{\partial x_{i}}\right)$
$\displaystyle=\frac{\partial m(\boldsymbol{\gamma}(t))}{\partial x_{i}},$
$\displaystyle\mathbb{E}\left(\int_{0}^{1}\frac{\partial
F(\boldsymbol{\gamma}(t))}{\partial x_{i}}dt\right)$
$\displaystyle=\int_{0}^{1}\frac{\partial
m(\tilde{\mathbf{x}}+t(\mathbf{x}-\tilde{\mathbf{x}}))}{\partial x_{i}}dt,$
$\displaystyle\mathbb{E}(attr_{i}(\mathbf{x}|F))$
$\displaystyle=\mathbb{E}\left((x_{i}-\tilde{x}_{i})\int_{0}^{1}\frac{\partial
F(\boldsymbol{\gamma}(t))}{\partial x_{i}}dt\right)$
$\displaystyle=(x_{i}-\tilde{x}_{i})\int_{0}^{1}\frac{\partial
m(\tilde{\mathbf{x}}+t(\mathbf{x}-\tilde{\mathbf{x}}))}{\partial x_{i}}dt.$

The resulting mean attribution function is defined since $m\in C^{1}$ implies
that $\partial m/\partial x_{i}$ is $C^{0}$. Since the derivative is
continuous, it is also integrable over the interval $[0,1]$. The calculation
for the covariance function is more tedious, but follows the same logic. Let
us define
$\boldsymbol{\beta}(\mathbf{x},t)=\tilde{\mathbf{x}}+t(\mathbf{x}-\tilde{\mathbf{x}})$
for conciseness. We find that

$\displaystyle\text{Cov}(F(\mathbf{x}),F(\mathbf{x}^{\prime}))$
$\displaystyle=k(\mathbf{x},\mathbf{x}^{\prime}),$
$\displaystyle\text{Cov}\left(\frac{\partial F(\mathbf{x})}{\partial
x_{i}},\frac{\partial F(\mathbf{x}^{\prime})}{\partial x_{i}}\right)$
$\displaystyle=\frac{\partial^{2}k(\mathbf{x},\mathbf{x}^{\prime})}{\partial
x_{i}\partial x_{i}^{\prime}},$ $\displaystyle\text{Cov}\left(\frac{\partial
F(\boldsymbol{\beta}(\mathbf{x},s))}{\partial x_{i}},\frac{\partial
F(\boldsymbol{\beta}(\mathbf{x}^{\prime},t))}{\partial x_{i}}\right)$
$\displaystyle=\frac{\partial^{2}k(\boldsymbol{\beta}(\mathbf{x},s),\boldsymbol{\beta}(\mathbf{x}^{\prime},t))}{\partial
x_{i}\partial x_{i}^{\prime}},$
$\displaystyle\text{Cov}\left(\int_{0}^{1}\frac{\partial
F(\boldsymbol{\beta}(\mathbf{x},s))}{\partial
x_{i}}ds,\int_{0}^{1}\frac{\partial
F(\boldsymbol{\beta}(\mathbf{x}^{\prime},t))}{\partial x_{i}}dt\right)$
$\displaystyle=\int_{0}^{1}\int_{0}^{1}\frac{\partial^{2}k(\boldsymbol{\beta}(\mathbf{x},s),\boldsymbol{\beta}(\mathbf{x}^{\prime},t))}{\partial
x_{i}\partial x_{i}^{\prime}}dsdt,$

and finally,

$\displaystyle\text{Cov}(attr_{i}(\mathbf{x}),attr_{i}(\mathbf{x}^{\prime}))$
$\displaystyle=\text{Cov}\left((x_{i}-\tilde{x}_{i})\int_{0}^{1}\frac{\partial
F(\boldsymbol{\beta}(\mathbf{x},s))}{\partial
x_{i}}ds,(x_{i}^{\prime}-\tilde{x}_{i})\int_{0}^{1}\frac{\partial
F(\boldsymbol{\beta}(\mathbf{x}^{\prime},t))}{\partial x_{i}}dt\right)$
$\displaystyle=(x_{i}-\tilde{x}_{i})(x_{i}^{\prime}-\tilde{x}_{i})\text{Cov}\left(\int_{0}^{1}\frac{\partial
F(\boldsymbol{\beta}(\mathbf{x},s))}{\partial
x_{i}}ds,\int_{0}^{1}\frac{\partial
F(\boldsymbol{\beta}(\mathbf{x}^{\prime},t))}{\partial x_{i}}dt\right)$
$\displaystyle=(x_{i}-\tilde{x}_{i})(x_{i}^{\prime}-\tilde{x}_{i})\int_{0}^{1}\int_{0}^{1}\frac{\partial^{2}k(\boldsymbol{\beta}(\mathbf{x},s),\boldsymbol{\beta}(\mathbf{x}^{\prime},t))}{\partial
x_{i}\partial x_{i}^{\prime}}dsdt.$

Again, to see that this expression is defined, we recall that $k\in C^{2}$,
and therefore $\partial^{2}k/\partial x_{i}\partial x_{i}^{\prime}\in C^{0}$,
and thus is integrable over the unit square.

To show completeness of the output GP, we need to know that almost-surely
every sample from a differentiable GP is differentiable. Let $f$ be any
differentiable sample from the GP. The attributions $attr_{i}(\mathbf{x}|f)$
satisfy completeness by definition, and thus
$\sum_{i}attr_{i}(\mathbf{x}|f)=f(\mathbf{x})-f(\tilde{\mathbf{x}})$. Since
this holds for almost-every $f$, the distributions of the left and right hand
sides coincide. We can see from either side that the distribution is Gaussian,
and we can compute that

$\mathbb{E}(F(\mathbf{x})-F(\tilde{\mathbf{x}}))=m(\mathbf{x})-m(\tilde{\mathbf{x}}),$

and

$\displaystyle\text{Cov}(F(\mathbf{x})-F(\tilde{\mathbf{x}}),F(\mathbf{x}^{\prime})-F(\tilde{\mathbf{x}}))$
$\displaystyle=\text{Cov}(F(\mathbf{x}),F(\mathbf{x}^{\prime}))+\text{Cov}(F(\tilde{\mathbf{x}}),F(\tilde{\mathbf{x}}))$
$\displaystyle\quad-\text{Cov}(F(\mathbf{x}),F(\tilde{\mathbf{x}}))-\text{Cov}(F(\tilde{\mathbf{x}}),F(\mathbf{x}^{\prime}))$
$\displaystyle=k(\mathbf{x},\mathbf{x}^{\prime})+k(\tilde{\mathbf{x}},\tilde{\mathbf{x}})$
$\displaystyle\quad-k(\mathbf{x},\tilde{\mathbf{x}})-k(\mathbf{x}^{\prime},\tilde{\mathbf{x}}).$

If $F\sim\mathcal{GP}(m,k)$ is the GPR posterior predictive distribution, then
$attr_{i}(\mathbf{x}|F)\sim\mathcal{GP}$ whenever $m\in C^{1}$ and $k\in
C^{2}$. If we satisfy the more refined assumption that $m\in C^{\infty}$ and
$k\in C^{\infty}$, then we have that $attr_{i}(\mathbf{x}|F)\in C^{\infty}$ as
well.

##  Appendix B Derivations for Specific Models

In this section, we derive the expressions for the mean and variance of the
attributions for specific models.

### ARD-SE Kernel

To derive the GPR attributions when using an ARD-SE kernel, it is instructive
to work with Corollary 1. From here, we observe that one must derive
expressions for $\mathbf{attr}_{i}(\mathbf{x}|\mathbf{k})$ and the double
integral in (30). In the main text, we defined

$A_{n,i}(\mathbf{x})\stackrel{{\scriptstyle\Delta}}{{=}}attr_{i}(\mathbf{x}|k(\cdot,\mathbf{x}_{n})),$

which are the entries of the vector
$\mathbf{attr}_{i}(\mathbf{x}|\mathbf{k})=[A_{n,i}]_{n=1:N}$. To compute
$A_{n,i}$, we must apply the $attr_{i}$ operator to the kernel function $k$,
which we know to be the ARD-SE function,

$k_{\text{ARD-
SE}}(\mathbf{x},\mathbf{x}^{\prime})=\sigma^{2}_{f}\exp\left(-\sum_{i=1}^{D}\frac{(x_{i}-x_{i}^{\prime})^{2}}{2\ell_{i}^{2}}\right).$

If we let $\mathbf{x}^{\prime}=\mathbf{x}_{n}$ and
$\mathbf{L}=\text{diag}([\ell_{1}\cdots\ell_{D}])$, then we can write out the
kernel function in matrix form:

$k(\mathbf{x},\mathbf{x}_{n})=\sigma^{2}_{f}\exp\left(-\frac{1}{2}(\mathbf{x}-\mathbf{x}_{n})^{\top}\mathbf{L}^{-2}(\mathbf{x}-\mathbf{x}_{n})\right).$

The partial derivative of this is

$\frac{\partial k(\mathbf{x},\mathbf{x}_{n})}{\partial
x_{i}}=\sigma^{2}_{f}\exp\left(-\frac{1}{2}(\mathbf{x}-\mathbf{x}_{n})^{\top}\mathbf{L}^{-2}(\mathbf{x}-\mathbf{x}_{n})\right)\left(-\frac{x_{i}-x_{n,i}}{\ell_{i}^{2}}\right).$

The expression begins to look hairy when we precompose with
$\boldsymbol{\gamma}(t)=\tilde{\mathbf{x}}+t(\mathbf{x}-\tilde{\mathbf{x}})$:

$\displaystyle\frac{\partial
k(\boldsymbol{\gamma}(t),\mathbf{x}_{n})}{\partial x_{i}}$
$\displaystyle=\sigma^{2}_{f}\exp{\left(-\frac{1}{2}(\tilde{\mathbf{x}}+t(\mathbf{x}-\tilde{\mathbf{x}})-\mathbf{x}_{n})^{\top}\mathbf{L}^{-2}(\tilde{\mathbf{x}}+t(\mathbf{x}-\tilde{\mathbf{x}})-\mathbf{x}_{n})\right)}$
$\displaystyle\times\left(-\frac{\tilde{x}_{i}+t(x_{i}-\tilde{x}_{i})-x_{n,i}}{\ell_{i}^{2}}\right).$

At this point, the most challenging task is to compute the integral from $t=0$
to $t=1$. The trick is to notice that the expression above is of the form

$\exp(-at^{2}-bt-c)(d\times t+f)$

as a function of $t$, where $a,b,c,d$ and $f$ are constants that we later plug
in. From here, one can produce closed-form expressions for the integral. It
can be checked that

$\int
t\exp(-at^{2}-bt-c)dt=-\frac{e^{-at^{2}-bt-c}}{2a}-\frac{\sqrt{\pi}e^{(b^{2}/4a)-c}\text{erf}\left(\frac{2at+b}{2\sqrt{a}}\right)}{4a^{3/2}},$

and

$\int\exp(-at^{2}-bt-c)dt=-\frac{\sqrt{\pi}e^{(b^{2}/4a)-c}\text{erf}\left(\frac{2at+b}{2\sqrt{a}}\right)}{2\sqrt{a}},$

where erf is the error function,

$\text{erf}(t)=\frac{2}{\sqrt{\pi}}\int_{0}^{z}e^{-t^{2}}dt.$

From the indefinite integral, we can compute the definite integral by
evaluating the expressions at $t=0$ and $t=1$. Based on the expression above,
the appropriate definitions for the constants $a,b,c,d$ and $f$ are

$\displaystyle a$
$\displaystyle=(\mathbf{x}-\tilde{\mathbf{x}})^{\top}\mathbf{L}^{-2}(\mathbf{x}-\tilde{\mathbf{x}}),$
$\displaystyle b$
$\displaystyle=2(\mathbf{x}-\tilde{\mathbf{x}})^{\top}\mathbf{L}^{-2}(\tilde{\mathbf{x}}-\mathbf{x}_{n}),$
$\displaystyle c$
$\displaystyle=(\tilde{\mathbf{x}}-\mathbf{x}_{n})^{\top}\mathbf{L}^{-2}(\tilde{\mathbf{x}}-\mathbf{x}_{n}),$
$\displaystyle d$
$\displaystyle=-\sigma^{2}_{0}\frac{(x_{i}-\tilde{x}_{i})^{2}}{\ell_{i}^{2}},$
$\displaystyle f$
$\displaystyle=-\sigma^{2}_{0}\frac{(x_{i}-\tilde{x}_{i})(\tilde{x}_{i}-x_{n,i})}{\ell_{i}^{2}}.$

After some algebra we deduce that

$A_{n,i}(\mathbf{x})=\frac{e^{\frac{-a-b-c}{2}}\left[\sqrt{2\pi}e^{\frac{(2a+b)^{2}}{8a}}\left(\text{erf}\left(\frac{b}{2\sqrt{2a}}\right)-\text{erf}\left(\frac{2a+b}{2\sqrt{2a}}\right)\right)(bd-2af)+4\sqrt{a}d\left(e^{\frac{a+b}{2}}-1\right)\right]}{4a^{3/2}},$

and thus

$\mathbb{E}(attr_{i}(\mathbf{x}))=\mathbf{attr}_{i}(\mathbf{x}|\mathbf{k})^{\top}\boldsymbol{\alpha}=\sum_{n=1}^{N}\alpha_{n}A_{n,i}(\mathbf{x}).$

Computation of the variance is also tractable, but requires some more work.
Let $\mathbf{\Lambda}=(\mathbf{K}+\sigma^{2}_{n}\mathbf{I})^{-1}$. Then we
observe that

$\text{Var}(attr_{i}(\mathbf{x}))=B_{i}(\mathbf{x})-\mathbf{attr}_{i}(\mathbf{x}|\mathbf{k})^{\top}\mathbf{\Lambda}\mathbf{attr}_{i}(\mathbf{x}|\mathbf{k}),$

where we have defined

$B_{i}(\mathbf{x})\stackrel{{\scriptstyle\Delta}}{{=}}(x_{i}-\tilde{x}_{i})^{2}\int_{0}^{1}\int_{0}^{1}\frac{\partial^{2}k(\boldsymbol{\gamma}(s),\boldsymbol{\gamma}(t))}{\partial
x_{i}\partial x_{i}^{\prime}}dsdt.$

We already know $\mathbf{attr}_{i}(\mathbf{x}|\mathbf{k})$ from computing the
mean attribution, and all that remains is to compute $B_{i}(\mathbf{x})$. Like
before, our goal is to convert $B_{i}(\mathbf{x})$ into a generic expression
and then use a known result to solve the integral. Firstly, we compute the
second derivative, which we may express as

$\displaystyle\frac{\partial^{2}k(\mathbf{x},\mathbf{x}^{\prime})}{\partial
x_{i}\partial x_{i}^{\prime}}$
$\displaystyle=\sigma^{2}_{0}\frac{\partial^{2}}{\partial x_{i}\partial
x_{i}^{\prime}}\exp\left(-\frac{1}{2}\sum_{i=1}^{D}\frac{(x_{i}-x_{i}^{\prime})^{2}}{\ell_{i}^{2}}\right)$
$\displaystyle=\sigma^{2}_{0}\exp\left(-\frac{1}{2}\sum_{i=1}^{D}\frac{(x_{i}-x_{i}^{\prime})^{2}}{\ell_{i}^{2}}\right)\left(\frac{1}{\ell_{i}^{2}}-\frac{(x_{i}-x_{i}^{\prime})^{2}}{\ell_{i}^{4}}\right)$
$\displaystyle=\sigma^{2}_{0}\exp\left(-\frac{(\mathbf{x}-\mathbf{x}^{\prime})^{\top}\mathbf{L}^{-2}(\mathbf{x}-\mathbf{x}^{\prime})}{2}\right)\left(\frac{1}{\ell_{i}^{2}}-\frac{(x_{i}-x_{i}^{\prime})^{2}}{\ell_{i}^{4}}\right).$

The next step is to precompose and integrate. When
$\mathbf{x}\neq\mathbf{x}^{\prime}$, the general integral is complicated. The
formulas simplify considerably if we constrain our attention to the case in
which $\mathbf{x}=\mathbf{x}^{\prime}$, meaning that we only compute the
variance of an attribution, and not the full covariance function. In this
case, we let
$\boldsymbol{\gamma}(t)=\tilde{\mathbf{x}}+t(\mathbf{x}-\tilde{\mathbf{x}})$,
and then precomposition yields

$\displaystyle\frac{\partial^{2}k(\boldsymbol{\gamma}(s),\boldsymbol{\gamma}(t))}{\partial
x_{i}\partial x_{i}^{\prime}}$
$\displaystyle=\sigma^{2}_{0}\exp\left(-\frac{(\boldsymbol{\gamma}(s)-\boldsymbol{\gamma}(t))^{\top}\mathbf{L}^{-2}(\boldsymbol{\gamma}(s)-\boldsymbol{\gamma}(t))}{2}\right)\left(\frac{1}{\ell_{i}^{2}}-\frac{(\gamma_{i}(s)-\gamma_{i}(t))^{2}}{\ell_{i}^{4}}\right)$
$\displaystyle=\sigma^{2}_{0}\exp\left(-\frac{(s-t)^{2}(\mathbf{x}-\tilde{\mathbf{x}})^{\top}\mathbf{L}^{-2}(\mathbf{x}-\tilde{\mathbf{x}})}{2}\right)\left(\frac{1}{\ell_{i}^{2}}-\frac{(s-t)^{2}(x_{i}-\tilde{x}_{i})^{2}}{\ell_{i}^{4}}\right)$
$\displaystyle=\exp\left(-\frac{a(s-t)^{2}}{2}\right)\left(w+v(s-t)^{2}\right),$

where
$a=(\mathbf{x}-\tilde{\mathbf{x}})^{\top}\mathbf{L}^{-2}(\mathbf{x}-\tilde{\mathbf{x}})$
is the same as in the calculation for the mean, and where we additionally
defined

$\displaystyle v$
$\displaystyle=-\sigma^{2}_{0}\frac{(x_{i}-\tilde{x}_{i})^{2}}{\ell_{i}^{4}},$
$\displaystyle w$ $\displaystyle=\frac{\sigma^{2}_{0}}{\ell_{i}^{2}}.$

In this form, the integral can be evaluated, yielding

$\int_{0}^{1}\int_{0}^{1}e^{-a(s-t)^{2}/2}(w+v(s-t)^{2})dsdt=\frac{\sqrt{2\pi}\text{erf}\left(\frac{\sqrt{a}}{\sqrt{2}}\right)(aw+v)}{a^{3/2}}+\frac{2e^{-\frac{a}{2}}\left(e^{a/2}-1\right)(aw+2v)}{a^{2}}.$

We may then write

$\displaystyle
B_{i}(\mathbf{x})=(x_{i}-\tilde{x}_{i})^{2}\left(\frac{\sqrt{2\pi}\text{erf}\left(\frac{\sqrt{a}}{\sqrt{2}}\right)(aw+v)}{a^{3/2}}-\frac{2e^{-\frac{a}{2}}\left(e^{a/2}-1\right)(aw+2v)}{a^{2}}\right).$

### Random Feature GPs

To derive the RFGP attributions, we can use Theorem 1 directly with (21) and
(22). From (21), we see that the RFGP mean can be expressed as a weighted sum.
If we define the vector
$\boldsymbol{\alpha}\stackrel{{\scriptstyle\Delta}}{{=}}\mathbf{A}^{-1}\boldsymbol{\Phi}\mathbf{y}$,
then we may write

$\boldsymbol{\phi}(\mathbf{x})^{\top}\mathbf{A}^{-1}\boldsymbol{\Phi}\mathbf{y}=\boldsymbol{\phi}(\mathbf{x})^{\top}\boldsymbol{\alpha}=\sum_{m=1}^{2M}\phi_{m}(\mathbf{x})\alpha_{m},$

where $\alpha_{m}$ are constants and

$\phi_{m}(\mathbf{x})=\begin{cases}\sin(\mathbf{x}^{\top}\mathbf{v}_{(m+1)/2}),&\text{if
$m$ odd}\\\
\cos(\mathbf{x}^{\top}\mathbf{v}_{m/2}),&\text{otherwise}\end{cases}.$

In either case, $\phi_{m}(\mathbf{x})$ is of the form
$\sin(\mathbf{x}^{\top}\mathbf{v})$ or $\cos(\mathbf{x}^{\top}\mathbf{v})$,
for some value of $\mathbf{v}$. By linearity, it suffices to compute the
attribution for this generic form, and then plug in to find the resulting
expression.

Since $\cos(x)$ can be represented as $\sin(x+b)$ where $b=\pi/2$, we show the
generic calculation for $\sin(\mathbf{x}^{\top}\mathbf{v}+b)$ that will
produce the desired result for both the sine and cosine cases. Let
$\phi(\mathbf{x})=\sin(\mathbf{x}^{\top}\mathbf{v}+b)$. Then

$\displaystyle attr_{i}(\mathbf{x}|\phi)$
$\displaystyle=(x_{i}-\tilde{x}_{i})\int_{0}^{1}\frac{\partial\phi(\tilde{\mathbf{x}}+t(\mathbf{x}-\tilde{\mathbf{x}}))}{\partial
x_{i}}dt$
$\displaystyle=(x_{i}-\tilde{x}_{i})\int_{0}^{1}\cos((\tilde{\mathbf{x}}+t(\mathbf{x}-\tilde{\mathbf{x}}))^{\top}\mathbf{v}+b)v_{i}dt$
$\displaystyle=(x_{i}-\tilde{x}_{i})\left[\sin((\tilde{\mathbf{x}}+t(\mathbf{x}-\tilde{\mathbf{x}}))^{\top}\mathbf{v}+b)\right|_{t=0}^{1}\frac{v_{i}}{(\mathbf{x}-\tilde{\mathbf{x}})^{\top}\mathbf{v}}$
$\displaystyle=(x_{i}-\tilde{x}_{i})\left(\sin(\mathbf{x}^{\top}\mathbf{v}+b)-\sin(\tilde{\mathbf{x}}^{\top}\mathbf{v}+b)\right)\frac{v_{i}}{(\mathbf{x}-\tilde{\mathbf{x}})^{\top}\mathbf{v}}.$

From this calculation, we can derive the expression for the vector of $\phi$
attributions:

$\mathbf{attr}_{i}(\mathbf{x}|\boldsymbol{\phi})=\begin{bmatrix}attr_{i}(\mathbf{x}|\phi_{1})\\\
\vdots\\\
attr_{i}(\mathbf{x}|\phi_{2M})\end{bmatrix}=(x_{i}-\tilde{x}_{i})\begin{bmatrix}\frac{v_{1,i}}{\mathbf{v}_{1}^{\top}(\mathbf{x}-\tilde{\mathbf{x}})}(\cos(\mathbf{v}_{1}^{\top}\mathbf{x})-\cos(\mathbf{v}_{1}^{\top}\tilde{\mathbf{x}}))\\\
\frac{v_{1,i}}{\mathbf{v}_{1}^{\top}(\mathbf{x}-\tilde{\mathbf{x}})}(\sin(\mathbf{v}_{1}^{\top}\mathbf{x})-\sin(\mathbf{v}_{1}^{\top}\tilde{\mathbf{x}}))\\\
\vdots\\\
\frac{v_{M,i}}{\mathbf{v}_{M}^{\top}(\mathbf{x}-\tilde{\mathbf{x}})}\left(\cos(\mathbf{v}_{M}^{\top}\mathbf{x})-\cos(\mathbf{v}_{M}^{\top}\tilde{\mathbf{x}})\right)\\\
\frac{v_{M,i}}{\mathbf{v}_{M}^{\top}(\mathbf{x}-\tilde{\mathbf{x}})}\left(\sin(\mathbf{v}_{M}^{\top}\mathbf{x})-\sin(\mathbf{v}_{M}^{\top}\tilde{\mathbf{x}})\right)\end{bmatrix}$

and thus

$\mathbb{E}(attr_{i}(\mathbf{x}))=\mathbf{attr}_{i}(\mathbf{x}|\boldsymbol{\phi})^{\top}\boldsymbol{\alpha}.$

In the main text, we use the notation

$\boldsymbol{\zeta}(\mathbf{x})=\begin{bmatrix}\frac{v_{1,i}}{\mathbf{v}_{1}^{\top}(\mathbf{x}-\tilde{\mathbf{x}})}(\cos(\mathbf{v}_{1}^{\top}\mathbf{x})-\cos(\mathbf{v}_{1}^{\top}\tilde{\mathbf{x}}))\\\
\frac{v_{1,i}}{\mathbf{v}_{1}^{\top}(\mathbf{x}-\tilde{\mathbf{x}})}(\sin(\mathbf{v}_{1}^{\top}\mathbf{x})-\sin(\mathbf{v}_{1}^{\top}\tilde{\mathbf{x}}))\\\
\vdots\\\
\frac{v_{M,i}}{\mathbf{v}_{M}^{\top}(\mathbf{x}-\tilde{\mathbf{x}})}\left(\cos(\mathbf{v}_{M}^{\top}\mathbf{x})-\cos(\mathbf{v}_{M}^{\top}\tilde{\mathbf{x}})\right)\\\
\frac{v_{M,i}}{\mathbf{v}_{M}^{\top}(\mathbf{x}-\tilde{\mathbf{x}})}\left(\sin(\mathbf{v}_{M}^{\top}\mathbf{x})-\sin(\mathbf{v}_{M}^{\top}\tilde{\mathbf{x}})\right)\end{bmatrix}$

so that we may compactly write

$\mathbb{E}(attr_{i}(\mathbf{x}))=(x_{i}-\tilde{x}_{i})\boldsymbol{\zeta}^{\top}\boldsymbol{\alpha}=(x_{i}-\tilde{x}_{i})\boldsymbol{\zeta}^{\top}\mathbf{A}^{-1}\boldsymbol{\Phi}\mathbf{y}.$

This notation is also useful for deriving the covariance, since when we write

$\text{Cov}(attr_{i}(\mathbf{x}),attr_{i}(\mathbf{x}^{\prime}))=(x_{i}-\tilde{x}_{i})(x_{i}^{\prime}-\tilde{x}_{i})\int_{0}^{1}\int_{0}^{1}\frac{\partial^{2}k(\tilde{\mathbf{x}}+s(\mathbf{x}-\tilde{\mathbf{x}}),\tilde{\mathbf{x}}+t(\mathbf{x}^{\prime}-\tilde{\mathbf{x}}))}{\partial
x_{i}\partial x_{i}^{\prime}}dtds,$

for
$k(\mathbf{x},\mathbf{x}^{\prime})=\sigma^{2}_{n}\boldsymbol{\phi}(\mathbf{x})^{\top}\mathbf{A}^{-1}\boldsymbol{\phi}(\mathbf{x}^{\prime})$,
we can split up the derivatives and integrals. Without writing everything out,
the desired result is

$\text{Cov}(attr_{i}(\mathbf{x}),attr_{i}(\mathbf{x}^{\prime}))=\sigma^{2}_{n}\mathbf{attr}_{i}(\mathbf{x}|\boldsymbol{\phi})^{\top}\mathbf{A}^{-1}\mathbf{attr}_{i}(\mathbf{x}^{\prime}|\boldsymbol{\phi}),$

which immediately yields

$\text{Cov}(attr_{i}(\mathbf{x}),attr_{i}(\mathbf{x}^{\prime}))=\sigma^{2}_{n}(x_{i}-\tilde{x}_{i})(x_{i}^{\prime}-\tilde{x}_{i})\boldsymbol{\zeta}(\mathbf{x})^{\top}\mathbf{A}^{-1}\boldsymbol{\zeta}(\mathbf{x}^{\prime}).$

##  Appendix C Additional Experiments

In this appendix, we include some additional experiments and validation that
we omit from the main text.

###  C.1 Wine: Quadrature Rules

In this section, we provide some additional figures to analyze the convergence
of the approximated integral. First, we visualize the convergence of the
approximated integral to the true IG attributions, for the Wine data set in
Figure 6, as only one feature was shown in Figure 4(a). Figure 6 shows that
the rate of convergence is similar across all input features, and does not
obviously correspond to the magnitude of the exact attribution. The ARD-SE
kernel also provides a measure of feature importance, where relevance of a
(normalized) feature is defined to be $r_{i}=1/\ell_{i}$. In this case, we
observe that for features with relevance, their attributions converge much
more quickly. This is expected since a low feature relevance indicates that
model training has preferentially underused a given feature in making
predictions.

Figure 6: Convergence of the attributions with approximated integration to the
exact attributions, in the Wine data set. The first 11 subfigures shown the
approximation convergence, as described in Figure 4(a). The last two
subfigures show the relative sizes of each attribution in the exact result,
and also the relative importance of each input feature (as measured by the ARD
kernel), respectively.

Second, we visualize the impact of different quadrature rules. The formula in
(36) uses the right-hand rule to approximate the integral of a function on the
interval $[0,1]$, but there exist other quadrature rules, or methods for
numerical integration, that may achieve noticeable improvements in performance
with only changes to the weights in the sum used to approximate the integral
Heath (2018). In Figure 7, we compare the application of the right-hand rule
with two other quadrature rules: the trapezoid rule and Simpson’s rule. Each
method approximates the integral using a partition of the unit interval:

$0=t_{0}<t_{1}<\cdots<t_{L-1}<t_{L}=1.$

For each $l=1,...,L$, each quadrature approximates a portion of the integral
along that segment. For a general function $g(t)$, the three quadratures use
the partition to approximate the integral differently:

$\int_{0}^{1}g(t)dt\approx\sum_{l=1}^{L}S_{l}$ Right-hand: $\displaystyle
S_{l}=g(t_{l}),$ Trapezoid: $\displaystyle
S_{l}=\frac{1}{2}\left(g(t_{l-1})+g(t_{l})\right),$ Simpson’s rule:
$\displaystyle
S_{l}=\frac{1}{4}\left(g(t_{l-1})+2g\left(\frac{t_{l-1}+t_{l}}{2}\right)+g(t_{l})\right).$

We note that Simpson’s rule includes an extra evaluation per partition. As a
result, Simpson’s rule should be compared to other methods relative to the
number of function calls, instead of $L$. In Figure 7, we compare the
convergence of each quadrature rule against the number of partitions $L$ as
well as the number of required function evaluations. We additionally visualize
the mean-squared error (MSE) of the means of each attribution GP. Both the
trapezoid rule and Simpson’s rule demonstrate a marked improvement over the
right-hand rule.

((a)) Attribution GPs vs the number of partitions $L$.

((b)) Attribution GPs vs the number of function evaluations.

((c)) MSE when estimating the mean attribution.

Figure 7: Comparison of the three quadrature rules applied to compute the
attributions of a GPR model with ARD-SE kernel. Each method was tasked with
computing the attribution to Alcohol, as seen in Figure 6.

## References

  * mis (2018) Real estate valuation data set.  UCI Machine Learning Repository, 2018.  DOI: https://doi.org/10.24432/C5J30W.
  * Adebayo et al. (2018) Julius Adebayo, Justin Gilmer, Michael Muelly, Ian Goodfellow, Moritz Hardt, and Been Kim.  Sanity checks for saliency maps.  _Advances in neural information processing systems_ , 31, 2018.
  * Adler and Taylor (2007) Robert J Adler and Jonathan E Taylor.  _Random Fields and Geometry_.  Springer, 2007.
  * Angelov et al. (2021) Plamen P Angelov, Eduardo A Soares, Richard Jiang, Nicholas I Arnold, and Peter M Atkinson.  Explainable artificial intelligence: an analytical review.  _Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery_ , 11(5):e1424, 2021.
  * Arrieta et al. (2020) Alejandro Barredo Arrieta, Natalia Díaz-Rodríguez, Javier Del Ser, Adrien Bennetot, Siham Tabik, Alberto Barbado, Salvador García, Sergio Gil-López, Daniel Molina, Richard Benjamins, et al.  Explainable artificial intelligence (XAI): Concepts, taxonomies, opportunities and challenges toward responsible AI.  _Information Fusion_ , 58:82–115, 2020.
  * Aumann and Shapley (2015) Robert J Aumann and Lloyd S Shapley.  _Values of Non-Atomic Games_.  Princeton University Press, 2015.
  * Bach et al. (2015) Sebastian Bach, Alexander Binder, Grégoire Montavon, Frederick Klauschen, Klaus-Robert Müller, and Wojciech Samek.  On pixel-wise explanations for non-linear classifier decisions by layer-wise relevance propagation.  _PLOS One_ , 10(7):e0130140, 2015.
  * Baehrens et al. (2010) David Baehrens, Timon Schroeter, Stefan Harmeling, Motoaki Kawanabe, Katja Hansen, and Klaus-Robert Müller.  How to explain individual classification decisions.  _The Journal of Machine Learning Research_ , 11:1803–1831, 2010.
  * Carlini et al. (2023) Nicolas Carlini, Jamie Hayes, Milad Nasr, Matthew Jagielski, Vikash Sehwag, Florian Tramer, Borja Balle, Daphne Ippolito, and Eric Wallace.  Extracting training data from diffusion models.  In _32nd USENIX Security Symposium (USENIX Security 23)_ , pages 5253–5270, 2023.
  * Chau et al. (2023) Siu Lun Chau, Krikamol Muandet, and Dino Sejdinovic.  Explaining the uncertain: Stochastic Shapley values for Gaussian process models.  _arXiv preprint arXiv:2305.15167_ , 2023.
  * Cortez et al. (2009a) Paolo Cortez, Antonio Luíz Cerdeira, Fernando Almeida, Telmo Matos, and José Reis.  Modeling wine preferences by data mining from physicochemical properties.  _Decis. Support Syst._ , 47:547–553, 2009a.  URL https://api.semanticscholar.org/CorpusID:2996254.
  * Cortez et al. (2009b) Paulo Cortez, A. Cerdeira, F. Almeida, T. Matos, and J. Reis.  Wine Quality.  UCI Machine Learning Repository, 2009b.  DOI: https://doi.org/10.24432/C56S3T.
  * Damianou and Lawrence (2013) Andreas Damianou and Neil D Lawrence.  Deep Gaussian processes.  In _Artificial Intelligence and Statistics_ , pages 207–215. PMLR, 2013.
  * de Bruijn et al. (2022) Hans de Bruijn, Martijn Warnier, and Marijn Janssen.  The perils and pitfalls of explainable AI: Strategies for explaining algorithmic decision-making.  _Government Information Quarterly_ , 39(2):101666, 2022.
  * Duvenaud et al. (2011) David K Duvenaud, Hannes Nickisch, and Carl Rasmussen.  Additive Gaussian processes.  _Advances in Neural Information Processing Systems_ , 24, 2011.
  * Friedman and Moulin (1999) Eric Friedman and Herve Moulin.  Three methods to share joint costs or surplus.  _Journal of Economic Theory_ , 87(2):275–312, 1999.
  * Friedman (2004) Eric J Friedman.  Paths and consistency in additive cost sharing.  _International Journal of Game Theory_ , 32:501–518, 2004\.
  * Ghassemi et al. (2021) Marzyeh Ghassemi, Luke Oakden-Rayner, and Andrew L Beam.  The false hope of current approaches to explainable artificial intelligence in health care.  _The Lancet Digital Health_ , 3(11):e745–e750, 2021.
  * Gibbs (1998) Mark N Gibbs.  _Bayesian Gaussian processes for regression and classification_.  PhD thesis, Citeseer, 1998.
  * Glorot et al. (2011) Xavier Glorot, Antoine Bordes, and Yoshua Bengio.  Deep sparse rectifier neural networks.  In _Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics_ , pages 315–323. JMLR Workshop and Conference Proceedings, 2011.  URL http://proceedings.mlr.press/v15/glorot11a/glorot11a.pdf.
  * Grömping (2007) Ulrike Grömping.  Estimators of relative importance in linear regression based on variance decomposition.  _The American Statistician_ , 61(2):139–147, 2007\.
  * Heath (2018) Michael T Heath.  _Scientific Computing: An Introductory Survey_.  SIAM, second edition, 2018.
  * Hensman et al. (2013) James Hensman, Nicolo Fusi, and Neil D Lawrence.  Gaussian processes for big data.  In _Uncertainty in Artificial Intelligence_ , page 282, 2013.
  * Hess and Ditzler (2023) Samuel Hess and Gregory Ditzler.  Protoshotxai: Using prototypical few-shot architecture for explainable AI.  _Journal of Machine Learning Research_ , 24(325):1–49, 2023.  URL http://jmlr.org/papers/v24/21-1261.html.
  * Huang and Wang (2022) Bin Huang and Jianhui Wang.  Applications of physics-informed neural networks in power systems-a review.  _IEEE Transactions on Power Systems_ , 38(1):572–588, 2022.
  * Jin et al. (2022) Di Jin, Elena Sergeeva, Wei-Hung Weng, Geeticka Chauhan, and Peter Szolovits.  Explainable deep learning in healthcare: A methodological survey from an attribution view.  _WIREs Mechanisms of Disease_ , 14(3):e1548, 2022\.
  * Krizhevsky et al. (2012) Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.  ImageNet classification with deep convolutional neural networks.  In _Advances in Neural Information Processing Systems_ , volume 25, 2012.
  * Lázaro-Gredilla et al. (2010) Miguel Lázaro-Gredilla, Joaquin Quinonero-Candela, Carl Edward Rasmussen, and Aníbal R Figueiras-Vidal.  Sparse spectrum Gaussian process regression.  _The Journal of Machine Learning Research_ , 11:1865–1881, 2010.
  * Le Gall (2016) Jean-François Le Gall.  _Brownian Motion, Martingales, and Stochastic Calculus_ , volume 274 of _Graduate Texts in Mathematics_.  Springer, 2016.
  * Lee (2012) John M Lee.  _Introduction to Smooth Manifolds_ , volume 218 of _Graduate Texts in Mathematics_.  Springer, 2012.
  * Lundberg and Lee (2017) Scott M Lundberg and Su-In Lee.  A unified approach to interpreting model predictions.  In _Advances in Neural Information Processing Systems_ , volume 30, 2017.  URL https://proceedings.neurips.cc/paper_files/paper/2017/file/8a20a8621978632d76c43dfd28b67767-Paper.pdf.
  * Lundstrom et al. (2022) Daniel D Lundstrom, Tianjian Huang, and Meisam Razaviyayn.  A rigorous study of integrated gradients method and extensions to internal neuron attributions.  In _International Conference on Machine Learning_ , pages 14485–14508. PMLR, 2022.
  * Mitchell et al. (2022) Rory Mitchell, Joshua Cooper, Eibe Frank, and Geoffrey Holmes.  Sampling permutations for Shapley value estimation.  _The Journal of Machine Learning Research_ , 23(1):2082–2127, 2022.
  * Murphy (2022) Kevin P. Murphy.  _Probabilistic Machine Learning: An Introduction_.  MIT Press, 2022.  URL probml.ai.
  * Neal (1996) Radford M Neal.  _Bayesian Learning for Neural Networks_ , volume 118 of _Lecture Notes in Statistics_.  Springer, 1996.
  * Peng et al. (2022) Xi Peng, Yunfan Li, Ivor W Tsang, Hongyuan Zhu, Jiancheng Lv, and Joey Tianyi Zhou.  XAI beyond classification: Interpretable neural clustering.  _The Journal of Machine Learning Research_ , 23(1):227–254, 2022.
  * Rasmussen and Williams (2006) Carl E. Rasmussen and Christopher K. Williams.  _Gaussian Processes for Machine Learning_.  MIT Press, 2006.  URL https://gaussianprocess.org/gpml/.
  * Ribeiro et al. (2016) Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin.  Model-agnostic interpretability of machine learning.  _arXiv preprint arXiv:1606.05386_ , 2016.
  * Riutort-Mayol et al. (2023) Gabriel Riutort-Mayol, Paul-Christian Bürkner, Michael R Andersen, Arno Solin, and Aki Vehtari.  Practical Hilbert space approximate Bayesian Gaussian processes for probabilistic programming.  _Statistics and Computing_ , 33(1):17, 2023.
  * Robnik-Šikonja and Kononenko (2008) Marko Robnik-Šikonja and Igor Kononenko.  Explaining classifications for individual instances.  _IEEE Transactions on Knowledge and Data Engineering_ , 20(5):589–600, 2008.
  * Rudin (2019) Cynthia Rudin.  Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead.  _Nature Machine Intelligence_ , 1(5):206–215, 2019.
  * Samek et al. (2021) Wojciech Samek, Grégoire Montavon, Sebastian Lapuschkin, Christopher J Anders, and Klaus-Robert Müller.  Explaining deep neural networks and beyond: A review of methods and applications.  _Proceedings of the IEEE_ , 109(3):247–278, 2021\.
  * Seitz (2022) Sarem Seitz.  Gradient-based explanations for Gaussian process regression and classification models.  _arXiv preprint arXiv:2205.12797_ , 2022.
  * Shapley (1953) Lloyd S Shapley.  A value for n-person games.  _Contributions to the Theory of Games_ , pages 307–317, 1953.
  * Silva et al. (2020) Andrew Silva, Matthew Gombolay, Taylor Killian, Ivan Jimenez, and Sung-Hyun Son.  Optimization methods for interpretable differentiable decision trees applied to reinforcement learning.  In _International conference on artificial intelligence and statistics_ , pages 1855–1865. PMLR, 2020.
  * Simonyan et al. (2014) Karen Simonyan, Andrea Vedaldi, and Andrew Zisserman.  Deep inside convolutional networks: visualising image classification models and saliency maps.  In _Proceedings of the International Conference on Learning Representations (ICLR)_. ICLR, 2014.
  * Solak et al. (2002) Ercan Solak, Roderick Murray-Smith, WE Leithead, D Leith, and Carl Rasmussen.  Derivative observations in Gaussian process models of dynamic systems.  _Advances in Neural Information Processing Systems_ , 15, 2002.
  * Springenberg et al. (2014) Jost Tobias Springenberg, Alexey Dosovitskiy, Thomas Brox, and Martin Riedmiller.  Striving for simplicity: The all convolutional net.  _arXiv preprint arXiv:1412.6806_ , 2014.
  * Stein (1999) Michael L Stein.  _Interpolation of Spatial Data: Some Theory for Kriging_.  Springer Science & Business Media, 1999.
  * Štrumbelj and Kononenko (2014) Erik Štrumbelj and Igor Kononenko.  Explaining prediction models and individual predictions with feature contributions.  _Knowledge and Information Systems_ , 41(3):647–665, 2014.
  * Štrumbelj et al. (2009) Erik Štrumbelj, Igor Kononenko, and M Robnik Šikonja.  Explaining instance classifications with interactions of subsets of feature values.  _Data & Knowledge Engineering_, 68(10):886–904, 2009.
  * Sturmfels et al. (2020) Pascal Sturmfels, Scott Lundberg, and Su-In Lee.  Visualizing the impact of feature attribution baselines.  _Distill_ , 5(1):e22, 2020.  URL https://distill.pub/2020/attribution-baselines/.
  * Sundararajan and Najmi (2020) Mukund Sundararajan and Amir Najmi.  The many Shapley values for model explanation.  In _International Conference on Machine Learning_ , pages 9269–9278. PMLR, 2020.
  * Sundararajan et al. (2017) Mukund Sundararajan, Ankur Taly, and Qiqi Yan.  Axiomatic attribution for deep networks.  In _International Conference on Machine Learning_ , pages 3319–3328. PMLR, 2017.
  * van der Waa et al. (2021) Jasper van der Waa, Elisabeth Nieuwburg, Anita Cremers, and Mark Neerincx.  Evaluating XAI: A comparison of rule-based and example-based explanations.  _Artificial Intelligence_ , 291:103404, 2021.
  * Voita et al. (2019) Elena Voita, David Talbot, Fedor Moiseev, Rico Sennrich, and Ivan Titov.  Analyzing multi-head self-attention: Specialized heads do the heavy lifting, the rest can be pruned.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_. Association for Computational Linguistics, 2019\.
  * Wahba (1990) Grace Wahba.  _Spline Models for Observational Data_.  SIAM, 1990.
  * Wilson et al. (2016) Andrew Gordon Wilson, Zhiting Hu, Ruslan Salakhutdinov, and Eric P Xing.  Deep kernel learning.  In _Artificial Intelligence and Statistics_ , pages 370–378. PMLR, 2016.
  * Wolberg et al. (1995) William Wolberg, W. Street, and Olvi Mangasarian.  Breast Cancer Wisconsin (Prognostic).  UCI Machine Learning Repository, 1995.  DOI: https://doi.org/10.24432/C5GK50.
  * Ying et al. (2019) Zhitao Ying, Dylan Bourgeois, Jiaxuan You, Marinka Zitnik, and Jure Leskovec.  GNNExplainer: Generating explanations for graph neural networks.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Yoshikawa and Iwata (2021) Yuya Yoshikawa and Tomoharu Iwata.  Gaussian process regression with interpretable sample-wise feature weights.  _IEEE Transactions on Neural Networks and Learning Systems_ , 2021\.
  * Zhang et al. (2018) Hao Zhang, Shinji Nakadai, and Kenji Fukumizu.  From black-box to white-box: Interpretable learning with kernel machines.  In _Machine Learning and Data Mining in Pattern Recognition: 14th International Conference, MLDM 2018, New York, NY, USA, July 15-19, 2018, Proceedings, Part I 14_ , pages 213–227. Springer, 2018.