arXiv CS open access

Learning Early-to-Final Solution Consistency for MILP Acceleration

2026 · arxiv_cs

arXiv CS · Papers · 2026

Conceptio Open Archive

Analyze & explore related research in Conceptio

Connect this record directly into your research workflow or AI agent. Explore citations across open repositories, query full text, and export structured proof bundles.

Explore in Search →Connect via MCP / API ↗

Learning Early-to-Final Solution Consistency for MILP Acceleration Guanlin Li1, 2∗ , Chengrui Gao1, 2∗ , Chenguang Wang1, 2 , Haopu Shang1, 2 , Zherong Zhang1, 2 , Ke Xue1, 2 , Jixiang Lu3 , Weiyong Yang3 , Chao Qian1, 2† 1

arXiv:2608.19953v1 [cs.AI] 20 Aug 2026

State Key Laboratory of Novel SoftwareTechnology, Nanjing University 2 School of Artificial Intelligence,Nanjing University 3 State Key Laboratory of Technology and Equipment for Defense Against Power System Operational Risks, Nari Technology Co., Ltd. Abstract Mixed-Integer Linear Programming (MILP) is a fundamental problem class in operations research and combinatorial optimization, with broad applications to industrial decisionmaking. Owing to their NP-hardness, however, modern solvers may struggle to find high-quality solutions for challenging MILP instances within practical time limits. Recent learningbased approaches seek to accelerate MILP solving by directly predicting high-quality solutions from static instancelevel features, such as variable-constraint bipartite graphs. Yet accurate solution prediction from instance features alone is difficult, and these methods largely overlook the information revealed during the solver’s search process. In this paper, we find that solutions produced at the early search stage of MILP solvers, which are computationally cheap to obtain, are often structurally close to the solutions found after full-budget search. Motivated by this observation, we propose a new solver-informed paradigm that shifts the learning target from variable assignment to early-to-final consistency: for each variable, we predict whether its early-stage assignment should persist in full-budget solutions. The predicted consistency naturally guides downstream search, for instance by fixing the assignments deemed consistent. At inference time, we further ensemble consistency predictions across multiple early-stage solutions to improve robustness. Experiments across four MILP benchmarks show our method improves prediction-guided search across diverse downstream pipelines. With Gurobi, our proposed method reduces the primal gap by 56.9% on average and closes it completely on combinatorial auction instances. Besides, we transferred the Gurobi-trained model zero-shot to SCIP without adaptation, achieving a 36.4% average gap reduction across benchmarks. Code is available at: https://github.com/lamda-bbo/EnCore.

Introduction

Mixed-Integer Linear Programming (MILP) is a fundamental modeling language for combinatorial optimization and decision making, with broad applications in scheduling (Floudas and Lin 2005), routing (Laporte 1992), network design (Crainic 2000) and production planning (Nam and Logendran 1992). Prevailing MILP solvers employ branch-andbound (Land and Doig 1960) frameworks to systematically ∗ †

These authors contributed equally. Corresponding author. Email: qianc@lamda.nju.edu.cn.

explore the combinatorial search space while guaranteeing optimality. Built upon this, the MILP community has developed a variety of sophisticated techniques (Padberg and Rinaldi 1991; Achterberg 2009) over the past decades, including presolve, cutting planes, primal heuristics, and branching strategies, which further reduce the search space and accelerate the solving process. However, solving MILP remains computationally expensive for modern solvers due to their inherent NP-hardness (Karp 1972). In real-world applications, most MILP instances are derived from similar problem classes with varying parameters, exhibiting shared structural patterns that can be naturally exploited by machine learning methods. As investigated by Bengio, Lodi, and Prouvost (2021), learning models can capture these patterns from instance characteristics and historical solving experiences, and infer effective optimization strategies. Such learned knowledge can be integrated into MILP solvers to accelerate the optimization process (Gasse et al. 2019; Li et al. 2023; Nair et al. 2020). Existing learningbased methods for MILP can be broadly classified into three types: (1) learning solver decisions, where machine learning models are employed to guide individual solver component, such as branching, node selection, and cutting plane selection (Khalil et al. 2016; Gasse et al. 2019; Gupta et al. 2020; Ferber et al. 2020; Tang, Agrawal, and Faenza 2020; Ling, Wang, and Wang 2024; Strang et al. 2026); (2) learning solver configurations, where models automatically tune algorithmic parameters or select suitable configurations for different instances (Hutter et al. 2009; Li et al. 2023); and (3) learning solution predictions, where models directly predict complete or partial solutions to provide high-quality initial solutions or reduce the search space (Nair et al. 2020). Among the solution prediction methods, recent advances have focused on two directions: improving learning-guided search mechanisms and enhancing prediction capacity. On the search side, early studies such as neural diving (Nair et al. 2020) fixed high-confidence variable assignments predicted by neural models, producing reduced MILP instances for efficient downstream optimization. Building upon this idea, Han et al. (2023) replaced hard variable fixing with a trust-region search strategy, reducing the risk of being trapped by inaccurate assignment. More recently, Liu et al. (2025) proposed to alternate prediction with trust-region correction, further improving robustness against prediction errors. On the predic-

tion side, Wang, Li, and Wang (2026) represented MILPs as geometric feature sequences and used multi-scale attention to overcome expressiveness limitations of graph-based models. Meanwhile, Li et al. (2026) introduced a flow-matching generator that jointly models integer and continuous variables to enhance solution generation quality, while Pu et al. (2026) exploited variable interactions through inter-variable contrastive learning and an intra-constraint competitive graph network. However, these methods still follow the conventional paradigm of predicting a complete solution directly from the MILP formulation. Under this paradigm, the model is required to determine all variable assignments, leading to a prediction task that approaches the complexity of solving the original MILP. In this paper, we begin by analyzing how primal solutions evolve during the search of modern MILP solvers. Across all evaluated benchmarks, we observe that the solution quality improves rapidly in the early stage of the search and then progresses much more slowly with only local corrections, indicating that early incumbents are highly informative about the final, full-budget assignments. On workload apportionment, for instance, an early incumbent already agrees with the final solution on 95.63% of the binary variables (as shown in Figure 2). Given such an informative starting point, the optimal solution need not be reconstructed from scratch. Conditioned on what the solver has already produced, the learning task reduces to deciding which of its assignments to trust. Based on this insight, we propose a new solver-informed paradigm that reformulates solution prediction as an early-tofinal consistency estimation problem. Given an MILP graph G, an early solution X ES , and a reference solution X ∗ obtained from a full-budget run, conventional solution prediction learns P (X ∗ | G) and thus requires the model to infer all variable assignments from static instance features alone. In contrast, we estimate P (x∗i = xES | G, X ES ), namely i the probability that each early assignment persists in the final solution. This target can substantially ease the learning task. Conditioned on the rich anchor X ES , the model only needs to identify a relatively small fraction of assignments that have yet to stabilize. Theoretically, we prove that conditioning on the early solution can strictly increases the best achievable accuracy of predicting the full-budget solution, under a mild posterior-crossing condition, and that this information gain provably persists under finite-sample model selection. Moreover, since the predicted consistent assignments are drawn from a feasible early solution, fixing them preserves feasibility and concentrates the subsequent search on assignments that the solver has left unresolved. We evaluate our method on multiple MILP benchmarks with Gurobi (Gurobi Optimization, LLC 2026) as the target solver and integrate it into existing prediction-based search frameworks. Experimental results show that our proposed refiner consistently reduces the primal gap in most evaluated settings. Specifically, it reduces the primal gap of Predict-and-Search (Han et al. 2023) by 56.9% on average and closes the gap on combinatorial auction instances. Furthermore, without retraining, the model trained with Gurobigenerated data can be directly transferred to SCIP (Achterberg 2009). These results demonstrate that the proposed

early-to-final consistency prediction paradigm is clearly effective and showcases transferability across different solvers.

Preliminaries

Mixed Integer Linear Programming A mixed-integer linear program (MILP) is an optimization problem with a linear objective and linear constraints, where a subset of decision variables is restricted to integer values. We formulate a MILP instance as minn c⊤ x x∈R

s.t.

Ax ≤ b, l ≤ x ≤ u, xi ∈ Z ∀i ∈ I,

(1)

where x = (x1 , . . . , xn ) is the vector of decision variables, c ∈ Rn is the objective coefficient vector, A ∈ Rm×n is the constraint matrix, and b ∈ Rm is the right-hand-side vector. The bounds l, u ∈ (R∪{±∞})n defines the variable domains, and I ⊆ {1, . . . , n} is the index set of variables required to take integer values. An important special case is that of binary variables, for which xi ∈ {0, 1}. General integer variables can be transformed into binary variables using standard preprocessing techniques.

Graph Representations for MILP Learning-based MILP methods commonly encode an instance as a bipartite graph (Salvagnin 2016; Gasse et al. 2019). For the MILP formulated in Eq. (1), the associated graph is defined as G = (V, C, E), where V = {v1 , . . . , vn } and C = {q1 , . . . , qm } denote the sets of variable nodes and constraint nodes, corresponding to decision variables and constraints, respectively. An edge (vi , qj ) ∈ E is introduced when the coefficient Aji is non-zero. Each node is initialized with features describing its role in the MILP (Han et al. 2023). A variable node vi contains the objective coefficient ci , the lower and upper bounds li and ui , and an indicator of whether xi is continuous or integer. A constraint node qj contains the right-hand-side value bj and an encoding of its constraint sense, such as ≤, =, or ≥. In addition, each edge carries the matrix coefficient Aji as its edge feature. Collectively, these node and edge attributes provide a description of the local algebraic structure of the MILP instance, which serves as the input representation for graph neural networks (GNN).

