# Rethinking Counterfactual Explanations as Local and Regional Counterfactual
Policies

Salim I. Amoukou
LaMME
University Paris Saclay
Stellantis Paris
&Nicolas J-B. Brunel
LaMME
ENSIIE, University Paris Saclay
Quantmetry Paris

###### Abstract

Among the challenges not yet resolved for Counterfactual Explanations (CE),
there are stability, synthesis of the various CE and the lack of
plausibility/sparsity guarantees. From a more practical point of view, recent
studies show that the prescribed counterfactual recourses are often not
implemented exactly by the individuals and demonstrate that most state-of-the-
art CE algorithms are very likely to fail in this noisy environment. To
address these issues, we propose a probabilistic framework that gives a sparse
local counterfactual rule for each observation: we provide rules that give a
range of values that can change the decision with a given high probability
instead of giving diverse CE. In addition, the recourses derived from these
rules are robust by construction. These local rules are aggregated into a
regional counterfactual rule to ensure the stability of the counterfactual
explanations across observations. Our local and regional rules guarantee that
the recourses are faithful to the data distribution because our rules use a
consistent estimator of the probabilities of changing the decision based on a
Random Forest. In addition, these probabilities give interpretable and sparse
rules as we select the smallest set of variables having a given probability of
changing the decision. Codes for computing our counterfactual rules are
available, and we compare their relevancy with standard CE and recent similar
attempts.

##  1 Introduction

In recent years, many explanations methods have been developed for explaining
machine learning models, with a strong focus on local analysis, i.e.,
generating explanations for individual prediction (see (Molnar, 2022) for a
survey). Among this plethora of methods, one of the most prominent and active
techniques are Counterfactual Explanations (Wachter et al., 2017b). Unlike
popular local attribution methods, e.g., SHAP (Lundberg et al., 2020) and LIME
(Ribeiro et al., 2016), which highlight the importance score of each feature,
Counterfactuals Explanations (CE) describe the smallest modification to the
feature values that changes the prediction to a desired target. Although CE
are intuitive and user-friendly by giving recourse in some scenarios (e.g.,
loan application), they have many shortcomings in practice. Indeed, most
counterfactual methods rely on a gradient-based algorithm or heuristics
approaches (Karimi et al., 2020b), thus can fail to identify the most natural
explanations and lack guarantees. Most algorithms either do not guarantee
sparse counterfactuals (changes in the smallest number of features) or do not
generate in-distribution samples (see (Verma et al., 2020; Chou et al., 2022)
for a survey on counterfactuals methods). Although some works (Parmentier and
Vidal, 2021; Poyiadzi et al., 2019; Looveren and Klaise, 2019) try to solve
the plausibility/sparsity problem, the suggested solutions are not entirely
satisfactory.
In another direction, many papers (Mothilal et al., 2020; Karimi et al.,
2020a; Russell, 2019) encourages the generation of diverse counterfactuals in
order to find actionable recourse (Ustun et al., 2019). Actionability is a
vital desideratum, as some features may be non-actionable, and generating many
counterfactuals increases the chance of getting actionable recourse. However,
the diversity of CE makes the explanations less intelligible, and the
synthesis of various CE or local explanations, in general, is yet to be
comprehensively solved (Lakkaraju et al., 2022). In addition, recently
Pawelczyk et al. (2022) highlights a new problem of local CE called: noisy
responses to prescribed recourses. Indeed, in real-world scenarios, some
individuals may not be able to implement exactly the prescribed recourses, and
they show that most CE methods fail in this noisy environment. Therefore, we
propose to reverse the usual way of explaining with counterfactual by
computing Counterfactual rules. We introduce a new line of counterfactuals: we
build interpretable policies for changing a decision with a given probability
that ensure the stability of the deduced recourse. These policies are optimal
(in sparsity) and faithful to the data distribution. Their computation comes
with statistical guarantees as they use a consistent estimator of the
conditional distribution. Our proposal is to find a general policy or rule
that permits changing the decision while fixing some features instead of
generating many counterfactual samples. One of the main challenges is to
identify the (minimal) set of features that provide the best promising
directions for changing the decision to the desired output. We also show this
approach can be extended for finding a collection of regional counterfactuals,
such that we have a global counterfactual policy for analyzing a model. An
example of the counterfactual rules that we introduce is given in figure 1.

Figure 1: Illustration of the local and regional Counterfactuals Rules that we
introduced on a dataset with 4 variables: Age, Salary, Sex, and HoursPerWeek.
The Counterfactual Rules define intervals on the minimal subset of features to
change the decision of a model prediction in the local counterfactual rule or
the decision of a rule that applies on a sub-population in the regional
counterfactual rule. In Blue, we have the proposed rules to change the
decision.

##  2 Motivation and Related works

Most of the methods that propose Counterfactuals Explanations are based on the
approach of the seminal work of Wachter et al. (2017a): the counterfactuals
are generated by optimizing a cost, but this procedure does not account
directly the plausibility of the counterfactual examples (see (Verma et al.,
2020) for classification of CE methods). Indeed, a major shortcoming is that
the adverse decision needed for obtaining the counterfactual is not designed
to be feasible or representative of the underlying data distribution. However,
some recent studies proposed ad-hoc plausibility constraint in the
optimization, using for instance an outlier score (Kanamori et al., 2020), an
Isolation Forest (Parmentier and Vidal, 2021) or a density-weighted metrics
(Poyiadzi et al., 2019) to generate in-distribution samples. In another
direction, Looveren and Klaise (2019) proposes to use an autoencoder that
penalizes out-of-distribution candidates. Instead of relying on ad-hoc
constraints, we propose CE that gives plausible explanations by design.
Indeed, for each observation, we identify the variables and associated ranges
of values that have the highest probability of changing the prediction. We can
compute this probability with a consistent estimator of the conditional
distribution $P(Y|\boldsymbol{X}_{S})$. As a consequence, the sparsity of the
counterfactuals is not encouraged indirectly by adding a penalty term
($\ell_{0}$ or $\ell_{1}$) as existing works (Mothilal et al., 2020). Our
approach is inspired by the concept of Same Decision Probability (SDP)
(introduced in (Chen et al., 2012)) that can be used for identifying the
smallest subset of features to guarantee (with a given probability) the
stability of a prediction. This minimal subset is called Sufficient
Explanations. In (Amoukou and Brunel, 2021), it has been shown that the SDP
and the Sufficient Explanations can be estimated and computed efficiently for
identifying important local variables in any classification and regression
models. For counterfactuals, we are interested in the dual set: we want the
minimal subset of features that have a high probability of changing the
decision (when the other features are fixed). Another limitation of the
current CE is their local nature and the multiplicity of the explanations
produced. While some papers (Mothilal et al., 2020; Karimi et al., 2020a;
Russell, 2019) promote the generation of diverse counterfactual samples to
ensure actionable recourse, such diverse explanations should be summarized to
be intelligible (Lakkaraju et al., 2022), but the compilation of local
explanations is often a very difficult problem. To address this problem, we do
not generate counterfactual samples, but we build a rule Counterfactual Rules
(CR) from which we can derive counterfactuals. Contrary to classic CE which
gives the nearest instances with a desired output, we find the most effective
rule for each observation (or group of similar observations) that changes the
prediction to the desired target. This local rule easily aggregates similar
counterfactuals. For example, if $\boldsymbol{x}=\\{\texttt{Age=20,
Salary=35k, HoursWeek=25h, Sex=M},\dots\\}$ with Loan=False, fixing the
variables Age and Sex and changing the Salary and HoursWeek change the
decision. Therefore, instead of given multiples combination of Salary and
HoursWeek (e.g. 35k and 40h or 40k and 55h, …) that result in many instances,
the counterfactual rule gives the range of values: IF HoursWeek
$\in\texttt{[35h, 50h]}$, Salary $\in$ [40k, 50k], and the remaining features
are fixed THEN Loan=True. It can be extended at a regional scale, e.g., given
a rule $\textbf{R}=\\{\texttt{IF Salary}\in\texttt{[35k, 20k],
Age}\in\texttt{[20, 80] THEN Loan=False}\\}$, the regional Counterfactual Rule
(CR) could be $\\{\texttt{{IF } Salary}\in\texttt{[40k, 50k],}\texttt{
HoursWeek}\in\texttt{[35h, 50h] and the}\textbf{ remaining rules are
fixed}\texttt{ THEN Loan=True}\\}$. The main difference between a local and a
global CR is that the Local-CR explain a single instance by fixing the
remaining feature values (not used in the CR) ; while a regional-CR is defined
by keeping the remaining variables in a given interval (not used in the
regional-CR). Moreover, by giving ranges of values that guarantee a high
probability of changing the decision, we partly answer the problem of noisy
responses to prescribed recourses (Pawelczyk et al., 2022) so long as the
perturbations are within our ranges.

Although the Local Counterfactual Rule is new, the Regional Counterfactual
Rule can be related to some recent works. Indeed, Rawal and Lakkaraju (2020)
proposed Actionable Recourse Summaries (AReS), a framework that constructs
global counterfactual recourses in order to have a global insight of the model
and detect unfair behavior. While AReS is similar to the Regional
Counterfactual Rule, we emphasize some significant differences. Our methods
can address regression problems and deal with continuous features. Indeed,
AReS needs to discretize the continuous features, inducing a trade-off between
speed and performance as noticed by (Ley et al., 2022). Thus, too few bins
result in unrealistic recourse, while too many bins result in excessive
computation time. In addition, AReS uses a greedy heuristic search approach to
find global recourse, which might produce sub-optimal recourse. As we have
already mentioned, the changes we provide overcome these two limitations
because the consistency of our counterfactual is controlled by an estimation
of the probability of changing the decision, and because we favor changes of a
minimum number of features. Another global CE framework has been introduced in
(Kanamori et al., 2022) to ensure transparency: the Counterfactual Explanation
Tree (CET) partitions the input space with a decision tree and assigns an
appropriate action for changing the decision of each subspace. Therefore, it
gives a unique action for changing the decision of multiple instances. In our
case, we offer more flexibility in the counterfactual explanations because we
provide a range of possible values that guarantee a change with a given
probability. In our approach, we do not make any assumption about the cost of
changing the feature nor the causal structure. If we have such information,
then we can add it as additional post-processing such that it can be made more
explicit and more transparent for the final user as required for trustworthy
AI.

##  3 Minimal Counterfactual Rules

We assume that we have an i.i.d sample
$\mathcal{D}_{n}=\\{(\boldsymbol{X}_{i},Y_{i})_{i=1,\dots,n}\\}$ such that
$(\boldsymbol{X},Y)\sim P_{(\boldsymbol{X},Y)}$ where
$\boldsymbol{X}\in\mathcal{X}$ (typically $\mathcal{X}=\mathbb{R}^{p}$) and
$Y\in\mathcal{Y}$. The output $\mathcal{Y}$ can be discrete or continuous. We
want to explain the predictor $f:\mathbb{R}^{p}\mapsto\mathcal{Y}$, that has
been learned with the dataset $\mathcal{D}_{n}$. We use uppercase letters for
random variables and lowercase letters for their value assignments. For a
given subset $S\subset[p]$, $\boldsymbol{X}_{S}=(X_{i})_{i\in S}$ denotes a
subgroup of features, and we write
$\boldsymbol{x}=(\boldsymbol{x}_{S},\boldsymbol{x}_{\bar{S}})$ (with some
abuse of notation).

