Primal Improvement Heuristics

Primal Improvement Heuristics#

Branch and Bound (B&B) certifies optimality by closing the gap between a dual bound (the best lower bound on the objective, from solving relaxations) and a primal bound (the objective of the best feasible incumbent found so far) [Belotti et al., 2013, Land and Doig, 1960]. The dual bound is what proves optimality, but the primal bound is what makes the search fast: every node whose relaxation is worse than the incumbent can be pruned without ever being explored.

This is why finding a good incumbent early matters so much. A strong incumbent discovered at the root front-loads pruning, tightens the optimality gap immediately, and powers reduced-cost variable fixing — so the tree the solver actually walks is dramatically smaller than the one it would walk if it stumbled onto good solutions only late in the search [Achterberg, 2007, Berthold, 2006]. Relaxation solvers do not, on their own, produce integer-feasible points; primal heuristics are the machinery that turns a fractional relaxation solution into a usable incumbent.

discopt ships four families of primal heuristics, all in discopt._relax.primal_heuristics:

  1. Feasibility pump — round the fractional integers, fix them, and re-solve the continuous part; perturb and retry if the rounding is infeasible [Fischetti et al., 2005].

  2. Diving — fix one fractional integer at a time and re-solve, descending the tree along a single branch (fractional or objective-guided rounding).

  3. RINS (Relaxation Induced Neighborhood Search) — fix the integers on which the incumbent and the relaxation agree, then dive on the rest, searching the neighbourhood “between” them [Danna et al., 2005].

  4. Local branching — search the Hamming-radius-\(k\) neighbourhood of a binary incumbent for something better [Fischetti and Lodi, 2003].

A defining property of every one of these heuristics in discopt is that they only ever propose feasible incumbents. They never touch the dual bound, so they cannot weaken the global optimality certificate — any candidate they return is re-verified for integer- and constraint-feasibility before injection, and the B&B tree keeps it only if it strictly improves the incumbent.

import numpy as np

import discopt.modeling as dm

Heuristics inside Model.solve()#

You do not normally call these functions yourself: discopt engages them automatically during B&B. After the root node is solved, a feasibility pump runs on the best relaxation solution; on nonconvex models an integer local search and a scheduled sub-NLP heuristic run at the root and periodically afterwards (controlled by subnlp_frequency / subnlp_max_calls), and a root binary-seed enumeration covers small binary blocks. Every injected point is re-verified before it becomes the incumbent.

Here is a small nonconvex MINLP with bilinear coupling and integer-heavy structure — the regime the sub-NLP heuristic exists for. The subnlp_* counters on the result report its activity, and on this instance they all read 0: the heuristic is gated on not _root_optimum_proven(), and the root solve already produced the optimal incumbent, so there was nothing for it to repair. Zero is the correct reading, not a missing feature — and the cell asserts the reason rather than leaving the counters to be read as failure. (Measured: they stay 0 on this model even at subnlp_frequency=1, subnlp_max_calls=50.)

m = dm.Model("nonconvex_minlp")
y = [m.integer(f"y{i}", lb=0, ub=5) for i in range(3)]
x = [m.continuous(f"x{i}", lb=0, ub=5) for i in range(3)]

# Nonconvex: bilinear product terms in objective and constraint.
m.minimize(
    sum((x[i] - 2.5) ** 2 for i in range(3))
    - x[0] * y[0]
    - x[1] * y[1]
    + sum(y[i] for i in range(3))
)
m.subject_to(x[0] * y[1] + x[1] * y[2] <= 12, name="bilinear")
m.subject_to(sum(y[i] for i in range(3)) >= 4, name="cover")

result = m.solve(time_limit=30)

print(f"Status:    {result.status}")
print(f"Objective: {result.objective:.4f}")
print(f"Nodes:     {result.node_count}")
print()
print(f"subnlp_calls:             {result.subnlp_calls}")
print(f"subnlp_feasible:          {result.subnlp_feasible}")
print(f"subnlp_incumbent_updates: {result.subnlp_incumbent_updates}")

