Flexible and Robust Counterfactual Explanations with Minimal Satisfiable Perturbations

Let us consider a pretrained model $f:\mathcal{X}\to\mathcal{Y}$, where $\mathcal{X}\subseteq\mathbb{R}^{d}$ denotes the feature space and $\mathcal{Y}$ is the prediction space. For simplicity, let $\mathcal{Y}=\{0,1\}$, where $0$/$1$ denotes unfavorable/favorable prediction, respectively. Given an input instance $\mathbf{x}\in\mathcal{X}$, which is predicted to be the unfavorable outcome ($f(\mathbf{x})=0$), a counterfactual explanation (CFE) $\mathbf{c}$ is a data point that leads to a favorable prediction, i.e., $f(\mathbf{c})=1$, with minimal perturbations of $\mathbf{x}$. Formally, a counterfactual explanation method $g:f\times\mathcal{X}\to\mathcal{X}$ can be mathematically defined as follows:

where $\text{cost}(\cdot,\cdot):\mathcal{X}\times\mathcal{X}\to\mathbb{R}^{+}$ is a distance or cost metric that quantifies the efforts required in order to change from an input $\mathbf{x}$ to its CFE $\mathbf{c}$. In practice, the commonly used cost function includes $L_{1}$/MAD (Wachter et al., 2017; Russell, 2019), total log-percentile shift (Ustun et al., 2019), and $L_{2}$ norm on latent space (Pawelczyk et al., 2020a). To optimize Eqn. (1), it can be further transformed to the Lagrangian form (Wachter et al., 2017), as shown below:

where $\ell(\cdot,\cdot)$ is a differential function to measure the gap between $f(\mathbf{c})$ and the favorable prediction $1$, and $\lambda$ is a positive trade-off factor. By optimizing the above objective, a CFE method $g(f,\mathbf{x})$ returns a single CFE or a set of CFEs for an input $\mathbf{x}$.

The definition in Eqn. (1) captures the most basic form of counterfactual explanations. Usually, additional constraints are often required to ensure that the produced CFEs are useful and actionable for specific applications (Verma et al., 2020).

Motivated by the formalization in (Artelt et al., 2021), we formally define the robustness of CFEs in more general cases that include slight input perturbations and model changes. However, before presenting technical details, one critical question that needs to answer is “Do we want CFEs to remain consistent after a series of slight changes, or should they vary to reflect such changes?”. The answer depends on practical scenarios. In certain applications, one may expect CFEs to be sensitive to such tiny changes. For example, in the study of the effects of climate change on sea turtles (Blechschmidt et al., 2020), one may expect CFEs to be sensitive to temperature changes. In this paper, we assume that such trivial changes are either irrelevant or less important to the generation of CFEs. As such, we aim to produce CFEs that are robust to trivial changes, such as inputs added with random noise, and model retraining on new data from the same distribution.

Let $\hat{\mathbf{x}}$ represent a slightly perturbed sample that is close to $\mathbf{x}$, meaning $\hat{\mathbf{x}}\sim p(\mathbf{x})$, where $p(\mathbf{x})$ is the density estimation of perturbed samples that yield the same prediction as the input $\mathbf{x}$. Similarly, let $\hat{f}\in\mathcal{F}$ denote a retrained model belonging to the class $\mathcal{F}$, which consists of potential models that perform equivalently well as the original one.

A lower value indicates higher robustness. By minimizing the above expectation, we can generate robust CFEs. However, in real life, $p(\mathbf{x})$ and $\mathcal{F}$ are typically unknown. Intuitively, they can be determined based on specific changes that users desire to be robust against. For instance, one can consider adding Gaussian noise to the input, masking certain features, or retraining the models on data from the same distribution, to decide $p(\mathbf{x})$ and $\mathcal{F}$.

Here, we explain the root causes of non-robustness. The total loss in Eqn. (2) form is usually non-convex due to the non-convex decision surface of probabilistic models and other constraints. As shown in Figure 1, multiple local minima can be found, but current methods often select a single or $k$ CFEs from them. Such selected CFEs can be different in each trial. Next, we discuss several influential factors that result in non-robust CFEs.

