Interpretable Counterfactual Explanations Guided by Prototypes

TL;DR

This paper introduces a fast, model-agnostic method for generating interpretable counterfactual explanations using class prototypes, improving speed and interpretability for classifier predictions on image and tabular data.

Contribution

The method leverages class prototypes to accelerate counterfactual search and introduces new metrics for local interpretability evaluation.

Findings

Significantly faster counterfactual generation using prototypes

More interpretable explanations compared to baseline methods

Effective on both image and tabular datasets

Abstract

We propose a fast, model agnostic method for finding interpretable counterfactual explanations of classifier predictions by using class prototypes. We show that class prototypes, obtained using either an encoder or through class specific k-d trees, significantly speed up the the search for counterfactual instances and result in more interpretable explanations. We introduce two novel metrics to quantitatively evaluate local interpretability at the instance level. We use these metrics to illustrate the effectiveness of our method on an image and tabular dataset, respectively MNIST and Breast Cancer Wisconsin (Diagnostic). The method also eliminates the computational bottleneck that arises because of numerical gradient evaluation for $\textit{black box}$ models.

Tables

Table 1. Table 1: Summary statistics with 95 95 95 % confidence bounds for each loss function for the MNIST experiment.

Method	Time (s)	Gradient steps	IM1	IM2 ($\times$10)
A	$13.06\pm 0.23$	$5158\pm 82$	$1.56\pm 0.03$	$1.65\pm 0.04$
B	$8.40\pm 0.38$	$2380\pm 113$	$1.36\pm 0.02$	$1.60\pm 0.03$
C	$3.06\pm$0.11	$835\pm 36$	$1.16\pm 0.02$	$1.09\pm 0.02$
D	$2.31\pm$0.04	$497\pm 10$	$1.21\pm 0.02$	$1.26\pm 0.03$
E	$1.93\pm$0.10	$777\pm 44$	$1.10\pm 0.02$	$1.10\pm 0.03$
F	$4.01\pm$0.05	$1116\pm 14$	$1.19\pm 0.02$	$1.27\pm 0.03$

Table 2. Table 2: Summary statistics with 95 95 95 % confidence bounds for each loss function for the Breast Cancer Wisconsin (Diagnostic) experiment.

Method	Time (s)	Gradient steps	IM1	IM2 ($\times$10)
A	$2.68\pm 0.20$	$2752\pm 203$	$2.07\pm 0.16$	$7.65\pm 0.79$
B	$0.35\pm 0.03$	$253\pm 33$	$0.94\pm 0.10$	$1.47\pm 0.15$
C	$0.22\pm 0.02$	$182\pm 30$	$0.88\pm 0.10$	$1.41\pm 0.15$

Table 3. Table 3: Mean time in seconds needed to compute 100 optimization steps for objective functions A ′ superscript 𝐴 ′ A^{\prime} , E 𝐸 E and F 𝐹 F with 95 95 95 % confidence bounds. A ′ superscript 𝐴 ′ A^{\prime} is the equivalent of A 𝐴 A without access to the model architecture. As a result, we can only query the prediction function and need to evaluate gradients numerically. One test instance is used for each class in MNIST.

Method	Time (s)
A’	$54.64\pm 1.28$
E	$0.53\pm 0.04$
F	$0.72\pm 0.01$

Table 4. Table 4: Total number of counterfactual MNIST instances generated for each loss function and per counterfactual predicted class. We highlight the lowest number of instances per class which defines the sample size for performing the per-class test.

Method	Total	0	1	2	3	4	5	6	7	8	9
A	5375	81	46	737	884	207	610	373	344	1398	695
B	5377	48	8	505	708	98	155	197	233	2819	606
C	5206	398	88	231	788	397	663	386	632	606	1017
D	4425	247	54	201	681	286	529	352	515	600	960
E	5232	399	83	231	788	394	664	387	659	609	1018
F	5120	379	90	226	762	390	627	392	633	604	1017

Equations

$$\underset{\delta }{\min }c\cdot f_{\kappa }(x_0,\delta )+f_{\text{dist}}(\delta ).$$

$$L_{\text{pred}}$$

$$=\max ([f_{\text{pred}}(x_0+\delta ){]}_{t_0}-\underset{i\ne t_0}{\max }[f_{\text{pred}}(x_0+\delta ){]}_i,-\kappa ),$$

$$f_{\text{dist}}(\delta )=\beta \cdot \Vert \delta {\Vert }_1+\Vert \delta {\Vert }_2^2=\beta \cdot L_1+L_2.$$

$$L_{\text{AE}}=\gamma \cdot \Vert x_0+\delta -\text{AE}(x_0+\delta ){\Vert }_2^2.$$

$$L=c\cdot L_{\text{pred}}+\beta \cdot L_1+L_2+L_{\text{AE}}.$$

$$L=c\cdot L_{\text{pred}}+\beta \cdot L_1+L_2+L_{\text{AE}}+L_{\text{proto}},$$

$${\mathrm{proto}}_i:=\frac{1}{K}\sum _{k=1}^K\text{ENC}(x_k^i)$$

$$j=\underset{i\ne t_0}{\mathrm{argmin}}\Vert \text{ENC}(x_0)-{\mathrm{proto}}_i{\Vert }_2.$$

$$L_{\text{proto}}=\theta \cdot \Vert \text{ENC}(x_0+\delta )-{\mathrm{proto}}_j{\Vert }_2^2,$$

$$L_{\text{proto}}=\theta \cdot \Vert x_0+\delta -{\mathrm{proto}}_j{\Vert }_2^2.$$

$$\frac{\partial L_{\text{pred}}}{\partial x}=\frac{\partial f_{\kappa }(x)}{\partial x}=\frac{\partial f_{\kappa }(x)}{\partial f_{\text{pred}}}\frac{\partial f_{\text{pred}}}{\partial x},$$

