# Supplementary Material for
GPEX, A Framework For Interpreting Artificial Neural Networks

Amir Akbarnejad, Gilbert Bigras, and Nilanjan Ray

##  S1 Deriving the Variational Lower-bound

In this section we derive the variational lower-bound introduced in Sec.3.3 of
the main article. We firstly introduce Lemmas 1 and 2 as they appear in our
derivations.

###### Lemma 1

The KL-divergence between two normal distributions
$\mathcal{N}_{1}(.\;;\;\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{1})$ and
$\mathcal{N}_{2}(.\;;\;\boldsymbol{\mu}_{2},\boldsymbol{\Sigma}_{2})$ can be
computed as follows:

$\begin{split}KL\Big{(}\mathcal{N}_{1}\;||\;\mathcal{N}_{2}\Big{)}=\frac{1}{2}&\Big{(}\log(\frac{|\boldsymbol{\Sigma}_{2}|}{|\boldsymbol{\Sigma}_{1}|})\;-\;D\;+\;trace\\{\boldsymbol{\Sigma}_{2}^{-1}\boldsymbol{\Sigma}_{1}\\}\;+\;\\\
&(\boldsymbol{\mu}_{2}-\boldsymbol{\mu}_{1})^{T}\boldsymbol{\Sigma}_{2}^{-1}(\boldsymbol{\mu}_{2}-\boldsymbol{\mu}_{1})\Big{)}.\blacksquare\end{split}$
(S1)

###### Lemma 2

Let $p_{1}$ and $p_{2}$ be two normal distributions:

$\begin{split}p_{1}(x)=\mathcal{N}\big{(}x\;;\;\mu_{1},\sigma^{2}_{1}\big{)},\\\
p_{2}(x)=\mathcal{N}\big{(}x\;;\;\mu_{2},\sigma^{2}_{2}\big{)}.\end{split}$

We have that

$\begin{split}\mathbb{E}_{x\sim
p_{2}}\big{[}\log\;p_{1}(x\;;\;\mu_{1},\sigma^{2}_{1})\big{]}&=\\\
-\frac{(\mu_{1}-\mu_{2})^{2}+\sigma_{2}^{2}}{2\sigma_{1}^{2}}-&\frac{1}{2}\log(\sigma_{1}^{2})-\frac{1}{2}\log(2\pi).\blacksquare\end{split}$
(S2)

Figure S1: The proposed framework as a probabilistic graphical model.

Fig.S1 illustrates the framework as a probabilistic graphical model. A general
feed-forward pipeline takes in a set of input(s) $\mathcal{X}$ and produces a
set of output(s) $\mathcal{Y}$. The general pipeline is required to have at
least one ANN as a submodule. The ANN submodule is required to take in only
one input $\boldsymbol{x}$ and to produce only one output $\boldsymbol{v}$,
where $\boldsymbol{x}$ and $\boldsymbol{v}$ are tensors of arbitrary sizes. As
illustrated in Fig.S1, the ANN’s input $\boldsymbol{x}$ can depend arbitrarily
on some other intermediate variables in the pipeline. This relation is modeled
by the conditional distribution
$p\big{(}\boldsymbol{x}_{n}|Parent(\boldsymbol{x}_{n})\big{)}$ where
$Parent(\boldsymbol{x}_{n})$ is the set of all variables which are connected
to $\boldsymbol{x}_{n}$. Similarly, as illustrated in Fig.S1 the pipeline’s
output $\mathcal{Y}$ can arbitrarily depend on some intermediate variables in
the pipeline. This relation is modeled by the conditional distribution
$p\big{(}\mathcal{Y}_{n}|Parent(\mathcal{Y}_{n})\big{)}$. In Fig.S1 the lower
boxes are the inducing points and other variables that determine the GPs’
posterior. More precisely, in Fig.S1
$\\{\tilde{\boldsymbol{x}}_{m}\\}_{m=1}^{M}$ are some inducing points (e.g.
some training images). Vectors in the kernel space are denoted by
$\tilde{\boldsymbol{u}}$ and $\boldsymbol{u}$. Moreover, the observed values
are denoted by $v$ and $\tilde{v}$. Informally, $\boldsymbol{u}$ and $v$
denote the input/output of the GPs. When referring to one of the $M$ inducing
points a ”tilde” is used (as $(\tilde{\boldsymbol{u}},\tilde{v})$), however
$(\boldsymbol{u},v)$ corresponds to a point that can be anywhere in the
kernel-space.

The inducing instances $\\{\tilde{\boldsymbol{x}}_{m}\\}_{m=1}^{M}$ are mapped
to the kernel-spaces by the kernel mappings $\\{f_{1}(.),...,f_{L}(.)\\}$. In
Fig.S1 the variables $\\{\tilde{\boldsymbol{u}}_{m}\\}_{m=1}^{M}$ are the
kernel-space representations of the inducing points
$\\{\tilde{\boldsymbol{x}}_{m}\\}_{m=1}^{M}$. Moreover,
$\\{\tilde{\boldsymbol{v}}_{m}\\}_{m=1}^{M}$ are the GP’s output values at the
inducing points. Given an instance $\boldsymbol{x}_{n}$, it is firstly fed to
the kernel mappings $\\{f_{1}(.),...,f_{L}(.)\\}$ and the kernel-space
representations $\boldsymbol{u}_{n}$ are obtained. Afterwards, the GPs’
outputs on $\boldsymbol{u}_{n}$ depend on $\boldsymbol{u}_{n}$ as well as all
other inducing points because the inducing points actually determine the GPs’
posterior on all kernel-space points including $\boldsymbol{u}_{n}$.
Therefore, in Fig.S1 the variable $\boldsymbol{v}_{n}$ is not only connected
to $\boldsymbol{u}_{n}$ but it is also connected to the box at the bottom
(i.e. all inducing points and other variables associated with them).

As usual, the variational lower-bound is equal to

$\mathcal{L}=\mathbb{E}_{\tiny{\sim q}}\big{[}\log p(\text{all
variables})\big{]}\;-\;\mathbb{E}_{\tiny{\sim q}}\big{[}\log q(\text{hidden
variables})\big{]}.$ (S3)

The likelihood of all variables in Eq.S3 factorizes as the product of
conditional distributions of each variable given its parents. Therefore

$\begin{split}p(\text{all
variables})=\prod_{variable\;t}p\big{(}t|Parent(t)\big{)}.\end{split}$ (S4)

