ElliCE: Efficient and Provably Robust Algorithmic Recourse via the Rashomon Sets

Thank you so much for your review and feedback.

W1: The motivation for robustness under the Rashomon effect is weak.

The robustness under the Rashomon effect is a critical and practical motivation for our work. Most counterfactual explanation methods operate on the assumption that the underlying predictive model is fixed. They generate explanations, often very close to the decision boundary, that are valid for one specific model. In practice, this assumption rarely holds. In many real-world applications, especially on tabular data, numerous distinct models can achieve similar, near-optimal performance. This is known as the Rashomon effect, or model multiplicity. Institutions may switch between these equally good models as part of regular maintenance, retraining cycles, or to adapt to minor shifts in data. This leads to a significant reliability problem. A counterfactual action guaranteed to work for one model may fail completely if another one is deployed. This undermines the entire purpose of algorithmic recourse, which is to provide reliable and actionable paths to a desired outcome.

Our work, ElliCE, directly confronts this failure point. Instead of generating an explanation for a single model, we produce recourse that is provably valid across an entire ellipsoidal approximation of the Rashomon set. This ensures the explanation given to individuals is a robust and reliable action, closing a gap in providing trustworthy algorithmic recourse under real-world model uncertainty.

Small comments

Thank you for such a detailed list of typos and suggestions. We will incorporate them, including fixing typos, notation, increasing fonts, improving margins, and so on. Next, we focus on clarifying questions and details.

In line 98, the parameter vector is defined as an element of the (strict) subset of R^d, but if the hypothesis space is a multi-layer perceptron, this is not correct. Also, why is this Theta a strict subset?

The notation is a general formalization. For MLP, d represents the total number of weights and biases across all layers. The strictness is not necessary. It is a common formality used for the convenience of theoretical results.

Line 116, R(epsilon) has a hat, while it was defined without one in line 109. Also line 296, or is this the approximated Rashomon set defined in line 278-289?

We used a hat to show that the Rashomon set is empirical, as constructed on the empirical data and not on the distribution. We will improve the consistency of the notation.

Line 148, for MLPs, is the CE given in terms of the original feature space, or of the output values of the last MLP layer?

The counterfactual explanations are given in the original feature space. While we construct the ellipsoidal approximation of the Rashomon set in the embedding space produced by the penultimate layer of the neural network, we map this information back to the original feature space. Our algorithms (both gradient-based and cutting plane methods) operate directly in the input space while using the worst-case model identification from the embedding space to guide the search for robust counterfactuals. This approach allows us to provide explanations in terms of the original features that users can understand and act upon.

Line 154, is this approximation an over or under approximation?

Using the Hessian matrix to define the set generally results in an under-approximation of the Rashomon set. The ellipsoidal set can be understood as a tractable geometric approximation based on the local curvature of the loss function.

Please note that for L2 regularized mean-squared error, the Rashomon set can be computed exactly as an ellipsoid, meaning it is not an approximation in that scenario.

Limitations and Line 159, 160, how likely is the assumption that the Hessian matrix H is positive definite?

The positive definiteness of the Hessian matrix is a reasonable assumption in our context. For linear models with convex losses, the Hessian is positive definite by definition. For neural networks, where the loss landscape may be non-convex, we ensure numerical stability by adding a small regularization term to the computed Hessian (equivalent to L2 regularization with a parameter). This practice is standard in optimization literature and guarantees that our ellipsoidal approximation remains well-defined. Our experiments confirm that this approach works effectively across diverse datasets and different model classes. Therefore, the positive definite Hessian assumption is a standard condition in the context of regularized models and is satisfied within our paper’s framework.

Line 198 states that the data-supported CE generation is used in section 6. Is the non-data supported CE generation also empirically tested? Or if not, why is it included in the text?

Yes, we provide the results for non-data supported counterfactual generation in Appendix C.7. We couldn’t include these results in the main paper due to space constraints, so we moved a lot of empirical findings to Appendix.

Line 251-254 mentions the trade-off between robustness and proximity of the CE. How, in practice, is the right value for epsilon chosen?

