Dantzig Wolfe Decomposition

Dantzig-Wolfe Decomposition

You are an expert in Dantzig-Wolfe (DW) decomposition for linear and mixed-integer programs. This skill covers detecting block-angular structure in a constraint matrix, reformulating the original (compact) model into a master problem with convexity constraints plus independent pricing subproblems, solving the reformulation by column generation, and interpreting the resulting bound against the LP relaxation, the Lagrangian dual, and the integer optimum. Use the framework below to decide whether a model is worth decomposing, to execute the reformulation correctly, and to verify the bound relationships on the user's instance.

Initial Assessment

Establish the following before reformulating anything:

• Structure. Which constraint rows couple otherwise-independent variable groups? Ask the user to name the natural blocks (plants, vehicles, machines, periods, scenarios). If they cannot, run structure detection on the constraint matrix (section below).
• Linking fraction. Count linking rows m0 versus total rows. DW pays off when m0 is small relative to the block rows — a useful rule of thumb is linking rows below 10-20% of all rows.
• Problem class. LP or MIP? If MIP, locate the integrality: integer variables inside blocks make the DW bound potentially stronger than the LP bound; integrality that lives only in the linking rows gains nothing from convexification.
• Block inventory. Number of blocks K, variables per block, and whether the blocks are identical (same costs, same constraint data). Identical blocks call for the aggregated master, which removes symmetry.
• Pricing tractability. What does one block look like in isolation? A knapsack, a shortest path, a small assignment, a small MIP? The whole method stands or falls on solving the pricing problem quickly and repeatedly.
• Boundedness. Are the block polyhedra bounded? If not, the implementation must handle extreme rays, not just extreme points.
• Goal. A tighter dual bound, a faster LP solve on a huge structured model, or an integer optimum? The first two end with column generation; the third requires branch-and-price (hand off to the column-generation skill).
• Baseline. Solve the compact model (or its LP relaxation) first. Record z_LP, the MIP gap, and the time. Without this baseline you cannot say whether DW helped.
• Solver access. Gurobi license for master and pricing? If not, plan for HiGHS as the master LP solver and GCG/SCIP for an automatic end-to-end alternative.
• Stopping policy. Run column generation to proven optimality, or stop early on a Lagrangian-bound gap? Agree on the tolerance up front.
• Time budget. Per-iteration cost is one RMP LP plus K pricing solves. Estimate iterations in the tens-to-hundreds range and check the budget supports that.

Decomposition Anatomy

Block-angular form

DW applies to problems whose constraint matrix is block-angular: a few linking rows across all variables, then independent diagonal blocks.

$$ z{\mathrm{IP}} \;=\; \min \;\sum{k=1}^{K} ck^{\top} xk \quad \text{s.t.} \quad \sum{k=1}^{K} Ak xk \;\{\le,=,\ge\}\; b \;\;[\pi], \qquad xk \in X_k, \;\; k = 1,\dots,K, $$

where $Xk = \{x \in \mathbb{Z}+^{nk} \text{ (or } \mathbb{R}+^{nk}\text{)} : Dk x \le dk\}$ collects the block-k constraints and the linking constraints have $m0$ rows with duals $\pi$.

Reformulation by convexification

By the Minkowski-Weyl theorem, every point of $\operatorname{conv}(Xk)$ is a convex combination of its extreme points $\{xk^p\}{p \in Pk}$ plus a conic combination of its extreme rays $\{rk^q\}{q \in Q_k}$. Substituting

$$ xk = \sum{p \in Pk} \lambda{kp}\, xk^p + \sum{q \in Qk} \mu{kq}\, rk^q, \qquad \sum{p \in Pk} \lambda{kp} = 1, \quad \lambda, \mu \ge 0 $$

into the linking constraints yields the DW master problem:

$$ z{\mathrm{DW}} = \min \sum{k,p} (ck^{\top} xk^p)\, \lambda{kp} + \sum{k,q} (ck^{\top} rk^q)\, \mu_{kq} $$