• •

  Input perturbations. Input instances can be perturbed by adding some noise or masking random features. Due to the local sensitivity of large models, such trivial perturbations can significantly influence model predictions, leading to different counterfactual explanations (CFEs) (Slack et al., 2021).

• •

  Model updates. The predictive model $f$ in Eqn. (1) is typically retrained periodically in deployment. The updated model $f^{\prime}$ may exhibit slightly different behavior compared to the previous model and thus may have a great impact on the cost of the desired prediction.

• •

  Random factors. Heuristic search methods $g(f,\mathbf{x})$ for Eq. (2) often involve random factors, e.g., random initial points in gradient descent (Wachter et al., 2017), random samples in Growing Sphere (Lash et al., 2017; Pawelczyk et al., 2020a), and random selection in genetic algorithms (Sharma et al., 2020). Algorithm randomness often causes different solutions in non-convex situations.

• •

  Algorithmic configurations. Different configurations of $g(f,\mathbf{x})$ also affect the search process, e.g., the trade-off factor $\lambda$, the step size in gradient descent, the stop condition, the guide point (Van Looveren and Klaise, 2019), the radius of the sphere (Lash et al., 2017; Pawelczyk et al., 2020a), variational autoencoders used to approximate data distribution (Pawelczyk et al., 2020a; Joshi et al., 2019).

Let $\mathcal{N}=\{1,2,\ldots,d\}$ be the indices of $d$-dimensional input features. $\mathbf{x}^{i},i\in\mathcal{N}$ denotes the $i$-th feature value of an input instance $\mathbf{x}$, and $[\mathbf{a}^{i},\mathbf{b}^{i}]$ denotes the normal range of $i$-th feature. We assume the normal range of each feature can be acquired from domain knowledge. For example, the normal range of BMI (CDC, 2022) is $[18.5,24.9]$, and the normal heart rate (Shmerling, 2020) is between 60 to 100. For a given input $\mathbf{x}$, then we partition $\mathcal{N}$ into two disjoint sets by checking whether each feature is in the normal range or not: $\mathcal{N}_{0}$ (with size $d_{0}$) for abnormal features and $\mathcal{N}_{1}$ (with size $d_{1}$) for features within the normal range.

We further assume that undesired predictions are attributed to abnormal features. Our intuition is to bring an abnormal feature within its normal range, which would increase the probability of the desired prediction. To ensure consistency, we replace an abnormal feature with the closest endpoint of the corresponding normal range. For instance, if the BMI of an obese person is $40$, and the suggested BMI is $24.9$ (BMI normal range $[18.5,24.9]$). We aim to convert a subset of abnormal features into their respective normal ranges, to achieve lower cost and higher sparsity. Let $\mathcal{A}\subseteq\mathcal{N}_{0}$ represent the subset of abnormal features to be replaced. We employ the function $\eta(\cdot,\cdot)\in\mathbb{R}^{d}$ to calculate the values after replacement for the abnormal features, where the $i$-th element is computed as follows:

In applications where normal ranges are not available (we set $\mathcal{N}_{0}=\mathcal{N}$), we replace the input feature values with the corresponding values of a guide point or prototype that has the desired prediction. Correspondingly, given a guide point or prototype $\hat{\mathbf{x}}$, $\eta(\cdot,\cdot)^{i}$ can be expressed as follows:

We use the binary vector $\mathbf{m}(\mathcal{A})\in\{0,1\}^{d}$ (abbreviated as $\mathbf{m}$) to represent the presence or absence of each feature in subset $\mathcal{A}$. Specifically, each element $\mathbf{m}^{i}$ indicates whether the $i$-th feature is included in $\mathcal{A}$, i.e., $\mathbf{m}^{i}=\mathds{1}_{i\in\mathcal{A}}$. The feature replacement operation for an input $\mathbf{x}$ on subset $\mathcal{A}$ is defined as follows:

where $\odot$ is the Hadamard (element-wise) product of vectors. If the $i$-th feature is not included in $\mathcal{A}$, we keep its original value, otherwise, we change it to the closest endpoint of its normal range.