For an observation $(\boldsymbol{x},y=f(\boldsymbol{x}))$, we have a target
set $\mathscr{Y}^{\star}\subset\mathcal{Y}$, such that
$y\notin\mathscr{Y}^{\star}$. For the simple case of classification problem,
$\mathscr{Y}^{\star}=\\{y^{\star}\\}$ is the standard singleton such that
$y^{\star}\in\mathcal{Y}$ is different of $y$. Contrary to standard
approaches, our definition of the counterfactual deals also with the
regression case by considering $\mathscr{Y}^{\star}=[a,b]\subset\mathbb{R}$;
our definitions and computations of counterfactuals are the same for both
classification and regression. We remind that the classic CE problem (defined
only for classification) is to find a function
$\boldsymbol{a}:\mathcal{X}\mapsto\mathcal{X}$, such that for all observations
$\boldsymbol{x}\in\mathcal{X}$, $f(\boldsymbol{x})\neq y^{\star}$, and we have
$f(\boldsymbol{a}(\boldsymbol{x}))=y^{\star}$. With standard CE, the function
is defined point-wise by solving an optimisation program. Most often
$\boldsymbol{a}(\cdot)$ is not a real function, as $\boldsymbol{a}(x)$ may be
in fact a collection of (random) values
$\\{\boldsymbol{x}_{1}^{\star},\dots,\boldsymbol{x}_{p}^{\star}\\}$. A more
recent point of view was proposed by Kanamori et al. (2022), and it defines
$\boldsymbol{a}$ as a decision tree, where in each leaf $L$, the best
perturbation $a_{L}$ is predicted and add it to all the instances
$\boldsymbol{x}\in L$.
Our approach is hybrid, because we do not propose a single action for each
subspace of $\mathcal{X}$ or sub-group of population, but we give sets of
possible perturbations. Indeed, a _Local Counterfactual Rule_ (Local-CR) for
$\mathscr{Y}^{\star}$ and observation $\boldsymbol{x}$ (with
$f(\boldsymbol{x})\notin\mathscr{Y}^{\star}$) is a rectangle
$C_{S}(\boldsymbol{x};\mathscr{Y}^{\star})=\prod_{i\in
S}[a_{i},b_{i}],a_{i},b_{i}\in\overline{\mathbb{R}}$ such that for all
perturbations of
$\boldsymbol{x}=\left(\boldsymbol{x}_{S},\boldsymbol{x}_{\bar{S}}\right)$
obtained as
$\boldsymbol{x}^{\star}=\left(\boldsymbol{z}_{S},\boldsymbol{x}_{\bar{S}}\right)$
with $\boldsymbol{z}_{S}\in C_{S}(\boldsymbol{x};\mathscr{Y}^{\star})$ and
$\boldsymbol{x}^{\star}$ an in-distribution sample, then
$f\left(\boldsymbol{x}^{\star}\right)$ is in $\mathscr{Y}^{\star}$ with a high
probability.
Similarly, a _Regional Counterfactual Rule_ (Regional-CR)
$C_{S}(\boldsymbol{R};\mathscr{Y}^{\star})$ is defined for
$\mathscr{Y}^{\star}$ and a rectangle
$\boldsymbol{R}=\prod_{i=1}^{d}[a_{i},b_{i}],a_{i},b_{i}\in\overline{\mathbb{R}}$,
if for all observations
$\boldsymbol{x}=(\boldsymbol{x}_{S},\boldsymbol{x}_{\bar{S}})\in\boldsymbol{R}$,
the perturbations obtained as
$\boldsymbol{x}^{\star}=(\boldsymbol{z}_{S},\boldsymbol{x}_{\bar{S}})$ with
$\boldsymbol{z}_{S}\in C_{S}(\boldsymbol{R},\mathscr{Y}^{\star})$ and
$\boldsymbol{x}^{\star}$ an in-distribution sample are such that
$f\left(\boldsymbol{x}^{\star}\right)$ is in $\mathscr{Y}^{\star}$ with high
probability.
We build such rectangles sequentially, first, we propose to find the best
directions $S\subset[p]$ that offers the best probability of change. Then, we
find the best intervals $[a_{i},b_{i}],i\in S$ that change the decision to the
desired target. A central tool in this approach is the Counterfactual Decision
Probability.

######  Definition 3.1.

Counterfactual Decision Probability (CDP). The Counterfactual Decision
Probability of the subset $S\subset\left\llbracket 1,p\right\rrbracket$, w.r.t
$\boldsymbol{x}=\left(\boldsymbol{x}_{S},\boldsymbol{x}_{\bar{S}}\right)$ and
the desired target $\mathscr{Y}^{\star}$ (s.t.
$f(\boldsymbol{x})\notin\mathscr{Y}^{\star})$ is

$CDP_{S}\left(\mathscr{Y}^{\star};\boldsymbol{x}\right)=P\left(f(\boldsymbol{X})\in\mathscr{Y}^{\star}\left|\boldsymbol{X}_{\bar{S}}=\boldsymbol{x}_{\bar{S}}\right.\right).$

The $CDP$ of the subset S is the probability that the decision changes to the
desired target $\mathscr{Y}^{\star}$ by sampling the features
$\boldsymbol{X}_{S}$ given
$\boldsymbol{X}_{\bar{S}}=\boldsymbol{x}_{\bar{S}}$. It is related to the Same
Decision Probability
$SDP_{S}(\mathscr{Y};\boldsymbol{x})=P\left(f(\boldsymbol{X})\in\mathscr{Y}|\boldsymbol{X}_{S}=\boldsymbol{x}_{S}\right)$
used in (Amoukou and Brunel, 2021) for solving the dual problem of selecting
the most local important variables for obtaining and maintaining the decision
$f(\boldsymbol{x})\in\mathscr{Y}$ (where
$f(\boldsymbol{x})\in\mathscr{Y}\subset\mathcal{Y}$). The set $S$ is called
the Minimal Sufficient Explanation. Indeed, we have
$CDP_{S}(\mathscr{Y}^{\star};\boldsymbol{x})=SDP_{\bar{S}}(\mathscr{Y}^{\star};\boldsymbol{x})$.
The computation of these probabilities is challenging and discussed in Section
4. We now focus on the minimal subset of features $S$ such that the model
makes the desired decision with a given probability $\pi$.

######  Definition 3.2.

( Minimal Divergent Explanations). Given an instance $\boldsymbol{x}$ and a
desired target $\mathscr{Y}^{\star}$, $S$ is a Divergent Explanation for
probability $\pi>0$, if
$CDP_{S}\left(\mathscr{Y}^{\star};\boldsymbol{x}\right)\geq\pi$, and no subset
$Z$ of $S$ satisfies
$CDP_{Z}\left(\mathscr{Y}^{\star};\boldsymbol{x}\right)\geq\pi$. Hence, a
Minimal Divergent Explanation is a Divergent Explanation with minimal size.

The set minimizing this probability is not unique, and we can have several
Minimal Divergent Explanations. Note that the probability $\pi$ represents the
minimum level required for a set to be chosen for generating counterfactuals,
and its value should be as high as possible and depends on the use case. We
have now enough material to define our main criterion for building a Local
Counterfactual Rule (Local-CR):

######  Definition 3.3.

(Local Counterfactual Rule). Given an instance $\boldsymbol{x}$, a desired
target $\mathscr{Y}^{\star}\not\owns f(\boldsymbol{x})$ , a Minimal Divergent
Explanation $S$, the rectangle
$C_{S}(\boldsymbol{x};\mathscr{Y}^{\star})=\prod_{i\in
S}[a_{i},b_{i}],a_{i},b_{i}\in\overline{\mathbb{R}}$ is a Local Counterfactual
Rule with probability $\pi_{C}$ if

$\displaystyle
CRP_{S}(\mathscr{Y}^{\star},\boldsymbol{x},C_{S}(\boldsymbol{x};y^{\star}))\triangleq
P(f(\boldsymbol{X})\in\mathscr{Y}^{\star}\;|\boldsymbol{X}_{S}\in
C_{S}(\boldsymbol{x};\mathscr{Y}^{\star}),\boldsymbol{X}_{\bar{S}}=\boldsymbol{x}_{\bar{S}})\geq\pi_{C}.$
(3.1)

The $CRP_{S}$ is the Counterfactual Rule Probability.

The higher the probability $\pi_{C}$ is, the better the relevance of the rule
$C_{S}(\boldsymbol{x};\mathscr{Y}^{\star})$ is, for this instance. Given a set
$S$, we seek for the maximal rectangle in the direction $S$ satisfying
Definition 3.1.

In practice, we can observe that the Local-CR
$C_{S}(\cdot;\mathscr{Y}^{\star})$ for neighbors
$\boldsymbol{x},\boldsymbol{x}^{\prime}$ are often quite close, because the
Minimal Divergent Explanations are similar and the corresponding rectangles
often overlaps. Hence, this motivates a generalisation of these Local-CR to
hyperrectangle
$\boldsymbol{R}=\prod_{i=1}^{d}[a_{i},b_{i}],a_{i},b_{i}\in\overline{\mathbb{R}}$
regrouping similar observations. We denote
$\text{supp}(\boldsymbol{R})=\\{i:[a_{i},b_{i}]\neq\overline{\mathbb{R}}\\}$
the support of the rectangle, and we extend the Local-CR to Regional
Counterfactual Rules (Regional-CR). In order to do it, we denote
$\boldsymbol{R}_{\bar{S}}=\prod_{i\in\bar{S}}[a_{i},b_{i}]$ as the rectangle
with intervals of $\boldsymbol{R}$ in $\text{supp}(\boldsymbol{R})\cap\bar{S}$
and we also defines the corresponding Counterfactual Decision Probability CDP
(Definition 3.1) for rule $\boldsymbol{R}$ and subset $S$ as
$CDP_{S}(\mathscr{Y}^{\star};\boldsymbol{R})=P\left(f(\boldsymbol{X})\in\mathscr{Y}^{\star}\left|\boldsymbol{X}_{\bar{S}}\in\boldsymbol{R}_{\bar{S}}\right.\right)$.
Therefore, we can also compute the Minimal Divergent Explanation for rule
$\boldsymbol{R}$ using Definition 3.2 with the CDP for rules.

######  Definition 3.4.

(Regional Counterfactual Rule). Given any rectangle $\boldsymbol{R}$, a
desired target $\mathscr{Y}^{\star}$, a Minimal Divergent Explanation $S$ of
$R$, the rectangle $C_{S}(\boldsymbol{R};y^{\star})=\prod_{i\in
S}[a_{i},b_{i}]$ is a Regional Counterfactual Rule with probability $\pi_{C}$
if

$\displaystyle
CRP_{S}(\mathscr{Y}^{\star};\boldsymbol{R},C_{S}(\boldsymbol{R},\mathscr{Y}^{\star}))\triangleq
P(f(\boldsymbol{X})\in\mathscr{Y}^{\star}\;|\boldsymbol{X}_{S}\in
C_{S}(\boldsymbol{R},\mathscr{Y}^{\star}),\boldsymbol{X}_{\bar{S}}\in\boldsymbol{R}_{\bar{S}})\geq\pi_{C}.$
(3.2)

$CRP_{S}(\mathscr{Y}^{\star};\boldsymbol{R},C_{S}(\boldsymbol{R}))$ is the
corresponding Counterfactual Rule Probability for rule $\boldsymbol{R}$.

#### Remarks:

Local-CR and regional-CR differ slightly: for local, we condition by
$\boldsymbol{X}_{\bar{S}}=\boldsymbol{x}_{\bar{S}}$ in Eq. 3.1, while for
regional, we condition by
$\boldsymbol{X}_{\bar{S}}\in\boldsymbol{R}_{\bar{S}}$. For computing regional-
CR, we can start for a rectangle generated by any method, such as (Wang et
al., 2017; Lin et al., 2020). The only condition is that it contains a
homogeneous group, i.e. with almost the same output. However, by default we
use as initial rules the Sufficient Rules derived in (Amoukou and Brunel,
2021) as it handles regression problem. The Sufficient Rules are minimal
support rectangles define for a given output $\mathscr{Y}$ as
$C_{S}(\mathscr{Y})=\Pi_{i\in S}[a_{i},b_{i}]$ such that
$\forall\boldsymbol{x}\in\mathcal{X},\boldsymbol{x}_{S}\in
C_{S}(\mathscr{Y})$,
$P(f(\boldsymbol{X})\in\mathscr{Y}|\boldsymbol{X}_{S}=\boldsymbol{x}_{S})\geq\pi$.

##  4 Estimation of the $CDP$ and $CRP$

In order to compute the probabilities $CDP_{S}$ and $CRP_{S}$ for any $S$, we
use a dedicated Random Forest (RF) $m_{k,n}$ that learns the model $f$ to
explain. Indeed, the conditional probabilities $CDP_{S}$ and $CRP_{S}$ can be
easily computed from a RF by combining the Projected Forest algorithm (Bénard
et al., 2021a) and the Quantile Regression Forest (Meinshausen and Ridgeway,
2006): hence we can estimate consistently the probabilities
$CDP_{S}(\mathscr{Y}^{\star};\boldsymbol{x})$. We adapt the approach used in
(Amoukou and Brunel, 2021) and remind for the sake of completeness, the
computation of the estimate of $SDP_{S}$.

###  4.1 Projected Forest and $CDP_{S}$

The estimator of the $SDP_{S}$ is built upon a learned Random Forest (Breiman
et al., 1984). A Random Forest (RF) is a predictor consisting of a collection
of $k$ randomized trees (see (Loh, 2011) for a detailed description of
decision tree). For each instance $\boldsymbol{x}$, the predicted value of the
$j$-th tree is denoted $m_{n}(\boldsymbol{x},\Theta_{j})$ where $\Theta_{j}$
represents the resampling data mechanism in the $j$-th tree and the successive
random splitting directions. The trees are then averaged to give the
prediction of the forest as:

$\displaystyle\small
m_{k,n}(\boldsymbol{x},\Theta_{1:k},\mathcal{D}_{n})=\frac{1}{k}\sum_{l=1}^{k}m_{n}(\boldsymbol{x};\Theta_{l},\mathcal{D}_{n})$
(4.1)

However, the RF can also be view as an adaptive nearest neighbor predictor.
For every instance $\boldsymbol{x}$, the observations in $\mathcal{D}_{n}$ are
weighted by $w_{n,i}(\boldsymbol{x};\Theta_{1:k},\mathcal{D}_{n})$,
$i=1,\dots,n$. Therefore, the prediction of RF can be rewritten as

$\small
m_{k,n}(\boldsymbol{x},\Theta_{1:k},\mathcal{D}_{n})=\sum_{i=1}^{n}w_{n,i}(\boldsymbol{x};\Theta_{1:k},\mathcal{D}_{n})Y_{i}.$

This emphasizes the central role played by the weights in the RF’s algorithm,
see (Meinshausen and Ridgeway, 2006; Amoukou and Brunel, 2021) for detailed
description of the weights. Therefore, it naturally gives estimators of other
quantities e.g., Cumulative hazard function (Ishwaran et al., 2008), Treatment
effect (Wager and Athey, 2017), conditional density (Du et al., 2021). For
instance, Meinshausen and Ridgeway (2006) showed that we can used the same
weights to estimate the Conditional Distribution Function with the following
estimator:

$\displaystyle\widehat{F}(y|\boldsymbol{X}=\boldsymbol{x},\Theta_{1:k},\mathcal{D}_{n})=\sum_{i=1}^{n}w_{n,i}(\boldsymbol{x};\Theta_{1:k},\mathcal{D}_{n})\mathds{1}_{Y_{i}\leq
y}$ (4.2)

In another direction, Bénard et al. (2021a) introduced the Projected Forest
algorithm (Bénard et al., 2021c, a) that aims to estimate
$E[Y|\boldsymbol{X}_{S}]$ by modifying the RF’s prediction algorithm.

#### Projected Forest:

To estimate $E[Y|\boldsymbol{X}_{S}=\boldsymbol{x}_{S}]$ instead of
$E[Y|\boldsymbol{X}=\boldsymbol{x}]$ using a RF, Bénard et al. (2021b)
suggests to simply ignore the splits based on the variables not contained in
$S$ from the tree predictions. More formally, it consists of projecting the
partition of each tree of the forest on the subspace spanned by the variables
in S. The authors also introduced an algorithmic trick that computes the
output of the Projected Forest efficiently without modifying the initial tree
structures. We drop the observations down in the initial trees, ignoring the
splits which use a variable not in $S$: when it encounters a split involving a
variable $i\notin S$, the observations are sent both to the left and right
children nodes. Therefore, each instance falls in multiple terminal leaves of
the tree. To compute the prediction of $\boldsymbol{x}_{S}$, we follow the
same procedure, and gather the set of terminal leaves where
$\boldsymbol{x}_{S}$ falls. Next, we collect the training observations which
belong to every terminal leaf of this collection, in other words, we keep only
the observations that fall in the intersection of the leaves where
$\boldsymbol{x}_{S}$ falls. Finally, we average their outputs $Y_{i}$ to
generate the estimation of $E[Y|\boldsymbol{X}_{S}=\boldsymbol{x}_{S}]$.
Notice that the author show that this algorithm converges asymptotically to
the true projected conditional expectation
$E[Y|\boldsymbol{X}_{S}=\boldsymbol{x}_{S}]$.

As the RF, the PRF gives also a weight to each observation. The associated PRF
is denoted
$m_{k,n}^{(\boldsymbol{x}_{S})}(\boldsymbol{x}_{S})=\sum_{i=1}^{n}w_{n,i}(\boldsymbol{x}_{S})Y_{i}$.
Therefore, as the weights of the original forest was used to estimate the CDF
in equation 4.2, Amoukou and Brunel (2021) used the weights of the Projected
Forest Algorithm to estimate the $SDP$ as
$\widehat{SDP}_{S}\left(\mathscr{Y};\boldsymbol{x}\right)=\sum_{i=1}^{n}w_{n,i}(\boldsymbol{x}_{S})\mathds{1}_{Y_{i}\in\mathscr{Y}}$.
The idea is essentially to replace $Y_{i}$ by
$\mathds{1}_{Y_{i}\in\mathscr{Y}}$ in the Projected Forest equation defined
above. The authors also show that this estimator converges asymptotically to
the true $SDP_{S}$. Therefore, we can estimate the $CDP$ with the following
estimator

$\widehat{CDP}_{S}\left(\mathscr{Y}^{\star};\boldsymbol{x}\right)=\sum_{i=1}^{n}w_{n,i}(\boldsymbol{x}_{\bar{S}})\mathds{1}_{Y_{i}\in\mathscr{Y}^{\star}}.$
(4.3)

#### Remarks:

Note that we only give the estimator of the $CDP_{S}$ of an instance
$\boldsymbol{x}$. The estimator of the $CDP_{S}$ of a rule $R$ will be
discussed in the next section as it is related to the estimator of the
$CRP_{S}$.

###  4.2 Regional RF and $CRP_{S}$

In this section, we focus on the estimation of the
$CRP_{S}(\mathscr{Y}^{\star},\boldsymbol{x},C_{S}(\boldsymbol{x};\mathscr{Y}^{\star}))=P(f(\boldsymbol{X})\in\mathscr{Y}^{\star}\;|\boldsymbol{X}_{S}\in
C_{S}(\boldsymbol{x};\mathscr{Y}^{\star}),\boldsymbol{X}_{\bar{S}}=\boldsymbol{x}_{\bar{S}})$
and
$CRP_{S}(\mathscr{Y}^{\star},\boldsymbol{R},C_{S}(\boldsymbol{R};\mathscr{Y}^{\star}))=P(f(\boldsymbol{X})\in\mathscr{Y}^{\star}\;|\boldsymbol{X}_{S}\in
C_{S}(\boldsymbol{R};\mathscr{Y}^{\star}),\boldsymbol{X}_{\bar{S}}\in\boldsymbol{R}_{\bar{S}})$.
For simplicity, we remove the dependency of the rectangles in
$\mathscr{Y}^{\star}$. Based on the previous Section, we already know that the
estimators using the RF will be in the form of
$\widehat{CRP}_{S}\left(\mathscr{Y}^{\star},\boldsymbol{x},C_{S}(\boldsymbol{x})\right)=\sum_{i=1}^{n}w_{n,i}(\boldsymbol{x})\mathds{1}_{Y_{i}\in\mathscr{Y}^{\star}}$,
thus we only need to find the right weighting. The main challenge is that we
have a condition based on a region, e.g., $\boldsymbol{X}_{S}\in
C_{S}(\boldsymbol{x})$ or
$\boldsymbol{X}_{\bar{S}}\in\boldsymbol{R}_{\bar{S}}$ (regional-based) instead
of condition of type $\boldsymbol{X}_{S}=\boldsymbol{x}_{S}$ (fixed value-
based) as usually. However, we introduced a natural generalization of the RF
algorithm to make predictions when the conditions are both regional-based and
fixed value-based. Thus, the case where there are only regional-based
conditions are naturally derived.

#### Regional RF to estimate
$CRP_{S}(\mathscr{Y}^{\star},\boldsymbol{x},C_{S}(\boldsymbol{x}))=P(f(\boldsymbol{X})\in\mathscr{Y}^{\star}\;|\boldsymbol{X}_{S}\in
C_{S}(\boldsymbol{x}),\boldsymbol{X}_{\bar{S}}=\boldsymbol{x}_{\bar{S}})$:

The algorithm is based on a slight modification of RF. Its works as follow: we
drop the observations in the initial trees, if a split used variable
$i\in\bar{S}$, i.e., fixed value-based condition, we use the classic rules of
RF, if $x_{i}\leq t$, the observations go to the left children, otherwise the
right children. However, if a split used variable $i\in S$, i.e, regional-
based condition, we use the rectangles
$C_{S}(\boldsymbol{x})=\prod_{i=1}^{|S|}[a_{i},b_{i}]$. The observations are
sent to the left children if $b_{i}\leq t$, right children if $a_{i}>t$ and if
$t\in[a_{i},b_{i}]$ the observations are sent both to the left and right
children. Therefore, we use the weights of the Regional RF algorithm to
estimate the $CRP_{S}$ as in equation 4.3, the estimator is
$\widehat{CRP}_{S}(y^{\star};\boldsymbol{x},C_{S}(\boldsymbol{x}))=\sum_{i=1}^{n}w_{n,i}(\boldsymbol{x})\mathds{1}_{Y_{i}=y^{\star}}$.
A more detailed version of the algorithm is provided and discussed in
Appendix.