$$\frac{\partial f_{\text{pred}}}{\partial x_k}\approx \frac{f_{\text{pred}}(x+{ϵ}_k)-f_{\text{pred}}(x-{ϵ}_k)}{2ϵ},$$

$$\text{IM1}({\text{AE}}_i,{\text{AE}}_{t_0},x_{\text{cf}}):=\frac{\Vert x_0+\delta -{\text{AE}}_i(x_0+\delta ){\Vert }_2^2}{\Vert x_0+\delta -{\text{AE}}_{t_0}(x_0+\delta ){\Vert }_2^2+ϵ}.$$

$$\text{IM2}({\text{AE}}_i,\text{AE},x_{\text{cf}}):=\frac{\Vert {\text{AE}}_i(x_0+\delta )-\text{AE}(x_0+\delta ){\Vert }_2^2}{\Vert x_0+\delta {\Vert }_1+ϵ}.$$

$$A$$

$$B$$

$$C$$

$$D$$

$$E$$

$$F$$

$$A$$

$$B$$

$$C$$

$$d(v_1,v_2)=\sum _{i=1}^n{\left|\frac{c_1^i}{c_1}-\frac{c_2^i}{c_2}\right|}^{\alpha },$$

$$d(v_1,v_2)=\sum _{j\ne i}\Psi \left(p(x_j|x_i=v_1),p(x_j|x_i=v_2))\right),$$

Full text

Machine Learning, ICML

1 Introduction

Humans often think about how they can alter the outcome of a situation. What do I need to change for the bank to approve my loan? or Which symptoms would lead to a different medical diagnosis? are common examples. This form of counterfactual reasoning comes natural to us and explains how to arrive at a desired outcome in an interpretable manner. Moreover, examples of counterfactual instances resulting in a different outcome can give powerful insights of what is important to the underlying decision process, making it a compelling method to explain predictions of machine learning models (Figure 1).

In the context of predictive models, given a test instance and the model’s prediction, a counterfactual instance describes the necessary change in input features that alter the prediction to a predefined output (Molnar, 2019). For classification models the predefined output can be any target class or prediction probability distribution. Counterfactual instances can then be found by iteratively perturbing the input features of the test instance until the desired prediction is reached. In practice, the counterfactual search is posed as an optimization problem—we want to minimize an objective function which encodes desirable properties of the counterfactual instance with respect to the perturbations. The key insight of this formulation is the need to design an objective function that allows us to generate high quality counterfactual instances. A counterfactual instance $x_{\text{cf}}$ should have the following desirable properties:

The model prediction on $x_{\text{cf}}$ needs to be close to the predefined output. 2. 2.

The perturbation $\delta$ changing the original instance $x_{0}$ into $x_{\text{cf}}=x_{0}+\delta$ should be sparse. 3. 3.

The counterfactual $x_{\text{cf}}$ needs to be interpretable. We consider an instance $x_{\text{cf}}$ interpretable if it lies close to the model’s training data distribution. This definition does not only apply to the overall data set, but importantly also to the training instances that belong to the counterfactual class. Let us illustrate this with an intuitive example. Assume we are predicting house prices with features including the square footage and the number of bedrooms. Our house is valued below £500,000 and we would like to know what needs to change about the house in order to increase the valuation above £500,000. By simply increasing the number of bedrooms and leaving the other features unchanged, the model predicts that our counterfactual house is now worth more than £500,000. This sparse counterfactual instance lies fairly close to the overall training distribution since only one feature value was changed. The counterfactual is however out-of-distribution with regards to the subset of houses in the training data valued above £500,000 because other relevant features like the square footage still resemble a typical house valued below £500,000. As a result, we do not consider this counterfactual to be very interpretable. We show in the experiments that there is often a trade-off between sparsity and interpretability. 4. 4.

The counterfactual instance $x_{\text{cf}}$ needs to be found fast enough to ensure it can be used in a real life setting.

An overly simplistic objective function may return instances which satisfy properties 1. and 2., but where the perturbations are not interpretable with respect to the counterfactual class.

In this paper we propose using class prototypes in the objective function to guide the perturbations quickly towards an interpretable counterfactual. The prototypes also allow us to remove computational bottlenecks from the optimization process which occur due to numerical gradient calculation for black box models. In addition, we propose two novel metrics to quantify interpretability which provide a principled benchmark for evaluating interpretability at the instance level. We show empirically that prototypes improve the quality of counterfactual instances on both image (MNIST) and tabular (Wisconsin Breast Cancer) datasets. Finally, we propose using pairwise distance measures between categories of categorical variables to define meaningful perturbations for such variables and illustrate the effectiveness of the method on the Adult (Census) dataset.

2 Related Work

Counterfactual instances—synthetic instances of data engineered from real instances to change the prediction of a machine learning model—have been suggested as a way of explaining individual predictions of a model as an alternative to feature attribution methods such as LIME (Ribeiro et al., 2016) or SHAP (Lundberg & Lee, 2017).

Wachter et al. (2018) generate counterfactuals by minimizing an objective function which sums the squared difference between the predictions on the perturbed instance and the desired outcome, and a scaled $L_{1}$ norm of the perturbations. Laugel et al. (2018) find counterfactuals through a heuristic search procedure by growing spheres around the instance to be explained. The above methods do not take local, class specific interpretability into account. Furthermore, for black box models the number of prediction calls during the search process grows proportionally to either the dimensionality of the feature space (Wachter et al., 2018) or the number of sampled observations (Laugel et al., 2018; Dhurandhar et al., 2019), which can result in a computational bottleneck. Dhurandhar et al. (2018, 2019) propose the framework of Contrastive Explanations which find the minimal number of features that need to be changed/unchanged to keep/change a prediction.