$$ \text{s.t.} \quad \sum{k,p} (Ak xk^p)\, \lambda{kp} + \sum{k,q} (Ak rk^q)\, \mu{kq} \;\{\le,=,\ge\}\; b \;\;[\pi], \qquad \sum{p} \lambda{kp} = 1 \;\;[\sigma_k], \;\; k=1,\dots,K. $$

The per-block equations $\sump \lambda{kp} = 1$ are the convexity constraints; their duals $\sigma_k$ price block membership. The master has one column per extreme point or ray — exponentially many — so it is solved by column generation over a restricted master problem (RMP): solve the RMP, read $(\pi, \sigma)$, and for each block solve the pricing problem

$$ \bar{c}k \;=\; \min{x \in Xk} \;(ck - Ak^{\top}\pi)^{\top} x \;-\; \sigmak . $$

If $\bar{c}_k list[tuple[np.ndarray, np.ndarray]]: """Split A into independent blocks after removing the linking rows.

Returns one (rowindices, colindices) pair per block, with row indices referring to the original matrix. Columns appearing only in linking rows form their own zero-row blocks and are reported too. """ m, n = A.shape keep = np.setdiff1d(np.arange(m), np.asarray(linkingrows, dtype=int)) pattern = (A[keep, :] != 0).astype(np.int8) # Bipartite graph: kept rows are nodes 0..len(keep)-1, columns follow. graph = sp.bmat([[None, pattern], [pattern.T, None]], format="csr") , labels = connected_components(graph, directed=False) blocks = [] for comp in np.unique(labels[len(keep):]): # components holding a column rows = keep[labels[: len(keep)] == comp] cols = np.where(labels[len(keep):] == comp)[0] blocks.append((rows, cols)) return blocks

def greedylinkingrows(A: sp.csrmatrix, maxlinking: int) -> np.ndarray: """Cheap detection heuristic: peel off the densest rows one at a time.

Dense rows are the most likely coupling rows. Keep the smallest removal set that maximizes the block count within the budget. Serious detection uses hypergraph partitioning (Bergner et al. 2015). """ density = np.diff(A.tocsr().indptr) # nonzeros per row order = np.argsort(-density) best = np.array([], dtype=int) bestcount = len(findblocks(A, best)) for take in range(1, maxlinking + 1): cand = order[:take] count = len(findblocks(A, cand)) if count > bestcount: best, bestcount = cand, count return np.sort(best)

if __name__ == "__main__": # 3 blocks of 2 vars each, 2 block rows per block, 2 dense linking rows on top. block = np.array([[3.0, 5.0], [4.0, 2.0]]) A = sp.csrmatrix(np.vstack([np.ones((2, 6)), sp.blockdiag([block] * 3).toarray()])) linking = greedylinkingrows(A, maxlinking=3) print("linking rows:", linking) for rows, cols in findblocks(A, linking): print("block rows", rows, "cols", cols) # Expected: linking rows [0 1]; three blocks with rows [2 3], [4 5], [6 7] # and columns [0 1], [2 3], [4 5] respectively.

If the heuristic finds no decomposition, the model may still decompose after assigning a small number of "ambiguous" rows to blocks by hand, or after a variable permutation suggested by domain knowledge (one block per machine, period, or scenario). When the matrix is born from a model you wrote, prefer annotating blocks at modeling time over rediscovering them numerically.

## Generic Dantzig-Wolfe Implementation

The driver below is problem-independent. Each block supplies its cost vector, its linking-row coefficients, and a pricing oracle; the driver owns the RMP, the artificial columns, the dual extraction, the ray handling, and the Lagrangian bound.

