Column generation

Deep dive

How it works

Parts:

• Master problem — the full LP min c·x, A x >= b, x >= 0, one column per pattern; never built.
• Restricted master (RMP) — same LP over a small column set P (initial feasible columns); solved with simplex (warm start) or an interior-point method.
• Duals y — one price per master row, read from the RMP optimum.
• Pricing subproblem — finds min over feasible patterns (c_j - y·a_j); a knapsack (cutting stock) or resource-constrained shortest path (routing, crews).

1. Start with columns that make the RMP feasible (one homogeneous pattern per item, or big-cost artificials).
2. Solve the RMP; take optimum z and duals y.
3. Price: maximise y·a over feasible patterns (cutting stock: unbounded knapsack by DP over the roll width). The reduced cost is c_j - y·a_j; for cutting stock c_j = 1.
4. If the best reduced cost is >= -eps, stop: z is the optimum of the full LP.
5. Otherwise append the column (often several) and go to 2.
6. At any iteration z / (best y·a) (Farley bound, unit costs) is a valid lower bound, so stop early when ceil(bound) equals the incumbent.
7. Integer solution: branch-and-price (Variants), or round and repair heuristically.

Resource-constrained shortest path (RCSP) for routing and crew problems: nodes are tasks or customers, arcs carry dual-adjusted costs (can be negative), and each label tracks resources (time, load, duty hours) with dominance, a Dijkstra-like label-setting, see Dijkstra and time-window feasibility. Elementary routes (no repeated customer) are NP-hard to price; the usual relaxation allows cycles or ng-routes.

   RMP (few columns) --duals y--> pricing (knapsack / RCSP)
        ^                               |
        |      column with c - y.a < 0  |
        +-------------------------------+
   no such column  =>  LP optimal  =>  branch-and-price for integers

Dantzig-Wolfe view: write the subproblem's polyhedron by its extreme points and rays; the master picks convex combinations and keeps the coupling rows. Its bound equals the Lagrangian relaxation of those rows, and the duals are the multipliers.

Code