A key contributions of this paper is the use of prototypes to guide the counterfactual search process. Kim et al. (2016); Gurumoorthy et al. (2017) use prototypes as example-based explanations to improve the interpretability of complex datasets. Besides improving interpretability, prototypes have a broad range of applications like clustering (Kaufmann & Rousseeuw, 1987), classification (Bien & Tibshirani, 2011; Takigawa et al., 2009), and few-shot learning (Snell et al., 2017). If we have access to an encoder (Rumelhart et al., 1986), we follow the approach of (Snell et al., 2017) who define a class prototype as the mean encoding of the instances which belong to that class. In the absence of an encoder, we find prototypes through class specific k-d trees (Bentley, 1975).

To judge the quality of the counterfactuals we introduce two novel metrics which focus on local interpretability with respect to the training data distribution. This is different from (Dhurandhar et al., 2017) who define an interpretability metric relative to a target model. Kim et al. (2016) on the other hand quantify interpretability through a human pilot study measuring the accuracy and efficiency of the humans on a predictive task. Luss et al. (2019) also highlight the importance of good local data representations in order to generate high quality explanations.

Another contribution of this paper is a principled approach to handling categorical variables during the counterfactual generation process. Some previously proposed solutions are either computationally expensive (Wachter et al., 2018) or do not take relationships between categories into account (Dhurandhar et al., 2019; Mothilal et al., 2020). We propose using pairwise distance measures to define embeddings of categorical variables into numerical space which allows us to define meaningful perturbations when generating counterfactuals.

3 Methodology

3.1 Background

The following section outlines how the prototype loss term is constructed and why it improves the convergence speed and interpretability. Finding a counterfactual instance $x_{\text{cf}}=x_{0}+\delta$, with both $x_{\text{cf}}$ and $x_{0}$ $\in$ $\mathcal{X}\subseteq\mathbb{R}^{D}$ where $\mathcal{X}$ represents the $D$-dimensional feature space, implies optimizing an objective function of the following form:

[TABLE]

$f_{\kappa}(x_{0},\delta)$ encourages the predicted class $i$ of the perturbed instance $x_{\text{cf}}$ to be different than the predicted class $t_{0}$ of the original instance $x_{0}$. Similar to (Dhurandhar et al., 2018), we define this loss term as:

[TABLE]

where $[f_{\text{pred}}(x_{0}+\delta)]_{i}$ is the $i$-th class prediction probability, and $\kappa\geq 0$ caps the divergence between $[f_{\text{pred}}(x_{0}+\delta)]_{t_{0}}$ and $[f_{\text{pred}}(x_{0}+\delta)]_{i}$. The term $f_{\text{dist}}(\delta)$ minimizes the distance between $x_{0}$ and $x_{\text{cf}}$ with the aim to generate sparse counterfactuals. We use an elastic net regularizer (Zou & Hastie, 2005):

[TABLE]

While the objective function (1) is able to generate counterfactual instances, it does not address a number of issues:

$x_{\text{cf}}$ does not necessarily respect the training data manifold, resulting in out-of-distribution counterfactual instances. Often a trade off needs to be made between sparsity and interpretability of $x_{\text{cf}}$. 2. 2.

The scaling parameter $c$ of $f_{\kappa}(x_{0},\delta)$ needs to be set within the appropriate range before a potential counterfactual instance is found. Finding a good range can be time consuming.

Dhurandhar et al. (2018) aim to address the first issue by adding in an additional loss term $L_{\text{AE}}$ which represents the $L_{2}$ reconstruction error of $x_{cf}$ evaluated by an autoencoder AE which is fit on the training set:

[TABLE]

The loss $L$ to be minimized now becomes:

[TABLE]

The autoencoder loss term $L_{\text{AE}}$ penalizes out-of-distribution counterfactual instances, but does not take the data distribution for each prediction class $i$ into account. This can lead to sparse but uninterpretable counterfactuals, as illustrated by Figure 2. The first row of Figure 2(b) shows a sparse counterfactual $3$ generated from the original $5$ using loss function (5). Both visual inspection and reconstruction of the counterfactual instance using AE in Figure 2(e) make clear however that the counterfactual lies closer to the distribution of a $5$ and is not interpretable as a $3$. The second row adds a prototype loss term to the objective function, leading to a less sparse but more interpretable counterfactual $6$.

The $L_{\text{AE}}$ loss term also does not consistently speed up the counterfactual search process since it imposes a penalty on the distance between the proposed $x_{\text{cf}}$ and its reconstruction by the autoencoder without explicitly guiding $x_{\text{cf}}$ towards an interpretable solution. We address these issues by introducing an additional loss term, $L_{\text{proto}}$.

3.2 Prototype loss term

By adding in a prototype loss term $L_{\text{proto}}$, we obtain the following objective function:

[TABLE]

where $L_{\text{AE}}$ becomes optional. The aim of $L_{\text{proto}}$ is twofold:

Guide the perturbations $\delta$ towards an interpretable counterfactual $x_{\text{cf}}$ which falls in the distribution of counterfactual class $i$. 2. 2.

Speed up the counterfactual search process without too much hyperparameter tuning.

To define the prototype for each class, we can reuse the encoder part of the autoencoder from $L_{\text{AE}}$. The encoder $\text{ENC}(x)$ projects $x\in\mathcal{X}$ onto an $E$-dimensional latent space $\mathbb{R}^{E}$. We also need a representative, unlabeled sample of the training dataset. First the predictive model is called to label the dataset with the classes predicted by the model. Then for each class $i$ we encode the instances belonging to that class and order them by increasing $L_{2}$ distance to $\text{ENC}(x_{0})$. Similar to (Snell et al., 2017), the class prototype is defined as the average encoding over the $K$ nearest instances in the latent space with the same class label:

[TABLE]

for the ordered $\{x_{k}^{i}\}_{k=1}^{K}$ in class $i$. It is important to note that the prototype is defined in the latent space, not the original feature space.

The Euclidean distance is part of a class of distance functions called Bregman divergences. If we consider that the encoded instances belonging to class $i$ define a cluster for $i$, then $\mathrm{proto}_{i}$ equals the cluster mean. For Bregman divergences the cluster mean yields the minimal distance to the points in the cluster (Banerjee et al., 2005). Since we use the Euclidean distance to find the closest class to $x_{0}$, $\mathrm{proto}_{i}$ is a suitable class representation in the latent space. When generating a counterfactual instance for $x_{0}$, we first find the nearest prototype $\mathrm{proto}_{j}$ of class $j\neq t_{0}$ to the encoding of $x_{0}$:

[TABLE]

The prototype loss $L_{\text{proto}}$ can now be defined as:

[TABLE]

where $\text{ENC}(x_{0}+\delta)$ is the encoding of the perturbed instance. As a result, $L_{\text{proto}}$ explicitly guides the perturbations towards the nearest prototype $\mathrm{proto}_{j\neq t_{0}}$, speeding up the counterfactual search process towards the average encoding of class $j$. This leads to more interpretable counterfactuals as illustrated by the experiments. Algorithm 1 summarizes this approach.

3.3 Using k-d trees as class representations

If we do not have a trained encoder available, we can build class representations using k-d trees (Bentley, 1975). After labeling the representative training set by calling the predictive model, we can represent each class $i$ by a separate k-d tree built using the instances with class label $i$. This approach is similar to (Jiang et al., 2018) who use class specific k-d trees to measure the agreement between a classifier and a modified nearest neighbour classifier on test instances. For each k-d tree $j\neq t_{0}$, we compute the Euclidean distance between $x_{0}$ and the $k$-nearest item in the tree $x_{j,k}$. The closest $x_{j,k}$ across all classes $j\neq t_{0}$ becomes the class prototype $\mathrm{proto}_{j}$. Note that we are now working in the original feature space. The loss term $L_{\text{proto}}$ is equal to:

[TABLE]

Algorithm 2 outlines the k-d trees approach.

3.4 Categorical variables

Creating meaningful perturbations for categorical data is not straightforward as the very concept of perturbing an input feature implies some notion of rank and distance between the values a variable can take. We approach this by inferring pairwise distances between categories of a categorical variable based on either model predictions (Modified Value Distance Metric) (Cost & Salzberg, 1993) or the context provided by the other variables in the dataset (Association-Based Distance Metric) (Le & Ho, 2005). We then apply multidimensional scaling (Borg & Groenen, 2005) to project the inferred distances into one-dimensional Euclidean space, which allows us to perform perturbations in this space. After applying a perturbation in this space, we map the resulting number back to the closest category before evaluating the classifier’s prediction.

3.5 Removing $L_{\text{pred}}$

In the absence of $L_{\text{proto}}$, only $L_{\text{pred}}$ encourages the perturbed instance to predict class $i\neq t_{0}$. In the case of black box models where we only have access to the model’s prediction function, $L_{\text{pred}}$ can become a computational bottleneck. This means that for neural networks, we can no longer take advantage of automatic differentiation and need to evaluate the gradients numerically. Let us express the gradient of $L_{\text{pred}}$ with respect to the input features $x$ as follows:

[TABLE]

where $f_{\text{pred}}$ represents the model’s prediction function. The numerical gradient approximation for $f_{\text{pred}}$ with respect to input feature $k$ can be written as:

[TABLE]

where $\epsilon_{k}$ is a perturbation with the same dimension as $x$ and taking value $\epsilon$ for feature $k$ and [math] otherwise. As a result, the prediction function needs to be evaluated twice for each feature per gradient step just to compute $\tfrac{\partial f_{\text{pred}}}{\partial x_{k}}$. For a $28\times 28$ MNIST image, this translates into a batch of $28\cdot 28\cdot 2=1568$ prediction function calls. Eliminating $L_{\text{pred}}$ would therefore speed up the counterfactual search process significantly. By using the prototypes to guide the counterfactuals, we can remove $L_{\text{pred}}$ and only call the prediction function once per gradient update on the perturbed instance to check whether the predicted class $i$ of $x_{0}+\delta$ is different from $t_{0}$. This eliminates the computational bottleneck while ensuring that the perturbed instance moves towards an interpretable counterfactual $x_{\text{cf}}$ of class $i\neq t_{0}$.

3.6 FISTA optimization

Like (Dhurandhar et al., 2018), we optimize our objective function by applying a fast iterative shrinkage-thresholding algorithm (FISTA) (Beck & Teboulle, 2009) where the solution space for the output $x_{\text{cf}}=x_{0}+\delta$ is restricted to $\mathcal{X}$. The optimization algorithm iteratively updates $\delta$ with momentum for $N$ optimization steps. It also strips out the $\beta\cdot L_{1}$ regularization term from the objective function and instead shrinks perturbations $|\delta_{k}|<\beta$ for feature $k$ to [math]. The optimal counterfactual is defined as $x_{\text{cf}}=x_{0}+\delta^{n^{*}}$ where $n^{*}=\operatorname*{argmin}_{n\in{1,...,N}}\beta\cdot\|\delta^{n}\|_{1}+\|\delta^{n}\|_{2}^{2}$ and the predicted class on $x_{\text{cf}}$ is $i\neq t_{0}$.

4 Experiments