Solution Predictions for MILP Solution prediction methods learn a mapping from an MILP instance to variable assignments, which can be exploited to reduce the search space for acceleration. Most existing predictors take the static bipartite graph G as input and employ a graph convolutional network (GCN) or a related messagepassing GNN to compute variable embeddings, from which the final predictions are decoded. (ℓ) Let hi be the embedding of variable node vi at layer ℓ, (ℓ) and let gj be the embedding of constraint node qj . A common bipartite message-passing layer first aggregates messages from variables to constraints and then sends messages

where rθ is a learnable readout function and σ is the sigmoid function. The output p̂i is the predicted probability that variable i takes value one in the reference solution. The predictions can then be used to construct a partial assignment, define a neighborhood around a candidate solution, or fix variables before invoking the solver on a reduced MILP.

Method

Motivation We begin with an observation about the dynamics of modern MILP solvers: the primal and dual gaps typically decrease sharply at the beginning of the solve, after which progress becomes much slower and is spread over a long period of search. Figure 1 illustrates this behavior with the evolution of the primal gap on example instances; the same qualitative pattern holds across solvers and problem classes (App. A). This behavior suggests a two-stage view of MILP solving. In the fast descent stage, the solver rapidly discovers useful early solutions and tightens bounds through presolve, heuristics, cutting planes, and early branch-and-bound decisions. In the long exploration stage, the solver spends most of its time searching for further improvements, proving bounds, and resolving difficult variables. This two-stage view motivates a new prediction target for learning-based MILP solving. In contrast to the prevailing paradigm, which constructs a full solution directly from the static MILP graph, we propose an Early-to-final solution Consistency prediction paradigm, termed EnCore. Specifically, given an early solution X ES collected after the fastdescent stage, we estimate the probability that the early assignment of each variable persists in the full-budget solution X ∗ . This target distinguishes assignments that have likely stabilized from those that merit further exploration, thereby directing both learning and search effort toward the decisions that remain open after the fast descent. The proposed paradigm rests on the premise that the early solution is informative of the final one. We now examine this premise empirically. Figure 2 compares each early solution X ES with the full-budget solution X ∗ by counting the flipped integer variables. Across benchmark classes, the

Primal Gap to BKS (%)

80 60 40 20 0 0

50 100 150 200 400 600 800 1000

100 80 60 40 20 0

Time (s)

(a) Gurobi Results

50 100 150 200 400 600 800 1000

Time (s)

(b) SCIP Results

Figure 1: Evolution of the primal gap over the solving process with Gurobi (left) and SCIP (right) on set covering instances. The time axis is scaled to highlight the early phase 0−200 s. Results on other problem classes are provided in App. A. Mean = 19.87% Median = 20.83%

Number of Instances

implemented as multilayer perceptrons with nonlinear activations. Each update depends on both neighboring node embeddings and the corresponding edge weight Aji . After L message-passing layers, the predictor applies a variable-level readout to each final variable embedding. For binary integer variables, a common formulation trains a predictor fθ on the graph G to estimate a target solution X ∗ : (L) p̂i = σ rθ (hi ) = fθ (G)i ≈ P (x∗i = 1 | G),

100

6 4 2 0

Flipped Binary Variables (%)

Mean = 4.37% Median = 1.80%

25 20 15 10 5 0

(a) Gurobi on CA Mean = 21.76% Median = 21.90%

10 8 6 4 2 0

Flipped Binary Variables (%)

(b) Gurobi on WA Number of Instances

(ℓ) (ℓ) (ℓ) where ψC and ψV are message functions, and ϕC and (ℓ) ϕV are node-update functions. These functions are usually

Primal Gap to BKS (%)

(ℓ+1) hi

(ℓ) (ℓ) P (ℓ) (ℓ) = ϕC gj , i:(vi ,qj )∈E ψC (hi , Aji ) , (ℓ) (ℓ) P (ℓ) (ℓ+1) = ϕV hi , j:(vi ,qj )∈E ψV (gj , Aji ) ,

Number of Instances

(ℓ+1)

Number of Instances

back to variables:

Flipped Binary Variables (%)

Mean = 4.21% Median = 2.50%

Flipped Binary Variables (%)

(d) SCIP on WA

Figure 2: Distributions of flipped integer variables between early and full-budget solutions. ‘CA’ and ‘WA’ stand for combinatorial auction and workload apportionment, respectively. early solutions exhibit agreement with the final ones, and the discrepancies tend to concentrate on a limited subset of variables rather than spreading arbitrarily across the instance. These observations indicate that the early solution carries rich information about the final solution, making it a valuable reference worth exploiting.

Consistency Prediction Given the MILP graph G and the collected early solution X ES , we design a consistency predictor that estimates, for each integer variable, the probability that its early assignment agrees with the full-budget solution X ∗ . The predictor must fuse two heterogeneous sources of information: the structural description of the instance, encoded in the bipartite variable–constraint graph, and the variable assignments in the early solution. We achieve this fusion at the input level, where each integer-variable node is augmented with its early assignment xES as an additional feature. Through rounds i of variable-to-constraint and constraint-to-variable message

passing, the network can thus assess each assignment jointly with the constraints it appears in and the assignments of neighboring variables, which provides the evidence needed to judge whether the assignment will persist. The backbone is a bipartite graph neural network with half-convolutions, a standard architecture for MILP representation learning (Han et al. 2023; Liu et al. 2025). A variable-level readout then produces the consistency probability p̂i for each integer variable. Architectural details are provided in Appendix B. To train the predictor, we fundamentally reshape the learning target, rather than imitating the complete solutions X ∗ . Since the early solution already encodes substantial information about the final one, the remaining task is not to predict every variable value from scratch, but to assess which early assignments are likely to persist in the final solution. We therefore define the binary consistency label yi = 1[x∗i = xES i ] for each integer variable i, and train the predictor fθ so that p̂i = fθ (G, X ES )i ≈ P (x∗i = xES | G, X ES ). i For each training instance G, we collect an early solution X ES and pair it with the full-budget solution X ∗ . The collection procedure is described later in this section. The predictor is then trained on the early-to-final consistency labels with the binary cross-entropy loss, X L(θ) = − [yi log p̂i + (1 − yi ) log(1 − p̂i )] . i∈I

This consistency target reshapes the learning task: rather than constructing the full solution from the static graph alone, the model refines an informative reference solution by assessing the stability of each assignment, which is a more focused and better-conditioned problem. Moreover, the conditioning solution X ES provides a rich source of instance-specific information that purely graph-based predictors cannot access. Early Solution Collection. We collect early solutions at the transition from fast descent to long exploration. During the solving process, we monitor feasible-solution improvement events and estimate the local dual-gap decay rate over a sliding window, i.e., the gap reduction within the window divided by the elapsed time. The run is terminated once this rate first falls below a threshold, indicating that rapid bound improvement has mostly ended, and the last explored solution is taken as the early solution X ES . In practice, we additionally impose a minimum probing time to avoid stopping before a meaningful dual bound, and a maximum time to restrict the runtime cost. Details are included in App. C. Inference-Time Augmentation. The early solution available for conditioning is not unique. Depending on when the collection run terminates, the solver may return slightly different solutions, and a variable predicted stable under one of them may still be predicted fluctuating under another. To mitigate this sensitivity, we employ the last few improving solutions for inference-time augmentation. Let {X1 , . . . , XK } denote the last K improving solutions found during the collection run, with XK = X ES . The trained predictor processes each Xk as the conditioning solution. Since a consistency score is always expressed relative to the value taken in

its own conditioning solution, the resulting predictions must be aligned to the common reference X ES . To this end, we flip the sign of each logit whose corresponding variable takes a different value in Xk than in X ES , as a high consistency score for such a variable argues against the reference value. The aligned logits are then averaged, ! K 1 X p̄i = σ τk,i p̂k,i , K k=1

where τk,i = +1 if xk,i = xES and −1 otherwise. The i ensembled probability p̄i replaces p̂i as the consistency score of xES used by the downstream search methods. i

Combine EnCore with Search Methods Our consistency predictor is independent of the downstream search strategy. In most learning-based MILP accelerators, the method must select a subset of variables whose predicted values are trusted, and the search space is then reduced accordingly. The central question is therefore which assignments can be trusted. Existing methods trust an assignment when the model is confident in its own prediction, so both the value and its evidence come from the model alone. In contrast, our early-to-final consistency design trusts values that the solver has already realized in a feasible early solution, and scores them by predicting whether they persist after the full solving process. Specifically, we rank variables by their consistency score p̄i and select the top-ranked ones. Each selected variable i is constrained to its early value xES i , while the remaining variables are left to the downstream search mechanism. This consistency-based design offers two main advantages: (1) Feasibility. The fixed values are drawn from a feasible solution, so fixing them never destroys feasibility, which model-generated assignments cannot guarantee. (2) Alignment with the search process. The consistency signals originate from the solver’s own search trajectory, and thus naturally guide the model toward the variables that the solver has yet to resolve.

Theoretical Analysis

Our method depends on the premise that conditioning on the early solution X ES makes predicting the full-budget solution X ∗ fundamentally easier than predicting from the instance graph alone. In this section, we formalize it by answering two questions: (1) Can the early solution improve the best achievable prediction accuracy, and under what condition is the gain strict (Theorem 1)? (2) Does this gain survive when the predictor must be selected from finite training instances (Theorem 2)? Complete proofs are in Appendix F. Let D be the distribution over triples (G, X ES , X ∗ ), and let J be a uniformly weighted variable index within an instance, so that Pr(b xJ = x∗J ) is the expected per-instance variable accuracy, i.e., the population variable-level evaluation criterion used in our analysis. A message-passing predictor at variable J considers only the instance features in its receptive MILP graph RJ , while EnCore additionally atES taches early assignments XR to construct the augmented J ES input UJ = (RJ , XRJ ).