As we state in Section 6, all multiplicity parameters ($\epsilon$ for ElliCE, $\delta$ for Delta Robustness, or $\tau$ for T:Rex) for each baseline are tuned via grid search on the validation set, basically maximizing validity. Given the robustness-proximity trade-off, this approach ensures that the counterfactual explanations are long enough to ensure the basic requirement, which is validity, but still satisfy the proximity requirement. We use grid search to provide each method the capacity to be robust, as the parameters of the baselines are not interpretable, and the papers of the baseline method often do not provide details on how to choose them. We provide guidelines on selecting a suitable epsilon in Appendix C4

Why are two different target epsilon values used in Fig. 2? How are these values chosen?

The target epsilon defines the allowable tolerance in empirical loss used to construct the evaluator baselines (Random Retrain, AWP, and Rashomon Dropout). For the analyses in Figure 2, the values of 0.02 and 0.04 were chosen arbitrarily from within a predefined range of [0, 0.1] as representative examples of possible loss increases due to model multiplicity.

Line 308-309, why is validity maximized in the grid search? Why not robustness or proximity or a trade-off?

We maximize validity during grid search because it's the fundamental prerequisite for counterfactual explanations. Optimizing directly for robustness would create circular reasoning since our evaluators depend on the same parameters we're tuning. Validity provides a fair comparison baseline across methods, after which the inherent robustness-proximity trade-offs of each method can be meaningfully compared.

Line 308-309, how can validity be used in the grid search for ElliCE if by line 206-209 validity is guaranteed, and line 322-323 also confirm 100% validity? Fig. 1, why is ElliCE’s robustness not always 100%?

For a given data point, a valid counterfactual might not exist if the robustness constraint ($\epsilon$) is too large, making the problem infeasible. The grid search tunes epsilon to maximize the number of data points for which a valid counterfactual can be found. ElliCE's robustness is not always 100% because it is guaranteed for its internal ellipsoidal approximation, but it is evaluated against different evaluators (e.g., Random Retrain, AWP, Dripout) . Since these approximation sets are not identical, a counterfactual proven robust for one set may not be valid for all models in the other.

Line 320-321, why does ElliCE operate on the features of the penultimate layer of a neural network? What guarantees does this provide?

ElliCE operates on the penultimate layer to treat the neural network's final layer as a simple linear model. This makes the problem tractable and allows the same convex optimization framework to be applied directly.

How are the datasets presented in Table 1 selected? And why are robustness results in Table 4 much lower for the baselines?

Table 1 or the main text is a subtable of Table 4 from the Appendix. The choice of datasets is random. The robustness results for the baselines are lower because ElliCE performs better.

Thank you very much for your review and feedback.

W1 and Q1: Can you comment on the plausibility performance of the proposed method?

Thank you for asking about plausibility. Our data-supported counterfactuals are inherently plausible, as they lie on the data manifold. As you mentioned, for non-data supported counterfactuals, there are possible cases when the resulting counterfactuals for continuous data may be in low-density (not plausible) areas. Very often, it happens when the resulting counterfactual is too close to the decision boundary. Note that due to our robustness requirements, we push counterfactuals away from the decision boundary. We can illustrate this on the synthetic examples, based on two Gaussian distributions, where each class $P(y|x)$ comes from a separate Gaussian. Visually, we observe that the resulting ElliCE counterfactuals are pushed towards the mean of the Gaussian of the target class. Visually, they also look in good density regions, where other points are also present. We will add this illustration on the synthetic data to the paper. Also, note that we observe a robustness-proximity trade-off, as we report in Theorem 5.

Counterfactuals might not be plausible when we push too far away from the decision boundary. In case of ElliCE, its hyperparameter epsilon allows us to control for this, as we seek for closest examples for some robustness level.

To provide more rigorous analysis, we evaluated the plausibility of the counterfactuals using LOF metric (as in [1]). LOF is used for outlier detection, and values close to 1 are preferred. Please see our results in the table below. For data-supported counterfactuals, we have small values comparable to baselines. For non-data supported counterfactuals, ElliCE performs slightly better than TREX and better than ROAR on the diabetes dataset. PROPLACE, which is designed for plausibility, achieves better results on LOF. However, please also note that ElliCE achieved better robustness results as compared to PROPLACE. We will add a more comprehensive study of the plausibility to the paper, if accepted.