The experiments are conducted on an image and tabular dataset. The first experiment on the MNIST handwritten digit dataset (LeCun & Cortes, 2010) makes use of an autoencoder to define and construct prototypes. The second experiment uses the Breast Cancer Wisconsin (Diagnostic) dataset (Dua & Graff, 2017). The latter dataset has lower dimensionality so we find the prototypes using k-d trees. Finally, we illustrate our approach for handling categorical data on the Adult (Census) dataset (Dua & Graff, 2017).

4.1 Evaluation

The counterfactuals are evaluated on their interpretability, sparsity and speed of the search process. The sparsity is evaluated using the elastic net loss term $\text{EN}(\delta)=\beta\cdot\|\delta\|_{1}+\|\delta\|_{2}^{2}$ while the speed is measured by the time and the number of gradient updates required until a satisfactory counterfactual $x_{\text{cf}}$ is found. We define a satisfactory counterfactual as the optimal counterfactual found using FISTA for a fixed value of $c$ for which counterfactual instances exist.

In order to evaluate interpretability, we introduce two interpretability metrics IM1 and IM2. Let $\text{AE}_{i}$ and $\text{AE}_{t_{0}}$ be autoencoders trained specifically on instances of classes $i$ and $t_{0}$, respectively. Then IM1 measures the ratio between the reconstruction errors of $x_{\text{cf}}$ using $\text{AE}_{i}$ and $\text{AE}_{t_{0}}$:

[TABLE]

A lower value for IM1 means that $x_{\text{cf}}$ can be better reconstructed by the autoencoder which has only seen instances of the counterfactual class $i$ than by the autoencoder trained on the original class $t_{0}$. This implies that $x_{\text{cf}}$ lies closer to the data manifold of counterfactual class $i$ compared to $t_{0}$, which is considered to be more interpretable.

The second metric IM2 compares how similar the reconstructed counterfactual instances are when using $\text{AE}_{i}$ and an autoencoder trained on all classes, AE. We scale IM2 by the $L_{1}$ norm of $x_{\text{cf}}$ to make the metric comparable across classes:

[TABLE]

A low value of IM2 means that the reconstructed instances of $x_{\text{cf}}$ are very similar when using either $\text{AE}_{i}$ or AE. As a result, the data distribution of the counterfactual class $i$ describes $x_{\text{cf}}$ as good as the distribution over all classes. This implies that the counterfactual is interpretable. Figure 2 illustrates the intuition behind IM1 and IM2.

The uninterpretable counterfactual $3$ ($x_{\text{cf,1}}$) in the first row of Figure 2(b) has an IM1 value of 1.81 compared to 1.04 for $x_{\text{cf,2}}$ in the second row because the reconstruction of $x_{\text{cf,1}}$ by $\text{AE}_{5}$ in Figure 2(d) is better than by $\text{AE}_{3}$ in Figure 2(c). The IM2 value of $x_{\text{cf,1}}$ is higher as well—0.15 compared to 0.12 for $x_{\text{cf,2}}$)—since the reconstruction by AE in Figure 2(e) yields a clear instance of the original class $5$.

Finally, for MNIST we apply a multiple model comparison test based on the maximum mean discrepancy (Lim et al., 2019) to evaluate the relative interpretability of counterfactuals generated by each method.

4.2 Handwritten digits

The first experiment is conducted on the MNIST dataset. The experiment analyzes the impact of $L_{\text{proto}}$ on the counterfactual search process with an encoder defining the prototypes for $K$ equal to 5. We further investigate the importance of the $L_{\text{AE}}$ and $L_{\text{pred}}$ loss terms in the presence of $L_{\text{proto}}$. We evaluate and compare counterfactuals obtained by using the following loss functions:

[TABLE]

For each of the ten classes, we randomly sample 50 numbers from the test set and find counterfactual instances for 3 different random seeds per sample. This brings the total number of counterfactuals to 1,500 per loss function.

The model used to classify the digits is a convolutional neural network with 2 convolution layers, each followed by a max-pooling layer. The output of the second pooling layer is flattened and fed into a fully connected layer followed by a softmax output layer over the 10 possible classes. For objective functions $B$ to $F$, the experiment also uses a trained autoencoder for the $L_{\text{AE}}$ and $L_{\text{proto}}$ loss terms. The autoencoder has 3 convolution layers in the encoder and 3 deconvolution layers in the decoder. Full details of the classifier and autoencoder, as well as the hyperparameter values used can be found in the supplementary material.

4.2.1 Results

Table 1 summarizes the findings for the speed and interpretability measures.

4.2.2 Speed

Figure 3(a) shows the mean time and number of gradient steps required to find a satisfactory counterfactual for each objective function. We also show the standard deviations to illustrate the variability between the different loss functions. For loss function $A$, the majority of the time is spent finding a good range for $c$ to find a balance between steering the perturbed instance away from the original class $t_{0}$ and the elastic net regularization. If $c$ is too small, the $L_{1}$ regularization term cancels out the perturbations, but if $c$ is too large, $x_{\text{cf}}$ is not sparse anymore.

The aim of $L_{\text{AE}}$ in loss function $B$ is not to speed up convergence towards a counterfactual instance, but to have $x_{\text{cf}}$ respect the training data distribution. This is backed up by the experiments. The average speed improvement and reduction in the number of gradient updates compared to $A$ of respectively 36% and 54% is significant but very inconsistent given the high standard deviation. The addition of $L_{\text{proto}}$ in $C$ however drastically reduces the time and iterations needed by respectively 77% and 84% compared to $A$. The combination of $L_{\text{AE}}$ and $L_{\text{proto}}$ in $D$ improves the time to find a counterfactual instance further: $x_{\text{cf}}$ is found 82% faster compared to $A$, with the number of iterations down by 90%.