A solution predictor h maps RJ to a prediction of x∗J , matching the conventional paradigm. A consistency predictor a observes UJ and decides whether to retain the early value, inducing ES xJ , a(UJ ) = 1, ha (UJ ) = 1 − xES , a(UJ ) = 0. J Given xES J , consistency decisions and final assignments determine each other (see Appendix F), so the two paradigms can be compared as predictors of x∗J . Information gain of the early solution. Define the best population accuracies attainable from each input, A∗sol := sup Pr h(RJ ) = x∗J , h A∗con := sup Pr ha (UJ ) = x∗J . a

If D were known, these would be achieved by the Bayesoptimal rules on each input (Devroye, Györfi, and Lugosi 1996). Their comparison isolates the informational value of the early solution from any modeling concern. Theorem 1. For any fixed early-solution collection procedure, A∗con ≥ A∗sol . Let ηJ := Pr(x∗J = 1 | UJ ). The inequality is strict if ! Pr(ηJ < 12 | RJ ) > 0 Pr > 0. and Pr(ηJ > 12 | RJ ) > 0 Proof sketch. Every solution predictor h defines ah (UJ ) = 1{h(RJ ) = xES J }, for which hah (UJ ) = h(RJ ). Hence the augmented input can reproduce every solution predictor. For strictness, conditional Jensen’s inequality (Kallenberg 2021) shows that refining the input strictly increases Bayes accuracy under the stated condition. The inequality itself is unsurprising since additional input can never hurt a Bayes-optimal predictor. The substantive content is about the strictness condition, which we call posterior crossing: for a non-negligible set of instance inputs, different early assignments compatible with the same instance structure favor different final values. In this regime, no predictor on instance features alone, however expressive, can resolve the ambiguity, while the early solution can. Figure 2 provides complementary empirical motivation by showing that early-to-final discrepancies are often concentrated on a small fraction of variables. A sparse-correction reading of the gain. Let a∗ be the Bayes rule on the augmented input, p = Pr(xES ̸= x∗J ) J ∗ represents the early-solution error rate, ρ = Pr(a (UJ ) = 0) denotes its flipping rate, and q is the precision on flipped variables, where q > 1/2 due to the Bayes optimality. Then we have 1 − A∗con = p − ρ(2q − 1), ∗ Acon − A∗sol = 1 − A∗sol − p + ρ(2q − 1). The gain is positive exactly when p < 1 − A∗sol + ρ(2q − 1). This identity captures the residual nature of EnCore: corrections with q > 1/2 remove more errors than they introduce,

and when early-to-final discrepancies are sparse, learning can focus on a targeted subset of unstable assignments. The condition can also hold when 1−p < A∗sol (see Appendix F), showing that the early solution is valuable not only as a feasible candidate but also as an informative conditioning signal. Finite-sample guarantee. Theorem 1 optimizes over all rules and thus characterizes only the information available in principle. However, the predictor is selected from finite data in practice, and one may worry whether the population gain can survive estimation errors. Following the standard finite-class empirical-risk-minimization analysis based on Bernstein’s inequality and a union bound (Boucheron, Lugosi, and Massart 2013), we fix a finite class A = {a1 , . . . , aN } of consistency rules before observing m independent training instances, and let b a minimize the average proportion of variable errors; only instances need be independent, matching how MILP training data are collected. The class advantage decomposes as ∆A = ∆B − αA , where ∆B := A∗con − A∗sol is the gain from Theorem 1 and αA := A∗con − maxa∈A Acc(ha ) is the approximation gap, with Acc(h) = Pr(hJ = x∗J ). Motivated by the sparse discrepancies in Figure 2, we assume the expected per-instance flip rate satisfies maxa∈A E[sa ] ≤ s̄. Theorem 2. Let b hES = hba , and let b hsol be any possibly data-dependent solution predictor. For every δ ∈ (0, 1), with probability at least 1 − δ, Acc(b hES ) − Acc(b hsol ) ≥ ∆B − αA − εm (δ), r 2s̄ 2N 8 2N εm (δ) = 2 ln + ln . m δ 3m δ Proof sketch. Use the rule that always keeps the early solution as a common reference. A candidate differs from it only on variables it flips, so its instance-level loss difference has variance at most s̄. Bernstein’s inequality with a union bound over the N rules (Boucheron, Lugosi, and Massart 2013) then controls all empirical loss differences simultaneously. Empirical risk minimization converts this into the penalty εm (δ), while Bayes optimality lower-bounds the risk of every instance-input solution predictor. Thus the learned consistency predictor provably outperforms any solution predictor, even one trained with unlimited data, whenever p ∆B > αA + εm (δ). p The dominant term of εm (δ) scales as s̄/m rather than 1/m: sparse correction rules are cheaper to select reliably. A numerical illustration in Appendix F shows that under representative parameters, fewer than 200 training instances suffice, which is similar to the typical sizes of training datasets. In summary, Theorem 1 establishes that the early solution provides a strict population-level information gain under a stated posterior-crossing condition, and the sparse-correction identity explains where the gain comes from: repairing a small, better-than-chance set of unstable assignments. Theorem 2 shows that this gain persists under finite-sample model selection, with a sample cost discounted by the very sparsity that motivates our design. Together, they justify shifting the learning target from full solution prediction to early-to-final consistency estimation.

CA ↑ (BKS 98627.99) SC ↓ (BKS 123.37) WA ↓ (BKS 706.86) IP ↓ (BKS 11.72) Method

Obj

Gap

Obj

Gap

Obj

Gap

Obj

Gap

Gurobi (3600s) Gurobi (1000s)

98448.84 97311.69

179.15 123.37 1316.30 123.64

0.00 706.86 0.27 707.36

0.00 11.72 0.50 13.77

0.00 2.05

ND PS Apollo

94340.63 97906.20 98083.79

4287.36 123.62 721.79 123.60 544.20 123.56

0.25 707.10 0.23 707.09 0.19 707.06

0.24 14.15 0.23 12.08 0.20 11.97

2.43 0.36 0.25

EnCore-ND EnCore-PS EnCore-Apollo

97847.92 98627.99 98491.40

780.70 123.58 0.00 123.50 136.59 123.57

0.21 707.03 0.13 706.98 0.20 707.03

0.17 13.47 0.12 11.95 0.17 11.82

1.75 0.23 0.10

Best Gap Reduction

721.79

100.0%

43.5%

0.10

0.11

47.8%

0.15

60.0%

Table 1: Main results with Gurobi. Obj is the final feasible objective and Gap is the absolute gap to the in-study BKS. Arrows indicate the preferred objective direction; bold marks the best learning-based result or a positive gap reduction. Gurobi EnCore-ND

Primal gap to BKS (%)

EnCore-PS

100 0

ND EnCore-Apollo

SC 103 102 101 100

101

200 400 600 800 1000 Solver time (s)

0.24 0.18 0.12 800

Apollo EnCore starts at 60 s

WA 10 900

1000

100

200 400 600 800 1000 Solver time (s)

0.09 0.06 0.03 800

103

102 900

1000

101 100

200 400 600 800 1000 Solver time (s)

Figure 3: Average primal gap to the BKS versus time under a 1,000-second time limit. EnCore starts after the early solution collection. Each curve is shown only after all test instances have obtained a feasible solution.

Experiments

Experimental Setup Benchmarks. Following Liu et al. (2025), we evaluate EnCore on Combinatorial Auctions (CA) (Leyton-Brown, Pearson, and Shoham 2000), Set Covering (SC) (Balas and Ho 1980), Workload Apportionment (WA), and Item Placement (IP) (Gasse et al. 2022). These benchmarks cover MILPs with diverse scales and structures. For zero-shot transfer, we additionally evaluate on the eleven-instance MIPLIB IIS benchmark considered by Liu et al. (2025), following the IIS setting introduced by Wang et al. (2024). The benchmark is a structurally related subset of MIPLIB centered on binary setcovering formulations, while exhibiting substantial variation in problem size, constraint structure, and sparsity. Dataset sources and instance statistics are reported in App. D. Baselines and variants. We compare against the native MILP solver and three prediction-guided search pipelines: Neural Diving (ND) (Nair et al. 2020), Predict-and-Search (PS) (Han et al. 2023), and Apollo-MILP (Apollo) (Liu et al. 2025). For each learning-based pipeline, we replace its original solution predictor with our consistency predictor while keeping other configurations; the resulting methods are denoted as EnCore-ND, EnCore-PS, and EnCore-Apollo. This design isolates the contribution of the proposed consistency prediction mechanism from downstream solver policies.

Evaluation protocol. All methods use a 1,000-second endto-end budget that includes early-solution collection, inference, and downstream solving. For a fair comparison, EnCore does not use the collected early solutions to warm-start the final search. For EnCore-Apollo, early solutions are collected once and reused across correction rounds. We report the final feasible objective, its absolute gap to the best-known solution (BKS), and the gap reduction relative to the corresponding solution predictor baseline. The BKS is the best mean final objective among all evaluated methods, including the 3,600-second Gurobi reference.

Main Results Table 1 reports the final objectives under Gurobi. Substituting EnCore for the original solution predictor reduces the final gap in 11 of the 12 settings (4 benchmarks × 3 search mechanisms), by 38.7%, 56.9%, and 36.2% on average for ND, PS, and Apollo, respectively. The sole exception is Apollo on SC, where the gap grows marginally from 0.19 to 0.20. On CA, EnCore-PS finds the best solution among all compared methods, improving on the 3,600-second Gurobi incumbent by 179.15 within a 1,000-second budget, and EnCore-Apollo also exceeds this reference (+42.56). On SC, WA, and IP, where the baselines have already been close to the BKS, EnCore still removes 43.5%, 47.8%, and 60.0% of the best baseline gap. Remarkably, EnCore remains effective even

