Bound tightening: FBBT, row activity, and OBBT

Bound tightening: FBBT, row activity, and OBBT#

Tightening variable bounds is one of the highest-leverage steps in global optimization: a smaller box yields tighter convex relaxations, fewer branch-and-bound nodes, and faster convergence [Puranik and Sahinidis, 2017]. discopt runs a layered bound-tightening stack:

FBBT (feasibility-based bound tightening) — propagate each constraint’s row activity under the current bounds and isolate each variable, iterated across constraints to a fixpoint. Cheap (interval arithmetic, no LP) and exact on the constraint graph.
Implied / row-activity propagation — the linear specialization of the above.
DBBT (duality-based) — one objective LP’s reduced costs tighten every variable at once when an incumbent cutoff is known.
OBBT (optimization-based) — solve min/max of each variable over the McCormick [McCormick, 1976] relaxation of the model, iterated as the box shrinks and the envelopes tighten [Tawarmalani and Sahinidis, 2005].

Every layer only ever derives valid bounds, so the tightened box always contains the whole feasible region — bound tightening never removes the optimum.

This notebook focuses on the sound, LP-free core — FBBT — exposed publicly as discopt.tightening.fbbt_box, and then shows the full stack at work inside a global solve.

FBBT row-activity propagation#

Consider \(a \ge b + 5\) and \(b \ge c + 3\) with \(a, b, c \in [0, 10]\). No single constraint pins \(a\), but propagating them together tightens the lower bounds: from the second, \(b \ge 3\); feeding that into the first, \(a \ge 8\). FBBT reaches this fixpoint automatically (and tightens the upper bounds the other way).

import os

import numpy as np
import discopt.modeling.core as dm
from discopt.tightening import fbbt_box

m = dm.Model("chain")
a = m.continuous("a", lb=0, ub=10)
b = m.continuous("b", lb=0, ub=10)
c = m.continuous("c", lb=0, ub=10)
m.minimize(a + b + c)
m.subject_to(a - b >= 5)
m.subject_to(b - c >= 3)

res = fbbt_box(m)
print("tightened lower bounds:", np.round(res.lb, 3))   # -> [8, 3, 0]
print("tightened upper bounds:", np.round(res.ub, 3))   # -> [10, 5, 2]
print("variables tightened   :", res.n_tightened)

tightened lower bounds: [8. 3. 0.]
tightened upper bounds: [10.  5.  2.]
variables tightened   : 3

Proving infeasibility#

If propagation collapses a variable’s interval to empty, the problem is infeasible — a result that downstream tools (e.g. conflict analysis) use to emit no-good cuts.

m2 = dm.Model("infeasible")
x = m2.continuous("x", lb=0, ub=10)
m2.minimize(x)
m2.subject_to(x >= 5)
m2.subject_to(x <= 3)

print("FBBT proves infeasible:", fbbt_box(m2).infeasible)

FBBT proves infeasible: True

The full stack inside a global solve#

On a nonconvex MINLP the solver applies the whole stack — FBBT at presolve, then DBBT and OBBT over the McCormick relaxation at the root — before and during spatial branch-and-bound. The bounds only ever shrink, so the certified global optimum is unchanged; the tightening just gets there with less search [Puranik and Sahinidis, 2017].

m3 = dm.Model("nonconvex")
xi = m3.integer("xi", lb=0, ub=4)
y = m3.continuous("y", lb=0, ub=4)
m3.minimize((xi - 1.7) ** 2 + (y - 2.3) ** 2)
m3.subject_to(xi * y <= 3.0)   # nonconvex bilinear

result = m3.solve()
print(f"status: {result.status}, objective: {float(result.objective):.4f}")

status: optimal, objective: 0.4900

Integrality-aware bound snapping#

FBBT is exact on the continuous relaxation of each constraint, but a tightened bound on an integer or binary variable can be sharpened further: a derived bound of n <= 2.5 means n <= 2 once integrality is enforced, and that floored bound can in turn propagate to other variables. discopt’s FBBT now snaps the tightened bounds of integer/binary variables to integer values as it propagates (rounding a lower bound up with ceil, an upper bound down with floor), so each sweep feeds sharper integer-aware intervals back into the fixpoint — the same binary-indicator-style inference that drives strong MIP presolve [Achterberg et al., 2020, Belotti et al., 2009].

Here a single linear constraint forces a fractional upper bound on an integer variable, which FBBT snaps to the next integer.