So far we have assumed access to the model architecture to take advantage of automatic differentiation during the counterfactual search process. $L_{\text{pred}}$ can however form a computational bottleneck for black box models because numerical gradient calculation results in a number of prediction function calls proportionate to the dimensionality of the input features. Consider $A^{\prime}$ the equivalent of loss function $A$ where we can only query the model’s prediction function. $E$ and $F$ remove $L_{\text{pred}}$ which results in approximately a 100x speed up of the counterfactual search process compared to $A^{\prime}$. The results can be found in the supplementary material.

4.2.3 Quantitative interpretability

IM1 peaks for loss function $A$ and improves by respectively 13% and 26% as $L_{\text{AE}}$ and $L_{\text{proto}}$ are added (Figure 3(b)). This implies that including $L_{\text{proto}}$ leads to more interpretable counterfactual instances than $L_{\text{AE}}$ which explicitly minimizes the reconstruction error using AE. Removing $L_{\text{pred}}$ in $E$ yields an improvement over $A$ of 29%. While $L_{\text{pred}}$ encourages the perturbed instance to predict a different class than $t_{0}$, it does not impose any restrictions on the data distribution of $x_{\text{cf}}$. $L_{\text{proto}}$ on the other hand implicitly encourages the perturbed instance to predict $i\neq t_{0}$ while minimizing the distance in latent space to a representative distribution of class $i$.

The picture for IM2 is similar. Adding in $L_{\text{proto}}$ brings IM2 down by 34% while the combination of $L_{\text{AE}}$ and $L_{\text{proto}}$ only reduces the metric by 24%. For large values of $K$ the prototypes are further from $\text{ENC}(x_{0})$ resulting in larger initial perturbations towards the counterfactual class. In this case, $L_{\text{AE}}$ ensures the overall distribution is respected which makes the reconstructed images of $\text{AE}_{i}$ and AE more similar and improves IM2. The impact of $K$ on IM1 and IM2 is illustrated in the supplementary material. The removal of $L_{\text{pred}}$ in $E$ and $F$ has little impact on IM2. This emphasizes that $L_{\text{proto}}$—optionally in combination with $L_{\text{AE}}$—is the dominant term with regards to interpretability.

Finally, performing kernel multiple model comparison tests (Lim et al., 2019) indicates that counterfactuals generated by methods not including the prototype term ($A$ and $B$) result in high rejection rates for faithfully modelling the predicted class distribution (see supplementary material).

4.2.4 Visual interpretability

Figure 4 shows counterfactual examples on the first row and their reconstructions using AE on the second row for different loss functions. The counterfactuals generated with $A$ or $B$ are sparse but uninterpretable and are still close to the manifold of a $2$. Including $L_{\text{proto}}$ in Figure 4(d) to (g) leads to a clear, interpretable [math] which is supported by the reconstructed counterfactuals on the second row. More examples can be found in the supplementary material.

4.2.5 Sparsity

The elastic net evaluation metric $\text{EN}(\delta)$ is also the only loss term present in $A$ besides $L_{\text{pred}}$. It is therefore not surprising that $A$ results in the most sparse counterfactuals (Figure 3(c)). The relative importance of sparsity in the objective function goes down as $L_{\text{AE}}$ and $L_{\text{proto}}$ are added. $L_{\text{proto}}$ leads to more sparse counterfactuals than $L_{\text{AE}}$ ($C$ and $E$), but this effect diminishes for large $K$.

4.3 Breast Cancer Wisconsin (Diagnostic) Dataset

The second experiment uses the Breast Cancer Wisconsin (Diagnostic) dataset which describes characteristics of cell nuclei in an image and labels them as malignant or benign. The real-valued features for the nuclei in the image are the mean, error and worst values for characteristics like the radius, texture or area of the nuclei. The dataset contains 569 instances with 30 features each. The first 550 instances are used for training, the last 19 to generate the counterfactuals. For each instance in the test set we generate 5 counterfactuals with different random seeds. Instead of an encoder we use k-d trees to find the prototypes. We evaluate and compare counterfactuals obtained by using the following loss functions:

[TABLE]

The model used to classify the instances is a 2 layer feedforward neural network with 40 neurons in each layer. More details can be found in the supplementary material.

4.3.1 Results

Table 2 summarizes the findings for the speed and interpretability measures.

4.3.2 Speed

$L_{\text{proto}}$ drastically reduces the time and iterations needed to find a satisfactory counterfactual. Loss function $B$ finds $x_{\text{cf}}$ in 13% of the time needed compared to $A$ while bringing the number of gradient updates down by 91%. Removing $L_{\text{pred}}$ and solely relying on the prototype to guide $x_{\text{cf}}$ reduces the search time by 92% and the number of iterations by 93%.

4.3.3 Quantitative interpretability

Including $L_{\text{proto}}$ in the loss function reduces IM1 and IM2 by respectively 55% and 81%. Removing $L_{\text{pred}}$ in $C$ results in similar improvements over $A$.

4.3.4 Sparsity

Loss function $A$ yields the most sparse counterfactuals. Sparsity and interpretability should however not be considered in isolation. The dataset has 10 attributes (e.g. radius or texture) with 3 values per attribute (mean, error and worst). $B$ and $C$ which include $L_{\text{proto}}$ perturb relatively more values of the same attribute than $A$ which makes intuitive sense. If for instance the worst radius increases, the mean should typically follow as well. The supplementary material supports this statement.

4.4 Adult (Census) Dataset

The Adult (Census) dataset consists of individuals described by a mixture of numerical and categorical features. The predictive task is to determine whether a person earns more than $50k/year. As the dataset contains categorical features, it is important to use a principled approach to define perturbations over these features. Figure 5 illustrates our approach using the association based distance metric (Le & Ho, 2005)(ABDM) to embed the feature “Education” into one dimensional numerical space over which perturbations can be defined. The resulting embedding defines a natural ordering of categories in agreement with common sense for this interpretable variable. By contrast, the frequency embedding method as proposed by (Dhurandhar et al., 2019) does not capture the underlying relation between categorical values.