To estimate the $CDP$ of a rule
$CDP_{S}\left(\mathscr{Y}^{\star};\boldsymbol{R}\right)=P\left(f(\boldsymbol{X})\in\mathscr{Y}^{\star}\left|\boldsymbol{X}_{\bar{S}}\in\boldsymbol{R}_{\bar{S}}\right.\right)$,
we just have to apply the projected Forest algorithm to the Regional RF, i.e.,
when a split involving a variable outside of $\bar{S}$ is met, the
observations are sent both to the left and right children nodes, otherwise we
use the Regional RF split rule, i.e., if an interval of
$\boldsymbol{R}_{\bar{S}}$ is below $t$, the observations go to the left
children, otherwise the right children and if $t$ is in the interval, the
observations go to the left and right children. The estimator of the
$CRP_{S}(\mathscr{Y}^{\star};\boldsymbol{R},C_{S}(\boldsymbol{R}))$ for rule
is also derived from the Regional RF. Indeed, it is a special case of the
Regional RF algorithm where there are only regional-based conditions.

##  5 Learning the Counterfactual Rules

We compute the Local and Regional CR using the estimators of the previous
section. First, we find the Minimal Divergent Explanation in the same way as
Minimal Sufficient Explanation can be found (Amoukou and Brunel, 2021). As the
exploration of all possible subsets is exponential, we search the Minimal
Divergent Subset among the $K=10$ most frequently selected variables in the RF
$m_{k,n}$ used to estimate the probabilities $CDP_{S},CRP_{S}$ ($K$ is an
hyper-parameter to select according to the use case and computational power).
We can also use any importance measure.

Given an instance $\boldsymbol{x}$ or rectangle $\boldsymbol{R}$ (and set
$\mathscr{Y}^{\star}$) and their corresponding Minimal Divergent Explanation
S, we want to find a rule $C_{S}(\boldsymbol{x})=\prod_{i\in S}[a_{i},b_{i}]$
s.t. given $\boldsymbol{X}_{\bar{S}}=\boldsymbol{x}_{\bar{S}}$ or
$\boldsymbol{X}_{\bar{S}}\in\boldsymbol{R}_{\bar{S}}$ and
$\boldsymbol{X}_{S}\in C_{S}(\boldsymbol{x})$, the probability that
$Y\in\mathscr{Y}^{\star}$ is high. More formally, we want:
$P(f(\boldsymbol{X})\in\mathscr{Y}^{\star}|\boldsymbol{X}_{S}\in
C_{S}(\boldsymbol{x}),\boldsymbol{X}_{\bar{S}}=\boldsymbol{x}_{\bar{S}})$ or
$P(f(\boldsymbol{X})\in\mathscr{Y}^{\star}|\boldsymbol{X}_{S}\in
C_{S}(\boldsymbol{x}),\boldsymbol{X}_{\bar{S}}\in\boldsymbol{R}_{\bar{S}})$
above $\pi_{C}$.

The computation of the rectangles $C_{S}(\boldsymbol{x})=\prod_{i\in
S|}[a_{i},b_{i}]$ relies heavily on our use of RF and on the algorithmic trick
of the projected RF. Indeed, the rectangles defining the rules arise naturally
from RF, while AReS (Rawal and Lakkaraju, 2020) relies on binned variables to
generate candidate rules and tests all these possible rules for choosing an
optimal one. We overcome the computational burden and the challenge of
choosing the number of bins.

Figure 2: The partition of the RF learned to classify the toy data (Green/Blue
stars). Its has 10 leaves. The explainee $\boldsymbol{x}$ is the Blue triangle
in leaf 5.

Figure 3: The partition of the projected Forest when we condition on $X_{0}$,
i.e., ignoring the splits based on $X_{1}$ (the dashed lines).

Figure 4: The optimal CR for $\boldsymbol{x}$ when we condition given
$X_{0}=x_{0}$ is the Green region, its corresponds to the union of leaf 3 and
4 of the forest

To illustrate the idea, we use a two-dimensional data $(X_{0},X_{1})$ with
label Y represented as Green/Blue stars in figure 4. We fit a Random Forest to
classify this dataset and show its partition in figure 4. The explainee
$\boldsymbol{x}$ is the Blue triangle observation.

By looking at the different cells/leaves of the RF, we can guess that the
Minimal Divergent Explanation of $\boldsymbol{x}$ is $S=X_{1}$. Indeed, in
figure 4, we observe the leaves of the Projected Forest when we do not
condition on $S=X_{1}$, thus projected the RF’s partition only on the subspace
$X_{0}$. Its consists of ignoring all the splits in the other directions (here
the $X_{1}$-axis), thus $\boldsymbol{x}$ falls in the projected leaf 2 (see
figure 4) and its $CDP$ is
$CDP_{X_{1}}(\text{Green};\boldsymbol{x})=\frac{10\text{ Green}}{10\text{
Green}+17\text{ Blue}}=0.58$.

Finally, the problem of finding the optimal rectangle
$C_{S}(\boldsymbol{x})=[a_{i},b_{i}]$ in the direction of $X_{1}$ s.t. the
decision changes can be easily solved by using the leaves of the RF. In fact,
by looking at the leaves of the RF (figure 4) of the observations that belong
in the projected RF leaf 2 (figure 4) where $\boldsymbol{x}$ falls, we see in
figure 4 that the optimal rectangle to change the decision given $X_{0}=x_{0}$
or being in the projected RF leaf 2 is the union of the intervals on $X_{1}$
of the leaf 3 and 4 of the RF (see the Green region of figure 4).

Given an instance $\boldsymbol{x}$ and its Minimal Divergent Explanation $S$,
the first step is the collect of the observations which belong to the leaf of
the Projected Forest given $\bar{S}$ where $\boldsymbol{x}$ falls. It
corresponds to the observations that has positive weights in the computation
of the
$CDP_{S}(\mathscr{Y}^{\star};\boldsymbol{x})=\sum_{i=1}^{n}w_{n,i}(\boldsymbol{x}_{\bar{S}})\mathds{1}_{Y_{i}\in\mathscr{Y}^{\star}}$,
i.e., $\\{\boldsymbol{x}_{i}:w_{n,i}(\boldsymbol{x}_{\bar{S}})>0\\}$. Then, we
used the partition of the original forest to find the possible leaves
$C_{S}(\boldsymbol{x})$ in the direction $S$. The possible leaves is among the
RF’s leaves of the collected observations
$\\{\boldsymbol{x}_{i}:w_{n,i}(\boldsymbol{x}_{\bar{s}})>0\\}$. Let denote
$L(\boldsymbol{x}_{i})$ the leaves of the observations $\boldsymbol{x}_{i}$
with $w_{n,i}(\boldsymbol{x}_{\bar{S}})>0$. A possible leaf is a leaf
$L(\boldsymbol{x}_{i})$ s.t.
$CRP_{S}(\mathscr{Y}^{\star},\boldsymbol{x},L(\boldsymbol{x}_{i})_{S})=P(f(\boldsymbol{X})\in\mathscr{Y}^{\star}|\boldsymbol{X}_{S}\in
L(\boldsymbol{x}_{i})_{S},\boldsymbol{X}_{\bar{S}}=\boldsymbol{x}_{\bar{S}})\geq\pi_{C}$.
Finally, we merge all the neighboring possible leaves to get the largest
rectangle, and this maximal rectangle is the counterfactual rule. Note that
the union of the possible leaves is not necessary a connected space, thus we
can have multiple counterfactual rules.

We apply the same idea to find the regional CR. Given a rule $\boldsymbol{R}$
and its Minimal Divergent Explanation $S$, we used the Projection given
$\boldsymbol{X}_{\bar{S}}\in\boldsymbol{R}_{\bar{S}}$ to find the compatible
observations and their leaves and combine the possible ones to obtain the
regional CR that has
$CRP_{S}(\mathscr{Y}^{\star},\boldsymbol{R},C_{S}(\boldsymbol{R}))\geq\pi_{C}$.
For example, if we consider the leaf 5 of the original forest as a rule: If
$\boldsymbol{X}\in$ Leaf 5, then predict Blue. Its Minimal Divergent
Explanation is also $S=X_{1}$. The R-CR would also be the Green region in
figure 4. Indeed, if we satisfy the $X_{0}$ condition of the leaf 5 and
$X_{1}$ condition of the leaf 3 and 4, then the decision change to Green.

##  6 Experiments

To demonstrate the performance of our framework, we conduct two experiments on
real-world datasets. The first consists of showing how we can use the Local
Counterfactual Rules for explaining a regression model. In the second
experiment, we compare our approaches with the 2 baselines methods in
classification problem: (1) CET (Kanamori et al., 2022), which partition the
input space using a decision tree and associate a vector perturbation for each
leaf, (2) AReS (Rawal and Lakkaraju, 2020) performs an exhaustive search for
finding global counterfactual rules, but we used the implementation of
Kanamori et al. (2022) that adapts the algorithm for returning counterfactuals
samples instead of rules. We compare the methods only in classification
problem as most prior works do not deal regression problem. In all
experiments, we split our dataset into train ($75\%$) - test ($25\%$), and we
learn a model $f$, a LightGBM (estimators=50, nb leaves=8), on the train set
that is the explainee. We learn $f$’s predictions on the train set with an
approximating RF $m_{nb,n}$ (estimators=20, max depth=10): that will be used
to generate the CR with $\pi=0.9$. The used parameters for AReS, CET are max
rules=8, bins=10 and max iterations=1000, max leaf=8, bins=10 respectively.
Due to page limitation, the detailed parameters of each method are provided in
Appendix.

#### Sampling CE using the Counterfactual Rules:

Notice that our approaches cannot be directly compare with the baseline
methods since they all return counterfactual samples while we give rules
(range of vector values) that permit to change the decision with high
probability. However, we adapt the CR to generate also counterfactual samples
using a generative model. For example, given an instance
$\boldsymbol{x}=(\boldsymbol{x}_{S},\boldsymbol{x}_{\bar{S}})$, target
$\mathscr{Y}^{\star}$ and its counterfactual rule
$C_{S}(\boldsymbol{x};\mathscr{Y}^{\star})$, we want to find a sample
$x^{\star}=(\boldsymbol{z}_{S},\boldsymbol{x}_{\bar{S}})$ with
$\boldsymbol{z}_{S}\in C_{S}(\boldsymbol{x},\mathscr{Y}^{\star})$ s.t
$\boldsymbol{x}^{\star}$ is an in-distribution sample and
$f(\boldsymbol{x}^{\star})\in\mathscr{Y}^{\star}$. Instead of using a complex
conditional generative model as (Xu et al., 2019; Patki et al., 2016) that can
be difficult to calibrate, we use an energy-based generative approach
(Grathwohl et al., 2020; Lecun et al., 2006). The core idea is to find
$\boldsymbol{z}_{S}\in C_{S}(\boldsymbol{x},y^{\star})$ s.t.
$\boldsymbol{x}^{\star}$ maximize a given energy score to ensure that it is an
in-distribution sample. As an example of an energy function, we use the
negative outlier score of an Isolation Forest (Liu et al., 2008). We use
Simulated Annealing (see (Guilmeau et al., 2021) for a review) to maximize the
negative outlier score using the information of the counterfactual rules
$C_{S}(\boldsymbol{x};\mathscr{Y}^{\star})$. In fact, the range values given
by the CR $C_{S}(\boldsymbol{x};\mathscr{Y}^{\star})$ reduce the search space
for $\boldsymbol{z}_{S}$ drastically. We used the training set
$\mathcal{D}_{n}$ to find the possible values i.e., we defined $P_{i}$,
$P_{S}$ as the list of values of the variable $i\in S$ found in
$\mathcal{D}_{n}$ and
$P_{S}=\\{\boldsymbol{z}_{S}=(z_{1},\dots,z_{S}):\boldsymbol{z}_{S}\in
C_{S}(\boldsymbol{x},y^{\star}),z_{i}\in P_{i}\\}$ the possible values of
$\boldsymbol{z}_{S}$ respectively. Then, we sample $\boldsymbol{z}_{S}$ in the
set $P_{S}$ and use Simulated Annealing to find a $\boldsymbol{x}^{\star}$
that maximizes the negative outlier score. Note that the algorithm is the same
for sampling CE with the Regional-CR. A more detailed version of the algorithm
is provided in Appendix.