Zero-shot cross-solver transfer to SCIP. We directly apply the Gurobi-trained checkpoint to SCIP-generated trajectories without retraining or model selection. As shown in Table 2, EnCore-ND consistently improves all four benchmarks, reducing the absolute gap by 53.6%, 22.0%, 33.1%, and 50.7% on CA, SC, WA, and IP, respectively. Moreover, EnCore-PS also improves performance on CA, WA, and IP. These results show that the learned consistency representation remains informative under a solver shift. CA ↑

SC ↓

WA ↓

IP ↓

SCIP(1000s)

94701.91

127.99

709.06

23.49

ND PS Apollo

94342.12 97097.46 97335.92

125.14 125.03 124.97

707.33 708.49 708.43

18.19 17.02 16.21

EnCore-ND EnCore-PS EnCore-Apollo

96641.43 97707.64 97293.51

124.75 125.04 125.38

707.24 708.22 707.91

14.91 16.13 17.39

53.6%

22.0%

33.1%

50.7%

Method

Best Gap Reduction

Table 2: Zero-shot transfer from Gurobi to SCIP under a 1,000-second total budget. Bold marks the best result. The BKS is presented in Table 1.

alone improves most benchmarks. However, the consistency target contributes substantially larger gains than the extra features themselves, indicating that reformulating the prediction objective is empirically effective. Finally, introducing earlysolution ensembling further improves CA, SC, and WA, delivering a smaller but complementary robustness benefit. CA ↑

Variant

Ablation Study Table 3 presents a cumulative ablation of the Predict-andSearch pipeline. Adding early solutions as input features

IP ↓

Table 3: Ablation study of key components under the Predictand-Search pipeline. Average objective values are reported. Trade-off Analysis of Early-Solution Collection Figure 4 varies the collection budget Tmax under the fixed 1,000second total budget. Across all benchmarks, performance peaks at a small Tmax : a short collection stage suffices to obtain informative early solutions, whereas larger Tmax leads to diminishing improvements in early-solution quality while consuming more time budget for downstream search. (a) CA ↑ 98640 98560 98480

707.6 707.4 707.2 707.0 0

124.8 124.4 124.0 123.6

500

(b) SC ↓

(d) IP ↓ 13.2 12.8 12.4 12.0 1000 0

500

1000

Early Solution collection-time cap Tmax (s)

Figure 4: Average final objective of EnCore-PS under different maximum early-solution collection times.

6 Zero-shot cross-family transfer to MIPLIB IIS. We train our EnCore models on the mixed SC, CA, WA, and IP instances, and evaluate them on the eleven unseen MIPLIB instances used by Liu et al. (2025) under the same 1,000second budget. The clearest advantage of EnCore is about feasibility. With ND, the original predictor finds feasible solutions on only 3 of 11 instances, while EnCore-ND succeeds on all 11. With PS and Apollo, where the baselines already achieve full feasibility, EnCore still yields consistent improvements in mean objectives. These results demonstrate that our proposed EnCore exhibits greater reliability and robustness than existing methods when faced with substantial distribution shifts. Details are provided in Appendix D.

SC ↓ WA ↓

Predict-and-Search 97906.20 123.60 707.09 12.08 + early solution as feature 98318.22 123.55 707.04 12.17 + consistency prediction target 98616.65 123.53 706.99 11.90 + early-solution ensemble 98627.99 123.50 706.98 11.95

Average objective

under weak easy-to-final consistency. On CA, where 19% of the variables flip between the early and final solutions, the consistency predictor still captures the learnable patterns underlying these flips, allowing EnCore to achieve a 100% average gap reduction. Figure 3 traces the primal gap over time. Although the best EnCore variant on each benchmark begins predictionguided search only after the 60-second collection stage, it catches up with its corresponding baseline at 71, 241, 164, and 93 seconds of end-to-end runtime on CA, SC, WA, and IP, respectively, and remains ahead thereafter. EnCore-PS overtakes all competitors early on CA and attains the lowest gap on SC and WA, with the late-stage separation shown in the insets. On IP, EnCore-Apollo exhibits the strongest latestage anytime performance and widens its lead over time.

Conclusion

In this paper, we presented EnCore, a solver-informed paradigm that accelerates MILP solving via early-to-final solution consistency prediction. Instead of constructing complete solutions from the static MILP formulation alone, EnCore conditions on feasible solutions collected during the early search stage and identifies the assignments likely to persist under a full solving budget. The resulting consistency scores integrate seamlessly with diverse prediction-guided search frameworks, and ensembling multiple early solutions further enhances robustness. Theoretically, we showed that conditioning on early solutions yields a strict gain in achievable prediction accuracy, and that this gain survives finite-sample model selection at a cost discounted by correction sparsity. Empirically, EnCore consistently improves prediction-guided solving across four benchmarks and generalizes zero-shot to unseen solvers and problem families. Future work will explore jointly learning consistency prediction and downstream search policies.

References Achterberg, T. 2009. SCIP: Solving Constraint Integer Programs. Mathematical Programming Computation, 1(1): 1– 41. Balas, E.; and Ho, A. 1980. Set Covering Algorithms Using Cutting Planes, Heuristics, and Subgradient Optimization: A Computational Study, volume 12, 37–60. Springer Berlin Heidelberg. Bengio, Y.; Lodi, A.; and Prouvost, A. 2021. Machine Learning for Combinatorial Optimization: A Methodological Tour d’Horizon. European Journal of Operational Research, 290(2): 405–421. Boucheron, S.; Lugosi, G.; and Massart, P. 2013. Concentration Inequalities: A Nonasymptotic Theory of Independence. Oxford: Oxford University Press. Crainic, T. G. 2000. Service network design in freight transportation. European Journal of Operational Research, 122(2): 272–288. Devroye, L.; Györfi, L.; and Lugosi, G. 1996. A Probabilistic Theory of Pattern Recognition, volume 31 of Stochastic Modelling and Applied Probability. New York: Springer. Ferber, A.; Wilder, B.; Dilkina, B.; and Tambe, M. 2020. MIPaaL: Mixed Integer Program as a Layer. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, 1504–1511. New York, NY, USA. Floudas, C. A.; and Lin, X. 2005. Mixed Integer Linear Programming in Process Scheduling: Modeling, Algorithms, and Applications. Annals of Operations Research, 139(1): 131–162. Gasse, M.; Bowly, S.; Cappart, Q.; Charfreitag, J.; Charlin, L.; Chételat, D.; Chmiela, A.; Dumouchelle, J.; Gleixner, A. M.; Kazachkov, A. M.; Khalil, E. B.; Lichocki, P.; Lodi, A.; Lubin, M.; Maddison, C. J.; Morris, C.; Papageorgiou, D. J.; Parjadis, A.; Pokutta, S.; Prouvost, A.; Scavuzzo, L.; Zarpellon, G.; et al. 2022. The Machine Learning for Combinatorial Optimization Competition (ML4CO): Results and Insights. In Proceedings of the NeurIPS 2021 Competitions and Demonstrations Track, volume 176 of Proceedings of Machine Learning Research, 220–231. Gasse, M.; Chételat, D.; Ferroni, N.; Charlin, L.; and Lodi, A. 2019. Exact Combinatorial Optimization with Graph Convolutional Neural Networks. In Advances in Neural Information Processing Systems 32, volume 32. Vancouver, BC, Canada. Gleixner, A.; Hendel, G.; Gamrath, G.; Achterberg, T.; Bastubbe, M.; Berthold, T.; Christophel, P. M.; Jarck, K.; Koch, T.; Linderoth, J.; Lübbecke, M.; Mittelmann, H. D.; Ozyurt, D.; Ralphs, T. K.; Salvagnin, D.; and Shinano, Y. 2021. MIPLIB 2017: Data-Driven Compilation of the 6th MixedInteger Programming Library. Mathematical Programming Computation, 13(3): 443–490. Gupta, P.; Gasse, M.; Khalil, E. B.; Mudigonda, P.; Lodi, A.; and Bengio, Y. 2020. Hybrid Models for Learning to Branch. In Advances in Neural Information Processing Systems 33, volume 33, 18087–18097. Virtual. Gurobi Optimization, LLC. 2026. Gurobi Optimizer Reference Manual.