Since ABDM infers distances from other variables by computing dissimilarity based on the K-L divergence, it can break down if there is independence between categories. In such cases one can use MVDM (Cost & Salzberg, 1993) which uses the difference between the conditional model prediction probabilities of each category. A counterfactual example changing categorical features is shown in Figure 1.

5 Discussion

In this paper we introduce a model agnostic counterfactual search process guided by class prototypes. We show that including a prototype loss term in the objective results in more interpretable counterfactual instances as measured by two novel interpretability metrics. We demonstrate that prototypes speed up the search process and remove the numerical gradient evaluation bottleneck for black box models thus making our method more appealing for practical applications. By fixing selected features to the original values during the search process we can also obtain actionable counterfactuals which describe concrete steps to take to change a model’s prediction. To facilitate the practical use of counterfactual explanations we provide an open source library with our implementation of the method (Klaise et al., ).

Appendix A Breast Cancer Wisconsin experiment details and results

The classification model used to classify the cell nuclei into the malignant or benign categories is a 2 layer feedforward neural network with 40 neurons and ReLU activations in each layer. The model is trained on standardized features with stochastic gradient descent for 500 epochs with batch size 128 and reaches 100% accuracy on the test set.

The class specific autoencoders used to evaluate IM1 and IM2 consist of 3 dense layers in the encoder with respectively 20, 10 and 6 neurons for each layer. The first 2 layers have ReLU activations whilst the last one has a linear activation. The dense layers in the decoder contain 10 and 20 neurons followed by a linear layer projecting the reconstructed instance back to the input feature space. The autoencoders are optimized with Adam and trained for 500 epochs on batches of 128 instances with the mean squared error between the original and reconstructed instance as the loss function.

Similar to the MNIST experiment, parameters $c$, $\kappa$ and $\beta$ are kept constant throughout the experiments at 1, 0 and 0.1. The results of the experiment are visualized in Figure 6.

We also study the impact of different values for hyperparameters $\theta$ and $k$ which is visualized in Figures 7 and 8. The figures show that a broad range of values for both $\theta$ and $k$ work well.

Finally we study the number of unique features changed by each method. The results are visualized in Figure 9.

All experiments were run on a Thinkpad T480 with an Intel Core i7-8550U Processor.

Appendix B MNIST experiment details

The classification model consists of 2 convolutional layers with respectively 64 and 32 $2\times 2$ filters and ReLU activations. Each convolutional layer is followed by a $2\times 2$ max-pooling layer. Dropout with fraction 30% is applied during training. The output of the second pooling layer is flattened and fed into a fully connected layer of size 256 with ReLU activation and 50% dropout. This dense layer is followed by a softmax output layer over the 10 classes. The model is trained with an Adam optimizer for 3 epochs with batch size 64 on MNIST images scaled to $[-0.5,0.5]$ and reaches a test accuracy of 98.6%.

The autoencoder used in objective functions $B$ to $F$ has 3 convolutional layers in the encoder. The first 2 contain 16 $3\times 3$ filters and ReLU activations and are followed by a $2\times 2$ max-pooling layer which feeds into a convolution layer with 1 $3\times 3$ filter and linear activation. The decoder takes the encoded instance as input and feeds it into a convolutional layer with 16 $3\times 3$ filters and ReLU activations, followed by a $2\times 2$ upsampling layer and again the same convolutional layer. The final convolutional is similar to the last layer in the encoder. All the convolutions have same padding. The autoencoder is trained with an Adam optimizer for 4 epochs with batch size 128 and uses the mean squared error between the original and reconstructed instance as the loss function.

The class specific autoencoders used to evaluate IM1 and IM2 consist of 3 convolutional layers with $3\times 3$ filters and ReLU activations in the encoder, each followed by $2\times 2$ max-pooling layers. The first one contains 16 filters while the others have 8 filters. The decoder follows the same architecture in reversed order and with upsampling instead of max-pooling. The autoencoder is trained with an Adam optimizer for 30 epochs and batch size 128.

Parameters $c$, $\kappa$, $\beta$ and $\gamma$ are kept constant throughout the experiments at 1, 0, 0.1 and 100. Both $L_{\text{AE}}$ and $L_{\text{proto}}$ are reconstruction errors, but $L_{\text{AE}}$ works on the full input feature space while $L_{\text{proto}}$ operates on the compressed latent space. $\theta$ is therefore set at 200 for loss functions $C$ and $E$, and 100 if used in combination with $L_{\text{AE}}$ in $D$ and $F$.

All experiments were run on a Thinkpad T480 with an Intel Core i7-8550U Processor.

Appendix C MNIST additional results

Table 3 shows the impact on the speed of counterfactual search if the model is only exposed as a black box requiring the calculation of numerical gradients.

Figure 10 shows the effect of the parameter $K$ on the speed, quality and sparsity of counterfactuals generated on the MNIST dataset.

Figures 12 and 13 show additional examples of counterfactual instances on the MNIST dataset generated with various objective functions.

Appendix D Categorical feature embedding details

Our goal is to embed every categorical feature into one-dimensional numerical space in which perturbations can be defined and thus counterfactual search performed. We want the embedding to reflect the underlying relation between categorical values.

We proceed in two steps. First, for each categorical variable with $k$ values we want to find a $k\times k$ matrix describing distances between the values. Second, we use multidimensional scaling to transform the pairwise distances into an embedding of feature values into one-dimensional space. We consider the Modified Value Distance Metric (MVDM) and the Association Based Distance Metric (ABDM) to find the required pairwise distance matrix for each categorical variable.

