Joint Distribution–Informed Shapley Values for Sparse Counterfactual Explanations
Authors: Lei You
Yijun Bian
Lele Cao
Venue: ICLR 2026 Poster (2026)
Keywords: Counterfactual Explanations
Shapley Values
Optimization
Explainable Machine Learning
Primary area: interpretability and explainable AI

Abstract: Counterfactual explanations (CE) aim to reveal how small input changes flip a model’s prediction, yet many methods modify more features than necessary, reducing clarity and actionability. We introduce COLA, a model- and generator-agnostic post-hoc framework that refines any given CE by computing a coupling via optimal transport (OT) between factual and counterfactual sets and using it to drive a Shapley-based attribution p-SHAP that selects a minimal set of edits while preserving the target effect. Theoretically, OT minimizes an upper bound on the $W_1$ divergence between factual and counterfactual outcomes and that, under mild conditions, refined counterfactuals are guaranteed not to move farther from the factuals than the originals. Empirically, across four datasets, twelve models, and five CE generators, COLA achieves the same target effects with only 26–45% of the original feature edits. On a small-scale benchmark, COLA shows near-optimality.

Published as a conference paper at ICLR 2026
JOINT DISTRIBUTION-INFORME SHAPLEY VALUES FOR SPARSE COUNTERFACTUAL EXPLANATIONS
Lei You"
Yijun Bian University of Copenhagen DK-2100 Copenhagen, Denmark yibi@ ku. dk
Technical University of Denmark DK-2750 Ballerup, Denmark le leiyoddtu.dk
Lele Cao
AI Labs, King (Microsoft Gaming) SE-11157 Stockholm, Sweden
l eleeacole cecicemmicomm wwwsoft.cccooomoooooo. com
ABSTRACT
Counterfactual explanations (CE) aim to reveal how small input changes flip a model's prediction, yet many methods modify more features than necessary, reducing clarity and actionability. We introduce COLA, model- and generator- agnostic post-hoc framework that refines any given CE by computing a coupling via optimal transport (OT) between factual and counterfactual sets and using it to drive a Shapley-
attribution (p-SHAP) that selects a minimal set of edits while preserving the target effect. Theoretically, OT minimizes an upper bound on the W1 divergence between factual and counterfactual outcomes and that, under mild conditions, refined counterfactuals are guaranteed not to move farther from the factuals than the originals. Empirically, across four datasets, twelve models, and five CE generators, COLA achieves the same target effects with only 26-45% of the original feature edits. On small-scale benchmark, COLA shows near-optimality.
Experiment code: https: / http//tthhuucccc http:///ooooaaoo/o//////o/////o/ Software (Zhu & You, 2026): https: http://ggg/poooott//alollollll
BACKGROUND
Explainable Artificial Intelligence (XAI) is essential for making artificial intelligence systems trans- parent and trustworthy (Arrieta et al., 2020; Das & Rad, 2020). Within this area, feature attributions (FA) methods, such as Shapley values (Sundararajan & Najmi, 2020; Lundberg & Lee, 2017), deter- mine how much each input feature contributes to a machine learning (ML) model's output. This helps simplify complex models by highlighting the most influential features. For example, in healthcare model, Shapley values can identify key factors like age and medical history, assisting clinicians in understanding the model's 's decisions (Ter-Minassian al., 2023; Nohara et al., 2022). Another technique counterfactual explanations (CE) (Wachter et al., 2017; Guidotti, 2022) show how small changes in input features can lead to different outcomes. While hundreds of CE algorithms have been proposed (Guidotti, 2022; Verma et al., 2020) to date, it is hardly practical to find one single CE algorithm that suits for all user cases, due to each of them is tailored particularly for their own different scenarios, goals, and tasks. For instance, the objective in one CE algorithm can be defined as finding a single counterfactual instance for each factual instance sometimes, while at othertimes, it could be treating multiple instances as group and seeking or more/multiple counterfactual instances for each factual observation. In some cases, the focus of a CE algorithm could be on the entire dataset, aiming identify global CE that indicate the direction to move the factual instances to achieve the desired model output (Rawal& & Lakkaraju, 2020: Ley et al., 2022; 2023; Carrizosa et al., 2024). Yet in other scenarios, the group of factual instances is viewed distribution, aiming to find a counterfactual distribution that remains similar in shape to the factual distribution (You et al.,
X Corresponding Author.

---

Published as a conference paper at ICLR 2026
2025), and ensuring comparable costs. Besides, some CE algorithms assume differentiable models, whereas others are designed specifically for tree-based or ensemble models.
Problem Description and Challenges We address a general and comprehensive problem this paper, which builds on the extensive foundations established in the literature (see Appendix A). To be specific, we seek to answer the following question: Although both and CE are vital for making Al models more interpretable and accountable, relying only on one of them will cause drawbacks. That is, FA alone may not provide tangible steps as explanations, and CE alone might include unnecessary feature changes that not practical. Therefore, address a general and comprehensive problem in this paper, which builds on the extensive foundations established in the literature (see Appendix A). To be specific, seek to answer the following question:
Given (group of) factual instance(s), how can we devise action plan that re- quires the least feature modifications to achieve a desired counterfactual outcome?
z'
z"
yz
R
8
R
200
5
No
250
Yes
200
No
Yes
150
3
No
150
No
150
Yes
Yes
100
2
No
350
Yes
100
No
Yes
150
6
No
150
No
350
Yes
Yes
Figure 1: [Example: User engagement on an e-commerce platform] A platform aims to increase user registrations. The platform has collected data user interactions, such the amount of money spent (LO),, the number of clicks ( ). and whether the user has registered R In the original data (x),
registered. Action plans z' and z" adjust user characteristics achieve the desired outcome (y*) ) of full registration. Both plans achieve a half counterfactual effect, but z" requires fewer modifications compared z'.
Three major challenges remain in addressing this problem. First, it is unrealistic to expect single CE algorithm to meet all the needs universally, as the problem is often task-specific. Second, the approach should not rely strong assumptions about the model (for example, requiring differentiability special structures) to ensure its applicability across wide range of models. Third, FA like feature importance can be misleading due to the lack of coherence between the FA scores and the changes for counterfactual effect. In other words, it is not effective to perform FA independently of CE to select the most important features to change. We will demonstrate later that this decoupling can result in counterproductive feature modifications (also referred to as actions from users/stakeholddrs) in Result II of Section 6, as the features deemed important generally may not align with the specific pathways to achieve the desired counterfactual outcomes.
Main Contributions We introduce COunterfactuals with Limited Actions (COLA), a general post-hoc framework that refines CE across models and CE generators by using an optimal transport ind)-idduced coupling between factual and counterfactual sets to guide Shapley attributions. This yields p-SHAP. which unifies other commonly used Shapley methods under appropriate couplings and offers modular interface for attribution and edit selection. We provide theoretical guarantees: OT minimizes an upper bound on the W2 divergence, and under mild conditions refined counterfactuals provably farther from factuals than the originals. Empirically, across four datasets, twelve models, and five CE methods, COLA with p-SHAP achieves the same target effects with only 26-45% of the original feature edits. In small-scale benchmark, COLA is near-optimal.
To our best knowledge, this is the first method proposed for systematically addressing the problem in Figure without specifying certain CE algorithms and models.
PROBLEM FORMULATION
We formally formulate the problem described in Figure 1. Denote R: - R as black-box ML model. Denote by x any observed (factual) data with n rows and d columns (x E xdd xd n 2 1, and d 2 1). Let y" be the target model output (y" E R m 2 1). The optimization is to look for
2

---

