Generating Likely Counterfactuals Using Sum-Product Networks

## Abstract

Abstract Explainability of decisions made by AI systems is driven by both recent regulation and user demand. These decisions are often explainable only post hoc , after the fact. In counterfactual explanations, one may ask what constitutes the best counterfactual explanation. Clearly, multiple criteria must be taken into account, although “distance from the sample” is a key criterion.
Recent methods that consider the plausibility of a counterfactual seem to sacrifice this original objective.
Here, we present a system that provides high-likelihood explanations that are, at the same time, close and sparse.
We show that the search for the most likely explanations satisfying many common desiderata for counterfactual explanations can be modeled using mixed-integer optimization (MIO).
In the process, we propose an MIO formulation of a Sum-Product Network (SPN) and use the SPN to estimate the likelihood of a counterfactual,
which can be of independent interest.

## 1 Introduction

There is no doubt that the use of artificial intelligence (AI) is increasing.
Consequently, a better understanding of the AI models deployed is needed, especially in high-risk scenarios [].
Understanding AI models is usually seen through the lens of trustworthy and explainable AI (XAI), which is concerned with techniques that help people understand, manage, and improve trust in AI models [, , ]. Explanations also serve an important role in debugging models to ensure that they do not rely on spurious correlations and traces of processing correlated with labels, such as timestamps.

In a *post-hoc* explanation, a vendor of an AI system provides an individual user with a personalized explanation of an individual decision made by the AI system,
improving the model’s trustworthiness [, ].
In this context, personalized explanations are often called local explanations because they explain the model’s decision locally, around a given sample, such as one person’s input.
Thus, local explanations provide information relevant to the user without revealing global information about the model,
regardless of whether the model is interpretable *a priori*.
Consider, for example, credit decision-making in financial services. The models utilized need to be interpretable *a priori*,
cf. Equal Credit Opportunity Act in the US []
and related regulation [, ]
in the European Union,
but an individual who is denied credit may still be interested in a personalized, local explanation.

A well-known example of local explanations is the counterfactual explanation (CE).
CE answers the question “How should a sample be changed to obtain a different result?” [].
In the example of credit decision-making, a denied client might ask what they should do to obtain the loan. The answer would take the form of a CE. For example, “Had you asked for half the amount of credit, your application would have been accepted”.
As illustrated, CE can be easily understood [, ].
However, their usefulness is influenced by many factors [], including validity, similarity, sparsity, diversity, actionability, causality, and plausibility.

Here, we focus on the plausibility of counterfactual explanations. In the above example, halving the credit amount might not be plausible because the purpose of the credit application might require a higher amount.
A closer CE with similar plausibility might be of greater use to the client.
At the same time, plausibility does not have a clear definition. The definition of suggests CE not being an outlier and measures it as the mean distance to data []. A Local Outlier Factor is often used [e.g., ]. Alternatively, define a “plausible region” as a convex hull of $k$ nearest neighbors of the factual [].
Another approach to the plausibility of CE uses (Conditional) Variational Auto Encoders [, , ] in likelihood estimation.
However, variational auto-encoders provide only a lower bound on likelihood and are intractable models of probability distributions [, ], which means that one cannot efficiently compute the exact likelihood or marginalize variables with polynomial-time complexity with respect to the size of the encoder [].
We, therefore, investigate how to maximize the likelihood when generating CE explanations.

##### Our Contribution

We propose a method for Likely Counterfactual Explanations (LiCE) using Sum-Product Networks (SPN) [] to estimate the likelihood of a counterfactual, thus plausibility, while satisfying most other common desiderata modeled within mixed-integer optimization (MIO), see Table 1 for comparison to existing methods.

MIO has been used to search for CEs before the term CE was first used []. MIO has the favorable properties of being model-agnostic (as long as the model can be formulated in MIO) and satisfying validity by design. [] considered generating multiple explanations using a custom MIO solver. One can also implement various constraints [e.g. ].
At the same time, this improves upon the early uses of tractable probabilistic models in CE generation, such as [], which estimate a Gaussian Mixture Model and which we refer to as PlaCE (Plausible CEs),
where PlaCE is limited in its inability
to handle categorical features and non-linear classifiers.

**Table 1: *Method comparison.* A check mark indicates that a given method claims to possess the given feature. The star symbol (*) means that the method is model-agnostic as long as the classifier can be expressed using MIO. Complex data means data with continuous, categorical, ordinal, and discrete contiguous values. Exogenous property holds if a method can generate unseen data samples as CEs.**
|  | LiCE | PROPLACE [] | C-CHVAE [] | FACE [] | DiCE [] | PlaCE [] | DACE [] |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Causality | ✓ |  |  |  |  |  |  |
| Sparsity | ✓ |  |  |  | ✓ | ✓ |  |
| Plausibility | ✓ | ✓ | ✓ | ✓ |  | ✓ | ✓ |
| Model-agnostic | ✓* |  | ✓ | ✓ |  |  | ✓* |
| Complex data | ✓ |  |  | ✓ |  |  |  |
| Exogenous | ✓ | ✓ | ✓ |  | ✓ | ✓ | ✓ |

The advantage can be illustrated on the German Credit dataset [].
See Figure 1, where CEs produced by a number of methods considering the diversity or plausibility of CE are compared against the factual (white cross) in the plane, where the horizontal axis represents the amount of credit and where the vertical axis is the duration.
The heatmap corresponds to the likelihood of the CE, evaluated using an SPN trained on the randomly split training set (80% of the dataset).
For example,
C-CHVAE []and FACE [] suggest approximately halving the credit amount.
The most plausible explanation produced by DiCE []
suggests decreasing the credit amount by almost a third while reducing the Duration to a single year.
VAE and PROPLACE [] suggest decreasing the credit amount even further to below 3000.
In contrast, MIO finds a counterfactual with the sought credit amount and suggests decreasing the Duration by two months. Because the visualization is a 2-dimensional projection, some changes are not visualized. LiCE, however, changes only one “hidden” feature (Installment rate as a percentage of disposable income). To compare, all other methods change at least six features (except MIO, which changes two), showing poor sparsity.

Figure: Figure 1: The heatmap shows the marginalized likelihood distribution of the German Credit dataset into a 2-dimensional space of Credit amount and Duration features. The factual (white cross) and CEs are also projected to the two dimensions. The factual is classified as being denied. Most CE methods choose distant points, sometimes with poor likelihood. The proposed method (LiCE) strikes a balance between likelihood and proximity. For methods capable of outputting multiple CEs (DiCE, VAE, MIO, LiCE), we show the most likely out of 10, based on the SPN.
Refer to caption: extracted/5871844/figures/motivation.png

In particular, we propose:

- •
Sum-Product Network encoded into an MIO formulation.
- •
Sum-Product Network as a measure of plausibility of CE, which allows the integration of the condition directly into the MIO formulation.
- •
LiCE method for the generation of CEs. An MIO model that can be constrained by or optimized with respect to the most common desiderata regarding CE generation.

##### Notation used

Throughout the paper, we consider a classification problem where the dataset $\mathcal{D}$ is a set of 2-tuples $(\mathbf{x},y)\in\mathcal{D}$. Each input vector $\mathbf{x}\in\mathcal{X}\subseteq{}^{P}$ consists of $P$ features and is taken from the input space $\mathcal{X}$ that can be smaller than $P$-dimensional real space (e.g., can contain categorical values). $x_{j}$ is the value of the $j$-th feature of the sample $\mathbf{x}$. We have $C$ classes and describe the set of classes $[C]=\{1,\ldots,C\}$. $y\in[C]$ is the true class of the sample $\mathbf{x}$.
Finally, we have a classifier $h(\mathbf{x})=\hat{y}\in[C]$ that predicts the class $\hat{y}$ for the sample $\mathbf{x}$. Our goal is to find a set of $M$ counterfactuals $\mathcal{C_{\mathbf{x}}}=\{\mathbf{x}^{\prime}_{(1)},\ldots,\mathbf{x}^{\prime
}_{(M)}\}$ for a given $\mathbf{x}$. Note that the bold $\mathbf{x}^{\prime}$ means that the subscripts refer to indices of counterfactual instances within $\mathcal{C_{\mathbf{x}}}$ rather than a feature index. More details on the notation are in Appendix A.1.

## 2 Related Work

There is a plethora of work on the search for CEs, as recently surveyed [, , , e.g.].
See [, cf. Table 1] for a one-page overview.
In our overview, we will focus on methods
that support both categorical and continuous-valued features,
with objectives related to the plausibility, actionability, or likelihood of the CE
and methods utilizing MIO.

Historically, a number of terms have been used for related approaches.
Focusing on approaches that use MIO, pioneering work [] considered classifiers based on additive tree models
and extracted an optimal plan to change a given input to a desired class at a minimum cost, calling the problem
“optimal action extraction” (OAE).
In parallel, similar approaches have been developed under the banner of “actionable recourse” [] or “algorithmic recourse” [, ].
Developing upon this, [] distinguish between contrasting explanations and consequential explanations, where actions are modeled explicitly
in a causal model.
We use the term counterfactual explanations (CEs), popularized by [, e.g.].

Given the many CE methods that generate various CEs, one must deal with the disagreement problem [], where a user could be misled by the explainer who selects the CEs that align with their interests. We argue for the use of exact methods (e.g., MIO) that return the same optimal CEs each time.

### 2.1 Counterfactuals

We define a counterfactual explanation in accordance with previous works as $\mathbf{x}^{\prime}\in\mathcal{X}$ such that $h(\mathbf{x})\neq h(\mathbf{x}^{\prime})$ and the distance between $\mathbf{x}^{\prime}$ and $\mathbf{x}$ is in some sense minimal [, ]. We refer to $\mathbf{x}$ as factual and $\mathbf{x}^{\prime}$ as counterfactual.
As mentioned above, there are many desiderata regarding the properties of CEs. Following [], the common desiderata include:

- •
*Validity.* $\mathbf{x}^{\prime}$ should be classified differently than $\mathbf{x}$
- •
*Similarity.* $\mathbf{x}^{\prime}$ should be similar (close) to $\mathbf{x}$
- •
*Sparsity.* $\mathbf{x}^{\prime}$ should change only a few features compared to $\mathbf{x}$, i.e., minimize $\lVert\mathbf{x}^{\prime}-\mathbf{x}\rVert_{0}$
- •
*Diversity.* Each $\mathbf{x}^{\prime}_{(m)}\in\mathcal{C}_{\mathbf{x}}$ should be as different as possible from any other CE in the set, ideally by proposing changes in different features. For example, one CE recommends increasing the income; another one should recommend decreasing the loan amount instead. An important example of a CE library aiming for diversity is DiCE [].
In MIO, this is usually achieved by adding constraints and resolving the formulation [, ].
- •
*Actionability.* A counterfactual should change only features that can be influenced in reality. It should also change them correctly, such as only increasing the age.
- •
*Causality.* Given that we know some causal relationships between the features, the generated CEs should follow them. For example, if $\mathbf{x}^{\prime}$ contains a decrease in the total loan amount, the number of payments or their amount should also decrease.
- •
*Plausibility.*
While computing the probability of a CE in a given probabilistic model may be NP-Hard,
CEs should not be outliers in the dataset $\mathcal{D}$ [, ].
A natural to measure normality using the likelihood of $\mathbf{x}^{\prime}$, which can be estimated for example using autoencoders [, , e.g.] (providing just lower bound), by
Gaussian Mixture Models [], or by probabilistic circuits, e.g., Sum-Product Networks (SPNs), used in this paper,
and promote CEs with the highest estimated likelihood. The probabilistic models usually scale better than nearest-neighbour based methods like the Local Outlier Factor used in [].