Han, Q.; Yang, L.; Chen, Q.; Zhou, X.; Zhang, D.; Wang, A.; Sun, R.; and Luo, X. 2023. A GNN-Guided Predict-andSearch Framework for Mixed-Integer Linear Programming. In The 11th International Conference on Learning Representations. Kigali, Rwanda. Hutter, F.; Hoos, H. H.; Leyton-Brown, K.; and Stützle, T. 2009. ParamILS: An Automatic Algorithm Configuration Framework. Journal of Artificial Intelligence Research, 36: 267–306. Kallenberg, O. 2021. Foundations of Modern Probability. Cham: Springer, 3 edition. Karp, R. M. 1972. Reducibility among combinatorial problems. In Proceedings of the Symposium on the Complexity of Computer Computations, 85–103. Yorktown Heights, New York, USA. Khalil, E. B.; Le Bodic, P.; Song, L.; Nemhauser, G.; and Dilkina, B. 2016. Learning to Branch in Mixed Integer Programming. In Proceedings of the Thirtieth AAAI Conference on Artificial Intelligence, volume 30, 724–731. Phoenix, AZ, USA. Land, A. H.; and Doig, A. G. 1960. An Automatic Method of Solving Discrete Programming Problems. Econometrica, 28(3): 497–520. Laporte, G. 1992. The vehicle routing problem: An overview of exact and approximate algorithms. European Journal of Operational Research, 59(3): 345–358. Leyton-Brown, K.; Pearson, M.; and Shoham, Y. 2000. Towards a Universal Test Suite for Combinatorial Auction Algorithms. In Proceedings of the 2nd ACM Conference on Electronic Commerce, 66–76. Minneapolis, MN, USA. Li, H.; Yuan, H.; Zhang, H.; Lin, J.; Ge, D.; Wang, M.; and Ye, Y. 2026. FMIP: Joint Continuous-Integer Flow for MixedInteger Linear Programming. In The 14th International Conference on Learning Representations. Rio de Janeiro, Brazil. Li, S.; Ouyang, W.; Paulus, M. B.; and Wu, C. 2023. Learning to Configure Separators in Branch-and-Cut. In Advances in Neural Information Processing Systems, volume 36, 60021– 60034. New Orleans, LA, USA. Ling, H.; Wang, Z.; and Wang, J. 2024. Learning to Stop Cut Generation for Efficient Mixed-Integer Linear Programming. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 38, 20759–20767. Vancouver, BC, Canada. Liu, H.; Wang, J.; Geng, Z.; Li, X.; Zong, Y.; Zhu, F.; Hao, J.; and Wu, F. 2025. Apollo-MILP: An Alternating PredictionCorrection Neural Solving Framework for Mixed-Integer Linear Programming. In The 13th International Conference on Learning Representations. Singapore. Nair, V.; Bartunov, S.; Gimeno, F.; von Glehn, I.; Lichocki, P.; Lobov, I.; O’Donoghue, B.; Sonnerat, N.; Tjandraatmadja, C.; Wang, P.; Addanki, R.; Hapuarachchi, T.; Keck, T.; Keeling, J.; Kohli, P.; Ktena, I.; Li, Y.; Vinyals, O.; and Zwols, Y. 2020. Solving Mixed Integer Programs Using Neural Networks. arXiv:2012.13349. Nam, S.; and Logendran, R. 1992. Aggregate production planning — A survey of models and methodologies. European Journal of Operational Research, 61(3): 255–272.

Padberg, M.; and Rinaldi, G. 1991. A Branch-and-Cut Algorithm for the Resolution of Large-Scale Symmetric Traveling Salesman Problems. SIAM Review, 33(1): 60–100. Pu, T.; Li, J.; Gao, Y.; Liu, S.; Geng, Z.; Liu, H.; Chen, C.; and Fan, C. 2026. CoCo-MILP: Inter-Variable Contrastive and Intra-Constraint Competitive MILP Solution Prediction. In Proceedings of the AAAI Conference on Artificial Intelligence, 24882–24890. Singapore. Salvagnin, D. 2016. Detecting Semantic Groups in MIP Models. In Integration of AI and OR Techniques in Constraint Programming, 329–341. Strang, P.; Alès, Z.; Bissuel, C.; Juan, O.; Kedad-Sidhoum, S.; and Rachelson, E. 2026. Planning in Branch-and-Bound: Model-Based Reinforcement Learning for Exact Combinatorial Optimization. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 40, 25627–25635. Tang, Y.; Agrawal, S.; and Faenza, Y. 2020. Reinforcement Learning for Integer Programming: Learning to Cut. In Proceedings of the 37th International Conference on Machine Learning, volume 119, 9367–9376. Vienna, Austria. Wang, H. P.; Liu, J.; Chen, X.; Wang, X.; Li, P.; and Yin, W. 2024. DIG-MILP: A Deep Instance Generator for MixedInteger Linear Programming with Feasibility Guarantee. Transactions on Machine Learning Research. Wang, R.; Li, X.; and Wang, M. 2026. MILPnet: A MultiScale Architecture with Geometric Feature Sequence Representations for Advancing MILP Problems. In The 14th International Conference on Learning Representations.

60 40 20 0

60 40 20 0 0

50 100 150 200 400 600 800 1000

Time (s)

20 0 0

50 100 150 200 400 600 800 1000

Time (s)

Number of Instances

Figure 5: Primal gap versus solving time for Gurobi.

2 0

Flipped Binary Variables (%)

Mean = 0.70% Median = 0.72%

Number of Instances

20 15 10 5 0

Flipped Binary Variables (%)

(b) SC Mean = 4.37% Median = 1.80%

Flipped Binary Variables (%)

Mean = 17.68% Median = 17.71%

40 20 0

50 100 150 200 400 600 800 1000

Flipped Binary Variables (%)

100 80 60 40 20 0

60 40 20 0 0

50 100 150 200 400 600 800 1000

Time (s)

50 100 150 200 400 600 800 1000

Time (s)

(b) SC Primal Gap to BKS (%)

Primal Gap to BKS (%)

(a) CA 80

Flipped Binary Variables (%)

10.0 7.5 5.0 2.5

Number of Instances

0.0

Mean = 15.24% Median = 15.81%

12.5

Flipped Binary Variables (%)

(d) IP

Mean Mean ± Std Knee = 60 s

100 80 60 40 20 0 0

3 2 1 0

50 100 150 200 400 600 800 1000

Mean = 4.12% Median = 1.13%

Time (s)

Flipped Binary Variables (%)

(b) Aggregate flip distribution.

50 Mean = 4.12%

40 30 20

17.74% 15.33%

10 4.67%

3 os am

ds pc

3.94%

2.02%

-sc

a gl

ou m ey

1.13%

-c

c s-h

1 10

0.17%

0.13%

pk sc

0 st

0.11%

ci 4s pj sc

0.10%

0.00%

pl sc

-g iis

-c

MIPLIB Instance

Time (s)

100

Flipped Binary Variables (%)

Figure 9: Solver behavior on the MIPLIB IIS dataset with Gurobi.

(d) IP

Primal Gap to BKS (%)

0 15

Figure 6: Early-to-final flip distributions for Gurobi; red and green dashed lines denote the mean and median. 100

(a) Primal gap over time.

(a) CA 25

(b) SC

(d) IP

Figure 8: Early-to-final flip distributions for SCIP; red and green dashed lines denote the mean and median.

Time (s)

100

50 100 150 200 400 600 800 1000

Mean = 19.87% Median = 20.83%

Flipped Binary Variables (%)

Mean = 4.21% Median = 2.50%

(b) SC Primal Gap to BKS (%)

Primal Gap to BKS (%)

(a) CA

Mean = 0.76% Median = 0.76%

(a) CA

Time (s)

100

50 100 150 200 400 600 800 1000

100

Primal Gap to Run Best (%)

Number of Instances

Mean = 21.76% Median = 21.90%

Number of Instances

100

Primal Gap to BKS (%)

Figures 5 and 7 show the primal gap versus solving time for Gurobi and SCIP on the four benchmarks. Figures 6 and 8 show the distributions of integer variables flipped between the early and full-budget solutions. Figure 9 reports the same analyses on MIPLIB.

Number of Instances

Solver Behavior

Number of Instances

100 80 60 40 20 0 0

50 100 150 200 400 600 800 1000

Time (s)

(d) IP

Figure 7: Primal gap versus solving time for SCIP.

Model Structure

Graph Representation. We follow the MILP representation of Han et al. (2023), adding only one dimension to each variable node to encode its early-solution value. Detailed feature definitions are provided in Table 4. Network architecture. We strictly follow the network architecture proposed by Han et al. (2023), modifying only the variable encoder to accommodate the additional earlysolution feature. Each variable node is represented by 19 input features, consisting of the original 18 features and the early-solution value introduced by our method. Constraint nodes and edges have 4 and 1 input features, respectively. All node representations are embedded into a 64-dimensional latent space. For an input feature vector of dimension d, the node encoder is defined as Encd (x) = ReLU(W2 ReLU(W1 LN(x))) , where W1 : Rd → R64 and W2 : R64 → R64 . Accordingly, the variable and constraint encoders use Enc19 and Enc4 ,

Index Feature Variable-node features 0 Objective 1 Variable coefficient 2 Variable degree 3 4 5

Maximum coefficient Minimum coefficient Variable type

6–17 Position embedding 18 Early-solution value

Description Normalized objective coefficient. Average variable coefficient across all constraints. Degree of the variable node in the bipartite graph. Maximum variable coefficient across all constraints. Minimum variable coefficient across all constraints. Indicator of whether the variable is integer. Binary encoding of the variable’s order among all variables. Value of the variable in the early solution collected as described in Appendix C.

Constraint-node features 0 Constraint Average of the nonzero coefficients in coefficient the constraint. 1 Constraint degree Degree of the constraint node in the bipartite graph. 2 Bias Normalized right-hand side of the constraint. 3 Sense Sense of the constraint. Edge features 0 Coefficient

Coefficient connecting the constraint and variable nodes.

Table 4: Graph features for input. respectively, while the scalar edge feature is normalized by LayerNorm(1). The model performs two rounds of bidirectional message passing between variable and constraint nodes: V → C → V → C → V. For a directed edge j → i, the message is computed as mij = WM ReLU[LN(WL hi + WE eij + WR hj )] , where hi and hj denote the target and source node representations, and eij is the encoded edge feature. Incoming messages are summed and normalized:   X m̄i = LN mij  . j∈N (i)

The target node representation is then updated by h′i = WU,2 ReLU(WU,1 [hi ∥ m̄i ]) , where ∥ denotes concatenation, WU,1 : R128 → R64 , and WU,2 : R64 → R64 . Finally, each variable representation is mapped to a scalar logit: ⊤ ℓv = wout ReLU(Wout hv ) , where the final output layer has no bias. The predicted probability σ(ℓv ) indicates whether the early-solution value of variable v is consistent with its value in the full-budget reference solution.

Early Solution Collection Details

We first quantify the progress of the dual bound using the average dual-gap decay rate, and then describe the earlysolution collection procedure in Algorithm 1. Specifically, for the k-th monitoring interval, the decay rate is defined as rk =