Published as a conference paper at ICLR 2026
(group of) counterfactual data instance(s) (z R=xd) subject to maximum number of allowed feature changes, C, to achieve model output(s) as close as possible to y". Let D denote a divergence function that measures the dissimilarity between two entities. The problem is below.
min 0,8
y(
(1a)
s.t.
(2,,x) Se,
(1b)
E k= Cik S C 2sh S Ta (1 ((c + Mea i == 1.......1. n1, 4: 1, d
(1c)
(1d)
2sh x. (1- ca) Me i i= 1, ".....,.. n, h 1, ,,...,d d.
(le)
The objective equation (Ia) and the constraint equation (Ib) formulate the typical CE optimization. Namely, 2 is expected to make f(z close to y* yet stays close to x. Then z can be used as a counterpart reference explain why f(x does not achieve y". We do not limit the function D to any specific type of divergence function, allowing it to stay general. Example functions of D can be Euclidean distance, OT, maximum mean discrepancy (MMD), or differences of the model outcome in mean or median. Then, equation (Ic)-equatio (1e) compose the CE optimization constrained by actions. On top of the counterfactual data z, we also optimize an indicator variable c, such that 24 is not allowed to change iff ca = O. Maximum C changes are allowed imposed by equation (Ic). Inspecting equation (Id) and equation (Ie), if Cjk O, Fikk equals 2 and no changes happen at (i, k). Otherwise, if Gı remark that M is sufficiently large constant such that 2 has good freedom to change.
To solve equation (1), we resort to FA to identify the most influential features to obtain the modification indicator variable c. The next section introduces commonly used Shapley value methods for FA, which, together with our later proposed one, are integrated into our algorithmic framework COLA, to obtain the refined counterfactual z. The problem is computationally difficult even when d' = 1 for linear models, see Appendix B.
PRELIMINARIES ON SHAPLEY VALUE
This section introduces commonly used Shapley value methods for FA, which, together with our later proposed one, are integrated into our algorithmic framework COLA. We first introduce the concept of Shapley value in game theory, followed by a discussion various commonly used Shapley methods for FA. Readers could find in (Sundararajan & Najmi, 2020) for a comprehensive overview of different Shapley value methods.
In cooperative game theory. a coalitional game is characterized by a finite set of players (in our context, features), denoted by F 21. d), and a characteristic function This function v 2 + R maps each subset S C F to a real number v(S), representing the total payoff that can be achieved by the members of S through cooperation, with the condition v0) O. The Shapley value, as defined based on e later, is a fundamental concept in this framework, providing an equitable method to distribute the overall payoff, v(F), of the grand coalition among the individual players.
Prior to defining the Shapley value, denote A(k, S) the incremental value produced by player (feature) & to a coalition S. aligning with the concept of marginal contribution. This is defined as the additional payoff realized by the coalition due to the inclusion of the player k, i.e.,
A S) = vSS U k) - uSS).
The Shapley value Φ, for a player i in cooperative game is defined in equation (3) below.
Φb = Σ
d (d-11 | (1s)
(kk,
SCF
The value Φe is computed by averaging the normalized marginal contributions A A(k, S) of feature k across all subsets of the feature set excluding k: The normalization factor is the inverse of the binomial coefficient ((1))) where is the total number of features. By summing these values for all subsets and dividing by d, the equation provides a fair and axiomatic distribution of payoffs among players, reflecting each features's contribution among the whole collection of features.
3

---

Published as a conference paper at ICLR 2026
For arbitrary data point x,, x to represent its part containing only the features in S, and F,,s s the other remaining part. The same rule applies for other notations as well. Below, we introduce commonly used Shapley value methods.
Baseline Shapley (B-SHAP). This approach relics on a specific baseline data point r as the reference for any x,. The set function is defined to be
! (S) = (x - f(r));
where feature's absence is modeled using its value in the reference baseline data point (Lundberg & Lee, 2017; Sun & Sundararajan, 2011; Merrick & Taly, 2020). This method assumes an exact alignment between and r.
Random Baseline Shapley (RB-SHAP). This approach is a variant of BShap and is implicitly used in (Lundberg & Lee, 2017) as well as the famous "SHAP" library. The equation is as follows.
RR
( (S) Ex [(x,ss:
[f
(5)
where D is the background distribution that is commonly selected to be the training dataset (Lundberg & Lee, 2017; Merrick & Taly, 2020).
Counterfactual Shapley (CF-SHAP). This is an extension from RB-SHAP, by taking the distribution D as the counterfactual distribution defined on every single factual instance x; (i 1,2,..... n2). (6)
c (S) = Er--((x) (x - E [f(r)]
This method has been adopted in (Albini al., 2022; Kommiya Mothilal al., 2021), demonstrated with advantages for contrastive FA. This method assumes probabilistic alignment between x and r.
4 THE PROPOSED ITS THEORETICAL ASPECTS
SHAPLEY (p-SHAP) AND
(Proposed) p-SHAP We generalize equation (4) equation (6) by integrating an algorithm Aprb that returns their joint probability.' Our p-SHAP is defined as follows,
e (S) = Er-p(rix.) [f(x (F(s) Er-w(r) [f(r)),
(7a)
s.t. p Aprb(xx r)
(7b)
Crucially, unlike prior methods that rely on random baselines, P-SHAP solves the alignment problem by identifying the optimal coupling via OT. This reformulates feature attribution not merely as prediction decomposition, but as a transport problem that minimizes the cost of explanation.
Remark that Apob is independent of the CE algorithm, but only depends on the generated counterfac- tual r. By focusing solely on these fixed components, p-SHAP ensures consistency FA without being influenced by the variability of different CE generation processes, which is major difference to CF-SHAP. Contrary to common expectations, we demonstrate that OT can be more effective than relying counterfactual distribution defined by CE generation mechanism done by Albini et al. (2022), in Result II of Section 6 later.
Especially, of the focus in this paper is to consider the OT problem (also the 2- Wassssssteiin divergence) defined below. And the transportation plan Por obtained by solving OT is used as the joint distribution of and in SHAP.
m ΣΣ
m Pij log
Por arg min
Piyllx |x; -
r2 11+++
+
Pij u4v
(8)
pET ==
Note that and v represent the marginal distributions of x and respectively, and I[ v) the set of joint distributions (i.e. all possible transport plans). The term Σ;_ Σ Pty log(Pt// (Har,) the entropic regularization with E 2 0 being the coefficient. Such regularization (E 0) helps accelerate the computation of OT.
First, p SHAP degrades to B- SHAP in equation (4) when Aper defines a joint distribution between x and r that indicates an i <+ j alignment of for any x,, Second, p-SHAP degrades to EB-SHAP in equation (5) when Apos is defined to be independent of CR but associates with arbitrary distribution D. Third, p-SHAP degrades CF-SHAP in equation (6) when Apos built upon known distribution of CE.
4

---

Published as a conference paper at ICLR 2026
Theoretical Aspects of p-SHAP OT minimizes the total feature modification cost (i.e. modifying x towards r) under its obtained alignment between factual x and counterfactual r. This directly corresponds to objective of finding feature modifications that achieve the counterfactual outcomes at minimal cost. We can further strengthen this connection theoretically under the Lipschitz continuity assumption of the predictive model f. In Theorem 4.1 below (proof in Appendix C), establish that the transportation plan Por used in p-SHAP is effective in minimizing an upper bound on the divergence between fx) and y". Specifically, the I-Wasserstein distance between f(x) and y", is bounded by the Lipschitz constant (assuming Lipschitz continuity of multiplied by the square root of the minimized expected cost of changing x towards r, i.e. Σ,, Pty ||x; ry where Pij = 2. ................................ m) quantify how the feature values of x; should be adjusted towards those of one or multiple Fz. Practically, this means that in p-SHAP, the OT plan por provides strategy to adjust the feature values of x towards those of in a way that minimizes the expected modification cost En, Piy ||; Piy]|x;
ryllă- Compared to other modification plans (p E IT), Por yields the minimal possible cost, which in turn provides the tightest upper bound on the violation of the counterfactual effect Wf(f(x), y*) in proportion this cost.
Theorem 4.1 (p-SHAP Towards Counterfactual Effect). Consider the 1-Wasserstein divergence Wi, ie. Wi(J(x),y") = mill > D5=1 Tej [xxx) - -Y]l- Suppose the counterfactual outcome y" is fully achieved by i.e. y; = f(r,) (j .., ..). Assume that the model -Rd IR is Lipschitz continuous with Lipschitz constant L. The expected absolute difference in model outputs between the factual and counterfactual instances, weighted by Por (with 0), is bounded by:
m
y")
P
r1
<
Piy xx
ryl12,
Vp E II.
Namely, Por minimizes the upper bound of (2(f(x). where the upper bound is based on the expected feature modification cost.
Although Theorem 4.1 bounds the transport cost, it serves as a convex proxy for the A'P-hard discrete sparsity (I-)) problem. By concentrating attribution mass the most efficient paths, the OT coupling effectively guides the greedy selection in Algorithm to discard non-essential features. In addition, P-SHAP is conceptually correct in attributing the causal behavior the modifications of the characteristics, stated in Theorem 4.2 below (proof in Appendix D).
Theorem 4.2 (Interventional Effect of p-SHAP). For any subset S C F and any x; (i = 11,2.... 22.....n) .n). (S) represents the causal effect of the difference between the expected value of fr) under the intervention on features S and the unconditional expected value of f(r). Mathematically, this is expressed as:
[[f(r)] [E()=)) ( (rs x1,s)l
Furthermore, we remark that P-SHAP preserves nice axioms of B-SHAP and RB-SHAP, which makes an effective tool for attributing features. We omit the proof but refer to (Sundararajan Najmi, 2020; Lundberg & Lee, 2017) reference for axioms of Shapley.
THE ALGORITHMIC FRAMEWORK COLA
Sketch The algorithmic framework COLA, stated in Algorithm 1 below, aims to solve equation (1) and is established on four categories of algorithms. We explain Algorithm 1 in details below, along with an illustration in Figure 2.
Line 1 (Applying CE algorithm to find counterfactual r). The CE algorithm ACE takes the model f. the factual x, the target outcome y", and the tolerance as input. The algorithm returns a counterfactual staying close with x, with y fr) and D(r,x). c.
Line 2 (Seeking a joint distribution of x and r). We use an algorithm Apob for this task, which takes x and as input, and outputs matrix representing the joint distribution of all and data points in x and r, respectively. The joint distribution represents an alignment relationship (or matching) between the factual and counterfactual rows, and it in Line 5 to compute
5

---

Published as a conference paper at ICLR 2026
Algorithm 1 COunterfactuals with Limited Actions (COLA)
Require: Model J. factual IRxd target y" E I c, and C
Ensure: Action plan E Rexd and correspondingly refined counterfactual 2 E RRadd
1: Use ACE(S, x, ,,,,)
obtain r E Rmmxd, with y* = f(r) and D(r, x)
2: Use Apro((X, r) to obtain the joint distribution matrix p E R
3: Use Astap r, P) to obtain the shapley values Φ E Raxd for each element of
4: Normalize the element-wise absolute values of Φ. i.e., Pik + / (4p E 12*d).
5: Use Avalue (r. p) to obtain matrix q E R=xd
6: For e E {0, 1 1 Cik f- 0 1......."" :, A l=....,d))
7: Sample C pairs (i, k) according to the probability matrix s, and let Cskc 1 for them.
8: Let z + x (z E [
9: for i 4- 1 to do
10: for k + to do
11:
if Cak 1 then
12:
2ik + qike
13:
end if
14:
end for
end for
16:
return and z
Obtain and compute
Use to compute Shapley Match to Fargmax, p(r, jx,
Get c and z by 40 and
r' f
Ace
0.5
P = 0.2
P12 = 0.4
T11
r12
0.1
[ 0.3
(222 0.1
F11
F12
F21
122
0.3
C 2
F21
T22
F31
F32.
4 Ashup
q
Value
4 (1, (2, 1)}
p(x, r) <- Aprob
112
F11
712
I11
912
P2z
P222
Pal
F32
921
T22
Lines 1-2
Lines
4
Line 5
Lines 6-
Figure 2: [An illustration of COLA] This figure shows how COLA gets c and z for equation (1). We use A mus for illustration in line 5 due to its simplicity. In lines 6-16, we assume C = and the sampling yields exactly two positions for modfications according to the probability matrix sp
the values that can be used for composing later As discussed in Section 3, Apro can be based on OT to compute a joint distribution that yields the smallest OT distance between x and r, if the alignment relationship between their joint distribution that accurately reflects the alignment between the rows in x and r.
unknown. Otherwise, it is recommended select
Lines 3-4 (p-SHAP FA). We apply equation (7) to compute the shapley value for x. The joint distribu- tion p can be used here (without being forced) to properly align each row of x with its corresponding counterfactual rows of r, such that the selected rows in reference for the row in x. Numerically, this alignment significantly influences our contrastive FA. Then, the shapley values (as matrix) of taken element- wisely with the absolute values and normalized (such that all values sum up one). The resulted matrix & forms our FA.
as the most representative contrastive
Line 5 (Computing feature values). The algorithm Avalue is used for this task, which takes the counterfactual r and the joint distribution p as input. For any row i in Avalue selects one or multiple rsw(s) in r, used as references for making changes in x. The algorithm returns a matrix q E IR nux where each element gok serves as counterfactual candidate for aas (Vi, k). Below, we introduce A and A respectively for the cases of selecting single row and selecting multiple rows.
For any X;, Value selects the row of r with the highest probability.
q
max Value (r. p), where qsk TT()). and T(i) arg max Pij
(9)
j=1,2....m

---

Published as a conference paper at ICLR 2026
The algorithm A computes lik convex combination weighted average of Pik, P2k:... Tmmk
172
q A A
p), where qik
Pij 2 Pij'
Tsk
(10)
j=1
Lines 6-16. Recall that the non-negative matrix is normalized to have its summation being one. We could hence treat as policy to select the positions in x for value replacement (i.c. Gik- = 1), as what line does. Then, for any and k with Cue- = 1, Zik gets modified to lik. and the modified matrix then returned as z together with с forming the optimized solutions of the problem in equation (1).
By Theorem 4.1, COLA is designed to minimize the dissimilarity between (z) and y" by modifying z based on feature attribution results, which identify the most important features to adjust to achieve the desired counterfactual effect. The theorem below (proof in Appendix E) demonstrates that the refined 2, produced by the COLA framework, satisfies the constraint equation (1b) in the typical scenario where
m, using the Frobenius supporting the general applicability of this conclusion can be found in Table 3.
as the distance measure. Empirical evidence
Theorem 5.1 (Counterfactual Proximity). Let q E Raxd be the aligned reference matrix derived from the counterfactual set In via the value assignment step q = Ava p), where p is the coupling matrix. For any dimensions m 2 1, the refined counterfactual z constructed by the COLA framework satisfies:
2 x 11---
This inequality guarantees that the refined counterfactual z. is strictly closer (or equal) to the factual input x compared to the aligned counterfactual proposal q, ensuring that the refinement process introduces no additional divergence.
Corollary 5.2. In the special case where = m and the OT plan p corresponds to a deterministic permutation (i.e., obtained without entropic regularization, ε 0), the matrix q becomes row-permuted version of (denoted as where To(s)). ). In this setting, the bound simplifies to:
x S [ro x||F,
recovering the intuition that the refined counterfactual at least as close to the factuals as the original counterfactual set reordered by a.
Complexity of COLA Let O(Mcr) be the algorithm complexity of AcE- For algorithm Ashup- consider using weighted linear regression to estimate Shapley values, and denote by Mshap the number d, and the regularization parameter
of sampled subsets. The complexity of COLA with respect e of entropic OT is O(McE) + Onn log(1/E)) + O(ndMstap) + N where = O(nm) + Ond))
A max is used and = O(nmd) Avalue is used. See Appendix for more details.
6 NUMERICAL RESULTS
This section evaluates the effectiveness of COLA in addressing the problem in equation (1), with ="= f"= f(r) where is the counterfactual obtained from a CE method ACE. We adopt four different divergence functions: OT evaluates the distance between entire distributions. MMD evaluates the divergence between the means of two distributions in a high-dimensional feature space. The absolute mean difference (MeanD) and absolute median difference (MedianD) evaluate the divergence between mean and median, respectively. The numerical results aim at showing: I) COLA's 's effectiveness for modification minimality. II) P-SHAP's superior performance than other Shapley methods towards modification minimality. III) COLA's near-optimat performance.
The experiments are conducted with 4 datasets for binary classification tasks, CE algorithms that are designed for diverse goals, and 12 classifiers, shown Table 1, where combination of dataset, Acr algorithm, and a model defines an "experiment scenario".
Experiment Setup
"The code is available on both ?
https: : / /a /ancnymous. 4open. wcwience/r//ontrraaeeeeeeee -etprrreeetteeeereeereeeeete Attribbtt Atw.rtttttoooorreat....o oF-DFB1 and the submitted supplementary files.
7

---

Published as a conference paper at ICLR 2026
Table [Experiment Scenarios Setup] Four datasets are used to benchmark five CE algorithms over 12 models. A "scenario" is defined to be a combination of dataset, AcE algorithm, and a model f.
Dataset
HELOC (FICO, 2018), German Credit (Hofmann, 1994), Hotel Bookings (Antonio
et al., 2019), COMPAS (Jeff Larson et al., 2016)
ACE
DiCE (Mothilal et al., 2020), AReS (Rawal & Lakkaraju, 2020), GlobeCE (Ley et al., 2023), KNN (Albini et al., 2022; Contardo et al.; Forel et al., 2023), Discount (You ct al., 2025)
Model f
Bagging, LightGBM, Support Vector Machine (SVM). Gaussian Process (GP), Radial Basis Function Network (RBF), XGBoost, Deep Neural Network (DNN), Random Forest (RndForest), AdaBoost, Gradient Boosting (GradBoost), Logistic Regression (LR), Quadratic Discriminant
Analysis (QDA)
Table 2: [Experiment Methods Setupl The table defines 6 methods for comparisons, colored to align with Figure 3. Each method is put in experiment scenario, as defined in Table 1, for benchmarking.
Method
The probability p used by Avalue and Ashap
The Shapley algorithm AShap RB-SHAP, D = trainset of f
RB-Puni
p ApjobiPij = 1/nm (7i,
RB-Por
p Aq. equation (8) (but not used in Astap)
RB-SHAP, = trainset of f
CF CF-Puai
P t ApsobiPij = 1/nm (Vi,j)
CF-SHAP, D(x:) p
CF CRod
P Any x; matched randomly to an ry
CF-SHAP, D D(x:) = p:
CF-Por
P f- Eq. equation (8)
p-SHAP with p
CF-PEet
p AprobiAny x; matched known counterpart r,
(CF or B)-SHAP, D(x.) - ry
We briefly introduce the many ACE in Table 1. DiCE and KNN are data-insttance-bsee CE methods, which yield counterfactual(s) respectively for each factual instance. AReS and GlobeCE are group- based CE methods, which find a collection of counterfactual instances for the whole factual data as a group. The algorithm Discount treats the factual instances as an empirical distribution and seeks a counterfactual distribution that stays in proximity to it.
Table 2 defines 6 methods, where CF-Por is the proposed p- SHAP and the others are baselines. Each is put in many experiment scenarios in Table 1, benchmarked comprehensively, Each method is determined by combination of the three algorithms Aprob. A and Aship For example, the first row RB-pum uses uniform distribution as the algorithm Aproh to compute p. and then the computed p is sequentially used in equation (9) or equation (10) of Avalue for computing q, and Ashap is RB-SHAP.
These methods are carefully designed for ablation studies. First, RB-PUni differs with CF-PUa in that the latter uses CE information whereas the former does not. Second, CF-PUmi CF-PRed, and CF-Por use CE in FA with different joint distributions, and we demonstrate later that the distribution computed by OT outperforms the others significantly. Third, want to make sure that OT is useful because it informs Ashap with respect to the factual-counterfactual alignment, not because of other factors, and hence comparison of CF-Por that uses such an alignment RB-por that does not. Finally, CF-P represents a special case that each counterfactual originates from a known source, which is used as an exact couttal-outttrfactuulll aligmment, making CF-SHAP also B-SHAP.
Result I: COLA achieves significant action reduction with a minor loss in counterfactual effect In Table 3, we set a goal for such that f(x) reaches 80% or 100% counterfacrtual effect of (rr) Observe that COLA is effective in achieving this goal, requiring significantly fewer actions in z (i.e., modifications of features in z). compared to the original CE r. Using COLA, one could expect to only resort to 13%-25% of the feature changes (calculated by ||z x||/||r - x||) to achieve the counterfactual effect of 80%. In particular, only COLA with SHAP (that is, CF-por) can reach the goal of the counterfactual effect of 100%, with only 45% of the feature changes in original r. Especially, PHAP leads to the best actional minimality, which is analyzed in details below.
NNoe that by definition, we have ((f(r), ) = O. which represents a 100% counterfactual effect since y f(r). To define counterfactual effect 80%, consider the proportion of divergence reduced by the refined CE z. That is, Counterfactual Effect 1 ()f(z), D(f(x), v²) 80%, with 1) heing OT.
8

---

Published as a conference paper at ICLR 2026
Table [COLA for Modification Minimality] This table shows the number of modified features in by each method, when (z) reaches 80% counterfactual effect of fr).. The result of each method is averaged by running in 4 randomly selected scenarios in Table 1, with 100 runs each scenario. The symbol "_" means the target counterfactual effect cannot be achieved.
Method 80% Counterfactual Effect
100% Counterfactual Effect
Dataset
Modified Features
--xl # Modified Features
r-x x
RB-Puni
German
RB-por
5.29
75.9%
Credit
CF-Pusi
[F]]
-R-PRad
CF-por
1.70(#0.02)
24.3%
8.13(=0.08)
44.9%
RB-Puni
7.10(+0.08)
24.5%
Hotel Bookings
RB-poT
8.55(+0.08)
50.2%
CF-PUni
7.01(++0.07)
41.1%
[F
29
CF-PRod
10.63(+0.08)
62.4%
CF-por
2.50(=0.03)
14.6%
(.4400.02
26.0%
RB-PUni
(.020)
82.7%
RB-POT
COMPAS
CF Puni
2.80(+0.04)
34.4%
F 15
CF-PRad
2.58 +0.04)
32.1%
CF-por
1.25 #0.03)
14.8%
2.45 00.03)
30.0%
RB-Puni
RB-Por
HELOC [F| = 23
CU-PUai
CF-PRdd
2.70(+0.04) 2.35(+0.03)
15.7%
CF-POT
13.4%
7.745(10.05)
44.7%
Result II: p-SHAP outperforms the other Shapley methods in achieving counterfactual effect We provide the evaluation of different Shapley methods in equation (4)-equation (7) in Figure 3, where the x-axis the number of allowed feature changes C and the y-axis is the term ())(2), in equation (I), First, RB-PUai and RB-PoT perform significantly worse than the others, indicating the importance of using CE information attributing features towards modification minimality in CE. Second, the result showing that CF-Por outperforms RB-Por demonstrates that the use of OT enhances the performance of p-SHAP specifically by providing effective Dectual-counteraccuu alignment, rather than being influenced by other factors, due to that they only differ in Ashap (see Table 2). Third, p-SHAP significantly outperforms CF-puni and CF-pRod- This shows the effectiveness of the joint distribution obtained in OT in using S-HH Namely, merely using the counterfactual information for FA (as done by CF-Puni and CF-PRad) is not enough, and a proper alignment (which does not necessarily mean the defined by the exact counterfactual generation mechanism, as revealed in Figure later) between factual and counterfactual must be considered.
We observed that COLA with S-SP possess good performance for many of the experiment scenarios defined in Table 1, as shown in Figure 3 as well as Table 3. massive amount of such results are further demonstrated in Appendix H. In conclusion, PHSP outperforms all the other Shapley methods for the sparsity in CE.
Result II:: COLA competitive with the mixed integer linear programming (MILP) optimum on tractable benchmarks This result demonstrates the effectiveness of p-SHAP in eliminating the influence of the CE generation process by replacing the CE algorithm-deependent knowledge" of D with the OT joint distribution between the factual and counterfactual data, shown in Figure 4. We benchmark the method CF-PEc using COLA, and focus on solving equation (1) with MeanD as known cettuall-ounterfacttu alignment,
the divergence function D. Since CF-PE>: relies benchmark the effectiveness of COLA for using this alginment. The theoretical optimality of COLA in this case can be obtained by solving MILP, Appendix G for how the MILP formulation is derived. Note that solving MILP is computationally heavy, hence only done for German Credit.
"This knowledge, i.c. exact alignment between factual and counterfactual, is available only by DiCE and KNN among all CE methods considered.

---

Published as a conference paper at ICLR 2026
LightGBM, Acr =KNN
RBF, Ace = -Diccount
XGBoost, AcE =DICE
GP, AcE =DICE
0.8
0.8
0.8
[German Credit)
0.6
0.6
6
0.4
0.4
0.5
0.2
0.2
0.2
0
100
200
100
200
0
100
200
100
200
DNN, ACE DisCount
AdaBoost, ACE -DiCE
LR. ACE KNN
GradBoost, ACE
KNN
(Hotel Bookings
0.8
MMD
0.6
0.4
0.5
0.5
0.5
0.2
0
0
0
200
400
200
400
0
200
400
200
400
Bagging, Acc =ReeS
SVM. Ace AReS
LightGBMM, Acr GGlobeCB
RBF, Ace Discount
0.8
0.8
0.8
MeanD
0.6
0.6
0.6
0.8
[COMPAS]
0.4
0.4
0.4
0.6
0.2
0.2
0.2
0.4
0
0.2
0
50
100
50
100
0
50
50
100
DNN, ACE Discount
RndForest, AcE KNNN XGBoost, AcE GGlobeCE
QDA, ACE GGlbbeE
1
RB-pun
MedianD
[HELOC)
RB-pur
0.5
0.5
CF-pus
0.5
0.5
CF-pess
CF-por
0
0 200 400 600 800
200 400 600 800
0 200 400 600 800
200
400 600
Figure 3: ((f(z) y") vs. allowed actions C. Experiments are with 100 runs. The shadows show the 99.9% confidence intervals. A is used for HELOC and COMPAS, and A is used for German
Credit and Hotel Bookings. The legend inside "RnDForest, ACE KNNN" applies all plots.
XGBoost, AcE= DiCE
LightGBM, ACE KNNN
DNN, ACE KNN
AdaBoost, ACE DiCE
CF-por
0.8
0.8
CF-Pea
0.6
0.5
Optimal
0.4
0.5
0.4
0.2
0.2
0
0
0
100
200
0
100
200
0
100
200
0
100
200
Figure 4: [German Credit) D(f(=), y²" vs. allowed actions C, with 1) being MeanD.
We that CF-Pra possess near-optimal performance using DiCE. Remark that for DiCE and KNN, the factual-counteefaatuaa pairs are independent to each others and hence we argue that COLA effective for instance-based CE, even though our formulation in equation (1) a generalization for group or distributional CE. It is interesting to note that CF-Por sometimes performs better than CF-Pra. because it utilizes more theoretically grounded approach identify the key features that require modification, whereas CF-PEa relies on CE algorithm-dependee solid justification on its effectiveness for FA. We notice that there is still a gap between CF-PFet. CF-Por and the optimal result in KNN. Finding the best alignment is still an open question.
knowledge, which lacks
CONCLUSIONS
This paper introduces a novel framework, COLA, for refining CE by joint-distributionnnniee Shapley values, ensuring the refined CE maintains the counterfactual effect with fewer actions.
10

---

Published as a conference paper at ICLR 2026
ETHICS STATEMENT
All authors have read and agree to abide by the ICLR Code Ethics. This paper proposes method (COLA with p-SHAP) for refining CE. Our work is methodological and empirical: it does not involve human subjects or interventions, and we did not collect, annotate, or release any new personal data.
All experiments use standard, de-identified, publicly available tabular datasets (c.g.. Adult/Census Income, German Credit, HELOC), We follow dataset licenses/terms and the common research practice of using these data solely for benchmarking. We do not attempt re-identification linkage, and we release no additional personal information. Our code will include options to exclude protected attributes and to enforce feasibllty/immuuabiiitty constraints on features.
CE can surface biases in underlying models but do not by themselves guarantee fairness. Our method reduce the number of recommended feature edits; it does not authorize edits to sensitive immutable attributes (c.g., sex, race, age) nor does it confer causal validity. We discourage use of CE "game" high-stakes systems (e.g., lending, hiring, healthcare) or pressure individuals into unrealistic or unethical behavior changes. Any deployment in consequential domains should include domain-exte -eaneeeeeeeettttiiieeee feasibility constraints, fairness audits (e.g.. disparate impact/benefit across groups), human oversight, and compliance with applicable data cnta-protection and Al regulations, We explicitly caution that our attribution distributional (not causal) and must not be interpreted as causal effect.
Risks include (i) misinterpretation attributions as causal, (ii) recommending infeasible or harmful edits, and (iii) enabling strategic manipulation. We mitigate these by (a) stating non-causal scope and limitations the paper, (b) supporting cost/feasibility masks and immutable features in the implementation, (c) reporting performance across subgroups where relevant, and (d) providing guidance for responsible use in documentation.
We provide anonymized implementation, detailed settings, and seeds to facilitate replication; report model types, hyperparameters, and compute budgets, and we avoid result cherry-picking by using fixed evaluation protocols. To knowledge, there are conflicts of interest sensitive sponsorships related to this submission; any funding disclosures will be provided after the double- blind review.
REPRODUCIBILITY STATEMENT
We take reproducibility seriously and provide pointers to all necessary components. The COLA framework and p-SHAP attribution are specified with pseudocode and notation in the main paper (Algorithm Method and Preliminaries sections), with all theoretical assumptions and complete proofs deferred to the appendix (theorems and lemmas with line-by-line proofs).
Implementation details hyperarammterrss,, optimization settings, environment specifications, and random seeds-are are documented in Appendix I, which also describes the datasets used and the evalua- tion protocol (including data splits and preprocessing). Baselines (models and CE generators) and coupling choices are enumerated in the experimental sections and ablations, with their configurations mirrored in the supplemental materials. A downloadable code repository
https: / / /ii thub. cww./oulll///o/ooooooooooo...... wowl/yolllll/e e 10/2//XA-COOLA
contains scripts to reproduce all results end-to-end, including a requirementtsccviroomeee file and seed-controlled runs; also report hardware and runtime in Appendix I to contextualize compute requirements, and provide ablation/sensitivity studies and the small-scale MILP benchmark setup in the appendix to support independent verification.
REFERENCES
OR-Tools URL https:t/////.. /ttp///tpp//eeeeeooorroooooooo/p
Gurobi Optimization, LLC. URL wtt::tt:::/www..grrot
11

---

Published as a conference paper at ICLR 2026
Machine Learning in
URL https: /ttpptpps///aattttlaaaoooo////aa
Emanuele Albini, Jason Long. Danial Dervovic, and Daniele Magazzeni. Counterfactual shapley additive explanations. Proceedings of the 2022 ACM Conference on Fairness, Accountability, and Transparency, pp. 1054-1070, 2022.
Nuno Antonio, Ana de Almeida, and Luis Nunes. Hotel booking demand datasets. Data in brief, 22: 41 41-49, 2019.
Alejandro Barredo Arrieta, Natalia Díaz-Rodriguez, Javier Del Ser, Adrien Bennetot, Siham Tabik, Alberto Barbado, Salvador Gareía, Sergio Gil-López, Daniel Molina, Richard Benjamins, et al. Explainable artificial intelligence (xai): Concepts, taxonomies, opportunities and challenges toward responsible ai. Information fusion, 58:82-115, 2020.
Emilio Carrizosa, Jasone Ramirez-Ayerhe, and Dolores Romero Morales. Mathematical optimization modelling for group counterfactual explanations. European Journal of Operational Research, 2024.
Tianqi Chen and Carlos Guestrin. Xgboost: A scalable tree boosting system. In Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining, pp. 785-794, 2016.
Claudio Contardo, Ricardo Fukasawa, Louis-Martin Rousseau, and Thibaut Vidal. Optimal coun- terfactual explanations for k-nearest neighbors using mathematical optimization and constraint programming. URL https: http///ptiiiiiiiiiiiiiiieeeetett line. www.liineeegg.............
Arun Das and Paul Rad. Opportunities and challenges in explainable artificial intelligence (xai): A survey. arXiv preprini arXiv:2006.1337/, 2020.
FICO. Fico explainable machine learning challenge, 2018. URL https: h/op://it.//m/ty..////// fico. machine Leannneeearnnnnnng challenge.
com/ s exppainaaaaee wawaaaallemmcciieeeeeee
Alexandre Forel, Axel Parmentier, and Thibaut Vidal. Explainable data-driven optimization: From context to decision and back again. In International Conference on Machine Learning, pp. 10170- 10187. PMLR, 2023.
Riccardo Guidotti. Counterfactual explanations and how to find them: literature review and bench- marking. Data Mining and Knowledge Discovery, pp. 1-55, 2022.
Hans Hofmann. Statlog (German Credit Data). UCI Machine Learning Repository, 1994. DOI: http://ddo..rr11024433///
Lauren Kirchner Jeff Larson, Surya Mattu, and J Angwin. How we analyzed the com-
pas recidivism algorithm, 2016. URL https: //wppp::/wwww.. wwwppplica.aac......p. wtoppplic...o////tttct..... org/article/ hatpaaaaaaaaaaat---ee--eeeeerdee algooiimaaaaaiiiammmmmmma
Kentaro Kanamori, Takuya Takagi. Ken Kobayashi, and Hiroki Arimura. Dace: Distribution-wwace counterfactual explanation by mixed-integer linear optimization. In IJCAI, pp. 2855-2862, 2020.
Kentaro Kanamori, Takuya Takagi, Ken Kobayashi, and Yuichi Ike. Counterfactual explanation trees:
Transparent and consistent actionable recourse with decision trees. In International Conference Artificial Intelligence and Statistics, pp. 1876-1870. PMLR, 2022.
Amir- Moir-Hosseiin Karimi, Bernhard Schölkopf. and Isabel Valera. Algorithmic recourse: from coun- terfactual explanations to interventions. In Proceedings of the 2021 ACM conference on fairness, accountability, and transparency, 353-362, 2021.
Guolin Ke, Qi Meng, Thomas Finley, Taifeng Wang, Wei Chen, Weidong Ma, Qiwei Ye, and Tie- Yan Liu. Lightghm: highly efficient gradient boosting decision tree. Advances neural information processing systems, 30, 2017.
Ramaravind Kommiya Mothilal, Divyat Mahajan, Chenhao Tan, and Amit Sharma. Towards uni- fying feature attribution and counterfactual explanations: Different means to the same end. In Proceedings of the 202/ AAAI/ACM Conference on AI, Ethics, and Society, pp. 652-663, 2021.
12

---

Published as a conference paper at ICLR 2026
Yongchan Kwon and James Y Zou. Weightedshap: analyzing and improving shapley based feature attributions. Advances in Neural Information Processing Systems, 35:34363-3453 2022.
Dan Ley, Saumitra Mishra, and Daniele Magazzeni. Global counterfactual explanations: Investiga- tions, implementations and improvements. arXiv preprint arXiv:2204.069) 06910, 2022.
Dan Ley, Saumitra Mishra, and Danicle Magazzeni. Globe-ce: translation based approach for global counterfactual explanations. In International Conference on Machine Learning, pp. 19315-19342. PMLR, 2023.
Scott M Lundberg and Su-In Lee. A unified approach to interpreting model predictions. Advances in neural information processing systems, 30, 2017.
Luke Merrick and Ankur Taly. The explanation game: Explaining machine learning models using shapley values. In Machine Learning and Knowledge Extraction: 4th IFIP TC 5, TC 12, WG 8.4, WG 8.9. WG 12.9 International Pooss-Domain Conference, CD-MAKE 2020, Dublin, Ireland, August 25-28, 2020, Proceedings 4, pp. 17-38. Springer, 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 Iransparency, pp. 607-617, 2020.
Yasunobu Nohara, Koutarou Matsumoto, Hidehisa Soejima, and Naoki Nakashima. Explanation of machine learning models using shapley additive explanation and application for real data in hospital. Computer Methods and Programs in Biomedicine, 214:106584, 2022.
Axel Parmentier and Thibaut Vidal. Optimal counterfactual explanations in tree ensembles. In International conference on machine learning, pp. 8422-8431. PMLR, 2021.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. Pytorch: An imperative style, high-performance deep learning library. Advances in neural information processing systems, 32, 2019.
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.
Amit Sharma, Hua Li, and Jian Jiao. The counteffactual-saaleey value: Attributing change system metrics. In NeurlPS 2022 Workshop on Causality for Real-world Impact.
Yi Sun and Mukund Sundararajan. Axiomatic attribution for multilinear functions. In Proceedings of the 12th ACM conference on Electronic commerce, pp. 177-178, 2011.
Mukund Sundararajan and Amir Najmi. The many shapley values for model explanation. In International conference on machine learning, pp. 9269-9278. PMLR, 2020.
Lucile Minasiina,, Sahra Ghalebikesabi, Karla Diaz-Ordaz, and Chris Holmes. Challenges and opportunities of shapley values in clinical context. arXiv preprint arXiv:2306. 14698, 2023.
Berk Ustun, Alexander Spangher, and Yang Liu. Actionable recourse in linear classification. In Proceedings of the conference on fairness, accountability, and transparency, pp. 10-19, 2019.
Sahil Verma, John Dickerson, and Keegan Hines. Counterfactual explanations for machine learning:
A review. arXiv preprint arXiv:2010. 10596, 2, 2020.
Sandra Wachter, Brent Mittelstadt, and Chris Russell. Counterfactual explanations without opening the black box: Automated decisions and the gdpr. Harv. JL & Tech., 31:841, 2017.
Lei You, Lele Cao, Mattias Nilsson, Bo Zhao, and Lei Lei. Distributional counterfactual explanations with optimal transport. In The 28th International Conference Artificial Intelligence and Statistics, 2025.
Lin Zhu and Lei You. xai-cola: A python library for sparsifying counterfactual explanations 2026. URL https: // //wrxivv..v////////v/////////vvv org/ absaabs/222 220.2 21845.
13

---

Published as a conference paper at ICLR 2026
APPENDICES
Comparison Existing Approaches
14
APPhhadneess the Problem in equation (1)
14
Proof of Theorem 4.1: y-SHAP Towards Counterfactual Effeet
15
Proof of Theorem 4.2: Interventional Effect
17
Proof of Theorem 5.1: Counterfactual Proximity
18
F Analysis of Computational Complexity of COLA
19
G An MILP Formulation of equation (1) With MeanD
19
H Extended Numerical Results
20
1 Experiments Reproducibility
21
A COMPARISON TO EXISTING APPROACHES
The authors in (Albini et al., 2022) proposed CF-SHAP, which uses counterfactual data points the background distribution for Shapley. Yet, it assumes the counterfactual data distribution is defined conditionally each single data instance, which implies that there is known probabilistic alignment between every factual and every counterfactual instance. Similar assumptions are made in (Merrick & Taly, 2020; Kommiya Mothilal et al., 2021). In many scenarios where global explanations expected, this assumption fails. We note that the setup in (Albini et al., 2022) for counterfactual data distribution is special case of ours. The authors (Kwon & Zou, 2022) proposed WeightedSHAP adding weights to features rather than treating them as contributing equally. Our proposed method weights the contributions for data points and can be straightforwardly extended to consider weighting both rows and columns. Literature (Kommiya Mothilal et al., 2021) establishes framework for utilizing both FA and CE for explainability. Yet, the CE-based FA have the same assumption as (Albini et al., 2022), making it difficult to generalize to group (Ley et al., 2023; 2022) or distributional CE (You et al., 2025) cases. Another relevant work is (Sharma et al.), a method that models counterfactuals using system's causal structure and time-series predictors to fairly decompose single observed change system-level metric into contributions from individual inputs. More importantly, the aforementioned literature does not address the minimal actions CE problem, which the focus of our paper.
The problems formulated in (Kanamori et al., 2022; Karimi al., 2021) quite close to the one investigated in this paper. In (Karimi et al., 2021) the authors minimize the cost of performing actions with assumptiosn of known structural causal model (SCM), which is rarely known in practice. The authors in (Kanamori et al., 2022) pointed out that it remains open whether existing CE methods can be used for solving that problem. An-MILP-solveee based approach proposed for linear classifiers, tree ensembles, and deep ReLU networks, built upon the works (Ustun et al., 2019; Kanamori et al., 2020; Parmentier Vidal, 2021). However, solving MILP is costly, which makes it difficult scale.
B NP-HARDNESS OF THE PROBLEM IN EQUATION (1)
Theorem B.1 below states that one does not expect scalable exact algorithm for solving it generally, unless P = N', and the hardness lies not only on the non-linearity of f. but also its combinatorial nature.
Theorem B.1. Problem equation (1) is generally N'P-hard for non-trivial divergences. More specifically, it is hard even when d for linear models.
Proof. Consider the Sparse Regression (SR) with Cardinality Constraint problem, defined as follows. Given a matrix W Rmxn a target vector y* E RT and a sparsity level K E N, the goal is to find a vector E R that minimizes the residual error while having most K non-zero
14

---

Published as a conference paper at ICLR 2026
elements:
min
[z y"
s.t. l S K,
where |zz|oo denotes the number of non-zero elements in This problem is known to be P'r-hard due to the combinatorial nature of selecting the subset of variables to include the model.
We will map this SR problem of features be d = 1 and the number of instances be n (the same as the dimension of the SR problem). Let the factual data be given in the SR problem). We define the model f as a linear function [(z) Wz.
problem equation (I) with the following settings. Let the number
0 R" (the zero vector), and the target model output be y" E R (as
For the model output, define D(f(z),y") [l y" 12, i.e. the Euclidean distance squared. For the instances, we define D x) 2 x11 1 since x = O. We set the maximum allowed feature changes C' to be equal k (the sparsity level from the SR problem). The large constant M can be any sufficiently large positive number, for example, M 2 max; |Za|-
Given that
0, constraints equation (Id) and equation (le) simplified to:
2 0 (1 c2) + Mc Me
2, 2 0- c,) Mc. = Mc
4 k = I....,
This means that if c = 0, then 2, I 0 and 2; 2 0, so 2, = 0. If c = 1, then 2, E { -MMM The constraint equation (Ic) becomes Σk=1 C k.
Our problem equation (1) thus becomes:
min 0,2
||z y 11
(11a)
s.t.
||a122 E,
(1lb)
Σl Gi S k, 2 = 0, 2, E [-M, M),
(11c)
if 04 0,
(11d)
if c; 1, Vk I....; n.
(1le)
4 E {0,
(11f)
In this formulation, the variables indicate whether the variable is allowed to change (c, = 1) or
not (c: = O). The constraints enforce that 2; = 0 when 6 = O, mirroring the sparsity constraint in the SR problem. The constraint Σk=1 Q S k ensures that at most k features can change, matching the sparsity level. The objective function is identical to that of SR.
Since our problem formulation directly mirrors the SR problem with cardinality constraint, which known to be Aharra solving Problem equation (I) is at least as hard as solving the SR problem. Therefore, Problem equation (1) is W'P-ha even when d = 1, the model f is linear, and the divergence functions are standard Euclidean distances.
C PROOF OF THEOREM 4.1: P-SHAP TOWARDS COUNTERFACTUAL EFFECT
Let (x+],_ Raxd be the of factual data points with associated probability weights /4; 2 such that Di #4 1, and let {r
1m E gmxd be the set of counterfactual data points with
associated probability weights vj 0 such that ΣØ= Uj 1. Let por E mxxm be the OT plan
between x and that minimizes the expected transportation cost:
Por arg min (E
ΣΣ
Pis |x 1
where 11(24, v) is the set of joint distributions satisfying the marginal constraints Σ Pi = P4 for all i and Σ Piy My for all j. The theorem below provides that feature attributions are aligned with the expected costs of feature modifications, leading to action plans that are cost-efficient in achieving counterfactual outcomes.
15

---

Published as a conference paper at ICLR 2026
Theorem C.1 (Theorem 4.1 in the main text). Consider the 1-Wasserstein divergence W1, i.e. (1( = wil Ei=1 -Ø=1 Tij fxx) - -j]. Suppose the counterfactual outcome y" is fully achieved by r, i.e. = (j 2,2... m.). Assume that the model f -2 + R is Lipschitz continuous with Lipschitz constant The expected absolute difference in model outpuss between the factual and counterfactual instances, weighted by por (with 0), is bounded by:
f(r,)
y)( ΣE
P
r3113 S L
Σ
P ryl1 EI
Namely, Por minimizes the upper bound of (1(f(xx) where the upper bound is based on the expected feature modification cost.
Proof. Since the model f is Lipschitz continuous with constant L, for any E R2 and ry E Rd. it holds that:
lf(x) f(rr)) S L||x ryl22.
Multiplying both sides of the inequality by Pij 2 0, we obtain:
Psj fxx:) - f(,) S Lp x; - 12l12-
Summing both sides over 1,..., and j = 1, ,,....... m, we have:
PM f((x)) - f(r,))
|x; r5512-
Let us denote Ey Ei=1 inequality then becomes:
0T j=1 Pij
[xx;) f(r,)) and Ea Ei Dj=1 Pij Xx
ryl12. The
Ej S LEa.
To further bound E,,, we apply the Cauchy-Schwarz, inequality. Observe that the weights pij non-negative and satisfy Σx E P.j = because por a probability distribution over the joint space of x and The Sauchy-Schward inequality states that for any real-valued functions a,j and bij.
Σ arybij
Σ i.j
Σ f.J
b2
1
Setting dij
V Pij and bij
v 11 11
r3ll2, have:
Ea
OT Xx; ryyl= Σ
V
Pij
[x
ry112
Σ
asybig.
i,
Applying the Cauchy- inequality:
E3
S
Σ
a
b
Σ
Pij
Σ ps
r2,12
2
dnd
Since Ees Pij = 1, this simplifies to:
x
11
Pij
OT
Xx rr,13
Taking the square root of both sides, we obtain:
n
Fat
S
Σ
OT Pij
X
-r
j=1
16

---

Published as a conference paper at ICLR 2026
Substituting back into the inequality for E we have:
B S
OO X; r
Note that the 1-Wasserstein divergence no more than E then,
12
Wi(f(x),y"
Ty [(xx) yj E
X x
Therefore, the 1- 1-Wasserstein divergence between f(x and y" is bounded by the Lipschitz constant times the square root of the expected transportation cost under pot- Since por minimizes the expected transportation cost Ci, Pijl Pijl|xi- -ryl2 over all feasible transport plans in II( v), we have
ryl1
Piy ||
ry r,,
vp M(,,
i=
==
Hence the conclusion follows.
D PROOF OF THEOREM 4.2: INTERVENTIONAL EFFECT
Theorem D.1 (Theorem 4.2 in the main text). For any subset S C F and any (i = 1.2.... n), (S) represents the causal effect of the difference between the expected value of f() under the intervention on features S and the unconditional expected value of f(r). Mathematically, this is expressed as:
x1,))l
Proof. Let p(x, r) be the joint probability of x and r obtained from APrb. Under the intervention do(rs x, the features are set to x,,s, s. and the features in F remain distributed according to their marginal distribution p(rr)s). Therefore, the expected value of f(r) under the intervention is:
B(f(r) do(rs- = xi.s)) S
Remark that by the definition of (1)
(i0 [ -E] ff JRms f(x.ssttsss)(((sss)ttt) dryws B[f(r))
such that
+ = [[f(r) do (rs = xi,)]]
Hence the conclusion.
The value function (S) captures the expected value of f(r) when we intervene by setting the features S to x;ss denoted as do(rs = x;,s). This intervention independent of any predefined joint probability distribution p(x, r). Therefore, the expression [r) + ((() (S) represents the combined effect of the base expected value of fr) and the additional causal impact of the attribution (0 (S).
Note that F.y = Ey fx.)) - f(r))) Because æ and por denote joint probability and however, makes the summation the minimum for the Wi term across all possible joint distributions, whereas por does not.
17

---

Published as a conference paper at ICLR 2026
E PROOF OF THEOREM 5.1: COUNTERFACTUAL PROXIMITY
Theorem E.1 (Theorem 5.1 in the main text), Let q E RRaxd be the aligned reference matrix derived from the counterfactual set mxxd via the value assignment step q Avalue (r,P). where p is the coupling matrix. For any dimensions n, m 2 1, the refined counterfactual z constructed by the COLA framework satisfies:
2 x S ||9 x|
This inequality guarantees that the refined counterfactual z is strictly closer (or equal) to the factual input x compared to the aligned counterfactual proposal q, ensuring that the refinement process introduces no additional divergence.
Proof. Let be the aligned reference matrix generated by Avalue p). For the elements where Cak (modified elements), 2ik = 9sk- For elements where Gik- = 0, 2ik = Iik- Therefore, we can write:
(qa xaa) if Cak = 1 10
(Zik xzx)2
(12)
if cas- = 0.
The squared Frobenius norms with respect to q and z are computed as follows:
d
|/q x
quh
d
||z
x|
2i
Iakh 3
im
k=1
And,
|q 1--x
x
x
(qtk Tikk )
(Zik Tskh
1=
i=1 k=
Mik x2x)
(yih Fitk )²
(k,))oraa=1
(gat -xi) 2 0,
(1.))o0a 0
where the equality (i) holds because of equation (12).
Since the difference ||q x
||z
| 2 0, follows that:
x
||q x|
Taking square roots, we get:
||z
x
||q x||F.
Hence the conclusion.
Corollary E.2 (Corollary 5.2 in the main text). In the special case where n = m and the OT plan corresponds a deterministic permutation (i.e., obtained without entropic regularization, the matrix q becomes a row-permuted version of (denoted as To, where qi To(s)). In this setting, the bound simplifies to:
0),
zz - x S -
recovering the intuition that the refined counterfuctual at least as close to the factuals as the original counterfactual set reordered by a.
Proof. Under the assumptions n = m and E = O. the optimal transport plan becomes a permutation matrix associated with a bijection a. Consequently, the value assignment algorithm Avatue (whether maximizing or averaging) assigns q: = ro(i) for every instance i. Substituting q with the permuted matrix r, into the general bound established in Theorem 5.1 (i.e., ||z x||r S ||q - x||r) directly yields the result.
18

---

Published as a conference paper at ICLR 2026
F ANALYSIS OF COMPUTATIONAL COMPLEXITY OF COLA
We perform analysis of the computational complexity of COLA as follows.
First, we analyze Aprah. If the alignment between x and r is known priori, then Apmb just constructs the matrix with the prior knowledge, which takes On x m). Otherwise, we consider solving OT to obtain P. which, by the Sinkhorn Knopp algorithm, takes On x m x log(1/=)).
Then analyze Ashap. For each subset, the model is evaluated on all n. data points, leading to ((n) evluations per subset. Incorporating baseline values from the reference data r involves replacing the values of certain features with their corresponding baseline values. This operation is (m) because it requires accessing the baseline values from the reference table r for each of the features. If assume that the reference values can be precomputed and accessed in constant time, then the complexity of incorporating these values can be considered as Odd). The number of Mshap subsets results in Mshap model evaluations. Combining the above steps, the complexity of Ashap is On x d x MShap).
The normalization in line of COLA takes O(n d).
avg To compare the complexities of the two algorithms, Value and Avalee we analyze each algorithm step-by-step. For A Value for each row x; in the data table, we need to (1) compute the probabilities Pay for all j E {1,2... 2....... ................................ m), which involves (m)) operations per row, (2) identify the row F in the reference data with the highest probability Pij. which involves O(m) operations per row,
and (3) assign Qik Tr(() & where r(i) arg maxj Pij. which involves ((d) operations per row.
Since there are rows in the data table, the total complexity for A is Onn x (m + m + d))
On x (2m + d)) = On x m + n x d) = On (m + d)).
For A for each row in the data table, we need to 1) compute the probabilities Ps for all j E which involves (m)) operations per row, 2) compute the sum -y= Pij', which in- volves ((m) operations per row, and 3) calculate the weighted average Oik: Σド= Sr Pay Tyk. which involves (m x d) operations per row. Since there n rows in the data table, the total complexity for A is On x (m++++m x d)) = O(n (2m x d)) = O(n (m+mx x On x + d)) On X x d).
For lines 6-16, the entire complexity is straightforwardly O(n x d) + (C) = Onn x d) due to the fact C C n d.
Therefore, the complexity of COLA using A Value equals
((Mc)) + Omm log( (1/)) + Ondd Shap)
+ ((dd) + O(n(m + d)) + O(n))
-M(cc)) + Onm log((11))) + OIddMssaap) (nm) + O(d)),
and the complexity of COLA using A equals
((Mc) + O(am log [I)))) + O(ndMsmap) + ((n)) + (nmdd) + O(nd)
(Mc)) + O(am log [1/=)) + O(ddMshap)) + ((nd)).
Hence the complexity of COLA with respect n, m, d, and the regularization parameter e of entropic OT is
((Mcs) + O logm//)))(/))) (ndd + N
where N = (mm) + (nd) if Value is used and N (amd) if A is used.
G AN MILP FORMULATION OF EQUATION (1) WITH MEAND
In this section, we provide a global optimality benchmark for using a known alignment between factual and counterfactual in solving equation (1) with D being MeanD, namely
1 Σ f(z) - 1
D(f(z),y)) =
19

