Year: 2021

Abstract:
Counterfactual explanations are attracting significant attention due to the flourishing applications of machine learning models in consequential domains. A counterfactual plan consists of multiple possibilities to modify a given instance so that the model’s prediction will be altered. As the predictive model can be updated subject to the future arrival of new data, a counterfactual plan may become ineffective, or infeasible, with respect to the future values of the model parameters. In this work, we study the counterfactual plans under model uncertainty, in which the distribution of the model parameters is partially prescribed using only the firstand second-moment information. First, we propose an uncertainty quantification tool to compute the lower and upper bounds of the probability of validity for any given counterfactual plan. We then provide corrective methods to adjust the counterfactual plan to improve the validity measure. The numerical experiments validate our bounds and demonstrate that our correction increases the robustness of the counterfactual plans in different real-world datasets.

text: Counterfactual explanations are attracting significant attention due to the flourishing applications of machine learning models in consequential domains. A counterfactual plan consists of multiple possibilities to modify a given instance so that the model’s prediction will be altered. As the predictive model can be updated subject to the future arrival of new data, a counterfactual plan may become ineffective, or infeasible, with respect to the future values of the model parameters. In this work, we study the counterfactual plans under model uncertainty, in which the distribution of the model parameters is partially prescribed using only the first- and second-moment information. First, we propose an uncertainty quantification tool to compute the lower and upper bounds of the probability of validity for any given counterfactual plan. We then provide corrective methods to adjust the counterfactual plan to improve the validity measure. The numerical experiments validate our bounds and demonstrate that our correction increases the robustness of the counterfactual plans in different real-world datasets.

heading: 1 INTRODUCTION
text: Machine learning models, thanks to their superior predictive performance, are blooming with increasing applications in consequential decision-making tasks. Along with the potential to help making better decisions, current machine learning models are also raising concerns about their explainability and transparency, especially in domains where humans are at stake. These domains span from loan approvals (Siddiqi, 2012), university admission (Waters & Miikkulainen, 2014) to job hiring (Ajunwa et al., 2016). In these applications, it is instructive to understand why a particular algorithmic decision is made, and counterfactual explanations act as a useful toolkit to comprehend (black-box) machine learning models (Wachter et al., 2017). Counterfactual explanation is also known in the field of interpretable machine learning as contrastive explanation (Miller, 2018; Karimi et al., 2020b) or recourse (Ustun et al., 2019). A counterfactual explanation suggests how an instance should be modified so as to receive an alternate algorithmic outcome. As such, it could be used as a suggestion for improvement purposes. For example, a student is rejected from graduate study, and the university can provide one or multiple counterfactuals to guide the applicant for admission in the following year. A concrete example may be of the form “get a GRE score of at least 325” or “get a 6-month research experience”.
In practice, providing a counterfactual plan consisting of multiple examples is highly desirable because a single counterfactual to every applicant with the same covariates may be unsatisfactory in practice (Wachter et al., 2017). Indeed, the covariates can barely capture the intrinsic behaviors, constraints and unrevealed preferences of the person they represent so that the users with the same features may have different preferences to modify their input. As a consequence, a pre-emptive design choice is to provide a “menu” of possible recourses, and let the applicant choose the recourse that fits them best. Viewed in this way, a counterfactual plan has the potential to increase satisfaction and build trust among the stakeholders of any machine learning applications.
Constructing a counterfactual plan, however, is not a straightforward task because of the many competing criteria in the design process. By definition, the plan should be valid: by committing to any counterfactual in the plan, the application should be able to flip his current unfavorable outcome to a favorable one. However, each possibility in the plan should be in the proximity of the covariates
of the applicant so that the modification is actionable. Further, the plan should consist of a diverse range of recourses to accommodate the different tastes and preferences of the population.
Russell (2019) propose a mixed integer programming method to generate a counterfactual plan for a linear classifier, in which the diversity is imposed using a rule-based approach. In Dandl et al. (2020), the authors propose a model-agnostic approach using a multi-objective evolutionary algorithm to construct a diverse counterfactual plan. More recently, Mothilal et al. (2020) use the determinantal point process to measure the diversity of a plan. The authors then formulate an optimization problem to find the counterfactual plan that minimizes the weighted sum of three terms representing validity, proximity, and diversity.
A critical drawback of the existing works is the assumption of an invariant predictive model, which often fails to hold in practical settings. In fact, during a turbulent pandemic time, it is difficult to assume that the demographic population of students applying for postgraduate studies remain unchanged. And even in the case that the demography remains unchanged, special pandemic conditions such as hybrid learning mode or travel bans may affect the applicants’ package, which in turn leads to fluctuations of the covariate distribution in the applicant pool.
These shifts in the data are channeled to the shift in the parameters of the predictive model: when the machine learning models are re-trained or re-calibrated with new data, their parameters also change accordingly (Venkatasubramanian & Alfano, 2020). This raises an emerging concern because the counterfactual plan is usually designed to be valid to only the current model, but that is not enough to guarantee any validity on the future models. Thus, the counterfactual plan carries a promise of a favorable future outcome, nevertheless, this promise is fragile.
It is hence reasonable to demand the counterfactual plan to be robust with respect to the shift of the parameters. Pawelczyk et al. (2020) study the sparsity of counterfactuals and its non-robustness under different fixed models (predictive multiplicity). Rawal et al. (2020) consider the counterfactual plan problem and describe several types of model shift related to the correction, temporal, and geospatial shift from data. They also study the trade-off between the recourse proximity and its validity regarding to the model updates. Most recently, Upadhyay et al. (2021) leverage robust optimization to generate a counterfactual that is robust to some constrained perturbations of the model’s parameters. However, both works consider only the single counterfactual settings.
Contributions. We study the many facets of the counterfactual plans with respect to random future model parameters. We focus on a linear classification setting and we prescribe the random model parameters only through the first- and second-moment information. We contribute concretely
1. a diagnostic tool to assess the validity of a counterfactual plan. It provides a lower and upper bound on the probability of joint validity of a given plan subject to uncertain model parameters.
2. a correction tool to improve the validity of a counterfactual plan, while keeping the modifications to each counterfactual at a minimal level. The corrections are intuitive and admit closed-form expression.
3. a COunterfactual Plan under Ambiguity (COPA) framework to construct a counterfactual plan which explicitly takes the model uncertainty into consideration. It minimizes the weighted sum of validity, proximity and diversity terms, and can be solved efficiently using gradient descents.
Each of our above contributions is exposed in Section 2, 3 and 4, respectively. In Section 5, we conduct experiments on both synthetic and real-world datasets to demonstrate the efficiency of our corrections and of our COPA framework. All proofs can be found in the appendix.
General setup. Consider a covariate space Rd and a linear binary classification setting. Each linear classifier can be parametrized by θ ∈ Rd with decision output Cθ(x) = 1 if θ>x ≥ 0, and
0 otherwise, where 0 represents an unfavorable outcome. Note that we omit the bias term to avoid clutter, taking the bias term into account can be achieved by extending the dimension of x and θ by an extra dimension. A counterfactual plan is a set of J counterfactual explanations {xj}j=1,...,J , and we denote {xj} for short. When J = 1, we have a single counterfactual explanation problem, which is the subject of recent works (Ustun et al., 2019; Karimi et al., 2020a; Upadhyay et al., 2021).
Next, we define the joint validity of a counterfactual plan. Definition 1.1 (Joint validity). A counterfactual plan {xj} is valid with respect to a realization θ if Cθ(xj) = 1 for all j = 1, . . . , J .
Notations. We use Sd++ (Sd+) to denote the space of symmetric positive (semi)definite matrices. For any A ∈ Rm×m, the trace operator is Tr [ A ] = ∑d i=1Aii. For any integer J , [J ] , {1, . . . , J}.

