Lagrangian Relaxation

Lagrangian relaxation [Geoffrion, 1974] dualizes the coupling (linking) constraints of a model — the rows that tie otherwise-independent blocks together. Moving them into the objective with multipliers \(\lambda \ge 0\) leaves a relaxed problem that separates into independent blocks and, for every \(\lambda\), yields a lower bound on the optimum:

\[L(\lambda) = \min_z\; c^\top z + \lambda^\top (A_c z - r_c) \quad\text{s.t. block constraints},\; z \in \text{bounds}.\]

Because the dualized term is \(\le 0\) at any feasible point, \(L(\lambda) \le z^\star\) always. The Lagrangian dual \(\max_{\lambda \ge 0} L(\lambda)\) tightens this bound; for problems whose blocks lack the integrality property it can dominate the LP relaxation [Fisher, 2004]. discopt maximizes the dual by a subgradient method or a bundle (cutting-plane) method, and recovers a feasible primal incumbent with a Lagrangian heuristic [Barahona and Anbil, 2000].

import discopt.modeling as dm

Two blocks tied by a conflict

Two independent blocks each must select one item; a coupling constraint forbids selecting item 0 and item 2 together. Marking that constraint with mark_coupling tells the solver which row to dualize.

Without the coupling the cheapest picks are items 0 and 2 (cost \(2 + 2 = 4\)); the conflict forces a swap, and the optimum is \(3 + 2 = 5\) (items 1 and 2).

m = dm.Model("two_blocks")
x = m.binary("x", shape=(4,))
m.minimize(2 * x[0] + 3 * x[1] + 2 * x[2] + 4 * x[3])
m.subject_to(x[0] + x[1] >= 1)  # block A: pick one
m.subject_to(x[2] + x[3] >= 1)  # block B: pick one

conflict = x[0] + x[2] <= 1  # coupling: not both
m.subject_to(conflict)
m.mark_coupling(conflict)

result = m.solve(decomposition="lagrangian")
print("status   :", result.status)
print("objective:", round(result.objective, 4))
print("bound    :", round(result.bound, 4))
print("x =", result.x["x"])

status   : optimal
objective: 5.0
bound    : 5.0
x = [-0.  1.  1. -0.]

The dual bound is a rigorous lower bound on the optimum; here it meets the recovered incumbent, so the gap is certified. The auto-detector also identifies the linking row, which detect_decomposition reports:

import discopt

print(discopt.detect_decomposition(m).summary())

DecompositionStructure (mixed)
  blocks: 1 (separable=False)
  complicating vars: 1
  coupling constraints: 1

Subgradient vs bundle

The dual can be maximized two ways. The subgradient method (default) takes a Polyak step along the constraint residual; the bundle method builds a piecewise-linear model of the concave dual and maximizes it over a trust box, typically converging in fewer iterations on harder duals. Both return the same rigorous bound.

for method in ("subgradient", "bundle"):
    r = m.solve(decomposition="lagrangian", lagrangian_method=method)
    print(f"{method:12s} -> objective {round(r.objective, 4)}, bound {round(r.bound, 4)}")

subgradient  -> objective 5.0, bound 5.0
bundle       -> objective 5.0, bound 5.0

Inside branch-and-bound

The Lagrangian dual is most powerful inside a MILP tree: at a node whose branching has tightened the variable bounds, the Lagrangian bound can dominate the node’s LP relaxation and prune the subtree. Passing lagrangian_bound=True to a MILP solve enables this hook — multipliers are fixed once at the root, and at each node a valid Lagrangian lower bound is combined with the LP bound. It is opt-in (default off), applies to linear minimization models with coupling structure, and no-ops cleanly otherwise.

# A small integer program with coupling structure, solved with the hook on.
mk = dm.Model("knap_blocks")
z = mk.binary("z", shape=(6,))
mk.minimize(sum((i + 1) * z[i] for i in range(6)))
mk.subject_to(z[0] + z[1] + z[2] >= 1)  # block A
mk.subject_to(z[3] + z[4] + z[5] >= 1)  # block B
link = z[0] + z[3] <= 1  # coupling
mk.subject_to(link)
mk.mark_coupling(link)

r = mk.solve(lagrangian_bound=True)
print("status:", r.status, " objective:", round(r.objective, 4))

status: optimal  objective: 6.0

When to use Lagrangian relaxation

Reach for it when a handful of constraints couple otherwise-separable blocks (shared resources, budgets, network flow balance) — dualizing them yields strong bounds and independent, parallelizable subproblems. The v1 implementation targets linear (mixed-integer) models; the reported bound is always a valid lower bound. Dantzig-Wolfe / column generation / branch-and-price are on the roadmap.