D.1 Pairwise distances

MVDM (Cost & Salzberg, 1993) takes as input the data $X$, the labels or predictions of a classifier $y$ and a parameter $\alpha$. In the following we assume that each column of $X$ corresponds either to a single numerical variable or a single categorical variable with ordinal encoding (one-hot encoded categorical variables can be trivially converted to ordinal encoding if necessary). Let $n$ be the number of classes in $y$ and let $v_{1},v_{2}$ be two categories of a categorical feature $v$. Then the distance between the two values is defined as

[TABLE]

where $c_{1}^{i}$ is the number of instances with the value $v_{1}$ that were classified into class $i$ and $c_{1}$ the total number of instances with the value $v_{1}$ (likewise for $c_{2}^{i}$ and $c_{2}$). We set $\alpha=1$. The interpretation of this metric is that we want two categories to be more similar if they appear with the same relative frequency across all prediction classes.

ABDM (Le & Ho, 2005) takes as input the data $X$. Let $x_{1},\dots,x_{m}$ denote each categorical feature in $X$. The distance between two values of a feature $x_{i}$ is defined as

[TABLE]

where $\Psi$ is a dissimilarity function between two probability distributions and $p(x_{j}|x_{i})$ is the conditional probability distribution of feature $x_{j}$ given feature $x_{i}$. Thus the distance between two values $v_{1},v_{2}$ of a categorical feature is directly proportional to the distances between the conditional probability distributions of other features given $v_{1},v_{2}$. Following Le & Ho (2005) we use the Kullback-Leibler divergence as the dissimilarity function $\Psi$. In practice, we also discretize every numerical feature in $X$, i.e. we calculate a histogram and map each original value to a bin. This allows the method to use both categorical and numerical features for inferring categorical distances.

D.2 Multidimensional scaling

After inferring the pairwise distances between categories we use multidimensional scaling to embed each category into 2-dimensional Euclidean space. We then use the norms of the embeddings as the one-dimensional numerical values of categories. We use the datapoint with the largest Frobenius norm in the embedded space as the origin. Finally we scale the numerical values for each category using either standard or min-max scaling to ensure the embedded categorical features are in the same range as pre-processed numerical features.

Appendix E Multiple model comparison test description and results on MNIST

We want to evaluate the generated counterfactual instances from each loss function on how well they model the counterfactual predicted class distribution. To do this we perform a kernel multiple model comparison test (Lim et al., 2019) using the open source library from the authors.†††https://github.com/wittawatj/model-comparison-test

The method compares $l>2$ models (for us $l=6$ corresponding to loss functions $A$ to $F$) on their relative fit to the data generating distribution (the test set of MNIST). The goal is to decide whether each candidate model is worse than the best one in the candidate list. Thus the test proceeds in two steps. In step one, a reference model is selected such that it is the model minimizing the discrepancy measure (we use maximum mean discrepancy or MMD) between samples generated from it and the data generating distribution (note that this is a random variable). In step two the actual test is performed with the null hypothesis that each model not selected as the reference model in step one is worse than the reference model given that the reference model was selected. Thus the hypothesis is conditional on the selection event of the reference model.

We run the tests on two flavours of the method presented in Lim et al. (2019). RelMulti partitions the sample into two disjoint sets so that step one and two for choosing the reference model and performing the test are done on independent sets. RelPSI uses the same sample for both steps. RelPSI is presented as an alternative for controlling the false positive rate which is lost during the RelMulti procedure (when the selection step is wrong, the test will give a lower true positive rate if using RelMulti).

To generate samples, for each loss function $A$ to $F$ we generate a counterfactual instance from real instances from the test set of MNIST and record the original predicted class and the counterfactual predicted class. The total number of samples generated per counterfactual prediction class for each loss is summarized in Table 4. The original instances used were the same for each loss function, the discrepancy between the total numbers is due to constraints on compute.

Note that because we can’t control which class the counterfactual instance will end up in some numbers are very low (e.g. it is very hard for any digit to become a counterfactual $1$). This has an impact on the sample size of the per-class tests, e.g. for the counterfactual class $1$ loss function $B$ only returned $8$ instances classified as $1$ by the model which defines the sample size for the test of class $1$ across all loss functions.

For a single test for each class we thus set the sample size to be the minimum number of instances found across all losses (highlighted in Table 4). For each loss that returned a higher number of instances we sample the same number without replacement. For the comparison with respect to real instances we sample the same number without replacement from the MNIST test set for that class.

Following Lim et al. (2019) we use the Inverse Multiquadratic kernel on 256 features extracted before the softmax layer of the CNN described in Appendix B. For RelMulti the sample is split in proportion 50:50 to perform the fitting and the test respectively.

We perform the tests 100 times on each class and following Lim et al. (2019) measure the selection and rejection rates of each loss function to gauge the quality of generated counterfactuals.

The results are shown in Figure 11. For each counterfactual predicted class we report the selection and rejection rate of each loss function across two methods—RelMulti and RelPSI. We can see that for classes with sample size $>200$, the non-zero rejection rates mostly correspond to loss functions $A$ and $B$ which do not include the prototype term.

Note that it is important to perform the test for each class of counterfactual instances separately. The reasoning is the same as for introducing the prototype loss term in the first place as we care about how close in distribution the generated counterfactual instances are to the same class of real instances. If, on the other hand, we run the test across all classes we will not get a reliable measure of this. For example, we noted that running the test across all classes resulted in small rejection rates for loss $A$. This is because the counterfactuals of loss $A$ are generated only using sparsity constraints, thus most counterfactuals resemble the original instance whilst giving a different prediction by the model (akin to adversarial examples). Thus aggregating classes would result in good test results for loss $A$.