# Flexible and Robust Counterfactual Explanations with Minimal Satisfiable
Perturbations

Yongjie Wang  0000-0003-4718-0742<EMAIL_ADDRESS>Nanyang
Technological UniversitySingapore ,  Hangwei Qian  0000-0003-4831-0748
<EMAIL_ADDRESS>A*STARSingapore ,  Yongjie Liu
0009-0007-7935-8601<EMAIL_ADDRESS>Shandong UniversityChina ,  Wei
Guo  0000-0002-8124-5186<EMAIL_ADDRESS>Shandong UniversityChina and
Chunyan Miao  0000-0002-0300-3448<EMAIL_ADDRESS>Nanyang Technological
UniversitySingapore

(2023)

###### Abstract.

Counterfactual explanations (CFEs) exemplify how to minimally modify a feature
vector to achieve a different prediction for an instance. CFEs can enhance
informational fairness and trustworthiness, and provide suggestions for users
who receive adverse predictions. However, recent research has shown that
multiple CFEs can be offered for the same instance or instances with slight
differences. Multiple CFEs provide flexible choices and cover diverse
desiderata for user selection. However, individual fairness and model
reliability will be damaged if unstable CFEs with different costs are
returned. Existing methods fail to exploit flexibility and address the
concerns of non-robustness simultaneously. To address these issues, we propose
a conceptually simple yet effective solution named Counterfactual Explanations
with Minimal Satisfiable Perturbations (CEMSP). Specifically, CEMSP constrains
changing values of abnormal features with the help of their semantically
meaningful normal ranges. For efficiency, we model the problem as a Boolean
satisfiability problem to modify as few features as possible. Additionally,
CEMSP is a general framework and can easily accommodate more practical
requirements, e.g., casualty and actionability. Compared to existing methods,
we conduct comprehensive experiments on both synthetic and real-world datasets
to demonstrate that our method provides more robust explanations while
preserving flexibility.

Counterfactual explanations, multiplicity, normal ranges, flexibility,
robustness

††journalyear: 2023††copyright: rightsretained††conference: Proceedings of the
32nd ACM International Conference on Information and Knowledge Management;
October 21–25, 2023; Birmingham, United Kingdom††booktitle: Proceedings of the
32nd ACM International Conference on Information and Knowledge Management
(CIKM ’23), October 21–25, 2023, Birmingham, United Kingdom††doi:
10.1145/3583780.3614885††isbn: 979-8-4007-0124-5/23/10††ccs: Computing
methodologies Philosophical/theoretical foundations of artificial
intelligence††ccs: Computing methodologies Symbolic and algebraic
manipulation††ccs: Theory of computation Logic and verification††ccs:
Computing methodologies Knowledge representation and reasoning

##  1\. introduction

Understanding the internal mechanisms behind model predictions is difficult
due to the large volume of parameters in machine learning models. This problem
is particularly significant in high-stakes domains such as healthcare and
finance, where incorrect predictions are disastrous (Daws, 2021).
Counterfactual explanation (CFE) (Wachter et al., 2017) aims to identify
minimal changes required to modify the input to achieve a desired prediction
and provides insights into why a model produces a certain prediction instead
of the desired one. CFEs can help understand the underlying logic of certain
predictions (Schoeffer et al., 2022), detect the inherent model bias for
fairness (Kasirzadeh and Smart, 2021), and provide suggestions to users who
receive adverse predictions (Karimi et al., 2021; Ustun et al., 2019).
Therefore, CFEs can be adopted in broad applications of healthcare, finance,
education, justice, and other domains.

Despite the valuable insights provided by CFEs, recent studies (Wachter et
al., 2017; Schoeffer et al., 2022; Sokol and Flach, 2019; Virgolin and
Fracaros, 2023; Slack et al., 2021; Laugel et al., 2019) have shown that
multiple CFEs can exist with equivalent evaluation metrics (e.g., validity,
proximity, sparsity, plausibility), yet significantly differ on feature values
for an input or seemingly indifferent inputs. For example, Wachter et al.
(Wachter et al., 2017) pointed out “multiple counterfactuals are possible, …”
and “multiple outcomes based on changes to multiple variables may be
possible”. Laugel et al. (Laugel et al., 2019) demonstrated that instances
that are close to each other can have different CFEs. Moreover, Virgolin and
Fracaros (Virgolin and Fracaros, 2023) have argued that CFEs lack robustness
to adverse perturbations if not deliberately designed. While multiple CFEs can
lead to the same desired prediction, each CFE tells a different story to reach
the target.

It is important to note that counterfactual multiplicity has both advantages
and disadvantages. On one hand, multiple CFEs can be beneficial because they
afford more flexibility and freedom to select user-friendly CFEs when a single
CFE may be overly restricted for users. Specifically, diverse CFEs (Russell,
2019; Mothilal et al., 2020) are offered to potentially cover broad user
preferences; interactive human-computer interfaces (Cheng et al., 2021; Wang
et al., 2021) are designed on multiple CFEs to obtain more satisfiable ones.
On the other hand, users who have the same feature values or seemingly
inconsequential differences may receive inconsistent CFEs (e.g., two different
diverse sets) as the CFE method itself does not store historical CFEs and
guarantee the optimal solutions either. Such inconsistency inevitably raises
fairness issues (Virgolin and Fracaros, 2023; Artelt et al., 2021) and
undermines users’ trust (Schoeffer et al., 2022) in CFEs. For example, two
financially similar individuals are rejected when they apply for a loan. Yet,
CFEs for two people are quite different-one needs to update the salary
slightly while the other is required to get a higher education degree and a
better job. Another negative example is when users make some efforts towards
previous CFEs but receive a significantly different CFE, rendering their
previous efforts futile. Therefore, it is crucial to take advantage of the
flexibility of multiple CFEs and maintain consistency for users having the
same feature values or slight differences.

Recent research (Upadhyay et al., 2021; Black et al., 2022; Pawelczyk et al.,
2020b; von Kügelgen et al., 2022; Dutta et al., 2022) on robustness mainly
studies CFEs with consistent predictions under slight model updates (by
restricting CFEs to preserve causal constraints, or follow data distribution,
etc.), rather than generates CFEs with consistent feature values. Therefore,
these studies fail to address fairness concerns and ignore the freedom of user
selection. Generally, models to explain are highly complex and non-convex,
e.g., DNN models. Even constrained by consistent prediction, heuristic search
strategies (Wachter et al., 2017; Laugel et al., 2017; Sharma et al., 2020;
Tolomei et al., 2017) can still converge to different non-optimal solutions
due to the huge search space. Meanwhile, these works do not explicitly exploit
flexibility to meet user preferences. As the number of possible CFEs can be
huge, existing methods on flexibility or robustness can be explained to be
different selection strategies from a solution pool (without optimizing
diversity and robustness in advance), i.e., selecting CFEs that are diverse
(Mothilal et al., 2020), follow data distribution (Pawelczyk et al., 2020b),
or withhold causal constraints (von Kügelgen et al., 2022). Motivated by this,
we target to design a novel method that obtains a diverse and robust set of
CFEs simultaneously.

To overcome the above limitations, we propose to incorporate task priors
(normal ranges, a.k.a. reference intervals) to stabilize valid search regions,
while ensuring that counterfactual explanations (CFEs) are diverse to meet
various user requirements. It should be noted that robustness measures the
differences between two sets of CFEs in different trials while diversity
measures the inherent discrepancy within a set of CFEs. Normal ranges in our
approach commonly exist in broad domains and are easy to obtain from prior
knowledge. For example, the normal range of heart rate per minute is between
60 and 100; the IELTS score should be greater than 6.5 and the minimal GPA is
3.5 for Ph.D. admissions. We assume that the undesired prediction results from
certain features outside of normal ranges and thus, we attempt to move
abnormal features into normal ranges to generate CFEs with the desired
prediction. Specifically, we replace an abnormal feature with the closest
endpoint of its normal range. As the endpoints are stationary, CFEs after
feature replacement tend to have the same feature value in different trials
for the same/similar input. In practice, it may be unnecessary to move all
abnormal features into normal ranges for the desired prediction. Therefore, we
aim to select minimal subsets of abnormal features to replace where each
subset corresponds to a CFE. CFEs determined by all minimal subsets are
diverse as an arbitrary minimal subset is not contained by another subset.