In Eq.S4 only some conditional distributions appear in our derivations which
are discussed at the following.

  * •

The variable $\boldsymbol{x}_{n}$: the ANN’s input $\boldsymbol{x}_{n}$ can
depend arbitrarily on some other intermediate variables in the pipeline. In
our derivations we leave this conditional distribution as
$p\big{(}\boldsymbol{x}_{n}|Parent(\boldsymbol{x}_{n})\big{)}$.

  * •

The variable $\boldsymbol{u}_{n}$: Given a training instance
$\boldsymbol{x}_{n}$, the kernel-space representations $\boldsymbol{u}_{n}$
are deterministically obtained by feeding the instance to the kernel-mappings
$[f_{1}(.),...,f_{L}(.)]$.

  * •

The variable $\boldsymbol{v}_{n}$: The ANN’s output is required to depend only
on the input, so

$p\big{(}\boldsymbol{v}_{n}|Parent(\boldsymbol{v}_{n})\big{)}=p\big{(}\boldsymbol{v}_{n}|\boldsymbol{u}_{n},\boldsymbol{x}_{n},\\{\tilde{\boldsymbol{x}}_{m},\tilde{\boldsymbol{u}}_{m},\tilde{v}_{m}\\}_{m=1}^{M}\big{)}.$
(S5)

The above distribution is actually the GPs’ posterior at $\boldsymbol{u}_{n}$
(i.e. the normal distribution of Eq.1 of the main article).

  * •

The variable $\tilde{\boldsymbol{x}}_{m}$: the inducing point
$\tilde{\boldsymbol{x}}_{m}$ can depend arbitrarily on some other intermediate
variables in the pipeline. In our derivations we leave this conditional
distribution as
$p\big{(}\tilde{\boldsymbol{x}}_{m}|Parent(\tilde{\boldsymbol{x}}_{m})\big{)}$.

  * •

The variable $\tilde{\boldsymbol{u}}_{m}$: Given an inducing point
$\tilde{\boldsymbol{x}}_{m}$, the kernel-space representations
$\tilde{\boldsymbol{u}}_{m}$ are deterministically obtained by feeding the
inducing point $\tilde{\boldsymbol{x}}_{m}$ to the kernel-mappings
$[f_{1}(.),...,f_{L}(.)]$.

  * •

The variables $\hat{\boldsymbol{v}}_{m}$: Given the kernel-space
representations $\\{\tilde{\boldsymbol{u}}^{(\ell)}_{m}\\}_{m=1}^{M}$, the
variables $\\{\tilde{v}^{(\ell)}_{1},...,\tilde{v}^{(\ell)}_{M}\\}$ follow a
$M$-dimensional Gaussian distribution with zero mean and a covariance matrix
determined by the GP prior covariance among the variables
$\\{\tilde{\boldsymbol{u}}^{(\ell)}_{m}\\}_{m=1}^{M}$.

  * •

The variable $\mathcal{Y}_{n}$: the pipeline’s output $\mathcal{Y}$ can
arbitrarily depend on some intermediate variables in the pipeline. In our
derivations we leave this conditional distribution as
$p\big{(}\mathcal{Y}_{n}|Parent(\mathcal{Y}_{n})\big{)}$.

According to Eq.S4, the likelihood of all variables factorizes as