heading: 2 VALIDITY BOUNDS OF COUNTERFACTUAL PLANS
text: In this section, we propose a diagnostic tool to benchmark the validity of a pre-computed counterfactual plan {xj}. We model the random model parameters θ̃ with a nominal distribution P̂. Instead of making a strong assumption on a specific parametric form of P̂ such as Gaussian distribution, we only assume that P̂ is known only up to the second moment. More specifically, we assume that under P̂, θ̃ has a nominal mean vector µ̂ and nominal covariance matrix Σ̂ ∈ Sd++. Definition 2.1 (Gelbrich distance). The Gelbrich distance between two pairs (µ1,Σ1) ∈ Rd × Sd+ and (µ2,Σ2) ∈ Rd × Sd+ is defined as
G ( (µ1,Σ1), (µ2,Σ2) ) , √ ‖µ1 − µ2‖22 + Tr [ Σ1 + Σ2 − 2 ( Σ 1 2 2 Σ1Σ 1 2 2 ) 1 2 ] .
The Gelbrich distance is closely related to the optimal transport distance between Gaussian distributions. Indeed, G ( (µ1,Σ1), (µ2,Σ2) ) is equal to the type-2 Wasserstein distance between two Gaussian distributions N (µ1,Σ1) and N (µ2,Σ2) (Gelbrich, 1990). It is thus trivial that G is a distance on Rd×Sd+, and as a consequence, it is symmetric andG ( (µ1,Σ1), (µ2,Σ2) ) = 0 if and only if (µ1,Σ1) = (µ2,Σ2). Using the Gelbrich distance to design the moment ambiguity set for distributionally robust optimization leads to many desirable properties such as computational tractability and performance guarantees (Kuhn et al., 2019). Motivated by this idea, we first construct the following uncertainty set
U , {(µ,Σ) ∈ Rd × Sd+ : G((µ,Σ), (µ̂, Σ̂)) ≤ ρ},
which is formally a ρ-neighborhood in the mean vector-covariance matrix space around the nominal moment (µ̂, Σ̂). The ambiguity set for the distributions of θ̃ is obtained by lifting U to generate a family of probability measures that satisfy the moment conditions