Finally, we compare the methods on unseen observations using three criteria.
Correctness is the average number of instances for which acting as prescribed
change to the desired prediction. Plausibility is the average number of inlier
(predict by an Isolation Forest) in the counterfactual samples. Sparsity is
the average number of features that have been changed, and especially for the
global counterfactual methods (AReS, Regional-CR) that do not ensure to cover
all the instances, we compute Coverage that corresponds to the average number
of unseen observations we cover.

#### Local counterfactual rules for regression:

We give recourse for the California House Price dataset (Kelley Pace and
Barry, 1997) derived from the 1990 U.S. census. We have information about each
district (demography, …), and the goal is to predict the median house value of
each district.

To illustrate the efficiency of the Local-CR, we select all the observations
in the test set having a price lower than $100k$ (1566 houses), and we aim to
find the recourse that permit to increase their price : we want the price $y$
to be in the interval $\mathscr{Y}^{\star}=[200k,250k]$. For each instance
$\boldsymbol{x}$, we compute the Minimal Divergent Explanation $S$, the Local-
CR $C_{S}(\boldsymbol{x};[200k,250k])$ and a CE using the Simulated Annealing
as described above. We succeed in changing the decision of all the
observations $(\textit{Correctness}=1)$ and most of them passed the outlier
test with $\textit{Plausibility}=0.92$. On top of that, our Local-CR have
sparse support ($\textit{Sparsity}=4.45$). For example, the Local-CR of the
instance $\boldsymbol{x}=$ (Longitude=-118.2, latitude=33.8, housing median
age=26, total rooms=703, total bedrooms=202, population=757, households=212,
median income=2.52) is $C_{S}(\boldsymbol{x},[200k,250k])=$ (total room
$\in[2132,3546],$ total bedrooms $\in[214,491]$). It means if total room and
total bedrooms satisfy the conditions in $C_{S}(\boldsymbol{x},[200k,250k])$
and the remaining features of $\boldsymbol{x}$ is fixed, then the probability
that the price is in $[200k,250k]$ is 0.97.

#### Comparisons of Local-CR and Regional-CR with baselines (AReS, CET):

We use 3 real-world datasets: Diabetes (Kaggle, 2016) contains diagnostic
measurements and aims to predict whether or not a patient has diabetes, Breast
Cancer Wisconsin (BCW) (Dua and Graff, 2017) consists of predicting if a tumor
is benign or not using the characteristic of the cell nuclei, and Compas
(Larson et al., 2016) was used to predict recidivism, and it contains
information about the criminal history, demographic attributes. During the
evaluation, we observe that AReS, CET are very sensitive to the number of bins
and the maximal number of rules or actions as noticed by (Ley et al., 2022). A
bad parameterization gives completely useless explanations. Moreover, a
different model needs to be trained for each class to be accurate, while we
only need to have a RF that has good precision.

In table 1, we notice that the Local and Regional-CR succeed in changing
decisions with a high accuracy in all datasets, outperforming AReS and CET
with a large margin on BCW, and Diabetes. Moreover, we notice that the
baselines struggle to change at the same time the positive and negative class,
(e.g. CET has Acc=1 in the positive class, and 0.21 for the negative class on
BCW) or when they have a good Acc, the CE are not plausible. For instance, CET
has Acc=0.98 and Psb=0 on Compas, meaning that all the CE are outlier.
Regarding the coverage of the global CE, CET covers all the instances as it
partitions the space, but we observe that AReS has a smaller
Coverage$=\\{0.43,0.44,0.81\\}$ than the Regional-CR which has $\\{1,0.7,1\\}$
for BCW, Diabetes, and Compas respectively. To sum up, the CR is easier to
train and provides more accurate and plausible rules than the baselines
methods.

Table 1: Results of the Correctness (Acc), Plausibility, and Sparsity (Sprs)
of the different methods. We compute each metric according to the positive
(Pos) and negative (Neg) class.

| COMPAS | BCW | Diabetes
---|---|---|---
| Acc | Psb | Sps | Acc | Psb | Sps | Acc | Psb | Sps
| Pos | Neg | Pos | Neg | Pos | Neg | Pos | Neg | Pos | Neg | Pos | Neg | Pos | Neg | Pos | Neg | Pos | Neg
L-CR | 1 | 0.9 | 0.87 | 0.73 | 2 | 4 | 1 | 1 | 0.96 | 1 | 9 | 7 | 0.97 | 1 | 0.99 | 0.8 | 3 | 4
R-CR | 0.9 | 0.98 | 0.74 | 0.93 | 2 | 3 | 0.89 | 0.9 | 0.94 | 0.93 | 9 | 9 | 0.99 | 0.99 | 0.9 | 0.87 | 3 | 4
AReS | 0.98 | 1 | 0.8 | 0.61 | 1 | 1 | 0.63 | 0.34 | 0.83 | 0.80 | 4 | 3 | 0.73 | 0.60 | 0.77 | 0.86 | 1 | 1
CET | 0.85 | 0.98 | 0.7 | 0 | 2 | 2 | 1 | 0.21 | 0.6 | 0.80 | 8 | 2 | 0.84 | 1 | 0.60 | 0.20 | 6 | 6

##  7 Conclusion

Most current works that generate CE are implicit through an optimization
process or a brunch of random samples, thus lacking guarantees. For this
reason, we rethink CE as Counterfactual Rules. For any individual or sub-
population, it gives the simplest policies that change the decision with high
probability. Our approach learns robust, plausible, and sparse adversarial
regions where the observations should be moved. We make central use of Random
Forests, which give consistent estimates of the interest probabilities and
naturally give the counterfactual rules we want to extract. In addition, it
permits us to deal with regression problems and continuous features.
Consequently, our methods are suitable for all datasets where tree-based model
performs well (e.g., tabular data). A prospective work is to evaluate the
robustness of our methods to noisy human responses, i.e., when the prescribed
recourse is not implemented exactly, and to refine the methodology for
selecting the threshold probabilities $\pi$ and $\pi_{C}$.