Table: Plausibility study of Ellice and baselines under different evaluation strategies.

To further help to control counterfactuals' actionability and plausibility, we did a series of updates to our algorithm. Besides having immutable features, we now allow:

1. Change of the feature in one direction (as if the value can only increase but not decrease, say if a person plans to increase income or spending amount when applying for a credit). This is implemented by blocking gradient components that would move features in the prohibited direction. For the data-supported method we do prefiltering.

2. Fixed range for the feature. This is achieved by restricting gradient flow to stay within bounds and pruning KD-tree searches to exclude out-of-range solutions.

3. sparsity of the counterfactual explanations. We additionally employ L0 and L1 regularizations to allow for more sparse explanations for the cases when the user wants to change a small number of features. For data-supported counterfactuals, this is implemented by defining a custom distance metric for Ball tree search that combines sparsity and proximity ($metric = C \times \ell_0 + \ell_1$, where a coefficient C prioritizes minimizing the number of changed features). We hope that these changes will help the users to find more plausible and actionable counterfactuals and overall contribute to the flexibility of our framework.

[1] Jiang, Junqi, et al. "Provably robust and plausible counterfactual explanations for neural networks via robust optimisation." Asian Conference on Machine Learning. PMLR, 2024.

W2: The method applies only to the linear classifier case (line 98 ) and to binary classification.

ElliCE works beyond linear models and works for non-linear hypothesis spaces that can be optimized with gradient descent. In our experimental section, we included results for the multi-layer perceptrons and show that our algorithm effectively works for this case. We describe in the Appendix in Algorithms 1 and 2 the method that we use to compute counterfactuals for the neural network space. Our theoretical section applies to the penultimate layer of the MLP. In this layer, we do the Rashomon set approximation and use our theory in Theorem 1 to compute the worst-case model for the Rashomon set. Please note that a lot of counterfactual generation methods known to us focus on binary classification (except PROPLACE), as it is very popular in tabular data settings. However, it is easy to extend ElliCE to the multiclass setting. We describe this extension below.

Q2: Does the proposed method also apply beyond the binary classification setting?

There is a simple way to extend our framework to handle one-versus-one multiclass classification. For one versus one, we will assume that we know the target class for the counterfactual. Our task succeeded if, for the resulting counterfactual, the target class’s $y$ logit is higher than each logit of some other class $c$ individually, under the worst-case model within the Rashomon set.

Formally, for a counterfactual point $x_c$:

$$\min_{\hat{\mathcal{R}}(\epsilon)} (f_{\theta}^{y} (x_c) - f_{\theta}^{c} (x_c)) \ge \tau, \qquad \forall \space c \neq y$$

where $\hat{\mathcal{R}}(\epsilon)$ is our ellipsoidal approximation of the Rashomon set, and $\tau$ is a margin hyperparameter that we can set to 0. This guarantees that the explanation remains valid across all models within the Rashomon set.

This way, we get $K−1$ second-order constraints (one for each rival class) and can use the same closed-form inner optimization as in the binary case. The overall optimization problem remains a convex second-order cone program (SOCP). The runtime grows roughly linearly in $K$.

For neural networks, we apply the same approach using a Rashomon set approximation in the embedding space and a gradient descent or cutting-plane approach, as in our binary case.

Q3: It appears in the code that baseline implementations are taken from the ROBUSTX library, which is not mentioned in the paper.

Thank you for pointing this out. You are correct, we used baseline implementations from the ROBUSTX library [2]. We will update the manuscript to include a proper citation and clarify that baselines implementation was done using ROBUSTX.

[2] Jiang, Junqi, et al. "RobustX: Robust Counterfactual Explanations Made Easy." arXiv preprint arXiv:2502.13751 (2025).

W3: The optimization approach in line 167 and Equation 1 is standard.