heading: B , {Q ∈ P : ∃(µ,Σ) ∈ U such that Q ∼ (µ,Σ)} ,
text: where P is a set of all probability measures supported on Rd and Q ∼ (µ,Σ) indicates that Q has mean vector µ and covariance matrix Σ.
The central question of this section is: If the distribution of θ̃ belongs to B, what is the probability that a given plan {xj} is valid? To answer this question, we define the event set Θ({xj}) that contains all model parameter values that renders {xj} jointly valid. Under the definition of a linear model, Θ({xj}) is an intersection of J open hyperplanes of the form
Θ({xj}) , { θ ∈ Rd : x>j θ ≥ 0 ∀j ∈ [J ] } . (1)
We name Θ({xj}) the set of favorable parameters. The probability of validity for a plan under a measure Q is Q(θ̃ ∈ Θ({xj})). We are interested in evaluating the lower and the upper bound probability that the plan {xj} is valid uniformly over all distributions Q ∈ B. This is equivalent to quantifying the following quantities
inf Q∈B Q(θ̃ ∈ Θ({xj})) and sup Q∈B Q(θ̃ ∈ Θ({xj})).
In the remainder of this section, we discuss how to evaluate the bounds for these terms. Lower bound. We denote by Θ◦ the interior of the set Θ, that is, Θ◦({xj}) ,{ θ ∈ Rd : x>j θ > 0 ∀j } . Note that all the inequalities defining Θ◦({xj}) are strict inequalities. By definition, we have Θ◦({xj}) ⊂ Θ({xj}), and hence infQ∈B Q(θ̃ ∈ Θ◦({xj})) ≤ infQ∈B Q(θ̃ ∈ Θ({xj})). Because Θ◦({xj}) is an open set, we can leverage the generalized Chebyshev lower bound to evaluate the minimum quantity of Q(θ̃ ∈ Θ◦({xj})) over all distributions with a given mean and covariance matrix (Vandenberghe et al., 2007). Adding moment uncertainty via the set U is obtained by rejoining two minimization layers. The next theorem presents this result.
Theorem 2.2 (Lower bound). For any ρ ∈ R+, µ̂ ∈ Rd and Σ̂ ∈ Sd+, let L? be the optimal value of the following semidefinite program
L? =  inf 1− ∑ j∈[J] λj s. t. µ ∈ Rd, Σ ∈ Sd+, C ∈ Rd×d, M ∈ Sd+ λj ∈ R, zj ∈ Rd, Zj ∈ Sd ∀j ∈ [J ] −x>j zj ≥ 0, [ Zj zj z>j λj ] 0 ∀j ∈ [J ]∑ j∈[J] [ Zj zj z>j λj ] [ M µ µ> 1 ] , [ Σ C C> Σ̂ ] 0, [ M − Σ µ µ> 1 ] 0
‖µ̂‖2 − 2µ̂>µ+ Tr [ M + Σ̂− 2C ] ≤ ρ2.
(2)
Then we have L? = inf Q∈B Q(θ̃ ∈ Θ◦({xj})) ≤ inf Q∈B Q(θ̃ ∈ Θ({xj})).
Upper bound. Because Θ({xj}) is a closed set, we can leverage a duality result to evaluate the maximum quantity of Q(θ̃ ∈ Θ({xj})) over all distributions with a given mean and covariance matrix (Isii, 1960). Adding moment uncertainty via the set U is obtained by invoking the support function of the moment set. This result is presented in the next theorem
Theorem 2.3 (Upper bound). For any ρ ∈ R+, µ̂ ∈ Rd and Σ̂ ∈ Sd+, let U? be the optimal value of the following semidefinite program
U? =  inf z0 + γ(ρ 2 − ‖µ̂‖22 − Tr [ Σ̂ ] ) + q + Tr [ Q ] s. t. γ ∈ R+, z0 ∈ R, z ∈ Rd, Z ∈ Sd+, q ∈ R+, Q ∈ Sd+, λ ∈ RJ+[ γI − Z γΣ̂ 12 γΣ̂ 1 2 Q ] 0, [ γI − Z γµ̂+ z γµ̂> + z> q ] 0[ Z z
z> z0
] 0, [ Z z
z> z0 − 1
] ∑ j∈[J] λj [ 0 12xj
1 2x > j 0
] .
Then we have supQ∈B Q(θ̃ ∈ Θ({xj})) ≤ U?.
Thanks to the choice of the Gelbrich distanceG, both optimization problems in Theorems 2.2 and 2.3 are linear semidefinite programs, and they can be solved efficiently by standard, off-the-shelf solvers such as MOSEK to high dimensions (MOSEK ApS, 2019). Other choices of distance (divergence) are also available: for example, one may opt for the Kullback-Leibler (KL) type divergence between Gaussian distribution to prescribe U as in Taskesen et al. (2021). Unfortunately, the KL divergence will entail a log-determinant term, and the resulting optimization problems are no longer linear programs and are no longer solvable using MOSEK. Equipped with L? and U?, we have the bounds
L? ≤ Q({xj} is a valid plan) ≤ U? ∀Q ∈ B
on the validity of the counterfactual plans {xj} under the distributional ambiguity set B. Complementary information. The previous results show that we can compute the lower bound L? and upper bound U? for the probability of validity by solving semidefinite programs. We now show that the two quantities L? and U? are complementary to each other in a specific sense. Proposition 2.4 (Complementary information). For any instance, either L? = 0 or U? = 1. More specifically, we have: (i) If µ̂ ∈ Θ({xj}), then U? = 1, and (ii) If µ̂ 6∈ Θ({xj}), then L? = 0.
Because L? and U? are bounds for a probability quantity, they are only informative when they are different from 0 and 1. Proposition 2.4 asserts that the upper bound U? is trivial when µ̂ ∈ Θ({xj}), while the lower bound L? becomes trival if µ̂ 6∈ Θ({xj}). Next, we leverage these insights to improve the validity of a given counterfactual plan.

heading: 3 COUNTERFACTUAL PLAN CORRECTIONS
text: Given a counterfactual plan {xj}, it may happen that {xj} have low probability of being valid under random realizations of the future model parameter θ̃. The diagnostic tools proposed in Section 2 indicates that {xj} have low validity when the bounds are low, and we are here interested in correcting this plan such that the lower bounds L? are increased. Indeed, increasing L? guarantees a higher confidence that the plan is valid, should the distribution of θ̃ belongs to the ambiguity set. At this point, one may be tempted to optimize L? directly with {xj} by first converting problem (2) into a maximization problem, and then jointly maximize with {xj} being decision variables. Unfortunately, this approach entails bilinear terms x>j zj in the constraints, and this approach is notoriously challenging to solve. We thus resort to heuristics for correction. Towards this end, the results from Proposition 2.4 suggest that there are two correction operations that we need to perform to improve the validity of the counterfactual plan:
(i) When µ̂ 6∈ Θ({xj}), then Proposition 2.4 suggests that we should modify the plan {xj} so that the resulting set of favorable parameters contains µ̂. We term this type of correction as a Requirement correction, and we consider one specific Requirement correction in Section 3.1.
(ii) When µ̂ ∈ Θ({xj}), we can also modify {xj} to as to increase the lower bound L?. This type of correction is termed an Improvement correction because its goal is to increase the validity of the counterfactual plans. We consider the Mahalanobis Improvement correction in Section 3.2.
We emphasize that the corrections of the plan {xj} are designed such that the modifications to each counterfactual xj should be minimal. This is achieved by two main criteria: the correction should modify as few counterfactuals as possible, and the modification to each counterfactual should also be as small as possible.

heading: 3.1 REQUIREMENT CORRECTION
text: We propose a Requirement correction with the goal of obtaining a corrected plan {x′j} from the given plan {xj} such that µ̂ lies inside (or strictly inside) the set Θ({x′j}). A simple Requirement correction is to construct the corrected plan {x′j} by
∀j ∈ [J ] : x′j = { xj if µ̂>xj ≥ , arg min{‖x− xj‖2 : µ̂>x ≥ } if µ̂>xj < ,
for some ≥ 0. Using this rule, x′j is the smallest modification of xj measured in the Euclidean distance such that x′j is valid with -margin with respect to the expected future parameter µ̂. The margin adds a layer of robustness: if > 0 then µ̂ lies in the interior of the set Θ({x′j}), while if = 0 then µ̂ lies on the boundary of the set Θ({x′j})). Moreover, it is easy to see that in the case µ̂>xj < , the resulting x′j is the Euclidean projection of xj onto the hyperplane µ̂
>xj = . The proposed Requirement correction admits thus the analytical form:
∀j ∈ [J ] : x′j = xj − min{0, µ̂>xj − }
‖µ̂‖22 µ̂.