As there could be a large number of subsets $\mathcal{A}$ that can achieve the desired prediction after feature replacement, we aim to find the minimal subsets $\mathcal{A}^{*}$. Note that a minimal subset does not necessarily indicate minimal cardinality due to the difficulty of comparing the costs of different features. For example, increases in “academic degree”, “credit score”, and “salary” cannot be measured concisely in a mathematical cost. As such, we treat the subset (on “academic degree”, “credit score”) and the subset (on “salary”) as two different solutions, although the latter subset has the smaller cardinality, same to (Wang et al., 2021). Based on the above discussion, we formulate this problem to generate nondominated sets of CFEs as follows:

Diversity Analysis: In (Mothilal et al., 2020), authors propose a diversity metric over a set of CFEs of size $k$, named count-diversity, to measure the inner discrepancy, written by,

Robustness Analysis: Let $\mathbf{z}=\mathbf{c}-\mathbf{x}$ represent the recommended actions for a user. In our method, $\mathbf{z}$ consistently applies to slightly perturbed instance $\hat{\mathbf{x}}$, except in the following two situations: (1) $f(\hat{\mathbf{x}}+\mathbf{z})$ is no longer valid, which occurs when slight perturbations have a negative impact on the desired prediction. For example, normal features may be turned into abnormal ones. We need more effort than $\mathbf{z}$ to achieve the desired prediction. (2) Changing fewer abnormal features is sufficient to achieve the desired prediction, indicating that slight perturbations are beneficial. In this case, $\mathbf{z}$ is omitted as there exist more cost-efficient solutions. As both model continuity and perturbation strategies can influence the $\mathbf{z}$, we leave the determination of the maximal bound of perturbation, to which our method remains robust, for future work.

The brute-force method that evaluates all possible subsets is exponentially complex with respect to the number of abnormal features. Next, we propose a technique Counterfactual Explanations with Minimal Satisfiable Perturbation (CEMSP) to boost the search process. Our method starts with finding the binary vectors $\mathbf{m}$ that satisfy the desired prediction after feature replacement. This can be converted into the Boolean satisfiability problem that checks whether there exists a Boolean value assignment on $d_{0}$ variables (features in abnormal ranges) such that the conjunction of Boolean formulas evaluates to $True$. For better efficiency, we introduce the following proposition from domain knowledge.

This proposition aligns with common sense in practical applications. It is important to note that the undesired prediction arises from specific abnormal features according to our assumption. Intuitively, moving an abnormal feature into the normal range should never decrease the desired probability. Additionally, we assume that the predictive model $f$ has learned the relationship between feature normal ranges and the predicted classes.

Based on the monotonicity of function $f(r(\mathbf{x},\cdot))$, we derive the following two theorems for any $\mathcal{A}\subseteq\mathcal{B}\subseteq\mathcal{N}_{0}$.

Theorem 3 illustrates that if we can achieve the desired prediction by replacing abnormal features in $\mathcal{A}$, there is no need to change more abnormal features. The theorem tends to produce sparser results at a lower cost. Theorem 4 demonstrates that if we cannot achieve the desired prediction by changing abnormal features in $\mathcal{B}$, there is no need to check the satisfiability of any subsets of $\mathcal{B}$.

To prune as many as subsets at a time with these two theorems, we need to find the minimal satisfiable subset (MSS) and the maximal unsatisfiable subset (MUS), shown as boxes filled with green/red background in Figure 2. Next, we introduce two algorithms to achieve this: Grow($\cdot$) in Algorithm 1 and Shrink($\cdot$) in Algorithm 2. The Grow($\cdot$) algorithm starts with an arbitrary unsatisfiable subset $\hat{\mathcal{A}}$ and iteratively attempts to change other abnormal features until a maximal unsatisfiable subset is found. The Shrink($\cdot$) algorithm starts with an arbitrary satisfiable subset and iteratively attempts to reserve some features until a minimal satisfiable subset is found. Note that Grow($\cdot$) and Shrink($\cdot$) algorithms serve two plugins in our method, which can be replaced by any advanced algorithm with the same purpose.