|Dka − Dkb | , tkb − tka

(2)

where ka and kb denote the observations at the beginning and end of the interval, respectively; Dka and Dkb are the corresponding dual gaps; and tka and tkb are their wall-clock times. Thus, rk measures the average absolute reduction in the dual gap per unit time during the k-th monitoring interval. Algorithm 1: Early Solution Collection Require: MILP instance G, minimum probing time Tmin , maximum probing time Tmax , window size w, decay threshold ϵ Ensure: Early solution X ES and probing time tprobe 1: Initialize an empty callback record list R. 2: Start a probing solver run on G. 3: while the solver is running do 4: if a MIPSOL callback occurs at time tk then 5: Record the explored solution Xk and relative dual gap Dk . 6: Append (tk , Xk , Dk ) to R. 7: if |R| ≥ w then 8: Let ka and kb be the first and last records in the latest w-record window. 9: Compute rk by Eq. (2). 10: if tk ≥ Tmin and rk < ϵ then 11: Stop the probing run. 12: end if 13: end if 14: end if 15: if the elapsed probing time reaches Tmax then 16: Stop the probing run. 17: end if 18: end while 19: Set X ES to the last explored solution during probing. 20: Set tprobe to the elapsed probing time. 21: return X ES and tprobe .

Experimental Details

Benchmark Details We evaluate on the four benchmark families used by ApolloMILP (Liu et al. 2025). Combinatorial Auctions (CA) instances follow the CATS generator (Leyton-Brown, Pearson, and Shoham 2000), Set Covering (SC) follows the classical construction of Balas and Ho (1980), and Workload Apportionment (WA) and Item Placement (IP) come from ML4CO (Gasse et al. 2022). CA is a maximization problem, whereas SC, WA, and IP are minimization problems. Table 5 reports average instance sizes. We use 240 training, 60 validation, and 100 testing instances, following the settings in Han et al.

Constraint Number 2590.33 3000 195 64306 Variable Number 1500 5000 1083 61000 Binary Variables Number 1500 5000 1050 1000 Continuous Variables Number 0 0 33 60000 Integer Variables Number 0 0 0 0 Table 5: Statistical information of the benchmark instances. (2023) and Liu et al. (2025). We train all models with 500 epochs and a learning rate of 0.001.

Protocols for Combining EnCore with Downstream Search Methods This section details how EnCore is integrated with the three downstream methods evaluated in Section 5. Following the evaluation protocol proposed by Han et al. (2023) and Liu et al. (2025), every run is given a total wall-clock budget of 1000s, and all overhead introduced by EnCore—including the time spent collecting the early solution—is counted against this budget. Let tES denote the time used to obtain the early solution xES on a given instance. Predict-and-Search and Neural Diving. Both methods consume the predictor output once, before the final solve, so the integration requires no change to their search procedures. Given an instance, we run the solver on the original MILP for tES seconds and record the early solution xES together with the associated solving-process features. The consistency predictor then assigns each binary variable a score p̄i , the estimated probability that xES agrees with the i reference solution. Since both methods specify their fixing budgets separately for the two binary values, selection is value-conditioned: we rank the variables with xES = 0 by i p̄i and retain the top k0 , and independently retain the top k1 among those with xES = 1. On CA under Predict-andi Search, for example, the setting (k0 , k1 , ∆) = (600, 0, 20) selects the 600 zero-valued variables with the highest consistency scores and no one-valued ones, and the resulting partial assignment {xi = 0} defines the center of a search neighborhood of radius ∆ = 20. Neural Diving consumes the selection in the same way except that the assignments are imposed as hard fixings and the reduced MILP is solved. In both cases the solver is restarted from scratch with a budget of 1000 − tES seconds, so that the end-to-end wall-clock time matches that of the baselines. Apollo-MILP. Apollo-MILP performs four prediction– correction rounds with per-round budgets of 100s, 100s, 200s, and 600s. Extending the protocol above directly would require collecting a fresh early solution at the start of every round, which is too costly: the collection phase alone would consume a non-negligible fraction of the 100s rounds. We instead adopt the following assumption. Since Apollo-MILP realizes its fixing step by imposing bound constraints on the selected variables rather than by restructuring the problem, we assume that the early search behavior of the bounded problem stays close to that of the original instance, so that

the early values realized on the original instance remain informative anchors throughout the iterative process. Under this assumption, we collect the early solution only once, on the original instance and in exactly the same way as above, and reuse xES and its consistency scores in the prediction step of every round. Each round applies the same value-conditioned selection, restricted to the variables not yet fixed in earlier rounds. To keep the total budget unchanged, the collection cost tES is charged to the final round, whose budget becomes 600 − tES seconds, while the first three rounds are left untouched; the overall wall-clock time thus remains 1000s.

Evaluation Protocol Hardware. All experiments are conducted on a Linux server (Ubuntu 20.04.6 LTS) equipped with an Intel Xeon Platinum 8124M CPU (72 logical cores, 3.50GHz), 256GB of RAM, and an NVIDIA GeForce RTX 3090 GPU (24GB). Neural network training and inference run on the GPU, while all solver calls run on the CPU. The solver configuration (e.g., thread count) is kept identical across all methods on each benchmark. Metrics. We report the final feasible objective and its absolute primal gap gapabs = |z − zBKS |,

(3)

where z is the mean final objective of a method and zBKS is the best-known solution (BKS) on the corresponding benchmark. The BKS is defined as the best mean final objective observed among all evaluated methods, including a 3,600second Gurobi run: Gurobi provides this reference on SC, WA, and IP, while EnCore-PS provides it on CA. The BKS is thus an in-study best-observed reference rather than an independently certified optimum. To quantify the benefit of our method over a downstream baseline, we additionally report the relative gap reduction Improvement =

gapbase − gapours × 100%. gapbase

(4)

Zero-Shot MIPLIB IIS Transfer Protocol The IIS setting traces to Wang et al. (2024), which formed a two-instance dataset from MIPLIB: iis-glass-cov was used for training and iis-hc-cov for testing. ApolloMILP later expanded this setting into an eleven-instance subset of MIPLIB 2017 (Gleixner et al. 2021; Liu et al. 2025). Specifically, it selected instances according to similarities measured with 100 human-designed structural features, then discarded instances whose presolving exceeded 300 seconds or whose inference exceeded GPU memory. Apollo-MILP used eight of the resulting instances for training and three— ramos3, scpj4scip, and scpl4—for testing. In contrast, we use all eleven instances exclusively as a zero-shot test set. The instances range from 214 to 200,000 variables and from 32,805 to 2,000,000 nonzero coefficients. In this experiment, both predictors are trained on the union of the CA, SC, WA, and IP training sets. They are validated on the union of the corresponding validation sets, and checkpoint selection uses only this combined validation set. The

selected checkpoints are frozen before evaluation on IIS. No MIPLIB instance is used for training, validation, fine-tuning, or checkpoint selection. Table 6 gives the overall performance while Table 7 gives the instance-level results. The PS improvement is mainly attributable to ramos3 and ex1010-pi, partially offset by scpk4. Apollo+Ours improves ex1010-pi, degrades on ramos3, and ties on the remaining instances. Predictor

Downstream

Mean Obj ↓

Feas.

Origin EnCore

Neural Diving Neural Diving

243.00† 173.45

3/11 11/11

Origin EnCore

Predict-and-Search Predict-and-Search

172.00 171.73

11/11 11/11

Origin EnCore

Apollo-MILP Apollo-MILP

172.91 172.82

11/11 11/11

Table 6: Zero-shot cross-family transfer to the MIPLIB IIS subset. Mean Obj is the mean of the per-instance final objectives. Feas. reports the number of instances with a finite feasible objective. † means computed only over its three feasible instances.

ND Instance

GCN

PS Ours

GCN

Apollo-MILP Ours

GCN

Ours

ex1010-pi – 242.00 237.00 236.00 fast0507 – 174.00 174.00 174.00 glass-sc – 23.00 23.00 23.00 iis-glass-cov – 21.00 21.00 21.00 iis-hc-cov – 17.00 17.00 17.00 ramos3 – 235.00 229.00 226.00 scpj4scip 132.00 132.00 132.00 132.00 scpk4 328.00 330.00 326.00 327.00 scpl4 269.00 269.00 269.00 269.00 seymour – 423.00 423.00 423.00 v150d30-2hopcds – 42.00 41.00 41.00

241.00 174.00 23.00 21.00 17.00 231.00 132.00 330.00 269.00 423.00 41.00

238.00 174.00 23.00 21.00 17.00 233.00 132.00 330.00 269.00 423.00 41.00

Table 7: Final objectives on the eleven MIPLIB IIS instances. All instances are minimization problems. Bold marks the lowest objective in each row; a dash indicates that no finite feasible solution was found.

Hyperparameter Settings

Early Solution Collection For early-solution collection, we set the sliding-window size to w = 5, the dual-gap decay threshold to ϵ = 0.01%/s, the minimum probing time to Tmin = 20 seconds, and the maximum probing time to Tmax = 60 seconds. We use an ensemble size of K = 3 for SC, CA and WA and use K = 2 for IP.

Neighborhood Search Parameters We report the hyper-parameters we used in our study. EnCore and baselines use the same hyper-parameters. Table 8 shows the hyperparameters of Predict-and-Search and

EnCore-PS, and the Neural Diving and EnCore-ND uses the same (k0 , k1 ). Table 9 shows the Apollo hyperparameters for baseline Apollo-MILP and EnCore-apollo. Benchmark

PS+Gurobi (600, 0, 20) (2000, 0, 100) (400, 5, 10) (0, 500, 10) PS+SCIP (400, 0, 20) (2000, 0, 100) (400, 5, 1) (0, 600, 5)

