Robust Counterfactual Explanations for Neural Networks With
  Probabilistic Guarantees
There is an emerging interest in generating robust counterfactual
explanations that would remain valid if the model is updated or changed even
slightly. Towards finding robust counterfactuals, existing literature often
assumes that the original model $m$ and the new model $M$ are bounded in the
parameter space, i.e., $\|\text{Params}(M){-}\text{Params}(m)\|{<}\Delta$.
However, models can often change significantly in the parameter space with
little to no change in their predictions or accuracy on the given dataset. In
this work, we introduce a mathematical abstraction termed
$\textit{naturally-occurring}$ model change, which allows for arbitrary changes
in the parameter space such that the change in predictions on points that lie
on the data manifold is limited. Next, we propose a measure -- that we call
$\textit{Stability}$ -- to quantify the robustness of counterfactuals to
potential model changes for differentiable models, e.g., neural networks. Our
main contribution is to show that counterfactuals with sufficiently high value
of $\textit{Stability}$ as defined by our measure will remain valid after
potential $\textit{naturally-occurring}$ model changes with high probability
(leveraging concentration bounds for Lipschitz function of independent
Gaussians). Since our quantification depends on the local Lipschitz constant
around a data point which is not always available, we also examine practical
relaxations of our proposed measure and demonstrate experimentally how they can
be incorporated to find robust counterfactuals for neural networks that are
close, realistic, and remain valid after potential model changes. This work
also has interesting connections with model multiplicity, also known as, the
Rashomon effect.
Faisal Hamman, Erfaun Noorani, Saumitra Mishra, Daniele Magazzeni, Sanghamitra Dutta
# Robust Counterfactual Explanations for Neural Networks

###### Abstract

There is an emerging interest in generating robust counterfactual explanations that would remain valid if the model is updated or changed even slightly. Towards finding robust counterfactuals, existing literature often assumes that the original model \(m\) and the new model \(M\) are bounded in the parameter space, i.e., \(\|\text{Params}(M)-\text{Params}(m)\|<\Delta\). However, models can often change significantly in the parameter space with little to no change in their predictions or accuracy on the given dataset. In this work, we introduce a mathematical abstraction termed _naturally-occurring_ model change, which allows for arbitrary changes in the parameter space such that the change in predictions on points that lie on the data manifold is limited. Next, we propose a measure - that we call _Stability_ - to quantify the robustness of counterfactuals to potential model changes for differentiable models, e.g., neural networks. Our main contribution is to show that counterfactuals with sufficiently high value of _Stability_ as defined by our measure will remain valid after potential "naturally-occurring" model changes with high probability (leveraging concentration bounds for Lipschitz function of independent Gaussians). Since our quantification depends on the local Lipschitz constant around a data point which is not always available, we also examine practical relaxations of our proposed measure and demonstrate experimentally how they can be incorporated to find robust counterfactuals for neural networks that are close, realistic, and remain valid after potential model changes. This work also has interesting connections with model multiplicity, also known as, the Rashomon effect.

## 1 Introduction

Counterfactual explanations (Wachter et al., 2017; Karimi et al., 2020; Barocas et al., 2020) have garnered significant interest in various high-stakes applications, such as lending, hiring, etc. Counterfactual explanations aim to guide an applicant on how they can change a model outcome by providing suggestions for improvement. Given an original data-point (e.g., an applicant who is denied a loan), the goal is to try to find a point on the other (desired) side of the decision boundary (a hypothetical applicant who is approved for the loan) which also satisfies several other preferred constraints, such as, (i) proximity to the original point; (ii) changes in as few features as possible; and (iii) conforming to the data manifold. Such a data-point that alters the model decision is widely referred to as a "counterfactual."

However, in several real-world scenarios, such as credit lending, the models making these high-stakes decisions have to be updated due to various reasons (Upadhyay et al., 2021; Black et al., 2021), e.g., to retrain on a few additional data points, change the hyper-parameters or seed, or transition to a different model class (Pawelczyk et al., 2020). Such model changes can often cause the counterfactuals to become invalid because typically they are quite close to the original data point, and hence, also quite close to the decision boundary. For instance, suppose the counterfactual explanation suggests an applicant to increase their income by \(10K\) to get approved for a loan and they actually act upon that, but now, due to updates to the original model, they are still denied by the updated model.

If counterfactuals become invalid due to model updates, this can lead to confusion and distrust in the use of algorithms in high-stakes applications altogether. Users would typically act on the suggested counterfactuals over a period of time, e.g., increase their income for credit lending, but only to find that it is no longer enough since the model has slightly changed (perhaps due to retraining with a new seed or hyperparameter). This cycle of invalidation and regenerating new counterfactuals can not only be frustrating and time-consuming for users but also potentially hurts an institution's reputation.
This motivates our primary question: _How do we provide theoretical guarantees on the robustness of counterfactuals to potential model changes?_

Towards addressing this question, in this work, we introduce the abstraction of "naturally-occurring" model change for differentiable models. Our abstraction allows for arbitrary changes in the parameter space such that the change in predictions on points that lie on the data manifold is limited. This abstraction motivates a measure of robustness for counterfactuals that arrives with provable probabilistic guarantees on their validity under naturally-occurring model change. We also introduce the notion of "targeted" model change and provide an impossibility result for such model change. Next, by leveraging a _computable relaxation_ of our proposed measure of robustness, we then design and implement algorithms to find robust counterfactuals for neural networks. Our experimental results validate our theoretical understanding and illustrate the efficacy of our proposed algorithms. We summarize our contributions here:

* **Abstraction of "naturally-occurring" model change for differentiable models:** Existing literature (Upadhyay et al., 2021; Black et al., 2021) on robust counterfactuals often assumes that the original model \(m\) and the new model \(M\) are bounded in the parameter space, i.e., \(\|\text{Params}(M)-\text{Params}(m)\|{<}\Delta\). Building on Dutta et al. (2022) for tree-based models, we note that models can often change significantly in the parameter space with little to no change on their predictions or accuracy on the given dataset. To capture this, we introduce an abstraction (see Definition 5), that we call _naturally-occurring_ model change, which instead allows for arbitrary changes in the parameter space such that the change in predictions on points that lie on the data manifold is limited. Our proposed abstraction of naturally-occurring model change also has interesting connections with predictive/model multiplicity, also known as, the Rashomon Effect (Breiman, 2001; Marx et al., 2020).
* that we call Stability
- to quantify the robustness of counterfactuals to potential model changes for differentiable models. Stability of a counterfactual \(x\in\mathbb{R}^{d}\) with respect to a model \(m(\cdot)\) is given by: \[R_{k,\sigma^{2}}(x,m)=\frac{1}{k}\sum_{x_{i}\in N_{x,k}}\left(m(x_{i})-\gamma_ {x}\|x-x_{i}\|\right),\] where \(N_{x,k}\) is a set of \(k\) points in \(\mathbb{R}^{d}\) drawn from the Gaussian distribution \(\mathcal{N}(x,\sigma^{2}\text{I}_{d})\) with \(\text{I}_{d}\) being the identity matrix, and \(\gamma_{x}\) is the local Lipschitz constant of the model \(m(\cdot)\) around \(x\) (Definition 6). Our main contribution in this work is to provide a theoretical guarantee (see Theorem 1) that counterfactuals with a sufficiently high value of Stability (as defined by our measure) will remain valid with high probability after "naturally-occurring" model change. Our result leverages concentration bounds for Lipschitz functions of independent Gaussian random variables (see Lemma 3). Since our proposed Stability measure depends on the local Lipschitz constant which is not always available, we also examine practical relaxations of our measure of the form: \[\hat{R}_{k,\sigma^{2}}(x,m)=\frac{1}{k}\sum_{x_{i}\in N_{x,k}}\left(m(x_{i})
- \left|m(x)-m(x_{i})\right|\right).\] The first term essentially captures the mean value of the model output in a region around it (higher mean is expected to be more robust and reliable). The second term captures the local variability of the model output in around it (lower variability is expected to be more reliable). This intuition is in alignment with the results in Dutta et al. (2022) for tree-based models.
* **Impossibility under targeted model change:** We also make a clear distinction between our proposed naturally-occurring and targeted model change. Under targeted model change, we provide an impossibility result (see Theorem 2) that given any counterfactual for a model, one can always design a new model that is quite similar to the original model and that renders that particular counterfactual invalid. However, in this work, our focus is on non-targeted model change such as retraining on a few additional data points, changing some hyperparameters or seed, etc., for which we have defined the abstraction of "naturally-occurring" model change (see Definition 5).
* **Experimental results:** We explore methods for incorporating our relaxed measure into generating robust counterfactuals for neural networks. We introduce T-Rex:I (Algorithm 1), which finds robust counterfactuals that are close to the original data point. T-Rex:I can be integrated into any base technique for generating counterfactuals to improve robustness. We also propose T-Rex:NN (Algorithm 2), which generates robust counterfactuals that are data-supported, making them more realistic (along the lines of Dutta et al. (2022) for tree-based models). Our experiments show that T-Rex:I can improve robustness for neural networks without significantly increasing the cost, and T-Rex:NN consistently generates counterfactuals that are similar to the data manifold, as measured using the local outlier factor (LOF).

### Related Works

Counterfactual explanations have seen growing interest in recent years (Verma et al., 2020; Karimi et al., 2020; Wachter et al., 2017). Regarding their robustness to model changes, Pawelczyk et al. (2020); Kanamori et al. (2020); Poyiadzi et al. (2020) argue that counterfactuals situated on the data
manifold are more likely to be more robust than the closest counterfactuals. Later, Dutta et al. (2022) demonstrate that generating counterfactuals on the data manifold may not be sufficient for robustness. While the importance of robustness in local explanation methods has been emphasized (Hancox-Li, 2020), the problem of specifically generating robust counterfactuals has been less explored, with the notable exceptions of some recent works (Upadhyay et al., 2021; Rawal et al., 2020; Black et al., 2021; Dutta et al., 2022; Jiang et al., 2022). In Upadhyay et al. (2021), the authors propose an algorithm called ROAR that uses min-max optimization to find the _closest_ counterfactuals that are also robust. In Rawal et al. (2020), the focus is on analytical trade-offs between validity and cost. Jiang et al. (2022) introduces a method for identifying close and robust counterfactuals based on a framework that utilizes interval neural networks. Black et al. (2021) propose that local Lipschitzness can be leveraged to generate consistent counterfactuals and propose an algorithm called Stable Neighbor Search to generate consistent counterfactuals for neural networks. Our research builds on this perspective and further performs Gaussian sampling around the counterfactual, leading to a novel estimator for which we are also able to provide probabilistic guarantees going beyond the bounded model change assumption. Furthermore, examining all three performance metrics, namely, cost, validity (robustness), and likeness to the data-manifold has received less attention with the notable exception of Dutta et al. (2022) but they focus only on tree-based models (non-differentiable). We also refer to Mishra et al. (2021) for a survey.

We note that Laugel et al. (2019); Alvarez-Melis and Jaakkola (2018) propose an alternate perspective of robustness in explanations (called \(L\)-stability in Alvarez-Melis and Jaakkola (2018)) which is built on similar individuals receiving similar explanations. Pawelczyk et al. (2022); Maragno et al. (2023); Dominguez-Olmedo et al. (2022) focus on finding counterfactuals that are robust to small input perturbations (noisy counterfactuals). In contrast, our focus is on counterfactuals remaining valid after some changes to the model, and providing theoretical guarantees thereof.

Our work also shares interesting conceptual connections with a body of work on model multiplicity or predictive multiplicity, also known as the Rashomon effect (Breiman, 2001; Marx et al., 2020; Black et al., 2022; Hsu and Calmon, 2022). Breiman (2001) suggested that models can be very different from each other but have almost similar performance on the data manifold. The term predictive multiplicity was suggested by Marx et al. (2020), defining it as the ability of a prediction problem to admit competing models with conflicting predictions and introduce formal measures to evaluate the severity of predictive multiplicity. Black et al. (2022) investigates ways to leverage model multiplicity beneficially in model selection processes while simultaneously addressing its concerning implications. Watson-Daniels et al. (2023) offered a framework for measuring predictive multiplicity in classification, introducing measures that encapsulate the variation in risk estimates over the ensemble of competing models. Hsu and Calmon (2022) unveiled a novel metric, Rashomon Capacity, for measuring predictive multiplicity in probabilistic classification. Our proposed abstraction of naturally-occurring model change, as explored in this work, can be viewed as a fresh perspective on model multiplicity that further emphasizes on the models that are more likely to occur.

## 2 Preliminaries

Let \(\mathcal{X}\subseteq\mathbb{R}^{d}\) denote the input space and let \(\mathcal{S}{=}\{x_{i}\in\mathcal{X}\}_{i=1}^{n}\) be a dataset consisting of \(n\) independent and identically distributed data points generated from a density \(q\) over \(\mathcal{X}\). We also let \(m(\cdot):\mathbb{R}^{d}\rightarrow[0,1]\) denote the original machine learning model that takes a \(d\)-dimensional input value and produces an output probability lying between \(0\) and \(1\). The final decision is denoted by \(\mathbb{1}(m(x)\geq 0.5)\) where \(\mathbb{1}(\cdot)\) denotes the indicator function.