m4 = dm.Model("int_snap")
n = m4.integer("n", lb=0, ub=10)
z = m4.continuous("z", lb=0, ub=10)
m4.minimize(n + z)
m4.subject_to(2 * n <= 5)   # continuous FBBT gives n <= 2.5 -> snapped to n <= 2

res_int = fbbt_box(m4)
print("integer n upper bound:", float(res_int.ub[0]))   # -> 2.0, not 2.5
print("variables tightened  :", res_int.n_tightened)

integer n upper bound: 2.0
variables tightened  : 1

Lifted-LP FBBT and relaxation equilibration#

The FBBT above propagates the original model’s rows. On a nonconvex MINLP the spatial branch-and-bound also builds a lifted McCormick relaxation [McCormick, 1976, Tawarmalani and Sahinidis, 2005]: every bilinear product w = x_i·x_j (and every monomial) gets its own auxiliary column z, and the McCormick envelopes become linear rows A_ub·z <= b_ub over those lifted columns. That relaxation carries structure the original constraint graph hides.

Per-node lifted-LP FBBT#

Purely-linear FBBT on the original model can miss bilinear-implied factor bounds. A classic case: a product column pinned by the rest of the model (say w = c·x_i·x_j forced to w >= 1) implies a lower bound on a factor x_i — but only through the bilinear relation, which linear FBBT never sees. The opt-in lifted-LP FBBT pass (discopt._relax.mccormick_lp._lifted_lp_fbbt, driving MccormickLPRelaxer._lifted_fbbt_rebuild) propagates the relaxation’s own rows — including McCormick facets such as w <= ub_j·x_i — to recover exactly those factor bounds, then rebuilds the relaxation on the tightened box so its envelopes no longer underestimate the product to zero in the box interior. This is what lets the global dual bound climb off a structural 0 on instances like ex1252 / ex1252a [Puranik and Sahinidis, 2017].

Every row used is a valid relaxation inequality, so each derived bound is implied by valid constraints plus the incoming box — the pass only ever tightens and never excludes a feasible point. It is off by default (it adds an FBBT sweep plus a conditional relaxation rebuild per node); enable it with the environment variable DISCOPT_LIFTED_FBBT=1.

Lifted-relaxation LP equilibration#

The lifted McCormick rows of a product over a wide variable box mix tiny constants (~1e-9) with large bound-derived coefficients (~1e7), so the coefficient spread can exceed 1e15 on boundary sub-boxes. An external LP backend (HiGHS) stalls on such an ill-conditioned matrix. Before the external solve, when the coefficient spread exceeds 1e6, discopt applies geometric-mean (Ruiz) row/column scaling [Ruiz, 2001], with the scale factors snapped to powers of two so the rescaling is exact in floating point (equilibrate_relaxation_lp in discopt._relax.milp_relaxation). Because it is an exact diagonal rescaling x = D·x', the LP optimum — and hence the dual bound — is unchanged; only the returned solution point is mapped back through the column scale. It conditions the solve, it never alters the result. (The pure-Rust simplex equilibrates internally, so the Python pre-scaling is applied only on the external-backend path.)

Below we solve a small bilinear model both ways. On a tiny, well-scaled box the lifted FBBT pass and the equilibration trigger have no visible effect on the already-correct objective — both runs return the same certified optimum, which is the whole point: these passes accelerate and condition the search without ever changing the answer. Their payoff appears on harder instances with wide boxes and deep bilinear coupling, as described above.

def build_bilinear():
    mb = dm.Model("bilinear")
    x = mb.continuous("x", lb=1.0, ub=5.0)
    y = mb.continuous("y", lb=1.0, ub=5.0)
    mb.minimize(x + y)
    mb.subject_to(x * y >= 4.0)   # nonconvex; optimum at x = y = 2, objective 4
    return mb

# Baseline solve (default: lifted-LP FBBT off).
os.environ.pop("DISCOPT_LIFTED_FBBT", None)
base = build_bilinear().solve()
print(f"default     : status={base.status}, objective={float(base.objective):.4f}")

# Same model with per-node lifted-LP FBBT enabled. The certified optimum is
# identical (the pass only tightens valid bounds); on this tiny box it conditions
# the search rather than changing the answer.
os.environ["DISCOPT_LIFTED_FBBT"] = "1"
lifted = build_bilinear().solve()
print(f"lifted FBBT : status={lifted.status}, objective={float(lifted.objective):.4f}")
_ = os.environ.pop("DISCOPT_LIFTED_FBBT", None)

default     : status=optimal, objective=4.0000
lifted FBBT : status=optimal, objective=4.0000