heading: 3.2 MAHALANOBIS IMPROVEMENT CORRECTION
text: Given a plan {xj} such that µ̂ ∈ Θ({xj}) and an integer K between 1 and J , the Mahalanobis Improvement correction aims to modify K out of J plans to obtain the corrected plan {x′j}. The goal of this correction is to increase the lower bound value L? associated with the plan {x′j}, while at the same time keeping the amount of modification as small as possible. To attain this goal, we first describe the geometric intuition behind the lower bound L? in (2), and then leverage this intuition to generate the correction.
Geometric intuition. We first analyze the distribution of the random vector θ̃ that attains the validity lower bound L?. To simplify the exposition, we assume that λ? > 0 and define λ?0 = 1− ∑ j λ ? j . Following the same argument as in Vandenberghe et al. (2007, §2.2), the distribution of θ̃ can be constructed as a mixture of J + 1 random vectors θ̃j satisfying:
∀j = 0, . . . , J : θ̃ = θ̃j with probability λ?j ,
where for each j = 1, . . . , J , we have E[θ̃j ] = z?j /λ ? j
and θ̃0 follows a properly chosen distribution. By the validity of z?, we can verify that the location z?j /λ ? j lies on the hyperplane {θ : x>j θ = 0}. Thus, we can think of λ?j as the marginal increase in the lower bound L
? if we slightly perturb xj so that the point z?j /λ ? j lies inside the set of favorable parameters. This observation underlies the Mahalanobis correction which we describe next.
Correction procedure. If we can adjust K out of J counterfactuals to improve the validity, then it is reasonable to modify the K counterfactuals associated with the K largest values of λ?j , where λ ? j is the optimal value of the variable λj in problem (2). Without any loss of generality, assume that λ?j have decreasing values, and in this case, our correction procedure will modify the counterfactuals xj for j = 1, . . . ,K. Further, to correct each individual action, we find x′j in a ∆-neighborhood of xj such that the Mahalanobis distance from µ̂ to the hyperplane {θ : θ>x′j = 0} is maximized, where the Mahalanobis distance is computed with the nominal covariance matrix Σ̂. This is equivalent to solving a max-min problem
x′j = arg max minθ:θ>x=0
√ (θ − µ̂)>Σ̂−1(θ − µ̂)
s. t. x ∈ Rd, ‖x− xj‖2 ≤ ∆. (3)
The next result indicates that x′j can be found by solving a conic optimization problem.
Theorem 3.1 (Mahalanobis Improvement correction). The Mahalanobis correction of xj is x′j = v?/t?, where (v?, t?) is the optimal solution of the following conic optimization problem
min { v>Σ̂v : v ∈ Rd, t ∈ R+, ‖v − txj‖2 ≤ ∆t, v>µ̂ = 1 } .
We have specifically modified x′j in (3) with respect to the nominal mean vector and covariance matrix (µ̂, Σ̂) of the random vector θ̃. Alternatively, we can also use (µ?,Σ?), where (µ?,Σ?) is the optimal solution in the variable (µ,Σ) of (2) to form the optimization problem. Theorem 3.1 holds with the corresponding parameters (µ?,Σ?). Similarly, equation (4) recovers the Euclidean projection if we use an identity matrix for weighting. The conic optimization problem in Theorem 3.1 can be solved using standard off-the-shelf solvers such as Mosek (MOSEK ApS, 2019).

heading: 4 COUNTERFACTUAL PLAN CONSTRUCTION UNDER AMBIGUITY
text: We propose in this section the COunterfactual Plan under Ambiguity (COPA) framework to devise a counterfactual plan that has high validity under random future model parameters. Given an input instance x0, COPA builds a plan {xj} of J ≥ 1 counterfactuals that balances between competing objectives including proximity, diversity and validity. We next describe each cost component.
Proximity. It is reasonable to ask that each counterfactual xj should be close to the input x0 so that xj is actionable. We suppose that the distance between x0 and xj can be measured using a function c : Rd × Rd → R. In general, the cost c is used to capture the ease of adopting the changes for a specific variable (e.g., one could barely change their height or race). The proximity of a plan {xj} is
simply the average distance from x0 to each counterfactual in the plan. More specifically, we have
Proximity({xj}, x0) , 1
J J∑ j=1 c(xj , x0). (4)
Diversity. We measure the diversity of a plan using the determinant point process (Kulesza, 2012) similar to the approach in Mothilal et al. (2020). The diversity is given by:
Diversity({xj}) , det(K), where Ki,j = (1 + c(xi, xj))−1 ∀1 ≤ i, j ≤ J. (5) Then, a plan with a larger value Diversity({xj}) is more diverse.
Validity. Given the moment information (µ̂, Σ̂), one potential approach to compute the validity of a plan {xj} is to compute the value L? in (2). However, for large covariate dimension d or high number of counterfactual J , the semidefinite program (2) becomes time-consuming to solve and is not practical. This entails us to derive the a computationally efficient proxy for the validity of {xj}. Towards this goal, we use the volume of the maximum-volume ellipsoid with center µ̂ and covariance Σ̂ that can be inscribed in Θ({xj}). Following Boyd & Vandenberghe (2004, §8.4.2), an ellipsoid with center µ̂, covariance matrix Σ̂ and radius r can be written in the parametric form as E(µ̂,Σ̂)(r) = {Σ̂ 1 2u+ µ̂ : ‖u‖2 ≤ r}. The validity of the plan {xj} is thus defined as
Validity({xj}, µ̂, Σ̂) , max{r : r ≥ 0, E(µ̂,Σ̂)(r) ⊆ Θ({xj})}. The next result asserts that the above validity measure can be re-expressed in closed form, which justifies its computational efficiency.
Lemma 4.1 (Validity value). If µ̂ ∈ Θ({xj}), then Validity({xj}, µ̂, Σ̂) = minj µ̂>xj/‖Σ̂ 1 2xj‖2.
Lemma 4.1 and the analysis in Proposition 2.4 also suggest that the counterfactual plan should satisfy µ̂ ∈ Θ({xj}) so as to improve the validity. Similar to Section 3.1, we will impose the constraints that µ̂>xj ≥ ∀j for some margin ≥ 0 for validity purposes. COPA framework. Our COPA framework finds the counterfactual plan that minimizes the weighted sum of the proximity, the diversity and the validity measure. More precisely, the COPA counterfactual plan is the minimizer of
min x1,...,xJ Proximity({xj}, x0) + λ1Validity({xj}, µ̂, Σ̂)− λ2Diversity({xj}) s. t. µ̂>xj ≥ ∀j
(6)
for some non-negative parameters λ1 and λ2. The COPA problem (6) can be solved efficiently under mild conditions using a projected (sub)gradient descent algorithm.
A projected gradient descent algorithm can be used to solve the COPA problem (6). The gradient of the objective function of (6) can be computed using auto-differentiation. We now discuss further the projection operator. Let X , {x ∈ Rd : µ̂>x ≥ }, then the feasible set of the COPA problem (6) is a product space X J . The projection operator ProjXJ on the product set X J is decomposable into simpler projections onto individual set X as ProjXJ ({x′j}) = {ProjX (x′1), . . . ,ProjX (x′J)}, where each individual projection is
ProjX (x ′ j) = arg min{‖x− x′j‖2 : µ̂>x ≥ } = x′j −min{0, µ̂>x′j − }µ̂/‖µ̂‖22.
Note that the second equality above follows from the analytical formula for the Euclidean projection onto a half-space, which was previously used in Section 3.1.