Further, we introduce how to solve it under Boolean satisfiability problem (Liffiton and Malik, 2013; Bendík et al., 2018). In particular, any subset can be converted to a satisfiable Boolean assignment under a set of propositional logic formulas in conjunctive normal form (CNF), i.e., $\mathcal{A}\iff\mathbf{m}:\text{CNF}=True$. For example, for a subset $\mathcal{A}=\{1,2\}$ in Figure 2, we can write the CNF in the following equation, and $[1,1,0,0]$ is the only solution.

By employing this approach, the explicit materialization of all subsets can be avoided, thereby mitigating the exponential space complexity. The crux lies in devising the appropriate propositional logic formulas.

Our complete algorithm is shown in Algorithm 3. Initially, in line $1$, we merely forbid the changes on normal features and any possible binary assignments on abnormal features can satisfy the CNF. $\text{getMask}(\text{CNF})$ uses an SAT solver to return a solution satisfying the CNF in line $3$. In our paper, we adopt the Z3 package¹¹1https://github.com/Z3Prover/z3. In line $7$, we convert the binary vector $\mathbf{m}$ to indices of subsets. Next, we check whether replacing features in this subset can achieve the desired prediction. If not satisfying the desired prediction, we call the Grow($\cdot$) algorithm to find the maximal unsatisfiable subset and then prune all subsets of it. The prune operation is achieved by the following propositional Boolean formulas which are conjugated into the existing CNF.

Similarly, if we find a subset satisfying the desired prediction, we call Shrink($\cdot$) function to return a minimal satisfiable subset, which can induce a CFE with minimal perturbations of features. Then, we prune all supersets of it with the following logic formula,

If no solution satisfies the CNF in the current iteration, we stop our algorithm in the line $4$ and $5$.

Correctness. In our algorithm, each subset is either evaluated to be an MSS/MUS or pruned by an MSS/MUS. Hence, our algorithm will return all minimal CFEs, providing the same solutions as the brute-force search.

Space Complexity. The space complexity depends on how many minimal CFEs are returned. The worst case is $O(\binom{d_{0}}{\lfloor d_{0}/2\rfloor})$, which corresponds to that feature replacements on arbitrary abnormal features of size $\lfloor d_{0}/2\rfloor$ are satisfiable.

Time Complexity. The runtime of our method primarily depends on two parts: (1) a solver (line 3) that takes a set of constraints and returns a mask $\mathbf{m}$. The solver is considerably faster than calling the pretrained model. (2) evaluating the prediction on a subset. Compared with brute-force search, which calls the deep model $2^{d_{0}}$ times to check all $2^{d_{0}}$ possible sets, our method reduces the number of calls on the model by pruning on certain subsets and supersets. Therefore, the empirical running time will decrease.

The major theme of recent research is to model various constraints into CFE generation. Here, we show how to write these constraints by propositional logic formulas that can be conjugated into the CNF in line $1$ of our Algorithm 3.

Immutable features. Considering some features are immutable (e.g., race, birthplace), CFEs should avoid perturbations on these features. To achieve this, we add the following Boolean logic formula for a set of immutable features $\mathcal{I}$,

Accordingly, these features should be ignored in Grow($\cdot$) algorithm when it searches for the maximal unsatisfiable set. Alternatively, we can directly treat immutable features as normal features to avoid any changes to them.

Conditional immutable features. These features must change in one direction, e.g., education degree. We can examine whether moving a feature value into its normal range follows the valid direction. If violating the valid direction, we treat this feature as an immutable feature, otherwise, we put no restriction on this feature.

Causality. In practice, changing one feature may cause a change in other features. Such causal relations among features are generally written by a set of triplets in the structural causal model (SCM) (Pearl, 2009), that is, $\mathcal{M}=\langle U,V,F\rangle$, where $U$ are exogenous features, $V$ are endogenous features, and $F:U\to V$ is a set of functions that describe how endogenous features are quantitatively affected by exogenous features. To adapt causality to our method, we only keep these triplets that normal exogenous features lead to normal endogenous features as our method merely considers discrete feature changes (from abnormal to normal). For example, feature $\mathbf{x}^{1}$ is an exogenous feature that affects two endogenous features $\mathbf{x}^{2}$ and $\mathbf{x}^{3}$, and $\mathbf{x}^{2}$ and $\mathbf{x}^{3}$ become normal in the consequence of its normal ancestor feature $\mathbf{x}^{1}$. In this example, we can add two material conditions in the following that restrict the feature change of CFEs to follow the causal relations.