**Definition 1** (\(\gamma-\)Lipschitz).: _A function \(m(\cdot)\) is said to be \(\gamma-\)Lipschitz if \(|m(x)-m(x^{\prime})|{\leq}\gamma\|x-x^{\prime}\|\) for all \(x,x^{\prime}{\in}\mathbb{R}^{d}\)._

Here \(\|\cdot\|\) denotes the Euclidean norm, i.e., for \(u\in\mathbb{R}^{d}\), we have \(\|u\|=\sqrt{u_{1}^{2}+u_{2}^{2}+\ldots+u_{d}^{2}}\). In Remark 2, we also discuss relaxations to local Lipschitz constants from global Lipschitz constants. We denote the updated or changed model as \(M(\cdot):\mathbb{R}^{d}\rightarrow[0,1]\) where \(M\) is a random entity. We mostly use capital letters to denote random entities, e.g., \(M\), \(X\), etc., and small letters to denote non-random entities, e.g., \(m\), \(x\), \(\gamma\), \(n\), etc.

### Background on Counterfactuals

**Definition 2** (Closest Counterfactual \(\mathcal{C}_{p}(x,m)\)).: _Given \(x\in\mathbb{R}^{d}\) such that \(m(x)<0.5\), its closest counterfactual (in terms of \(l_{p}\)-norm) with respect to the model \(m(\cdot)\) is defined as a point \(x^{\prime}\in\mathbb{R}^{d}\) that minimizes the \(l_{p}\) norm \(\|x-x^{\prime}\|_{p}\) such that \(m(x^{\prime})\geq 0.5\)._

\[\mathcal{C}_{p}(x,m)=\arg\min_{x^{\prime}\in\mathbb{R}^{d}}\|x-x^{\prime}\|_{ p}\text{ such that }m(x^{\prime})\geq 0.5.\]

When one is interested in finding counterfactuals by changing as few features as possible, the \(l_{1}\) norm is used (enforcing a sparsity constraint). These counterfactuals are also called _sparse_ counterfactuals (Pawelczyk et al., 2020).

Closest counterfactuals often fall too far from the data manifold, resulting in unrealistic and anomalous instances, as noted in Poyiadzi et al. (2020); Pawelczyk et al. (2020); Kanamori et al. (2020); Verma et al. (2020); Karimi et al. (2020); Albini et al. (2022). This highlights the need for generating counterfactuals that lie on the data manifold.
**Definition 3** (Closest Data-Manifold Counterfactual \(\mathcal{C}_{p,\mathcal{X}}(x,m)\)).: _Given \(x\in\mathbb{R}^{d}\) such that \(m(x)<0.5\), its closest data-manifold counterfactual \(\mathcal{C}_{p,\mathcal{X}}(x,m)\) with respect to the model \(m(\cdot)\) and data manifold \(\mathcal{X}\) is defined as a point \(x^{\prime}\in\mathcal{X}\) that minimizes the \(l_{p}\) norm \(\|x-x^{\prime}\|_{p}\) such that \(m(x^{\prime})\geq 0.5\)._

\[\mathcal{C}_{p,\mathcal{X}}(x,m)=\arg\min_{x^{\prime}\in\mathcal{X}}\|x-x^{ \prime}\|_{p}\text{ such that }m(x^{\prime})\geq 0.5.\]

In order to assess the similarity or anomalous nature of a point concerning the given dataset \(\mathcal{S}\subseteq\mathcal{X}\), various metrics can be employed, e.g., K-nearest neighbors, Mahalanobis distance, Kernel density, LOF. These metrics play a crucial role in understanding the quality of counterfactual explanations generated by a model. One widely used metric in the literature on counterfactual explanations (Pawelczyk et al., 2020; Kanamori et al., 2020; Dutta et al., 2022) is the Local Outlier Factor (LOF).

**Definition 4** (Local Outlier Factor (Breunig et al., 2000)).: _For \(x\in\mathcal{S}\), let \(L_{k}(x)\) be its \(k\)-Nearest Neighbors (\(k\)-NN) in \(\mathcal{S}\). The \(k\)-reachability distance \(rd_{k}\) of \(x\) with respect to \(x^{\prime}\) is defined by \(rd_{k}(x,x^{\prime})=\max\{\delta(x,x^{\prime}),d_{k}(x^{\prime})\}\), where \(d_{k}(x^{\prime})\) is the distance \(\delta\) between \(x^{\prime}\) and its the \(k\)-th nearest instance on \(\mathcal{S}\). The \(k\)-local reachability density of \(x\) is defined by \(lrd_{k}(x)=|L_{k}(x)|(\sum_{x^{\prime}\in L_{k}(x)}rd_{k}(x,x^{\prime}))^{-1}\). Then, the k-LOF of \(x\) on \(\mathcal{S}\) is defined as follows:_

\[LOF_{k,\mathcal{S}}(x)=\frac{1}{|L_{k}(x)|}\sum_{x^{\prime}\in L_{k}(x)}\frac {lrd_{k}(x^{\prime})}{lrd_{k}(x)}.\]

_Here, \(\delta(x,x^{\prime})\) is the distance between two \(d\)-dimensional feature vectors. The LOF Predicts \(-1\) for anomalous points and \(+1\) for inlier points._

**Goals:** In this work, our main goal is to provide _probabilistic guarantees_ on the robustness of counterfactuals to potential model changes for differential models such as neural networks. Towards achieving this goal, our objective involves: (i) introducing an abstraction that rigorously defines the class of model changes that we are interested in; and (ii) establishing a measure, denoted as \(R_{\Phi}(x,m)\), for a counterfactual \(x\) and a given model \(m(\cdot)\), that quantifies its robustness to potential model changes. Here, \(\Phi\) represents the hyperparameters of the robustness measure. Ideally, we desire that the measure \(R_{\Phi}(x,m)\) should be high if the counterfactual \(x\) is less likely to be invalidated by potential model changes. We seek to provide: (i) theoretical guarantees on the validity of counterfactuals with sufficiently high value of \(R_{\Phi}(x,m)\); and also (ii) incorporate our measure \(R_{\Phi}(x,m)\) into an algorithmic framework for generating robust counterfactuals which also meet other requirements, such as, low cost or likeness to the data manifold.

## 3 Main Theoretical Contributions