We actually think it is a good thing that the definition and high-level optimization problem for the counterfactual explanation is standardized in the literature and community. Note that the optimization problems that you mentioned are quite general: minimize the distance with respect to some constraints. And here is where ElliCE shines, as we provide new way of handling constraints: using the Rashomon set that we approximate with an ellipse. Our experiments show that this approximation of the Rashomon set is effective when we use multiplicity evaluation baselines such as AWP and random retrain.

Thank you very much for the review and feedback.

Quality W1. The ellipsoidal approximation of the Rashomon set is a simplification that may fail in real-world scenarios

The ellipsoidal approximation is a simplification of the Rashomon set that can be hard to find and fully approximate for a lot of datasets. However, our experiments are not limited to convex geometry only and effectively demonstrate that the approximation works and delivers robust results. More specifically, in Figure 1b, Table 1 and in Appendix C Figure 6 (Figure 11 for non data-supported) we use neural networks with cross-entropy loss and do not expect the loss to be convex. Therefore, we provide empirical evaluation in non-convex scenarios.

Q2 and Quality W2. Experiments focus on linear models and neural networks with linearized final layer, leaving scalability to deep architectures untested.

Our methodology is designed to apply to deep architectures by focusing its theoretical guarantees on the network's final linear layer. This layer operates on the embedding created by the preceding deep layers. This allows us to keep the robustness calculation tractable and fast.

Once a worst-case model is computed in the embedding space using Theorem 1, the framework uses an optimization technique, such as the cutting-plane or gradient descent (we support both), to find the corresponding counterfactual in the original input space. This two-stage process allows the method to scale to deep networks, as our experiments confirm. We will explain this in more detail in the paper.

Clarity W1: The Background and Notation section would benefit from better structure.

The background and notation section is only one page long, however, we will apply a modular structure and emphasize changes in the definitions better. We will add equation numbers to improve the paper's clarity.

Clarity W2: The transition from theoretical analysis to empirical validation in Section 6 is abrupt.

Thank you, we will add better smoothing text when transitioning from Section 5 to 6.

Clarity W3: Technical intuition behind the Hessian-based approximation of the Rashomon set and how it captures model uncertainty is underdeveloped, leaving non-experts struggling to grasp the geometric reasoning.

The core idea is to view the set of good models (the Rashomon set) as a "valley" in the loss landscape. Our method uses the Hessian matrix to mathematically describe the curvature or shape of this valley at its lowest point. A flat valley signifies high model uncertainty (many different models are good), while a steep valley indicates low uncertainty. The ellipsoid is a direct, geometric map of this curvature, providing a tractable representation of the model uncertainty that we use to find a provably robust counterfactual. We will add this more detailed explanation to the paper.

Please note that Hessian-based approximations are a common technique used in various settings (see [1]-[3]).

[1] Dinh, Laurent, et al. "Sharp minima can generalize for deep nets." International Conference on Machine Learning. PMLR, 2017.

[2] Zhong, Chudi, et al. "Exploring and interacting with the set of good sparse generalized additive models." Advances in neural information processing systems 36 (2023): 56673-56699.

[3] Sagun, Levent, et al. "Empirical analysis of the hessian of over-parametrized neural networks." arXiv preprint arXiv:1706.04454 (2017).

Limitation2 and Significance W1: While the paper emphasizes technical robustness, it lacks discussion on real-world deployment challenges, such as ensuring that suggested feature changes are actionable for users.

Thank you for mentioning this. For the specific case when the users can not change some feature or they are immutable, such as for example, date of birth, we allow the users to select a set of features that they do not want to change (please see Appendix C.5.1 for more details). Technically, to obtain such counterfactuals, we block the gradient when optimizing with gradient descent for non-data-supported counterfactuals and search for counterfactuals that do not change immutable features for the data-supported version.

Because we are focused on providing a robust and actionable recourse, we also introduced a series of changes in order to make counterfactuals more actionable.

More specifically, we now additionally allow users to

1. Change of the feature in one direction (as if the value can only increase but not decrease, say if a person plans to increase income or spending amount when applying for a credit).

2. Fixed range for the feature.