## References

  * Amoukou and Brunel (2021) Salim I Amoukou and Nicolas JB Brunel.  Consistent sufficient explanations and minimal local rules for explaining regression and classification models.  _arXiv preprint arXiv:2111.04658_ , 2021.
  * Bénard et al. (2021a) Clément Bénard, Gérard Biau, Sébastien Da Veiga, and Erwan Scornet.  Shaff: Fast and consistent shapley effect estimates via random forests.  _arXiv preprint arXiv:2105.11724_ , 2021a.
  * Bénard et al. (2021b) Clément Bénard, Gérard Biau, Sébastien Veiga, and Erwan Scornet.  Interpretable random forests via rule extraction.  In _International Conference on Artificial Intelligence and Statistics_ , pages 937–945. PMLR, 2021b.
  * Bénard et al. (2021c) Clément Bénard, Sébastien Da Veiga, and Erwan Scornet.  Mda for random forests: inconsistency, and a practical solution via the sobol-mda.  _arXiv preprint arXiv:2102.13347_ , 2021c.
  * Breiman et al. (1984) Leo Breiman, Jerome Friedman, Richard Olshen, and Charles Stone.  Classification and regression trees. wadsworth int.  _Group_ , 37(15):237–251, 1984.
  * Chen et al. (2012) S. Chen, Arthur Choi, and Adnan Darwiche.  The same-decision probability: A new tool for decision making.  2012\.
  * Chou et al. (2022) Yu-Liang Chou, Catarina Moreira, Peter Bruza, Chun Ouyang, and Joaquim Jorge.  Counterfactuals and causability in explainable artificial intelligence: Theory, algorithms, and applications.  _Information Fusion_ , 81:59–83, 2022.  ISSN 1566-2535.  doi: https://doi.org/10.1016/j.inffus.2021.11.003.  URL https://www.sciencedirect.com/science/article/pii/S1566253521002281.
  * Du et al. (2021) Qiming Du, Gérard Biau, François Petit, and Raphaël Porcher.  Wasserstein random forests and applications in heterogeneous treatment effects.  In _International Conference on Artificial Intelligence and Statistics_ , pages 1729–1737. PMLR, 2021.
  * Dua and Graff (2017) Dheeru Dua and Casey Graff.  UCI machine learning repository, 2017.  URL http://archive.ics.uci.edu/ml.
  * Grathwohl et al. (2020) Will Grathwohl, Kuan-Chieh Wang, Joern-Henrik Jacobsen, David Duvenaud, Mohammad Norouzi, and Kevin Swersky.  Your classifier is secretly an energy based model and you should treat it like one.  In _International Conference on Learning Representations_ , 2020.
  * Guilmeau et al. (2021) Thomas Guilmeau, Emilie Chouzenoux, and Víctor Elvira.  Simulated annealing: a review and a new scheme.  pages 101–105, 07 2021.  doi: 10.1109/SSP49050.2021.9513782.
  * Ishwaran et al. (2008) Hemant Ishwaran, Udaya B Kogalur, Eugene H Blackstone, and Michael S Lauer.  Random survival forests.  _The annals of applied statistics_ , 2(3):841–860, 2008.
  * Kaggle (2016) Kaggle.  Pima indians diabetes database, 2016.  URL https://www.kaggle.com/datasets/uciml/pima-indians-diabetes-database.
  * Kanamori et al. (2020) Kentaro Kanamori, Takuya Takagi, Ken Kobayashi, and Hiroki Arimura.  Dace: Distribution-aware counterfactual explanation by mixed-integer linear optimization.  In _IJCAI_ , 2020.
  * Kanamori et al. (2022) Kentaro Kanamori, Takuya Takagi, Ken Kobayashi, and Yuichi Ike.  Counterfactual explanation trees: Transparent and consistent actionable recourse with decision trees.  In _Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, PMLR 151:1846-1870_ , 2022.
  * Karimi et al. (2020a) Amir-Hossein Karimi, Gilles Barthe, Borja Balle, and Isabel Valera.  Model-agnostic counterfactual explanations for consequential decisions.  _ArXiv_ , abs/1905.11190, 2020a.
  * Karimi et al. (2020b) Amir-Hossein Karimi, Gilles Barthe, Bernhard Schölkopf, and Isabel Valera.  A survey of algorithmic recourse: definitions, formulations, solutions, and prospects.  _CoRR_ , abs/2010.04050, 2020b.  URL https://arxiv.org/abs/2010.04050.
  * Kelley Pace and Barry (1997) R. Kelley Pace and Ronald Barry.  Sparse spatial autoregressions.  _Statistics, Probability Letters_ , 33(3):291–297, 1997.  ISSN 0167-7152.  doi: https://doi.org/10.1016/S0167-7152(96)00140-X.  URL https://www.sciencedirect.com/science/article/pii/S016771529600140X.
  * Lakkaraju et al. (2022) Himabindu Lakkaraju, Dylan Slack, Yuxin Chen, Chenhao Tan, and Sameer Singh.  Rethinking explainability as a dialogue: A practitioner’s perspective.  _CoRR_ , abs/2202.01875, 2022.  URL https://arxiv.org/abs/2202.01875.
  * Larson et al. (2016) Jeff Larson, Surya Mattu, Lauren Kirchner, , and Julia Angwin.  How we analyzed the compas recidivism algorithm, 2016.  URL https://www.propublica.org/article/how-we-analyzed-the-compas-recidivism-algorithm.
  * Lecun et al. (2006) Yann Lecun, Sumit Chopra, and Raia Hadsell.  _A tutorial on energy-based learning_.  01 2006.
  * Ley et al. (2022) Dan Ley, Saumitra Mishra, and Daniele Magazzeni.  Global counterfactual explanations: Investigations, implementations and improvements, 2022.  URL https://arxiv.org/abs/2204.06917.
  * Lin et al. (2020) Jimmy J. Lin, Chudi Zhong, Diane Hu, Cynthia Rudin, and Margo I. Seltzer.  Generalized optimal sparse decision trees.  _ArXiv_ , abs/2006.08690, 2020.
  * Liu et al. (2008) Fei Tony Liu, Kai Ming Ting, and Zhi-Hua Zhou.  Isolation forest.  In _2008 eighth ieee international conference on data mining_ , pages 413–422. IEEE, 2008.
  * Loh (2011) Wei-Yin Loh.  Classification and regression trees.  _Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery_ , 1, 2011.
  * Looveren and Klaise (2019) Arnaud Van Looveren and Janis Klaise.  Interpretable counterfactual explanations guided by prototypes.  _CoRR_ , abs/1907.02584, 2019.  URL http://arxiv.org/abs/1907.02584.
  * Lundberg et al. (2020) Scott M. Lundberg, Gabriel Erion, Hugh Chen, Alex DeGrave, Jordan M. Prutkin, Bala Nair, Ronit Katz, Jonathan Himmelfarb, Nisha Bansal, and Su-In Lee.  From local explanations to global understanding with explainable ai for trees.  _Nature Machine Intelligence_ , 2(1):2522–5839, 2020.
  * Meinshausen and Ridgeway (2006) Nicolai Meinshausen and Greg Ridgeway.  Quantile regression forests.  _Journal of Machine Learning Research_ , 7(6), 2006.
  * Molnar (2022) Christoph Molnar.  _Interpretable Machine Learning_.  2 edition, 2022.  URL https://christophm.github.io/interpretable-ml-book.
  * Mothilal et al. (2020) Ramaravind K. Mothilal, Amit Sharma, and Chenhao Tan.  Explaining machine learning classifiers through diverse counterfactual explanations.  In _Proceedings of the 2020 Conference on Fairness, Accountability, and Transparency_ , FAT* ’20, page 607–617, New York, NY, USA, 2020. Association for Computing Machinery.  ISBN 9781450369367.  doi: 10.1145/3351095.3372850.  URL https://doi.org/10.1145/3351095.3372850.
  * Parmentier and Vidal (2021) Axel Parmentier and Thibaut Vidal.  Optimal counterfactual explanations in tree ensembles.  _CoRR_ , abs/2106.06631, 2021.  URL https://arxiv.org/abs/2106.06631.
  * Patki et al. (2016) N. Patki, R. Wedge, and K. Veeramachaneni.  The synthetic data vault.  In _2016 IEEE International Conference on Data Science and Advanced Analytics (DSAA)_ , pages 399–410, Oct 2016.  doi: 10.1109/DSAA.2016.49.
  * Pawelczyk et al. (2022) Martin Pawelczyk, Teresa Datta, Johannes van-den Heuvel, Gjergji Kasneci, and Himabindu Lakkaraju.  Algorithmic recourse in the face of noisy human responses, 2022.  URL https://arxiv.org/abs/2203.06768.
  * Poyiadzi et al. (2019) Rafael Poyiadzi, Kacper Sokol, Raúl Santos-Rodriguez, Tijl De Bie, and Peter A. Flach.  FACE: feasible and actionable counterfactual explanations.  _CoRR_ , abs/1909.09369, 2019.  URL http://arxiv.org/abs/1909.09369.
  * Rawal and Lakkaraju (2020) Kaivalya Rawal and Himabindu Lakkaraju.  Beyond individualized recourse: Interpretable and interactive summaries of actionable recourses.  _Advances in Neural Information Processing Systems_ , 33:12187–12198, 2020.
  * Ribeiro et al. (2016) Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin.  " why should i trust you?" explaining the predictions of any classifier.  In _Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining_ , pages 1135–1144, 2016.
  * Russell (2019) Chris Russell.  Efficient search for diverse coherent explanations.  In _Proceedings of the Conference on Fairness, Accountability, and Transparency_ , FAT* ’19, page 20–28, New York, NY, USA, 2019. Association for Computing Machinery.  ISBN 9781450361255.  doi: 10.1145/3287560.3287569.  URL https://doi.org/10.1145/3287560.3287569.
  * Ustun et al. (2019) Berk Ustun, Alexander Spangher, and Yang Liu.  Actionable recourse in linear classification.  _Proceedings of the Conference on Fairness, Accountability, and Transparency_ , 2019.
  * Verma et al. (2020) Sahil Verma, John P. Dickerson, and Keegan Hines.  Counterfactual explanations for machine learning: A review.  _CoRR_ , abs/2010.10596, 2020.  URL https://arxiv.org/abs/2010.10596.
  * Wachter et al. (2017a) Sandra Wachter, Brent Mittelstadt, and Chris Russell.  Counterfactual explanations without opening the black box: Automated decisions and the gdpr.  _Harv. JL & Tech._, 31:841, 2017a.
  * Wachter et al. (2017b) Sandra Wachter, Brent Daniel Mittelstadt, and Chris Russell.  Counterfactual explanations without opening the black box: Automated decisions and the gdpr.  _Cybersecurity_ , 2017b.
  * Wager and Athey (2017) Stefan Wager and Susan Athey.  Estimation and inference of heterogeneous treatment effects using random forests, 2017.
  * Wang et al. (2017) Tong Wang, Cynthia Rudin, Finale Doshi-Velez, Yimin Liu, Erica Klampfl, and Perry MacNeille.  A bayesian framework for learning rule sets for interpretable classification.  _J. Mach. Learn. Res._ , 18:70:1–70:37, 2017.
  * Xu et al. (2019) Lei Xu, Maria Skoularidou, Alfredo Cuesta-Infante, and Kalyan Veeramachaneni.  Modeling tabular data using conditional gan.  In _NeurIPS_ , 2019.


#  On fractional-order maps and their synchronization

Prashant M. Gade1,∗ and Sachin B. Bhalekar2,3
1 Department of Physics,
Rashtrasant Tukadoji Maharaj Nagpur University, Nagpur, India
Email<EMAIL_ADDRESS>(*Corresponding author)
2Department of Mathematics,
Shivaji University, Kolhapur, India
3School of Mathematics and Statistics,
University of Hyderabad, Hyderabad, India
Email<EMAIL_ADDRESS>
###### Abstract

We study the stability of linear fractional order maps. We show that in the
stable region, the evolution is described by Mittag-Leffler functions and a
well defined effective Lyapunov exponent can be obtained in these cases. For
one-dimensional systems, this exponent can be related to the corresponding
fractional differential equation. A fractional equivalent of map $f(x)=ax$ is
stable for $a_{c}(\alpha)<a<1$ where $\alpha$ is a fractional order parameter
and $a_{c}(\alpha)\approx-\alpha$. For coupled linear fractional maps, we can
obtain ‘normal modes’ and reduce the evolution to effectively one-dimensional
system. If the eigenvalues are real the stability of the coupled system is
dictated by the stability of effectively one-dimensional normal modes. For
complex eigenvalues, we obtain a much richer picture. However, in the stable
region, the evolution of modulus is dictated by Mittag-Leffler function and
the effective Lyapunov exponent is determined by modulus of eigenvalues. We
extend these studies to synchronized fixed points of fractional nonlinear
maps.

##  1 Introduction

After May’s influential paper[1], the difference equations gained wide
popularity in all walks of science. They were studied from the viewpoint of
nonlinear dynamics and chaos which was ubiquitous in disciplines ranging from
ecology to economics. Mathematicians studied it as an interesting object in
its right[2]. The studies on maps are computationally less intensive. In
certain cases, it was easier to track them analytically. Feigenbaum’s study of
period-doubling bifurcations in logistic maps is an example[3, 4]. Lower order
difference equations can show the same phenomena as higher-order differential
equations. For example, we can observe chaos in a logistic map that has a
single variable while we need at least three first-order autonomous
differential equations to observe chaos. Insights gained from studies in
difference equations were often (though not always) useful in studies in
differential equations and vice versa. From control problems to
synchronization schemes, many theoretical ideas have found applications in
both difference and differential equations[5]. In this work, we study
fractional linear difference equations and coupled fractional linear
difference equations. We show that the analysis of linear fractional
difference equations has remarkable similarities with corresponding fractional
differential equations. For coupled difference equations, the Jacobian, its
eigenvalues and eigenvectors play a central role. We show that similar
concepts can be very useful for coupled fractional difference equations.

Synchronization of dynamical systems has been extensively studied in the past
three decades both theoretically as well as experimentally. Exact
synchronization for master-slave type systems as well as synchronization for
mutually coupled systems has been investigated extensively. Apart from exact
synchronization, several other types such as anti-synchronization, anticipated
synchronization, lag synchronization, phase synchronization, generalized
synchronization, etc. have been studied. For spatially extended systems, the
connectivity matrix of underlying topology plays an important role in
synchronization. For exact synchronization, we have a clear mathematical
formalism for finding necessary and sufficient conditions for synchronization.

Ordinary differential equations have successfully described a variety of
physical system and found countless applications since Newton and Leibniz
introduced them. However, memory plays an important role in many physical
systems. Apart from a mathematical curiosity, such systems are modeled by
fractional-order differential equations from viewpoint of applications. For
fractional-order systems, fractional order maps have been introduced recently.
These systems are not very well investigated yet. In a previous work[6], a
coupled map lattice model of fractional order maps was studied. It was found
that it is possible to have synchronization even in the thermodynamic limit in
these systems. However, the error reduced as a power-law. The exponent is the
same as the fractional order of maps. We give certain pointers to understand
these findings analytically in this work.

As in differential equations, we start by studying the stability of a fixed
point in linear systems. We study two coupled fractional maps and derive
conditions for synchronization analytically. For simple linear maps the bounds
of synchronization and its relation to Mittag-Leffler function can be shown
analytically. Asymptotically, Mittag-Leffler function behaves as a power-law
and this could be the reason for the observation in previous work. For a
linear function with constant slope, this relation can be shown very clearly.
We will be studying symmetric coupling. However, most of the studies can be
easily extended to asymmetric coupling. We give conditions for the stability
of the synchronized state and demonstrate it with certain examples. We propose
that the usual definition of Lyapunov exponent using logarithm (which is
inverse of exponential) should be appropriately modified to obtain an accurate
quantifier that describes the convergence of trajectories in a stable regime.
(Mittag-Leffler function is a power-law asymptotically which is slower decay
than exponential. Thus the Lyapunov exponent will always be zero if we fit it
with an exponential.) However, the divergence of trajectories in the unstable
region is still exponential and the usual definition of Lyapunov exponent may
hold in this case.

There are several definitions of fractional differential equations and the
same is true for fractional difference equations. We will study the definition
obtained by Gejji and Deshpande [7]. The evolution depends on the value at all
previous time-steps. The weight of previous values decays as a power-law and
there is a long term memory built-in in the system. It is not surprising that
the fluctuations also decay extremely slowly and decay can be approximated by
power-law with power related to the order of fractional difference equation.
This is in turn related to the properties of Mittag-Leffler function which
plays an important role in fractional differential equations and plays a
similar role here. In difference equations, the concept of coupled maps was
introduced by Kapral, Kuznetsov, and Kaneko. Kaneko can be credited for making
it popular.