Table 8: The partial solution size parameters (k0 , k1 ) and neighborhood parameter ∆.

Iteration 1 (400, 0, 60) (1000, 0, 200) (100, 20, 50) (20, 200, 100) Iteration 2 (200, 0, 30) (500, 0, 100) (40, 15, 20) (10, 100, 50) Iteration 3 (100, 0, 15) (250, 0, 50) (20, 15, 10) (10, 5, 5) Iteration 4 (50, 0, 10) (10, 0, 5) (5, 50, 30) (1, 10, 5) (i)

(i)

Table 9: Hyperparameters (k0 , k1 , ∆(i) ) for Apollo-MILP.

Details and Proofs for the Theoretical Analysis

Population Setting and Notation Let D be the probability distribution over all triples consisting of an MILP graph, the early solution produced by the fixed collection procedure, and the full-budget solution. Draw Z = (G, X ES , X ∗ ) ∼ D. Let B(G) be the set of binary-variable indices in G. Assume |B(G)| ≥ 1 almost surely. To represent the per-instance variable average compactly, conditional on Z, choose J from B(G), giving every index probability 1/|B(G)|. We use J for this random index in population quantities and r for a deterministic summation index. Fix the message-passing depth L. For r ∈ B(G), let Rr denote the part of the static graph and its features within L ES message-passing steps of variable r, and let XR denote the r early assignments attached to the variable nodes in this local input. For any assignment predictor h, define X 1 1{hr ̸= x∗r }, ℓh (Z) := |B(G)| r∈B(G)

R(h) := E[ℓh (Z)] = Pr(hJ ̸= x∗J ), Acc(h) := 1 − R(h). Here a direct prediction hr is a function of Rr , whereas ES an induced consistency prediction may also use XR . The r corresponding posteriors are ηG := Pr(x∗J = 1 | RJ ), ES ηES := Pr(x∗J = 1 | RJ , XR ). J

Their Bayes errors are bG := E[min{ηG , 1 − ηG }], bES := E[min{ηES , 1 − ηES }].

Thus A∗sol = 1−bG is the best population accuracy attainable from the static local input, and 1 − bES is the best assignment accuracy attainable when the early assignments are also observed.

Proof of Theorem 1

Consistency–Assignment Equivalence

The inequality is strict if there exists an event E ∈ σ(RJ ) with Pr(E) > 0 such that, almost surely on E,

Lemma 1 (label–assignment equivalence). For r ∈ B(G), let yr = 1{xES = x∗r }. Given a binary consistency r ES prediction ybr , retain xr when ybr = 1 and use its binary complement otherwise. The induced assignment x br satisfies 1{b xr ̸= x∗r } = 1{b yr ̸= yr }. Equivalently, for a consistency rule a, define

(5)

r∈B(G)

Proof. For binary values, x∗r = xES ⊕ (1 − yr ). XOR r by the same value preserves equality, so ha (G, X ES )r ̸= x∗r exactly when a(G, X ES )r ̸= yr . Averaging over B(G) proves Eq. (5). □ The induced assignment ha is used only to compare variable-level prediction errors. In the downstream search method, assignments predicted inconsistent are left free rather than forcibly complemented, so these identities characterize the variable-level prediction task, while the downstream objective and runtime effects are evaluated empirically. Conversely, every direct solution predictor h defines a consistency rule ah (G, X ES )r := 1{h(Rr ) = xES r }. The induced predictor satisfies hah (G, X ES )r = h(Rr ). Thus consistency prediction with the augmented input can reproduce every direct predictor based on the static local input. ∗ Let YJ = 1{xES J = xJ } and ES ξJ := Pr(YJ = 1 | RJ , XR ). J

For a fixed variable r, ξr denotes the corresponding conditional posterior evaluated at r. Since xES J is observed in the augmented local input, ES ξJ = xES J ηES + (1 − xJ )(1 − ηES ).

Consequently, min{ξJ , 1 − ξJ } = min{ηES , 1 − ηES }, so the Bayes consistency error equals bES . Under variablelevel 0–1 prediction loss, the optimal rule predicts 1 when ξJ ≥ 1/2 and 0 otherwise. By Lemma 1, its induced assignment accuracy is A∗con = 1 − bES . Because the label transformation is bijective after xES J is observed, this is also the Bayes accuracy of a direct assignment predictor using the same augmented input. Theorem 1 therefore compares the information in the augmented input with that in the static input, rather than attributing a Bayes advantage to the consistency relabeling itself.

(6)

Pr(ηES < 12 | RJ ) > 0, Pr(ηES > 21 | RJ ) > 0. In words, the same static local input can be paired with early assignments that make either final value more likely. The tower property gives

ηG = E[ηES | RJ ].

Then X 1 1{a(G, X ES )r ̸= yr }. |B(G)|

A∗con − A∗sol = bG − bES ≥ 0.

Proof of Theorem 1.

ha (G, X ES )r := xES ⊕ (1 − a(G, X ES )r ). r ℓha (Z) =

The population accuracy difference can be written as

(7)

Using min{t, 1 − t} = 21 − |t − 12 |, the Bayes accuracy gap

A∗con − A∗sol = E[|ηES − 21 |]

− E[|ηG − 12 |]. Conditional Jensen’s inequality and Eq. (7) imply E[|ηES − 12 | | RJ ] ≥ |ηG − 12 |. Taking expectations proves Eq. (6). If the refined posterior lies on both sides of 1/2 with positive conditional probability, the conditional absolute-value inequality is strict. The stated positive-probability condition therefore gives A∗con > A∗sol . □ Scope of the comparison. The two optima compared in Theorem 1 range over functions of the same receptive field: A∗sol is the Bayes accuracy given the static local input RJ , and A∗con is the Bayes accuracy given the augmented local ES input (RJ , XR ), at the same fixed depth L. Statements in J the main text about “any solution predictor”, or about what “no predictor on instance features alone” can achieve, therefore quantify over predictors whose decision at a variable is a measurable function of its static local input, however large that function class is; they do not compare against predictors with a strictly larger input scope. In particular, if X ∗ were almost surely determined by the full graph G—for example, a unique optimal solution returned by a deterministic solver—then the whole-graph posterior would be {0, 1}valued, posterior crossing could not occur at the whole-graph level, and a whole-graph predictor would attain perfect accuracy in principle, leaving no information-theoretic gain for the early solution. The strict-gain content of Theorem 1 is thus per locality level: for every fixed message-passing depth, augmenting the local input strictly improves the best achievable accuracy whenever posterior crossing occurs at that depth. This is the operative comparison for the models considered in this line of work, including ours and the baselines, which are all message-passing GNNs whose receptive field is fixed by the architecture and is typically far smaller than the benchmark graphs. Moreover, on practical benchmarks X ∗ retains genuine randomness given G: instances frequently admit multiple optimal or near-optimal solutions, and the incumbent returned under a time budget depends on tie-breaking, thread timing, and budget truncation rather than

on G alone. Under such residual randomness the whole-graph posterior is not {0, 1}-valued either, so posterior crossing— and hence a strict gain from the early solution—can occur at any input scope, including the full graph.

Sparse-Correction Identity The Bayes consistency rule is a∗ (G, X ES )r := 1{ξr ≥ 12 }. Define the event C ∗ := {ξJ < 12 } = {a∗ (G, X ES )J = 0},

Recall ∆B := A∗con − A∗sol , and define the finite-class approximation gap αA := A∗con − A∗A . Then ∆A = ∆B −αA , separating the population information gain from the approximation cost of the finite class. For an instance Z, let sa (Z) be the proportion of binary variables that rule a would flip. Assume every candidate flips at most a proportion s̄ ∈ [0, 1] on average: X 1 sa (Z) := (1 − a(G, X ES )r ), |B(G)| r∈B(G)

max E[sa (Z)] ≤ s̄.

and define

a∈A

p := Pr(YJ = 0),

ρ := Pr(C ∗ ).

If ρ > 0, also define q := Pr(YJ = 0 | C ∗ ). Then q > 1/2, and a∗ satisfies bES = R(ha∗ ) = p − ρ(2q − 1), ∆B := A∗con − A∗sol = bG − p + ρ(2q − 1).

(8)

Thus ∆B > 0 exactly when p < bG + ρ(2q − 1).

(9)

For p ≥ bG and ρ > 0, this is equivalent to q>

Thus ∆A is the population accuracy difference between the best rule available in A and the best predictor based only on the static local input. It can be nonpositive if the finite class does not contain a sufficiently accurate consistency rule. For a predictor that depends on the realized training sample and any training randomness, R(h) and Acc(h) below denote conditional population quantities evaluated on a fresh test triple independent of that training information. Let akeep (G, X ES )r ≡ 1. Then hakeep (G, X ES ) = X ES and R(hakeep ) = p. Define the loss change relative to this baseline by Da (Z) := ℓha (Z) − ℓhakeep (Z),

1 p − bG + . 2 2ρ

b a := 1 D m

If ρ = 0, then bES = p and ∆B = bG − p. Derivation. For any consistency rule a, retaining an assignment has error 1 − YJ , whereas using its binary complement has error YJ . Hence R(ha ) = p + E[(2YJ − 1)(1 − a(G, X ES )J )].

(10)

∗

On C , the conditional probability of YJ = 0 is q. Substituting a∗ into Eq. (10) gives R(ha∗ ) = p + ρ(1 − q) − ρq = p − ρ(2q − 1). Moreover,

q = 1 − E[ξJ | C ∗ ] > 12 . The rule a∗ is the pointwise Bayes classifier for YJ , so Lemma 1 gives R(ha∗ ) = bES . Subtracting this risk from bG proves Eq. (8); the remaining statements follow by rearrangement. If ρ = 0, a∗ retains the early solution almost surely, and its error is p. □

