The Automatic Decomposition Advisor#
Many optimization models have structure a solver can exploit: fixing a few complicating variables disconnects the problem into independent recourse blocks (Benders decomposition [Benders, 1962], and its nonlinear generalization GBD [Geoffrion, 1972]); a few coupling constraints link otherwise-independent blocks (Lagrangian relaxation [Geoffrion, 1974], Dantzig–Wolfe [Dantzig and Wolfe, 1960]).
Historically the modeler had to recognize the structure and rewrite the model
by hand. The Decomposition Advisor automates that: it analyzes a model,
discovers exploitable structure as a graph problem, scores candidate
decompositions, explains its recommendation, and produces a solvable reformulation
— a compiler-style transformation pass on the same footing as presolve. See the
design spec in docs/design/decomposition-advisor.md; for a broader survey of
Benders methods see Rahmaniani et al. [2017].
This notebook walks the advisor’s analysis surface end to end. It is pure analysis — no solve is triggered — so it runs without the compiled solver backend; running the chosen decomposition is covered in the Benders and Lagrangian tutorials.
import discopt.modeling as dm
from discopt.decomposition import RecordStore, analyze_decomposition, record_outcome
from discopt.decomposition.learning import ObservedPerformance, Outcome
A two-stage facility model#
Three candidate facilities: binary build[i] decides whether to open facility
i (first stage), continuous x[i] is its production (recourse). Production is
gated by the build decision, and a budget limits how many we open — the budget
row couples the binaries, so the model is not already block-diagonal. But
once the binaries are fixed, the recourse splits into three independent blocks:
the textbook Benders structure.
m = dm.Model("facility")
# Scalar variables so each facility is a distinct node in the structure graph.
b = [m.binary(f"b{i}") for i in range(3)] # first-stage: open i?
x = [m.continuous(f"x{i}", lb=0, ub=10) for i in range(3)] # recourse: production
for i in range(3):
m.subject_to(x[i] <= 10 * b[i]) # production gated by build
m.subject_to(b[0] + b[1] + b[2] <= 2) # budget couples the binaries
m.minimize((x[0] + x[1] + x[2]) - 3 * (b[0] + b[1] + b[2]))
print(m)
Model: facility
Variables: 6 (3 continuous, 3 integer/binary)
Constraints: 4
Objective: minimize (((x0 + x1) + x2) - (3 * ((b0 + b1) + b2)))
Parameters: 0
Analyze the structure#
model.analyze_decomposition() returns a DecompositionAdvisor. Its
StructureReport summarizes the exploitable structure in one pass — sizes,
integer localization (how many recourse blocks appear once the integers are
fixed), coupling density, and more.
adv = m.analyze_decomposition()
print(adv.structure().summary())
StructureReport
size: 6 vars, 4 cons, 9 incidences
blocks: 1 (block-diagonal=False)
integer: 3 (50%); projection blocks=3
coupling: 1 constraints (density=25.0%)
bridges: 1; articulation vars=3
nonlinear: False; annotated=False
Candidate decompositions#
The advisor discovers candidate decompositions with a cheap-first cascade of
generators. Each candidate reuses the existing DecompositionStructure, so it
can drive the shipping solvers directly. The monolithic no-decomposition
baseline is always included.
for c in adv.candidates():
print("-", c.summary())
- No decomposition [proven_equivalent] — baseline
- Classical Benders [proven_equivalent, 3 blocks] — integer-variable projection
- Lagrangian relaxation [relaxation, 3 blocks] — bridge-constraint scan
Score and rank#
Scoring is method-aware and measured against the no-decomposition baseline (anchored at 0): a decomposition must beat solving the model as written. The score carries an analytic performance estimate (an Amdahl-style speedup capped by load balance) and a confidence. Unsound candidates are vetoed, not down-weighted.
for cand, sv in adv.scores():
print(f"{cand.method.label:22s} {sv.summary()}")
Classical Benders score=+1.52 (confidence 0.75) — ~2.0x speedup, 100% parallel efficiency over 3 blocks, ~11 iters [confidence 0.75]
Lagrangian relaxation score=+1.32 (confidence 0.54) — ~2.0x speedup, 100% parallel efficiency over 3 blocks, ~50 iters [confidence 0.60]
No decomposition score=+0.00 (confidence 1.00)
The recommendation, explained#
Every recommendation is explainable — the reasons are backed by the same metrics the scorer used, with concerns and the runner-up alternatives. The explanation renders as text, markdown, or JSON (the JSON feeds tooling / the LLM layer).
print(adv.explain())
Recommended: Classical Benders
────────────────────────────────────────
Why:
• Fixing 3 integer variable(s) disconnects the model into 3 recourse blocks.
evidence: integer_localization=1.0, components_after_projection=3
• Recourse subproblems are linear (Benders).
Estimated:
• ~2.0x speedup, 100% parallel efficiency over 3 blocks, ~11 iters [confidence 0.75]
Alternatives considered:
• Lagrangian relaxation (score +1.32) — feasible but relaxation; ranked lower
• No decomposition (score +0.00) — baseline (solve as written)
print(adv.explain("json"))
{
"recommendation": "benders",
"headline": "Recommended: Classical Benders",
"confidence": 0.75,
"reasons": [
{
"claim": "Fixing 3 integer variable(s) disconnects the model into 3 recourse blocks.",
"evidence": {
"integer_localization": 1.0,
"components_after_projection": 3
}
},
{
"claim": "Recourse subproblems are linear (Benders).",
"evidence": {}
}
],
"concerns": [],
"score": {
"num_blocks": 3.0,
"block_size_cv": 0.0,
"coupling_density": 0.25,
"parallel_efficiency": 1.0,
"estimated_speedup": 2.0,
"integer_localization": 1.0,
"expected_cut_strength": 0.9,
"est_communication_cost": 0.5,
"warm_start_potential": 1.0
},
"performance": {
"estimated_speedup": 2.0,
"parallel_efficiency": 1.0,
"num_blocks": 3,
"est_iterations": 11
},
"alternatives": [
{
"method": "lagrangian",
"aggregate": 1.325,
"reason": "feasible but relaxation; ranked lower"
},
{
"method": "none",
"aggregate": 0.0,
"reason": "baseline (solve as written)"
}
]
}
Counterfactuals: “why not method X?”#
The advisor also explains why a method was not chosen — including the missing structural precondition. Here the budget row couples the binaries, so the model is not already block-diagonal and an independent-block decomposition does not apply. (Lagrangian relaxation, by contrast, is a viable alternative — it appears as a runner-up above, since the budget row can be dualized.)
print(adv.explain(method="independent_blocks"))
Independent blocks was not applicable: the model is not already block-diagonal (removing no constraints, it stays connected)
Reformulate: build a solvable decomposition#
decompose() turns the recommended candidate into a DecomposedModel — a
master/subproblem partition with a soundness certificate and an
original-space variable mapping. Its solve() (not called here) dispatches to
the shipping Benders / GBD / Lagrangian drivers.
dcmp = adv.decompose()
print(dcmp.summary())
print()
print("variable roles:")
for name in ["b0", "b1", "b2", "x0", "x1", "x2"]:
print(f" {name}: {dcmp.var_map.role(name)}")
DecomposedModel: Classical Benders
Classical Benders: proven_equivalent — classical Benders reproduces the monolithic optimum for linear recourse
master (Classical Benders): 3 complicating vars, 1 coupling constraint(s)
3 subproblem block(s)
subproblem[0]: 1 vars, 0 cons
subproblem[1]: 1 vars, 0 cons
subproblem[2]: 1 vars, 0 cons
variable roles:
b0: master
b1: master
b2: master
x0: block0
x1: block1
x2: block2
Parallel execution plan#
The independent recourse blocks run on a pluggable communication backend
(sequential or thread pool today), ordered biggest-first to avoid stragglers. The
SchedulingGraph reports the load-balance speedup cap; map_subproblems runs a
per-block function and reduces results back into block order, so the answer is
deterministic regardless of backend. Here we map a trivial pure-Python function
(the block size) rather than a real solve.
print(dcmp.schedule().summary())
block_sizes = dcmp.map_subproblems(lambda sp: sp.size, backend="threads")
print("recourse block sizes (block order):", block_sizes)
SchedulingGraph: 3 blocks, total_cost=3, critical_path=1, ideal_speedup=3.00x, master=True
recourse block sizes (block order): [1, 1, 1]
Learn from experience#
Every solve can be recorded — the instance’s features, the chosen method, what
was predicted, and what was observed — so the advisor’s own error is
measurable and future recommendations can improve. Records persist as local JSON
lines; an InstanceBasedPolicy retrieves nearest past instances to bias
selection (deferring safely to the rule-based policy until enough data accrues).
store = RecordStore() # in-memory; pass a path to persist as JSONL
record_outcome(
adv,
store,
observed=ObservedPerformance(wall_clock_s=1.8, iterations=6),
outcome=Outcome.OPTIMAL,
timestamp=0.0,
)
rec = store.all()[0]
print("chosen :", rec.chosen)
print("predicted speedup:", rec.predicted_speedup)
print("features :", rec.features.to_dict())
chosen : benders
predicted speedup: 2.0
features : {'num_vars': 6, 'num_constraints': 4, 'num_incidences': 9, 'integer_fraction': 0.5, 'coupling_density': 0.25, 'num_blocks': 1, 'blocks_after_integer_projection': 3, 'nonlinear': False}
Summary#
The advisor took a model and, without any manual rewriting, discovered its
two-stage structure, scored the candidates against the no-decomposition baseline,
recommended Benders with an evidence-backed rationale, and produced a
solvable DecomposedModel with a soundness certificate and a parallel execution
plan — then recorded the outcome for future learning.
Future directions (design spec §14) include automatic hypergraph partitioning for block-angular models [Karypis and Kumar, 1998], learned decomposition policies, and graph-neural-network structure prediction.