### 2.2 Mixed-Integer Optimiation

Mixed-integer optimization (MIO, []) is a powerful framework for modeling and solving optimization problems, where some decision variables take values from a discrete set
while others are continuously valued. Non-trivially, the problem is in NP [] and is NP-Hard, in general.
There has been fascinating progress in the field in the past half-century [].
State-of-the-art solvers based on the branch-and-bound-and-cut approach can often find global, certified optima for instances with millions of binary variables
within hours, while there are pathological instances on under a thousand variables whose global optima are still unknown.
Naturally, MIO is widely used in areas of machine learning, where both discrete and continuous decision variables need to be optimized jointly [e.g., ].

In counterfactual reasoning, where we often need to consider both categorial and continuous-valued features,
MIO has been used since the early days [, , , , , , , ].
A crucial advance has been ’s mixed polytope formulation [], which neatly combines categorical and continuous values. A feature $j$ takes a continuous value from the $[L_{j},U_{j}]$ range or one of $K_{j}$ distinct categorical values. This is useful for modeling data with missing features, especially when there is a description of why the value is missing [].
To model the mixed polytope [] of a counterfactual for the feature $j$, we create a one-hot encoding for $K_{j}$ discrete values into binary variables $d_{j,k}$ and a continuous variable $c_{j}$ with a binary indicator variable $d_{j}^{\mathrm{~{}cont}}$ equal to 1 when the feature takes a continuous value. In summary:

$$ $\displaystyle\sum_{k=1}^{K_{j}}{d_{j,k}}+d_{j}^{\mathrm{~{}cont}}=1$ (1) $\displaystyle c_{j}=F_{j}d_{j}^{\mathrm{~{}cont}}-l_{j}+u_{j}$ (2) $\displaystyle 0\leq l_{j}\leq(F_{j}-L_{j})d_{j}^{\mathrm{~{}cont}}$ (3) $\displaystyle 0\leq u_{j}\leq(U_{j}-F_{j})d_{j}^{\mathrm{~{}cont}}$ (4) $\displaystyle d_{j}^{\mathrm{~{}cont}},d_{j,k}\in\{0,1\}\qquad\quad\forall k \in[K_{j}],$ (5) $$

where $F_{j}$ is either the original value $x_{j}$ or the median value if the factual $x_{j}$ has a categorical value instead. The computation of counterfactual $c_{j}$ in (2) uses two non-negative variables, $l_{j}$ and $u_{j}$, representing the decrease and increase in the continuous value, respectively. This construction facilitates the computation of the absolute difference from the factual. Since we minimize their (weighted) sum, at least one of them will always equal 0. []

### 2.3 Sum-Product Networks

Probabilistic circuits (PCs) [] are tractable probabilistic models (or rather, computational graphs)
that support exact probabilistic inference and marginalization in time linear w.r.t. their representation size. Probabilistic circuits are defined by a tuple $(\mathcal{G},\psi,\theta),$ where $\mathcal{G}=(\mathcal{V},\mathcal{E})$ is a Directed Acyclic Graph (DAG) defining the computation model, a scope-function $\psi:\mathcal{V}\rightarrow 2^{[P]}$ defines a subset of features over which the node defines its distribution,
and a set of parameters $\theta$. To simplify the notation,
we define a function $\mathrm{pred}({n})$, giving a set of children (predecessors) of an inner node $n$.

An important subclass of PCs are Sum-Product Networks (SPNs), which restrict PCs such that the inner (non-leaf) nodes can be only sums and products. We denote the set of sum nodes $\mathcal{V}^{\mathrm{\Sigma}}$ and product nodes $\mathcal{V}^{\mathrm{\Pi}}$.

Leaf node $n^{\mathrm{L}}\in\mathcal{V}^{\mathrm{L}}=\{n\,|\,\mathrm{pred}({n})=\emptyset\}$ within SPNs takes a value $O_{n^{\mathrm{L}}}$ from a (tractable) distribution over its scope $\psi(n^{\mathrm{L}})$ parametrized by $\theta_{n^{\mathrm{L}}}$.

Product node $n^{\mathrm{\Pi}}\in\mathcal{V}^{\mathrm{\Pi}}$ performs a product of probability distributions defined by its children

$$ $O_{n^{\mathrm{\Pi}}}(\mathbf{x})=\prod_{a\in\mathrm{pred}({n^{\mathrm{\Pi}}})} O_{a}(x_{\psi(a)})$ (6) $$

. The scope of product nodes must satisfy decomposability,
meaning that the scopes of its children are disjoint, i.e., $\bigcap_{a\in\mathrm{pred}({n^{\mathrm{\Pi}}}){}}\psi(a)=\emptyset$, but complete $\bigcup_{a\in\mathrm{pred}({n^{\mathrm{\Pi}}})}\psi(a)=\psi(n^{\mathrm{\Pi}}).$

Sum node $n^{\mathrm{\Sigma}}\in\mathcal{V}^{\mathrm{\Sigma}}$ has its value defined as

$$ $O_{n^{\mathrm{\Sigma}}}(\mathbf{x})=\sum_{a\in\mathrm{pred}({n^{\mathrm{\Sigma }}})}w_{a,n^{\mathrm{\Sigma}}}\cdot O_{a}(\mathbf{x}),$ (7) $$

where weights $w_{a,n^{\mathrm{\Sigma}}}\geq 0$ and
$\sum_{a\in\mathrm{pred}({n^{\mathrm{\Sigma}}})}w_{a,n^{\mathrm{\Sigma}}}=1$.
The value of a sum node is thus a mixture of distributions defined by its children. The scope of each sum node must satisfy completeness (smoothness), i.e., it must hold that $\psi(a_{1})=\psi(a_{2})\;\forall a_{1},a_{2}\in\mathrm{pred}({n^{\mathrm{
\Sigma}}})$.

The root node $n^{\mathrm{root}}$ (a node without parents) has the scope-function over all features, i.e., $\psi(n^{\mathrm{root}})=[P]$.
There are a variety of methods for training SPN [], out of which we use a variant of LearnSPN [] with k-means clustering, as implemented in the SPFlow library [].

## 3 SPN encoding into MIO

To utilize the computation of likelihood provided by SPNs to ensure plausible counterfactuals generated using MIO, we formulate them using integer linear constraints. Specifically, we propose an MIO formulation for a bounded approximation of an SPN.
Since the multiplication of variables in the product nodes is not desirable for linear formulations and numerical stability is an issue for MIO solvers, we perform all computations in the log space.

##### Leaf nodes

Leaves are represented by probability distributions over a single feature. These, especially in the case of histograms, can be expressed as (or approximated by) piece-wise linear functions. The modeling of piece-wise linear within MIO is well-studied [e.g., ], and various implementations are available out-of-the-box in the Pyomo modeling library [].
We suggest a possible alternative formulation of a histogram in Appendix A.2.1. We don’t need to map categorical values to histograms; we can use the $d_{j,k}$ variables directly by computing a dot product between probability values and the one-hot encoding.

##### Product nodes

Having obtained the output of the leaf nodes, we move to the inner nodes. Each node $n$ produces an output $o_{n}$ combining the outputs $o_{a},\;a\in\mathrm{pred}({n})$
of its predecessors. Consider now a product node $n\in\mathcal{V}^{\mathrm{\Pi}}$, with output defined as
a product of predecessor outputs. Since we consider all computations in log space, this translates to

$$ $o_{n}=\sum_{a\in\mathrm{pred}({n})}o_{a}\quad\forall n\in\mathcal{V}^{\mathrm{ \Pi}}.$ (8) $$

##### Sum nodes

So far, we produced the exact output values. However, this is no longer possible for $\mathrm{sum}$ nodes. Sum node $n\in\mathcal{V}^{\mathrm{\Sigma}}$ is defined as a weighted sum of predecessor $a$ outputs. In log space, the sum would translate to $o_{n}^{*}=\log\sum_{a\in\mathrm{pred}({n})}w_{a,n}\exp(o_{a})$, which we cannot model easily. To approximate it, we notice that $w_{a,n}\exp(o_{a})=\exp(o_{a}+\log w_{a,n})$, and that we can approximate $\log\sum\exp(z)$ by $\max z$. Specifically, let $z_{a}=o_{a}+\log w_{a,n}$ and we bound

$$ $\begin{split}\max_{a\in\mathrm{pred}({n})}z_{a}&=\log\exp(\max_{a\in\mathrm{ pred}({n})}z_{a})\\ &\leq\log\sum_{a\in\mathrm{pred}({n})}\exp(z_{a})=o_{n}^{*}\\ &\leq\log\left(|\mathrm{pred}({n})|\exp(\max_{a\in\mathrm{pred}({n})}z_{a}) \right)=\log(|\mathrm{pred}({n})|)+\max_{a\in\mathrm{pred}({n})}z_{a}.\end{split}$ $$

That is: we can bound the true $o_{n}^{*}$ from below, and the error in the estimate will be at most logarithm of the number of predecessors.
Putting this into our formulation, we once again utilize the general max constraint, linearized by the solver

$$ $o_{n}=\max_{a\in\mathrm{pred}({n})}(o_{a}+\log w_{a,n})\quad\forall n\in \mathcal{V}^{\mathrm{\Sigma}},$ (9) $$

where $o_{n}$ is in a conservative estimate of the true $o_{n}^{*}$ as explained above. If we wanted to use the overestimation variant, we could easily add $\log|\mathrm{pred}({n})|$ to the estimate $o_{n}$.
When this approximation of $\mathrm{sum}$ nodes by $\max$ nodes is used during training, it is called hard EM and was recommended by in the original paper introducing the SPN architecture []. Using hard EM during training, therefore, would make the approximation exact.

## 4 Our Method

We now present our model for Likely Counterfactual Explanations (LiCE). Combining the strictness and global properties of MIO with the likelihood modeling capabilities of SPNs, we create a model capable of satisfying most of the common desiderata for CEs.
We assume that all continuous values are normalized to the range $[0,1]$.

##### Input encoding

To encode the input vector, we extend the mixed polytope formulation [], as explained in Eqs. 1–5 on page 1.
The mixed polytope encoding works for purely continuous values by setting $K_{j}=0$. For purely categorical features, the original implementation contains an issue. The first categorical value (represented by zero) is mapped to the continuous variable. This seems to work fine for the logarithmic regression [], but it failed on non-monotone neural networks, leading to non-binary outputs. This was corrected by replacing the continuous variable $c_{j}$ with another binary decision variable, making it a standard one-hot encoding.

As an input to the neural network, we concatenate all variables $c_{j}$ and $d_{j,k}$ (but not $d_{j}^{\mathrm{~{}cont}}$) into a single vector. With some abuse of the notation, we denote this vector $\mathbf{x}^{\prime}$. We further use $\mathcal{X}$ for the space of encoded inputs and $P$ for the number of encoded variables
when there is no risk of confusion.

##### Model encoding

To encode the classification model, we use the OMLT library [].
It enables simple encoding of various machine-learning models, including neural networks and gradient-boosted trees. In our experiments, we used a simple neural network with ReLU activations. Linear combinations in layers are modeled directly, while ReLUs are modeled using general max constraints, which are mixed-integer representable using linear and SOS1 [] constraints [].

We now show that we can model various desiderata as constraints in the MIO formulation. Given that they are constraints, we can be certain that each generated counterfactual satisfies them. The potential of MIO to formulate such constraints is well discussed in the literature [e.g., , , , ], though they rarely provide a concrete formulation. We discuss MIO formulations of all considered desiderata in Appendix A.2.