3. Model for sparsity of the counterfactual explanations. We additionally employ L0 and L1 regularizations to allow for more sparse explanations for the cases when the user wants to change a small number of features.

We believe that these changes should significantly improve the actionability of our explanations as well as user-experience while working with ElliCE.

Limitation2 and Significance W2: The societal implications of counterfactuals, particularly in cases where they might reinforce biases or suggest unethical modifications, are not addressed.

While our work focuses on the technical challenge of providing reliable recourse, the resulting counterfactuals should be carefully audited in deployment to ensure they do not reinforce systemic biases. We will add this to our paper discussion and conclusions.

Originality W1: While the framework is original in its formulation, some ideas like parameter perturbation robustness [1], are inspired by prior work.

Thank you for your comment. We carefully studied $Delta$-robustness in our paper both empirically and theoretically. Note that $Delta$-robustness works for neural networks only and introduces the perturbation to the whole neural network, while ElliCE uses ellipsoidal approximation in the embedding space for neural networks. In Appendix D (Proposition 1 and Theorem 6) , we theoretically show the limitations of the approach that $Delta$-robustness takes and show that it has limited ability to capture meaningful model multiplicity. Empirically, Ellice achieves better robustness guarantees, please see Figures 4, 6, 10, 11 in Appendix C.

Originality W2: How the method advances over existing methods that rely on ensembles or MILP [2], rather than geometry-based approximations.

We do not rely on ensembles or MIP. Instead, we propose a new method that relies on ellipsoidal approximation for the model multiplicity, gradient optimization for non-data supported counterfactuals and KD-tree and Ball-tree search for data-supported counterfactuals. Our method is much faster (Table 2) and more robust (Figures 4, 6, 10, 11 and Table 4 in Appendix C).

Limitation1 and Q1. How does your method perform when the true Rashomon set has non-convex or disjoint regions like multiple local optima?

While the ellipsoidal approximation simplifies complex Rashomon sets, our empirical results show it remains effective even for non-convex cases. Our experiments with neural networks using cross-entropy loss (Figure 1b, Table 1, Appendix C Figures 6 and 11) demonstrate that ElliCE delivers robust counterfactuals despite non-convex loss landscapes.

In cases with disjoint regions of good models, the ellipsoidal approximation doesn't fail but provides a local approximation around one optimum . In the worst case, this locality can lead to more conservative counterfactuals (due to the "price of robustness") or overly optimistic ones (which are not robust to the true, larger Rashomon set) . To account for multiple local minima, one could consider an ensemble of models trained with different random initializations. Note that we compared ElliCE against such a "Random Retrain" evaluation baseline and found ElliCE’s performance to be consistently strong, demonstrating its practical effectiveness even in the mentioned case.

Q3: Section 4 mentions "feature-level constraints". How are these constraints integrated into the optimization?

In our framework, we allow for mutable and immutable features to handle feature-level constraints. More specifically, the user can tell the algorithms to fix some features and not change their values. We described this in Appendix C7. Implementation-based, we fix gradient flow for the gradient-based method and select counterfactuals that do not change the feature for tree-search optimization.

After submission, we further improved our algorithms, and now we also handle:

1. change of the feature in one direction. This is implemented by blocking gradient components that would move features in the prohibited direction. For the data-supported method, we do prefiltering.

2. fixed range for the feature. This is achieved by restricting gradient flow to stay within bounds and pruning KD-tree searches to exclude out-of-range solutions.

3. sparsity of the counterfactual explanations. For data-supported counterfactuals, this is implemented by defining a custom distance metric for Ball tree search that combines sparsity and proximity ($metric = C \times \ell_0 + \ell_1$, where a coefficient C prioritizes minimizing the number of changed features).

Q4: What happens when a positive definite Hessian assumption fails (e.g. for for non-convex losses or high-dimensional data).

The positive definite Hessian assumption is satisfied by design through regularization. If the Hessian is not positive definite, for instance, due to non-convexity or numerical issues, we perform additional regularization by adding a small multiple of the identity matrix, $\alpha*I$, to the computed Hessian. This standard technique guarantees the matrix becomes positive definite, ensuring numerical stability. This approach is effective regardless of data dimensionality or the non-convexity of the loss. We will make these details more explicit in the manuscript.