---

Published as a conference paper at ICLR 2026
with 3* = 1. 1= Vj. Since COLA is used, we have D(r,x) and z stays closer to x than
-ー
hence equation (1b) is dropped. The formulation of equation (1) then becomes:
min c22
E ((4) my
(13a)
n Σ
S.L.
Cik SC,
(13b)
1=1
2sk Cack + xx -x), 1= k=
(13c)
Note that the original constraints equation (Id) and equation (Ie) merge to be equation (13c), because CF-PEa is imposed to be used. That is, for any x,, there is exact ry serves its reference in Ashap and Avalue such that 2ik (k = 1, ......... d) either stays unchanged or can be changed to ryk- Therefore Zik Cik Tza (1 cas) of which the value depends on the binary variable cAk-
q is determined (also, both Amux and Values return the same q). For any data point and any feature set S S F, let z,s denotes the solution z, where we have all features k E S changed to lik. and the other features h E F stays 2ih. Hence the set of 2 (S C F) composes the domain of all possible values of z. Define a corresponding scaler variable for any Z,s:
Due to the known alignment between any x; and its corresponding
gis [(2s) y"
Then, for any z, equation (13), the value of the term Σ f(2) - ny" can be represented by binary variable a;s together with the scaler 9is, namely,
114)) - ny
ΣΣ
Iisas.
SC
The optimization problem equation (13) hence reformulated as mixed integer programming below.
min 8.7)
7
(14a)
s.t.
Giso:s S 17.
(14b)
gis0:s 2 -1]
(14c)
SCF
Σ ais = i= ==
(14d)
SCF
ΣΣ
Slais S C. $1
(14e)
Minimizing 7 under the two constraints equation (14b) and equation (14c) is equivalent to minimizing the objective function of equation (13). The constraints in equation (14d) guarantees that each data point is subject to one and only one feature modification plan 0gs for specific S (S C F). The constraint equation (14e) corresponds equation (13b). Solving equation (14) yields the theoretical optimality of COLA using a known alignment between factual and counterfactual, demonstrated in Figure in Section 6.
H EXTENDED NUMERICAL RESULTS
Observing Figures 5-8, CF-Pumi generally performs better than RB-Puni Second, consider RB-Por and CF-por that also differ only in Asuup the latter consistently outperforms the former. Hence RB-SHAP is not suitable for FA in CE.
We analyze how different Shapley methods affect FA, corresponding to lines 3-4 in COLA. The shapley methods be classified into two categories: First, consider RB-pum and CF-puai that differ
20