The fractional difference equations have been studied only recently. In 1989,
Miller and Ross began this investigations[8]. Some of the studies in
fractional difference equations are due to Atici and coworkers [9, 10],
Holm[11], and others [12, 13]. We extend this definition to the coupled
fractional difference equations. While exploratory works on fractional
equivalents of known nonlinear maps can yield useful insights, we focus on
linear systems in this work. Linearizations of nonlinear systems are a
standard tool in difference equtions as well as coupled difference equations.
Most of analytic work in these systems is dependent on linearization. Thus
understanding linear systems and coupled linear systems is extremely important
in nonlinear dynamics. We believe that studies in linear systems can be
equally useful in fractional difference equations.

##  2 Effective Lyapunov Exponent

In this section, we introduce an alternative viz. Effective Lyapunov Exponent
(ELE) for the classical Lyapunov exponent. We will follow the notation and
definition used by Deshpande and Gejji [7]. They define the fractional
equivalent for the $x(n+1)=f(x(n))$ in the following manner. They construct an
discrete Caputo-type fractional difference operator and define
$u(t)=u_{0}+{\frac{1}{\Gamma(\alpha)}}\sum_{j=1}^{t}{\frac{\Gamma(t-j+\alpha)}{\Gamma(t-j+1)}}\left(f(j-1,u(j-1))-u(j-1)\right)$
in general. We assume that the function $f$ does not depend on time. We define
$g_{\alpha}(k)={\frac{\Gamma(k+\alpha)}{\Gamma(k+1)}}$ and alternatively write
above expression as
$u(t)=u_{0}+{\frac{1}{\Gamma(\alpha)}}\sum_{j=1}^{t}g_{\alpha}(t-j)\left(f(u(j-1))-u(j-1)\right)$.

Few one-dimensional maps such as the Gaussian map and Bernoulli map have been
studied in this work. Liu has numerically investigated coupled fractional
Henon map[14]. Henon map is a two-dimensional map and Liu introduces memory in
only one of the variables. We will call such systems ‘fractional difference
equations of inhomogeneous order’. On the other hand, we will investigate maps
of homogeneous order.

Consider $f(x)=rx$ where $r\in{\cal R}$

$u(t)=u_{0}+{\frac{1}{\Gamma(\alpha)}}\sum_{j=1}^{t}g_{\alpha}(t-j)\left(r-1\right)u(j-1)$
(1)

can be identified with the continuous-time system

$D^{\alpha}u(t)=(r-1)u(t)$ (2)

for sufficiently small values of coefficient.

The continuous-time system (2) has exact solution in terms of Mittag-Leffler
function as below

$u(t)=u(0)E_{\alpha}\left((r-1)t^{\alpha}\right).$ (3)

Fig. 1 shows the Mittag-Leffler decay for various values of $\alpha$.

Figure 1: Mittag-Leffler decay.

This can be alternatively written as

$u(t)=u(0)E_{\alpha}\left(\lambda_{e}t^{\alpha}\right).$ (4)

where $\lambda_{e}$ is effective Lyapunov exponent. For a linear first order
ordinary differential equation $x^{\prime}(t)=\lambda x(t)$, the solution
would be $x(t)=x(0)\exp(\lambda t)$ where $\lambda$ is Lyapunov exponent. When
it is negative, the system goes to absorbing state. Our formulation could be
considered as generalization

$\lambda_{e}=\lim_{t\longrightarrow\infty}t^{-\alpha}E_{\alpha}^{-1}\left({\frac{u(t)}{u(0)}}\right).$
(5)

This formulation is very similar to standard definition of Lyapunov exponent
for $\alpha=1$ where $E_{\alpha}(x)=\exp(x)$. We will demonstrate that this
quantity is a well-defined quantity which indeed converges in the stable
regime. On the other hand, if we insist on using the definition of Lyapunov
exponent used for ordinary differential equations, it will lead to zero value.
The reason is that Mittag-Leffler function is a power-law asymptotically which
is slower than exponential. Like Lyapunov exponent, the effective Lyapunov
exponent is negative in the absorbing state. In Fig. 2, we sketch the
numerically computed $\lambda_{e}$ for various values of $r$. It is clear that
$\lambda_{e}=r-1$.

The system is unstable for $r>1$ for any $\alpha$. However, lower bound
$a_{c}(\alpha)$ depends on $\alpha$. Lower bound $a_{c}(\alpha)\rightarrow-1$
as $\alpha\rightarrow 1$. This is an expected limit for integer order
difference equation.

Figure 2: Computation of $\lambda_{e}$ for few values of $r$. The value
obtained is $r-1$.

##  3 Fractional order coupled maps

We define two coupled maps in this setting and define.

$\displaystyle x(t)$ $\displaystyle=$ $\displaystyle
x_{0}+{\frac{1}{\Gamma(\alpha)}}\sum_{j=1}^{t}g_{\alpha}(t-j)G(x(j-1),y(t-j)),$
$\displaystyle y(t)$ $\displaystyle=$ $\displaystyle
y_{0}+{\frac{1}{\Gamma(\alpha)}}\sum_{j=1}^{t}g_{\alpha}(t-j)G(y(j-1),x(t-j))$
(6)

where $G(a,b)=\delta f(a)+qf(b)-a$.

Case 1: Real ‘Normal Modes’

First we consider the case $f(x)=x$ for which the coefficient matrix of the
system has real eigenvalues viz. $\delta+q$ and $\delta-q$. Now we consider
two new variables $u(t)=x(t)+y(t)$ and $v(t)=x(t)-y(t)$ and obtain

$\displaystyle u(t)$ $\displaystyle=$ $\displaystyle
u_{0}+{\frac{1}{\Gamma(\alpha)}}\sum_{j=1}^{t}g_{\alpha}(t-j)((\delta+q-1)u(j-1)),$
(7) $\displaystyle v(t)$ $\displaystyle=$ $\displaystyle
v_{0}+{\frac{1}{\Gamma(\alpha)}}\sum_{j=1}^{t}g_{\alpha}(t-j)((\delta-q-1)v(j-1)).$
(8)

This is a simple decoupled system of linear difference equations.

Again, the discrete-time equations above can be identified with the
continuous-time system

$D^{\alpha}u(t)=(\delta+q-1)u(t)$ (9)

$D^{\alpha}v(t)=(\delta-q-1)v(t)$ (10)

for sufficiently small values of coefficient and exact solution in terms of
Mittag-Leffler function is given as

$u(t)=u(0)E_{\alpha}\left((\delta+q-1)t^{\alpha}\right).$ (11)

$v(t)=v(0)E_{\alpha}\left((\delta-q-1)t^{\alpha}\right).$ (12)

Thus the effective Lyapunov exponents are given by $\delta+q-1$ and
$\delta-q-1$ for the above decoupled system. We compute the system for
$T=8\times 10^{4}-10^{5}$ time-steps. The decay is extremely slow for smaller
$\alpha$. After discarding first 2000 time-steps, we check if the distance of
$(x(t),y(t))$ from origin every 100 time-steps and the trajectory is stable if
the distance does not increase till time $T$.

Figure 3: Stability regime for symmetrically coupled maps.

The stability regime for the above system is plotted as a function of $\delta$
and $q$ in Fig. 3. The stability region is a rhombus bounded by lines parallel
to $\delta=q$ and $\delta=-q$. For all values of $\alpha$, the system is
unstable for $\delta+q-1>0$ and $\delta-q-1>0$. The other two bounding lines
change with $\alpha$. For $q=0$, we have an effectively one-dimensional
system. The bounds for $q=0$ are $\delta=1$ and $\delta=a_{c}(\alpha)$. As
$\alpha\rightarrow 1$, $a_{c}(\alpha)\rightarrow-1$. The lines enclosing
stability region are given by $\delta\pm q=1$ and $\delta\pm
q=a_{c}(\alpha)$.‘ It is clear from $a_{c}(\alpha)$ are close to $-\alpha$.
This explains the rhombus structure in Fig. 3. (For $f(x)=rx$, the phase
diagram will be similar except that the values of $\delta$ and $q$ will be
scaled to ${\frac{\delta}{r}}$ and ${\frac{q}{r}}$.)

Synchronization and antisynchronization
We consider the linear system (6) with real normal modes. As shown in (7)–(8),
the sum and difference variables $u(t)=x(t)+y(t)$ and $v(t)=x(t)-y(t)$ are
effectively decoupled. We can also say the $u$ corresponds to $(1,1)$ mode and
$v$ corresonds to $(1,-1)$ mode. If the effective Lyapunov exponent
corresponding to variable $u$ is in the unstable region while the one
corresponding to $v$ is in stable region, we will find that $v(t)\rightarrow
0$. Thus $x(t)\rightarrow y(t)$ implying synchronization. On the other hand if
$\lambda_{e}$ corresponding to $v(t)$ is in unstable region and corresponding
to $u(t)$ is in the stable region we will observe that $v(t)\rightarrow 0$.
This phenomenon is termed as antisynchronization. Fig. 4 shows both these
phenomena for different values of parameters. We also find that the decaying
mode decays slower than exponential making it necessary to give a new
definition for Lyapunov exponent. However, growing mode increases
exponentially. The reason may be that approximating fractional difference
equation by fractional differential equation may not be valid for large
values. In coupled map lattices, it is known that the condition for chaotic
synchronization is that all modes except one corresponding to the mean, i.e.
eigenmode $(1,1,\ldots 1)$ should decay[15]. We can obtain normal modes in
fractional system in a similar manner and find conditions for chaotic
synchronization.

Figure 4: Chaotic synchronization and antisynchronization for symmetrically
coupled maps. $\alpha=0.9$,a) $\delta$=0.85, $q=-0.2$ and
b)$\delta=0.85,q=0.2$

Case 2: Complex ‘Normal Modes’

Now let us consider another coupled system

$\displaystyle x(t)$ $\displaystyle=$ $\displaystyle
x_{0}+{\frac{1}{\Gamma(\alpha)}}\sum_{j=1}^{t}g_{\alpha}(t-j)(\delta
x(j-1)+qy(t-j)-x(j-1)),$ (13) $\displaystyle y(t)$ $\displaystyle=$
$\displaystyle
y_{0}+{\frac{1}{\Gamma(\alpha)}}\sum_{j=1}^{t}g_{\alpha}(t-j)(\delta
y(j-1)-qx(t-j)-y(j-1)).$ (14)

The eigenvalues of coefficient matrix of the system (13)–(14) are
$\delta\pm\iota q$.

Figure 5: Stability region for coupled maps with antisymmetric coupling. It is
clear that the stability region approaches unit circle as $\alpha\rightarrow
1$ while it is significantly different for small $\alpha$. The stability
region is shown for different values of $\alpha$ in ascending order in a) and
in descending order in b).

We have plotted stability region in $\delta-q$ space for various values of
$\alpha$ in Fig. 5. This structure is far richer than one obtained in Fig. 3.
As $\alpha\rightarrow 1$, the stability region tends to the unit circle in the
complex plane which is a stability region for integer-order difference maps in
two dimensions.

As $\alpha$ decreases, the stability region gradually deforms from an unit
circle to a non-convex shape.

Consider $z(t)=x(t)+\iota y(t)$.