$\begin{split}p(&\text{all
variables})=\prod_{variable\;t}p\big{(}t|Parent(t)\big{)}\\\
&=\big{(}\prod_{n}p(\boldsymbol{x}_{n}|Parent(\boldsymbol{x}_{n}))\big{)}\times\big{(}\prod_{n}p(\boldsymbol{u}_{n}|\boldsymbol{x}_{n})\big{)}\times\\\
&\big{(}\prod_{n}\prod_{\ell}p(v^{(\ell)}_{n}|\boldsymbol{u}_{n},\boldsymbol{x}_{n},\\{\tilde{\boldsymbol{x}}_{m},\tilde{\boldsymbol{u}}_{m},\tilde{v}_{m}\\}_{m=1}^{M})\big{)}\times\\\
&\big{(}\prod_{m}p(\tilde{\boldsymbol{x}}_{m}|Parent(\tilde{\boldsymbol{x}}_{m})\big{)}\times\big{(}\prod_{m}p(\tilde{\boldsymbol{u}}_{m}|\tilde{\boldsymbol{x}}_{m})\big{)}\times\\\
&\big{(}\prod_{\ell}p(\tilde{\boldsymbol{v}}^{(\ell)}_{1:M}|\mathbf{0},\mathcal{K}_{prior}(\tilde{\boldsymbol{u}}^{(\ell)}_{1:M},\tilde{\boldsymbol{u}}^{(\ell)}_{1:M}))\big{)}\times\\\
&\big{(}\prod_{n}p(\mathcal{Y}_{n}|Parent(\mathcal{Y}_{n}))\big{)}\times\big{(}\prod_{\text{other
vars $t$}}p\big{(}t|Parent(t)\big{)}\big{)}.\end{split}$ (S6)

Now we derive the lower-bound $\mathcal{L}$ with respect to each parameter
separately.

###  S1.1 Deriving the Lower-bound With Respect to the Kernel-mappings

In the right-hand-side of Eq.S6 only the following terms are dependant on the
kernel-mappings $[f_{1}(.),...,f_{L}(.)]$:

$\begin{split}\big{[}\prod_{m}p(\tilde{\boldsymbol{u}}_{m}|\tilde{\boldsymbol{x}}_{m})\times\prod_{\ell}p(\tilde{v}^{(\ell)}_{m}|\mathbf{0},\mathcal{K}_{prior}(\tilde{\boldsymbol{u}}^{(\ell)}_{1:M},\tilde{\boldsymbol{u}}^{(\ell)}_{1:M}))\big{]}\times\\\
\big{[}\prod_{n}p\big{(}\boldsymbol{u}_{n}|\boldsymbol{x}_{n}\big{)}\times\prod_{\ell}p\big{(}v^{(\ell)}_{n}|\boldsymbol{u}_{n},\boldsymbol{x}_{n},\\{\tilde{\boldsymbol{x}}_{m},\tilde{\boldsymbol{u}}_{m},\tilde{v}_{m}\\}_{m=1}^{M}\big{)}\big{]}.\end{split}$
(S7)

Please note that in the above equation the terms
$p(\tilde{\boldsymbol{u}}_{m}|\tilde{\boldsymbol{x}}_{m})$ and
$p(\boldsymbol{u}_{n}|\boldsymbol{x}_{n})$ are equal to 1 because
$\tilde{\boldsymbol{u}}_{m}$ and $\boldsymbol{u}_{n}$ are deterministically
obtained from $\tilde{\boldsymbol{x}}_{m}$ and $\boldsymbol{x}_{n}$.
Therefore, in Eq.S3 the terms containing the kernel mappings
$[f_{1}(.),...,f_{L}(.)]$ are as follows:

$\begin{split}\mathcal{L}_{f}&=\mathbb{E}_{\tiny{\sim
q}}\big{[}\sum_{\ell}\log
p(v^{(\ell)}|\boldsymbol{u},\boldsymbol{x},\\{\tilde{\boldsymbol{x}}_{m},\tilde{\boldsymbol{u}}_{m},\tilde{v}_{m}\\}_{m=1}^{M})\big{]}+\\\
&\;\;\;\;\sum_{\ell}\mathbb{E}_{\tiny{\sim q}}\big{[}\log
p(\tilde{\boldsymbol{v}}^{(\ell)}_{1:M}|\mathbf{0},\mathcal{K}_{prior}(\tilde{\boldsymbol{u}}^{(\ell)}_{1:M},\tilde{\boldsymbol{u}}^{(\ell)}_{1:M}))\big{]}-\\\
&\;\;\;\;\sum_{\ell}\mathbb{E}_{\tiny{\sim q}}\big{[}\log
q_{2}(\tilde{\boldsymbol{v}}^{(\ell)}_{1:M})\big{]}\\\
&=\mathbb{E}_{\tiny{\sim q}}\big{[}\sum_{\ell}\log
p(v^{(\ell)}|\boldsymbol{u},\boldsymbol{x},\\{\tilde{\boldsymbol{x}}_{m},\tilde{\boldsymbol{u}}_{m},\tilde{v}_{m}\\}_{m=1}^{M})\big{]}-\\\
&\sum_{\ell=1}\mathbb{E}_{\tiny{\sim
q}}\big{[}KL\Big{(}q_{2}(\tilde{\boldsymbol{v}}^{(\ell)}_{1:M})\;\;||\;\;p(\tilde{\boldsymbol{v}}^{(\ell)}_{1:M}|\mathbf{0},\mathcal{K}_{prior}(\tilde{\boldsymbol{u}}^{(\ell)}_{1:M},\tilde{\boldsymbol{u}}^{(\ell)}_{1:M}))\Big{)}\big{]}.\end{split}$
(S8)

We simplify the two terms on the right-hand-side of Eq.S8. The first term is
the expected log-likelihood of a Gaussian distribution (i.e. the conditional
log-likelihood of $\tilde{v}^{\ell}$ as in Eq.1 of the main article). Also the
variational distribution $q(.)$ is Gaussian. Therefore, we can use Lemma.2 to
simplify the first term:

$\begin{split}&\mathbb{E}_{\tiny{\sim q}}\big{[}\sum_{\ell=1}^{L}\log
p(v^{(\ell)}|\boldsymbol{u},\boldsymbol{x},\\{\tilde{\boldsymbol{x}}_{m},\tilde{\boldsymbol{u}}_{m},\tilde{v}_{m}\\}_{m=1}^{M})\big{]}=\\\
&\sum_{\ell=1}^{L}\mathbb{E}_{\tiny{\sim q}}\big{[}\log
p(v^{(\ell)}|\boldsymbol{u},\boldsymbol{x},\\{\tilde{\boldsymbol{x}}_{m},\tilde{\boldsymbol{u}}_{m},\tilde{v}_{m}\\}_{m=1}^{M})\big{]}=\\\
&\sum_{\ell=1}^{L}\Big{[}-\frac{\big{(}\mu_{v}(\boldsymbol{u}^{(\ell)},\tilde{\mathbf{u}}^{(\ell)}_{1:M},\tilde{v}^{(\ell)}_{1:M})-g_{\ell}(\boldsymbol{x})\big{)}^{2}+\sigma_{g}^{2}}{cov_{v}(\boldsymbol{u}^{(\ell)},\tilde{\mathbf{u}}^{(\ell)}_{1:M},\tilde{v}^{(\ell)}_{1:M})}\\\
&\;\;\;\;\;\;\;\;-\frac{1}{2}\log\big{(}cov_{v}(\boldsymbol{u}^{(\ell)},\tilde{\mathbf{u}}^{(\ell)}_{1:M},\tilde{v}^{(\ell)}_{1:M})\big{)}\\\
&\;\;\;\;\;\;\;\;-\frac{1}{2}\log(2\pi)\Big{]}.\end{split}$ (S9)

Please note that the two terms of Eq.S9 are the two terms which were presented
and discussed in Eq.5 of the main article.

Now we simplify the KL-term on the right-hand-side of Eq.S8. According to
Lemma.1 we have that

$\begin{split}KL\Big{(}q_{2}(\tilde{\boldsymbol{v}}^{(\ell)}_{1:M})\;\;&||\;\;p(\tilde{\boldsymbol{v}}^{(\ell)}_{1:M}|\mathbf{0},\mathcal{K}_{prior}(\tilde{\boldsymbol{u}}^{(\ell)}_{1:M},\tilde{\boldsymbol{u}}^{(\ell)}_{1:M}))\Big{)}=\\\
&+0.5\big{(}\log(\frac{\sigma_{gp}^{2}}{\sigma_{\varphi}^{2}})\big{)}\\\
&-0.5M\\\ &+\frac{\sigma_{\varphi}^{2}}{\sigma_{gp}^{2}}\\\
&+\frac{\boldsymbol{\varphi}^{(\ell)^{T}}_{1:M}\boldsymbol{\varphi}^{(\ell)}_{1:M}}{\sigma_{gp}^{2}},\end{split}$
(S10)

where $\varphi$ are the variational parameters of $q_{2}(.)$ as in Eq.4 of the
main article. Therefore, the KL-term of Eq.S8 is a constant with respect to
the kernel mappings $[f_{1}(.),...,f_{L}(.)]$ and can be discarded. All in
all, the lower-bound for optimizing the kernel-mappings is equal to the right-
hand-side of Eq.S9 which was introduced and discussed in Sec.3.3. of the main
article.

###  S1.2 Deriving the Lower-bound With Respect to the ANN Parameters

According to Eq.4 of the main article, in our formulation the ANN’s parameters
appear as some variational parameters. Therefore, the likelihood of all
variables (Eq.S6) does not generally depend on the ANN’s parameters. But
according to the general ELBO formulation in Eq.S3 the ELBO $\mathcal{L}$
depends on ANN’s parameters, because when computing the expectation the
variables are drawn from the variational distribution $q(.)$. We estimated the
ELBO of Eq.S3 by the average over few samples. More precisely, given a
training instance $\boldsymbol{x}$, we firstly computed the kernel-space
representations as:

$\begin{split}\boldsymbol{u}^{(\ell)}=f_{\ell}(\boldsymbol{x}),\;\;1\leq\ell\leq
L.\end{split}$ (S11)

Afterwards, we used the reparametrization trick for Eq.1 of the main article
to draw a sample for $\boldsymbol{v}^{(\ell)}$ as follows:

$\begin{split}&z_{q2}^{(\ell)}\;\sim\;\mathcal{N}(0,1),\\\
&v^{(\ell)}\;\sim\;\mu_{v}(\boldsymbol{u}^{(\ell)},\tilde{\mathbf{u}}^{(\ell)}_{1:M},\tilde{v}^{(\ell)}_{1:M})+z_{q2}^{(\ell)}cov_{v}(\boldsymbol{u}^{(\ell)},\tilde{\mathbf{u}}^{(\ell)}_{1:M},\tilde{v}^{(\ell)}_{1:M}),\end{split}$
(S12)

where $\mu_{v}(.,.,.)$ and $cov_{v}(.,.,.)$ are defined in Eqs.2 and 3 of the
main article. Moreover, we continue the forward pass of the original pipeline
to get a sample $\mathcal{Y}$. Having drawn $\boldsymbol{x}$,
$\boldsymbol{u}$, $v$, and $\mathcal{Y}$ from the variational distribution, we
estimate the ELBO of Eq.S3 by these samples.

$\begin{split}\mathcal{L}&=\\\ &\mathbb{E}_{\tiny{\sim q}}\big{[}\log
p(\text{all variables})\big{]}-\mathbb{E}_{\tiny{\sim q}}\big{[}\log
q(\text{hidden variables})\big{]}\\\ &\approx\log p(\text{all
variables})\Big{|}_{\boldsymbol{x},\boldsymbol{u},v,\mathcal{Y}}-\sum_{m}^{M}\sum_{\ell}^{L}\mathbb{E}_{\tiny{\sim
q_{2}}}\big{[}\log q_{2}(\tilde{v}^{(\ell)}_{m})\big{]}\end{split}$ (S13)

In the above equation, the second term on the right-hand-side is the entropy
of a normal distribution and it only depends on the variance of the $q_{2}$
distribution. As we let the variance of $q_{2}$ be fixed ($\sigma_{g}^{2}$ in
Eq.4 of the main article), the second term is a constant. Therefore,

$\begin{split}\mathcal{L}\approx\log p(\text{all
variables})\Big{|}_{\boldsymbol{x},\boldsymbol{u},v,\mathcal{Y}}.\end{split}$
(S14)

Among the likelihood term on the right-hand-side of Eq.S6 the conditional
distribution of all variables before $\boldsymbol{u}_{n}$ (e.g.
$\boldsymbol{x}_{n}$ and $\mathcal{X}_{n}$) are independent of the ANN’s
parameters (i.e. the parameters of the function $g(.)$). On the other hand,
for all variables that appear after $\boldsymbol{u}_{n}$, the conditional
distribution depends on the ANN’s parameters. Indeed, according to Eq.S14

$\begin{split}\mathcal{L}_{ann}&\approx\big{[}\sum_{\ell=1}^{L}\log
p(v^{(\ell)}|\boldsymbol{u},\boldsymbol{x},\\{\tilde{\boldsymbol{x}}_{m},\tilde{\boldsymbol{u}}_{m},\tilde{v}_{m}\\}_{m=1}^{M})\big{]}\Big{|}_{\boldsymbol{x},\boldsymbol{v}}+\\\
&\;\;\;\;\log
p(\mathcal{Y}|Parent(\mathcal{Y}))\Big{|}_{\boldsymbol{x},\boldsymbol{v},\mathcal{Y}}+\\\
&\;\;\;\;\big{(}\sum_{\text{other vars after $\boldsymbol{u}_{n}$}}\log
p(t|Parent(t))\big{)}\Big{|}_{\boldsymbol{x},\boldsymbol{v},\mathcal{Y}}.\end{split}$
(S15)

In the above equation, the first term on the right-hand-side is the log-
likelihood of the normal distribution of Eq.1:

$\begin{split}&\log
p(v^{(\ell)}|\boldsymbol{u},\boldsymbol{x},\\{\tilde{\boldsymbol{x}}_{m},\tilde{\boldsymbol{u}}_{m},\tilde{v}_{m}\\}_{m=1}^{M})=\\\
&\;\;\;-\frac{1}{2}\;\big{[}\sum_{\ell=1}^{L}\frac{(\mu_{v}(\boldsymbol{u}^{(\ell)},\tilde{\mathbf{u}}^{(\ell)}_{1:M},\tilde{v}^{(\ell)}_{1:M})-g_{\ell}(\boldsymbol{x}))^{2}}{cov_{v}(\boldsymbol{u}^{(\ell)},\tilde{\mathbf{u}}^{(\ell)}_{1:M},\tilde{v}^{(\ell)}_{1:M})}\big{]}\\\
&\;\;\;+\big{(}\text{some terms independent from $g(.)$}\big{)}.\end{split}$
(S16)

In Eq.S15 the term $p\big{(}\mathcal{Y}|Parent(\mathcal{Y})\big{)}$ is the
likelihood of the output(s) of the whole pipeline as illustrated by Fig.1 of
the main article, given the ANN’s output and all other intermediate variables
on which the final output $\mathcal{Y}$ depends. This likelihood turns out to
be equivalent to commonly-used losses like the cross-entropy loss or the mean-
squared loss. Here we elaborate upon how this happens. Let the task be a
classification, and let $\hat{\mathcal{Y}}\in\mathbb{R}^{L}$ be the pipeline’s
output. The final model prediction $\mathcal{Y}$ is done as follows:

$\mathcal{Y}\sim\mathcal{C}ategorical(\hat{\mathcal{Y}}_{K},...,\hat{\mathcal{Y}}_{K})$
(S17)

Therefore we have that

$p\big{(}\mathcal{Y}|Parent(\mathcal{Y})\big{)}=(\hat{\mathcal{Y}}_{1})^{I[\mathcal{Y}==1]}\times...\times(\hat{\mathcal{Y}}_{K})^{I[\mathcal{Y}==K]},$
(S18)

where $I[.]$ is the indicator function. So, we have that

$\begin{split}&\log p\big{(}\mathcal{Y}|Parent(\mathcal{Y})\big{)}=\\\
&\;\;\;{I[\mathcal{Y}==1]}\log(\hat{\mathcal{Y}}_{1})+...+I[\mathcal{Y}==K]\log(\hat{\mathcal{Y}}_{K}).\end{split}$
(S19)

Therefore, when the pipeline is for classification, $\log
p(\mathcal{Y}|\mathbf{v},\;etc.)$ will be equal to the cross-entropy loss.
This conclusion was introduced and discussed in Eq.6 of the main article. We
can draw similar conclusions when the pipeline is for other tasks like
regression, or even a combination of tasks.

In the general pipeline of Fig.S1, if all stages after $\boldsymbol{v}$ are
deterministic (of course except the final stage which is probabilistic like
Eq.S17), the third term on the right-hand-side of Eq.S15 becomes 1. Therefore,
the right-hand-side of Eq.S15 is equal to Eq.6 of the main article. As we
discussed in Sec.3.3 of the main article, $\mathcal{L}_{ann}$ has two terms:
the first terms encourages the GP-ANN analogy and the second term seeks to
lower the task-loss.

###  S1.3 Deriving the Lower-bound With Respect to $q_{2}(.)$ Parameters

In Eq.4 of the main article we considered the variational parameters
$\\{\varphi^{(\ell)}_{m}\\}_{m=1}^{M}$ for the hidden variables
$\\{\tilde{v}^{(\ell)}_{m}\\}_{m=1}^{M}$. The ELBO of Eq.S3 can be optimized
with respect to $\\{\varphi^{(\ell)}_{m}\\}_{m=1}^{M}$ as well. But we noticed
that optimizing $\\{\varphi^{(\ell)}_{m}\\}_{m=1}^{M}$ is computationally
unstable. Therefore, we set $\\{\varphi^{(\ell)}_{m}\\}_{m=1}^{M}$ according
to the following rule:

$\begin{split}&\varphi_{m}^{(\ell)}=g_{\ell}(\tilde{\boldsymbol{x}}_{m}),\\\
&\;\;\;\;1\leq m\leq M,\;\;1\leq\ell\leq L.\end{split}$ (S20)

We set $\\{\varphi^{(\ell)}_{m}\\}_{m=1}^{M}$ as above because
$\tilde{v}^{(\ell)}_{m}$ is simply the $\ell$-th GP posterior mean at the
inducing point $\tilde{\boldsymbol{x}}_{m}$. To make the GP’s posterior mean
equal to the ANN’s output, $\tilde{v}^{(m)}_{\ell}$ should be equal to the
ANN’s (i.e. $g(.)$’s) output at the $m$-inducing point.

##  S2 Efficiently Computing Gaussian Process Posterior

Let $\boldsymbol{A}$ be an arbitrary $M\times D$ matrix where $M>>D$.
Moreover, let $\boldsymbol{b}$ be a $M$-dimensional vector and let $\sigma$ be
a real number. The computational techniques [30] allow us to efficiently
compute:

$(\mathbf{A}\mathbf{A}^{T}+\sigma^{2}\boldsymbol{I}_{M\times
M})^{-1}\;\boldsymbol{b}.$

The idea is that $\mathbf{A}\mathbf{A}^{T}$ and therefore its inverse are of
rank $D$. Therefore, $(\mathbf{A}\mathbf{A}^{T})^{-1}$ has $D$ non-zero
eigenvalues like $\\{\lambda_{1},...,\lambda_{D}\\}$ and the rest of its
eigenvalues are zero. Let the corresponding eigenvectors be
$\\{\boldsymbol{e_{1}},...,\boldsymbol{e_{D}}\\}$. To compute
$(\mathbf{A}\mathbf{A}^{T})^{-1}\boldsymbol{b}$ we can simply project
$\boldsymbol{b}$ to the $D$-dimensional space of the eigenvectors. By doing
so, we avoid the $\mathcal{O}(M^{3})$ computational complexity. Let
$\\{\lambda_{1},...,\lambda_{D}\\}$ be the non-zero eigenvalues of
$\mathbf{A}\mathbf{A}^{T}$ and let
$\\{\boldsymbol{e}_{1},...,\boldsymbol{e}_{D}\\}$ be the corresponding
eigenvectors. From linear algebra, it follows that for
$\mathbf{A}\mathbf{A}^{T}+\sigma^{2}\mathbf{I}_{M\times M}$ the eigenvalues
and the eigenvectors are
${\\{\lambda_{1}+\sigma^{2},...,\lambda_{D}+\sigma^{2},\sigma^{2},...,\sigma^{2}\\}}$
and $\\{\boldsymbol{e_{1}},...,\boldsymbol{e_{D}}\\}$, respectively. Please
note that $M-D$ eigenvectors are added all of which are equal to $\sigma^{2}$.
Similarly, from linear algebra it follows that for the inverse of
${\mathbf{A}\mathbf{A}^{T}+\sigma^{2}\mathbf{I}_{M\times M}}$ the eigenvalues
and eigenvectors are
$\\{\frac{1}{\lambda_{1}+\sigma^{2}},...,\frac{1}{\lambda_{D}+\sigma^{2}},\frac{1}{\sigma^{2}},...,\frac{1}{\sigma^{2}}\\}$
and
$\\{\boldsymbol{e_{1}},...,\boldsymbol{e_{D}},\boldsymbol{e_{D+1}},...,\boldsymbol{e_{M}}\\}$
respectively. Please note that although there are $M$ eigenvectors, only the
first $D$ eigenvectors appear in our computations. More precisely, let
$\mathbf{E}\in\mathbb{R}^{M\times D}$ be a matrix whose columns are
$\\{\boldsymbol{e_{1}},...,\boldsymbol{e_{D}}\\}$. Let $\boldsymbol{\Lambda}$
be a diagonal matrix whose diagonal is formed by
$\\{\frac{1}{\lambda_{1}+\sigma^{2}},...,\frac{1}{\lambda_{D}+\sigma^{2}}\\}$.
In the space of the $D$ eigenvectors the linear transformation on any vector
like $\boldsymbol{b}$ is equal to
$\mathbf{E}\boldsymbol{\Lambda}\mathbf{E}^{T}\boldsymbol{b}$, meaning that
multiplication by $\mathbf{E}^{T}$ transforms $\boldsymbol{b}$ to the space of
the $D$ eigenvectors, multiplication by $\boldsymbol{\Lambda}$ performs the
transformation in that space, and multiplication by $\mathbf{E}$ transforms
the result back to the original space. The $(M-D)$ eigenvalues that correspond
to the rest of the eigenvectors are all the same and are equal to
$\frac{1}{\sigma^{2}}$. Therefore, there is no need to project
$\boldsymbol{b}$ to the space of the $(M-D)$ eigenvectors because the linear
transformation in that space is simply a scaling by $\frac{1}{\sigma^{2}}$.
All in all, we have that

$\big{(}\mathbf{A}\mathbf{A}^{T}+\sigma^{2}\mathbf{I}_{M\times
M}\big{)}^{-1}\boldsymbol{b}\;=\;\mathbf{E}\boldsymbol{\Lambda}\mathbf{E}^{T}\boldsymbol{b}+\frac{1}{\sigma^{2}}(\boldsymbol{b}-\mathbf{E}\mathbf{E}^{T}\boldsymbol{b}).$
(S21)

Complexity of computing the right-hand-side of Eq.S21 is way lower than the
$\mathcal{O}(M^{3})$ requirement of the standard matrix inversion. We borrowed
more computational ideas from the work on fast spectral clustering [30]. To
compute the first $D$ eigenvlaues and eigenvectors of
$\mathbf{A}\mathbf{A}^{T}$, we worked with the $D$-by-$D$ matrix
$\mathbf{A}^{T}\mathbf{A}$ rather than the $M$-by-$M$ matrix
$\mathbf{A}\mathbf{A}^{T}$ (recall that $D<<M$), because given the eigenvalues
and eigenvectors of $\mathbf{A}^{T}\mathbf{A}$, those of
$\mathbf{A}\mathbf{A}^{T}$ are easily computable [30]. The procedure is
explained in Alg.5. In Alg.5, lines 1-3 compute the eigenvalues/vectors of the
matrix $\mathbf{A}^{T}\mathbf{A}$. Afterwards, lines 4 and 5 compute the first
$D$ eigenvalues/vectors of $\mathbf{A}\mathbf{A}^{T}$ using those of
$\mathbf{A}^{T}\mathbf{A}$. Finally, line 8 computes
$(\mathbf{A}\mathbf{A}^{T}+\sigma^{2}\mathbf{I})^{-1}\boldsymbol{b}$ according
to the right-hand-side of Eq.S21. To make the computations faster, we made use
of the following equation
$\mathbf{A}\mathbf{A}^{T}=\sum_{m}\mathbf{A}[m,:]\mathbf{A}[m,:]^{T}$, where
${\mathbf{A}[m,:]}$ is the $m$-th row of the matrix $\mathbf{A}$. Thanks to
this equation, we compute $\mathbf{A}\mathbf{A}^{T}$ only once at the
beginning of the training. Afterwards, as each mini-batch alters only some
rows of $\mathbf{A}$, we update the previously computed
$\mathbf{A}\mathbf{A}^{T}$ by considering only the effect of the modified
rows.

##  S3 Computing Pixel Contributions to the Similarity

Let the kernel mapping $f(.)$ be a convolutional neural network that produces
a volumetric map of size $C\times H\times W$ followed by a spatial average
pooling that produces the $C$-dimensional vector in the kernel-space. In this
case, $\mathcal{K}(\boldsymbol{x}_{1},\boldsymbol{x}_{2})$ is as follows:

$\begin{split}\mathcal{K}(\boldsymbol{x}_{1},\boldsymbol{x}_{2})&=f(\boldsymbol{x}_{1})^{T}\;f(\boldsymbol{x}_{2})\\\
&=\big{(}\sum_{i=1}^{H}\sum_{j=1}^{W}\boldsymbol{z}^{(1)}_{ij}\big{)}^{T}\big{(}\sum_{k=1}^{H}\sum_{\ell=1}^{W}\boldsymbol{z}^{(2)}_{k\ell}\big{)}\\\
&=\sum_{i=1}^{H}\sum_{j=1}^{W}\sum_{k=1}^{H}\sum_{\ell=1}^{W}\big{(}\boldsymbol{z}^{(1)^{T}}_{ij}\boldsymbol{z}^{(2)}_{k\ell}\big{)},\end{split}$
(S22)

where $\boldsymbol{z}^{(1)}$ and $\boldsymbol{z}^{(2)}$ are the volumetric
maps of size ${C\times H\times W}$ and the indices $(i,j)$ and $(k,\ell)$
index the spatial locations over the volumetric maps. The last term in Eq.S22
shows that the total similarity
$\mathcal{K}(\boldsymbol{x}_{1},\boldsymbol{x}_{2})$ is the sum of the
contributions from each pair of positions $(i,j)$ on $\boldsymbol{x}_{1}$ and
$(k,\ell)$ on $\boldsymbol{x}_{2}$. To compute the contribution of a specific
location like $(i,j)$ on $\boldsymbol{x}_{1}$, we sum up the contributions of
$(i,j)$ on $\boldsymbol{x}_{1}$ and all possible locations
${\\{(k,\ell)\\}_{k=1}^{H}{}_{\ell=1}^{W}}$ on $\boldsymbol{x}_{2}$.

The kernel-mappings that we used have a slightly different architecture than a
volumetric map followed by spatial average pooling. Our kernel mappings
produce a volumetric map of size $C\times H\times W$ followed by a spatial
average pooling that produces a $C$-dimensional vector. Afterwards, the
resulting vector is divided by its $\ell_{2}$-norm to produce a vector of norm
1. Consequently, this vector of norm 1 is fed to a leaky ReLU layer that
produces the final kernel-space representation $f(\boldsymbol{x})$. For this
architecture the pixel contributions can be computed according to an equation
similar to Eq.S22 as follows. Our kernel mappings produce the volumetric map
$\boldsymbol{z}$ of size $C\times H\times W$ followed by a spatial average
pooling that produces the $C$-dimensional vector $\boldsymbol{a}$:

$\boldsymbol{a}=\sum_{i=1}^{H}\sum_{j=1}^{W}\boldsymbol{z}_{ij}.$ (S23)

Afterwards, the resulting vector is divided by its $\ell_{2}$-norm to produce
the vector $\boldsymbol{b}$ of norm 1:

$\boldsymbol{b}=[\frac{a_{1}}{||\boldsymbol{a}||_{2}},\;...\;,\frac{a_{C}}{||\boldsymbol{a}||_{2}}].$
(S24)

Consequently, this vector of norm 1 is fed to a leaky ReLU layer that produces
the final kernel-space representation $f(\boldsymbol{x})$:

$f(\boldsymbol{x})=leakyReLU(\boldsymbol{b}).$ (S25)

We begin with simplifying Eq.S25. The leaky ReLU activation function
multiplies the input by a constant and this constant depends on the sign of
the input. Therefore, applying the leaky ReLU activation is equivalent to
multiplication by a diagonal matrix $\boldsymbol{\Lambda}$. Therefore,

$f(\boldsymbol{x})=\boldsymbol{\Lambda}\boldsymbol{b}.$ (S26)

Let $\boldsymbol{x}_{1}$ and $\boldsymbol{x}_{2}$ be two images, and
$\boldsymbol{z}^{(1)}$ and $\boldsymbol{z}^{(2)}$ be the corresponding
volumetric maps. We have that

$\begin{split}\boldsymbol{a}^{(1)}=\sum_{i=1}^{H}\sum_{j=1}^{W}\boldsymbol{z}^{(1)}_{ij},\\\
\boldsymbol{a}^{(2)}=\sum_{k=1}^{H}\sum_{\ell=1}^{W}\boldsymbol{z}^{(2)}_{k\ell}.\end{split}$
(S27)

And

$\begin{split}\boldsymbol{b}^{(1)}=[\frac{a^{(1)}_{1}}{||\boldsymbol{a}^{(1)}||_{2}},\;...\;,\frac{a^{(1)}_{C}}{||\boldsymbol{a}^{(1)}||_{2}}],\\\
\boldsymbol{b}^{(2)}=[\frac{a^{(2)}_{1}}{||\boldsymbol{a}^{(2)}||_{2}},\;...\;,\frac{a^{(2)}_{C}}{||\boldsymbol{a}^{(2)}||_{2}}].\end{split}$
(S28)

And

$\begin{split}f(\boldsymbol{x}^{(1)})=\boldsymbol{\Lambda}^{(1)}\boldsymbol{b}^{(1)},\\\
f(\boldsymbol{x}^{(2)})=\boldsymbol{\Lambda}^{(2)}\boldsymbol{b}^{(2)}.\end{split}$
(S29)

Now we simplify the similarity
$\mathcal{K}(\boldsymbol{x}_{1},\boldsymbol{x}_{2})$:

$\begin{split}\mathcal{K}(\boldsymbol{x}_{1},\boldsymbol{x}_{2})&=\big{(}\boldsymbol{\Lambda}^{(1)}\boldsymbol{b}^{(1)}\big{)}^{T}\big{(}\boldsymbol{\Lambda}^{(2)}\boldsymbol{b}^{(2)}\big{)}\\\
&=\big{(}\boldsymbol{\Lambda}^{(1)^{T}}\boldsymbol{\Lambda}^{(2)}\big{)}\big{(}\boldsymbol{b}^{(1)^{T}}\boldsymbol{b}^{(2)}\big{)}\\\
&=\frac{\big{(}\boldsymbol{\Lambda}^{(1)^{T}}\boldsymbol{\Lambda}^{(2)}\big{)}}{||\boldsymbol{a}^{(1)}||_{2}\;\;||\boldsymbol{a}^{(2)}||_{2}}\;\big{(}\sum_{i=1}^{H}\sum_{j=1}^{W}\boldsymbol{z}^{(1)}_{ij}\big{)}^{T}\big{(}\sum_{k=1}^{H}\sum_{\ell=1}^{W}\boldsymbol{z}^{(2)}_{k\ell}\big{)}\\\
&=\frac{\big{(}\boldsymbol{\Lambda}^{(1)^{T}}\boldsymbol{\Lambda}^{(2)}\big{)}}{||\boldsymbol{a}^{(1)}||_{2}\;\;||\boldsymbol{a}^{(2)}||_{2}}\;\sum_{i=1}^{H}\sum_{j=1}^{W}\sum_{k=1}^{H}\sum_{\ell=1}^{W}\big{(}\boldsymbol{z}^{(1)^{T}}_{ij}\boldsymbol{z}^{(2)}_{k\ell}\big{)}.\end{split}$
(S30)

Indeed, as the used architecture for kernel-mappings is slightly different
than producing a volumetric map followed by spatial average pooling, instead
of Eq.S22, we used Eq.S30 that we derived above.

##  S4 Practical Details and Parameter Settings

In this section we discuss some practical details which we have not yet
discussed in this paper. Moreover, we provide the exact parameter settings
that we used throughout our experiments. As explained in Sec.4, there are $L$
kernel mappings that we denoted by $[f_{1}(.),f_{2}(.),...,f_{L}(.)]$. One can
implement this kernel mappings by, e.g., considering $L$ independent CNNs.
However, doing so dramatically increases the computation cost. Therefore, we
modeled the $L$ mappings by a common ResNet-50 [38] backbone. After the common
backbone, we placed $L$ branches. Each branch has two convolutional layers
followed by global spatial average pooling that produce a vector. Each branch
ends with an L2 normalizer layer (that sets the L2-norm of the vector to 1)
followed by a leaky-ReLU layer. During our experiments, we noticed that the
L2-normalization layer and the final leaky-ReLU layer are essential. Without
the L2 normalization layer, the vectors in the kernel-space can have
arbitrarily-small or arbitrarily-big elements, and this makes the training
unstable. We included the last leaky-ReLU layer, because according to GP
posterior mean formula, vectors in the kernel-space go through a linear
transformation. Therefore, without the last leaky-ReLU layer, the pipeline
would have two consequtive linear layers. Throughout our experiments, we set
the output of each branch (i.e. vectors in the kernel-space of each GP) to be
20-dimensional.

As illustrated by Fig.10, we need to make the inducing dataset as large as
possible. Therefore, throughout our experiments we selected the whole training
dataset as the inducing dataset. Unlike training instances, we didn’t apply
data-augmentation on inducing instances. By doing so, the training dataset and
the inducing dataset will have very similar instances. This causes a
difficulty that we are going to discuss in this part. The kernel-mappings
$[f_{1}(.),...,f_{L}(.)]$ are trained according to Alg.4. After selecting an
instance like $\boldsymbol{x}$ from the training dataset, $\boldsymbol{x}$ is
actually the augmented version of an inducing instance like
$\tilde{\boldsymbol{x}}_{m}$. Indeed, we have that
${\boldsymbol{x}=DataAug(\tilde{\boldsymbol{x}}_{m})}$. Because
$\boldsymbol{x}$ and $\tilde{\boldsymbol{x}}_{m}$ are very similar, they will
be very close to one another regardless of what parameters
$[f_{1}(.),...,f_{L}(.)]$ have. Therefore, regardless of the kernel-mappings,
the GP-mean will match the ANN value at $\boldsymbol{x}$, and there will be no
training signal for the kernel-mappings $[f_{1}(.),...,f_{L}(.)]$. Please note
that in this case the GPs match the ANNs only on training instances and not
testing instances. To avoid this issue, we optimized the GP-ANN analogy (i.e.
the objective in Eq.5 of the main article) on instances like
$\lambda\boldsymbol{x}_{i}+(1-\lambda)\boldsymbol{x}_{j}$, where
$\boldsymbol{x}_{i}$ and $\boldsymbol{x}_{j}$ are two instances randomly
selected from the training set and $\lambda$ is a scalar uniformly selected
from $[-1,2]$.

When applying our proposed GPEX we used Adam optimizer [40]. Although the
AMSGrad version [42] of this optimizer is often recommended, for our proposed
GPEX we noticed the Adam optimizer [40] without AMSGrad works the best. For
explaining classifier ANNs, we used a learning-rate of 0.0001 while for
explaining the attention submodules we used a learning rate of 0.00001.

Figure S2: Scatters for Cifar10 (classifier).

Figure S3: Scatters for Cifar10 (attention).

Figure S4: Scatters for MNIST (classifier).

Figure S5: Scatters for MNIST (attention).

Figure S6: Scatters for Kather dataset (classifier).

Figure S7: Scatters for Kather dataset (attention).

Figure S8: Explanations for MNIST (set 1).

Figure S9: Explanations for MNIST (set 2).

Figure S10: Explanations for MNIST (set 3).

Figure S11: Explanations for MNIST (set 4).

Figure S12: Explanations for MNIST (set 5).

Figure S13: Explanations for MNIST (set 6).

Figure S14: Explanations for MNIST (set 7).

Figure S15: Explanations for MNIST (set 8).

Figure S16: Explanations for MNIST (set 9).

Figure S17: Explanations for MNIST (set 10).

Figure S18: Explanations for MNIST (set 11).

Figure S19: Explanations for MNIST (set 12).

Figure S20: Explanations for Kather dataset (set 1).

Figure S21: Explanations for Kather dataset (set 2).

Figure S22: Explanations for Kather dataset (set 3).

Figure S23: Explanations for Kather dataset (set 4).

Figure S24: Explanations for Kather dataset (set 5).

Figure S25: Explanations for DogsWolves (set 1).

Figure S26: Explanations for DogsWolves (set 2).

Figure S27: Explanations for DogsWolves (set 3).

Figure S28: Explanations for DogsWolves (set 4).

Figure S29: Explanations for DogsWolves (set 5).

Figure S30: Explanations for DogsWolves (set 6).

Figure S31: Explanations for DogsWolves (set 7).

Figure S32: Explanations for DogsWolves (set 8).

Figure S33: Explanations for DogsWolves (set 9).

Figure S34: Explanations for DogsWolves (set 10).

Figure S35: Explanations for DogsWolves (set 11).

Figure S36: Explanations for Cifar10 (set 1).

Figure S37: Explanations for Cifar10 (set 2).

Figure S38: Explanations for Cifar10 (set 3).

Figure S39: Explanations for Cifar10 (set 4).

Figure S40: Explanations for Cifar10 (set 5).

Figure S41: Explanations for Cifar10 (set 6).

Figure S42: Explanations for Cifar10 (set 7).

Figure S43: Explanations for Cifar10 (set 8).

Figure S44: Explanations for Cifar10 (set 9).

Figure S45: Explanations for Cifar10 (set 10).

Figure S46: Explanations for Cifar10 (set 11).

Figure S47: Explanations for Cifar10 (set 12).

Figure S48: Explanations for Cifar10 (set 13).

Figure S49: Explanations for Cifar10 (set 14).

Figure S50: Comparing the proposed GPEX with representer point selection [12]
(set 1).

Figure S51: Comparing the proposed GPEX with representer point selection [12]
(set 2).

Figure S52: Comparing the proposed GPEX with representer point selection [12]
(set 3).

Figure S53: Comparing the proposed GPEX with representer point selection [12]
(set 4).

Figure S54: Comparing the proposed GPEX with representer point selection [12]
(set 5).

Figure S55: Comparing the proposed GPEX with representer point selection [12]
(set 6).

Figure S56: Comparing the proposed GPEX with representer point selection [12]
(set 7).

Figure S57: Comparing the proposed GPEX with representer point selection [12]
(set 8).

head 1

head 2

head 3

head 4

head 5

head 6

head 7

head 8

head 9

head 10

Figure S58: Parameter analysis of Sec.6 of the main article.

| Cifar10 [33] | MNIST [32] | Kather [34] | DogsWolves [35]
---|---|---|---|---
ANN accuracy | 95.43 | 99.56 | 96.80 | 80.50
GPs accuracy | 92.26 | 99.41 | 93.60 | 78.75
Table S1: Accuracies of ANN classifiers versus the accuracies of the explainer
GPs on four datasets.

batch 1

batch 2

batch 3

batch 4

Figure S59: Comparing GP and ANN outputs for four batches of Cifar10 dataset
[33]. The red rectangles highlight the instnaces for which the predictions of
GP and ANN (i.e. the class with maximum score) are different.

batch 1

batch 2

batch 3

batch 4

Figure S60: Comparing GP and ANN outputs for four batches of MNIST dataset
[32]. The red rectangles highlight the instnaces for which the predictions of
GP and ANN (i.e. the class with maximum score) are different.

batch 1

batch 2

batch 3

batch 4

Figure S61: Comparing GP and ANN outputs for four batches of Kather dataset
[34]. The red rectangles highlight the instnaces for which the predictions of
GP and ANN (i.e. the class with maximum score) are different.

batch 1

batch 2

batch 3

batch 4

Figure S62: Comparing GP and ANN outputs for four batches of DogsWolves
dataset [35]. The red rectangles highlight the instnaces for which the
predictions of GP and ANN (i.e. the class with maximum score) are different.