---

Published as a conference paper at ICLR 2026
only in Astap- Observing Figures 5-8, CF-pum generally performs better than RB-Puai Second, consider RB-por and CF-por that also differ only in Ashap, the latter consistently outperforms the former. Hence RB-SHAP is not suitable for FA in CE.
Besides FA, the other equally important step of COLA is line 5, i.. using the joint probability p(x, r) which direction change its features towards the target model outcome. We observe in Figures 5-8 that CF-Por consistently outperforms all other methods throughout all experiments. Note that all the three methods CT-pusi. CF-pRad, and CF-por provide solution's 's for the joint probability when the exact alignment between factuals and counterfactuals are unknown. Yet, their performance differ significantly. Simply knowing the CE (and its marginal distribution) is insufficient.
to compose the matrix q. telling the factual
move
OT proves be exceptionally useful when the alignment information between factual and counterfac- tual instances is missing or inaccurate. Even when the CE algorithm explicitly matches each factual instance corresponding counterfactual, it is challenging to justify that the known alignment optimizes performance. This is supported by Figure 4 in Section 6.
Note that por does not need to be the true joint distribution of x and r from a data generation perspective. Instead, it should guide COLA to treat x and r together for both FA and CE. Furthermore, the QDA column in Figure shows stableness of -b-sd methods, while others diverge significantly from the target. We emphasize that COLA, however, is not limited to using OT as Aprob. As indicated by Figure 4, any known best p still has non-negligible gap to the global optimality, Devising a better Aprob algorithm is hence an interesting topic worth exploration.
I EXPERIMENTS REPRODUCIBILITY
The experiments are conducted on high performance computing (HPC) cluster, running with four nodes (for the four datasets) in parallel, with each node equipped with two Intel Xeon Processor 2660v3 (10 core, 2.60GHz) and 128 GB memory. The experiment runs approximately 5-10 hours in each node, dependent on the size of the dataset. It is also possible to reproduce the experiment on a laptop. while it costs more computational time generally than using an HPC cluster.
For the four datasets, the numerical features are standardized, and the categorical features follow either label-encoding or one-hot encoding. Practically, we did not observe remarkable difference between the two encoding methods in terms of COLA'S performance. The train-test split follows 3.
The optimality baseline as shown in Figure 4 is solved by Gurobi 11.0.2 (gur). In order to reproduce the optimality baseline, a license of Gurobi is required. Otherwise, we can resort to open-source source operations research libraries such Google-OR tools (goo). We remark that solving the MILP in Appendix G is computationally expensive, such that it may only apply to small scale datasets such as German Credit. If one wants compute the optimality baseline for other datasets, then the number of used features needs to be reduced.
The hyperparameters of the models used in the experiment are specified follows. The models Bagging, GP, RBF, RndForest, AdaBoost, GradBoost, and QDA are scikit-learn models (skl), where all hyper-parameters are kept default. The models DNN, SVM. RBF, and LR are implemented by PyTorch (Paszke et al., 2019). The DNN has three layers. The SVM uses the linear kernel. The models XGBoost (Chen& & Guestrin, 2016) and LightGBM (Ke al., 2017) used by their scikit-leam interface, with all hyper-parameter kept default.
THE USE OF LARGE LANGUAGE MODELS (LLMS)
We used large language models (LLMs) only as general-purposs assist tools. Their role was limited grammar, wording, and light copy-editing of author-written text. LLMs did not contribute research ideation, modeling choices, experimental design, or results. All algorithms, proofs, datasets, and analyses were created and verified by the authors. Any LLM-suggested phrasing was reviewed and edited before inclusion in the paper. The authors take full responsibility for all content; LLMs are not authors or contributors. This disclosure complies with the ICLR policy on LLM usage.
21