As mentioned earlier, the problem of finding CFEs boils down to selecting
minimal subsets of abnormal features to replace, which can be formulated as
either the Maximally Satisfiable Subsets (MSS) or Minimal Unsatisfiable
Subsets (MUS) problem (Liffiton and Malik, 2013; Liffiton and Sakallah, 2005).
However, finding all minimal sets for satisfiable CFEs is an NP-Complete
problem as an exponential number of subsets should be checked. To enhance
efficiency, we covert the enumeration of minimal subsets to the Boolean
satisfiability problem (SAT) that finds satisfiable Boolean assignments over a
series of Boolean logic formulas, which can be solved with efficient modern
solvers. As for commonly mentioned constraints (e.g., actionability,
correlation, and causality), we can conveniently write them in Boolean logic
formulas, which can be conjugated into current clauses in conjunctive normal
form (CNF). Therefore, our framework is flexible to provide feasible
counterfactual recommendations.

The main contributions of this paper are summarized as follows.

  * •

We reformulate the counterfactual explanation problem to satisfy both
flexibility and robustness by replacing a minimal subset of abnormal features
with the closest endpoints of normal ranges.

  * •

We convert this problem by checking the satisfiability of Boolean logic
formulas for a Boolean assignment, which can be solved by modern SAT solvers
efficiently. In addition, common constraints can be easily incorporated into
current Boolean logic formulas and solved together.

  * •

We conduct intensive experiments on both synthetic and real-world datasets to
demonstrate that our approach produces more consistent and diverse CFEs than
state-of-the-art methods.

##  2\. Related Work

Counterfactual explanations (Wachter et al., 2017) refer to perturbed
instances with the minimum cost that result in a different prediction from a
pre-trained model given an input instance. These explanations provide ways to
comprehend the model’s prediction logic and offer advice to users receiving
adverse predictions. Most existing algorithms focus on modeling practical
requirements and user preferences with proper constraints. Typical constraints
include actionability (Ustun et al., 2019), which freezes immutable features
such as race, gender, etc; plausibility (Joshi et al., 2019; Van Looveren and
Klaise, 2019), which requires CFEs to follow the data distribution; diversity
(Mothilal et al., 2020; Russell, 2019), which generates a diverse set of
explanations at a time; sparsity (Dhurandhar et al., 2018), which favors fewer
features changed; causality (Joshi et al., 2019; Karimi et al., 2020), which
restricts CFEs to meet specific causal relations. However, recent studies
(Wachter et al., 2017; Laugel et al., 2019; Russell, 2019; Schoeffer et al.,
2022; Sokol and Flach, 2019; Virgolin and Fracaros, 2023) have revealed that
there often exist multiple CFEs with equivalent performance but different
feature values for an input or seemingly indifferent inputs. Next, we review
research that takes advantage of counterfactual multiplicity and addresses
concerns regarding non-robustness.

Multiple CFEs provide users with more flexibility to prioritize their
preferences without compromising the validity and proximity of CFEs. When a
single CFE is inadequate to meet users’ requirements, employing a diverse set
of CFEs is an effective and straightforward strategy to overcome this
limitation. For example, Wachter et al. (Wachter et al., 2017) generate a
diverse set by running multiple times with different initializations; Russell
(Russell, 2019) prohibits the transition to previous CFEs in each run; while
Mothilal et al. (Mothilal et al., 2020) add a DPP (Determinantal Point
Processes) term to ensure that the CFEs are far apart from each other. In
addition, multiple CFEs also enable researchers to develop Human-Computer
Interaction (HCI) tools for interactively satisfying user requirements (Cheng
et al., 2021; Wang et al., 2021). However, such a diverse set can be
inconsistent for two inputs with no or slight differences. In our paper, we
aim to generate a consistent and diverse set of CFEs for an input or seemingly
different inputs, to enhance the robustness and reliability of CFEs.

The non-robustness issue of CFEs has garnered significant attention recently.
As introduced in (Slack et al., 2021; Laugel et al., 2019), even a slight
perturbation to the input can result in drastically different CFEs. To verify
this phenomenon for neural network models, Slack et al. (Slack et al., 2021)
train an adversarial model that is sensitive to trivial input changes. Some
relevant works (Upadhyay et al., 2021; Black et al., 2022; Pawelczyk et al.,
2020b; Dutta et al., 2022) propose to generate CFEs that yield consistent
predictions when the model is retrained. For example, (Pawelczyk et al.,
2020b) proves that adhering to the data manifold ensures stable predictions
for CFEs; (Upadhyay et al., 2021) incorporates adversarial training to produce
robust models for generating explanations; Black et al. (2022) states that
closeness to the data manifold is insufficient to indicate counterfactual
stability, and they propose Stable Neighbor Search (SNS) to find an
explanation with the lower model Lipschitz continuity and higher confidence.
However, constraining consistent predictions of CFEs does not necessarily
ensure CFEs with the same or similar feature values, and still fails to
address the unfairness issue. Moreover, these robustness methods lack the
flexibility to meet user requirements while our work considers both
flexibility and robustness simultaneously.

##  3\. Preliminary

###  3.1. Counterfactual Explanations

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:

(1) $\displaystyle\mathop{\arg\min}_{\mathbf{c}}\ $ $\displaystyle\
\text{cost}(\mathbf{x},\mathbf{c})$ $\displaystyle s.t.\ $ $\displaystyle\
f(\mathbf{c})=1$

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:

(2)
$\displaystyle\mathcal{L}(\mathbf{c},\lambda)=\text{cost}(\mathbf{x},\mathbf{c})+\lambda\ell(f(\mathbf{c}),1)$

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).

###  3.2. Robustness of Counterfactual Explanations

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.

###### Definition 0 (Robustness of Counterfactual Explanations).

Given a function $d(\cdot,\cdot)$ computing the similarity between two sets of
CFEs, we quantify the robustness of the explanations $g(f,\mathbf{x})$ by
assessing the expected similarity between the current set of CFEs, and a new
set of CFEs after potential input perturbations or model changes.

(3) $\displaystyle\mathop{\mathbb{E}}_{\mathbf{x}^{\prime}\sim
p(\mathbf{x}),\atop\hat{f}\in\mathcal{F}}[d(g(f,\mathbf{x}),g(\hat{f},\mathbf{x}^{\prime}))]$

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}$.

Figure 1. The toy example to illustrate the non-robustness issue caused by
input perturbations for an input $(0.3,0)$. The figures in the top row depict
the decision surface, $L_{2}$ cost, and total loss in the Lagrangian form
($\lambda=5$ for simplicity). The bottom row shows the contour lines
corresponding to each figure above. It is evident that even a slight
perturbation on this input can result in different local minima.

###  3.3. Causes of Non-robustness

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).

##  4\. Proposed Method

In this section, we propose a method offering a diverse set of CFEs that
remain stable when the input or model being explained is changed subtly by
incorporating domain knowledge into the search.

###  4.1. Problem Formulation

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:

(4)
$\displaystyle\eta(\mathbf{x},\mathcal{A})^{i}=\begin{cases}\mathbf{a}^{i},&i\in\mathcal{A}\text{
and }\mathbf{x}^{i}<\mathbf{a}^{i}\\\ \mathbf{b}^{i},&i\in\mathcal{A}\text{
and }\mathbf{x}^{i}>\mathbf{b}^{i}\\\
\mathbf{x}^{i},&\text{otherwise}\end{cases}$

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:

(5)
$\displaystyle\eta(\mathbf{x},\mathcal{A})^{i}=\begin{cases}\hat{\mathbf{x}}^{i},&i\in\mathcal{A}\\\
\mathbf{x}^{i},&\text{otherwise}\end{cases}$

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:

(6) $\displaystyle
r(\mathbf{x},\mathcal{A})=\mathbf{x}\odot(1-\mathbf{m})+\eta(\mathbf{x},\mathcal{A})\odot\mathbf{m}$

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:

###### Definition 0 (Minimal Satisfiable Counterfactual Explanations).

Given a pre-trained model $f$, input $\mathbf{x}$ and feature normal ranges,
the goal is to find all minimal satisfiable subsets, where for each subset
$\mathcal{A}^{*}\subseteq\mathcal{N}_{0}$, $f(r(\mathbf{x},\mathcal{A}^{*}))$
can achieve the desired prediction and $r(\mathbf{x},\mathcal{A}^{*})$
represents a counterfactual explanation $\mathbf{c}$.

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,

(7) $\textit{count-
diversity}=\frac{2}{k(k-1)d}\sum_{i=1}^{k-1}\sum_{j=i+1}^{k}\sum_{l=1}^{d}\mathds{1}_{[\mathbf{c}_{i}^{l}\neq\mathbf{c}_{j}^{l}]}$

###### Theorem 2.

The lower bound of count-diversity defined in Eqn. (7) over solutions of our
problem definition is $\frac{2}{d}$.

###### Proof.

Let’s denote the inner loop term
$\sum_{l=1}^{d}\mathds{1}_{[\mathbf{c}_{i}^{l}\neq\mathbf{c}_{j}^{l}]}$ as
$D(\mathbf{c}_{i},\mathbf{c}_{j})$ for brevity. We aim to prove that for any
two arbitrary CFEs $\mathbf{c}_{i}$ and $\mathbf{c}_{j}$, the pairwise
distance $D(\mathbf{c}_{i},\mathbf{c}_{j})$ is always greater than or equal to
$2$. To demonstrate this, we consider two contradictive cases by assuming
$0\leq D(\mathbf{c}_{i},\mathbf{c}_{j})<2$.

######  Case 1 ($D(\mathbf{c}_{i},\mathbf{c}_{j})=0$):

This case implies $\mathbf{c}_{i}=\mathbf{c}_{j}$, which contradicts the fact
that two distinct CFEs are from the solution set.

######  Case 2 ($D(\mathbf{c}_{i},\mathbf{c}_{j})=1$):

In this case, there is only one feature difference. Let $\mathcal{A}_{i}$ and
$\mathcal{A}_{j}$ denote the indices of abnormal features of two solutions.
Then, $\mathcal{A}_{i}\subseteq\mathcal{A}_{j}$ or
$\mathcal{A}_{j}\subseteq\mathcal{A}_{i}$ must exist, which also contradicts
the minimal set to return in our problem definition. One CFE should be
excluded because it costs more than the other.

The above contradictive proof shows that $D(\mathbf{c}_{i},\mathbf{c}_{j})\geq
2$ holds. Summing up all $\frac{k(k-1)}{2}$ pair-wise distance
$D(\mathbf{c}_{i},\mathbf{c}_{j})$, we can obtain the lower bound
$\frac{2}{d}$. ∎

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.

###  4.2. Problem Solving

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.

Figure 2. The figure shows all subsets of a toy example with $4$ abnormal
features. The bitvectors denote the binary vector $\mathbf{m}$. Boxes with
red/green borders represent the unsatisfiable/satisfiable subsets
respectively. The minimal subsets for CFEs are filled with green background
and the maximal unsatisfiable subsets are filled with red background.

###### Proposition 1 (Monotoncity of $f(r(\mathbf{x},\cdot))$).

The function $f(r(\mathbf{x},\cdot))$ is monotone, that is,
$f(r(\mathbf{x},\mathcal{A}))\leq f(r(\mathbf{x},\mathcal{B}))$ holds for all
$\mathcal{A}\subseteq\mathcal{B}\subseteq\mathcal{N}_{0}$.

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.

If $f(r(\mathbf{x},\mathcal{A}))$ can achieve the desired target,
$f(r(\mathbf{x},\mathcal{B}))$ can also achieve the desired target for any
superset $\mathcal{B}$ of $\mathcal{A}$.

###### Theorem 4.

If $f(r(\mathbf{x},\mathcal{B}))$ cannot achieve the desired target,
$f(r(\mathbf{x},\mathcal{A}))$ cannot satisfy the desired target either for
any subset $\mathcal{A}$ of $\mathcal{B}$.

###### Proof.

We first prove theorem 3. If $f(r(\mathbf{x},\mathcal{A}))\geq\delta$, we can
induce that $f(r(\mathbf{x},\mathcal{B}))\geq\delta$ as
$f(r(\mathbf{x},\mathcal{B}))\geq f(r(\mathbf{x},\mathcal{A}))$ holds for
$\mathcal{A}\subseteq\mathcal{B}$, where $\delta$ is the confidence threshold
of desired prediction. Similarly, theorem 4 can be proved. ∎

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.

Algorithm 1 Grow($\mathcal{A}$)

1:An unsatisfiable subset $\mathcal{A}$.

2:A maximal unsatisfiable subset $\hat{\mathcal{A}}$.

3:$\hat{\mathcal{A}}=\mathcal{A}$

4:for $i\in\mathcal{N}_{0}\backslash\mathcal{A}$ do

5: if $f(r(\mathbf{x},\hat{\mathcal{A}}\bigcup\\{i\\}))==0$ then

6: $\hat{\mathcal{A}}=\hat{\mathcal{A}}\bigcup\\{i\\}$

7: end if

8:end for

9:return $\hat{\mathcal{A}}$

Algorithm 2 Shrink($\mathcal{A}$)

1:A subset $\mathcal{A}$ satisfying the desired prediction.

2:A minimal subset $\mathcal{A}^{*}$ for a CFE.

3:$\mathcal{A}^{*}=\mathcal{A}$

4:for $i\in\mathcal{A}$ do

5: if $f(r(\mathbf{x},\mathcal{A}^{*}\backslash\\{i\\}))==1$ then

6: $\mathcal{A}^{*}=\mathcal{A}^{*}\backslash\\{i\\}$

7: end if

8:end for

9:return $\mathcal{A}^{*}$

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.

(8)
$\text{CNF}=\mathbf{m}^{1}\land\mathbf{m}^{2}\land\lnot\mathbf{m}^{3}\land\lnot\mathbf{m}^{4}$

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
package111https://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.

(9)
$\text{PruneSubSet}(\hat{\mathcal{A}})=\vee_{i:i\in\mathcal{N}_{0}\backslash\hat{\mathcal{A}}}\
\mathbf{m}^{i}$

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,

(10) $\text{PruneSuperSet}(\mathcal{A}^{*})=\vee_{i:i\in\mathcal{A}^{*}}\
\lnot\mathbf{m}^{i}$

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.

Algorithm 3 CEMSP($\mathbf{x},f$)

1:An input $\mathbf{x}$, a pretrained model $f$.

2:All minimal subsets $\mathcal{A}^{*}$ for CFEs.

3:$\text{CNF}=\land_{i:i\in\mathcal{N}_{1}}\ \lnot\mathbf{m}^{i}$

4:while True do

5: $\mathbf{m}\leftarrow\text{getMask}(\text{CNF})$

6: if not $\mathbf{m}$ then $\triangleright$ No assignment $\mathbf{m}$
returned.

7: Break

8: end if

9: $\mathcal{A}\leftarrow\\{i\in\mathcal{N}_{0}:\mathbf{m}^{i}=1\\}$

10: if $f(r(\mathbf{x},\mathcal{A}))==0$ then

11: $\hat{\mathcal{A}}\leftarrow\text{Grow}(\mathcal{A})$