In this section, we first introduce our proposed abstraction of _naturally-occurring_ model change and then propose a novel measure - that we call _Stability_ - to quantify the robustness of counterfactuals to potential model changes. We derive a theoretical guarantee that counterfactuals that have a sufficiently high value of _Stability_ will remain valid after potential _naturally-occurring_ model change with high probability. But since our quantification would depend on the local Lipschitz constant around a data point, which is not always known, we also examine a practical relaxation of our proposed measure and demonstrate its applicability.

### Naturally-Occurring Model Change

A popular assumption in existing literature (Upadhyay et al., 2021; Black et al., 2021) to quantify potential model changes is to assume that the model changes are bounded in the parameter space, i.e.,

\[\|\text{Params}(M)-\text{Params}(m)\|<\Delta\text{ for a constant }\Delta.\]

Here, \(\text{Params}(M)\) denote the parameters of the model \(M\), e.g., weights of a neural network. However, we note that

Figure 1: Models can often change drastically in the parameter space causing little to no change in the actual decisions on the points on the data manifold.

Figure 2: Illustrates our proposed abstraction of naturally-occurring model change: The distribution of the changed model outputs \(M(x)\) (stochastic) is centered around the original model output \(m(x)\). The points specifically lying on the data-manifold acting as anchors without much change as they exhibit lower variance in model outputs compared to points outside the manifold. This visualization also connects with the Rashomon effect, encapsulating the diverse yet similarly accurate models that can be learned from a given dataset.
models can often change drastically in the parameter space causing little to no change on the actual decisions on the points on the data manifold (see Fig. 1 for an example). In this work, we relax the bounded-model-change assumption, and instead introduce the notion of a naturally-occurring model change as defined in Definition 5. Our abstraction allows for arbitrary model changes such that the change in predictions on points that lie on the data manifold is limited (see Fig. 2 for an illustration).

**Definition 5** (Naturally-Occurring Model Change).: _A naturally-occurring model change is defined as follows:_

1. \(\mathbb{E}\left[M(X)|X=x\right]=\mathbb{E}\left[M(x)\right]=m(x)\) _where the expectation is over the randomness of_ \(M\) _given a fixed value of_ \(X=x\in\mathbb{R}^{d}\)_._
2. _Whenever_ \(m(x)\) _is_ \(\gamma_{m}\)_-Lipschitz, any updated model_ \(M(x)\) _is also_ \(\gamma-\)_Lipschitz for some constant_ \(\gamma\)_. Note that, this constant_ \(\gamma\) _does not depend on_ \(M\) _since we may define_ \(\gamma\) _to be an upper bound on the Lipschitz constants for all possible_ \(M\) _as well as_ \(m\)_._
3. \(\mathrm{Var}\left[M(X)|X=x\right]=\mathrm{Var}\left[M(x)\right]=\nu_{x}\) _which depends on the fixed value of_ \(X=x\in\mathbb{R}^{d}\)_. Furthermore, whenever_ \(x\) _lies on the data manifold_ \(\mathcal{X}\)_, we have_ \(\nu_{x}\leq\nu\) _for a small constant_ \(\nu\)_._

Closely connected to naturally-occurring model change is the idea of the Rashomon effect, alternatively known as predictive or model multiplicity. (Breiman, 2001; Pawelczyk et al., 2020; Marx et al., 2020; Hsu and Calmon, 2022) which suggests that models can be very different from each other but have almost similar performance on the data manifold, e.g., \(\frac{1}{n}\sum_{i=1}^{n}|M(x_{i})-m(x_{i})|\) is small when the points \(x_{i}\) lie on the data manifold. Under the naturally-occurring model change, this holds in expectation as follows:

**Lemma 1** (Connection to Roshomon Effect).: _For points \(x_{1},\ldots,x_{n}\in\mathcal{X}\) (lying on the data-manifold) under naturally-occurring model change, the following holds:_

\[\mathbb{E}\left[\frac{1}{n}\sum_{i=1}^{n}|M(x_{i})-m(x_{i})|\right]\leq\sqrt{ \nu}. \tag{1}\]

Thus, Definition 5 might be better suited over boundedness in the parameter space. Proof of Lemma 1 is in Appendix B.

**Remark 1** (Targeted Model Change).: _In contrast with naturally-occurring model change, we also introduce the notion of targeted model change (adversarial, worst-case) which essentially refers to a model change that is more deliberately targeted to make a particular counterfactual invalid. For example, one could have a new model \(M(x)=m(x)\) almost everywhere except at or around the targeted point \(x^{\prime}\), i.e., \(M(x^{\prime})=1-m(x^{\prime})\). See Section 3.3 for more details._

### A Measure of Robustness With Probabilistic Guarantees on Validity

#### 3.2.1 Proposed Measure: Stability

**Definition 6** (Stability).: _The stability of a counterfactual \(x\in\mathbb{R}^{d}\) is defined as follows:_

\[R_{k,\sigma^{2}}(x,m)=\frac{1}{k}\sum_{x_{i}\in N_{x,k}}\left(m(x_{i})-\gamma \|x-x_{i}\|\right), \tag{2}\]

_where \(N_{x,k}\) is a set of \(k\) points drawn from the Gaussian distribution \(\mathcal{N}(x,\sigma^{2}\mathrm{I}_{d})\) with \(\mathrm{I}_{d}\) being the identity matrix, and \(\gamma\) is an upper bound on the Lipschitz constant for all models \(M(\cdot)\) under naturally-occurring change._

**Remark 2** (Relaxations to local Lipschitz).: _While we prove our theoretical result (Theorem 1) with the global Lipschitz constant \(\gamma\), we can relax this to local Lipschitz constants \(\gamma_{x}\), around a given point \(x\). This is because we sample from a Gaussian centered around the point \(x\) and hence mainly capture the variability around \(x\). So most points will be very close to \(x\) but a few points can still lie far away. Potential extensions of our guarantees could apply to truncated Gaussian and uniform sampling methods, given their sub-Gaussian properties. This is because Lipschitz concentration inherently extends to sub-Gaussian random variables (Baraniuk et al., 2008)._

#### 3.2.2 Probabilistic Guarantee

**Theorem 1** (Probabilistic Guarantee).: _Let \(X_{1},X_{2},\ldots,X_{k}\) be \(k\) iid random variables with distribution \(\mathcal{N}(x,\sigma^{2}I_{d})\) and \(Z=\frac{1}{k}\sum_{i=1}^{k}(m(X_{i})-M(X_{i}))\). Suppose \(\left|\mathbb{E}\left[Z|M\right]-\mathbb{E}\left[Z\right]\right|<\epsilon^{\prime}\). Then, for any \(\epsilon>2\epsilon^{\prime}\), a counterfactual \(x\in\mathcal{X}\) under naturally-occurring model change satisfies:_