---

Published as a conference paper at ICLR 2026
DNN, AcE Discount
RndForest. AcE
KNN
XGBoost, AcE GlobeCE
QDA, ACE E
1
0.6
RB-pumi
0.8
0.8
RB- RB-por por
0.6
0.4
0.5
0.4
0.6
0.2
0.2
0.4
0
0
0
200
400
600 800
0
200
400
600 800
0
200 400 600 800
0
200
400
600
1
CH-pusi
0.4
MMD
CR-prad
0.5
0.5
0.2
0.5
0
0
0
200
400
600
800
0
200
400
600
800
0
200
400
600 800
0
200 400
600
1
1
0.6
CF-por
0.8
0.8
MeanD
0.6
0.4
0.5
0.4
0.6
0.2
0.2
0.4
0
200
400
600 800
0
200
400
600 800
0
200
400
600 800
0
200 400 600
MedianD
0.5
0.5
0.5
0.5
0
0
0
0 200 400
600
800
0
200 400 600 800
0 200 400 600 800
0 200
400
600
Figure 5: [HELOC] D(f(z), y²) vs. allowed actions C. Experiments are with 100 runs. The shadows
show the 99.9% confidence intervals. The legends apply to all plots A is used.
Value
LightGBM, ACE KKNN
RBF, Ace -Discount
XGBoost, AcE DDICE
GP, Ace DICE
0.8
0.8
RB-puai
RB-por
0.8
0.6
0.6
0.6
5
0.4
0.4
0.5
0.4
0.2
0.2
0.2
0
0
0
100
200
100
200
0
100
200
0
100
200
0.8
CF-Puii
1
-F-Per
0.8
MMD
0.6
CF-PRed
0.6
0.5
0.4
0.5
0.4
0.2
0.2
0
0
0
100
200
100
200
0
100
200
100
200
0.8
0.8
0.8
0.6
0.6
0.6
MeanD
0.4
0.4
0.5
0.4
0.2
0.2
0.2
0
0
0
100
200
100
200
0
100
200
100
200
1
MedianD
0.5
0.5
0.5
0.5
0
0
100
200
100
200
0
100
200
100
200
Figure 6: [German Credit] D(f(z) vs. allowed actions C. Experiments are with 100 runs. The shadows show the 99.9% confidence intervals. The legends apply to all plots. Aavg used.
22

---

Published as a conference paper at ICLR 2026
DNN, Ace. DiiCCount
AdaBoost, AcE =DCC
LR, Ace =KNN
GradBoost, ACE =DICE
0.8
RB-puai
0.6
RB- RB-por
5 0.4
0.5
0.5
0.5
0.2
0
200
400
0
200
400
0
200
400
0
200
400
0.8
Cl-puui
1
1
MMD
0.6
CF-pRaD
0.4
CF- Por
0.5
0.5
0.5
0.2
0
0
0
200
400
200
400
0
200
400
0
200
400
0.8
1
1
0.6
MeanD
0.4
0.5
0.5
0.5
0.2
0
0
0
200
400
200
400
200
400
0
200
400
1
1
MedianD
0.5
0.5
0.5
0.5
0
200
400
200
400
200
400
200
400
Figure 7: [Hotel Bookings] D(f(z), y*) vs. allowed actions C. Experiments are with 100 runs. The shadows show the 99.9% confidence intervals. The legends apply to all plots. Amax is used.
Bagging, ACE =ReSS
SVM, Ace =AReS
LightGBM. ACE
-GlobeCE
RBF, Ace = Discount
0.8
0.8
0.8
0.6
0.6
0.8
0.6
5
0.4
0.4
0.4
0.6
0.2
0.2
0.2
0.4
0
0
0
0.2
50
100
0
50
100
0
50
100
0
50
100
0.8
MMD
0.6
0.6
0.5
0.4
0.4
0.5
0.2
0.2
RB -UBPU
0
0
50
100
0
50
100
0
50
100
0
50
100
0.8
0.8
0.8
0.6
0.6
0.8
McanD
0.6
0.4
0.4
0.4
0.6
0.2
0.2
0.2
0.4
0
0
0.2
RB-PoT
50
100
50
100
0
50
100
0
50
100
McdianD
0.5
0.5
0.5
0.5
CF-pui
CF-pRod
0
0
0
CF-por
50
100
50
100
0
50
100
0
50
100
Figure 8: [COMPAS] D(J(z), vs. allowed actions C. Experiments are with 100 runs. The shadows show the 99.9% confidence intervals. The legends apply to all plots. Amax is used.
Value
23


Reviews:
[{"review_id": "e2MEtNHf3E", "reviewer": "Reviewer_Puot", "summary": "The paper presents a method for obtaining sparse counterfactual explanations titled COunterfactuals with Limited Actions (COLA). The main novelty resides on the one side of it being a post-hoc method that refines existing counterfactuals to reduce their cardinality while minimizing the impact on the counterfactual's performance, and secondly to apply an Optimal Transport setting for guiding the refinement of the counterfactuals using Shapley attributions.", "questions": "- In the problem formulation (1a)-(1e), the authors state their model without but decide not to include the flipping of the target as an additional constraint. In fact, model (1a)-(1e) would by itself not produce any counterfactuals but merely minimum-distance synthetic samples. I think that should be included in the model.\n- In Section 4, it remains unclear to me what the authors mean with an i <-> j alignment for any x_i, r_j (line 145). I guess it refers to the matching step in Algorithm 1, but the paper would benefit from a clear definition of what alignment means in this case. In the subsequent line, it is also unclear what it means for an algorithm to be \"independent of CE\". Can the authors clarify this?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper's main strength is its combinaton of theory-guided framework with its extensive experimental setup and computational experiments. Another major strength of the work is the source code provided by the authors, which allows for reproducibility of the results end-to-end. The appendices are not used merely as extension of the paper but provide extensive theoretical and practical background about the study.", "weaknesses": "At times, the theoretical exposition remains a bit obscure, mainly due to the fact that the authors frame the problem around the group setting (generating a group of counterfactuals for a given group of factual instances) instead of in my opinion much common framing of the problem where a set of counterfactuals is generated for a specific factual instance. This makes the exposition harder for non-specialists. It is worth asking if the benefit (more generalizability) overweights the cost (clarity of exposition).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ufgcePI1FH", "reviewer": "Reviewer_Myko", "summary": "This paper introduces COLA, a post-hoc framework for refining counterfactual explanations by leveraging optimal transport (OT)-induced joint distributions to guide a novel p-SHAP attribution method. The key innovation lies in using OT to model the alignment between factual and counterfactual instances, which then informs Shapley-based feature attribution to select the minimal set of actionable edits that preserve the desired counterfactual outcome. Theoretically, the authors prove that OT minimizes an upper bound on the W₁ divergence between factual and counterfactual predictions under Lipschitz continuity, and that refined counterfactuals remain no farther from the original facts than the original CEs. Empirically, COLA achieves 26–45% fewer feature edits than baseline CEs across 4 datasets, 12 models, and 5 CE generators, with p-SHAP consistently outperforming other Shapley variants.", "questions": "Is the OT coupling used only to define the joint distribution p for p-SHAP, or does it also influence the actual value of the Shapley contributions beyond the reference distribution? Could you provide a small example where $p_{OT}$ leads to different attribution than $p_{Uni}$ or $p_{Rnd}$?\n\nHow does COLA perform when the initial CE is poor (high divergence from factuals, or incorrect predictions)? Could the refinement process amplify errors?\n\nWhat are the practical limits of the OT step in terms of n (number of instances) and d (features)? Have you tested COLA on larger datasets (>10k instances)?\n\nWhy not use the OT plan directly for edit selection? Instead of using $p_{OT}$ to compute p-SHAP, could one directly use the OT plan to select which features to modify by summing $p_{ij}$ over $j$)? How does this compare to p-SHAP?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "## Novelty and Theoretical Depth:\nThe integration of OT as a coupling mechanism to inform Shapley values is a fresh and well-justified approach. It moves beyond ad-hoc or model-specific alignments and provides a principled way to define meaningful contrastive references.\nThe p-SHAP framework generalizes B-SHAP, RB-SHAP, and CF-SHAP under a unified probabilistic coupling, offering a modular and flexible interface for attribution.\n\n## Strong Theoretical Guarantees:\nTheorem 4.1 is compelling: it establishes that OT minimizes an upper bound on the $W_1$ divergence between $f(x)$ and $y*$, linking the cost of feature modification directly to prediction fidelity. This is a non-trivial theoretical contribution.\nTheorem 5.1 provides a provable guarantee that refined CEs do not move farther from the factuals than the originals (in Frobenius norm).\n\n## Comprehensive Empirical Evaluation:\nThe experiments are rigorous and well-designed: 4 datasets, 12 models, 5 CE generators, multiple divergence metrics (OT, MMD, MeanD, MedianD), and ablations across different Shapley methods. The use of CF-$p_{OT}$ as the proposed method is clearly contrasted with baselines, demonstrating that alignment matters more than mere use of counterfactuals. The MILP-based optimality benchmark adds significant credibility, showing COLA’s near-optimality in a controlled setting.\n\n## Clear Motivation and Problem Framing:\nThe paper correctly identifies the de-coupling problem: standard FA methods (like Shapley) ignore the specific path to counterfactual outcomes, leading to suboptimal or even counterproductive edits. COLA addresses this by grounding attribution in the factual-counterfactual alignment.", "weaknesses": "## Ambiguity in the Role of OT in p-SHAP:\nThe paper frames OT as a means to define a joint distribution p, which is then used in p-SHAP. However, OT is not directly used in the Shapley computation.\n\n## Assumption of Known Counterfactuals:\nThe framework assumes that a counterfactual set r is already available via some CE method. While this is standard in post-hoc refinement, the paper does not discuss how errors in the initial CE propagate into the final refinement. A sensitivity analysis to noisy or suboptimal CEs would strengthen the claim.\n\n## Limited Discussion on Scalability:\nWhile computational complexity is analyzed, authors do not discuss practical bottlenecks: OT with entropic regularization is expensive for high dimensionality. Solving OT on 1,000 × 1,000 matrices might be feasible but slow. Authors assume OT is \"cheap\" via Sinkhorn, but does not report scalability limits.\n\n## Averaging Across Datasets/Models May Mask Heterogeneity:\nThe results are averaged across scenarios. While this is standard, it could obscure cases where COLA fails (like with highly non-linear models or sparse features). A per-scenario analysis (e.g., worst-case performance) would help assess robustness.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "1dAOklLITA", "reviewer": "Reviewer_2t4Y", "summary": "This paper addresses the problem of sparse counterfactual explanations: how to achieve a desired model prediction with minimal feature modifications. The authors propose COLA, a general post-hoc framework applicable to any model and counterfactual generator. The method introduces a novel joint distribution–informed Shapley attribution (p-SHAP), which uses an optimal transport coupling between factual and counterfactual samples to guide feature selection. Theoretical results show that OT minimizes an upper bound on the Wasserstein distance between model outputs and targets.", "questions": "1. Can the authors clearly articulate the new theoretical or methodological insight that distinguishes p-SHAP from prior Shapley-based counterfactual attribution methods?\n2. How does Theorem 4.1 translate into sparse modification guarantees in practice? Please provide an intuitive or empirical justification of this link.\n3. What is the computational complexity and runtime overhead of COLA compared to direct CE methods？\n4. Can COLA handle mixed data types (categorical + continuous) and feasibility constraints in real applications？", "limitations": "", "rating": 4, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "1. The paper tackles a meaningful and practical issue in explainable AI: generating concise and actionable counterfactual explanations.\n2. The problem and algorithmic design are mathematically well-structured.\n3. Experiments across several datasets and models show consistent improvements in sparsity and performance.\n4. The writing is generally clear and organized; the figures and tables support understanding of the method.", "weaknesses": "1. The proposed p-SHAP is largely a recombination of existing ideas (Optimal Transport + Shapley values). The conceptual contribution beyond combining these two paradigms remains limited. It is unclear what fundamentally new insight p-SHAP provides.\n2. The theorems mainly show that OT minimizes a Wasserstein upper bound, but this does not directly imply minimal feature edits in the discrete action space. The theoretical link between Theorem 4.1/5.1 and the claimed “minimal action” property is not convincingly established.\n3. All experiments are on small-scale tabular datasets. There is no evidence that the approach scales to more complex or high-dimensional data, nor comparison to recent counterfactual sparsity benchmarks.\n4. Since COLA may select arbitrary feature edits, how does it ensure the modifications are feasible or ethically valid (e.g., immutable features like age or gender)? The authors briefly mention this in the ethics statement but do not model it in the algorithm.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "VRYRt8rzty", "reviewer": "Reviewer_2QKe", "summary": "This paper introduces COLA (COunterfactuals with Limited Actions), a post-hoc framework that refines existing counterfactual explanations to use fewer feature modifications while maintaining the desired outcome. The key innovation is p-SHAP, a generalized Shapley value method that uses optimal transport (OT) to compute a coupling between factual and counterfactual instances, then leverages this coupling to identify which features are most critical to modify. Existing counterfactual explanation methods often modify more features than necessary to flip a model's prediction.", "questions": "The paper's theoretical contributions rest on restrictive assumptions that may not hold in practice. Theorem 4.1 requires Lipschitz continuity of the model f, yet many of the 12 tested models (particularly DNNs with ReLU activations) are not globally Lipschitz continuous. The authors do not verify which models satisfy this assumption or report how tight the theoretical bound is empirically. Theorem 5.1 only applies when n=m and ε=0 (deterministic matching), yet most experimental scenarios violate these conditions. Furthermore, while Theorem B.1 proves the problem is NP-hard even for simple cases, the paper provides no approximation guarantees or worst-case bounds for the proposed heuristic solution, leaving a significant gap between the hardness result and the practical algorithm's performance guarantees.\n\nSeveral core design choices in Algorithm 1 lack adequate justification. The probabilistic sampling mechanism in Line 7 introduces unnecessary randomness when deterministic top-C selection would be simpler and potentially more stable—no ablation study compares these approaches. The specific formulation of p-SHAP in Equation 7 is presented without explaining why this particular way of incorporating the coupling p into Shapley values is optimal compared to alternatives like weighted Shapley methods. Most puzzling is the observation in Figure 4 where CF-pOT sometimes outperforms CF-pEct despite the latter using the \"true\" alignment from the CE generator—this counterintuitive result demands deeper investigation as it questions fundamental assumptions about what constitutes the \"correct\" factual-counterfactual pairing. Additionally, the choice between A^max_Value and A^avg_Value appears arbitrary across datasets without clear selection criteria.\n\nThe experimental evaluation, while broad in scope (4 datasets, 12 models, 5 CE methods), lacks critical depth in several dimensions. Most notably, there is no computational cost analysis despite complexity bounds being provided in Appendix F—readers cannot assess whether COLA's sparsity improvements justify its computational overhead. The \"near-optimal\" claim relies solely on the German Credit dataset, where Figure 4 actually shows 20-40% gaps to optimality undermining this assertion. No statistical significance tests are provided, making it unclear whether observed improvements are meaningful or due to random variation. The extremely small standard deviations in Table 3 (±0.02-0.09) are suspicious given the sampling process and warrant verification. Crucially, the paper lacks comparisons to methods explicitly designed for sparse counterfactual generation, such as L1-regularized CE methods, sparse regression approaches, or greedy feature selection baselines—the current baselines only compare different Shapley variants within COLA.", "limitations": "", "rating": 4, "confidence": 5, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The paper provides a rigorous theoretical foundation for its approach with three well-formulated theorems. Theorem 4.1 establishes an important connection between optimal transport and counterfactual effect, proving that the OT-derived coupling minimizes an upper bound on the W₁ divergence between model outputs. This provides principled justification for why OT is the right tool for aligning factual and counterfactual instances. Theorem 4.2 demonstrates that p-SHAP captures genuine interventional effects, showing that the value function represents the causal impact of intervening on feature subsets. Theorem 5.1 offers proximity guarantees, proving that under certain conditions, the refined counterfactuals stay at least as close to the factuals as the original counterfactuals. These theoretical contributions are properly formalized with complete proofs in the appendices, demonstrating mathematical rigor and providing confidence in the method's foundations beyond purely empirical validation.", "weaknesses": "Weaknesses:\n\nTheory-practice gap: The main theoretical results assume Lipschitz continuity (Theorem 4.1) and deterministic OT matching with n=m (Theorem 5.1). These conditions are restrictive and don't hold generally. The paper doesn't adequately discuss when these assumptions are violated in practice.\nNP-hardness result (Theorem B.1) undermines the approach: If the problem is NP-hard even for d=1 linear models, why should we expect good approximations from the proposed heuristic? No approximation guarantees are provided.\nLimited theoretical justification for OT choice: While Theorem 4.1 shows OT minimizes an upper bound, it doesn't prove this leads to optimal feature selection for the discrete action problem in Eq. 1.\nExperimental methodology concerns:\n- Standard deviations in Table 3 are very small (±0.02-0.09), suggesting possibly limited diversity in runs\n- The \"near-optimal\" claim in Result III is based only on German Credit dataset - too limited to generalize\n- Figure 4 shows significant gaps between CF-pOT/CF-pEct and optimal, contradicting \"near-optimal\" claims", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)