Optimal Counterfactual Search in Tree Ensembles: A Study Across Modeling and Solution Paradigms
arXiv:2605.06561v1 [cs.LG] 7 May 2026
Awa Khouna Polytechnique Montréal [email protected]
Youssouf Emine Ivey Business School [email protected]
Julien Ferry Polytechnique Montréal [email protected]
Thibaut Vidal Polytechnique Montréal [email protected]
Abstract Trust in counterfactual explanations depends critically on whether their recommended changes are truly minimal: suboptimal explanations may vastly overshoot the actual changes needed to alter a decision, and heuristic errors can affect individuals unevenly, giving some users relevant recourse while assigning others unnecessarily costly recommendations. Consequently, we study the problem of computing optimal counterfactual explanations for tree ensembles under plausibility and actionability constraints. This is a combinatorial problem: for a fixed model, counterfactual search boils down to selecting consistent branching decisions and threshold-defined regions under a distance objective. We exploit this structure through C PCF, a constraint programming (CP) formulation in which numerical features are encoded as interval domains induced by split thresholds, while discrete features retain native finite-domain representations. This yields a compact finitedomain formulation that supports multiple distance objectives without continuous split-boundary search. We then place C PCF in a broader comparison across mathematical programming paradigms: we extend a maximum Boolean satisfiability (MaxSAT) formulation, originally designed for hard-voting random forests, to softvoting ensembles, and compare against the current state-of-the-art mixed-integer linear programming (MILP) optimal approach. Across ten datasets and three types of tree ensembles, we analyze scalability, anytime performance, and sensitivity to distance metrics. We observe that CP achieves the best overall performance. More importantly, our results identify regimes in which the specific strengths of each paradigm make it best suited: CP is most versatile overall, MaxSAT handles hardvoting ensembles particularly well, and MILP remains competitive in amortized inference settings with a moderate number of split levels.
1
Introduction
Counterfactual explanations have become a widely used form of local explanation for machine learning models [Karimi et al., 2022, Verma et al., 2024]. Given an input instance, a counterfactual specifies how its features should change to obtain a desired prediction, thereby providing interpretable recourse in domains such as credit, healthcare, and public decision-making. However, generating high-quality counterfactuals can be challenging for many models of interest: their decision boundaries may be non-convex and non-differentiable, their feature spaces can combine different feature types, and realistic explanations must satisfy plausibility and feasibility constraints [Joshi et al., 2019, Poyiadzi et al., 2020]. As a result, many counterfactual approaches rely on heuristics or generative models [Looveren and Klaise, 2020, Nemirovsky et al., 2022], without guarantees on the quality of the returned solution. This limitation has severe consequences. Since counterfactual explanations are Preprint.
Soft
Performance metric 67.2%
only solver time total tim
Mean #split levels per feature 34.3%
≤ 7.42
> 7. 42
e
Voting 100.0% Hard
<2
#Ordinal features 32.8%
≥2
MILP 13.5%
CPCF 58.8% #Estimators 6.8%
> 100
≤ 100
MAXSAT 27.7%
Figure 1: Decision diagram summarizing the best-performing optimal counterfactual formulation for tree ensembles, for each setting (see Appendix A for construction details). Percentages indicate the share of experimental configurations reaching each node. The proposed C PCF performs best overall, M AX SAT Raevskaya and Lehtonen [2025] excels on moderate-size hard-voting random forests or datasets with many ordinal features, while MILP [Parmentier and Vidal, 2021] is competitive in soft-voting settings with moderately many split levels per feature, especially if only solver time matters (e.g., if optimization model-building time is amortized over many queries).
often interpreted as minimal changes required to alter a decision, suboptimal outputs can overstate the required effort and suggest unnecessarily costly or impractical forms of recourse [Parmentier and Vidal, 2021]. Errors may also be unevenly distributed across the input space: some individuals may receive near-minimal recommendations, while others are assigned substantially more burdensome changes than necessary. This motivates the study of optimal counterfactual explanations, which provably minimize a prescribed notion of change under prediction and feasibility constraints. Tree ensembles are a central model class in this context. They remain among the strongest predictors for many tabular-data tasks that require explanations and recourse, yet their non-convex, nondifferentiable decision functions pose substantial challenges for counterfactual generation. At the same time, their structure exposes a combinatorial problem that can be leveraged by optimization algorithms: for any fixed ensemble, counterfactual search reduces to selecting consistent split decisions and leaf regions. Several heuristics exploiting this structure have been proposed. Early methods such as actionable feature tweaking and related recourse procedures search for favorable leaf regions or candidate modifications without fully solving the underlying optimization problem [Tolomei et al., 2017, Cui et al., 2015, Fernández et al., 2020]. Other approaches prioritize speed, diversity, or realism, often through approximate search or data-driven restrictions of the feasible region [Carreira-Perpinan and Hada, 2023, Zhang and Zhong, 2025]. Differentiable approximations of the decision boundaries of the ensembles have also been explored [Lucic et al., 2022]. However, these heuristic algorithms regularly lead to vastly suboptimal explanations, e.g., by factors of over 10× as seen in [Parmentier and Vidal, 2021, Khouna et al., 2026]. To yield explanations with guarantees, recent work has explored mathematical programming formulations, based on mixed-integer linear programming (MILP) [Parmentier and Vidal, 2021], satisfiability-based methods [Karimi et al., 2020], and weighted maximum Boolean satisfiability (MaxSAT) [Raevskaya and Lehtonen, 2025]. These approaches differ in how they represent and explore the search space, and on which mathematical programming paradigm they are grounded to prune the search space and ensure optimality. In this work, we first introduce C PCF, a constraint programming formulation for optimal counterfactual search in tree ensembles. Our starting point is that the finite-domain primitives of CP align naturally with the split-induced regions of tree ensembles: numerical features can be represented through threshold intervals, discrete features retain their finite-domain structure, and tree paths are enforced through logical constraints. This yields a compact formulation that directly exploits the problem’s combinatorial structure and supports multiple distance objectives, as well as plausibility and actionability requirements. Compared to MILP formulations, C PCF avoids reasoning over continuous split boundaries during search and does not require linearization; compared to Boolean encodings, it retains a compact representation of numerical features. We then place C PCF in a broader empirical comparison of exact mathematical programming paradigms for counterfactual search. Across ten datasets, multiple ensemble types, and several distance objectives, C PCF achieves the best overall performance. More importantly, our analysis shows that no single paradigm dominates across all regimes: CP is the most versatile overall, MaxSAT is particularly effective for hard-voting random forests, and MILP is competitive when model-building time can be amortized in moderate split-level regimes. This provides practical guidance for selecting an exact counterfactual solver, summarized by 2
the decision diagram in Figure 1, which relates the preferred formulation to ensemble characteristics, dataset structure, and the performance metric of interest. Our main contributions are as follows: 1. We introduce C PCF, the first constraint programming formulation for optimal counterfactual explanations in tree ensembles. It provides a compact finite-domain model of counterfactual search and supports common tree ensembles, heterogeneous feature types, distance objectives, plausibility and actionability requirements. 2. We provide a unified implementation of exact mathematical programming approaches for this problem, spanning CP, MILP, SAT, and MaxSAT. As part of this framework, we extend the existing MaxSAT formulation of Raevskaya and Lehtonen [2025] beyond hard-voting random forests to support soft-voting ensembles and boosted trees. A user-friendly library incorporating all formulations will be released upon acceptance, under an MIT license. 3. We conduct a large-scale empirical comparison across ten datasets with varied characteristics, multiple ensemble types, and several distance objectives. Beyond aggregate performance, our analysis characterizes when each paradigm is preferable, providing practical guidance for selecting an exact counterfactual solver.
2
Background
Let T be a tree ensemble (e.g., random forest or XGBoost), with classification function hT : X 7→ Y. For any given input sample x ∈ X , its predicted label is the class with maximum confidence score: X hT (x) = arg max sy (x) where sy (x) = by + wt st,y (x), y∈Y
t∈T
where by is a constant bias for class y, wt denotes the weight associated with tree t ∈ T , and st,y (x) denotes the score assigned by tree t to class y. Two voting schemes are commonly used. Hard voting occurs when each tree casts a vote for a single class, i.e., wt = 1 and st,y (x) ∈ {0, 1}, whereas soft voting generalizes this setting by allowing each tree to contribute class-specific confidence scores, with wt ∈ [0, 1] and st,y (x) ∈ R , for each tree t ∈ T and class y ∈ Y. Let x̂ ∈ X denote a query sample and let y ⋆ ∈ Y denote a desired target class, such that hT (x̂) ̸= y ⋆ . An optimal counterfactual explanation for x̂ is a sample x ∈ X that is classified as the target class y ⋆ , lies as close as possible to x̂ under a prescribed cost function, and satisfies optional plausibility and actionability constraints. We therefore consider the optimization problem: min cost(x, x̂)
x∈X
s.t.
sy⋆ (x) ≥ sy (x)+εc ,
x ∈ Xplausible ∩Xactionable ,
∀y ∈ Y\{y ⋆ }, (1)
where cost(x, x̂) measures the cost of changing the query instance into the counterfactual, and εc > 0 is a small margin used to enforce a strict target prediction. The set Xplausible restricts counterfactuals to regions well supported by the empirical data distribution, excluding low-density or out-of-distribution instances, while Xactionable encodes the values attainable from x̂ under structural, immutability, and monotonicity constraints. This formulation is shared by all methods considered in this paper. The differences between C PCF and the baselines therefore do not lie in the high-level optimization problem, but rather in how the feasible set, tree structure, and prediction constraints are encoded within different mathematical programming paradigms. Several mathematical programming formulations have been proposed to model and solve Problem (1). First, the seminal work of Karimi et al. [2020], known as M ACE, maps the nearest-counterfactual problem into a sequence of satisfiability problems, encoding the predictive model, the distance function, and feasibility requirements as logical formulae. Subsequently, O CEAN [Parmentier and Vidal, 2021] proposed a MILP formulation for tree ensembles, in which feature changes, tree paths, and class scores are represented through linear constraints and both continuous and binary decision variables. More recently, Raevskaya and Lehtonen [2025] leveraged a partial weighted Maximum Satisfiability (M AX SAT) formulation, where validity and feasibility conditions are encoded as hard clauses, while the counterfactual cost is represented through weighted soft clauses. This formulation is restricted to hard-voting ensembles, since incorporating numerical weights and confidence scores requires pseudo-Boolean encodings and pairwise class-score comparisons, which can be computationally costly in a purely propositional encoding. For comparison purposes, we 3
extend it to handle the general case of soft-voting. To this end, we encode each pairwise target-class constraint as a hard pseudo-Boolean inequality over leaf-selection variables, and translate these inequalities into weighted CNF using a standard pseudo-Boolean encoding [Een and Sörensson, 2006, Philipp and Steinke, 2015]. This extends the original M AX SAT formulation beyond hard-voting while preserving a purely propositional optimization model (see details in Appendix B). Our proposed C PCF formulation instead leverages constraint programming (CP) [Rossi et al., 2006], a mathematical programming paradigm that models combinatorial problems through decision variables, finite domains, and constraints. A key strength of CP lies in its ability to exploit global constraints, which provide compact, high-level representations of structured subproblems and dedicated filtering algorithms. These algorithms propagate constraints by removing inconsistent values from variable domains, thereby pruning large portions of the search space before and during the search. In contrast to MILP formulations, CP does not require all constraints to be linearized, and can instead preserve and exploit the discrete and logical structure of the counterfactual search problem in tree ensembles. Another distinctive feature of C PCF is that numerical and ordinal features are represented through interval domains induced directly by the split thresholds appearing in the ensemble.
3
Methodology
We now describe how C PCF encodes and solves Problem (1). For each tree t ∈ T , we use Lt to denote its set of leaves. Each leaf ℓ ∈ Lt is associated to a confidence score pt,ℓ,y for each possible class y ∈ Y, such that st,y (x) = pt,ℓ,y if example x falls into leaf ℓ of tree t. Modeling the counterfactual features. For each feature f ∈ {1, . . . , d}, we define one decision variable encoding the counterfactual’s value, according to the feature type: • Each binary feature is modeled through a binary variable xf ∈ {0, 1}. • Each categorical feature is represented through a one-hot encoding group Gk ⊆ {1, . . . , d} of binary variables xf ∈ {0, 1}. To ensure that exactly one category is selected in each group, P we enforce f ∈Gk xf = 1 via a global ExactlyOne constraint. • Numerical features are modeled through interval encoding. Let Af = {τf,1 , . . . , τf,kf } be the sorted set of split thresholds appearing in T for feature f . Together with global lower and upper bounds lb(f ) and ub(f ), these thresholds define the ordered partition If = lb(f ), τf,1 , . . . , τf,kf , ub(f ) . Instead of reasoning on the feature values directly, C PCF introduces an interval-index integer variable xf ∈ {0, . . . , kf } for each feature f , where xf = m indicates that its value belongs to the interval (If [m], If [m + 1]]1 . Note that this discretization remains exact with respect to the ensemble: two values in the same interval are indistinguishable to the model, as they induce identical tree routing decisions. Modeling the counterfactual paths. We use a binary variable zt,ℓ ∈ {0, 1} to indicate whether the counterfactual x falls into leaf ℓ ∈ Lt of tree t ∈ T . Since the counterfactual falls into exactly one leaf of each tree, we leverage global ExactlyOne constraints to impose X zt,ℓ = 1 ∀t ∈ T . (2) ℓ∈Lt
Each leaf corresponds to a conjunction of split conditions along its root-to-leaf path. We therefore enforce consistency between leaf assignments and counterfactual feature values: • For numerical features, split conditions are enforced through the corresponding interval index. If the path to leaf ℓ in tree t contains the condition xf ≤ τ , then zt,ℓ = 1 ⇒ xf ≤ mf (τ ), (3) where mf (τ ) is the index of threshold τ in the ordered partition If . Likewise, if the path contains xf > τ , then zt,ℓ = 1 ⇒ xf ≥ mf (τ ) + 1. (4) 1 For XGBoost ensembles, we use [I [m], I [m + 1)). f f
4
• For binary and one-hot encoded categorical features, path conditions are enforced directly on the corresponding binary variable. If the path to leaf ℓ requires feature f to take value v ∈ {0, 1}, then zt,ℓ = 1 ⇒ xf = v. (5) Implications (3) to (5) are implemented using efficient reification constraints. Since the split conditions within each tree induce a partition of the feature space, these implications implicitly enforce the leaf-selection constraints in (2). However, explicitly imposing them through global constraints yields stronger and faster propagation. The class score for each class y ∈ Y is computed as XX sy = by + wt pt,ℓ,y zt,ℓ ,
Enforcing target class prediction.
(6)
t∈T ℓ∈Lt
where wt is the weight of tree t, by is an optional class-dependent base score, and pt,ℓ,y denotes the P contribution of leaf ℓ in tree t to class y, so that st,y (x) = ℓ∈Lt pt,ℓ,y zt,ℓ . We then ensure that the counterfactual is predicted as the target class y ⋆ by imposing ∀y ∈ Y \ {y ⋆ }.
sy⋆ ≥ sy + εc
(7)
Objective value. As in Problem (1), the objective is to minimize the distance, or cost, between the query instance x̂ and the counterfactual x. In C PCF, this objective is represented featurewise, with per-feature actionability costs. For separable distances, C PCF optimizes min cost(x, x̂) = x,z
d X
αf δf (xf , x̂f ),
(8)
f =1
where αf is a user-specified feature-specific cost coefficient. For binary and categorical features, δf (xf , x̂f ) is defined directly from the change in feature value. For numerical features, the distance depends on the selected interval. For instance, using the Lp -type separable cost with p ≥ 0, we define if x̂f ∈ [If [m], If [m + 1]), 0, p δf (m, x̂f ) = (If [m] − x̂f ) , (9) if x̂f < If [m], p (x̂f − If [m + 1]) , if x̂f ≥ If [m + 1], which corresponds to the minimum displacement required to move the query value into the selected interval. In practice, the cost associated with each interval is retrieved from its index, leveraging efficient Element constraints. More elaborate objectives, such as nonlinear actionability costs or direction-dependent costs, can be incorporated by modifying the per-interval costs computed in (9). Actionability and plausibility. The plausibility and actionability restrictions appearing in Problem (1) can be imposed directly on the feature variables. In particular, immutable features can be fixed to their original values, directional constraints can enforce monotonic changes, and admissibledomain restrictions can be added independently of the ensemble structure. This modularity is one of the practical advantages of the CP formulation, whose large and generic catalog of constraints allows additional desiderata to be expressed efficiently. Furthermore, similar to Parmentier and Vidal [2021], plausibility can be enforced in C PCF through a user-provided isolation forest [Liu et al., 2008], a popular outlier-detection model that computes anomaly scores from the average path length of an example across a set of randomly built trees. As detailed in Appendix C, the path of the counterfactual through each isolation tree is then modeled using constraints similar to Constraints (2)–(5), and a hard constraint is imposed on the resulting anomaly score.
4
Numerical Experiments
This section evaluates C PCF alongside existing mathematical programming formulations on benchmark cases spanning heterogeneous feature types, ensemble types, and controlled variations in ensemble complexity. We assess the overall efficiency of C PCF and characterize how each paradigm behaves across tree depth, ensemble size, cost function, and plausibility constraints. 5
Table 1: Summary of the datasets used in our experiments. Columns N , O, C, and B report the number of numerical, ordinal, categorical, and binary features, respectively. Categorical features are one-hot encoded during preprocessing, resulting in a total of #X (OHE) features and #Y classes. Dataset AD: Adult BC: Breast Cancer Wisconsin CP: COMPAS CC: Credit Card GC: German Credit ON: Online News Popularity PH: Phishing SE: Seeds SP: Spambase ST: Students Performance
4.1
# samples
# features
N
O
C
B
#X (OHE)
#Y
45,222 683 5,278 29,623 1,000 39,644 11,055 210 4,601 395
11 9 5 14 9 47 30 7 57 30
2 0 0 0 1 37 0 7 57 0
3 9 2 11 5 6 8 0 0 13
4 0 0 0 3 2 0 0 0 4
2 0 3 3 0 2 22 0 0 13
41 9 5 14 19 59 30 7 57 43
2 2 2 2 2 2 2 3 2 2
Experimental Setup
Baselines. We evaluate C PCF alongside three methods spanning the different mathematical programming paradigms for optimal counterfactual generation in tree-based models: O CEAN, current state-of-the-art MILP formulation [Parmentier and Vidal, 2021]; M ACE, a satisfiability-based approach [Karimi et al., 2020]; and the M AX SAT [Raevskaya and Lehtonen, 2025] formulation. For O CEAN, we use the original formulation of Parmentier and Vidal [2021], with minor solver-parameter adjustments to preserve fidelity to the underlying scikit-learn random forest implementation. For M ACE, we set ε = 10−3 . C PCF is solved using the OR-Tools CP-SAT solver v9.14.6206 [Perron and Didier], O CEAN using Gurobi v12.0, M AX SAT using the CaDiCaL SAT solver [Biere et al., 2024], and M ACE using Z3 [De Moura and Bjørner, 2008], all through their Python bindings. Further details on numerical precision and method-specific implementation choices are provided in Appendix D. Datasets. We use ten tabular classification datasets with heterogeneous feature types, summarized in Table 1. For AD, CC, CP, GC, ON, PH, SP, and ST, we follow the benchmark introduced by Parmentier and Vidal [2021], in order to maintain direct comparability with prior exact approaches for tree-ensemble counterfactual explanations. We do not impose dataset-specific actionability restrictions, and complement this benchmark with two UCI datasets, namely BC [Wolberg, 1990] and SE [Charytanowicz et al., 2010]. BC adds a small-scale binary classification task with predominantly ordinal features, while SE introduces a multiclass setting that is absent from the original benchmark. For BC, we remove incomplete observations, retaining 683 complete samples. Models, configurations, and compute budget. For each dataset, we train tree ensemble classifiers and evaluate all counterfactual methods under a common protocol. We consider three ensemble types: hard-voting and soft-voting random forests [Breiman, 2001], as implemented in the scikit-learn library [Pedregosa et al., 2011], and XGBoost classifiers [Chen and Guestrin, 2016] from the eponymous library, which inherently use soft-voting. Since M ACE does not support boosted ensembles, it is not benchmarked on XGBoost. Unless stated otherwise, the default configuration uses 100 estimators with maximum depth 5. To isolate sources of combinatorial growth, we use two scaling protocols: depth scaling, with 100 estimators and maximum depth in {3, 4, 5, 6, 7, 8}, and ensemble scaling, with maximum depth 5 and number of estimators in {10, 20, 50, 100, 200, 500}. For each dataset and configuration, we sample 50 query instances and repeat experiments over five random seeds, yielding 250 counterfactual computations. All experiments run on homogeneous cluster nodes equipped with AMD EPYC 9654 (Zen 4) @ 2.40GHz CPUs, 64 GB RAM, and 8 threads per run, with a 900-second time limit per counterfactual search. The source code required to reproduce our experiments will be publicly released upon acceptance. Evaluation metrics. We focus on metrics that distinguish exact counterfactual methods. Since all methods solve the same optimization problem under identical prediction and feasibility constraints, standard counterfactual-quality measures such as validity, sparsity, and final distance are not informative for comparison: once optimality is certified, validity is guaranteed and the optimal objective value is the same across formulations. The central question is therefore computational: how quickly each method finds optimal counterfactuals and how they scale with the ensemble’s characteristics. We report final performance using median total generation time and optimality status. Total generation 6
Fraction proven optimal
RF soft voting
1.0 0.8 0.6 0.4 0.2 0.0
RF hard voting
XGB
CPCF OCEAN MACE MaxSAT
100
101
102
103
101 102 Time (s)
100
103
100
101
102
103
Figure 2: Cactus plots for the default ensemble configuration, aggregated over all ten datasets. Each curve reports the fraction of queries solved (out of a total of 2,500) to proven optimality (higher is better) as a function of total time, including model construction.
Mean normalized error
RF soft voting
RF hard voting
XGB
CPCF OCEAN
10 1 10 3 10 5 10 7 10 1
100
101
102
103
10 1
100
101 Time (s)
102
103
10 1
100
101
102
103
Figure 3: Anytime solution quality for C PCF and O CEAN for the default ensemble configuration, aggregated over all ten datasets and restricted to queries solved to optimality by both methods (which covers more than 99% of the queries for all three ensemble types). The curve reports the mean normalized cost of the best incumbent objective; lower is better. time includes both model-building time and solver time, where model-building time corresponds to instantiating the optimization model from the trained ensemble. This distinction matters because encoding overhead can differ substantially across paradigms; it also differs from prior work, which typically reports solver time only [Parmentier and Vidal, 2021, Raevskaya and Lehtonen, 2025]. We also evaluate anytime performance by tracking the best incumbent objective value over time, thereby measuring solution quality before optimality is certified. Note that only O CEAN and C PCF support such anytime behavior. Finally, although all main experiments use the L1 -norm, the choice of recourse cost can substantially affect the structure of the optimization problem. We therefore assess cost sensitivity by comparing the two formulations that support all considered objectives, C PCF and O CEAN, under L0 , L1 , and L2 -type costs. We exclude M AX SAT and M ACE from this analysis because they did not solve enough instances within the time limit to support a reliable comparison, and because M AX SAT does not support the L2 objective. Further details regarding the experimental setup are provided in Appendix E. 4.2
Results
We organize the empirical analysis around five findings, covering default benchmark performance, anytime behavior, scaling with ensemble complexity, sensitivity to the objective, and the effect of plausibility constraints. Result 1. The best-performing formulation depends on the setup. Figure 2 shows the fraction of counterfactuals proven optimal as a function of total running time, aggregated over datasets and random seeds, for hard-voting random forests, soft-voting random forests, and XGBoost. The relative performance of the different methods is strongly regime-dependent, as summarized in Figure 1. For hard-voting random forests, M AX SAT is most effective, solving the largest fraction of instances the fastest, consistent with the target-class condition remaining close to a Boolean majority constraint. In soft-voting setups, however, it is no longer efficient, as it requires substantially larger Boolean encodings to model numerical confidence scores. In these settings, C PCF performs best, proving optimality for almost all instances while achieving the fastest certification times. When model-build time is excluded, visible as the x-axis offset before search begins, O CEAN and C PCF exhibit very similar solver-time performance. M ACE is markedly less competitive in all tested configurations (and does not support XGBoost ensembles). 7
RF soft voting
Median time (s)
103
XGB
CPCF OCEAN MACE MaxSAT
102 101 100 10 1
3
4
5
6
7
8
3
4
5 6 Max depth
7
8
3
4
5
(a) Scaling with the estimators’ depth. RF hard voting
RF soft voting
103 Median time (s)
RF hard voting
6
7
8
300
400
500
XGB
102 101
CPCF OCEAN MACE MaxSAT
100 10 1 0
100
200
300
400
500 0
100
200 300 400 Number of trees
500 0
100
200
(b) Scaling with the number of estimators.
Figure 4: Median total time to optimality (error bars indicate the first and third quartiles) with 100 trees of varying depth, and with varying numbers of trees of maximum depth 5. Both panels aggregate all ten datasets. Missing or non-optimal runs are counted at the 900-second time limit. Result 2. C PCF reaches high-quality incumbents earlier than O CEAN. Beyond certifying optimality, mathematical programming formulations such as CP and MILP can also be valuable because they provide anytime incumbent counterfactuals, together with lower bounds on the optimal cost, yielding a certifiable approximation gap for the current best solution. Figure 3 compares C PCF and O CEAN, the only benchmarked methods supporting anytime results, through the best incumbent objective found over time. Across all three regimes, C PCF improves the incumbent more rapidly in the early search phase. The difference is most pronounced for random forests, where the normalized counterfactual cost, with 0 denoting the optimal cost, decreases earlier and more sharply, indicating that C PCF reaches near-optimal counterfactuals sooner. The same pattern remains visible for XGBoost, although with a smaller gap than for soft-voting random forests. Thus, when CP and MILP are the relevant broadly applicable paradigms, C PCF is preferable not only because it certifies optimality faster, but also because it provides useful explanations earlier during the search. Result 3. Solution methodologies differ in their scalability. We now study the impact of two characteristics of the ensembles that significantly impact the size of the search space: the number of trees and their depth. Figure 4 reports the median total solution time as these parameters vary. For soft-voting random forests and XGBoost ensembles, C PCF scales substantially better than O CEAN across both the number of trees and their depth, while both remain far more efficient than the Boolean baselines. For hard-voting random forests, M AX SAT retains the best performance across much of the range, while C PCF remains competitive and scales better with the number of trees. These experiments reinforce our previous observations: M AX SAT is preferable for hard-voting random forests, but C PCF is most suitable once the model departs from that purely Boolean regime through numerical confidence scores. Result 4. C PCF remains efficient across different objective choices. Figure 5 compares the total solution times of C PCF and O CEAN under L0 , L1 , and L2 on the default ensemble configuration. We observe that C PCF remains faster than O CEAN across all datasets and norms. This gap is especially pronounced under L2 , for datasets with many numerical features. This is coherent with the formulations: in C PCF, changing the norm mostly amounts to changing the coefficients in the objective, in particular through (9). By contrast, the MILP formulation is more sensitive to the objective: under the L2 norm, it becomes a quadratic program. Result 5. C PCF efficiently accommodates plausibility constraints. Table 2 compares C PCF and O CEAN, the two approaches supporting plausibility constraints, with and without such constraints under the default soft-voting random forest configuration. Plausibility is enforced using an isolation 8
Runtime (s)
103
CPCF OCEAN
L0
L1
L2
102 101 100 AD
BC
CP
CC
GC
ON
PH
SE
SP
ST
Figure 5: Distribution of total times to optimality for C PCF and O CEAN under L0 , L1 , and L2 -type cost functions, for the default configuration (100 depth-5 trees) of soft-voting random forests. Table 2: Effect of isolation forest plausibility constraints on O CEAN and C PCF for the default configuration of soft-voting random forests. For each dataset, P reports the proportion of plausible explanations among all queries, and T reports the median total time to optimality in seconds. Method
Metric
AD
BC
CP
CC
GC
ON
PH
SE
SP
ST
C PCF-noIF
P (%) T (s)
96 4.62
22 1.87
82 2.02
70 3.06
68 3.36
82 6.28
80 2.07
18 1.62
50 11.58
90 1.43
O CEAN-noIF
P (%) 96 T (s) 10.31
26 84 70 72 7.16 11.96 15.74 13.92
82 41.88
78 18 6.84 20.19
50 61.92
88 5.82
C PCF-IF
P (%) 100 100 100 100 100 T (s) 10.92 16.51 17.78 12.85 12.23
100 21.60
100 100 9.29 18.50
100 43.44
100 9.51
O CEAN-IF
P (%) 100 100 100 100 100 100 100 100 100 100 T (s) 27.00 38.02 48.41 39.86 50.73 132.48 28.03 78.45 106.68 25.81
forest of 100 trees with a contamination level (proportion of training examples flagged as outliers) of 10%. For each dataset, we report the proportion of plausible explanations among all queries and the median total runtime. Two conclusions emerge. First, the constraint behaves as intended: without it, the proportion of plausible explanations varies substantially across datasets, whereas both formulations return fully plausible counterfactuals once it is enforced. Second, this plausibility layer increases the solution time for both methods, as expected, since it makes the optimization problem more computationally demanding. However, it does not change their relative ranking. Across all the datasets, C PCF remains consistently faster than O CEAN, both with and without the isolation forest. The gap is especially marked on larger or more structured datasets, where the additional plausibility constraints amplify the benefit of stronger global propagation in CP.
5
Conclusion
We introduced C PCF, a constraint programming formulation for computing optimal counterfactual explanations in tree ensembles. Across a broad empirical evaluation, C PCF achieves state-of-the-art total time to proven optimality, strong anytime performance, and robust scalability across ensemble size, tree depth, cost functions, and plausibility constraints. Our results also show that the relative strengths of different mathematical programming paradigms depend on the specific setting. Future work could extend C PCF in several directions. First, more complex cost functions could be considered, for instance when per-feature costs are interdependent or when feature modifications must comply with a given causal model. Efficient automata representations in CP solvers could, for instance, be leveraged to encode causal dependencies between admissible feature modifications. Second, integrating additional desiderata, such as sparsity or robustness, would naturally lead to multi-objective counterfactual formulations. Providing precise and certifiable characterizations of the trade-offs between these desiderata would be an interesting direction. Finally, adapting C PCF beyond tree ensembles (that are axis-aligned piecewise-constant models) is another promising direction, although it would require moving beyond the current per-feature interval representation. More broadly, our results suggest that the efficiency of each mathematical programming paradigm depends on how well its modeling primitives match the structure of the problem at hand. This 9
perspective extends beyond counterfactual explanations for tree ensembles to other important problems in trustworthy machine learning, such as verification, robustness analysis, fairness-constrained learning, and sparse explanations, which combine logical, continuous, and combinatorial components. Identifying which paradigm best captures these structures, or how to combine paradigms when they coexist, is an important direction for future work.
References Armin Biere, Tobias Faller, Katalin Fazekas, Mathias Fleury, Nils Froleyks, and Florian Pollitt. CaDiCaL 2.0. In Arie Gurfinkel and Vijay Ganesh, editors, Computer Aided Verification, volume 14681, pages 133–152, 2024. Leo Breiman. Random forests. Machine Learning, 45(1):5–32, 2001. Miguel A. Carreira-Perpinan and Suryabhan Singh Hada. Very fast, approximate counterfactual explanations for decision forests. Proceedings of the Thirty-Seventh AAAI Conference on Artificial Intelligence, 37(6):6935–6943, 2023. Magorzata Charytanowicz, Jerzy Niewczas, Piotr Kulczycki, Piotr Kowalski, and Szymon Lukasik. Seeds. UCI Machine Learning Repository, 2010. DOI: https://doi.org/10.24432/C5H30K. 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, pages 785–794, 2016. Zhicheng Cui, Wenlin Chen, Yujie He, and Yixin Chen. Optimal action extraction for random forests and boosted trees. In Proceedings of the 21th ACM SIGKDD international conference on knowledge discovery and data mining, pages 179–188, 2015. Leonardo De Moura and Nikolaj Bjørner. Z3: an efficient SMT solver. In 14th International Conference on Tools and Algorithms for the Construction and Analysis of Systems, page 337–340, 2008. Niklas Een and Niklas Sörensson. Translating pseudo-boolean constraints into SAT. JSAT, 2:1–26, 2006. Rubén R Fernández, Isaac Martín De Diego, Víctor Aceña, Alberto Fernández-Isabel, and Javier M Moguerza. Random forest explainability using counterfactual sets. Information Fusion, 63: 196–207, 2020. Alexandre M. Florio, Pedro Martins, Maximilian Schiffer, Thiago Serra, and Thibaut Vidal. Optimal decision diagrams for classification. Proceedings of the Thirty-Seventh AAAI Conference on Artificial Intelligence, 37(6):7577–7585, 2023. Shalmali Joshi, Oluwasanmi Koyejo, Warut Vijitbenjaronk, Been Kim, and Joydeep Ghosh. Towards realistic individual recourse and actionable explanations in black-box decision making systems, 2019. URL https://arxiv.org/abs/1907.09615. Amir-Hossein Karimi, Gilles Barthe, Borja Balle, and Isabel Valera. Model-agnostic counterfactual explanations for consequential decisions. In International conference on artificial intelligence and statistics, pages 895–905, 2020. Amir-Hossein Karimi, Gilles Barthe, Bernhard Schölkopf, and Isabel Valera. A survey of algorithmic recourse: contrastive explanations and consequential recommendations. ACM Computing Surveys, 55(5):1–29, 2022. Awa Khouna, Julien Ferry, and Thibaut Vidal. Counterfactual maps: What they are and how to find them, 2026. URL https://arxiv.org/abs/2602.09128. Fei Tony Liu, Kai Ming Ting, and Zhi-Hua Zhou. Isolation forest. In 2008 Eighth IEEE international conference on data mining, pages 413–422, 2008. Arnaud Van Looveren and Janis Klaise. Interpretable counterfactual explanations guided by prototypes, 2020. URL https://arxiv.org/abs/1907.02584. 10
Ana Lucic, Harrie Oosterhuis, Hinda Haned, and Maarten De Rijke. Focus: Flexible optimizable counterfactual explanations for tree ensembles. Proceedings of the Thirty-Sixth AAAI conference on artificial intelligence, 36(5):5313–5322, 2022. Daniel Nemirovsky, Nicolas Thiebaut, Ye Xu, and Abhishek Gupta. Countergan: Generating counterfactuals for real-time recourse and interpretability using residual gans. In Uncertainty in Artificial Intelligence, pages 1488–1497, 2022. Axel Parmentier and Thibaut Vidal. Optimal counterfactual explanations in tree ensembles. In Proceedings of the 38th International Conference on Machine Learning, pages 8422–8431, 2021. Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, Jake Vanderplas, Alexandre Passos, David Cournapeau, Matthieu Brucher, Matthieu Perrot, and Édouard Duchesnay. Scikit-learn: Machine learning in python. Journal of Machine Learning Research, 12(85):2825– 2830, 2011. URL http://jmlr.org/papers/v12/pedregosa11a.html. Laurent Perron and Frédéric Didier. optimization/cp/cp_solver/.
Cp-sat.
URL https://developers.google.com/
Tobias Philipp and Peter Steinke. PBLib – A Library for Encoding Pseudo-Boolean Constraints into CNF. In Theory and Applications of Satisfiability Testing – SAT 2015, volume 9340, pages 9–16. 2015. ISBN 978-3-319-24317-7. Rafael Poyiadzi, Kacper Sokol, Raul Santos-Rodriguez, Tijl De Bie, and Peter Flach. Face: Feasible and actionable counterfactual explanations. In Proceedings of the AAAI/ACM Conference on AI, Ethics, and Society, page 344–350, 2020. Alesya Raevskaya and Tuomo Lehtonen. Optimal Counterfactual Explanations for Random Forests with MaxSAT. Proceedings of the 28th European Conference on Artificial Intelligence, pages 903–910, 2025. Francesca Rossi, Peter Van Beek, and Toby Walsh. Handbook of constraint programming. Elsevier, 2006. Gabriele Tolomei, Fabrizio Silvestri, Andrew Haines, and Mounia Lalmas. Interpretable predictions of tree-based ensembles via actionable feature tweaking. In Proceedings of the 23rd ACM SIGKDD international conference on knowledge discovery and data mining, page 465–474, 2017. Sahil Verma, Varich Boonsanong, Minh Hoang, Keegan Hines, John Dickerson, and Chirag Shah. Counterfactual explanations and algorithmic recourses for machine learning: A review. ACM Computing Surveys, 56(12):1–42, 2024. WIlliam Wolberg. Breast Cancer Wisconsin (Original). UCI Machine Learning Repository, 1990. DOI: https://doi.org/10.24432/C5HP4Z. Haifei Zhang and Jinfeng Zhong. Efficient and effective counterfactual explanations for random forests. Expert Systems with Applications, 293:128661, 2025. ISSN 0957-4174.
11
A
Details on the Construction of Figure 1
Figure 1 is obtained by fitting a small decision-diagram meta-classifier to our experimental benchmark results, using the optimal decision-diagram learning method introduced by Florio et al. [2023] with a 600 second time limit. Each training row corresponds to one configuration–metric pair and is described by simple dataset and ensemble descriptors: the numbers of numerical, ordinal, categorical, and binary features; the ensemble type and voting regime; the number of estimators; the maximum depth; the total, per-tree mean, and per-tree maximum numbers of tree nodes; the total, mean, and maximum numbers of split levels associated with each feature; the objective norm; whether isolation forest plausibility is enabled; and the performance metric of interest, either median total running time or median solver-only time, for which model-building time is ignored. The resulting training dataset contains 3,380 rows, corresponding to 2 × 1,690: each of the 1,690 unique experimental configurations contributes one row for median total time and one row for median solver-only time. Each row is labeled with the best-performing method among M ACE, O CEAN, M AX SAT, and C PCF for the corresponding configuration and metric. M ACE is the best-performing method on only one row in the full 3,380-row table. Since the purpose of the decision diagram is to summarize the practically relevant regimes rather than isolate a single outlier configuration, we remove that row and train the final diagram on the remaining 3,379 rows using only the three dominant exact paradigms: C PCF, O CEAN, and M AX SAT. The percentages displayed in the nodes report the empirical fraction of retained rows routed through each node. The resulting decision diagram (1) achieves an accuracy of 87.79% on the full dataset.
B
M AX SAT Formulation and Soft-Voting Extension
Raevskaya and Lehtonen [2025] formulate hard-voting random forest counterfactual search as a partial weighted M AX SAT instance. Using our notation, the formulation is naturally expressed through a set of hard clauses Φhard and a set of weighted soft clauses ΦL1 . Note that in a partial weighted M AX SAT instance, hard clauses must be satisfied, whereas soft clauses may be violated at a cost. Each soft clause is associated with a nonnegative weight, and the objective is to minimize the total weight of violated soft clauses. Here, the soft clauses encode the distance from the original point x̂: preserving the original value of a feature satisfies the corresponding soft clauses, while moving away from it violates clauses whose weights add up to the L1 cost of that move. Variables. For each tree t ∈ T and leaf ℓ ∈ Lt , let zt,ℓ indicate that the counterfactual reaches leaf ℓ. For each categorical feature f ∈ FC and category j ∈ {1, . . . , kf }, let νf,j be the corresponding onehot literal. Binary features FB are special cases of categorical features and can be represented with a single Boolean literal; in what follows, we treat them together with categorical features. For each numerical or ordinal feature f ∈ FN ∪ FO and threshold τf,m , let tf,m (such that m ∈ {1, . . . , kf }) denote the literal stating that the counterfactual lies on the left side of τf,m . Hard clauses.
Feature validity is encoded by Φdom = Φcat ∧ Φthr ,
where Φcat =
^
kf _
f ∈FC
j=1
! ^
νf,j ∧
(¬νf,j ∨ ¬νf,j ′ ) ,
1≤j<j ′ ≤kf
which means that for all categorical features, exactly one category is selected, and kf −1
Φthr =
^
^
(¬tf,m ∨ tf,m+1 ),
f ∈FN ∪FO m=1
which enforces consistency among threshold literals: crossing a threshold implies crossing all smaller thresholds, or equivalently, lying on the left side of a threshold implies lying on the left side of all larger thresholds. 12
For each t ∈ T , exactly one leaf must be active: Φleaf =
^_ zt,ℓ ∧
t∈T
ℓ∈Lt
^ ℓ,ℓ′ ∈Lt ℓ̸=ℓ′
(¬zt,ℓ ∨ ¬zt,ℓ′ ) .
Let L(t, ℓ) denote the set of literals induced by the root-to-leaf path of leaf ℓ, namely literals of the form tf,m or ¬tf,m for numerical and ordinal splits, and νf,j or ¬νf,j for categorical splits. Path consistency is then encoded by ^ ^ ^ Φpath = (¬zt,ℓ ∨ λ). t∈T ℓ∈Lt λ∈L(t,ℓ)
If y(ℓ) is the class label of leaf ℓ, the hard-voting target-class condition is encoded by ^ X X X X Φclass = Card zt,ℓ ≥ zt,ℓ + εc . y∈Y\{y ⋆ }
t∈T ℓ∈Lt : y(ℓ)=y ⋆
t∈T ℓ∈Lt : y(ℓ)=y
The full hard part is therefore Φhard = Φdom ∧ Φleaf ∧ Φpath ∧ Φclass . L1 soft clauses. To facilitate comparison with the other formulations, we retain only the L1 part of the original objective. For a binary or categorical feature f , let ρf be the literal preserving the current value of the query x̂. The corresponding soft clauses are ^ Φcat (ρf , wt = αf ). L1 = f ∈FB ∪FC
where αf > 0 denotes the weight assigned to feature f in the weighted L1 objective. For a numerical feature f , let x̂f ∈ (τf,j , τf,j+1 ], and define τ̄f,m =
τf,m − lb(f ) , ub(f ) − lb(f )
x̄f =
x̂f − lb(f ) . ub(f ) − lb(f )
The exact incremental L1 encoding is Φnum L1 (f ) =
j−1 ^
(¬tf,k , wt = αf (τ̄f,k+1 − τ̄f,k ))
k=1
∧(¬tf,j , wt = αf (x̄f − τ̄f,j )) ∧ (tf,j+1 , wt = αf (τ̄f,j+1 − x̄f )) ∧
kf ^
(tf,k , wt = αf (τ̄f,k − τ̄f,k−1 )).
k=j+2
This encoding exploits the monotonicity of the threshold literals. Since tf,m denotes the condition xf ≤ τf,m , the original value x̂f ∈ (τf,j , τf,j+1 ] satisfies ¬tf,k for all thresholds below x̂f , and tf,k for all thresholds above x̂f . The soft clauses therefore reward keeping the counterfactual on the same side of each threshold as x̂f . If the counterfactual crosses a threshold, the corresponding soft clause is violated and its weight contributes to the objective. The weights are chosen as normalized interval lengths, so the sum of violated weights is exactly the weighted normalized L1 displacement along feature f . For ordinal features, the same construction is used after replacing each threshold by the nearest admissible value on the appropriate side. The full soft-clause family is then ^ ΦL1 = Φcat Φnum L1 ∧ L1 (f ). f ∈FN ∪FO
13
Soft-voting extension. The original formulation is specific to hard-voting tree ensembles. In our benchmark, we extend it to support soft-voting tree ensembles. We keep Φdom , Φleaf , and Φpath unchanged, and replace Φclass by ! ^ X X (y⋆ ,y) soft Φclass = PB p̄t,ℓ zt,ℓ ≥ εc , y∈Y\{y ⋆ }
t∈T ℓ∈Lt
⋆
(y ,y)
where p̄t,ℓ is an integer scaling of pt,ℓ,y⋆ − pt,ℓ,y , and PB(·) denotes a CNF encoding of the corresponding pseudo-Boolean inequality [Philipp and Steinke, 2015, Een and Sörensson, 2006].
C
Integrating Isolation Forest Plausibility Constraints into C PCF
For each target class y ⋆ , we train a separate isolation forest on the training examples labeled y ⋆ . This model estimates the training-data distribution within the target class, yielding a plausibility criterion aligned with the target region rather than with the full data distribution. (y ⋆ )
Let TIF denote the isolation forest associated with target class y ⋆ . As in the explained tree ensemble, each isolation tree contributes exactly one active leaf, and each active leaf is assigned a corrected path length. The path of the counterfactual through the isolation forest is encoded in the same way as for the explained tree ensemble, using leaf-selection variables and path-consistency constraints (2)–(5). (y ⋆ )
For a leaf ℓ of an isolation tree t ∈ TIF , let deptht (ℓ) denote its depth, i.e., the length of the path from the root of t to ℓ, and let mt,ℓ be the number of training samples that reached ℓ during isolation forest construction. The corrected path length is then defined as ( deptht (ℓ), if mt,ℓ ≤ 1, Lt,ℓ = deptht (ℓ) + c(mt,ℓ ), if mt,ℓ > 1, where c(m) is the standard isolation forest correction term, 0, c(m) = 1, 2(m − 1) 2 ln(m − 1) + 2γ − , m and γ ≈ 0.57721 is the Euler–Mascheroni constant.
if m ≤ 1, if m = 2, if m > 2,
Let zt,ℓ ∈ {0, 1} denote the leaf-selection variable for leaf ℓ in isolation tree t. The average corrected path length of the counterfactual through the target-class isolation forest is then computed with: X X 1 H(x) = (y⋆ ) Lt,ℓ zt,ℓ . (10) |TIF | (y ⋆ ) ℓ∈LIF t∈TIF
t
The corresponding anomaly score is score(x) = −2−H(x)/cmax , where mmax is the isolation forest subsample size and cmax = c(mmax ). The isolation forest decision function is then defined as decision(x) = score(x) − offset, where offset is a constant chosen from the fitted isolation forest to match the prescribed contamination level, namely the proportion of training examples classified as outliers. A counterfactual is considered plausible if and only if decision(x) ≥ 0. Because the exponential form cannot be represented directly in the CP model, we reformulate this condition as an equivalent linear lower bound on the average path length. The final plausibility constraint is therefore enforced as: H(x) ≥ Hmin , (11) 14
where Hmin = −cmax log2 (−offset). Intuitively, examples isolated early or falling into leaves with few training samples have shorter corrected path lengths. They therefore have lower values of H(x) and are more likely to be classified as outliers by the isolation forest. Both C PCF and O CEAN use this corrected path-length formulation. This slightly refines the original O CEAN isolation forest constraint, which was based on the selected leaf depths and a fixed averagedepth threshold. In Table 2, the “noIF” variants do not enforce Equation (11); their plausibility rates are therefore computed a posteriori using the same target-class-specific isolation model. In contrast, the “IF” variants enforce this constraint during optimization, and therefore produce 100% plausible counterfactuals by construction.
D
Implementation Details
D.1
Numerical Confidence Scores Handling
The mathematical programming formulations considered in this work all solve the same optimization problem. However, while O CEAN can model confidence scores using continuous variables, as discussed in Appendix D.2, the other paradigms require a discrete numerical representation. In C PCF, confidence scores used to enforce the target-class prediction are represented as integers. We therefore multiply these scores by 109 and round them to the nearest integer. In soft-voting M AX SAT, pairwise score differences are similarly scaled before the pseudo-Boolean encoding, so that the targetclass constraints can be represented in CNF form. For M ACE, we retain the numerical treatment of the authors’ implementation: the random forest encoding uses real-valued auxiliary variables for tree-level class probabilities and enforces the forest prediction by comparing the corresponding sums. These transformations do not change the intended target-class comparison in Equation (7); they only adapt its numerical representation to the requirements of the underlying solver. In C PCF and soft-voting M AX SAT, confidence scores are converted to integer quantities by scaling them by 109 and rounding. This introduces at most 5 × 10−10 absolute error per scaled score, which is below the precision at which the confidence scores are exposed by the tree-ensemble implementations and below the numerical thresholds used in our experiments. Thus, the discrete encodings preserve the same effective target-class comparisons as the original floating-point scores. For consistency, we also tighten the numerical tolerances of the MILP solver used by O CEAN. This avoids comparing O CEAN under a looser default feasibility or integrality tolerance than the discrete encodings used by C PCF and M AX SAT. Consequently, all formulations are evaluated under the same effective precision regime, and observed differences in performance are not attributable to different numerical tolerances across solver backends. D.2
O CEAN-Specific Implementation Details
To align the numerical precision of O CEAN with that of the other baselines, we make a few adaptations to the implementation. In particular, we set the Gurobi parameters FeasibilityTol=1e-9, IntFeasTol=1e-9, and IntegralityFocus=1, and reduce the prediction margin from 10−4 to 10−7 , matching the value used in C PCF. These changes do not modify the formulation itself, but make the returned counterfactuals more faithful to the reference tree-ensemble predictions. A second difference concerns the target class. While Parmentier and Vidal [2021] report experiments in which counterfactuals are generated only toward target class 1, we generate counterfactual explanations for any admissible target class y ⋆ ∈ Y \ {hT (x̂)}. This choice is consistent with the multiclass setting considered in our benchmark and with the unified problem definition used throughout the paper.
E
Additional Experimental Details
This appendix provides additional details on the experimental protocol used in Section 4. 15
E.1
Model Training and Query Generation
For each configuration, dataset, and seed, we train the ensemble on the full processed dataset and then sample 50 query instances without replacement from the same dataset. For binary tasks, the target class is set to the opposite label. For multiclass tasks, the target class is selected uniformly at random from Y \ {hT (x̂)}. Soft-voting random forests use the class scores returned by the fitted scikit-learn leaves. Hardvoting random forests are obtained from the same trained forests by replacing each leaf score vector with a one-hot vote for its majority class, while keeping the tree structure unchanged. XGBoost is evaluated only in the soft-voting regime. For the eight datasets inherited from Parmentier and Vidal [2021], we keep the same preprocessing pipeline, but do not enforce the dataset-specific actionability restrictions used in that paper. This choice avoids introducing dataset-dependent feasibility constraints that could interact differently with the competing formulations. Instead, all methods are evaluated on the same unrestricted processed feature space, corresponding to the largest common search space induced by the preprocessing. E.2
Evaluation and Plotting Conventions
All reported runtimes are total times, defined as model-construction time plus solving time, unless explicitly stated otherwise. Indeed, before an optimization formulation can be solved, the corresponding model must first be instantiated through the solver API: variables must be created, domains specified, and constraints posted. This model-construction phase, which we refer to as build time, can be non-negligible and may vary substantially across optimization paradigms and solver interfaces. Reporting total time, therefore, captures the end-to-end computational cost incurred when solving a query. In the scaling plots, runs that do not prove optimality within the time limit are all censored at 900 seconds before computing medians. The anytime curves compare only C PCF and O CEAN, since these are the only two approaches in our benchmark that expose intermediate incumbents. They are restricted to instances solved to optimality by both methods, which covers a wide majority of the queries (more than 99% of the queries in the worst case in our experiments). For each such instance, the normalized error at time t is defined as costt − cost⋆ , costmax − cost⋆ where costt is the best incumbent objective available at time t, cost⋆ is the final optimal objective, and costmax is the worst objective value observed across the two incumbent traces for that instance.
16