\[\Pr(M(x)\geq R_{k,\sigma^{2}}(x,m)-\epsilon)\geq 1-\exp\left(\frac{-k \epsilon^{2}}{8(\gamma_{m}+\gamma)^{2}\sigma^{2}}\right).\]

_The probability is over the randomness of both \(M\) and \(X_{i}^{\prime}\)s._

**Intuition Behind Our Result:** This stability metric (Definition 6) is a way to measure the robustness of counterfactuals that are subject to natural model changes (see Definition 5). The first term in the metric, represented by \(\frac{1}{k}\sum_{i=1}^{k}m(X_{i})\), captures the average model outputs for a group of points centered around the counterfactual \(x\). The second term, represented by \(\gamma\|x-X_{i}\|\), is an upper bound on the potential difference in outputs of any new model on the points \(x\) and \(X_{i}\) (Recall the Lipschitz property of \(M\) around the point \(x\)). Using our measure, the guarantee in Theorem 1 can be rewritten as:

\[\Pr\left(\frac{1}{k}\sum_{i=1}^{k}m(X_{i})- M(x)\leq\frac{\gamma}{k}\sum_{i=1}^{k}\|x-X_{i}\|+\epsilon\right)\\ \geq 1-\exp\left(\frac{-k\epsilon^{2}}{8(\gamma+\gamma_{m})^{2} \sigma^{2}}\right). \tag{3}\]
This form of the inequality allows for the following interpretation of Theorem 1: The distance between the output of the new model on an input \(x\), i.e., \(M(x)\), and the average prediction of the neighborhood of the given input by the old model, i.e., \(\frac{1}{k}\sum m(X_{i})\) is upper bounded by \(\epsilon\)-corrected, \(\gamma\) multiplied average distance of the datapoints within the neighborhood of the input \(x\), i.e., \(\frac{1}{k}\sum\|x-X_{i}\|\).

Proof Sketch:.: The complete proof of Theorem 1 is provided in Appendix C.1. Here, we include a proof sketch. Notice that, using the Lipschitz property of \(M(\cdot)\) around \(x\), we have \(M(x)\geq M(X_{i})-\gamma\|x-X_{i}\|\) for all \(X_{i}\). Thus,

\[M(x) \geq\frac{1}{k}\sum_{i=1}^{k}(M(X_{i})-\gamma\|x-X_{i}\|) \tag{4}\] \[\stackrel{{(a)}}{{\geq}}\frac{1}{k}\sum_{i=1}^{k}( m(X_{i})-\gamma\|x-X_{i}\|)-\epsilon, \tag{5}\]

where (a) holds from Lemma 2 with probability at least \(1-\exp\Big{(}\frac{-k\epsilon^{2}}{8(\gamma+\gamma_{m})^{2}\sigma^{2}}\Big{)}\).

**Lemma 2** (Deviation Bound).: _Let \(X_{1},X_{2},\ldots,X_{k}\sim\mathcal{N}(x,\sigma^{2}I_{d})\) and \(Z{=}\frac{1}{k}\sum_{i=1}^{k}(m(X_{i}){-}M(X_{i}))\). Suppose \(|\mathbb{E}\left[Z|M\right]-\mathbb{E}\left[Z\right]|<\epsilon^{\prime}\). Then, under naturally-occurring model change, we have \(\mathbb{E}\left[Z\right]{=}0\). Moreover, for any \(\epsilon{>}2\epsilon^{\prime},\)_

\[\Pr(Z\geq\epsilon)\leq\exp\bigg{(}\frac{-k\epsilon^{2}}{8(\gamma+\gamma_{m}) ^{2}\sigma^{2}}\bigg{)}. \tag{6}\]

Proof Sketch:.: The proof of Lemma 2 leverages concentration bounds for Lipschitz functions of independent Gaussian random variables (see Lemma 3). The complete proof of Lemma 2 is provided in Appendix C. 

**Lemma 3** (Gaussian Concentration Inequality).: _Let \(W=(W_{1},W_{2},\ldots,W_{n})\) consist of n i.i.d. random variables belonging to \(\mathcal{N}(0,\sigma^{2})\), and \(Z=f(W)\) be a \(\gamma\)-Lipschitz function, i.e., \(|f(W)-f(W^{\prime})|\leq\gamma\|W-W^{\prime}\|.\) Then, we have,_

\[\Pr(Z-\mathbb{E}\left[Z\right]\geq\epsilon)\leq\exp\bigg{(}\frac{-\epsilon^{ 2}}{2\gamma^{2}\sigma^{2}}\bigg{)}\text{ for all }\epsilon>0. \tag{7}\]

For the proof of Lemma 3 refer to Boucheron et al. (2013) in p.125. Our robustness guarantee (Theorem 1) essentially states that \(\Pr(M(x)\leq R_{k,\sigma^{2}}(x,m)-\epsilon)\leq\exp\big{(}\frac{-k\epsilon^{2 }}{8(\gamma+\gamma_{m})^{2}\sigma^{2}}\big{)}\) under naturally-occurring model change. For instance, if we find a counterfactual \(x\) such that \(R_{k,\sigma^{2}}(x,m)-\epsilon\) is greater or equal to \(0.5\), then \(M(x)\) would also be greater than \(0.5\) with high probability. The term \(\exp\big{(}\frac{-k\epsilon^{2}}{8(\gamma+\gamma_{m})^{2}\sigma^{2}}\big{)}\) decays with \(k\).

#### 3.2.3 Practical Relaxation of Stability and Its Properties

While our proposed measure Stability (Definition 6) has probabilistic guarantees, we note that it relies on the Lipschitz constant \(\gamma\) (or the local Lipschitz constant \(\gamma_{x}\) around the point \(x\)), which is often unknown. Therefore, we next propose a practical relaxation of the measure as follows:

**Definition 7** (Stability (Relaxed)).: _The stability (relaxed) of a counterfactual \(x\in\mathbb{R}^{d}\) is defined as follows:_

\[\hat{R}_{k,\sigma^{2}}(x,m)=\frac{1}{k}\sum_{x_{i}\in N_{x,k}}(m(x_{i})-|m(x)- m(x_{i})|),\]

_where \(N_{x,k}\) is a set of \(k\) points drawn from the Gaussian distribution \(\mathcal{N}(x,\sigma^{2}\mathrm{I}_{d})\) with \(\mathrm{I}_{d}\) being the identity matrix._