# The counters are 0 because the root solve already found the optimum; pin both
# halves of that statement, so a build where the heuristic DOES fire here fails
# this cell instead of silently changing what the prose describes.
assert result.status == "optimal"
assert result.subnlp_calls == 0, (
    f"subnlp fired {result.subnlp_calls} time(s) -- the prose below says it did not"
)

Status:    optimal
Objective: -20.7400
Nodes:     15

subnlp_calls:             0
subnlp_feasible:          0
subnlp_incumbent_updates: 0

A non-zero subnlp_incumbent_updates would mean a primal heuristic had found a feasible point beating what plain node solving produced, and injected it as the incumbent. Here all three counters are zero and the tree still closes in a handful of nodes — the root relaxation’s own solution was already optimal, which is the case the heuristic schedule deliberately skips.

Calling a heuristic directly#

The heuristics are also exported as standalone functions for experimentation. They share a common shape: given the model and an NLP relaxation solution x_relax (and, for RINS / local branching, an incumbent), they return either a verified (x, objective) pair or None. They take an optional prebuilt NLPEvaluator so the evaluator’s tape is built once and reused across calls. (That tape is POUNCE’s, in the Rust core — JAX is not on the solve path.)

Let’s set up a small knapsack-style MILP and obtain its continuous relaxation solution. The NLPEvaluator treats integer variables as continuous within their bounds, so a single NLP solve from the box midpoint gives us a fractional relaxation point to feed the heuristics.

from discopt._relax.nlp_evaluator import NLPEvaluator
from discopt.solvers.nlp_backend import get_nlp_solver

# A 0/1 knapsack: maximize value subject to a weight budget.
values = [8, 5, 11, 6, 9]
weights = [5, 3, 7, 4, 6]
capacity = 12

m = dm.Model("knapsack")
items = [m.binary(f"x{i}") for i in range(len(values))]
m.maximize(sum(values[i] * items[i] for i in range(len(values))))
m.subject_to(sum(weights[i] * items[i] for i in range(len(values))) <= capacity, name="cap")

# Continuous relaxation solution at the root.
evaluator = NLPEvaluator(m)
lb, ub = evaluator.variable_bounds
backend = get_nlp_solver("auto")
relax = backend(evaluator, 0.5 * (lb + ub), options={"print_level": 0})
x_relax = np.asarray(relax.x)

print("Relaxation solution (fractional):", np.round(x_relax, 3))

Relaxation solution (fractional): [ 1.     1.     0.571 -0.    -0.   ]

Fractional diving#

fractional_diving (a thin wrapper over diving(..., mode="fractional")) repeatedly fixes the most fractional unfixed integer — the one closest to 0.5 — to its nearest integer and re-solves, until every integer is integral or a sub-NLP turns infeasible. It returns a verified (x, objective); note the objective is in minimization sense (discopt minimizes the negated objective for a maximize model).

from discopt._relax.primal_heuristics import fractional_diving

dive = fractional_diving(m, x_relax, evaluator=evaluator)
x_dive, obj_dive = dive
print("Diving incumbent:", np.round(x_dive, 0))
print(f"Objective (min sense): {obj_dive:.1f}  ->  knapsack value: {-obj_dive:.0f}")

Diving incumbent: [0. 1. 1. 0. 0.]
Objective (min sense): -16.0  ->  knapsack value: 16

RINS#

rins needs both an incumbent and a relaxation point. It fixes every integer on which the two agree, then dives on the variables where they disagree — restricting the search to the neighbourhood “between” the incumbent and the relaxation, which is often where better solutions hide. Here we start from a mediocre incumbent (items 0 and 3, value 14) and let RINS find something better.

from discopt._relax.primal_heuristics import rins

x_incumbent = np.array([1.0, 0.0, 0.0, 1.0, 0.0])  # items 0 and 3, value 14
improved = rins(m, x_incumbent, x_relax, evaluator=evaluator)

if improved is not None:
    x_rins, obj_rins = improved
    print("RINS incumbent:", np.round(x_rins, 0))
    print(f"Objective (min sense): {obj_rins:.1f}  ->  knapsack value: {-obj_rins:.0f}")
else:
    print("RINS found no improving point in the agreement neighbourhood.")