$z(t)={\frac{1}{\Gamma(\alpha)}}\sum_{j=1}^{t}g_{\alpha}(t-j)\left(((\delta-1)-\iota
q)x(j-1)+\iota(\delta-1-\iota q)y(j-1)\right).\\\ $ (15)

Thus

$z(t)={\frac{1}{\Gamma(\alpha)}}\sum_{j=1}^{t}g_{\alpha}(t-j)(\delta-1)-\iota
q)z(j-1).$ (16)

Similarly, if we define $\bar{z}(t)=x(t)-\iota y(t)$, we get

$\bar{z}(t)={\frac{1}{\Gamma(\alpha)}}\sum_{j=1}^{t}g_{\alpha}(t-j)(\delta-1)+\iota
q)\bar{z}(j-1).$ (17)

The variables $z$ and $\bar{z}$ can be associated to ordered pairs $(x,y)$ and
$(x,-y)$ describing a complex number and its conjugate in complex plane.

The equivalent continuous-time system of (13)–(14) is given by

$\displaystyle D^{\alpha}x$ $\displaystyle=$ $\displaystyle(\delta-1)x+qy$
(18) $\displaystyle D^{\alpha}y$ $\displaystyle=$ $\displaystyle-
qx+(\delta-1)y.$ (19)

The general solution of system (18)–(19) is

$x(t)+\iota y(t)=E_{\alpha}\left((\delta-1-\iota
q)t^{\alpha}\right)\left(x(0)+\iota y(0)\right).$ (20)

If we define $z(t)=x(t)+\iota y(t)$ and $\lambda^{\prime}_{e}=(\delta-1-\iota
q)$, we have

$z(t)=E_{\alpha}\left(\lambda^{\prime}_{e}t^{\alpha}\right)z(0).$ (21)

This motivates us to define the effective Lyapunov exponent
$\lambda^{\prime}_{e}$ of the above system as

$\lambda^{\prime}_{e}=\lim_{t\longrightarrow\infty}t^{-\alpha}E_{\alpha}^{-1}\left(\frac{z(t)}{z(0)}\right).$
(22)

This is a well-defined quantity and we always find
$\lambda^{\prime}_{e}=\delta-1-\iota q$ numerically. This is a complex number.
On the other hand, Lyapunov exponents for dynamical systems have real Lyapunov
exponents. We believe that a single real number can determine stability for
difference equation because the stability condition for integer-order
differential equation is that the real part of eigenvalues is less than zero.
This is a single condition. On the other, for fractional order differential
equations, stability condition is given by two lines. Thus a single number may
not be enough to determine stability of fractional equations. However, if we
insist that the effective Lyapunov exponents should be real, we can take
$-|\lambda^{\prime}_{e}|$ as an effective Lyapunov exponent. Alternatively, we
may define the effective Lyapunov exponent $\lambda_{e}$ of system (18)–(19)
as

$\lambda_{e}=\lim_{t\longrightarrow\infty}\left|t^{-\alpha}E_{\alpha}^{-1}\left(\frac{x(t)+\iota
y(t)}{x(0)+\iota y(0)}\right)\right|=\sqrt{(\delta-1)^{2}+q^{2}}.$ (23)

Note that,
$\frac{\sqrt{x(t)^{2}+y(t)^{2}}}{\sqrt{x(0)^{2}+y(0)^{2}}}=|E_{\alpha}\left((\delta-1-\iota
q)t^{\alpha}\right)|$. In general, it is not possible to find the inverse of
composite function on right side. However, for $\alpha\in[0.5,1]$ and
$(\delta,q)$ in the stability region we observed that the quantity
$\frac{\sqrt{x(t)^{2}+y(t)^{2}}}{\sqrt{x(0)^{2}+y(0)^{2}}}$ decays with
$\lambda_{e}$ and we have
$\lambda_{e}\approx\lim_{t\longrightarrow\infty}t^{-\alpha}E_{\alpha}^{-1}\left(\frac{\sqrt{x(t)^{2}+y(t)^{2}}}{\sqrt{x(0)^{2}+y(0)^{2}}}\right)$.

In Fig. 6 we have plotted effective Lyapunov found numerically for various
values of $\delta$, $q$ and $\alpha$ as a function of above quantity and it is
clear that there is an excellent match.

Figure 6: Effective Lyapunov exponent $\lambda_{e}$ is plotted as function of
expression $\sqrt{(\delta-1)^{2}+q^{2}}$.

##  4 Generalization to nonlinear case

Let us try to extend the analysis to nonlinear systems. Let us consider
stability of the fixed point of a nonlinear system. Nonlinear system can be
linearized by $f(x)=f(x^{*})+f^{\prime}(x^{*})(x-x^{*})+\ldots$. We define
$y=x-x^{*}$. Now $x_{n+1}=f(x_{n})$ is equivalent to
$y_{n+1}=x_{n+1}-x^{*}=f(x_{n})-x^{*}=f^{\prime}(x^{*})y_{n}=ry_{n}$ where
$r=f^{\prime}(x^{*})$. Thus we can conjecture that the stability regime for
the fixed point is given by $a_{c}(\alpha)<f^{\prime}(x^{*})<1$.
Interestingly, this conjecture works. Consider the cases of Bernoulli map and
Gauss map considered by Deshpande and Gejji [7].

For Bernouli map $f(x)=rx|_{mod1},r>0$, we can guess that for $r<1$, the fixed
point is stable for any value of $\alpha$. This is precisely what Deshpande
and Gejji [7] find. In the case of Gauss map given by
$f(x)=\exp(-7.5x^{2})+\beta$, the fixed point can be found numerically using
bisection method or other methods. The fixed point is close to $\beta$. The
slope is negative and the stability regime is dependent on value of $\alpha$.
One expects this fixed point to be more stable if $\alpha$ increases. It can
be checked that $a_{c}(0.4)=-0.3132$, $a_{c}(0.6)=-0.5085$ and
$a_{c}(0.8)=-0.7328$. The critical points corespond to values of $\beta$ at
which the $f^{\prime}(x^{*})$ matches with these values. Thus we expect the
fixed point to be stable for $\beta>0.49$ for $\alpha=0.8$, for $\beta>0.57$
for $\alpha=0.6$ and $\beta>0.65$ for $\alpha=0.4$ (cf. Fig. 7).

This can be confirmed from the above paper [7] as well as independent
simulations. We start with $\beta$ as an initial condition.

Figure 7: a)Bifurcation diagram for Gauss map as a function of $\beta$ for
various values of $\alpha$. b) Local slope of fixed point $x^{*}$ is plotted
as a function of $\beta$. Vertical lines show the values of $\beta$ at which
$f^{\prime}(x^{*})=a_{c}(\alpha)$ for these values of $\alpha$. The fixed
point is stable in the expected region even for a nonlinear function.

##  5 Discussion and Conclusions

Analytic studies in nonlinear dynamics or coupled nonlinear systems are often
based on local linearization of dynamics. Thus linear systems serve as a basis
which helps undertanding dynamics in these systems. We have studied the
fractional equivalent of linear maps, which are not studied before to our
knowledge. Our studies indicate that the results obtained are useful in
studies of nonlinear systems as well. For coupled systems, we have studied the
stability of the fixed point. We also find that the conditions for chaotic
synchronization and antisynchronization and find that the conditions are very
similar to those obtained for coupled integer order maps.

For one-dimensional $f(x)=rx$, we find that the stability regime is given by
$a_{c}(\alpha)<r<1$. In the stable regime the dynamics is governed by Mittag-
Leffler function. We also define the effective Lyapunov exponent and find that
$r-1$ is effective Lyapunov exponent in this case.

For two coupled linear systems, the behavior is different for symmetric and
antisymmetric coupling. The analysis is motivated by study of linear
difference equations which is essentially the theory of matrices. We can
reduce dynamics to ‘normal modes’ which could be real or complex. We can find
effective Lyapunov exponents in these cases as well. When the normal modes are
real, the stability condition is the same as the condition for a single linear
map for each mode. A much richer picture is observed for complex normal modes.
The effective Lyapunov exponents, in this case, are complex which is not
entirely unexpected.

The stability region of a continuous-time fractional-order dynamical system is
a superset of that of classical integer-order one. As we increase the
fractional order to 1, the cone-like stability region of the fractional case
gets contracted to the left-half complex plane and we get the usual region of
stability of classical case. In this article, we showed that the stability
properties of discrete-time systems are different. The cardioid-like stability
region of fractional order system gets deformed and converted to the unit disc
as we increase the order to 1.

We also extend this work to fixed points of nonlinear maps and confirm that a
similar criterion holds. This work can be extended in many directions. We can
try to find the stability of periodic orbits of the higher period and find
routes to chaos in low-dimensional fractional difference equations.

## Acknowledgements

PMG thanks DST-SERB for financial assistance (Ref. EMR/2016/006686). SBB
acknowledges the Science and Engineering Research Board (SERB), New Delhi,
India for the Research Grant (Ref. MTR/2017/000068) under Mathematical
Research Impact Centric Support (MATRICS) Scheme.

## References

  * [1] R. M. May, “Biological populations obeying difference equations: stable points, stable cycles, and chaos,” Journal of Theoretical Biology, vol. 51, no. 2, pp. 511–524, 1975.
  * [2] S. H. Strogatz, Nonlinear dynamics and chaos with student solutions manual: With applications to physics, biology, chemistry, and engineering.  CRC press, 2018.
  * [3] M. J. Feigenbaum, “The universal metric properties of nonlinear transformations,” Journal of Statistical Physics, vol. 21, no. 6, pp. 669–706, 1979.
  * [4] M. J. Feigenbaum, “Quantitative universality for a class of nonlinear transformations,” Journal of statistical physics, vol. 19, no. 1, pp. 25–52, 1978.
  * [5] E. Ott, Chaos in dynamical systems.  Cambridge university press, 2002.
  * [6] S. S. Pakhare, V. Daftardar-Gejji, D. S. Badwaik, A. Deshpande, and P. M. Gade, “Emergence of order in dynamical phases in coupled fractional gauss map,” Chaos, Solitons & Fractals, vol. 135, p. 109770, 2020.
  * [7] A. Deshpande and V. Daftardar-Gejji, “Chaos in discrete fractional difference equations,” Pramana, vol. 87, no. 4, p. 49, 2016.
  * [8] K. S. Miller and B. Ross, “Fractional difference calculus.,” in Proceedings of the International Symposium on Univalent Functions, Fractional Calculus and Their Applications, Nihon University, Koriyama, Japan,, p. 139, 1989\.
  * [9] F. M. Atıcı and S. Şengül, “Modeling with fractional difference equations,” Journal of Mathematical Analysis and Applications, vol. 369, no. 1, pp. 1–9, 2010.
  * [10] F. M. Atici and P. W. Eloe, “A transform method in discrete fractional calculus,” International Journal of Difference Equations, vol. 2, no. 2, 2007.
  * [11] M. T. Holm, “The theory of discrete fractional calculus: Development and application,” 2011.
  * [12] G.-C. Wu and D. Baleanu, “Chaos synchronization of the discrete fractional logistic map,” Signal processing, vol. 102, pp. 96–99, 2014.
  * [13] J. J. Mohan and G. Deekshitulu, “Fractional order difference equations,” International Journal of Differential Equations, vol. 2012, 2012.
  * [14] Y. Liu, “Discrete chaos in fractional hiénon maps,” International Journal of Nonlinear Science, vol. 2018, pp. 170–175, 2014.
  * [15] P. M. Gade, “Synchronization of oscillators with random nonlocal connectivity,” Physical Review E, vol. 54, no. 1, p. 64, 1996.