Formatting issues.

Thank you for pointing out the formatting issues. We will fix them.

W1: Multiplicity computing approach

The ellipsoid-based Rashomon sets will contain models that occur due to slight changes in the training data that T:Rex, Argumentative Ensembling are considering, especially for convex losses, such as logistic loss. This is because this Rashomon set contains different coefficient vectors (and thus feature weighting), allowing for diversity of the models. Therefore, if the data shift occurred, there are other models to choose from the ellipsoidal Rashomon set to account for the changes. To support this, we run a simple experiment to show the stability of the ellipsoid-based Rashomon set (see table below). More specifically, we add noise to the data that comes from a Gaussian distribution with 0 mean and variance $\sigma$ that ranges from 0.1 to 2.

We consider the hypothesis space of linear models and the German Credit dataset. First, we compute the Rashomon set on the original data (we use 10% of the optimal loss as epsilon and standardize the data). For each variance value, for 100 times, we do the following: 1) sample noise, 2) create a new dataset with added noise, 3) train logistic regression on this new data, 4) check if the Rashomon set contains models with similar performance as this optimal model. Finally, we compute the percentage of cases when the Rashomon set contained a similar performing model on this new data as the optimal one in the table below. We find that the Rashomon set is stable, even for higher levels of noise. Therefore, we expect the ellipsoid to contain models that perform similarly to optimal models on data folds (if not directly containing these optimal models).

Table. Fraction of cases when the Rashomon set trained on original data contains optimal models trained on noise data, with noise $\sim \mathcal{N}(0, \sigma)$.

Q1. How does the proposed method compare with the baselines w.r.t. the robustness to changes in training data or re-training on initially unseen data?

To further support that the Rashomon set handles re-training scenarios, we followed a similar setup that T:Rex and Argumentative Ensembling are using. More specifically, we split the data in 2 parts. We use first data part to train our Rashomon set (ElliCE) and train baselines (STCE, RNCE). Second data part we use to train evaluators (random retraining, adversarial weight perturbation, and Rashomon dropout). This, similar to T:Rex and Argumentative Ensembling, considers different data for the methods we evaluate and baselines with which we do the evaluation. For this new evaluation, we still found that ElliCE outperforms or performs on par with the baselines. Since we can not report images, please see the tables below. Here we use Target Epsilon 0.02.

Table: Re-train evaluation of ElliCE for the hypothesis space of linear models

Table: Re-train evaluation of ElliCE for the hypothesis space of multi-layer perceptrons.

W2 and Q2: The applicability of the results to the non-linear case.

Thank you for pointing this out. “Without loss of generality” in the paragraph on line 146 applies to the second sentence, as our results apply to the embedding space of the neural networks (penultimate layer). In this layer, we do the Rashomon set approximation and use our theory from Theorem 1 to compute the worst-case model for our Rashomon set. Then we can use 2 methods to find counterfactuals in the feature space:

1. gradient-descent method that is described in Appendix C7. On each step of the algorithm, we compute the weighted sum of the worst base model (based on the formula in Theorem and the optimal model, then do the gradient step

2. cutting-place method that is described on line 187.

We found the cutting plane method to work very well and better than baselines for smaller Rashomon set epsilons. This method tends to converge faster than gradient descent for smaller epsilons as well. For larger epsilon, we found that the gradient-based method is more stable.

Both of these methods work for neural networks and allow for robust counterfactuals as our results suggest in the experimental section and appendix. For both, we use our theoretical results from Section 4. So no, our theory does not apply only to linear models.

We will explain the methods and how we use our theoretical results more clearly in the paper. We will also remove “Without loss of generality” at the start of line 146 to not confuse the readers, as it refers to the second sentence. Our new paragraph will say: “We focus our theoretical analysis on linear predictors of the form $f_{\theta}(x)=\theta^{\top}x$. However, the methodology presented here is also applicable to the embedding space of the multilayer perceptrons.”

We will also fix minor formatting concerns. Thank you again for your review.