RINS incumbent: [ 1.  0.  1. -0.  0.]
Objective (min sense): -19.0  ->  knapsack value: 19

Local branching#

local_branching searches the Hamming-radius-\(k\) neighbourhood of a binary incumbent: it enumerates every flip of up to \(k\) binaries, fixes them, and solves the continuous sub-NLP for each, keeping the best feasible result. This realises the classic local-branching neighbourhood [Fischetti and Lodi, 2003] directly — without a recursive sub-MIP — which keeps it cheap and self-contained for the small/medium binary blocks where it pays off.

from discopt._relax.primal_heuristics import local_branching

x_incumbent = np.array([1.0, 1.0, 0.0, 0.0, 0.0])  # items 0 and 1, value 13
best = local_branching(m, x_incumbent, k=2, evaluator=evaluator)

x_lb, obj_lb = best
print("Starting incumbent value: 13")
print("Local-branching incumbent:", np.round(x_lb, 0))
print(f"Objective (min sense): {obj_lb:.1f}  ->  knapsack value: {-obj_lb:.0f}")

Starting incumbent value: 13
Local-branching incumbent: [1. 1. 0. 1. 0.]
Objective (min sense): -19.0  ->  knapsack value: 19

How each heuristic works, and when discopt triggers it#

Feasibility pump (Fischetti et al. [2005]). Given a relaxation solution with fractional integers, round the integers to their nearest values, fix them (tighten their bounds to that single value), and re-solve the continuous NLP. If the result is integer- and constraint-feasible, return it; otherwise randomly flip a fraction of the integers and retry for a few rounds. The rounding-and-repair loop is what converts a relaxation point into a genuine incumbent. discopt runs the feasibility pump once, right after the root node, on the best relaxation solution in the batch (feasibility_pump(model, x_relax, max_rounds=5, ...) in solver.py).

Diving (fractional_diving, objective_diving). Rather than rounding everything at once, diving fixes a single fractional integer per step and re-solves, descending a single branch of the tree. The "fractional" rule fixes the most fractional variable (closest to 0.5) to its nearest integer; the "objective" rule rounds in the direction the objective gradient prefers (toward the cheaper neighbour). The dive ends when all integers are integral (success) or a sub-NLP becomes infeasible (the dive is abandoned). discopt uses a root LP dive (_root_dive) to manufacture an early incumbent before the main loop begins.

RINS (Danna et al. [2005]). Relaxation Induced Neighborhood Search fixes every integer on which the current incumbent and the current relaxation agree (and which the relaxation reports as integral), then runs a fractional dive on the remaining free integers. The intuition: where the incumbent and the relaxation concur, they are probably right, so the interesting search lives in the disagreement set — a small, focused sub-problem. With full agreement there is no neighbourhood to explore (it returns None); with full disagreement it degenerates to a plain dive. discopt invokes it as an improvement heuristic once an incumbent and a fresh relaxation point are both available.

Local branching (Fischetti and Lodi [2003]). Classic local branching adds a constraint bounding the Hamming distance from a binary incumbent to at most \(k\) and re-solves a sub-MIP. discopt realises the same neighbourhood by direct enumeration: every flip of up to \(k\) binaries is fixed and the continuous sub-NLP is solved, and the best feasible point is kept. This is exact for the neighbourhood and needs no recursive solve, so it stays cheap for the small/medium binary blocks (capped by max_binaries) where local branching is worthwhile. It runs as an improvement step once a binary incumbent exists.

Two further root-level heuristics round out the inventory for nonconvex MINLP. integer_local_search descends the integer lattice directly (1-opt and 2-opt moves on the true-constraint violation), repairing the continuous variables with subnlp at each local minimum — essential when the relaxation’s integers satisfy the relaxed but not the true constraints. enumerate_binary_seeds_subnlp enumerates all 0/1 assignments over a capped binary set and solves a fixed-integer sub-NLP from each, removing a platform-sensitive nondeterminism in which disjunct a nonconvex relaxation locks onto. Both feed only subnlp-verified points to the tree, so — like every heuristic here — they can only ever strengthen the incumbent, never the certificate.