12: $\text{CNF}\to\text{CNF}\wedge\text{PruneSubset}(\hat{\mathcal{A}})$
$\triangleright$ Prune any subset of $\hat{\mathcal{A}}$

13: else

14: $\mathcal{A}^{*}\leftarrow\text{Shrink}(\mathcal{A})$

15: $\text{{yield} }\mathcal{A}^{*}$

16: $\text{CNF}\to\text{CNF}\wedge\text{PruneSuperset}(\mathcal{A}^{*})$
$\triangleright$ Prune any superset of $\mathcal{A}^{*}$

17: end if

18:end while

###  4.3. Compatibility with Other Constraints

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}$,

(11) $\displaystyle\text{Actionability}(\mathcal{I})=\land_{i\in\mathcal{I}}\
\lnot\mathbf{m}^{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.

(12)
$\displaystyle\text{Causality}=(\lnot\mathbf{m}^{1}\vee\mathbf{m}^{2})\land(\lnot\mathbf{m}^{1}\vee\mathbf{m}^{3})$

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,

(13)
$\displaystyle\text{Correlation}=(\lnot\mathbf{m}^{1}\vee\mathbf{m}^{2})\land(\lnot\mathbf{m}^{2}\vee\mathbf{m}^{1})$

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.

##  5\. Experiments

In this section, we undertake a quantitative comparison between our proposed
method CEMSP and state-of-the-art approaches. Additionally, we demonstrate
empirical examples of counterfactual explanations that effectively integrate
practical constraints. The source code is available at the GitHub
repository222https://github.com/wangyongjie-ntu/CEMSP.

#### Datasets.

We conducted a comprehensive series of experiments involving a synthetic
dataset and two real-world UCI medical datasets. Notably, the medical datasets
encompass diagnostic features with well-defined and clinically significant
normal ranges.

  * •

Synthetic Dataset is a binary class dataset consisting of $20,000$ samples
with $4$ features. Each feature is sampled from the normal distribution
independently. Regarding label balance, the binary label $y$ is assigned a
value of $1$ when the following equation is satisfied; otherwise, $y$ is set
to 0:

$\displaystyle(\mathbf{x}^{1}>0.5)\lor(\mathbf{x}^{2}>0.4\land\mathbf{x}^{3}>0)\lor(\mathbf{x}^{2}>0.4\land\mathbf{x}^{3}>0.5)$

We set the lower values of normal ranges of four features as
$[0.55,0.45,0.05,0.55]$ for a higher confidence prediction.

  * •

UCI HCV Dataset (Diaconis and Efron, 1983). This dataset contains $615$
instances. Following (BHATT, 2021), we convert $5$ categories of diagnosis
into binary classes. After label conversion, the dataset consists of $75$
individuals diagnosed with HCV and $540$ individuals labeled as healthy. Next,
we remove the “Age” and “Sex” and keep the other $10$ medical features with
normal ranges. We adopt the tight normal ranges from laboratory tests in
(Hoffmann et al., 2016) as certain normal ranges depend on “Sex” and we remove
the “Sex” attribute in preprocessing.

  * •

UCI Thyroid Dataset (Quinlan et al., 1987). The raw dataset contains $3,772$
instances where each instance is described by $15$ features and labeled as
either hypothyroid or normal class. We retain the most discriminative features
“FTI”, “TSH”, “T3”, “TT4” that have meaningful normal ranges and remove other
features. Subsequently, we drop certain rows with missing values. The final
dataset consists of $223$ patients and $2530$ healthy users. Normal ranges of
“TSH”, “T3”, “TT4” are from laboratory tests in (Joshi et al., 2011). As we do
not find the normal range of “FTI” that matches the “FTI” values in this
dataset, we simply choose the 1-sigma interval of “FTI” of the normal group.

#### Evaluation Metrics.

To comprehensively compare over CFEs across various approaches, we employ the
following evaluation metrics.

  * •

Inconsistency. We propose to adopt a modified Hausdorff distance (Dubuisson
and Jain, 1994; Jesorsky et al., 2001) to measure the inconsistency between
two sets of CFEs $\mathcal{C}$ and $\mathcal{C}^{\prime}$,

(14) $\displaystyle
H(\mathcal{C},\mathcal{C}^{\prime})=\max(h_{mod}(\mathcal{C},\mathcal{C}^{\prime}),h_{mod}(\mathcal{C}^{\prime},\mathcal{C}))$

where
$h_{mod}(\mathcal{C},\mathcal{C}^{\prime})=\frac{1}{|\mathcal{C}|}\sum_{\mathbf{c}\in\mathcal{C}}\min_{\mathbf{c}^{\prime}\in\mathcal{C}^{\prime}}||\mathbf{c}-\mathbf{c}^{\prime}||_{2}$
and the lower $H$ is better.

  * •

Average Percentile Shift (APS) (Pawelczyk et al., 2020a) measures the relative
cost of perturbations of CFEs,

(15)
$\displaystyle\text{APS}(\mathbf{x},\mathcal{C})=\frac{1}{d*|\mathcal{C}|}\sum_{\mathbf{c}\in\mathcal{C}}\sum_{i=1}^{d}|Q^{i}(\mathbf{c}^{i})-Q^{i}(\mathbf{x}^{i})|$

where $Q^{i}(\cdot)$ denotes the percentile of the $i$-th feature value
relative to all values of the feature in the whole data set. A lower score is
favored.

  * •

Sparsity. It measures the percentage of features that remain unchanged and we
prefer higher sparsity,

(16)
$\displaystyle\text{Sparsity}(\mathbf{x},\mathcal{C})=\frac{1}{d*|\mathcal{C}|}\sum_{\mathbf{c}\in\mathcal{C}}\sum_{i=1}^{d}\mathds{1}_{\mathbf{c}^{i}=\mathbf{x}^{i}}.$

  * •

Diversity. We consider two diversity metrics, named Diversity, which is
introduced in (Mothilal et al., 2020), and count-diversity (named C-Diversity
for abbreviation), which is defined in Eqn (7) in Section 4.1, to measure the
discrepancy within returned solutions.

(17)
$\displaystyle\text{Diversity}=\frac{2}{k(k-1)}\sum_{i=1}^{k-1}\sum_{j=i+1}^{k}\text{dist}(\mathbf{c}_{i},\mathbf{c}_{j})$

where $\text{dist}(\cdot,\cdot)$ represents the $L_{1}/MAD$ and $k$ is the
number of CFEs.

#### Baselines.

We compare our method with the following baseline methods.

  * •

GrowingSphere (GS) (Lash et al., 2017). This algorithm searches for CFEs from
random samples in the sphere neighborhood of the input. The radius of the
sphere grows until a CFE is found. It adopts the postprocessing on returned
CFEs to make sparser solutions.

  * •

PlainCF (Wachter et al., 2017). It minimizes the objective in Eqn. (2) with
gradient descent. We run this algorithm from a random initial point and stop
when an iteration threshold is reached or the loss difference is below a
specified threshold.

  * •

CFProto (Van Looveren and Klaise, 2019). It adds a prototype term to restrict
that CFEs should resemble the prototype of the desired class. In our
experiment, we set the prototype as the closest endpoints of normal ranges of
these abnormal features, that is, $r(\mathbf{x},\mathcal{N}_{0})$.

  * •

DiCE (Mothilal et al., 2020). Compared with PlainCF, it considers the
diversity constraint that is modeled by a $dpp(\cdot)$ term over a set of
CFEs.

  * •

SNS (Black et al., 2022). It finds a CFE with higher confidence and lower
Lipschitz constant in the neighborhood of a given CFE, to produce consistent
prediction under model update.

#### Experiment configurations.

We first randomly split the datasets into train/test sets at the ratio of
$7:3$ and normalize all features by a standard scaler on two UCI datasets (no
feature normalization on the synthetic dataset). Then, we train a 3-layer
Multilayer perceptron (MLP) model $f$ with Adam optimizer. The test accuracies
are 99%, 96%, and 98% on three datasets correspondingly. As we intend to
convert unhealthy patients to healthy ones, we produce CFEs for all correctly
classified patients in test sets for two UCI datasets. For saving time, we
only produce CFEs for 100 random true negative samples in the synthetic
dataset.

Our method produces a set of CFEs of varied size $k$, while GS, PlainCF,
CFProto, and SNS generate one at a time. We run GS, PlainCF, and CFProto $k$
times for fair comparisons to generate the same number of CFEs as ours. For
DiCE, we directly keep the size of the diverse set as ours. We evaluate all
CFE methods $g(f,\mathbf{x})$ under the following two kinds of slight updates
to measure the algorithm robustness.

  * •

Inputs are fixed, and we produce two sets of CFEs from two models that are
trained on the same dataset with different initializations.

  * •

Model is fixed, and we produce two sets of CFEs from an input $\mathbf{x}$ and
its perturbed instance $\mathbf{x}^{\prime}$, where
$\mathbf{x}^{\prime}=\mathbf{x}+\alpha$ and $\alpha$ is the random noise
sampling from a Gaussian distribution $\mathcal{N}(0,\sigma)$. In our
experiments, $\sigma\in\\{0.0001,0.001,0.01,0.1\\}$.

Figure 3. Evaluation of sparsity, APS, Diversity, and C-Diversity over three
datasets. The $\uparrow$/$\downarrow$ means the higher/lower score is better.
Diversity and C-Diversity of the Thyroid dataset are missing as our method
CEMSP only produces a single counterfactual explanation.

Figure 4. Evaluation of inconsistency score of model retraining.

Figure 5. Evaluation of inconsistency score of input perturbations. $x$-axis
represents the standard deviation $\sigma$ of added Gaussian noise.

###  5.1. Quantitative Evaluation

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.

###  5.2. Use-Case Evaluation

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.

Figure 6. Use-case evaluation of a patient in the HCV dataset. The top table
presents the original feature values of the patient, while the shaded features
represent the altered features of CFEs and their closest endpoints. The tables
below display CFEs before/after incorporating constraints. To incorporate
actionability and correlation, we introduced the expressions $\lnot m^{5}$ and
$(\lnot\mathbf{m}^{3}\land\mathbf{m}^{4})\vee(\mathbf{m}^{3}\land\lnot\mathbf{m}^{4})$,
respectively.

##  6\. Conclusion

Lacking robustness in counterfactual explanations can undermine both
individual fairness and model reliability. In this work, we present a novel
framework to generate robust and diverse counterfactual explanations (CFEs).
Our work leverages the feature normal ranges from domain knowledge and
generates CFEs that replace the minimal number of abnormal features to the
closest endpoints of their normal ranges. We convert this problem into the
Boolean satisfiability problem and solve it with modern SAT solvers.
Experiments on both synthetic and real-life datasets demonstrate that our
generated CFEs are more consistent than baselines while preserving flexibility
for user preferences.

##  7\. Limitations and Future Work

While our work offers the potential to address the non-robustness issue
through the utilization of domain knowledge, such as the normal ranges in
healthcare and finance, some limitations hinder its applicability in broader
contexts. Firstly, the scalability of the proposed method is underestimated.
SAT solvers exhibit exponential complexity in the worst-case scenario. When
dealing with a substantial number of features, the time required to find a
binary mask from an SAT solver may surpass that of a forward pass in the DNN
model. This concern can be addressed through empirical comparisons of the
execution time between the SAT solver and DNN model revoking. Secondly, our
approach is not directly applicable to scenarios where a portion of normal
ranges is unknown. It might be necessary to incorporate additional information
to determine the appropriate replacement values for these features. Thirdly,
our study is established on binary classification tasks. However, the direct
adaptation of our method to multi-class classification or regression tasks
remains challenging. Normal ranges are typically contingent upon the target
prediction. In the context of multi-class classification or regression, the
target predictions can become intricate, rendering the normal ranges
unattainable.

In future work, our ultimate goal is to investigate robust and flexible
counterfactual explanations in more general situations without any assumption
about normal ranges. In addition, we intend to develop a sustainable system
that offers users actionable recommendations and gathers valuable feedback to
nourish continual enhancements.

###### Acknowledgements.

This research is supported, in part, by Alibaba Group through Alibaba
Innovative Research (AIR) Program and Alibaba-NTU Singapore Joint Research
Institute (JRI), Nanyang Technological University, Singapore. This research is
also supported, in part, by the National Research Foundation, Prime Minister’s
Office, Singapore under its NRF Investigatorship Programme (NRFI Award No.
NRF-NRFI05-2019-0002). Any opinions, findings and conclusions or
recommendations expressed in this material are those of the authors and do not
reflect the views of National Research Foundation, Singapore.

## References

  * (1)
  * Artelt et al. (2021) André Artelt, Valerie Vaquet, Riza Velioglu, Fabian Hinder, Johannes Brinkrolf, Malte Schilling, and Barbara Hammer. 2021.  Evaluating robustness of counterfactual explanations. In _2021 IEEE Symposium Series on Computational Intelligence (SSCI)_. IEEE, 01–09.
  * Bendík et al. (2018) Jaroslav Bendík, Ivana Černá, and Nikola Beneš. 2018.  Recursive online enumeration of all minimal unsatisfiable subsets. In _International symposium on automated technology for verification and analysis_. Springer, 143–159.
  * BHATT (2021) KRISHNA BHATT. 2021\.  Hepatitis C Predictions.  https://www.kaggle.com/code/krishnabhatt4/hepatitis-c-predictions
  * Black et al. (2022) Emily Black, Zifan Wang, and Matt Fredrikson. 2022.  Consistent Counterfactuals for Deep Models. In _International Conference on Learning Representations_.  https://openreview.net/forum?id=St6eyiTEHnG
  * Blechschmidt et al. (2020) Jana Blechschmidt, Meike J Wittmann, and Chantal Blüml. 2020\.  Climate change and green sea turtle sex ratio—preventing possible extinction.  _Genes_ 11, 5 (2020), 588.
  * CDC (2022) USA CDC. 2022.  Assessing Your Weight.  https://www.cdc.gov/healthyweight/assessing/index.html
  * Cheng et al. (2021) F. Cheng, Y. Ming, and H. Qu. 2021.  DECE: Decision Explorer with Counterfactual Explanations for Machine Learning Models.  _IEEE Transactions on Visualization & Computer Graphics_ 27, 02 (feb 2021), 1438–1447.  https://doi.org/10.1109/TVCG.2020.3030342
  * Daws (2021) Ryan Daws. 2021\.  Medical chatbot using OpenAI’s GPT-3 told a fake patient to kill themselves.  _AI News_ (11 2021).  https://www.artificialintelligence-news.com/2020/10/28/medical-chatbot-openai-gpt3-patient-kill-themselves/
  * Dhurandhar et al. (2018) Amit Dhurandhar, Pin-Yu Chen, Ronny Luss, Chun-Chen Tu, Paishun Ting, Karthikeyan Shanmugam, and Payel Das. 2018.  Explanations based on the missing: Towards contrastive explanations with pertinent negatives. In _Advances in Neural Information Processing Systems_. 592–603.
  * Diaconis and Efron (1983) Persi Diaconis and Bradley Efron. 1983.  Computer-intensive methods in statistics.  _Scientific American_ 248, 5 (1983), 116–131.
  * Dubuisson and Jain (1994) M-P Dubuisson and Anil K Jain. 1994.  A modified Hausdorff distance for object matching. In _Proceedings of 12th international conference on pattern recognition_ , Vol. 1. IEEE, 566–568.
  * Dutta et al. (2022) Sanghamitra Dutta, Jason Long, Saumitra Mishra, Cecilia Tilli, and Daniele Magazzeni. 2022.  Robust counterfactual explanations for tree-based ensembles. In _International Conference on Machine Learning_. PMLR, 5742–5756.
  * Hoffmann et al. (2016) Georg Hoffmann, Ralf Lichtinghagen, and Werner Wosniok. 2016\.  Simple estimation of reference intervals from routine laboratory data.  _LaboratoriumsMedizin_ 39, s1 (2016).
  * Jesorsky et al. (2001) Oliver Jesorsky, Klaus J Kirchberg, and Robert W Frischholz. 2001\.  Robust face detection using the hausdorff distance. In _International conference on audio-and video-based biometric person authentication_. Springer, 90–95.
  * Joshi et al. (2019) Shalmali Joshi, Oluwasanmi Koyejo, Warut Vijitbenjaronk, Been Kim, and Joydeep Ghosh. 2019.  Towards realistic individual recourse and actionable explanations in black-box decision making systems.  _arXiv preprint arXiv:1907.09615_ (2019).
  * Joshi et al. (2011) Shashank R Joshi et al. 2011\.  Laboratory evaluation of thyroid function.  _J Assoc Physicians India_ 59, Suppl (2011), 14–20.
  * Karimi et al. (2021) Amir-Hossein Karimi, Bernhard Schölkopf, and Isabel Valera. 2021.  Algorithmic Recourse: from Counterfactual Explanations to Interventions. In _FAccT ’21: 2021 ACM Conference on Fairness, Accountability, and Transparency, Virtual Event / Toronto, Canada, March 3-10, 2021_ , Madeleine Clare Elish, William Isaac, and Richard S. Zemel (Eds.). ACM, 353–362.  https://doi.org/10.1145/3442188.3445899
  * Karimi et al. (2020) Amir-Hossein Karimi, Julius von Kügelgen, Bernhard Schölkopf, and Isabel Valera. 2020\.  Algorithmic recourse under imperfect causal knowledge: a probabilistic approach.  _Advances in Neural Information Processing Systems_ 33 (2020).
  * Kasirzadeh and Smart (2021) Atoosa Kasirzadeh and Andrew Smart. 2021.  The Use and Misuse of Counterfactuals in Ethical Machine Learning. In _FAccT ’21: 2021 ACM Conference on Fairness, Accountability, and Transparency, Virtual Event / Toronto, Canada, March 3-10, 2021_ , Madeleine Clare Elish, William Isaac, and Richard S. Zemel (Eds.). ACM, 228–236.  https://doi.org/10.1145/3442188.3445886
  * Lash et al. (2017) Michael T Lash, Qihang Lin, Nick Street, Jennifer G Robinson, and Jeffrey Ohlmann. 2017.  Generalized inverse classification. In _Proceedings of the 2017 SIAM International Conference on Data Mining_. SIAM, 162–170.
  * Laugel et al. (2019) Thibault Laugel, Marie-Jeanne Lesot, Christophe Marsala, and Marcin Detyniecki. 2019.  Issues with post-hoc counterfactual explanations: a discussion.  _arXiv preprint arXiv:1906.04774_ (2019).
  * Laugel et al. (2017) Thibault Laugel, Marie-Jeanne Lesot, Christophe Marsala, Xavier Renard, and Marcin Detyniecki. 2017\.  Inverse classification for comparison-based interpretability in machine learning.  _arXiv preprint arXiv:1712.08443_ (2017).
  * Liffiton and Malik (2013) Mark H Liffiton and Ammar Malik. 2013.  Enumerating infeasibility: Finding multiple MUSes quickly. In _International Conference on Integration of Constraint Programming, Artificial Intelligence, and Operations Research_. Springer, 160–175.
  * Liffiton and Sakallah (2005) Mark H Liffiton and Karem A Sakallah. 2005.  On finding all minimally unsatisfiable subformulas. In _International conference on theory and applications of satisfiability testing_. Springer, 173–186.
  * Mothilal et al. (2020) Ramaravind K Mothilal, Amit Sharma, and Chenhao Tan. 2020\.  Explaining machine learning classifiers through diverse counterfactual explanations. In _Proceedings of the 2020 Conference on Fairness, Accountability, and Transparency_. 607–617.
  * Pawelczyk et al. (2020a) Martin Pawelczyk, Klaus Broelemann, and Gjergji Kasneci. 2020a.  Learning Model-Agnostic Counterfactual Explanations for Tabular Data. In _Proceedings of The Web Conference 2020_. 3126–3132.
  * Pawelczyk et al. (2020b) Martin Pawelczyk, Klaus Broelemann, and Gjergji Kasneci. 2020b.  On counterfactual explanations under predictive multiplicity. In _Conference on Uncertainty in Artificial Intelligence_. PMLR, 809–818.
  * Pearl (2009) Judea Pearl. 2009\.  _Causality_ (2 ed.).  Cambridge University Press.  https://doi.org/10.1017/CBO9780511803161
  * Quinlan et al. (1987) J. R. Quinlan, P. J. Compton, K. A. Horn, and L. Lazarus. 1987\.  Inductive Knowledge Acquisition: A Case Study. In _Proceedings of the Second Australian Conference on Applications of Expert Systems_. Addison-Wesley Longman Publishing Co., Inc., USA, 137–156.
  * Russell (2019) Chris Russell. 2019\.  Efficient search for diverse coherent explanations. In _Proceedings of the Conference on Fairness, Accountability, and Transparency_. 20–28.
  * Schoeffer et al. (2022) Jakob Schoeffer, Niklas Kuehl, and Yvette Machowski. 2022\.  “There Is Not Enough Information”: On the Effects of Explanations on Perceptions of Informational Fairness and Trustworthiness in Automated Decision-Making. In _2022 ACM Conference on Fairness, Accountability, and Transparency_ (Seoul, Republic of Korea) _(FAccT ’22)_. Association for Computing Machinery, New York, NY, USA, 1616–1628.  https://doi.org/10.1145/3531146.3533218
  * Sharma et al. (2020) Shubham Sharma, Jette Henderson, and Joydeep Ghosh. 2020\.  _CERTIFAI: A Common Framework to Provide Explanations and Analyse the Fairness and Robustness of Black-Box Models_.  Association for Computing Machinery, New York, NY, USA, 166–172.  https://doi.org/10.1145/3375627.3375812
  * Shen et al. (2015) Jianying Shen, Jingying Zhang, Jing Wen, Qiang Ming, Ji Zhang, and Yawei Xu. 2015\.  Correlation of serum alanine aminotransferase and aspartate aminotransferase with coronary heart disease.  _International journal of clinical and experimental medicine_ 8, 3 (2015), 4399.
  * Shmerling (2020) Robert H. Shmerling, MD. 2020\.  How’s your heart rate and why it matters?  _Harvard Health_ (3 2020).  https://www.health.harvard.edu/heart-health/hows-your-heart-rate-and-why-it-matters
  * Slack et al. (2021) Dylan Slack, Anna Hilgard, Himabindu Lakkaraju, and Sameer Singh. 2021.  Counterfactual explanations can be manipulated.  _Advances in Neural Information Processing Systems_ 34 (2021), 62–75.
  * Sokol and Flach (2019) Kacper Sokol and Peter A Flach. 2019.  Counterfactual explanations of machine learning predictions: opportunities and challenges for AI safety. In _SafeAI@ AAAI_.
  * Tolomei et al. (2017) Gabriele Tolomei, Fabrizio Silvestri, Andrew Haines, and Mounia Lalmas. 2017.  Interpretable predictions of tree-based ensembles via actionable feature tweaking. In _Proceedings of the 23rd ACM SIGKDD international conference on knowledge discovery and data mining_. 465–474.
  * Upadhyay et al. (2021) Sohini Upadhyay, Shalmali Joshi, and Himabindu Lakkaraju. 2021\.  Towards robust and reliable algorithmic recourse.  _Advances in Neural Information Processing Systems_ 34 (2021), 16926–16937.
  * Ustun et al. (2019) Berk Ustun, Alexander Spangher, and Yang Liu. 2019\.  Actionable recourse in linear classification. In _Proceedings of the Conference on Fairness, Accountability, and Transparency_. 10–19.
  * Van Looveren and Klaise (2019) Arnaud Van Looveren and Janis Klaise. 2019.  Interpretable counterfactual explanations guided by prototypes.  _arXiv preprint arXiv:1907.02584_ (2019).
  * Verma et al. (2020) Sahil Verma, John Dickerson, and Keegan Hines. 2020\.  Counterfactual Explanations for Machine Learning: A Review.  _arXiv preprint arXiv:2010.10596_ (2020).
  * Virgolin and Fracaros (2023) Marco Virgolin and Saverio Fracaros. 2023.  On the Robustness of Sparse Counterfactual Explanations to Adverse Perturbations.  _Artif. Intell._ 316, C (feb 2023), 30 pages.  https://doi.org/10.1016/j.artint.2022.103840
  * von Kügelgen et al. (2022) Julius von Kügelgen, Amir-Hossein Karimi, Umang Bhatt, Isabel Valera, Adrian Weller, and Bernhard Schölkopf. 2022.  On the fairness of causal algorithmic recourse. In _Proceedings of the AAAI Conference on Artificial Intelligence_ , Vol. 36. 9584–9594.
  * Wachter et al. (2017) Sandra Wachter, Brent Mittelstadt, and Chris Russell. 2017\.  Counterfactual Explanations without Opening the Black Box: Automated Decisions and the GDPR.  arXiv:1711.00399 [cs.AI]
  * Wang et al. (2021) Yongjie Wang, Qinxu Ding, Ke Wang, Yue Liu, Xingyu Wu, Jinglong Wang, Yong Liu, and Chunyan Miao. 2021.  _The Skyline of Counterfactual Explanations for Machine Learning Decision Models_.  Association for Computing Machinery, New York, NY, USA, 2030–2039.  https://doi.org/10.1145/3459637.3482397


# Simple Quantum State Encodings
for Hybrid Programming of Quantum Simulators

Thomas Gabor, Marian Lingsch Rosenfeld, Claudia Linnhoff-Popien  LMU Munich

###### Abstract

Especially sparse quantum states can be efficiently encoded with simple
classical data structures. We show the admissibility of using a classical
database to encode quantum states for a few practical examples and argue in
favor of further optimizations for quantum simulation targeting simpler, only
“semi-quantum” circuits.

###### Index Terms:

quantum computing, quantum algorithms, quantum state

This is a pre-print of an article published in QSA 2022.

##  I Introduction

As of today, most sufficiently complex quantum algorithms are not actually run
on quantum hardware. Instead, quantum simulators have been known to offer more
qubit space and even faster computation times per financial resource spent –
and most obviously: fully controllable noise [7]. As literature proposes many
different approaches and suggestions on how to integrate future quantum
hardware with classical computation infrastructure, integrating inherently
classical quantum simulators with classical systems should be much easier.
However, it is clear why the developers of quantum simulators may opt to more
closely mimic the behavior of actual quantum computers in any aspect. Still,
in this paper we want to consider using the classical crutch to our advantage
while we still depend on it anyway. Effectively, we want to achieve two co-
dependent goals:

  1. (A)

Write quantum circuits that incorporate classical operations to a larger
extend than would be viable on a real quantum computer.

  2. (B)

Simulate circuits producing very simple, near-classical quantum states
efficiently within a quantum simulator that runs on classical hardware anyway

In this paper, we show almost trivial examples for (A) and showcase how the
simple approach of encoding quantum states as classical data structures as
long as possible enables us to realize (B) in this context. As such, this
paper should be read as a motivation to explore quantum algorithms beyond
their current physically realizable platform. Why should that be useful? The
success of quantum hardware will depend on quantum software and effective
quantum software is in desperate need of new algorithms and new methods that
yield algorithms with some kind of quantum advantage [2]. We reckon that
finding such algorithms is made easier by expanding the playground and
providing programming paradigms closer to what most non-quantum programmers
are used to. As of today, whether we like it or not, the question which
hardware platform our quantum algorithms will eventually be run on is
secondary.

##  II Quantum State Encodings

Qiskit is one of the most popular toolkits used to develop algorithms for
gate-based quantum computers [7]. As such, it comes with an integrated
simulator for quantum circuits that allows to test them on classical hardware.
Our argument is motivated by the observation that for certain kind of quantum
states very simple quantum state encodings perform much better at simulation
than the state-of-the-art Qiskit simulator. We quantify this performance in
Section III.

What we mean by “very simple state encodings” is using the most obvious
classical data structures to encode quantum states. For the purpose of this
paper we compare the following approaches empirically:

  1. (a)

Array. We encode the entire state as a standard Numpy array containing all
classical states occurring within the superposition and their respective
amplitudes. Note that we only include states with non-zero amplitudes this
way.

  2. (b)

Database. Similar to Array, we save all classical states with non-zero
amplitude within the quantum state. In this case, we use an sqlite3 database
to do so.

  3. (c)

Qiskit. As a benchmark, we perform all experiments using the state-of-the-art
quantum simulator provided with Qiskit.

  4. (d)

Mixed. We provide a simple heuristic to decide between using the Qiskit
simulator and the Database state encoding on a case-by-case basis, i.e., if
less than $\frac{2}{3}$ of the qubits have a Hadamard gate applied to them
somewhere within the circuit, the Database encoding is used, and if not, the
standard Qiskit simulator is used.

The main idea behind Array and Database encoding is enabling efficient
representation of very sparse quantum states. As most quantum algorithms try
to avoid these states (as they do not exploit the capabilities of quantum
hardware to the fullest extent), it is not surprising that standard quantum
simulators like the one provided by Qiskit are not entirely optimized for
these scenarios. In fact, using special encodings for rather sparse quantum
states is not a new idea [6]. However, we show that a state-of-the-art
simulator like Qiskit (as future work should naturally also consider more
quantum simulators) has a hard time recognizing near-classical states and
optimizing its treatment of them. While at first it may seem paradoxical to
use a quantum simulator to simulate almost or even fully classical algorithms,
we argue that having an all-in-one simulation framework opens up a new process
for hybrid programming, i.e., using quantum and classical parts within the
same algorithm, where we can transfer classical programs identically to a
quantum simulation framework to then embed quantum programming parts wherever
they seem beneficial. Thus, in contrast to most work on encoding quantum
states for quantum simulators (cf. [3, 1, 4], e.g.) we do not aim to improve
classical encodings for standard quantum algorithms (although that is a
valuable challenge, of course) but to reduce the overhead for already quasi-
classical states within the simulator.

In the following experiments, we also present another twist on Database
encoding called _state drop_ : We limit our database size to $1000$ entries by
simply removing the quantum states with the smallest amplitude after each
quantum (i.e., in our case Hadamard) operation. This approach was inspired by
the technique of approximate approximation [5], which works surprisingly well
for quantum annealing. Our results are more expected and the error rate
deteriorates quickly as our empirical results show. Still, for some
experiments we see a small middle ground where computation times can be
improved by admitting this approximation, perhaps making it worthwhile for
further investigation.

##  III Experiments

To test the quantum state encodings described in the previous section, we
execute various quantum algorithms. All code can be found in the accompanying
repository111https://github.com/marian-lingsch/qc-simulators.

###  III-A Generate Superposition

A Hadamard gate is applied to qubits $q_{i}$ for $i\in\\{1,...,r\\}$, where
for benchmarking purposes $r$ varied between $0$ and $n$ with $n$ being the
total amount of qubits. See Figure III-A for a depiction of the circuit.

Figure 1: Quantum circuit for the “generate superposition” benchmark.

Figure 2: Performance analysis for the “generate superposition” benchmark.
Identical 3D plot from two perspectives.

Figure 3: State drop performance for the “generate superposition” benchmark

We execute this benchmark using all quantum state encodings described in the
previous section. Note that we vary the overall size of the benchmark by
simply introducing more qubits in total, i.e., increasing $n$. We also vary
the “quantum-ness” of the circuit by applying the Hadamard gate only to some
(i.e., an amount of $r$) of the involved qubits (shown on the “Nondet qubits”
axis in Figure 2) as only these qubits’ values cannot be determined
classically. The Z axis shows the execution time in seconds on a log scale.
When running the experiments on our test machine was no longer feasible in
regard to computational resources, we cut off the experiments at the largest
possible number of total qubits. Note that with respect to the base surface
(i.e., “Total qubits”$\times$“Nondet qubits”) all plotted surfaces occupy the
upper triangle, i.e., “Nondet qubits” can never be larger than “Total qubits”.

It can be seen that the runtime of the database simulator is only dependent on
the amount of Nondet Qubits or $r$, while the runtime of the Qiskit simulator
depends on the Total Qubits or $n$, even when these add no real complexity to
the problem. We observe that for few Nondet Qubits the Database simulator
performs better than the Qiskit simulator, while the Qiskit simulator is
better than the Database simulator when the amount of Nondet Qubits is close
to the amount of Total Qubits; in fact Qiskit is the only computationally
feasible encoding for more than $27$ “Nondet qubits”. The Mixed simulator
allows us to have the best of both worlds in all instances. It can be seen
that the Array simulator performs much worse than the Database simulator,
except for very small inputs. For this reason, it was not considered in
further experiments.

Figure 3 shows the results for the state drop technique. We fixed the total
number of qubits to $20$ and compared computation times and produced error for
both the standard Database simulator as well as Database with state drop. Note
that the pure Database error rate is always $0$ so it is not plotted
separately. In this scenario, we see that computation times of state drop
compared to the full Database approach only improve in a regime where error
rates have already reached infeasible heights.

###  III-B Addition

We add two input variables with $k$ bits each. For benchmarking purposes, we
apply Hadamard gates to $r$, $r\in\\{0,...,2k-1\\}$, of these qubits. Qubits
$q_{0}$ through $q_{2k-1}$ represent the input variables and $q_{2k}$ through
$q_{3k-1}$ represent the output variable while $q_{3k}$ through $q_{3k+4}$ are
ancillary qubits. The total amount of qubits is $3k+5$. Figure III-B depicts
the whole circuit.

Figure 4: Quantum circuit for the “addition” benchmark.

Figure 5: Performance analysis for the “addition” benchmark. Identical 3D plot
from two perspectives.

Figure 6: State drop performance for the “addition” benchmark

The analysis in Figure 5 shows that the Qiskit simulator performs worse than
the Database simulator most of the time. This occurs because there are much
fewer states with non-zero probability than the total possible amount, given
the mostly classical nature of the circuit. The Mixed simulator always chooses
to use the Database simulator and thus their runtimes are very similar. For
more than $29$ total qubits, Qiskit could no longer be run practically on our
machine while Database (and thus Mixed) still work up to $30$ total qubits.

Figure 6 shows the state drop trade-off for $26$ total qubits. At $11$ “Nondet
qubits” we already have a huge normalized error and errors quickly climb to
maximum beyond that. However, compared to Figure 3, there exists a trade-off
at least, when with $11$, $12$, or $13$ “Nondet qubits” we save computation
time (albeit by introducing great amounts of error), which could be further
investigated in future work.

###  III-C Grover’s Search

We implement Grover’s algorithm using an oracle for $\mathbf{0}$, i.e., the
state where all the numbers are zero. The circuit used to implement this
algorithm is given in Figure III-C, where $U$ is the oracle, $A$ is the
addition operator as described in Subsection III-B, and $D$ is the diffusion
operator.

Figure 7: Quantum circuit for the “Grover’s search” benchmark

For the Qiskit simulator the diffusion operator is implemented as a series of
Hadamard and CNOT gates, while for the Database simulator we took a shortcut
by implementing the effect of the diffusion operator within a single SQL query
to the database to showcase the potential of classical approximation for
practical simulation.

Figure 8: Performance analysis for the “Grover’s search” benchmark.

Figure 9: State drop performance for the “Grover’s search” benchmark

Since all qubits in the search range need to be non-deterministic for Grover’s
algorithm, Figure 8 shows a two-dimensional plot this time. As the oracle
implements addition, the Mixed simulator’s behavior carries over from the
previous benchmark (cf. Subsection III-B) and it always resorts to the
Database simulator, which performs better than the Qiskit simulator. It is
interesting to observe that behavior even when the whole framework of Grover’s
search, one of the primary examples for quantum algorithms, is put around the
benchmark of addition.

Figure 9 shows a discouraging image where state drop errors rise early on and
we again only see a small regime at $11$ “Nondet qubits” (i.e., where the
“Nondet qubits” span over $2^{11}=2048$ possible states) where we can achieve
a computational benefit by accepting relatively large error rates.

##  IV Conclusion

We have analyzed four basic techniques for classically encoding quantum states
and have seen that even very simple mechanisms can beat the state-of-the-art
Qiskit simulator, when the circuit considered do not need to construct too
complex quantum states. We have shown a very simple “Mixed” heuristic that can
switch between these scenarios dynamically. Future work needs to improve that
heuristic for more complex real-world scenarios and perhaps integrate many
more means of quantum state simulation. Furthermore, the “state drop”
technique was first evaluated with limited success. However, we reckon that
more advanced heuristics might still be able to more precisely hit a sweet
spot between computational resources and acceptable error that might just push
classical quantum simulation a step further.

Overall, we think it is big challenge for the field to come up with a fully
integrated simulator that employs various techniques to recognize the circuit
it is given and that can be easily parametrized to yield an arbitrary
error/time trade. We should also consider how much more complex algorithms
might be split up and be simulated using a combination of various techniques.
We suspect that such a generic platform for quantum software development will
enable more complex applications for quantum computers by allowing developers
to start out with classical algorithms as they would use on classical hardware
and migrate gradually towards more usage of quantum states and thus quantum
hardware. Such tools could lead to a truly hybrid workflow for programming
quantum computers to be used within hybrid setups for real-world challenges.

## References

  * [1] Itai Arad and Zeph Landau.  Quantum computation and the evaluation of tensor networks.  SIAM Journal on Computing, 39(7):3089–3121, 2010.
  * [2] Thomas Gabor, Leo Sünkel, Fabian Ritz, Thomy Phan, Lenz Belzner, Christoph Roch, Sebastian Feld, and Claudia Linnhoff-Popien.  The holy grail of quantum artificial intelligence: Major challenges in accelerating the machine learning pipeline.  In Proceedings of the IEEE/ACM 42nd International Conference on Software Engineering Workshops, pages 456–461, 2020.
  * [3] Richard Jozsa.  On the simulation of quantum circuits.  arXiv preprint quant-ph/0603163, 2006.
  * [4] Nader Khammassi, Imran Ashraf, Xiang Fu, Carmen G Almudever, and Koen Bertels.  QX: A high-performance quantum computer simulation platform.  In Design, Automation & Test in Europe Conference & Exhibition (DATE), 2017, pages 464–469. IEEE, 2017.
  * [5] Irmi Sax, Sebastian Feld, Sebastian Zielinski, Thomas Gabor, Claudia Linnhoff-Popien, and Wolfgang Mauerer.  Approximate approximation on a quantum annealer.  In Proceedings of the 17th ACM International Conference on Computing Frontiers, pages 108–117, 2020.
  * [6] Guifré Vidal.  Efficient classical simulation of slightly entangled quantum computations.  Physical review letters, 91(14):147902, 2003.
  * [7] Robert Wille, Rod Van Meter, and Yehuda Naveh.  IBM’s Qiskit tool chain: Working with and developing for real quantum computers.  In 2019 Design, Automation & Test in Europe Conference & Exhibition (DATE), pages 1234–1240. IEEE, 2019.