Finite-Sample Analysis The population result ranges over all possible rules. For the finite-sample result, fix before observing the training instances a finite class A of N := |A| ≥ 1 consistency rules whose variable-level decisions use the augmented local inputs defined above. Let Z1 , . . . , Zm be i.i.d. training instances from D. Suppose training and test instances follow D. Define A∗A := max Acc(ha ), a∈A

∆A := A∗A − A∗sol = bG − min R(ha ). a∈A

m X

Da (Zk ).

k=1

Lemma 2 (uniform concentration). Following the finiteclass Bernstein–union-bound argument (Boucheron, Lugosi, and Massart 2013), we apply concentration to instance-level loss differences relative to the always-retain rule. For every δ ∈ (0, 1), with probability at least 1 − δ, b a − E[Da ]| ≤ βES , sup |D a∈A

4 2N 2s̄ 2N ln + ln . m δ 3m δ Proof. A rule differs from the always-retain baseline only on assignments that it marks inconsistent, so |Da (Z)| ≤ sa (Z) ≤ 1, βES :=

Var(Da (Z)) ≤ E[Da (Z)2 ] ≤ E[sa (Z)] ≤ s̄. Thus, when candidate corrections modify only a small proportion of early assignments, their loss differences have lower variance; this is the finite-class sparse-correction effect captured by s̄. Also, |Da (Z) − E[Da ]| ≤ 2. (The sharper bound |Da (Z) − E[Da ]| ≤ sa (Z) + s̄ ≤ 1 + s̄ is available, since |Da (Z)| ≤ sa (Z) and |E[Da ]| ≤ s̄; we retain the simpler constant 2, which only inflates the second-order term and keeps the penalty identical to the main-text statement.) To make the additive constant explicit, start from the two-sided Bernstein inequality in exponential form (Boucheron, Lugosi, and Massart 2013, Theorem 2.10): for a variance proxy σ 2 and centered range c, mϵ2 b Pr |Da − E[Da ]| ≥ ϵ ≤ 2 exp − 2 . 2σ + 2cϵ/3

Bounding the right-hand side by 2e−t gives the quadratic condition mϵ2 ≥ 2σ 2 t+2ctϵ/3, whose positive root satisfies r r ct 2σ 2 t c2 t2 2σ 2 t 2ct ≤ ϵ= + + + , 3m m 9m2 m 3m √ √ √ where the inequality uses a + b ≤ a + b. Substituting 2 the variance proxy p σ = s̄ and the range c = 2 shows that the deviation 2s̄t/m + 4t/(3m) is sufficient; the additive form is conservative relative to the exponential form because of this square-root relaxation. That is, for every fixed a and t > 0, ! r 2s̄t 4t b Pr |Da − E[Da ]| > ≤ 2e−t . + m 3m Setting t = ln(2N/δ) and applying a union bound over A proves the claim. □ Theorem 2 (restated). Let b a minimize the empirical 0–1 risk over A, and write b hES := hba . Let b hsol be any possibly data-dependent direct predictor determined without access to the fresh test triple and whose test-time decisions use only the static local inputs. Define r 2s̄ 2N 8 2N εm (δ) := 2 ln + ln . m δ 3m δ With probability at least 1 − δ, Acc(b hES )−Acc(b hsol ) ≥ ∆A −εm (δ) = ∆B −αA −εm (δ). On the same event, a∈A

Proof. Subtracting the same empirical loss of akeep from every candidate preserves the empirical minimizer. On the event in Lemma 2, choose a∗A ∈ arg min R(ha ). a∈A

Then b ba + βES ≤ D b a∗ + βES E[Dba ] ≤ D A

R(b hES ) ≤ R(ha∗A ) + 2βES (11)

Condition on the realized training sample. Any resulting direct predictor whose test decision uses only the static local input has, by Bayes optimality,

Acc(b hES ) − Acc(b hsol ) = R(b hsol ) − R(b hES ) ≥ ∆A − 2βES . Since εm (δ) = 2βES , the theorem follows. □

Let Λ := ln(2N/δ). The two bounds in Eq. (13) r 2s̄Λ ∆A 8Λ ∆A 2 < , < . m 2 3m 2 Adding the inequalities gives εm (δ) < ∆A ; Theorem 2 then yields the claim. □ Proof. imply

Illustrative Calculations Population gain. The sparse-correction condition can hold even when the early solution itself is less accurate than direct prediction. For example, take p = 0.21,

A∗sol = 0.82,

ρ = 0.12,

q = 0.70.

These numerical values are purely illustrative. The righthand side of Eq. (9) is 1 − A∗sol + ρ(2q − 1) = 0.18 + 0.12(0.40) = 0.228. Thus p = 0.21 < 0.228. The early solution alone has accuracy 0.79, which is below the direct-prediction accuracy 0.82, whereas identifying the inconsistent assignments gives A∗con = 1 − [0.21 − 0.12(0.40)] = 0.838.

Finite-sample bound. For a numerical evaluation of the bound, consider a fixed class of N = 4 correction rules constructed to satisfy E[sa ] ≤ 0.0180. The value s̄ = 0.0180 uses the 1.80% median WA early-to-final difference in Figure 2 as an illustrative correction-budget scale, with E[sa ] ≤ s̄ retained as the stated rule-class assumption. Set the confidence level to 1−δ = 0.9 and take a deliberately conservative illustrative class advantage ∆A = 0.12. Here ln

Because R(ha ) = p + E[Da ], it follows that

Combining Eqs. (11) and (12) gives

then εm (δ) < ∆A . Consequently, with probability at least 1 − δ, Acc(b hES ) > Acc(b hsol ).

2N = ln 80 ≈ 4.382. δ The two terms on the right-hand side of Eq. (13) are

≤ E[Da∗A ] + 2βES .

R(b hsol ) ≥ bG .

(13)

The early-solution input therefore raises the accuracy from 0.82 to 0.838, an improvement of 1.8 percentage points.

R(b hES ) − min R(ha ) ≤ εm (δ).

= 1 − A∗A + 2βES = bG − ∆A + 2βES .

A sufficient sample size. If ∆A > 0 and 32s̄ ln(2N/δ) 16 ln(2N/δ) , , m > max ∆2A 3∆A

(12)

32s̄ ln(2N/δ) ≈ 175.3, ∆2A 16 ln(2N/δ) ≈ 194.8. 3∆A Thus the convenient sufficient condition in Eq. (13) would require m = 195. It is conservative because it separately bounds each of the two terms in εm by ∆A /2. Directly solving εm (0.1) < 0.12, namely r 2(0.0180) 8 2 ln 80 + ln 80 < 0.12, m 3m gives m > 188.19, so the smallest integer sample size satisfying the exact inequality is m = 189.

Consistency requirement on the parameters. The illustrative values above are not jointly free. A rule that flips a proportion sa (Z) of assignments changes the induced assignment on at most that fraction of variables, so its accuracy exceeds the always-retain accuracy by at most the expected flip rate: Acc(ha ) ≤ (1 − p) + E[sa ] ≤ (1 − p) + s̄. Consequently, ∆A = A∗A − A∗sol ≤ bG − p + s̄, and the choice s̄ = 0.0180, ∆A = 0.12 implicitly requires bG − p ≥ 0.102: the Bayes error of static local prediction must exceed the early-solution error rate by roughly ten percentage points or more. On WA, where the 1.80% median early-to-final difference suggests p ≈ 0.018, this amounts to bG ≳ 0.12; that is, even the best predictor using only the static local input must mislabel at least about 12% of the binary variables. This is a substantive assumption about the difficulty of static prediction on the benchmark, not a consequence of the theorems, and we state it so that the example is read as a self-consistent operating point rather than a measured one. Smaller, still positive class advantages remain valid and simply rescale the sufficient sample size, which grows like ∆−2 A in Eq. (13). Each MILP instance is one independent observation, so labels within an instance may be arbitrarily dependent. Accordingly, the penalty depends on m, the number of instances, rather than the total number of variable labels. The result assumes a fixed finite rule class selected by empirical variablelevel 0–1 prediction loss. Its formal scope is finite-rule model selection, providing an abstraction that isolates instance-level estimation and the sparse-correction variance effect. In practice the predictor is a GNN trained by stochastic optimization with validation-based checkpoint selection, a procedure better described as comparing a small, data-dependent set of candidates—the checkpoints and configurations actually evaluated—than as exact empirical risk minimization over a class fixed in advance. Theorem 2 should therefore be read as an idealized account of that selection stage: it isolates the instance-level estimation cost and explains why sparse correction classes are cheap to select among, but it is not a uniform-convergence guarantee for the full GNN function class. Statements in the main text about the gain persisting under finite-sample model selection refer to this abstraction layer. Relation to the deployed algorithm. Two design choices of the deployed system fall outside the formal statements above and are worth flagging explicitly. First, the analysis scores rules by variable-level 0–1 accuracy—equivalently, thresholding the posterior at 1/2 over all binary variables— whereas the deployed method ranks variables by the predicted consistency probability and fixes only the top-ranked ones under a budget chosen by the downstream method. The theorems thus quantify the information available in the scores; how the scores are consumed by a ranking-based fixing rule is evaluated empirically. Second, the posteriors in this appendix condition on a single early solution, whereas at inference time the system averages, per variable, the logits computed from the last K improving solutions after aligning them to the reference solution X ES . Because conditioning

on the tuple of the last K solutions only enlarges the input, Theorem 1 applies verbatim to the ensembled variant, while the specific alignment-and-averaging estimator is a practical variance-reduction heuristic whose effect is evaluated empirically.

Learning Early-to-Final Solution Consistency for MILP Acceleration

Analyze & explore related research in Conceptio

Related documents