At the same time, Grow($\cdot$) and Shrink($\cdot$) algorithms should be updated to satisfy such causal relations when they attempt to add/remove a feature. This can be easily implemented by storing these causal relations by an inverted index where an entry is an exogenous feature and the inverted list contains all its endogenous features.

Correlation. Correlation can be regarded as bidirectional causal relations. For example, if features $\mathbf{x}^{1}$ and $\mathbf{x}^{2}$ are correlated and in the normal range simultaneously, we can write the correlation between $\mathbf{x}^{1}$ and $\mathbf{x}^{2}$ as,

The great advantage of our framework is that it allows us to insert these constraints gradually and flexibly, as the complete relation graphs (e.g., full causal graph) are often difficult to derive in the beginning.

We first report the quantitative comparison of sparsity, APS, Diversity, and C-Diversity in Figure 5, where the $x$-axis denotes a counterfactual explanation generation method, the $y$-axis represents the average score of each metric of all evaluated instances. We can see that our method achieves the competitive sparsity as GS. GS achieves sparsity through post-preprocessing techniques, whereas our method CEMSP focuses on making minimal modifications to subsets of abnormal features. In contrast, other methods do not explicitly optimize for sparsity and consequently fall behind in this aspect. For the APS, CEMSP is slightly better. This result is grounded in two key considerations: firstly, we substitute an abnormal feature with its closest endpoints and secondly, we aim to change the minimal number of abnormal features. It is worth noting that although PlainCF minimizes the $L_{1}/MAD$ distance in its objective, this does not equate to minimizing the APS since APS is a density-aware metric among the population. Our method achieves at least $\frac{2}{d}$ C-Diversity as expected. In contrast to sparsity, C-Diversity sums up the fraction of features that are different between any two CFEs. Therefore, the method with a higher sparsity often associates with a lower C-Diversity. As a result, our CEMSP appears to have less competitive C-Diversity than methods that make simultaneous changes on many numerous features. However, our CEMSP has a competitive Diversity score defined in (Mothilal et al., 2020).

Figure 5 and 5 report inconsistency scores under model retraining and input perturbations. Our CEMSP exhibits superior performance compared to other baseline methods. Specifically, GS, plainCF, and DiCE yield the poorest results when consistency restrictions are not enforced. CFProto, which incorporates a prototype term, achieves a better inconsistency score than PlainCF by directing all CFEs towards the prototype. As discussed earlier, our findings demonstrate that SNS does not perform well in generating CFEs with consistent feature values, despite having CFEs that yield consistent model predictions. In summary, our CEMSP outperforms the baseline methods across the aforementioned metrics, establishing its overall superiority.

Next, we use the use-case evaluation in Figure 6 to present the compatibility of our method. The input instance is a patient in the HCV dataset who has the undesired prediction. Without any constraint, our method generates 4 CFEs, as illustrated in the table located at the bottom left. By introducing additional constraints, our model can effortlessly generate new CFEs that meet the desired criteria. For example, if we want to keep the original value of $BIL$, we can easily incorporate the CNF ($\lnot\mathbf{m}^{5}$), leading to CFEs that solely modify the remaining features. Furthermore, domain knowledge reveals that $ALT$ and $AST$ are correlated (Shen et al., 2015). Consequently, we can incorporate correlation constraints that limit simultaneous changes to both features. This can be achieved by including the CNF $(\lnot\mathbf{m}^{3}\land\mathbf{m}^{4})\vee(\mathbf{m}^{3}\land\lnot\mathbf{m}^{4})$ in our method, effectively enforcing the desired correlation constraint. Although this use-case evaluation may not yet be fully applicable to real-life scenarios, it offers valuable insights and demonstrates the potential to accommodate more practical considerations.