Cutting stock (Gilmore-Gomory). The RMP is solved through its dual, max d·y s.t. pattern·y <= 1, whose slack basis is feasible, by a dense tableau simplex (Bland's rule). Pricing is an unbounded knapsack DP.

use std::time::Instant;

/// Restricted master: min sum x_j  s.t.  sum_j a_ij x_j >= d_i.
/// Solved through its dual  max d.y  s.t. pat_j . y <= 1  (slack basis is feasible),
/// by a dense tableau simplex with Bland's rule. Returns (bound, y = item duals, x = pattern use).
fn solve_rmp(d: &[f64], pats: &[Vec<usize>]) -> (f64, Vec<f64>, Vec<f64>) {
    let (n, m) = (d.len(), pats.len());
    let w = n + m + 1;
    let mut t = vec![vec![0.0; w]; m + 1]; // rows 0..m constraints, row m objective
    for (j, p) in pats.iter().enumerate() {
        for i in 0..n { t[j][i] = p[i] as f64; }
        t[j][n + j] = 1.0;
        t[j][w - 1] = 1.0;
    }
    for i in 0..n { t[m][i] = -d[i]; }
    let mut basis: Vec<usize> = (0..m).map(|j| n + j).collect();
    while let Some(c) = (0..n + m).find(|&c| t[m][c] < -1e-9) {
        let mut r = usize::MAX;
        for k in 0..m {
            if t[k][c] > 1e-9 {
                let better = r == usize::MAX || {
                    let (a, b) = (t[k][w - 1] / t[k][c], t[r][w - 1] / t[r][c]);
                    a < b - 1e-12 || (a < b + 1e-12 && basis[k] < basis[r])
                };
                if better { r = k; }
            }
        }
        assert!(r != usize::MAX, "dual unbounded: some item fits no pattern");
        let pv = t[r][c];
        for v in t[r].iter_mut() { *v /= pv; }
        let pr = t[r].clone();
        for k in 0..=m {
            if k != r {
                let f = t[k][c];
                if f != 0.0 { for (v, p) in t[k].iter_mut().zip(&pr) { *v -= f * p; } }
            }
        }
        basis[r] = c;
    }
    let mut y = vec![0.0; n];
    for k in 0..m { if basis[k] < n { y[basis[k]] = t[k][w - 1]; } }
    let x = (0..m).map(|j| t[m][n + j]).collect();
    (t[m][w - 1], y, x)
}

/// Pricing: unbounded knapsack max sum y_i a_i  s.t.  sum w_i a_i <= cap  (integer DP).
fn price(y: &[f64], w: &[usize], cap: usize) -> (f64, Vec<usize>) {
    let mut best = vec![0.0; cap + 1];
    let mut pick = vec![usize::MAX; cap + 1];
    for c in 1..=cap {
        best[c] = best[c - 1];
        pick[c] = usize::MAX; // MAX = leave a unit unused
        for i in 0..w.len() {
            if w[i] <= c && best[c - w[i]] + y[i] > best[c] + 1e-12 {
                best[c] = best[c - w[i]] + y[i];
                pick[c] = i;
            }
        }
    }
    let (mut pat, mut c) = (vec![0; w.len()], cap);
    while c > 0 {
        match pick[c] { usize::MAX => c -= 1, i => { pat[i] += 1; c -= w[i]; } }
    }
    (best[cap], pat)
}

fn column_generation(w: &[usize], dem: &[usize], cap: usize) -> (f64, Vec<Vec<usize>>, Vec<f64>) {
    let d: Vec<f64> = dem.iter().map(|&v| v as f64).collect();
    let mut pats: Vec<Vec<usize>> = (0..w.len())
        .map(|i| { let mut p = vec![0; w.len()]; p[i] = cap / w[i]; p }).collect();
    loop {
        let (z, y, x) = solve_rmp(&d, &pats);
        let (val, pat) = price(&y, w, cap);
        if val <= 1.0 + 1e-9 { return (z, pats, x); } // reduced cost 1 - val >= 0 for all columns
        pats.push(pat);
    }
}

fn main() {
    let (w, dem, cap) = ([45, 36, 31, 14], [97, 610, 395, 211], 100);
    let t0 = Instant::now();
    let (z, pats, x) = column_generation(&w, &dem, cap);
    let el = t0.elapsed();
    assert_eq!(column_generation(&w, &dem, cap).0.to_bits(), z.to_bits());

    // no improving column remains
    let d: Vec<f64> = dem.iter().map(|&v| v as f64).collect();
    let (_, y, _) = solve_rmp(&d, &pats);
    assert!(price(&y, &w, cap).0 <= 1.0 + 1e-7);

    let rolls: usize = x.iter().map(|v| (v - 1e-9).ceil() as usize).sum();
    assert!(rolls as f64 >= z - 1e-9);
    println!("patterns {}  LP bound {:.4}  rolls {}  time {:?}", pats.len(), z, rolls, el);
}

Run on widths 45/36/31/14, demands 97/610/395/211, roll 100: 6 patterns, LP bound 452.25, 454 rolls after rounding every count up (the bound proves at least 453), microseconds in release (timing: order of). Asserts: bound <= rolls, no improving column left, bit-identical second run.

Complexity and constants

• Per iteration: one LP re-solve plus one pricing. Cutting stock pricing is O(n * W) for n item types and width W; RCSP is pseudo-polynomial in resources, exponential when routes must be elementary.
• The column pool only grows: purge columns with large positive reduced cost now and then.

Efficiency

• Cost: the pricing solves dominate for routing and crew problems (tens of percent to most of the time, instance dependent, order of); for cutting stock the LP re-solve is comparable. Bound: compute, not memory; the pool is small.
• Levers:
  • Several columns per iteration — fewer RMP solves — when pricing yields many distinct columns cheaply.
  • Heuristic pricing first, exact only when it finds nothing — most iterations skip the expensive solver — when a good heuristic exists.
  • Dual stabilization (smoothing, boxes) — fewer oscillating iterations (surveyed in Lübbecke & Desrosiers 2005) — for tailing-off and degenerate masters.
  • Warm-start from the previous basis — a few pivots instead of a full solve — with simplex.
  • Stop by the Lagrangian bound — skips the tail — when the incumbent equals ceil(bound).
• Naive → standard → state of the art: enumerate all patterns · exponential in items; restricted master plus exact pricing · iterations x (LP + pricing); stabilised, heuristic-first pricing inside branch-and-price · fewest iterations, tail controlled.

Implementation notes

• Use a reduced-cost epsilon (such as 1e-9 at unit cost scale), or the loop cycles on a duplicate column.
• If duals jump, smooth: y_used = a * y_prev + (1 - a) * y_new (Wentges-style smoothing), and re-price at the true duals before declaring optimality (the mis-price check).
• Cost drivers are pricing runs and LP solves, not bytes: batch columns, warm-start, price heuristically first.

Pitfalls

• Tailing-off and dual oscillation in degenerate masters: stabilise, or stop by the bound.
• Branching on pattern variables breaks pricing (the banned column is generated again). Branch on arcs or Ryan-Foster pairs so pricing can encode it.

Variants and alternatives

In production

• Paper, film and steel mills: cutting stock with column generation is the standard bound and a common production method.
• Airlines: crew pairing and rostering are built as set-partitioning masters with RCSP pricing.
• Delivery and transport: exact vehicle routing solvers use branch-price-and-cut; nurse and shift rostering at hospitals; telecom and logistics path-based network design and multicommodity flow.
• Generic MIP-decomposition solvers: GCG detects structure and applies Dantzig-Wolfe automatically on top of SCIP.

Implementations

Sources

• Gilmore, Gomory, A Linear Programming Approach to the Cutting-Stock Problem, Operations Research 9(6), 1961, 849-859. paper
• Dantzig, Wolfe, Decomposition Principle for Linear Programs, Operations Research 8(1), 1960, 101-111. paper
• Lübbecke, Desrosiers, Selected Topics in Column Generation, Operations Research 53(6), 2005, 1007-1023. paper