##### Validity

Let $h^{\mathrm{raw}}:\mathcal{X}\rightarrow\mathcal{Z}$ be the neural network model $h(\cdot)$ without activation at the output layer. Let $h^{\mathrm{raw}}(\mathbf{x}^{\prime})$ be the result that we obtain from the model implementation.
Assuming that we have a binary classification task ($C=2$), a neural network typically has a single output neuron ($\mathcal{Z}=\real$). A sample $\mathbf{x}$ is classified based on whether the raw output is above or below 0, i.e., $h(\mathbf{x})=\mathds{1}\{h^{\mathrm{raw}}(\mathbf{x})\geq 0\}$.
Therefore, we can model the decision using the desired sign $\sigma\in\{-1,1\}$ as

$$ $\sigma h^{\mathrm{raw}}(\mathbf{x}^{\prime})\geq\tau,$ (10) $$

where $\tau\geq 0$ is a margin that can be set to ensure a higher certainty of the decision, improving the reliability of the CE. We present further specifications of validity for $C>2$ in Appendix A.2.2.

##### Similarity and Sparsity

To ensure similarity of the counterfactual, we follow [, ] and use
the somewhat non-standard $\lVert\cdot\rVert_{1,\textrm{MAD}}$ norm, weighed by inverse Median Absolute Deviation (MAD)

$$ $\displaystyle\lVert\mathbf{x}\rVert_{1,\textrm{MAD}}$ $\displaystyle=\sum_{j=1}^{P}{\left\lvert\frac{x_{j}{}}{\mathrm{MAD}_{j}}\right\rvert}$ (11) $\displaystyle\mathrm{MAD}_{j}$ $\displaystyle=\mathrm{median}_{(\mathbf{x},\cdot)\in\mathcal{D}}\left(|x_{j}{} -\mathrm{median}_{(\mathbf{x},\cdot)\in\mathcal{D}}(x_{j})|\right).$ $$

This metric also improves sparsity and adds scale invariance that is robust to outliers [].

##### Actionability and Causality

Both actionability and causality depend on prior knowledge of the data. For immutability, the constraint is simply $x_{j}{}=x_{j}^{\prime}$ for all immutable features $j$. We can also set the input value as a parameter instead of a variable, omitting the mixed polytope encoding.
In causality, the constraints are in the form of implications [].
We describe a way to model causal constraints where if one value changes in a certain direction, then another feature must change accordingly. This can model, among others, monotonicity, i.e., that a given value cannot decrease/increase. We call a CE *actionable* if it satisfies monotonicity and immutability constraints. Details are provided in Appendix A.2.4.

##### Diversity and Robustness

The diversity of CEs generated by MIO is discussed in the literature [, ].
Although their approach can be applied to our model too, here we simply generate a set of top-$M$ counterfactuals closest to the global optimum. We can optionally limit the maximal distance relative to the optimal CE; see the Appendix A.2.6.
Regarding the robustness of the counterfactuals,
show that finding plausible CEs indirectly improves the robustness [].
Thus, we do not add any further constraints to the model despite this being a viable option [e.g., , ].

##### Plausibility

Finally, we address plausibility. are the first to use a tractable model to approximate the likelihood of a counterfactual [].
Their formulation approximates a GMM with a quadratic term and uses a general convex optimization solver.
Instead, we propose to use a trained SPN (see Section 3) that better captures categorical data, and our formulation is linear.
For the SPN computation, we add the counterfactual class $\hat{y}^{\prime}$ as an additional input feature $x_{P+1}^{\prime}=\hat{y}^{\prime}$.
Obtaining the likelihood estimate by SPN means taking the output value of the root ($o_{\mathrm{root}}$).
We propose two variants of utilizing this value. We can add the value to the minimization objective with some multiplicative coefficient $\alpha>0$. We subtract this value since we wish to maximize likelihood. In experiments we use $\alpha=0.1$ since features are normalized to [0 0,$1$] and log-likelihood often takes values in the [$-100$, $-10$] range.

Alternatively, we set $\alpha=0$ and instead introduce a constraint forcing all generated CEs to have likelihood at or above a certain threshold $\delta^{\mathrm{SPN}}$. Such constraint is simply

$$ $o_{n^{\mathrm{root}}}\geq\delta^{\mathrm{SPN}},$ (12) $$

where $\delta^{\mathrm{SPN}}$ is a hyperparameter of our method. Similarly to [], we choose $\delta^{\mathrm{SPN}}$ as the median log-likelihood of the training data, estimated by our SPN.

##### Full LiCE model

In summary, our method optimizes the following problem:

$$ $\displaystyle\arg\min_{\mathbf{l},\mathbf{u},\mathbf{d}}\;$ $\displaystyle(\mathbf{l}+\mathbf{u})^{\mathsf{T}}\mathbf{v}^{\mathrm{cont}}+( \mathbf{d}-\mathbf{d}^{\mathrm{fact}})^{\mathsf{T}}\mathbf{v}^{\mathrm{bin}}- \alpha\cdot o_{n^{\mathrm{root}}}$ (13) s.t. mixed-polytope conditions ( 1 – 5 ) hold validity constraint (e.g., 10 ) holds SPN constraints ( 8 – 9 ) hold plausibility constraint ( 12 ) holds $\displaystyle\textrm{data-specific desiderata constraints hold},$ $$

where $\mathbf{l},\mathbf{u}$ and $\mathbf{d}$ represent the vectors obtained by concatenation of the parameters in Eqs. 1–5. The vector $\mathbf{d}^{\mathrm{fact}}$ is the vector of binary variables of the encoded factual $\mathbf{x}$. $\mathbf{v}^{\mathrm{cont}}$ and $\mathbf{v}^{\mathrm{bin}}$ represent weights for continuous and binary variables, respectively. In our case, the weights for feature $j$ are $1/\mathrm{MAD}_{j}$ and Eq. 13 corresponds Eq. 11. Details about the data-specific constraints are in Appendix A.4.1.

## 5 Experiments

We aim to compare the plausibility, similarity, and sparsity of counterfactuals generated by a number of methods. To quantify these criteria, we use a (negative) log-likelihood estimate using the SPN, $\lVert\cdot\rVert_{1,\textrm{MAD}}$, and the number of modified features, respectively.

To test the methods, we find CEs for 100 factuals. We take the factuals from both classes and look for the opposite one. For methods that can output more CEs (MIO, LiCE, DiCE, VAE), we find (at most) 10. We test for validity and actionability and measure their quality using an SPN (we use the SPFlow library []) trained on the same training set. Then, we select the CE with the highest estimated likelihood that is valid. If a method requires any prior training, we use the default hyperparameters (or some reasonable values for parameters without defaults, details in Appendix A.4.2) and train it on the same training set. If a given method can take into account actionability and causality constraints, we enforce them. The specific constraints for each dataset can be found in Appendix A.4.1.

We test finding counterfactuals for a simple neural network model consisting of 2 hidden layers with ReLU activations.

##### Data

We performed tests on the Give Me Some Credit (GMSC) dataset [], Adult dataset [] and German Credit (referred to as Credit) dataset []. We dropped some outlier data and some less informative features (details in Appendix A.3) and performed all experiments in a 5-fold cross-validation setting, training each model 5 times per dataset.

##### Methods

We compare the following CE methods:

**Table 2: The proportion of factual instances for which a given method was able to generate a *valid* (or *actionable*) counterfactual. LiCE with a fixed likelihood bound sometimes does not generate a feasible counterfactual within the 2-minute time limit. For the other methods, it is a similar mix of inability to generate a CE or to satisfy the conditions. Overall, these results show that MIO methods have a high success rate unless the constraints are too tight. Notably, once a CE is generated, it is guaranteed to be both valid and actionable.**
|  | GMSC [] | Adult [] | Credit [] |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| Method | Valid | Actionable | Valid | Actionable | Valid | Actionable |
| DiCE | $100.0$% | $100.0$% | $99.8$% | $56.2$% | $98.4$% | $3.400\,000\,000\,000\,000\,4$% |
| VAE | $1.400\,000\,000\,000\,000\,1$% | $0.2$% | $75.4$% | $10.2$% | $27.200\,000\,000\,000\,003$% | $0.0$% |
| C-CHVAE | $98.6$% | $21.6$% | $16.8$% | $8.6$% | $11.0$% | $8.799\,999\,999\,999\,999$% |
| FACE ($\epsilon$) | $98.6$% | $13.200\,000\,000\,000\,001$% | $62.0$% | $19.6$% | $27.200\,000\,000\,000\,003$% | $10.0$% |
| FACE (knn) | $98.6$% | $16.2$% | $79.4$% | $28.4$% | $27.200\,000\,000\,000\,003$% | $8.799\,999\,999\,999\,999$% |
| PROPLACE | $98.6$% | $6.600\,000\,000\,000\,000\,5$% | $79.4$% | $6.600\,000\,000\,000\,000\,5$% | $27.200\,000\,000\,000\,003$% | $11.799\,999\,999\,999\,999$% |
| MIO | $100.0$% | $100.0$% | $100.0$% | $100.0$% | $100.0$% | $100.0$% |
| LiCE (optimize) | $100.0$% | $100.0$% | $100.0$% | $100.0$% | $100.0$% | $100.0$% |
| LiCE (median) | $53.6$% | $53.6$% | $91.4$% | $91.4$% | $100.0$% | $100.0$% |

- •
MIO represents our method without the SPN model directly formulated. We use all constraints described in Section 4, including actionability and causality constraints.
- •
LiCE is our proposed model. It is similar to MIP, with the additional plausibility constraint (Eq. 12) using the SPN. We show two variants, one with bounded log-likelihood at the median log-likelihood value on training data and another where we optimize a combination of distance and likelihood with $\alpha=0.1$. We solve for up to 2 minutes, after which we recover (up to) 10 best solutions.
- •
DiCE is a well-known method focusing on generating a diverse set of counterfactuals
[].
- •
VAE is a method using VAE. We use the implementation available in version 0.4 of the DiCE library, based on [].
- •
C-CHVAE uses a Conditional VAE to search for plausible (they use the term faithful) CEs without the need for a metric in the original space [].
- •
FACE focuses on selecting a CE from the training set $\mathcal{D}$, rather than generating it from $\mathcal{X}$. It works by navigating a graph of the samples $\mathbf{x}$, where an edge exists between 2 samples if they are close ($\epsilon$ variant) or by connecting k nearest neighbors (knn variant)
[].
- •
PROPLACE is an MIO-based method for finding robust CEs within a “plausible region” constructed as a convex hull of the factual and its (robust) nearest neighbors. [].