```text
DANTZIG-WOLFE COLUMN GENERATION (minimization)

Input: blocks k = 1..K with costs c_k, linking coefficients A_k, sets X_k;
       linking RHS b with senses; tolerance tol.

1   build RMP: linking rows, one convexity row per block,
    big-M artificial columns so the RMP starts feasible
2   repeat
3       solve RMP (LP)  ->  z_RMP, duals pi (linking), sigma_k (convexity)
4       for each block k:
5           solve pricing  v_k = min { (c_k - A_k' pi)' x : x in X_k }
6           if unbounded with extreme ray r:  rc = (c_k - A_k' pi)' r
7           else:                             rc = v_k - sigma_k
8           if rc  Callable[[np.ndarray], PricingResult]:
    """Build a pricing oracle for X_k = {x >= 0 : D x  PricingResult:
        """Minimize price @ x over the block polyhedron."""
        model.setObjective(price @ x, GRB.MINIMIZE)
        model.optimize()
        if model.Status == GRB.OPTIMAL:
            return model.ObjVal, x.X.copy(), False
        if model.Status == GRB.UNBOUNDED and vtype == GRB.CONTINUOUS:
            ray = np.array(model.getAttr("UnbdRay", model.getVars()))
            ray /= np.abs(ray).max()           # scale for numerical sanity
            return float(price @ ray), ray, True
        raise RuntimeError(f"pricing status {model.Status}: integer blocks must be bounded")

    return pricing

class DantzigWolfe:
    """Dantzig-Wolfe master with one convexity constraint per block."""

    def __init__(self, blocks: list[Block], b: np.ndarray, senses: str,
                 big_m: float = 1e6) -> None:
        """senses is one character per linking row, each of ''."""
        self.blocks = blocks
        self.b = np.asarray(b, dtype=float)
        self.master = gp.Model("rmp")
        self.master.Params.OutputFlag = 0
        smap = {"": GRB.GREATER_EQUAL}
        self.link = [self.master.addLConstr(gp.LinExpr(), smap[s], rhs, name=f"link[{i}]")
                     for i, (s, rhs) in enumerate(zip(senses, self.b))]
        self.conv = [self.master.addLConstr(gp.LinExpr(), GRB.EQUAL, 1.0, name=f"conv[{k}]")
                     for k in range(len(blocks))]
        # Artificial variables keep the RMP feasible before real columns exist.
        self.artificials: list[gp.Var] = []
        for i, s in enumerate(senses):
            if s in ("=", ">"):
                self.artificials.append(self.master.addVar(
                    obj=big_m, name=f"art_up[{i}]", column=gp.Column([1.0], [self.link[i]])))
            if s in ("=", " None:
        """Append one point or ray column for block k to the RMP."""
        blk = self.blocks[k]
        coeffs = [float(v) for v in blk.A @ x]
        constrs: list[gp.Constr] = list(self.link)
        if not is_ray:
            coeffs.append(1.0)               # convexity coefficient (points only)
            constrs.append(self.conv[k])
        var = self.master.addVar(obj=float(blk.c @ x), lb=0.0,
                                 name=f"col[{k},{len(self.columns)}]",
                                 column=gp.Column(coeffs, constrs))
        self.columns.append((k, x.copy(), is_ray, var))

    def solve(self, tol: float = 1e-6, max_iters: int = 200) -> dict:
        """Run column generation; return bounds, recovered x_k, and the iteration log."""
        best_lb = -math.inf
        z_rmp = math.inf
        for it in range(1, max_iters + 1):
            self.master.optimize()
            if self.master.Status != GRB.OPTIMAL:
                raise RuntimeError(f"RMP status {self.master.Status}")
            z_rmp = self.master.ObjVal
            pi = np.array([c.Pi for c in self.link])
            sigma = np.array([c.Pi for c in self.conv])
            new_cols, rc_sum, ray_seen = 0, 0.0, False
            for k, blk in enumerate(self.blocks):
                value, x, is_ray = blk.pricing(blk.c - blk.A.T @ pi)
                rc = value if is_ray else value - sigma[k]
                ray_seen = ray_seen or is_ray
                rc_sum += 0.0 if is_ray else rc
                if rc  1e-6 for a in self.artificials):
            raise RuntimeError("artificials positive at the end: instance infeasible "
                               "or big_m too small")
        x_blocks = [np.zeros(len(blk.c)) for blk in self.blocks]
        for k, x, _, var in self.columns:
            x_blocks[k] += var.X * x
        return {"z_dw": z_rmp, "lower_bound": best_lb,
                "iterations": len(self.history), "x_blocks": x_blocks}

Design notes. The pricing oracle is a closure over a persistent Gurobi model, so each call only swaps the objective — no model rebuild. Reduced-cost sums feed the Lagrangian bound $L(\pi) = z{\mathrm{RMP}} + \sumk \bar{c}k$ every iteration, which is what allows the early-exit gap test on line `zrmp - bestlb dict: """Three plants, two products, two resources per plant. Integer batches.""" return { "cost": np.array([[4.0, 7.0], [5.0, 6.0], [6.0, 5.0]]), # (K, J) "use": [np.array([[3.0, 5.0], [4.0, 2.0]]), # Rk: (2, J) np.array([[2.0, 4.0], [5.0, 3.0]]), np.array([[4.0, 3.0], [3.0, 4.0]])], "cap": [np.array([25.0, 22.0]), np.array([20.0, 24.0]), np.array([26.0, 23.0])], "demand": np.array([10.0, 8.0]), }

def solve_compact(data: dict, relax: bool) -> float: """Solve the compact model; relax=True drops integrality on the batches.""" K, J = data["cost"].shape model = gp.Model("multiplant") model.Params.OutputFlag = 0 vtype = GRB.CONTINUOUS if relax else GRB.INTEGER x = model.addMVar((K, J), lb=0.0, vtype=vtype, name="x") for k in range(K): model.addConstr(data["use"][k] @ x[k, :] = data["demand"], name="demand") model.setObjective((data["cost"] * x).sum(), GRB.MINIMIZE) model.optimize() assert model.Status == GRB.OPTIMAL return model.ObjVal

if __name__ == "__main__": data = makeinstance() print(f"compact LP relaxation: {solvecompact(data, relax=True):.4f}") print(f"compact IP optimum: {solve_compact(data, relax=False):.4f}") # Expected: compact LP relaxation 90.0909, compact IP optimum 92.0000

Now decompose: one `Block` per plant with integer pricing over its own capacity polytope, identity linking coefficients, and `>` senses on the two demand rows.

```python
"""Dantzig-Wolfe on the multi-plant instance: convexified integer blocks."""
import numpy as np
from gurobipy import GRB

from dw_framework import Block, DantzigWolfe, make_block_pricing   # framework above
from multiplant import make_instance, solve_compact                # compact model above

def run() -> None:
    """Compare z_LP " * J)
    result = dw.solve()
    z_lp = solve_compact(data, relax=True)
    z_ip = solve_compact(data, relax=False)
    print(f"z_LP = {z_lp:.4f}  z_DW = {result['z_dw']:.4f}  z_IP = {z_ip:.4f}")
    print(f"iterations = {result['iterations']}, columns = {len(dw.columns)}")
    for k, xk in enumerate(result["x_blocks"]):
        print(f"plant {k}: x = {np.round(xk, 4)}")
    for row in dw.history:
        print(row)

if __name__ == "__main__":
    run()
    # Expected: z_LP = 90.0909  [5, 1], plant 1 -> [3.2, 2.6], plant 2 -> [1.8, 4.4];
    # plants 1 and 2 are fractional convex combinations of integer pricing points.

Reading the run: the first two RMP values are dominated by big-M artificials; once real columns cover the demand rows, zRMP drops to the 90s while the Lagrangian bound climbs (73.2 → 86.0 → 91.0 → 91.4) until the two meet at the DW optimum. Every column generated is an integer production plan for one plant, yet the master mixes them fractionally — the bound improved from 90.09 to 91.4 exactly because conv(Xk) cuts off fractional plans like x = (8.33, 0) that the compact LP relaxation allows.

Worked Example 2: Cutting Stock — Compact vs. DW Reformulation

The one-dimensional cutting-stock problem is the canonical demonstration that the same problem under two formulations gives very different bounds. The compact (Kantorovich) model has one identical block per stock roll $k$: $yk \in \{0,1\}$ (roll used) and $x{ik} \in \mathbb{Z}+$ (copies of item i cut from roll k), with block constraint $\sumi li x{ik} \le W yk$ and linking constraints $\sumk x{ik} \ge di$. Its LP relaxation is

…

Source & license

This open-source skill is cataloged on AgentStack and links to its original source — we do not rehost the code.

• Author: hajibabaie
• Source: hajibabaie/combinatorial-optimization-skills
• License: MIT

Install and usage instructions live in the source repository linked above.