To arrive at this relaxation, we utilize the Lipschitz property to approximate the aspect that involves the Lipschitz constant, specifically, by approximating \(\gamma_{x}||x-x_{i}||\) with \(|m(x)-m(x_{i})|\). Another possibility is to consider an estimate of \(\gamma_{x}\) given by:

\[\hat{\gamma}_{x}=\max_{x_{i}\in N_{x,k}}\frac{|m(x)-m(x_{i})|}{\|x-x_{i}\|}. \tag{8}\]

We observed that the experimental results with both these stability estimates are in the same ballpark.

To gain a deeper understanding of the relaxed stability measure, we now consider some desirable properties of counterfactuals that make them more robust and then demonstrate that our proposed relaxation of Stability (Definition 7) satisfies those desirable properties. These properties are inspired from Dutta et al. (2022) which proposed these properties for tree-based ensembles. The first property is based on the fact that the output of a model \(m(x)\in[0,1]\) is expected to be higher if the model has more confidence in that prediction.

**Property 1**.: _For any \(x\in\mathbb{R}^{d}\), a higher value of \(m(x)\) makes it less likely to be invalidated due to model changes._

Having a high \(m(x)\) alone does not guarantee robustness, as local variability around \(x\) can make predictions less reliable. E.g., points with high \(m(x)\) near the decision boundary are also vulnerable to invalidation with model changes.

**Property 2**.: _An \(x{\in}\mathbb{R}^{d}\) is less likely to be invalidated if several points close to \(x\) (denoted by \(x^{\prime}\)) have a high value of \(m(x^{\prime})\)._

Counterfactuals may also be more likely to be invalidated if it lies in a highly variable region of the model output function. This is because the confidence of the model predictions in that region may be less reliable.

**Property 3**.: _An \(x\in\mathbb{R}^{d}\) is less likely to be invalidated if model outputs around \(x\) have low variability._
Our stability measure aligns with these desired properties. Given a point \(x\in\mathbb{R}^{d}\), it generates a set of \(k\) points centered around \(x\). The first term \(\frac{1}{k}\sum_{x^{\prime}\in N_{x,k}}m\left(x^{\prime}\right)\) is expected to be high if the model output value \(m(x)\) is high for \(x\) as well as several points close to \(x\). But the mean value of \(m(x^{\prime})\) around a point \(x\) may not always capture the variability in that region, hence, the second term of our stability measure, i.e., \(\frac{1}{k}\sum_{x^{\prime}\in N_{x,k}}|m(x)-m(x^{\prime})|\). This term captures the variability of the model output values in a region around \(x\).

It is worth noting that the variability term is only useful in conjunction with the mean term. This is because even points on the opposite side of the decision boundary can have varying levels of variability, regardless of whether \(m(x^{\prime})\) is less or greater than \(0.5\).

In Fig. 3, we provide an example on the synthetic moon dataset to observe the effect of our stability measure on naturally-changed models. Note that these changed models were realized from actual experiments by retraining with different weight initializations.

### Impossibility Under Targeted Model Change

In this work, we make a key distinction between naturally-occurring and targeted model changes. While we are able to provide probabilistic guarantees for naturally-occurring model change, we also demonstrate an impossibility result for targeted model change.

**Theorem 2** (Impossibility Under Targeted Change).: _Given a model and a counterfactual, one can design another similar model such that the particular targeted counterfactual can be invalidated._

What this result essentially demonstrates is that for a given model, one can design another similar model such that any particular targeted counterfactual can be invalidated. The proof relies on the possibility that one could have a new model \(M(x)=m(x)\) almost everywhere except at or around the targeted point \(x^{\prime}\), i.e., \(M(x^{\prime})=1-m(x^{\prime})\).

## 4 Generating Robust Counterfactuals using Our Proposed Measure: Stability

In this section, we examine two techniques of incorporating our proposed measure, Stability (relaxed; see Definition 7), for generating robust counterfactuals for neural networks.

To begin, along the lines of Dutta et al. (2022), we first define a counterfactual robustness test.

**Definition 8** (Counterfactual Robustness Test).: _A counterfactual \(x\in\mathbb{R}^{d}\) satisfies the robustness test if:_

\[\hat{R}_{k,\sigma^{2}}(x,m)\geq\tau. \tag{9}\]

Now, we would like to find a reasonable point \(x^{\prime}\) that receives a positive prediction from the model (essentially \(m(x^{\prime})\geq 0.5\)), while also satisfying the robustness test, \(\hat{R}_{k,\sigma^{2}}(x^{\prime},m)\geq\tau\). The threshold value of \(\tau\) can be adjusted based on the desired effective validity (recall Theorem 1). Hence, a larger threshold would likely ensure that the new model, \(M\), remains valid with high probability. In trying to find a reasonable point \(x^{\prime}\), one may strive to generate robust counterfactuals that are as close as possible to the original point. One might also want the generated counterfactuals to be as realistic as possible, i.e., lie on the data manifold. Toward that end, we propose two algorithms.

We propose Algorithm 1, T-Rex:I, which incorporates our measure to find robust counterfactuals that are close to the

Figure 3: Effect of stability measure on naturally-occurring model changes: (a) corresponds to the original data distribution and the trained model. (b)-(e) demonstrate some examples of changed models obtained on retraining with different weight initializations. One may notice that the model decision boundary is changing a lot in the sparse regions of the data-manifold (few data-points), possibly violating the bounded-parameter change assumption but the predictions on the dense regions of the data-manifold do not change much (in alignment with Rashomon effect). This motivates our proposed abstraction of naturally-occurring model change which allows for arbitrary changes in the parameter space with little change in the actual predictions on the dense regions of the data manifold. (f) demonstrates our proposed measure of stability \(\hat{R}_{k,\sigma^{2}}(x,m)\) (high mean model output, low variability, _almost_ like a Gaussian filter) for which we derive probabilistic guarantees on validity. In essence, we show that under the abstraction of naturally-occurring model change, the stability measure captures the reliable intersecting region of changed models with high probability. In the original model, we observe that certain non-robust regions (i.e., those caused by overfitting to certain data points in the original model) have higher local Lipschitz values and variability. Counterfactuals assigned to these regions (even if \(m(x)\) is high) would be invalidated in the changed models. The stability measure, which samples around a region, penalizes these higher local Lipschitz values.
original data point. T-Rex:I works with any preferred base method for generating counterfactuals. It evaluates the stability of the generated counterfactual and, if necessary, iteratively updates the generated counterfactual through a gradient ascent process until a robust counterfactual that meets the desired criteria is obtained.