We use the implementations of FACE and C-CHVAE provided in the CARLA library [].
MIO and LiCE are implemented using the open-source Pyomo modeling library [] that allows for simple use of (almost) any MIO solver. We use the Gurobi solver [] under an academic license. The entire implementation, together with data, is open source(^1^11Code is available at [https://github.com/Epanemu/LiCE](https://github.com/Epanemu/LiCE)). Details about hyperparameters and other configurations of experiments are in Appendix A.4.

##### Results

**Table 3: Mean negative log-likelihood (NLL), $\lVert\cdot\rVert_{1,\textrm{MAD}}$ distance, and the number of changed features, measured on *valid* generated counterfactuals for 500 testing samples from each dataset, with information about standard deviation. The log-likelihood is estimated by the SPN trained for a given fold of the 5-fold cross-validation. The number of valid counterfactuals generated by a given method varies (see Table 2), so the direct comparison between methods is non-trivial. The (+spn) means that the given method generates multiple CEs from which we choose the likeliest valid counterfactual based on the SPN. For all measures, a lower value is better.**
|  | GMSC [] | Adult [] | Credit [] |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Method | NLL $\downarrow$ | Similarity $\downarrow$ | Sparsity $\downarrow$ | NLL $\downarrow$ | Similarity $\downarrow$ | Sparsity $\downarrow$ | NLL $\downarrow$ | Similarity $\downarrow$ | Sparsity $\downarrow$ |
| DiCE (+spn) | $30.548\,357\,877\,371\,81$ $\pm$$3.744\,231\,256\,468\,958\,6$ | $20.729\,892\,468\,623\,28$ $\pm$$5.155\,386\,617\,387\,275$ | $6.378$ $\pm$$0.991\,522\,062\,286\,059\,2$ | $20.520\,055\,630\,843\,24$ $\pm$$3.020\,977\,294\,836\,486\,6$ | $23.498\,686\,273\,670\,65$ $\pm$$9.885\,384\,640\,171\,358$ | $4.605\,210\,420\,841\,684$ $\pm$$1.707\,713\,523\,037\,501$ | $51.642\,936\,129\,476\,11$ $\pm$$17.864\,706\,511\,740\,18$ | $28.139\,732\,404\,738\,78$ $\pm$$7.857\,014\,851\,632\,647$ | $8.658\,536\,585\,365\,853$ $\pm$$2.186\,727\,399\,610\,254\,5$ |
| VAE (+spn) | $23.133\,582\,347\,502\,05$ $\pm$$12.620\,280\,330\,907\,162$ | $15.360\,480\,876\,614\,137$ $\pm$$4.409\,715\,811\,918\,999$ | $9.142\,857\,142\,857\,142$ $\pm$$0.832\,993\,127\,835\,042\,9$ | $17.128\,829\,221\,598\,37$ $\pm$$2.990\,694\,423\,314\,944\,4$ | $33.317\,975\,395\,768\,76$ $\pm$$10.897\,937\,523\,366\,313$ | $5.511\,936\,339\,522\,546\,5$ $\pm$$1.453\,246\,878\,309\,409\,7$ | $48.999\,691\,246\,059\,31$ $\pm$$17.800\,686\,227\,752\,045$ | $28.836\,044\,522\,219\,06$ $\pm$$7.826\,043\,034\,825\,539$ | $10.860\,294\,117\,647\,058$ $\pm$$1.823\,751\,765\,380\,942\,5$ |
| C-CHVAE | $24.932\,351\,617\,002\,016$ $\pm$$2.440\,382\,556\,178\,696$ | $17.402\,834\,471\,526\,116$ $\pm$$4.700\,986\,135\,057\,483$ | $9.0$ $\pm$$0.0$ | $17.915\,342\,664\,707\,63$ $\pm$$3.233\,913\,972\,340\,564\,3$ | $8.566\,660\,930\,837\,79$ $\pm$$6.301\,532\,683\,865\,802$ | $3.0$ $\pm$$1.011\,834\,731\,470\,275\,1$ | $33.981\,052\,018\,647\,73$ $\pm$$6.517\,189\,401\,668\,139$ | $13.487\,000\,933\,770\,872$ $\pm$$4.934\,575\,257\,783\,774\,5$ | $7.581\,818\,181\,818\,182$ $\pm$$1.021\,423\,411\,286\,41$ |
| FACE ($\epsilon$) | $30.039\,852\,083\,916\,834$ $\pm$$9.041\,413\,205\,917\,035$ | $14.776\,436\,466\,364\,785$ $\pm$$3.789\,359\,455\,344\,945\,4$ | $8.470\,588\,235\,294\,118$ $\pm$$1.079\,801\,207\,581\,305\,7$ | $16.084\,158\,891\,829\,894$ $\pm$$3.021\,999\,820\,160\,709\,7$ | $14.612\,012\,040\,621\,806$ $\pm$$8.435\,818\,793\,669\,833$ | $3.722\,580\,645\,161\,290\,4$ $\pm$$1.226\,247\,798\,306\,360\,7$ | $46.376\,131\,908\,440\,3$ $\pm$$17.536\,288\,786\,253\,447$ | $18.059\,275\,679\,867\,06$ $\pm$$6.197\,235\,913\,596\,171$ | $6.963\,235\,294\,117\,647$ $\pm$$1.233\,169\,566\,576\,631$ |
| FACE (knn) | $29.541\,225\,265\,087\,704$ $\pm$$8.151\,337\,623\,332\,541$ | $14.584\,563\,027\,819\,925$ $\pm$$3.833\,834\,963\,338\,297\,5$ | $8.350\,912\,778\,904\,664$ $\pm$$1.101\,623\,446\,348\,438\,5$ | $15.626\,674\,184\,587\,388$ $\pm$$3.103\,101\,752\,671\,255$ | $14.108\,044\,390\,855\,772$ $\pm$$8.003\,764\,882\,445\,008$ | $3.687\,657\,430\,730\,478\,6$ $\pm$$1.220\,386\,775\,840\,469\,2$ | $44.439\,356\,830\,326\,85$ $\pm$$17.756\,711\,749\,612\,613$ | $18.546\,651\,052\,835\,877$ $\pm$$6.226\,349\,605\,762\,76$ | $7.117\,647\,058\,823\,529$ $\pm$$1.306\,425\,430\,265\,459\,4$ |
| PROPLACE | $27.318\,085\,896\,373\,727$ $\pm$$5.677\,743\,674\,768\,005$ | $12.137\,345\,451\,690\,734$ $\pm$$3.070\,203\,819\,539\,59$ | $7.367\,139\,959\,432\,048$ $\pm$$1.087\,031\,862\,675\,703\,4$ | $17.608\,996\,451\,752\,688$ $\pm$$3.335\,729\,262\,575\,719\,6$ | $22.040\,165\,945\,059\,343$ $\pm$$7.985\,934\,399\,464\,67$ | $4.848\,866\,498\,740\,554$ $\pm$$1.246\,589\,291\,334\,805\,6$ | $41.182\,498\,195\,113\,37$ $\pm$$17.016\,016\,655\,122\,04$ | $24.544\,981\,879\,687\,196$ $\pm$$6.283\,479\,472\,851\,551$ | $8.705\,882\,352\,941\,176$ $\pm$$1.295\,453\,855\,620\,308\,3$ |
| MIO (+spn) | $27.155\,729\,933\,087\,436$ $\pm$ $5.870\,217\,430\,532\,42$ | $6.112\,013\,439\,812\,737$ $\pm$ $1.773\,348\,813\,553\,345\,3$ | $2.654$ $\pm$ $0.879\,934\,088\,440\,719\,5$ | $18.016\,831\,844\,772$ $\pm$ $3.752\,305\,083\,335\,073\,3$ | $5.670\,272\,054\,723\,149$ $\pm$ $3.628\,840\,478\,792\,115$ | $2.164$ $\pm$ $0.934\,400\,342\,465\,690\,7$ | $47.541\,812\,289\,555\,13$ $\pm$ $18.212\,031\,751\,254\,422$ | $4.440\,776\,499\,678\,725$ $\pm$ $2.693\,063\,367\,539\,072$ | $2.252$ $\pm$ $1.045\,225\,334\,557\,099\,2$ |
| LiCE (opt) | $25.285\,510\,136\,463\,927$ $\pm$ $4.133\,463\,045\,002\,99$ | $6.937\,643\,384\,208\,569$ $\pm$ $2.485\,880\,918\,958\,983\,4$ | $3.026$ $\pm$ $1.041\,788\,846\,167\,974\,7$ | $18.102\,806\,617\,672\,538$ $\pm$ $3.769\,154\,829\,556\,925$ | $5.545\,350\,398\,862\,173$ $\pm$ $3.639\,405\,244\,177\,759$ | $2.022$ $\pm$ $0.966\,186\,317\,435\,721$ | $30.257\,270\,846\,631\,485$ $\pm$ $3.770\,407\,420\,093\,172\,3$ | $6.665\,132\,151\,715\,382$ $\pm$ $3.588\,694\,969\,754\,312\,2$ | $2.792$ $\pm$ $1.035\,729\,694\,466\,659\,6$ |
| LiCE (med) | $18.040\,588\,642\,162\,838$ $\pm$ $1.848\,664\,298\,701\,172\,8$ | $11.337\,135\,935\,754\,526$ $\pm$ $3.966\,806\,822\,095\,020\,6$ | $4.731\,343\,283\,582\,089$ $\pm$ $1.150\,367\,748\,497\,517\,4$ | $12.879\,907\,003\,114\,257$ $\pm$ $0.976\,996\,425\,856\,518$ | $9.717\,846\,465\,736\,828$ $\pm$ $6.541\,191\,462\,385\,036$ | $2.943\,107\,221\,006\,564\,5$ $\pm$ $1.344\,834\,340\,823\,649\,4$ | $30.502\,830\,993\,516\,877$ $\pm$ $3.246\,978\,526\,832\,697\,4$ | $6.795\,640\,167\,430\,78$ $\pm$ $3.853\,215\,969\,802\,196$ | $2.758$ $\pm$ $1.098\,833\,927\,397\,584\,7$ |

The comparison of the methods is non-trivial since factuals for which a given method generated a valid counterfactual are not the same for all methods. See Table 2 for details on the success rate of the presented methods. The lower rates are caused partly by a failure to create a counterfactual candidate at all or a failure (numerical or otherwise) to follow the validity/actionability criteria. Looking at Table 3, we notice that the proposed methods, LiCE and the two-stage MIO (generate CEs with MIO, then select using SPN), perform best not only with regards to Likelihood, but also in terms of similarity and sparsity, measured by the $\lVert\cdot\rVert_{1,\textrm{MAD}}$ and number of changed features, respectively. MIO-based methods are also able to generate CEs for more factuals, that is, even the more difficult ones.

DiCE, focusing on diversity, generates the least performant counterfactuals in terms of likelihood and similarity despite the use of SPN to select the most likely of the 10 CEs.
VAE and C-CHVAE both succeed in likelihood, although this is contrasted by lower success rates, bordering with unusability.
FACE variants have reasonable success rates for valid CEs as well as likelihood values. This is expected, since they return samples directly from the training data. Their similarity and sparsity, however, leave much to be desired.
Lastly, PROPLACE seems to perform quite well, lagging behind LiCE and MIO in success rates (esp. lacking in actionability), but also in similarity and sparsity measures.
Further comparisons and discussion of the results are in Appendices A.5 and A.6.

##### Limitations

Our method shares the limitations of all MIO methods. The additional SPN encoding leads to some computational overhead, especially when using the likelihood threshold, as exemplified in the results on the GMSC dataset in Table 2.

## 6 Discussion and Conclusions

We have presented a comprehensive method for generating counterfactual explanations called LiCE. In Section 4, we show that our method satisfies the most common desiderata–namely validity, similarity, sparsity, actionability, causality, and most importantly, plausibility.

Our method shows promising performance at the intersection of plausibility, similarity, and sparsity. It also reliably generates high-quality, valid, and actionable CEs. However, time concerns are relevant once the full SPN is formulated within the model.

Last but not least, the Sum-Product Network encoding to MIO opens up new avenues for research.

### Acknowledgement

This work has received funding from the European Union’s Horizon Europe research and innovation program under grant agreement No. 101070568.

## Appendix A Appendix

### A.1 Notation

Generally, notation follows these rules:

- •
Capital letters typically refer to amounts of something, as in classes, features, bins, etc. Exceptions are $U,L$, and $F$, which are taken from the original work [].
- •
Caligraphic capital letters denote sets or continuous spaces.
- •
Small Latin letters are used as indices, variables, or parameters of the MIP formulation.
- •
Small Greek letters refer to hyperparameters of the LiCE formulation or parameters of the SPN (scope $\psi$, parameters $\theta$).
- •
Subscript is used to specify the position of a scalar value in a matrix or a vector. When in parentheses, it specifies the index of a vector within a set.
- •
Superscript letters refer to a specification of a symbol with otherwise intuitively similar meaning. Except for ^P, where $P$ has the standard meaning of $P$-dimensional.
- •
A hat ($\hat{~{}}$) symbol above a variable means that the variable is the output of the Neural Network $h(\cdot)$.
- •
A prime (^′) symbol as a superscript of a variable means that the variable is a part of (or the output of) the counterfactual.
- •
In bold font are only vectors. When we work with a scalar value, the symbol is in regular font.

The specific meanings of symbols used in the article are shown in Tables 4 to 8. The symbols are divided into groups.

- •
Functions non-specific to our task (Table 4)
- •
Used indices (Table 5)
- •
LiCE (hyper)parameters that can be tuned (Table 6)
- •
Classification task and SPN symbols (Table 7)
- •
MIO formulation parameters and variables (Table 8)

**Table 4: General functions used**
| General function symbols |  |
| --- | --- |
| $|\cdot|$ | Absolute value (if scalar) or size of the set |
| $[\cdot]$ | Set of integers, $[N]=\{1,2,\ldots,N\}$ |
| $\mathds{1}\{\cdot\}$ | Equal 1 if input is true, 0 otherwise |
| $\|\cdot\|_{0}$ | $\ell_{0}$ norm, number of non-zero elements |
| $2^{[P]}$ | Set of all subsets of $[P]$ |

**Table 5: Symbols used as indices**
| Indices |  |
| --- | --- |
| $j$ | Index of features, typically $j\in[P]$ |
| $(m)$ | Index of counterfactuals within a set $\mathcal{C}_{\mathbf{x}}$, typically $m\in[M]$ |
| $n$ | A node of the SPN, $n\in\mathcal{V}$ |
| $i$ | Index of bins of a histogram in a leaf node ($n$), typically $i\in[B_{n}]$ |
| $a$ | A predecessor node (of node $n$) in the SPN, usually $a\in\mathrm{pred}({n})$ |
| $k$ | A class ($k\in[C]$) or categorical value ($k\in[K_{j}]$) index |
| $e$ | Index of the feature that is changed as an effect of causal relation R |

**Table 6: Input parameters into the LiCE formulation**
| LiCE (hyper)parameters |  |
| --- | --- |
| $\sigma$ | Desired sign of counterfactual for binary classification. $\sigma=-1$ corresponds to class 0 and $\sigma=1$ to class 1. |
| $\tau$ | The minimal difference between counterfactual class ($h^{\mathrm{raw}}(\mathbf{x}^{\prime})_{\hat{y}^{\prime}}$) and factual class ($h^{\mathrm{raw}}(\mathbf{x}^{\prime})_{\hat{y}}$) NN output value. Alternatively, for binary classification, it is the requirement for a minimal absolute value of the NN output before sigmoid activation ($h^{\mathrm{raw}}(\mathbf{x}^{\prime})$). |
| $\rho$ | Limit for the relative difference of values of the objective function within the set of closest counterfactuals $\mathcal{C}_{\mathbf{x}}$. |
| $\alpha$ | Weight of negative log-likelihood in the objective function |
| $\epsilon_{j}$ | Minimal change in continuous value $c_{j}$ of $j$-th feature. The absolute difference between $x_{j}^{\prime}$ and $x_{j}$ is either 0, or at least $\epsilon_{j}$. |
| $\delta^{\mathrm{SPN}}$ | Lower bound on the estimated value of likelihood of the generated counterfactual. |

**Table 7: Symbols of the classification task, CE search, and SPNs**
| Classification task symbols |  |
| --- | --- |
| $P$ | Number of features |
| $C$ | Number of classes |
| $\mathcal{D}$ | The dataset, set of 2-tuples $(\mathbf{x},y)\in\mathcal{D}$ |
| $\mathcal{X}$ | Input space $\mathcal{X}\subseteq{}^{P}$ |
| $\mathbf{x}$ | A (factual) sample $\mathbf{x}\in\mathcal{X}$ |
| $x_{j}$ | A $j$-th feature of sample $\mathbf{x}$ |
| $y$ | Ground truth of sample $\mathbf{x}$, $y\in[C]$ |
| $h(\cdot)$ | Classifier we are explaining $h:\mathcal{X}\rightarrow[C]$ |
| $\hat{y}$ | Classifier-predicted class<br>$h(\mathbf{x})=\hat{y}\in[C]$ |
| $h^{\mathrm{raw}}(\cdot)$ | NN classifier output without activation $h^{\mathrm{raw}}:\mathcal{X}\rightarrow\mathcal{Z}$ |
| $\mathcal{Z}$ | Output space of the NN classifier, without sigmoid/softmax activation |
| Counterfactual generation symbols |  |
| $\lVert\cdot\rVert_{1,\textrm{MAD}}$ | Counterfactual distance function (see Eq. 11) |
| $\mathcal{C}_{\mathbf{x}}$ | Set of generated counterfactuals for factual $\mathbf{x}$ |
| $M$ | Number of sought counterfactuals, $M\geq|\mathcal{C}_{\mathbf{x}}|$ |
| $\mathbf{x}^{\prime}$ | Counterfactual explanation of $\mathbf{x}$, $\mathbf{x}^{\prime}\in\mathcal{C}_{\mathbf{x}}$ |
| $\mathbf{x}^{\prime*}$ | Optimal (closest) counterfactual |
| $\mathbf{x}^{\prime}_{(m)}$ | $m$-th counterfactual explanation of factual $\mathbf{x}$ |
| $x_{j}^{\prime}$ | A value of $j$-th feature of the counterfactual |
| $\hat{y}^{\prime}$ | Predicted class of the counterfactual (can be a parameter of LiCE) |
| Sum Product Network symbols |  |
| $\mathcal{V}$ | Set of nodes of the SPN |
| $\mathcal{V}^{\mathrm{L}}$ | Set of leaf nodes |
| $\mathcal{V}^{\mathrm{\Sigma}}$ | Set of sum nodes |
| $\mathcal{V}^{\mathrm{\Pi}}$ | Set of product nodes |
| $\mathrm{pred}({\cdot})$ | Function returning children (predecessors) of a node |
| $\psi(\cdot)$ | Scope function mapping nodes to their input features $\psi:\mathcal{V}\rightarrow 2^{[P]}$ |
| $\theta$ | Parameters of the SPN |
| $O_{n}$ | Output value of a node $n\in\mathcal{V}$ |
| $w_{a,n}$ | Weight of output value of predecessor node $a$ in computing the value of sum node $n$. |
| $n^{\mathrm{root}}$ | Root node, its value is the value of the SPN |

**Table 8: Used variables and parameters in the MIO formulation**
| MIO formulation variables |  |
| --- | --- |
| $l_{j}$ | Decrease in continuous value of $j$-th feature. |
| $\mathbf{l}$ | Concatenated vector of all $l_{j}$. |
| $u_{j}$ | Increase in continuous value of $j$-th feature. |
| $\mathbf{u}$ | Concatenated vector of all $u_{j}$. |
| $c_{j}$ | Continuous value of $j$-th CE feature. |
| $d_{j,k}$ | 1 iff $x_{j}^{\prime}$ takes $k$-th categorical value $k\in K_{j}$. |
| $\mathbf{d}$ | All variables $d_{j,k}$<br>concatenated into a vector. |
| $d_{j}^{\mathrm{~{}cont}}$ | 1 iff $x_{j}^{\prime}$ takes continuous value $c_{j}$. |
| $h^{\mathrm{raw}}(\cdot)_{k}$ | Value of $h^{\mathrm{raw}}$, corresponding to class $k\in[C]$. |
| $g_{k}$ | 1 iff class $k\in[C]$ has higher $h^{\mathrm{raw}}$ value than the factual class. |
| $s_{j}$ | 1 iff $j$-the feature changed, i.e., $x_{j}\neq x_{j}^{\prime}$. |
| $r$ | 1 iff causal relation $R$ is activated, i.e., cause is satisfied and effect is enforced. |
| $\bar{b}_{n,i}$ | 1 iff $x_{j}^{\prime}$ does not belong to the $i$-th bin ($i\in[B_{n}]$), assuming $j$-th feature corresponds to node $n$, i.e., $\psi(n)=\{j\}$. |
| $o_{n}$ | Estimated output value of SPN node $n\in\mathcal{V}$. |
| MIO formulation parameters |  |
| $L_{j}$ | Lower bound on continuous values of $j$-th feature. In our implementation, equal to 0. |
| $U_{j}$ | Upper bound on continuous values of $j$-th feature. In our implementation, equal to 1. |
| $F_{j}$ | Default continuous value of $j$-th feature, equal to the value of the factual $x_{j}$, if it has continuous value. Otherwise equal to the median. |
| $K_{j}$ | Number of categorical values of $j$-th feature. |
| $f_{j}$ | Equal to $x_{j}$, if it has categorical value. If $x_{j}$ is continuous, $f_{j}$ is removed, and so are all constraints containing it. |
| $S$ | Maximal number of feature value changes of $\mathbf{x}^{\prime}$ compared to $\mathbf{x}$. Sparsity limit. |
| $R$ | Example causal relation: if $j$-th feature increases, $e$-th feature must decrease. |
| $B_{n}$ | Number of bins in the histogram of leaf node $n$. |
| $t_{n,i}$ | Threshold between $i-1$-th and $i$-th bin in histogram of leaf node $n$ |
| $q_{n,i}$ | Likelihood value of $i$-th bin of node $n$. |
| $\mathbf{v}^{\mathrm{bin}}$ | Vector of respective $\lVert\cdot\rVert_{1,\textrm{MAD}}$ weights for binary one-hot encodings. |
| $\mathbf{v}^{\mathrm{cont}}$ | Vector of respective $\lVert\cdot\rVert_{1,\textrm{MAD}}$ weights for continuous values. |
| $\mathbf{d}^{\mathrm{fact}}$ | One-hot encoded vector of factual categorical values corresponding to $\mathbf{d}$. |

### A.2 MIO formulations

The following MIO formulations of the desiderata are novel in that we came up with them, and, to the best of our knowledge, they were not formalized before. They are not too complex, but we formulate them for completeness.

A reasonable value of $\alpha$ is $0.1$ since features are normalized to [0 0,$1$] and log-likelihood often takes values in the [$-100$, $-10$] range.

#### A.2.1 SPN histogram formulation

In practice, the probability distribution of a leaf $n\in\mathcal{V}^{\mathrm{L}}$ trained on data is a histogram on a single feature $j$, i.e., $\psi(n)=\{j\}$. The interval of possible values of $x_{j}^{\prime}$ is split into $B_{n}$ bins, delimited by $B_{n}+1$ breakpoints
denoted $t_{n,i},\;i\in[B_{n}+1]$.

Because modeling that a variable belongs to a union of intervals is simpler than an intersection, we consider variables $\bar{b}_{n,i}$ that equal 1 if and only if the value $x_{j}^{\prime}$ does *not* belong to the interval $[t_{n,i},t_{n,i+1})$. This leads to a set of constraints

$$ $\displaystyle\bar{b}_{n,i}$ $\displaystyle\geq t_{n,i}-x_{j}^{\prime}$ $\displaystyle\forall n\in\mathcal{V}^{\mathrm{L}},\forall i\in[B_{n}]$ (14) $\displaystyle\bar{b}_{n,i}$ $\displaystyle\geq x_{j}^{\prime}+\epsilon_{j}-t_{n,i+1}$ $\displaystyle\forall n\in\mathcal{V}^{\mathrm{L}},\forall i\in[B_{n}]$ (15) $\displaystyle\sum_{i=1}^{B_{n}}\bar{b}_{n,i}$ $\displaystyle=B_{n}-1$ $\displaystyle\forall n\in\mathcal{V}^{\mathrm{L}}$ (16) $\displaystyle o_{n}$ $\displaystyle=\sum_{i=1}^{B_{n}}(1-\bar{b}_{n,i})\log{q_{n,i}}$ $\displaystyle\forall n\in\mathcal{V}^{\mathrm{L}}$ (17) $\displaystyle\bar{b}_{n,i}$ $\displaystyle\in\{0,1\}$ $\displaystyle\forall n\in\mathcal{V}^{\mathrm{L}},\forall i\in[B_{n}],$ (18) $$

where $q_{n,i}$ is the likelihood value in a bin $i$ and $o_{n}$ is the output value of the leaf node $n$. $\epsilon_{j}$ is again the minimal change in the feature $j$ and ensures that we consider an open interval on one side. We use the fact that all values $x_{j}$ (thus also $t_{n,i}$) are in the interval $[0,1]$. Eq. 14 sets $\bar{b}_{n,i}=1$ if $x_{j}^{\prime}<t_{n,i}$ and Eq. 15 sets $\bar{b}_{n,i}=1$ for values on the other side of the bin $x_{j}^{\prime}\geq t_{n,i+1}$. Eq. 16 ensures that a single bin is chosen and Eq. 17 sets the output value to the log value of the bin that $\mathbf{x}^{\prime}$ belongs to. This implementation of bin splitting is inspired by the formulation of interval splitting in piecewise function fitting [].

We assume that the bins cover the entire space, which we can ensure by adding at most 2 bins on both sides of the interval.

#### A.2.2 Validity

For $C>2$ classes, the raw output has $C$ dimensions ($\mathcal{Z}={}^{C}$), and the classifier assigns the class equal to the index of the highest value, i.e., $h(\mathbf{x})=\arg\max_{k\in[C]}h^{\mathrm{raw}}(\mathbf{x})_{k}$. Let $\hat{y}^{\prime}$ be the desired counterfactual class. The validity constraint, given that we specify the counterfactual class prior, is then

$$ $h^{\mathrm{raw}}(\mathbf{x}^{\prime})_{\hat{y}^{\prime}}-h^{\mathrm{raw}}( \mathbf{x}^{\prime})_{k}\geq\tau\quad\forall k\in[C]\setminus\{\hat{y}^{\prime }\}.$ (19) $$

Note that we can also implement a version where we do not care about the counterfactual class $\hat{y}^{\prime}$ in advance by the following

$$ $\begin{split}g_{k}=1\implies h^{\mathrm{raw}}(\mathbf{x}^{\prime})_{k}-h^{ \mathrm{raw}}(\mathbf{x}^{\prime})_{\hat{y}}&\geq\tau\quad\forall k\in[C] \setminus\{\hat{y}\}\\ g_{k}=0\implies h^{\mathrm{raw}}(\mathbf{x}^{\prime})_{k}-h^{\mathrm{raw}}( \mathbf{x}^{\prime})_{\hat{y}}&\leq\tau\quad\forall k\in[C]\setminus\{\hat{y} \}\\ \sum_{k\in[C]\setminus\{\hat{y}\}}g_{k}&\geq 1,\end{split}$ (20) $$

where $\implies$ can be seen either as an indicator constraint or as an implication
[], $g_{k}$ is equal to 1 if and only if class $k$ has a higher value than the factual class $\hat{y}$ in the raw output. The sum then ensures that at least one other class has a higher value.

A wide variety of constraints ensuring validity are possible. For example, we can ensure that the factual class has the lowest score by setting $\sum_{k\in[C]\setminus\{\hat{y}\}}g_{k}\geq C-1$, or we could enforce a custom order of classes.

#### A.2.3 Sparsity

To constrain the sparsity further, we can set an upper bound $S$ on the number of features changed

$$ $\displaystyle\sum_{j}s_{j}$ $\displaystyle\leq S$ $\displaystyle s_{j}$ $\displaystyle\geq 1-d_{j,f_{j}}$ $\displaystyle\forall j\in[P]$ $\displaystyle s_{j}$ $\displaystyle\geq d_{j,k}$ $\displaystyle\forall j\in[P],\;\forall k\in[K_{j}]\setminus\{f_{j}\}$ (21) $\displaystyle s_{j}$ $\displaystyle\geq l_{j}+u_{j}$ $\displaystyle\forall j\in[P]$ $\displaystyle s_{j}$ $\displaystyle\in\{0,1\}$ $\displaystyle\forall j\in[P],$ $$

where we use the binary value $s_{j}$ that equals 1 if the $j$-th feature changed, the $f_{j}$ is the categorical value of attribute $j$ of the factual (if applicable).

Neither LiCE nor MIO use this constraint.

#### A.2.4 Causality

Consider the following example of a causal relation $R$. If feature $j$ increases its value, another feature $e$ must decrease. For continuous ranges, this is formulated as

$$ $\begin{split}r&\geq u_{j}-l_{j}\\ l_{e}&\geq r\epsilon_{e}\\ u_{e}&\leq 1-r\\ r&\in\{0,1\},\end{split}$ (22) $$

where $\epsilon_{e}$ is a minimal change in the value of feature $e$ and $r$ equals 1 if the relation $R$ is active.
In the case when the features are ordinal, we can assume that their values are just categorical one-hot encoded variables ordered by indices and use:

$$ $\begin{split}r&\geq\sum_{k=f_{j}+1}^{K_{j}}d_{j,k}\\ r&\leq\sum_{k=1}^{f_{e}}d_{e,k}\\ r&\in\{0,1\},\end{split}$ (23) $$

where $f_{j}$ is the categorical value of the factual in feature $j$.
Naturally, one can see that we can use any combination of increasing/decreasing values in continuous and categorical feature spaces. With these formulations, we can also model monotone values, such as age or education. We simply replace the variable $r$ by 1.

#### A.2.5 Complex data

We use the umbrella term “Complex data” to refer to tabular data with non-real continuous values. This includes categorical (e.g., race), binary (e.g., migrant status), ordinal (e.g., education), and discrete contiguous (e.g., number of children) values.

For binary, we use a simple 0-1 encoding; categorical data is encoded into one-hot vectors; and discrete variables are discretized by fixing their value to an integer variable within the formulation. Since we normalize all values to the $[0,1]$ range, we use a proxy integer variable $m$:

$$ $\displaystyle(F_{j}-l_{j}+u_{j})*\mathrm{scale}_{j}+\mathrm{shift}_{j}$ $\displaystyle=m$ $\displaystyle m$ $\displaystyle\in\mathbb{Z}$ $$

For ordinal variables, we use the same encoding as categorical values, with the addition of the one-hot encoding being sorted by value rank to allow for the causality/monotonicity to be enforced.

#### A.2.6 Diversity

Instead of a single counterfactual, the solver returns (up to) $M$ counterfactuals closest to the global optimum, optionally within some distance range. This range is defined in terms of the objective function, which is the distance of a counterfactual in our case.
In other words, we search for a set $\mathcal{C}_{\mathbf{x}}=\{\mathbf{x}^{\prime}_{(1)},\ldots,\mathbf{x}^{\prime
}_{(M)}\}$ of counterfactuals that have a similar distance to the factual.

Let $\mathbf{x}^{\prime*}$ be the closest CE satisfying all other constraints; we can set a parameter $\rho$ that represents the relative distance of all CEs to the $\mathbf{x}^{\prime*}$ leading to the generation of set

$$ $\mathcal{C}_{\mathbf{x}}=\{\mathbf{x}^{\prime}\mid\lVert\mathbf{x}-\mathbf{x}^ {\prime}\rVert_{1,\textrm{MAD}}\leq(1+\rho)\cdot\lVert\mathbf{x}-\mathbf{x}^{ \prime*}\rVert_{1,\textrm{MAD}}\}.$ $$

Nevertheless, we disregard the relative distance parameter and search for the $M$ closest CEs. Later, we sift through the set $\mathcal{C}$ of top-$M$ counterfactuals, looking for the most likely CEs. Here, one could perform any filtering.

### A.3 Data modifications

We always remove samples with missing values. Optionally, we also remove some outlier data or uninformative features.

##### GMSC

We do not remove any feature in GMSC, but we keep only data with reasonable values to avoid numerical issues within MIO. The thresholds for keeping the sample are as follows

- •
MonthlyIncome $<$ 50000
- •
RevolvingUtilizationOfUnsecuredLines $<$ 1
- •
NumberOfTime30-59DaysPastDueNotWorse $<$ 10
- •
DebtRatio $<$ 2
- •
NumberOfOpenCreditLinesAndLoans $<$ 40
- •
NumberOfTimes90DaysLate $<$ 10
- •
NumberRealEstateLoansOrLines $<$ 10
- •
NumberOfTime60-89DaysPastDueNotWorse $<$ 10
- •
NumberOfDependents $<$ 10

this removes around 5.5% of data after data with missing values was removed. We could combat the same issues by taking a $\log$ of some of the features.
In our “pruned” GMSC dataset, there are 113,595 samples and 10 features, none of which are categorical, 7 are discrete contiguous, and the remaining 3 are real continuous. Further details are in the preprocessing code.

##### Adult

In the Adult dataset, we remove 5 features

- •
fnlwgt which equals the estimated number of people the data sample represents in the census, and is thus not actionable and difficult to obtain for new data, making it less useful for predictions,
- •
education-num because it can be substituted by ordinal feature education,
- •
native-country because it is again not actionable, less informative, and also heavily imbalanced,
- •
capital-gain and capital-loss because they contain few non-zero values.

It is not uncommon to remove the features we did, as some of them also have many missing values. We remove only about 2% of the data by removing samples with missing values.
We are left with 47,876 samples and 9 features, 5 of which are categorical, 1 is binary, 1 ordinal, and the remaining 2 are discrete contiguous. Further details are in the preprocessing code.

##### Credit

We do not remove any samples or features for the Credit dataset. The dataset contains 1,000 samples and 20 features, 10 of which are categorical, 2 are binary, 1 ordinal, 5 are discrete contiguous, and the remaining 2 are real continuous. Further details are in the preprocessing code.

All code used for the data preprocessing is in the supplementary material.

### A.4 Experiment setup

Here, we describe the details of our experiments.
The code used for experiments with examples is in the supplementary material and will be made available once the paper is accepted.

#### A.4.1 Additional Constraints

In addition to data type constraints described in Section A.3, we also constrain some features for immutability and causality.

##### GMSC

- •
Immutable: NumberOfDependents
- •
Monotone: age cannot decrease
- •
Causal: no constraints

##### Adult

- •
Immutable: race and sex
- •
Monotone: age cannot decrease and education cannot decrease
- •
Causal: education increases $\implies$ age increases

##### Credit

- •
Immutable: Number of people being liable to provide maintenance for, Personal status and sex, and foreign worker
- •
Monotone: Age cannot decrease
- •
Causal: Present residence since increases $\implies$ Age increases and Present employment since increases $\implies$ Age increases

#### A.4.2 Hyperparameter setup

The entire configuration can be found in the code, but we also present (most of) it here.

##### Neural Network

We compare methods on a neural network with four layers, first with a size equal to the length of the encoded input, then 20 and 10 for hidden layers, and a single neuron as output. It trained with batch size 64 for 50 epochs. We compare all methods on this neural network architecture, trained separately five times for each training set (from the five folds).

##### SPN

To create fewer nodes in the SPN (i.e., to not overtrain it), we set the min_instances_slice parameter to the number of samples divided by 20.

##### CE methods

We used default parameters for most methods. In cases when there were no default values set, we used the following:

- •
DiCE: we use the gradient method of searching for CEs.
- •
VAE: we set the size of the model to copy the predictor model. We parametrize the hinge loss with a margin of 0.1 and multiply the validity loss by 10 to promote validity. We use learning rate 1e-3 and batch size 64. We use weight decay of 1e-4 and train for 20 epochs (200 for the Credit data since the dataset is small).
- •
FACE: we only configure the fraction of the dataset used to search for the CE, increasing it to 0.5 for the Credit dataset due to its size.
- •
C-CHVAE: we set the size of the model to copy the predictor model. For the Credit dataset, we increase the number of training epochs to 50.
- •
PROPLACE: We create the retrained NN models to reflect the same architecture and train them for 15 epochs. We set up 1 instance of PROPLACE per class and set its delta by starting at 0.025 and decreasing by 0.005 until we are able to recover enough samples.
- •
LiCE + MIO: For our methods, we configure a time limit of 2 minutes for MIO solving. These are well enough for MIO, but LiCE tends to struggle with increasing likelihood requirements. We generate 10 closest CEs, not using the relative distance parameter. We set the decision margin $\tau=10^{-4}$ and we use one $\epsilon_{j}=10^{-4}$ for all features $j$ because they are normalized. We choose $\delta^{\mathrm{SPN}}$ equal to the median (or lower quartile) of likelihood on the dataset.

#### A.4.3 Computational resources

Most experiments ran on a personal laptop with 32GB of RAM and 16 CPUs AMD Ryzen 7 PRO 6850U, but since the proposed methods had undergone wider experimentation, their experiments were run on an internal cluster with assigned 32GB of RAM and 16 CPUs, some AMD EPYC 7543 and some Intel Xeon Scalable Gold 6146, based on their availability.

Regarding computational time, it is non-trivial to estimate. The time varies greatly for some methods since, for example, VAE retries generating a CE until a valid is found or a limit on tries is reached. Most methods we compared took a few hours for the 500 samples, including the method training.
The MIO method takes, on average, a few seconds to generate an optimal counterfactual, while LiCE often reaches the 2-minute time limit.

Considering the tests presented in this paper, we estimate 200 hours of real-time was spent generating them, meaning approximately 3,200 CPU hours.
If we include all preliminary testing, the compute time is estimated at around 15,000 CPU hours, though these are all inaccurate rough estimates, given that the hours were not tracked.

### A.5 Further comparisons

In this section, we would like to discuss some results that could not fit into the article’s main body.

#### A.5.1 LiCE variants

We tested multiple versions of using the SPN within LiCE. In Table 9, we show results

The results show that selecting the most likely CE out of 10 local optima given by MIO is surprisingly strong.
[SHOW WITHOUT THE SELECTION]
A two-stage setup can be quite performant. The results on similarity show that some of the MIO CEs are not globally optimal. This happens likely because the SPN in the second phase selects some of the locally optimal CEs.

**Table 9: Comparison of LiCE variants. (optimize) means that we optimize the likelihood together with the distance, with coefficient $\alpha=0.1$. (quartile) means that we constrain the CE to have the likelihood greater or equal to the lower quartile likelihood of training data. Finally, (median) is the same as (quartile), but we take the median instead of the quartile.**
|  | GMSC | Adult | Credit |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Method | NLL | Similarity | Sparsity | NLL | Similarity | Sparsity | NLL | Similarity | Sparsity |
| MIO (+spn) | $27.155\,729\,933\,087\,436$ $\pm$ $5.870\,217\,430\,532\,42$ | $6.112\,013\,439\,812\,737$ $\pm$ $1.773\,348\,813\,553\,345\,3$ | $2.654$ $\pm$ $0.879\,934\,088\,440\,719\,5$ | $18.016\,831\,844\,772$ $\pm$ $3.752\,305\,083\,335\,073\,3$ | $5.670\,272\,054\,723\,149$ $\pm$ $3.628\,840\,478\,792\,115$ | $2.164$ $\pm$ $0.934\,400\,342\,465\,690\,7$ | $47.541\,812\,289\,555\,13$ $\pm$ $18.212\,031\,751\,254\,422$ | $4.440\,776\,499\,678\,725$ $\pm$ $2.693\,063\,367\,539\,072$ | $2.252$ $\pm$ $1.045\,225\,334\,557\,099\,2$ |
| LiCE (optimize) | $25.285\,510\,136\,463\,927$ $\pm$ $4.133\,463\,045\,002\,99$ | $6.937\,643\,384\,208\,569$ $\pm$ $2.485\,880\,918\,958\,983\,4$ | $3.026$ $\pm$ $1.041\,788\,846\,167\,974\,7$ | $18.102\,806\,617\,672\,538$ $\pm$ $3.769\,154\,829\,556\,925$ | $5.545\,350\,398\,862\,173$ $\pm$ $3.639\,405\,244\,177\,759$ | $2.022$ $\pm$ $0.966\,186\,317\,435\,721$ | $30.257\,270\,846\,631\,485$ $\pm$ $3.770\,407\,420\,093\,172\,3$ | $6.665\,132\,151\,715\,382$ $\pm$ $3.588\,694\,969\,754\,312\,2$ | $2.792$ $\pm$ $1.035\,729\,694\,466\,659\,6$ |
| LiCE (quartile) | $26.303\,556\,070\,816\,796$ $\pm$ $2.917\,918\,334\,569\,642$ | $6.868\,610\,448\,936\,876\,4$ $\pm$ $2.472\,729\,061\,864\,144\,7$ | $2.458$ $\pm$ $0.899\,019\,465\,862\,669\,7$ | $18.408\,668\,095\,804\,902$ $\pm$ $3.498\,790\,186\,208\,691$ | $5.613\,073\,352\,113\,758$ $\pm$ $3.807\,060\,888\,195\,73$ | $2.038$ $\pm$ $0.996\,271\,047\,456\,464$ | $31.098\,135\,820\,088\,153$ $\pm$ $4.170\,588\,862\,049\,285$ | $6.578\,050\,644\,998\,177$ $\pm$ $3.576\,612\,541\,126\,014$ | $2.59$ $\pm$ $0.966\,385\,016\,440\,135$ |
| LiCE (median) | $18.040\,588\,642\,162\,838$ $\pm$ $1.848\,664\,298\,701\,172\,8$ | $11.337\,135\,935\,754\,526$ $\pm$ $3.966\,806\,822\,095\,020\,6$ | $4.731\,343\,283\,582\,089$ $\pm$ $1.150\,367\,748\,497\,517\,4$ | $12.879\,907\,003\,114\,257$ $\pm$ $0.976\,996\,425\,856\,518$ | $9.717\,846\,465\,736\,828$ $\pm$ $6.541\,191\,462\,385\,036$ | $2.943\,107\,221\,006\,564\,5$ $\pm$ $1.344\,834\,340\,823\,649\,4$ | $30.502\,830\,993\,516\,877$ $\pm$ $3.246\,978\,526\,832\,697\,4$ | $6.795\,640\,167\,430\,78$ $\pm$ $3.853\,215\,969\,802\,196$ | $2.758$ $\pm$ $1.098\,833\,927\,397\,584\,7$ |

#### A.5.2 Time concerns

Table 10 shows the median number of seconds required to generate (or fail to generate) a CE. We see that there are stark differences between methods and also between datasets. For methods using MIO, we constrain the maximal solution time to 120 seconds.

LiCE seems to be comparable on Adult as well as Credit datasets. Seeing that MIO seems to have better time complexity, we suggest that the main portion of the overhead is caused by solving the SPN formulation. Please note that the optimizing variant of LiCE takes long mainly to prove optimality. A (non-optimal) solution could likely be obtained even with a tighter time limit. [SAVE LOGS OF THE OPTIMIZATION and find out the numbers for this ]

**Table 10: Median time spent on the computation of a single CE. The values significantly above 120 in the LiCE computation on GMSC are partly caused by computational overhead in formulating the SPN and partly possibly by a mistake in the computation of the values. Time spent optimizing was limited to 120 seconds, as reflected in the LiCE duration on the Credit dataset.**
| Method | GMSC | Adult | Credit |
| --- | --- | --- | --- |
| DiCE (+spn) | $25.265\,180\,200\,018\,22$s | $17.043\,989\,650\,002\,914$s | $147.353\,215\,249\,997\,45$s |
| VAE (+spn) | $0.635\,165\,050\,014\,620\,6$s | $0.967\,617\,150\,003\,206\,9$s | $0.626\,853\,800\,007\,666$s |
| C-CHVAE | $0.430\,347\,100\,000\,005\,87$s | $0.629\,137\,349\,999\,979\,3$s | $0.586\,687\,150\,000\,017\,3$s |
| FACE ($\epsilon$) | $7.876\,486\,049\,999\,954$s | $7.457\,379\,199\,999\,991$s | $5.116\,963\,949\,999\,985$s |
| FACE (knn) | $5.654\,202\,399\,999\,974$s | $7.409\,264\,799\,999\,96$s | $5.278\,127\,100\,000\,006$s |
| PROPLACE | $0.379\,812\,350\,000\,065\,63$s | $0.238\,859\,649\,999\,994\,87$s | $0.203\,372\,999\,997\,824\,88$s |
| MIO | $0.252\,333\,779\,004\,402\,46$s | $0.359\,805\,232\,961\,662\,1$s | $0.410\,225\,712\,956\,162\,17$s |
| LiCE (optimize) | $122.894\,237\,936\,940\,04$s | $28.158\,922\,566\,450\,2$s | $0.868\,628\,578\,493\,371\,6$s |
| LiCE (median) | $122.273\,014\,989\,332\,3$s | $8.075\,433\,300\,458\,826$s | $0.807\,712\,459\,005\,415\,4$s |

#### A.5.3 Other plausibility metrics

We considered using Kernel Density Estimator (KDE) for the evaluation (similarly to []), but the KDE does not perform well on categorical data, so we decided against it.

### A.6 Further comments

Given the limited size of the Credit dataset, it is unsurprising to see so many failures of some methods. There is not much data for some methods to support the training. This might be behind the low success rate of computing a valid CE.

Regarding the other results, it is possible that the VAE method has been misconfigured for GMSC, returning very few results.

It is unclear why there is a drop in the success rate for Valid and Actionable results for the MIO-based methods. It could be caused by numerical inaccuracies, but mainly, it is important to note that despite a drop in success rate, the difference compared to other methods is significant.

The main disadvantage of LiCE is the time complexity of CE generation. This will be studied further, hopefully leading to improvement. We argue, however, that for some use cases, the user might be willing to wait to obtain a high-quality CE. We leave this decision of usefulness to the user.

#### A.6.1 Omitted methods

It is not in our power to test against all CE methods, which is why we looked for a selection of methods that consider the plausibility of generated CEs. Two methods were, however, not tested for the following reasons:

- •
PlaCE [] does not allow for an implementation for Neural Networks. It also cannot model categorical features well.
- •
DACE [] does not have a public implementation that would allow for Neural Networks as models. It might also struggle with the size of datasets used here since they are an order of magnitude larger, and DACE computes the Local Outlier Factor, meaning that the formulation size increases linearly with the increase in the number of samples.

## Peer Reviews

confidence: 4
contribution: 4
initial_rating: 8
presentation: 3
questions: - Would making structural assumptions (e.g., determinism) about SPNs make it easier to encode exactly? 
- Can you elaborate on the setup and significance of the VAE baseline from Mahajan et al. (2020)?
- The results mention that LICE (median) times out for some of the cases, and the appendix mentions that the time limit for each run was 2 minutes. Are these results sensitive to the time limit duration? Would they change drastically for a small increase in the time limit?
soundness: 4
strengths: - LiCE is a principled method of inferring *plausible* counterfactual explanations that satisfy several desiderata. The MIO approach is flexible enough to accommodate additional criteria.
- The MIO formulation for SPN inference is itself a significant contribution. It opens up space for work on SPN inference tasks such as finding the entire Maximum a Posterori (MAP) set possible. Existing work has focused on finding single solutions (e.g., Poon and Domingos, 2011 and Arya et al., 2024).

Arya, Shivvrat, Tahrima Rahman, and Vibhav Gogate. "Neural Network Approximators for Marginal MAP in Probabilistic Circuits." AAAI 2024.
summary: This paper deals with the task of finding counterfactual explanations (x’) for a multi-class classification model (h(x)). Specifically, the authors aim to find a counterfactual set (C) satisfying desiderata specified in Guidotti (2022). i.e. apart from being a valid counterfactual, each x’ \in C must be similar to the original example x, involve changing as few features as possible (sparse), comply with domain constraints such as monotonicity of age (actionable), and not be an outlier (plausibility). Additionally, C must be as diverse as possible and follow causal domain knowledge. The authors approach this problem from a mixed integer optimization perspective and use sum-product networks to model the data distribution (P(X)). To this effect, they translate the desiderata into constraints and solve the optimization problem using the OMLT library (Ceccon et al., 2022). While the trained SPN cannot be encoded exactly in log-space, the authors develop an approximate encoding and bound the likelihood. The authors call their proposed counterfactual generation approach Likely Counterfactual Explanations (LiCE) and define two variants of LiCE, one based on likelihood (upper) threshold of train-set median (LiCE (median)) and the other on minimizing a linear combination of distance and likelihood (LiCE (optimize)). They evaluate the two variants by comparing them against a non-SPN variant (MIO) and on prior work including DiCE (Mothilal et al., 2020), C-CHVAE (Pawelczyk et al., 2020), FACE (Poyiadzi et al., 2020) and PROPLACE (Jiang et al., 2024) on 3 financial data sets focusing on plausibility (as measured by log-likelihood), similarity (counterfactual distance), and sparsity (number of modified features). Their experiments show that both LiCE variants outperform the baselines, and while LiCE (median) excels at plausibility of generated counterfactual sets, LiCE (optimize) is the best at similarity even beating MIO.
weaknesses: The Accuracy of likelihood assessment is limited by the expressivity of SPNs. This might be at least partially resolved by using PFCs (Sidheekh et al., 2023) to improve performance on high-dimensional domains,

Sidheekh, Sahil, Kristian Kersting, and Sriraam Natarajan. "Probabilistic flow circuits: towards unified deep models for tractable probabilistic inference." UAI 2023.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: Content-related: 

- L53: What is the definition of an outlier in your context? Are you suggesting that an epsilon-ball around the factual could contain outliers due to factor interactions and their joint distribution? Please expand on this. 

- L80: Table 1 is not explained. Why are the "other desiderata" important? 

- L97-107: A definition of plausibility is missing in this section. You claim that DiCE is the most plausible—what criteria support this claim? Additionally, LiCE changes installment/disposable income while keeping amount and duration constant, suggesting a change in income. Is this more "realizable" than borrowing less? MIO changes two features (duration and income); could you provide more details on how these changes are measured?  

- Later in the paper, it becomes clearer that your definition of plausibility stems from the joint probability distribution extracted from SPNs. However, this raises concerns about whether this plausibility measure should be evaluated independently. In your experiments, you assume that this probability distribution provides the correct measure of plausibility without empirical validation. Could you justify this assumption further, or better yet, provide an evaluation of plausibility as a standalone concept? This is crucial, as your entire method relies on the assumption that SPN-based plausibility is the most accurate or appropriate for counterfactuals. 

- L184: Why do you take the median for categorical values? Could you use their actual values instead? Taking the median for categorical variables may cause all of them to start from the same point. 

- L188: Why will at least one always be 0? Please provide intuition. Also, the citation is placed after the full stop—should be corrected. 

- L237: Outliers are mentioned again without a definition. Also, plausibility is referenced but not defined earlier. An intuition would help the reader follow the motivation more clearly. 

- L260: Why do you work in log-space? Is it to simplify handling products? Please provide a rationale for this choice. 

- L332: Taking a threshold of 0 on the raw output of a neural network implies a threshold probability of 0.5 after applying a logit. Are there dependencies on this, and can it be changed? 

- L461: You select CEs based on SPN output, which establishes plausibility in your results. Does this give your method an advantage over baselines? For example, on L484, you claim "unparalleled" plausibility. 

- L471: You mention that the failure to terminate worsens the results. How pronounced is this effect? 


Minor points: 

- L82: "This work combines the tradition ..."—what tradition are you referring to? Please include a reference to clarify. 

- Fig2: The heatmap legend is missing. Are the most likely points represented by the most yellow areas, particularly near 0 amount and 12 months? 

- L252: "SPNs are a strict generalization"—this statement is missing a reference. 

- L293: Typo: "the it" 

- L376: Why do you need $v^{cont}$ and $v^{bin}$? Could you elaborate on the correspondence between these values and their respective roles in your approach? Does this correspondence offer advantages or properties? 

- L468: What do you mean by "considering the differences between methods" when referring to finding the error acceptable? 

- L497: What do you mean by "the limitations of all MIO methods"?
soundness: 2
strengths: - The use of SPNs to model plausibility in counterfactual explanations is a novel and creative approach, addressing limitations in handling mixed data that many existing methods struggle with. 

- The integration of SPNs and MIO to optimize multiple constraints simultaneously offers flexibility and the formulation of bounds for the output of SPNs that can be integrated as constraints into a MIO formulation is (as the authors mention) quite interesting on its own.
summary: The paper presents a novel method to model counterfactual explanation (CE) search that maximizes likelihood, closeness, and sparseness. The method leverages sum-product networks (SPNs) to estimate point likelihoods (or plausibility) and use mixed integer optimization (MIO) to optimize the different objectives simultaneously. The work includes empirical comparisons with state-of-the-art CE algorithms.
weaknesses: - The key concept of plausibility, derived from the output of SPNs, is assumed to be correct without evaluation. This is a critical assumption, as it underpins the entire method, and its validity should be examined and justified. The lack of empirical evaluation of this concept weakens the overall contribution. 

- The experiments appear to be biased in their evaluation of plausibility as the SPN distribution is taken as the true one, giving the proposed method an unfair advantage over the baselines 

- The clarity of the paper could be improved, particularly around the explanation of key terms and concepts like plausibility and outliers, which are not sufficiently defined when first introduced. 

- The method appears computationally expensive, which could be a practical concern in real-world applications. This issue is not directly addressed in the paper, but the experiments suggest that the optimization process may be time-consuming, particularly for more complex scenarios. The paper would benefit from a discussion of the computational cost and potential ways to mitigate it, such as pruning strategies, approximations, or discussing scalability.

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 2
questions: 1. Could the proposed method apply to other generative models measuring the likelihood, e.g., Gaussian processes?
2. For all variants of the proposed methods, LiCE (optimize) seems to perform the best. However, its implementation is somewhat ambiguous. The authors stated, "we optimize a combination of distance and likelihood with α = 0.1 and relax the plausibility constraint (Eq. 13)." Could the authors further elaborate on this statement? For instance, what is the "combination of distance" and how to "relax the plausibility constraint?" A pseudo-code description would be appreciated.
soundness: 3
strengths: - The writing is generally clear. The paper is well-organized. While no main theorem is presented in the paper, the derivation of the mixed integer program is technically sound.
- Comprehensive simulations were performed. Simulation results support the proposed MIO formulation.
- The authors have clearly stated the limitations of the proposed methods, including the computational challenges of solving MIO programs.
summary: This paper studies the problem of generating counterfactual explanations for a fixed classification function, say $f(x)$. The explanation here refers to different input $x'$ to the classifier $f(x)$ (with a factual baseline $x$) such that they result in different classification outcomes $f(x') \neq f(x)$. In general, there are many such realizations of input $x'$. This paper studies the problem of finding the most likely counterfactual explanations $x'$ satisfying a set of selection desiderata, including validity, similarity, sparsity, and actionability. The authors evaluate the likelihood of an input $x'$ by fitting a sum-product network (SPN). Using the formalism of SPNs, the selection desiderata could be translated into a series of equivalent polynomial constraints. The explanation generation is then reduced to an equivalent mixed integer optimization (MIO) problem.
weaknesses: - The SPN models might be limited. First, it can only account for discrete features, while most prediction tasks in practice involve high-dimensional, continuous feature inputs. Also, SPNs require structural constraints. It could be challenging to learn SPNs that best fit with the observed data.
- The explanation generation problem is reduced to an MIO program, which is generally NP-hard to solve. This presents challenges in generalizing the proposed methods to more complex domains.
- The reduction of SPN to MIO programs is not surprising. It has been known that constrained optimization in Bayesian networks is generally equivalent to linear integer programming. This paper would be most improved by proposing an efficient approximation algorithm to solve the MIO program while leveraging the graphical structure of the learned SPN.

---

confidence: 5
contribution: 3
initial_rating: 6
presentation: 3
questions: See weaknesses.
soundness: 3
strengths: The paper addresses an important problem and proposes a novel and interesting approach. Its strength is its ability to handle categorical features (often encountered in real-world data sets).
summary: The paper proposed a method for computing plausible counterfactual explanations by using sum-product networks.
weaknesses: Disclaimer: I have reviewed the paper before (Neurips 2024)—I hope that this does not disqualify me as a reviewer. In my opinion, the paper has improved significantly. In particular, the empirical evaluation is now done fairly, the notation and formalization are clear now, and the structure of the entire paper is now more accessible and more "scientific" compared to the previous version that I reviewed.

However, I still think that some more information on SPNs could be added: How powerful are SPNs? What types of distribution can they model? What are the assumptions? How to estimate SPN's parameters from data? Given the large number of existing methods (for computing plausible counterfactuals), this information would help practitioners select one of those existing methods. There is information given in the appendix but it might be also good to (briefly) answer some of those questions in the main text.

I understand that there is not enough space to fully introduce SPNs but I still think that it would be important to highlight what they can model and what they can not -- this would help practitioners to select an appropriate method for their particular use case. Currently, different methods exist for computing plausible counterfactuals based on Gaussian Mixture Models, Kernel Density Estimators, Neural Autoencoders, etc., and each of those methods has its advantages and disadvantages. However, I am happy to discuss this and step back if the authors or the other reviewers can convince me that this is not necessary or is somehow already in the paper. Thanks :)

Meta-review: This paper is truly in the borderline where two reviewers view this positively while two are others are a bit negative about this. All the reviewers agree that the integration of SPNs and MIO to optimize multiple constraints is quite novel. The idea of learning SPNs and then finding a counterfactual set is quite interesting. There were some concerns about the computational complexity, comparison to baselines and  better experimentation could improve the paper.

Decision: Accept (Poster)