Many decision trees fit the same data almost equally well, yet they can route a query point to different leaves and induce different local empirical distributions. We study decisions that meet prescribed cost, shortage or risk targets despite this predictive multiplicity. We propose the joint Rashomon and robustness optimization framework for optimal decision trees. It jointly selects an operational decision and the largest Rashomon set of trees, so that the targets hold under the local empirical distribution that every tree in this set induces at the query point. We specialize the framework to the regression setting, and we show that a tree affects the decision only through the training observations sharing the query leaf, which we call its query neighborhood. As a result, the robust problem involves only finitely many distinct constraints, which can be examined in order of increasing estimation loss. We develop a constraint generation algorithm that combines query-path pricing with dynamic programming to identify violating neighborhoods without enumerating trees. On synthetic newsvendor instances, the algorithm typically needs few neighborhoods and runs substantially faster than full neighborhood enumeration. On restaurant demand data, the robust orders increase the mean tolerated excess estimation loss by 17.6% and reduce the empirical conditional value-at-risk of the worst 10% of realized costs by 8.3% relative to the sample average approximation orders of the optimal tree, while the mean cost difference is not statistically significant. An interpretability analysis further shows how the retained neighborhoods explain the decision and its robustness limit.
Figures & tables
Figure 1: The monotone scan of Proposition 1 , drawn for a scalar decision. The neighborhood S5 empties the feasible set at η5 , which is the supremum radius r⋆ , and the decision z⋆ is selected from Fp⋆ by tie-breaking.
Figure 2: Query path, query region, and sibling regions for D=2 , dx=2 and Nmin=3 . The path π=((a1,b1),(a2,b2)) retains at each level the half-space Hj(π) containing x0 , here H1(π)={a1⊤x>b1} and H2(π)={a2⊤x≤b2} . Their intersection is the query region R(π) , whose observations form the query neighborhood ST(x0)=N(R(π)) , while the complementary half-spaces yield the sibling regions Q1(π) and Q2(π) .
dx=1
dx=2
dx=4
D
N
Solved
Median (s)
90th (s)
Solved
Median (s)
90th (s)
Solved
Median (s)
90th (s)
2
50
100
< 0.01
0.02
100
< 0.01
0.03
100
0.02
0.07
100
100
0.08
0.18
100
0.02
0.23
100
0.05
0.51
150
100
0.19
0.49
100
0.04
0.57
100
0.08
1.21
200
100
0.38
1.33
100
0.06
1.26
100
0.12
2.62
3
50
100
0.04
0.12
100
0.10
0.20
100
0.42
0.78
Table 1: Computational performance of Algorithm 1 by maximum depth D , sample size N and number of features dx . Each cell aggregates 3600 runs. Solved is the percentage of runs that returned a proven outcome up to numerical tolerance (certified radius, attained cap, or infeasibility at r=0 ). All remaining runs stopped at the limit of 106 cache entries. We report the median and 90th percentile wall-clock times in seconds over solved runs.
Runs (%)
Outcome (%)
Median
z⋆=z0
Separations
States
Time share (%)
D
γ
Solved
Stop in ρ⋆
Inf.
Brk.
Cap
r⋆/ρ⋆
(%)
Med.
Max
( ×103 )
DP
Path
Master
2
0.10
100
–
20.4
79.6
0
0.17
22.3
2
9
1.6
67.1
29.0
2.0
0.25
100
–
10.7
89.3
0
0.23
32.0
2
13
2.1
59.8
37.3
1.4
0.50
100
–
9.2
90.7
< 0.1
0.36
35.4
2
21
2.8
55.4
43.3
0.7
0.75
100
–
8.8
89.2
2.0
0.50
42.2
2
33
3.2
52.9
45.6
0.7
1.00
100
–
8.1
85.1
6.8
0.65
40.7
2
49
3.7
53.8
44.8
0.6
Table 2: Computational profile of Algorithm 1 by maximum depth D and target slack γ , with 8640 runs per row. Solved is the percentage of runs returning a proven outcome up to numerical tolerance. Stop in ρ⋆ is the percentage of unsolved runs that reached the cache limit while computing the nominal optimum ρ⋆ . Outcome splits solved runs into infeasibility at r=0 (Inf.), an excluded breakdown threshold (Brk.), and an attained cap rˉ (Cap). The median normalized value r⋆/ρ⋆ and the percentage of decisions differing from the nominal decision z0 are computed over solved Brk. and Cap outcomes only. Separations reports the median and maximum number of calls to the separation problem, and States reports the median number of memoized subtree values Ψh(A) , in thousands. Both are computed over all solved runs, including infeasible outcomes. Time share splits the solved runtime into evaluation of the dynamic program (DP), path search excluding DP (Path), and the master problem (Master). The remainder is fixed overhead such as model setup and result assembly.
N
Method
Total (s)
Generated neighborhoods
Distinct DP subtree states
50
Constraint generation (Alg. 1 )
0.0130
2.88
925
Constraint generation without pruning
0.0841
2.88
7819
Full neighborhood enumeration
0.0963
6215
6446
100
Constraint generation (Alg. 1 )
0.0359
2.96
2962
Constraint generation without pruning
0.5674
2.96
42482
Full neighborhood enumeration
0.7784
31573
30939
Table 3: Computational comparison with complete-tree enumeration, full query-neighborhood enumeration, and the pruning ablation. Entries are means over five training seeds and ten queries per seed, with depth D=2 , dx=4 , Nmin=5 , and λ=0.1 , on instances generated with σ=2 and solved with γ=1 . The upper block reports cached dynamic programming (DP) states, while the tree enumeration reports distinct recursive enumeration subproblems. All reported runs were solved, and all paired comparisons with Algorithm 1 agreed with absolute and relative comparison tolerances of 10−7 .
Realized cost
Method
Mean order
Mean (SD)
Mean difference: 95% CI
CVaR0.1
CVaR0.1 difference: 95% CI
Stockout days (%)
Operational success (%)
Methods returning a certificate
ODT-JRRO
36.9
13.8 (17.4)
–
51.3
–
21.6
81.0
JERO-OLS
28.5
18.2 (23.9)
[2.2,7.4]
73.2
[13.7,29.8]
63.4
67.3
Methods without a certificate
ODT-SAA
34.9
13.9 (19.4)
[−0.5,0.9]
55.9
[1.8,8.0]
28.8
81.0
Table 4: Held-out comparison on chicken demand from the YAZ dataset over the 153 consecutive days of the held-out period, with D=3 , Nmin=20 , λ=0.01 and γ=0 ( τ≈18.07 ). Differences are computed as method minus ODT-JRRO, so positive values favor ODT-JRRO. The pointwise 95% intervals are percentile intervals from a paired moving-block bootstrap with blocks of 7 consecutive days and 50000 resamples. Operational success is the percentage of days on which a method issues an order with realized cost at most τ . All methods issue orders on every test date.
Figure 3: Expansion of the admissible neighborhood family (upper panel) and contraction of the feasible-order interval (lower panel). The full family contains 26939 neighborhoods strictly below the breakdown threshold and 26940 at r⋆ , marked by the vertical dashed line. The shaded interval contains orders satisfying all admitted constraints, and the horizontal dashed line marks q⋆ . The feasible set is empty at r⋆ .
Figure 4: Membership of chronologically ordered training observations in the five retained neighborhoods. The right-hand columns report neighborhood size, entry level, and feasible-order interval.
Figure 5: Empirical mean cost as a function of the order for the nominal and breakdown neighborhoods. The horizontal dashed line is the target and the vertical dotted line is q⋆ . The minimum cost on the breakdown neighborhood exceeds the target, so no order can preserve feasibility when this neighborhood enters.
Figure 6: Nominal tree (left) and a minimum-loss tree inducing the breakdown neighborhood S5 (right). The highlighted paths route the query to its terminal leaf. Each leaf reports its sample size, mean demand, and empirical newsvendor order.
Figure 7: Sensitivity to the minimum leaf size. Left: finite breakdown thresholds (circles) and robustness over the entire admissible tree class (triangles). Full-class cases have infinite uncapped breakdown thresholds and are displayed at the common ceiling Rcommon=Vnorm+λ(2D−1)≈0.2335 , where Vnorm is the empirical variance of normalized training demand. Right: empirical upper-tail CVaR of the newsvendor cost at the returned order q⋆ , its minimum over q∈Z , and the target τ , using tail fraction Nmin/N . The vertical dotted lines mark Nmin/N=1/2 , above which the only feasible tree consists of a single leaf at the root.
The Rashomon effect describes the phenomenon that many models can achieve nearly equivalent performance on the same learning task, with significant ramifications for robustness, feature importance, and customizability. These use cases motivate the computation of Rashomon sets: the set of all models whose regularized loss is near-optimal. Decision trees are one of the few model classes for which Rashomon sets can be fully enumerated, but this computation has always been conditional on a binarization of the original data, either restricting which splits each tree is allowed to make or substantially increasing the complexity of an already difficult combinatorial problem. We introduce the first algorithm that exactly enumerates decision-tree Rashomon sets while exploiting the ordered structure of continuous features. We further develop a relaxation for approximate enumeration and an anytime algorithm that progressively refines the set of candidate thresholds, producing increasingly detailed approximations that converge to the continuous-feature Rashomon set. Experiments show that coarse binarization can miss many trees, important features, and predictive multiplicity; our algorithms achieve orders-of-magnitude speedups over existing enumeration methods, with approximations providing further speedups while maintaining near-perfect recall.
Zakk Heile, Hayden McTavish, Margo Seltzer +1
Department of Computer Science, Duke University, Durham, USA · Department of Computer Science, University of British Columbia, Vancouver, Canada
Standard machine learning pipelines often admit many near-optimal models. These "Rashomon sets" pose a range of challenges and opportunities for uncertainty-aware, robust decision making. They allow users to incorporate domain knowledge and preferences that would otherwise be difficult to specify directly in an objective, and they quantify diversity among valid models for a given training dataset and objective function. However, computation of Rashomon sets, even for simple, interpretable model classes such as sparse decision trees, continues to require immense memory and runtime resources. We present PRAXIS, an algorithm to approximate this Rashomon set with orders of magnitude improvement in runtime and memory usage. We validate that PRAXIS regularly recovers almost all of the full Rashomon set. PRAXIS allows researchers and practitioners to scalably model the Rashomon set for real-world datasets. Code for PRAXIS is available at https://github.com/zakk-h/PRAXIS
Zakk Heile, Hayden McTavish, Varun Babbar +2
Department of Computer Science, Duke University, Durham, USA · Department of Computer Science, University of British Columbia, Vancouver, Canada.
Robust optimization (RO) provides a principled framework for decision-making under uncertainty, but its performance critically depends on the choice of the uncertainty set. While large sets ensure reliability, they often lead to overly conservative decisions, whereas small sets risk excluding the true outcome. Recent data-driven approaches, particularly conformal prediction, offer finite-sample validity guarantees but remain largely task-agnostic, ignoring the downstream decision structure. In this paper, we propose a decision-aware conformal framework that learns uncertainty sets tailored to robust optimization objectives. Our approach parameterizes a flexible family of polyhedral sets via data-driven hyperplanes and learns their geometry by directly minimizing the induced robust loss, while preserving statistical validity through conformal calibration. To correct for data-dependent selection, we incorporate a re-calibration step on an independent dataset to restore coverage. The resulting sets capture directional and anisotropic uncertainty aligned with the decision objective while remaining computationally tractable. We provide finite-sample coverage guarantees and bounds on the sub-optimality gap to an oracle decision. This work bridges the gap between statistical validity and decision optimality, providing a principled framework for data-driven robust optimization.