**Remark 3** (Gradient of Stability).: _In Algorithm 1, we compute the gradient of \(R(x,m)\) with respect to \(x\) (not model parameters \(m\)). Such gradients w.r.t. \(x\) instead of m are also computed commonly in adversarial machine learning and also in feature-attributions for explainability. We use TensorFlow tf.GradientTape for automatic differentiation, which allows for the computation of gradients with respect to certain inputs._

To ensure that the counterfactuals are as realistic as possible, we also define the T-Rex:NN Counterfactual, which considers counterfactuals that lie within a dataset to avoid any unrealistic or anomalous results (see Algorithm 2).

**Definition 9** (Robust Nearest Neighbor Counterfactual).: _Given \(x\in\mathbb{R}^{d}\) such that \(m(x)<0.5\), its robust nearest neighbor counterfactual \(\mathcal{C}_{p,\mathcal{S}}^{(\tau)}(x,m)\) with respect to the model \(m(\cdot)\) and dataset \(\mathcal{S}\) is defined as another point \(x^{\prime}\in\mathcal{S}\) that minimizes the \(l_{p}\) norm \(\left\|x-x^{\prime}\right\|_{p}\) such that \(m\left(x^{\prime}\right)\geq 0.5\) and \(\hat{R}_{k,\sigma^{2}}(x^{\prime},m)\geq\tau\)._

The closest data-supported counterfactual serves as a reliable reference, as it inherently has a high Local Outlier Factor (LOF). The T-Rex:I algorithm may find counterfactuals with lower costs, but they may compromise on the LOF and result in unrealistic samples.

To address this, we propose Algorithm 2, T-Rex:NN, for finding data-supported counterfactuals. The algorithm first finds \(k\) nearest neighbors counterfactuals of \(x\) in dataset \(\mathcal{S}\), checks through each of them to see if they satisfy the robustness test, \(\hat{R}_{k,\sigma^{2}}(x^{\prime},m)\geq\tau\), and terminates once such a counterfactual is found.

```
Input: Model \(m(\cdot)\), Datapoint \(x\) with \(m(x)\)\(<\)\(0.5\), Dataset \(\mathcal{S}\), Algorithm parameters (\(K,\sigma^{2},k,\tau\)).  Find \(K\) nearest neighbor counterfactuals \(x^{\prime}_{i}\in\mathcal{S}\) to \(x\) with respect to model \(m(\cdot)\), i.e., \(\mathrm{NN}_{x}=(x^{\prime}_{1},x^{\prime}_{2},\ldots,x^{\prime}_{K})\). for\(x^{\prime}_{i}\in\mathrm{NN}_{x}\)do  Perform counterfactual robustness test on \(x^{\prime}_{i}\):  Check if \(\hat{R}_{k,\sigma^{2}}(x^{\prime}_{i},m)\geq\tau\) if counterfactual robustness test is satisfied: then  Output \(x^{\prime}_{i}\) and exit endif endfor  Output: No robust counterfactual found and exit
```

**Algorithm 2** T-Rex:N: Theoretically Robust EXplanations: Nearest Neighbor Version

## 5 Experiments

In this section, we present experimental results to demonstrate the effectiveness of our proposed measure in capturing robustness and then generating robust counterfactuals that remain valid after potential model changes. We illustrate how our proposed Algorithm 1 & 2 utilizes our stability measure to effectively generate robust counterfactuals.

**Datasets:** We conduct experiments on several benchmark datasets, namely, HELOC (FICO, 2018), German Credit, Cardiotocography (CTG), Adult (Dua & Graff, 2017), and Taiwanese Credit (Yeh & hui Lien, 2009). These have two classes, with one class representing the most favorable outcome, and the other representing the least desirable outcome for which we aim to generate counterfactuals. For simplicity, we normalize the features to lie between \([0,1]\).

**Performance Metrics:** Our metrics of interest are:

* Cost: Average \(l_{1}\) or \(l_{2}\) distance between counterfactuals \(x^{\prime}\) and original points \(x\).
* Validity (%): Percentage of counterfactuals that remain valid under the new model \(M\).
* LOF: Predicts \(-1\) for anomalous points, and \(+1\) for inliers. A high average LOF essentially suggests the points lie on the data manifold and hence more realistic, i.e., _higher is better_ (see Definition 4). We use an existing implementation from Scikit-Learn to compute the LOF.

**Methodology:** We begin by training a baseline neural network model and aim to find counterfactuals for data points with true negative predictions. To test the robustness of these counterfactual examples, we then train \(50\) new models (\(M\)) and evaluate the \(validity\) of the counterfactuals under different model change scenarios, which include: (i) Weight Initialization (WI): Retraining new models using the same hyperparameters but with different weight initialization by using different random seeds for each new model; and (ii)
Leave Out (LO): Retraining new models by randomly removing a small portion (\(1\%\)) of the training data each time (with replacement) as well as different weight initialization.

**Hyperparameter selection:** Our theoretical findings indicate that higher \(k\) improves robustness, but comes at the cost of increased computational cost. We determine \(k=1000\) was sufficient. The value of \(\sigma^{2}\) was determined by analyzing the variance of the features. In the dataset with the features between \([0,1]\), we found that a value of \(\sigma^{2}=0.01\) produced good results. The threshold \(\tau\) is a critical aspect of our method and can be adjusted based on the desired effective validity. A higher \(\tau\) value improves validity at the expense of \(l_{1}\) or \(l_{2}\) cost. See Appendix D for more details.

**Baseline:** We compare our approaches with established baselines. First, we find the min Cost (\(l_{1}\) and \(l_{2}\)) counterfactual (Wachter et al., 2017) and use it as our base method for generating counterfactuals. We then compare T-Rex:I to the Stable Neighbor Search (SNS) (Black et al., 2021) and Robust Algorithmic Recourse (ROAR) (Upadhyay et al., 2021). We evaluate the performance of our Robust Nearest Neighbor (Algorithm 2:T-Rex:NN) against the Nearest Neighbor (NN) counterfactuals (closest data-support robust counterfactual in Definition 9). We choose a value of \(\tau\) to get high validity and compare cost and LOF with baselines.

**Results:** Results for HELOC, German Credit, and CTG datasets are summarized in Table 1. Observe that the min Cost counterfactual is not robust to variations in the training data or weight initialization as expected. ROAR generates counterfactuals with high validity, albeit at the expense of a higher cost. Our proposed method, T-Rex:I, significantly improves the validity of the counterfactuals compared to the minimum cost. The T-Rex:I algorithm achieves comparable validity results to the SNS method for both types of model changes, and often accomplishes this with lower costs and higher LOF. This can be observed across all three datasets for both \(l_{1}\) and \(l_{2}\) cost metrics. The T-Rex:NN algorithm also significantly improves the validity of the counterfactuals compared to the traditional Nearest Neighbor (NN) method and maintains a high LOF. It comes at a price of increased cost, but the counterfactuals are guaranteed to be realistic since they are data supported. Refer to Appendix D for additional results for Adult and Taiwanese credit datasets.