heading: 5 NUMERICAL EXPERIMENTS
text: In this section, we evaluate the correctness of our validity bounds and the performance of our corrections and our COPA framework on both synthetic and real-world datasets. Our baseline for comparison is the counterfactual plan constructed from the state-of-the-art DiCE framework (Mothilal et al., 2020). Throughout the experiments, we set the number of counterfactuals to J = 5. For DiCE, we use the default parameters recommended in the DiCE source code. The Mahalanobis correction will use the counterfactual plan obtained by DiCE method with K = 3 and the perturbation limit ∆ is 0.1. In our COPA framework, we use Adam optimizer to implement Projected Gradient Descent and `2-distance to compute perturbation cost between inputs.

heading: 5.1 SYNTHETIC DATASET
text: We first generate 1000 samples with two-dimensional features from two Gaussian distributions N (µ0,Σ0) and N (µ1,Σ1) to create a synthetic dataset. Each instance is labelled as 0 or 1 corresponding to the distribution that generated it. For the Gaussian distributions, we use similar parameters as in Upadhyay et al. (2021). where µ0 = [−2,−2]>, µ1 = [2, 2]>, Σ0 = Σ1 = 0.5I with I being the identity matrix. This dataset is then used to train a logistic classifier with the present parameter θ0. This classifier Cθ0 is fixed for the experiments that follow. The impact of Gelbrich radius on the validity. Given a counterfactual plan generated by DiCE on the classifier Cθ0 , we consider two scenarios: µ̂ ∈ Θ({xj}) and µ̂ /∈ Θ({xj}). We choose µ̂ = θ0 for the case µ̂ ∈ Θ({xj}) and µ̂ = −θ0, otherwise. We also set Σ̂ = 0.5I . We then compute the lower and upper validity bound of this plan with respect to different Gelbrich bounds ρ ∈ [0, 1]. To evaluate the empirical validity of this plan, we simulate 1000 futures for θ̃. For each future, we generate µg and Σg randomly so that G((µ̂, Σ̂), (µg,Σg)) ≤ 1, and then we sample 106 values of θ̃ ∼ Ng(µg,Σg). The empirical validity of the plan for each future is the fraction of parameter samples from the future that the prescribed plan is valid. We plot the 1000 empirical validity of the plan in Figure 3. This result is consistent with our guarantees that the validity is between the two bounds. We also observe that increasing ρ loosens the validity bounds.
The impact of degree of distribution shift on validity of a plan. We explore the case µ̂ ∈ Θ({xj}), where µ̂ = θ0 and Σ̂ = 0.5I , to assess the impact of distribution shift to three algorithms DiCE, MahalanobisCrr, and COPA. In this experiment, we run our COPA framework with λ1 = 2.0, λ2 = 200.0. To assess the performance of three algorithms, we parameterize the ground truth distribution of the future parameters θ̃ ∼ N (µg,Σg) as follows: µg = µ̂+α[0,−1, 0]>,Σg = (1 + β)I . Here, we simulate three types of distributional shift of the parameters θ̃: (1) mean shift (α ∈ [0, 1], β = 0), (2) covariance shift (α = 0, β ∈ [0, 3]), and (3) mean and covariance shift ((α, β) ∈ [0, 1]× [0, 3]). For each shift’s type, we generate 100 counterfactual plans corresponding to 100 original inputs x0 and compute the empirical validity as previously described. The average and confidence range of the empirical validity are plotted in Figure 4. This result shows the tendency of decreasing validity measure of all algorithms when increasing the Gelbrich distance between estimate and ground truth distribution. However, COPA shows stability and robustness for
all shift types. The validity of DiCE deteriorates when the ground truth distribution is far from θ0. Meanwhile, MahalanobisCrr increases the robustness of the plans obtained by DiCE significantly.

heading: 5.2 REAL-WORLD DATASETS
text: In this experiment, we evaluate the robustness of the counterfactual plans obtained by three frameworks on the real datasets. We use three real-world datasets: German Credit (Dua & Graff, 2017; Groemping, 2019), Small Bussiness Administration (SBA) (Li et al., 2018), and Student performance (Cortez & Silva, 2008). Each dataset contains two datasets (the present data -D1 and the shifted data D2). The shifted dataset D2 could capture the correction shift (German credit), the temporal shift (SBA), or the geospatial shift (Student). More details for each dataset are provided in Appendix.
Experimental settings. For each present datasetD1, we train a logistic classifier Cθ0 with parameter θ0 on 80% of instances of the dataset and fix this classifier to construct counterfactual plans in whole experiment. We generate 100 counterfactual plans for 100 original inputs and report the average of values of our evaluation metrics. To estimate µ̂ and Σ̂, we train 1000 different classifiers from the present dataset D1 (each is trained on a random set containing 50% instances of D1), then use the empirical mean and covariance matrix of the parameter. We set Gelbrich radius ρ = 0.01.
Metrics. To compute the empirical validity in the shift dataset D2, we sample 50% instances of D2 1000 times to train 1000 different logistic classifiers. We then report the empirical validity of a plan as the fraction of the classifiers with respect to which the plan is valid. We also use the lower validity bound as a metric for evaluating the robustness of a plan. We use the formula in (4) and (5) to measure the proximity and diversity of a counterfactual plan.
Results. The results in Table 1 show that our COPA framework achieves highest empirical validity, L∗, and diversity (especially when increasing λ2) in all evaluated datasets. Comparing DiCE and Mahalanobis correction, we can observe that the trade-off of proximity and diversity of Mahalanobis correction is relatively small as compared to its improvement in terms of validity.

heading: 6 CONCLUSION
text: This paper studies the problem of generating counterfactual plans under the distributional shift of the classifier’s parameters given the fact that the classification model is usually updated upon the arrival of new data. We propose an uncertainty quantification tool to compute the bounds of the probability of validity for a given counterfactual plan, subject to uncertain model parameters. Further, we introduce a correction tool to increase the validity of the given plan. We also propose a COPA framework to construct a counterfactual plan by taking the model uncertainty into consideration. The experiments demonstrate the efficiency of our methods on both synthetic and real-world datasets. Further extensions, notably to incorporate nonlinearities, are presented in the appendix.

heading: A PROOFS

heading: A.1 PROOFS OF SECTION 2
text: Proof of Theorem 2.2. For any (µ,Σ) ∈ Rd × Sd+, let
P(µ,Σ) , {Q : EQ[θ̃] = µ, EQ[θ̃θ̃>] = µµ> + Σ}
denote the set of probability measures under which the random vector θ̃ has mean µ and covariance matrix Σ. The infimum probability can be decomposed as
inf Q∈B Q(θ̃ ∈ Θ◦({xj})) = inf (µ,Σ)∈U inf Q∈P(µ,Σ) Q(θ̃ ∈ Θ◦({xj}))
=  inf (µ,Σ)∈U inf 1− ∑ j∈[J] λj s. t. λj ∈ R, zj ∈ Rd, Zj ∈ Sd ∀j ∈ [J ] −x>j zj ≥ 0, [ Zj zj z>j λj ] 0 ∀j ∈ [J ]∑
j∈[J]
[ Zj zj z>j λj ] [ Σ + µµ> µ µ> 1 ] ,
where the second equality follows from Vandenberghe et al. (2007, §2). By Malagò et al. (2018, Proposition 2), we have
G2((µ,Σ), (µ̂, Σ̂)) =  min C∈Rd×d ‖µ− µ̂‖2 + Tr [ Σ + Σ̂− 2C ] s. t. [ Σ C
C> Σ̂
] 0.
=  min C∈Rd×d ‖µ̂‖2 − 2µ̂>µ+ Tr [ Σ + µµ> + Σ̂− 2C ] s. t. [ Σ C
C> Σ̂
] 0.
Hence, by combing two infimum operators, we have
inf Q∈B
Q(θ̃ ∈ Θ◦({xj})) =  inf 1− ∑ j∈[J] λj s. t. µ ∈ Rd, Σ ∈ Sd+, C ∈ Rd×d λj ∈ R, zj ∈ Rd, Zj ∈ Sd ∀j ∈ [J ] −x>j zj ≥ 0, [ Zj zj z>j λj ] 0 ∀j ∈ [J ]∑ j∈[J] [ Zj zj z>j λj ] [ Σ + µµ> µ µ> 1 ] ‖µ̂‖2 − 2µ̂>µ+ Tr [ Σ + µµ> + Σ̂− 2C ] ≤ ρ2, [ Σ C
C> Σ̂
] 0.
In the last step, we add an auxiliary variable M ∈ Sd+ with the constraint M = Σ + µµ>. Note that this constraint can be replaced by M Σ + µµ> without affecting the optimal value of the optimization problem. Using the Schur complement, this constraint is equivalent to[
M − Σ µ µ> 1
] 0.
This completes the proof.
Proof of Theorem 2.3. Let 1Θ(θ) be the indicator function of the set Θ({xj}), that is,
1Θ(θ) = { 1 if θ ∈ Θ({xj}), 0 otherwise.
By defining the loss function `(θ) = 1Θ(θ) and let Z be the convex feasible set defined by Z , { z0 ∈ R, z ∈ Rd, Z ∈ Sd : z0 + 2z>θ + 〈 Z, θθ> 〉 ≥ 1Θ(θ) ∀θ ∈ Rd } .
Notice that Z is a closed and convex set because it is an intersection of uncountably many closed and covex sets. Denote the following set V of mean - second moment matrices that is induced by U by V , {(µ,M) ∈ Rd × Sd+ : ∃(µ,Σ) ∈ U such that (µ,M) = (µ,Σ + µµ>)}. The support function δ∗V of the set V is defined as
δ∗V(z, Z) = sup{z>µ+ Tr [ ZM ] : (µ,M) ∈ V}.
Using these notations, we now have
sup Q∈B Q(θ ∈ Θ({xj})) = sup (µ,Σ)∈U sup Q∈P(µ,Σ) EQ[1Θ(θ̃)] (7a)
≤ sup (µ,Σ)∈U inf (z0,z,Z)∈Z z0 + 2µ >z + Tr
[ (Σ + µµ>)Z ] (7b)
= sup (µ,M)∈V inf (z0,z,Z)∈Z
z0 + 2µ >z + Tr [ MZ ] = inf
(z0,z,Z)∈Z sup (µ,M)∈V z0 + 2µ
>z + Tr [ MZ ] (7c)
= inf (z0,z,Z)∈Z
z0 + δ ∗ V(2z, Z),
where equality (7a) is from the two layer decomposition of the ambiguity set B, and inequality (7b) is from the Isii’s duality result Isii (1960). Equality (7c) follows from the Sion’s minimax theorem Sion (1958) which holds because the objective function is linear in each variable and because V is compact by the compactness of U (Nguyen et al., 2021, Lemma A.6). We thus have
sup Q∈B
Q(θ ∈ Θ({xj})) =  inf z0 + γ(ρ 2 − ‖µ̂‖22 − Tr [ Σ̂ ] ) + q + Tr [ Q ] s. t. γ ∈ R+, z0 ∈ R, z ∈ Rd, Z ∈ Sd, q ∈ R+, Q ∈ Sd+[ γI − Z γΣ̂ 12 γΣ̂ 1 2 Q ] 0, [ γI − Z γµ̂+ z γµ̂> + z> q ] 0
z0 + 2z >θ + 〈 Z, θθ> 〉 ≥ 1Θ(θ) ∀θ ∈ Rd,
where the equality follows by substituting the support function of V in Kuhn et al. (2019, Lemma 2). Consider now the last constraint of the above optimization problem, it is easy to see that it is equivalent to {
z0 + 2z >θ + 〈 Z, θθ> 〉 ≥ 0 ∀θ ∈ Rd, z0 + 2z >θ + 〈 Z, θθ> 〉 ≥ 1 ∀θ ∈ Θ({xj}).
The first semi-infinite constraint is equivalent to the semidefinite constraints Z 0, [ Z z
z> z0
] 0.
A sufficient condition for the second semi-infinite constraint is that ∃λ ∈ RJ+ : [ Z z
z> z0 − 1 ] ∑ j∈[J] λj [ 0 12xj 1 2x > j 0 ] ,
which holds thanks to the S-lemma Pólik & Terlaky (2007). Adding these above constraints into the optimization problem leads to the desired upper bound. This completes the proof.
The proof of Proposition 2.4 relies on the following result on the multivariate Chebyshev inequalities, which can be found in Marshall & Olkin (1960) and Bertsimas & Popescu (2005). Theorem A.1 (Multivariate Chebyshev inequality). Let S be a convex set, then
sup Q∼(µ,Σ) Q(θ̃ ∈ S) = 1 1 + κ , κ = inf θ∈S (θ − µ)>Σ−1(θ − µ).
We are now ready to prove Proposition 2.4.
Proof of Proposition 2.4. If µ̂ ∈ Θ({xj}), then it is clear that
inf θ∈Θ({xj})
(θ − µ̂)>Σ̂−1(θ − µ̂) = 0,
thus we have U? ≥ sup Q∼(µ̂,Σ̂) Q(θ̃ ∈ Θ({xj})) = 1 by Theorem A.1. This leads to U? = 1.
Consider the case when µ̂ 6∈ Θ({xj}). By the hyperplane separation theorem, there exists a vector x̄ ∈ Rd such that x̄>θ ≥ 0 for all θ ∈ Θ({xj}) and x̄>µ̂ < 0. Let T , {θ : x̄>θ ≥ 0}, then it is trivial that Θ({xj}) ⊆ T and that T is a convex set. We now have
inf Q∈B Q(θ̃ ∈ Θ({xj})) = 1− sup Q∈B Q(θ̃ 6∈ Θ({xj})) ≤ 1− sup Q∈B Q(θ̃ 6∈ T) = 1− 1 = 0,
where the penultimate equality follows from Theorem A.1. This leads to L? = 0

heading: A.2 PROOFS OF SECTION 3
text: Proof of Theorem 3.1. Notice that the optimal solution in x should satisfy x>µ̂ > 0. Fix any value of x 6= 0. Consider first the inner minimization problem of (3), and associate with the equality constraint a Lagrangian dual variable ζ ∈ R, we have
min θ:θ>x=0 (θ − µ̂)>Σ̂−1(θ − µ̂) = min θ∈Rd max ζ∈R (θ − µ̂)>Σ̂−1(θ − µ̂) + 2ζθ>x
= max ζ∈R min θ∈Rd
(θ − µ̂)>Σ̂−1(θ − µ̂) + 2ζθ>x
= max ζ∈R
− ζ2x>Σ̂x+ 2ζµ̂>x
= (x>µ̂)2
x>Σ̂x ,
where the second equality follows from convex duality result. The third equality follows from the fact that for every value of ζ, the optimal solution in the variable θ is
θ?(ζ) = µ̂− ζΣ̂x.
Moreover, the last equality follows from the optimality condition in ζ which gives ζ? = µ̂>x/(x>Σ̂x). Because the optimal solution in x should satisfy x>µ̂ > 0, problem (3) is hence equivalent to
max x>µ̂√ x>Σ̂x s. t. ‖x− xk‖ ≤ ∆.
Adding now two auxiliary variables t ∈ R+ and v ∈ Rd with the constraints:
1
x>µ̂ = t, v = tx,
the claim in the statement of the theorem now follows by a simple substitution to get
max 1√ v>Σ̂v s. t. v ∈ Rd, t ∈ R+, ‖v − txk‖2 ≤ ∆t, v>µ̂ = 1.
Swapping the maximum operator to a minimum operator completes the proof.

heading: A.3 PROOFS OF SECTION 4
text: Proof of Lemma 4.1. From the definition of the set Θ({xj}), we have
max{r : Er ⊆ Θ({xj})} =
{ max r
s. t. sup ‖u‖2≤r
−(Σ̂ 12u+ µ̂)>xj ≤ 0 ∀j
=
{ max r
s. t. −µ̂>xj + sup ‖u‖2≤r −x>j Σ̂ 1 2u ≤ 0 ∀j
=
{ max r
s. t. −µ̂>xj + r‖Σ̂ 1 2xj‖2 ≤ 0 ∀j,
where the last equality follows from the dual norm property. The proof now follows by finding the maximum value of r so that the problem is feasible.

heading: B EXPERIMENTS

heading: B.1 EXPERIMENTAL DETAIL
text: Real-world datasets Here, we provide more detail about three real-world datasets we used. Source code can be found at https://anonymous.4open.science/r/COPA-ICLR.
i German Credit (Dua & Graff, 2017). The dataset contains the information (e.g. age, gender, financial status,...) of 1000 customers who took a loan from a bank. The classification task is to determine the risk (good or bad) of an individual. There is another version of this dataset regarding to corrections of coding error (Groemping, 2019). We use the corrected version of this dataset as a shifted data to capture correction shift. The features we used in this dataset include ‘duration’, ‘amount’, ‘personal status sex’, and ‘age’.
ii Small Bussiness Administration (SBA) (Li et al., 2018). This data includes 2,102 observations with historical data of small business loan approvals from 1987 to 2014. We divide this dataset into two datasets (one is instances from 1989 - 2006 and one is instances from 2006 - 2014) to capture temporal shift. We use the following features: selected, ‘Term’, ‘NoEmp’, ‘CreateJob’, ‘RetainedJob’, ‘UrbanRural’, ‘ChgOffPrinGr’, ‘GrAppv’, ‘SBA Appv’, ‘New’, ‘RealEstate’, ‘Portion’, ‘Recession’.
iii Student performance (Cortez & Silva, 2008). This data includes the performance records of 649 students in two schools: Gabriel Pereira (GP) and Mousinho da Silveira (MS). The classification task is to determine if their final score is above average or not. We split this dataset into two sets in two schools to capture geospatial shift. The features we used are: ‘age’, ‘Medu’, ‘Fedu’, ‘studytime’, ‘famsup’, ‘higher’, ‘internet’, ‘romantic’, ‘freetime’, ‘goout’, ‘health’, ‘absences’, ‘G1’, ‘G2’.
Classifier Throughout this paper, we use a Logistic Regression for linear classifier and a threelayer MLP with 20, 50, 20 nodes and ReLU activation in each consecutive layer as nonlinear classifier. We use one-hot encoding for categorical features in the datasets to convert it to a vector of [0, 1]. We use min-max normalization to scale the numerical features to [0, 1]. We report the performance of the classifiers in three real-world datasets in Table 2

heading: B.2 ADDITIONAL EXPERIMENTS
text: The impact of degree of distribution shift on validity of a plan. We provide an additional experiment in different covariance shift Σg = (1 + β)A,A 0. In this experiment, we choose A as:
A = ( 1 −1 0 −1 1 1 0 1 1 ) .
The matrix A introduces both positive and negative correlations between classifier’s parameters. Other settings are set the same as the experiment in Section 5.1.
Mahalanobis correction on real-world dataset. In this experiment, we evaluate the Mahalanobis correction on different number of corrections K and different perturbation limit ∆. We set ρ = 0.01, = 0.1,K ∈ {0, . . . , J}, J = 5,∆ ∈ [0.05, 0.35]. (µ̂, Σ̂) is estimated using similar manner as in Section 5.2 of the main paper. The results in shown in Figure 6.
Counterfactual explanations for real-world datasets. To illustrate the use case of the counterfactual explanations, we provide some examples on the German dataset with J = 3 (Table 3) in which we consider the “personal status and sex” feature as immutable. Here, we can observe that three algorithms could provide diverse sets of counterfactuals that the users may prefer. However, by providing better empirical validity, the plans generated by MahalanobisCrr and COPA are more robust with distribution shift than DiCE (generated without considering the shift).

heading: C EXTENSION TO NONLINEAR CLASSIFIERS
text: In the main paper, our analysis is based on the linearity in both features and model’s parameters. We now discuss two extensions of our COPA framework to the nonlinear settings.

heading: C.1 NONLINEARITY IN INPUT FEATURES
text: This section extends to any linear classifier Cθ(x) = 1 if θ>φ(x) ≥ 0, and 0 otherwise, where φ : X → Rd is a (possibly nonlinear) feature mapping that maps input features to a latent representation in a covariate space Rd. Note that our bounds in Section 2 still hold in latent space Rd: for a concrete example, Theorem 2.2 holds with xj being replaced by φ(xj).
The COPA framework is also extendable to incorporate the feature map φ. Assuming that φ is differentiable, the COPA framework solves the following optimization problem:
min x1,...,xJ Proximity({xj}, x0) + λ1Validity({φ(xj)}, µ̂, Σ̂)− λ2Diversity({xj}). s. t. µ̂>xj ≥ ∀j
(8)
The proximity and diversity are measured in the input space and the validity term is now measured in latent space space instead. This optimization problem can be solved efficiently by projected gradient descent similar to Section 4.

heading: C.2 NONLINEARITY IN MODEL’S PARAMETERS
text: Similar to the prior works (Ustun et al., 2019; Rawal & Lakkaraju, 2020; Upadhyay et al., 2021), our work can adapt to nonlinear classifiers Cnl using a local surrogate models such as LIME (Ribeiro et al., 2016). LIME (Ribeiro et al., 2016) is a popular technique for explaining predictions of blackbox machine learning models. The main idea of LIME is to train a local surrogate model Cx0θ on a perturbed samples around a given input instance x0 to approximate the local decision boundary of the black-box models. We thus model the uncertainty of parameters θ in the surrogate model Cx0θ for x0 instead of the parameters of Cnl. For the experiment, we first generate a local linear model Cx0θ using LIME method with 5000 perturbed samples. We then choose (µ̂, Σ̂) = (θ, 0.05I), where I is identity matrix, to model the distributional uncertainty of the parameters. Similar to Section 5.2, we set Gelbrich radius ρ is to 0.01, J = 5, K = 3.
We report the performance of three algorithms on the MLP classifier in the real-world datasets in Table 4. The result is promising since the proposed COPA can increase the empirical validity significantly. However, the infidelity of LIME could lead to invalid counterfactual explanations, represented by a lower current validity value. The low current validity is also observed in the literature, see Upadhyay et al. (2021). For further investigation, one can use another local surrogate model that provides a better approximation of the decision boundary (e.g., BayLIME (Zhao et al., 2020)). Another direction is to use a mixture linear regression model to approximate the decision
boundary as in Guo et al. (2018). However, advocating for the mixture of linear models requires further analysis.

Review summaries:
This paper studies the effect of uncertainly in the parameters of an ML model used for consequential decision making on the validity, proximity, and diversity of sets of counterfactual explanations (CFE), namely counterfactual plans. In particular, the paper assumes the model parameters are sampled from a distribution would known mean and covariance, and studies the aforementioned metrics of counterfactual plans when the parameter distribution changes in a bounded manner (according to the Gelbrich distance). Beyond a diagnostic tool (providing bounds on how much the parametere distribution can change), the paper presents a method for counterfactual plan correction after having generated a plan (ex-post), and an approach to generate robust counterfactual plans before the fact (ex-ante). Experiments are presented for the proposed methods.

The paper considers the problem of how to provide a collection of counterfactual explanations for a binary linear classifier such that an informed choice can be made on how best to actualise a new input based on individual preferences such that a different classification is likely to be made, even under potentially changing model parameters. It introduces the Counterfactual Plan under Ambiguity (COPA) framework, consisting of a probabilistic validity assessment of a predetermined plan under a given parameter distribution of the classifier, followed by a correction of the predetermined plan to improve validity under that distribution. These are combined in the COPA framework to produce a plan that optimises a weighted combination of the desiderata ‘proximity’, ‘diversity’ and ‘validity’. Experiments show the resulting output satisfies the derived upper and lower bounds on the validity, with good performance on the desiderata.