**Ablation:** To evaluate the efficacy of our proposed stability measure, we conduct an ablation study on the German credit dataset. We first evaluate a robustness measure that solely relies on the model's prediction of the counterfactual, denoted as \(r(x^{\prime},m)=m(x^{\prime})\). We then examine a measure that only incorporates the mean, the average predictions for \(k\) points sampled from the distribution \(N(x^{\prime},\sigma^{2}I_{d})\), denoted as \(r_{k,\sigma^{2}}(x^{\prime},m)=\frac{1}{k}\sum_{x^{\prime}_{i}\in N_{x^{\prime },k}}m(x^{\prime}_{i})\). We compare these with our proposed robustness measure \(\hat{R}_{k,\sigma^{2}}(x^{\prime},m)\), which also takes into account the variability around the counterfactual. The results of the ablation study, for various \(\tau\) thresholds, are summarized in Table 7 in Appendix D.

\begin{table}
\begin{tabular}{l l c c c c c c c c} \hline \hline  & & \multicolumn{4}{c}{\(l_{1}\)_based_} & \multicolumn{4}{c}{\(l_{2}\)_based_} \\ \cline{3-10}  & Method & COST & LOF & WI VAL & LO VAL & COST & LOF & WI VAL & LO VAL \\ \hline \multirow{4}{*}{\begin{tabular}{} \end{tabular} } & min Cost & 0.40 & 0.49 & 38.8\% & 35.2\% & 0.11 & 0.75 & 13.5\% & 13.5\% \\  & min Cost+T-Rex:I (Ours) & 1.02 & 0.38 & 98.2\% & 98.1\% & 0.29 & 0.68 & 98.5\% & 98.2\% \\  & min Cost+SNS & 1.20 & 0.30 & 98.0\% & 97.8\% & 0.31 & 0.64 & 97.9\% & 97.0\% \\  & ROAR & 1.69 & 0.41 & 92.6\% & 91.2\% & 1.91 & 0.43 & 86.3 \% & 84.8\% \\ \cline{2-10}  & NN & 1.91 & 0.80 & 51.1\% & 50.3\% & 0.56 & 0.80 & 51.1\% & 50.3\% \\  & T-Rex:NN (Ours) & 2.50 & 0.92 & 84.0\% & 84.0\% & 0.77 & 0.92 & 84.0\% & 84.0\% \\ \hline \multirow{4}{*}{\begin{tabular}{} \end{tabular} } & min Cost & 1.42 & 0.77 & 58.8\% & 56.7\% & 0.48 & 0.81 & 26.6\% & 26.6\% \\  & min Cost+T-Rex:I (Ours) & 4.81 & 0.72 & 98.0\% & 96.5\% & 1.20 & 0.75 & 99.2\% & 98.7\% \\  & min Cost+SNS & 5.71 & 0.67 & 97.5\% & 98.1\% & 1.44 & 0.68 & 99.9\% & 98.9\% \\  & ROAR & 7.63 & 0.54 & 96.3\% & 92.3\% & 6.81 & 0.58 & 87.8\% & 85.2\% \\ \cline{2-10}  & NN & 7.05 & 1.00 & 95.3\% & 95.4\% & 2.50 & 1.00 & 95.3\% & 95.3\% \\  & T-Rex:NN (Ours) & 10.13 & 1.00 & 100\% & 100\% & 3.04 & 1.00 & 100\% & 100\% \\ \hline \multirow{4}{*}{
\begin{tabular}{} \end{tabular} } & min Cost & 0.21 & 0.94 & 74.6\% & 70.2\% & 0.08 & 1.00 & 19.7\% & 14.1\% \\  & min Cost+T-Rex:I (Ours) & 1.11 & 0.83 & 100\% & 98.8\% & 0.42 & 0.94 & 100\% & 99.7\% \\  & min Cost+SNS & 3.34 & -1.00 & 100\% & 98.2\% & 1.07 & -1.00 & 100\% & 99.3\% \\  & ROAR & 3.68 & 0.64 & 98.7\% & 96.4\% & 1.35 & 0.59 & 98.9\% & 97.2\% \\ \cline{2-10}  & NN & 0.39 & 1.00 & 70.5\% & 67.5\% & 0.15 & 1.00 & 70.5\% & 67.5\% \\  & T-Rex:NN (Ours) & 2.22 & -0.33 & 100\% & 100\% & 1.00 & -0.33 & 100\% & 100\% \\ \hline \hline \end{tabular}
\end{table}
Table 1: Experimental results.
## 6 Discussion

We introduce an abstraction called naturally-occurring model change and propose a measure, Stability, to quantify the robustness of counterfactuals with probabilistic guarantees. We show that counterfactuals with high Stability will remain valid after potential model changes with high probability. We investigate various techniques for incorporating stability in generating robust counterfactuals and introduce the T-Rex:I and T-Rex:NN algorithms. We also make a novel conceptual connection with the body of work on model multiplicity, further emphasizing on the models that are more likely to occur.

**Limitations and Broader Impact:** The naturally-occurring model changes rest on assumptions that may not apply to all models or datasets. Our relaxed stability, although practically implementable, lacks the same theoretical guarantees as the initial stability measure. Estimating the Lipschitz constant around a counterfactual can be computationally demanding, particularly when leveraging gradient descent to optimize stability. Though generating robust counterfactuals is a key step towards trustworthy AI, it can fall short of other important factors such as fairness (Sharma et al., 2019; Gupta et al., 2019; Ley et al., 2022; Raman et al., 2023; Ehyaei et al., 2023). Future research could explore links between robustness and fairness, improving the estimation of stability, or integrating Stability into training-time-based approaches for generating robust counterfactuals.

DisclaimerThis paper was prepared for informational purposes in part by the Artificial Intelligence Research group of JPMorgan Chase & Co. and its affiliates ("JP Morgan"), and is not a product of the Research Department of JP Morgan. JP Morgan makes no representation and warranty whatsoever and disclaims all liability, for the completeness, accuracy, or reliability of the information contained herein. This document is not intended as investment research or investment advice, or a recommendation, offer or solicitation for the purchase or sale of any security, financial instrument, financial product or service, or to be used in any way for evaluating the merits of participating in any transaction, and shall not constitute a solicitation under any jurisdiction or to any person, if such solicitation under such jurisdiction or to such person would be unlawful.
