© 2018–2026 by Alfred Galichon. Past and present support from NSF grant DMS-1716489, ERC grant CoG-866274 are acknowledged, as well as inputs from contributors listed here.
If you reuse material from this masterclass, please cite as:
Alfred Galichon, 'math+econ+code' masterclass series. https://www.math-econ-code.org/
import numpy as np
import gurobipy as grb
When a linear program is too large to be solved efficiently, we can sometimes exploit a particular structure. The Dantzig–Wolfe decomposition does this for LPs with a suitable block-diagonal form: it splits the problem into several independent subproblems and a coordinating master problem. The solving algorithm then alternates between them.
The key benefit is scalability: the method starts with a small set of variables and adds more only when they improve the solution, instead of optimizing over all variables from the start.
We will consider the following problem as motivation. A profit-maximizing planner allocates production across $K$ factories indexed by $k\in[K]$. Each factory can produce multiple products (say, different car models), and the set of products may be different across factories. Each factory operates under both local constraints and global constraints on its inputs. Local inputs are for instance work hours or machine time. Global inputs are for instance the different materials required to build the cars, that can be allocated to the different factories by the planner.
Decisions.
Let $G_k$ be the number of car models that can be produced by factory $k$.
The planner chooses for each factory $k$ the quantity $x_g^k$ of model-$g$ cars to produce for all $g \in [G_k]$,
\begin{equation}
x^k = (x_g^k)_{g \in [G_k]} \geq 0.
\end{equation}
Profit.
Each car model $g \in [G_k]$ produced by factory $k$ is sold for profit $\psi_g^k \in \mathbb R$.
The planner wants to maximize the total profit,
\begin{equation}
\sum_{k \in [K]} \sum_{g \in [G_k]} \psi_g^k \, x_g^k = \sum_{k\in[K]} (\psi^k)^\top x^k.
\end{equation}
Local constraints on inputs.
In each factory $k$, production uses a number of local inputs $I_k$ indexed by $i \in [I_k]$ (labor-hours, components, machine time, etc).
Producing one model-$g$ car in factory $k$ requires $a_{ig}^k$ units of input $i$ (different factories may have different production technologies).
The quantity of local input $i$ available to factory $k$ is $b_i^k$.
Thus, for each input $i \in [I_k]$ and factory $k$,
\begin{equation}
\sum_{g \in [G_k]} a_{ig}^k \, x_g^k \leq b_i^k,
\qquad \text{or in matrix form,} \qquad
A^k \, x^k \leq b^k
\end{equation}
where $A^k = (a_{ig}^k)$ is an $I_k \times G_k$ matrix.
Global constraints on inputs. The planner must allocate global resources across the different factories (e.g., raw materials needed to build cars). There is a number $J_0$ of these global inputs, indexed by $j \in [J_0]$. Producing one model-$g$ car in factory $k$ consumes $d_{jg}^k$ units of input $j$ (again, factories may have different production technologies). The total amount of input $j$ available is $b_j^0$. Thus for each global input $j \in [J_0]$, \begin{equation} \sum_{k\in[K]} \sum_{g \in [G_k]} d_{jg}^k \, x_g^k \leq b_j^0, \qquad \text{or in matrix form,} \qquad \sum_{k\in[K]} D^k \, x^k \leq b^0 \end{equation} where $D^k = (d_{jg}^k)$ is a $J_0 \times G_k$ matrix.
It is clear that the above problem is a linear program. However, because some constraints are local while others are global, it admits a special structure.
To see this, stack all decisions vectors $x^k$ into the same vector, and similarly for profits $\psi^k$: \begin{equation} x = \begin{pmatrix} x^1 \\ \vdots \\ x^K \end{pmatrix}, \qquad \psi = \begin{pmatrix} \psi^1 \\ \vdots \\ \psi^K \end{pmatrix}. \end{equation} Then the central planner's problem is the LP \begin{align} \max_{x \geq 0} ~ & \psi^\top x \\ \text{s.t.} ~ & A x \leq b \end{align} where \begin{equation} A = \begin{pmatrix} A^1 & 0 & \cdots & 0 \\ 0 & A^2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & A^K \\ D^1 & D^2 & \cdots & D^K \end{pmatrix} \qquad \text{and} \qquad b = \begin{pmatrix} b^1 \\ b^2 \\ \vdots\\ b^K \\ b^0 \end{pmatrix}. \end{equation}
This is the canonical structure for Dantzig–Wolfe: many local blocks plus a small set of global linking constraints.
We encode this block structure in a DantzigWolfe Python class as follows:
A_k: list of subproblem constraint matrices $\left[A^1, \ldots, A^K\right]$, where $A^k x^k \leq b^k \quad \forall k \in [K]$b_k: list of budget constraint vectors $\left[b^1, \ldots, b^K, b^0\right]$, the first $K$ entries are the subproblem constraints, and the last entry $b^0$ is the linking constraintD_k: list of linking constraint matrices $\left[D^1, \ldots, D^K\right]$, where $\sum_k D^k x^k \leq b^0$psi_k: list of objective vectors $\left[\psi^1, \ldots, \psi^K\right]$, where the objective is $\sum_k (\psi^k)^\top x^k$class DantzigWolfe():
def __init__(self, A_k, D_k, b_k, psi_k):
self.K = len(A_k)
self.I_k = [ A_k[k].shape[0] for k in range(self.K) ]
self.G_k = [ A_k[k].shape[1] for k in range(self.K) ]
if len(D_k) != self.K or len(b_k) != self.K+1 or len(psi_k) != self.K:
raise ValueError(f"Dimensions of inputs do not match.")
self.A_k, self.D_k, self.b_k, self.psi_k = A_k, D_k, b_k, psi_k
A_block = np.block([ [ A_k[k] if k == l else np.zeros((A_k[k].shape[0], A_k[l].shape[1])) for l in range(self.K) ]
for k in range(self.K) ])
self.A = np.vstack((A_block, np.hstack(D_k)))
self.b, self.psi = np.hstack(b_k), np.hstack(psi_k)
We generate a random example.
def generate_example(seed=0):
np.random.seed(seed)
K = 5
I_k = [np.random.randint(1, 5) for _ in range(K)] # numbers of local inputs per factory
J0 = np.random.randint(1, 5) # number of global inputs
G_k = [np.random.randint(1, 11) for _ in range(K)] # numbers of car models per factory
A_k = [np.random.randint(1, 6, size=(I_k[k], G_k[k])) for k in range(K)]
D_k = [np.random.randint(1, 6, size=(J0, G_k[k])) for k in range(K)]
b_k = [np.random.randint(50, 1000, size=I_k[k]) for k in range(K)] \
+ [np.random.randint(50, 1000, size=J0)]
psi_k = [np.random.randint(1, 21, size=G_k[k]) for k in range(K)]
return DantzigWolfe(A_k, D_k, b_k, psi_k)
ex = generate_example()
def DantzigWolfe_gurobi_solve(self, verbose=0):
m = grb.Model()
x = m.addMVar(np.sum(self.G_k), lb=0.0)
m.setObjective(self.psi.T @ x, grb.GRB.MAXIMIZE)
m.addConstr(self.A @ x <= self.b)
m.optimize()
return m.ObjVal, np.array(m.X)
DantzigWolfe.gurobi_solve = DantzigWolfe_gurobi_solve
ex.gurobi_solve()
Gurobi Optimizer version 12.0.0 build v12.0.0rc1 (mac64[arm] - Darwin 24.6.0 24G624)
CPU model: Apple M1
Thread count: 8 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 16 rows, 32 columns and 212 nonzeros
Model fingerprint: 0x76468854
Coefficient statistics:
Matrix range [1e+00, 5e+00]
Objective range [2e+00, 2e+01]
Bounds range [0e+00, 0e+00]
RHS range [8e+01, 1e+03]
Presolve removed 5 rows and 0 columns
Presolve time: 0.00s
Presolved: 11 rows, 32 columns, 180 nonzeros
Iteration Objective Primal Inf. Dual Inf. Time
0 1.0124271e+04 5.369962e+02 0.000000e+00 0s
7 1.9866667e+03 0.000000e+00 0.000000e+00 0s
Solved in 7 iterations and 0.00 seconds (0.00 work units)
Optimal objective 1.986666667e+03
(1986.6666666666667,
array([ 0. , 0. , 0. , 0. , 0. ,
0. , 0. , 4. , 0. , 0. ,
0. , 67. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 28.66666667,
0. , 0. , 0. , 0. , 0. ,
0. , 15.16666667, 0. , 0. , 0. ,
0. , 0. ]))
The Dantzig–Wolfe decomposition is a method to solve LPs of the above form at very large scales. To do this, it reduces the number of variables of the problem and iteratively adds new ones if this is needed.
For each $k$, consider the set of non-zero extreme points $\{x^{k,p}\}_{p\in[P_k]}$ (where $P_k$ is the number of such points) of the polytope defined by the subproblem: \begin{equation} A^k x^k \leq b^k \quad \text{with} \quad x^k \geq 0. \end{equation}
Our previous linear program can then be recast as one where we choose how much weight $\lambda^{k,p}$ to put on each extreme point $p$ of subproblem $k$. Namely:
\begin{align} \max_{\lambda^{k,p} \geq 0} ~ & \sum_{k\in[K]} \sum_{p\in[P_k]} \lambda^{k,p} \big( (\psi^k)^\top x^{k,p} \big) \tag{M} \\ \text{s.t.} ~ & \sum_{k\in[K]} \sum_{p\in[P_k]} \lambda^{k,p} (D^k \, x^{k,p}) \leq b^0 \\ & \sum_{p\in[P_k]} \lambda^{k,p} \leq 1 \quad \forall k \in [K]. \end{align}We call this the full master problem.
Initially, we know that the vector $0_{G_k}$ is an extreme point for subproblem $k$, but we do not know what the other extreme points are, i.e. we do not have an explicit list of the points $x^{k,p}$. However, we shall discover these points as needed.
To this end, we code a function which initializes the master problem while keeping all $\lambda^{k,p}$ equal to 0 (that is, without any nontrivial extreme points).
def DantzigWolfe_initialize_master(self):
m = grb.Model()
m.Params.OutputFlag = 0
m.setObjective(0, grb.GRB.MAXIMIZE) # objective is zero because initially all lambdas are zero
self.linking_j = [ m.addLConstr(0.0, '<', self.b_k[-1][j]) for j in range(len(self.b_k[-1])) ] # initialize 0 <= b^0
self.convexity_k = [ m.addLConstr(0.0, '<', 1.0) for k in range(self.K) ] # initialize 0 <= 1
self.lambda_vars = {}
self.master = m
DantzigWolfe.initialize_master = DantzigWolfe_initialize_master
ex.initialize_master()
To find a first nontrivial extreme point for each subproblem $k$, we can simply solve that subproblem: \begin{align} \max_{x^k \geq 0} ~ & (\psi^k)^\top x^k \\ \text{s.t.} ~ & A^k x^k \leq b^k. \end{align}
We find these first points with initalize_extreme_points and add the corresponding variable symbols $\lambda^{k,p}$, self.lambda_vars[k,self.p_k[k]], to our master problem with update_master.
Note that this is still simply setting up the problem, not numerically solving for $\lambda^{k,p}$ yet.
This corresponds to activating a new column (objective coefficient and constraint matrix) in our problem (i.e. using a new variable), hence the name column generation.
def DantzigWolfe_initalize_extreme_points(self):
self.x_k_p = [ [] for _ in range(self.K) ]
self.p_k = [ 0 for _ in range(self.K) ] # indices for last discovered extreme points
new_x_k = [ False for _ in range(self.K) ]
for k in range(self.K):
m = grb.Model()
x = m.addMVar(self.G_k[k], lb=0.0)
m.setObjective(self.psi_k[k].T @ x, grb.GRB.MAXIMIZE)
m.addConstr(self.A_k[k] @ x <= self.b_k[k])
m.Params.OutputFlag = 0
m.optimize()
new_x_k[k] = m.X
self.x_k_p[k].append(m.X)
return new_x_k
DantzigWolfe.initalize_extreme_points = DantzigWolfe_initalize_extreme_points
def DantzigWolfe_update_master(self, new_x_k):
for k in range(self.K):
if new_x_k[k]:
col = grb.Column()
col.addTerms((self.D_k[k] @ new_x_k[k]).tolist(), self.linking_j)
col.addTerms(1.0, self.convexity_k[k])
self.lambda_vars[k,self.p_k[k]] = self.master.addVar( obj=float(self.psi_k[k] @ new_x_k[k]), lb=0.0, column=col)
DantzigWolfe.update_master = DantzigWolfe_update_master
new_x_k = ex.initalize_extreme_points()
ex.update_master(new_x_k)
Now that our problem is updated with new extreme points, we can consider the restricted version of the master problem which only uses these points. We call this the restricted master problem, and we solve it.
ex.master.Params.OutputFlag = 1
ex.master.optimize()
ex.master.ObjVal
Set parameter OutputFlag to value 1
Gurobi Optimizer version 12.0.0 build v12.0.0rc1 (mac64[arm] - Darwin 24.6.0 24G624)
CPU model: Apple M1
Thread count: 8 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 9 rows, 5 columns and 25 nonzeros
Model fingerprint: 0xa671fe77
Coefficient statistics:
Matrix range [1e+00, 4e+03]
Objective range [1e+03, 2e+04]
Bounds range [0e+00, 0e+00]
RHS range [1e+00, 1e+03]
Presolve removed 6 rows and 0 columns
Presolve time: 0.00s
Presolved: 3 rows, 5 columns, 15 nonzeros
Iteration Objective Primal Inf. Dual Inf. Time
0 5.4720000e+03 4.093726e+01 0.000000e+00 0s
3 1.8720000e+03 0.000000e+00 0.000000e+00 0s
Solved in 3 iterations and 0.00 seconds (0.00 work units)
Optimal objective 1.872000000e+03
1872.0
The value of the objective is lower than the optimal one found above. This is normal: we haven't considered the full set of extreme points yet.
The main interest of considering a restricted master problem is that its size is smaller than that of the full master problem; in particular when the numbers $G_k$ are large, then the restricted problem is typically much smaller than the full one, hence much easier to solve.
Our goal is to construct a restricted master problem whose solution coincides with the solution to the full master problem.
The question of interest is therefore:
When is a solution to the restricted master problem also optimal for the full master problem?
Let's first derive the dual of the full master problem:
\begin{align} & \max_{\lambda^{k,p} \geq 0} ~ \bigg\{ \sum_{k\in[K]} \sum_{p\in[P_k]} \lambda^{k,p} \big( (\psi^k)^\top x^{k,p} \big) \quad \text{s.t.} ~ \sum_{k\in[K]} \sum_{p\in[P_k]} \lambda^{k,p} (D^k \, x^{k,p}) \leq b^0, \quad \sum_{p\in[P_k]} \lambda^{k,p} \leq 1 \quad (\forall k) \bigg\} \\ = & \max_{\lambda^{k,p} \geq 0} ~ \min_{\sigma^k, \pi_j \geq 0} ~ \sum_{k\in[K]} \sum_{p\in[P_k]} \lambda^{k,p} \big( (\psi^k)^\top x^{k,p} \big) + \pi^\top \Big( b^0 - \sum_{k\in[K]} \sum_{p\in[P_k]} \lambda^{k,p} (D^k \, x^{k,p}) \Big) + \sum_{k\in[K]} \sigma^k \Big( 1 - \sum_{p\in[P_k]} \lambda^{k,p} \Big) \\ = & \min_{\sigma^k, \pi_j \geq 0} ~ \max_{\lambda^{k,p} \geq 0} ~ \pi^\top b^0 + \sum_{k\in[K]} \sigma^k + \sum_{k\in[K]} \sum_{p\in[P_k]} \lambda^{k,p} \big( (\psi^k)^\top x^{k,p} - \pi^\top D^k \, x^{k,p} - \sigma^k \big) \end{align}which finally yields the dual problem to the full master problem $(\text{M})$:
\begin{align} \min_{\pi, \sigma \geq 0} ~& \pi^\top b^0 + \sum_{k\in[K]} \sigma^k \\ \text{s.t.} ~ & \pi^\top D^k \, x^{k,p} + \sigma^k \geq (\psi^k)^\top x^{k,p} \quad \forall k \in [K], \forall p\in[P_k]. \end{align}Recall the optimality conditions for a linear program: (1) primal feasibility, (2) dual feasibility, and (3) complementary slackness. If we have a primal-dual solution to a restricted master problem:
We are thus left to check whether dual feasibility holds.
The difficulty is that we cannot check the dual constraints one by one since we do not know the full set of extreme points $\{x^{k,p}\}_{p\in[P_k]}$ for all $k$. The key observation is that we actually do not need to know this set in order to determine whether at least one constraint is violated: we just need to study the following linear program (known as the reduced-cost problem).
Proposition (Column generation). Let $(\pi, \sigma^k)$ be a dual solution to a restricted master problem. For all $k$, the dual constraint $\pi^\top D^k \, x^{k,p} + \sigma^k \geq (\psi^k)^\top x^{k,p}$ is satisfied for all $p\in[P_k]$ if and only if the value of the linear program \begin{align} \max_{x^k \geq 0} ~ & \big( (\psi^k)^\top - \pi^\top D^k \big) x^k - \sigma^k \\ \text{s.t.} ~ & A^k x^k \leq b^k \end{align} is nonpositive.
If for some $k$, the program above takes a strictly positive value, it means that its solution is a new extreme point $x^{k,p}$ that we have to consider.
We thus add it to our restricted master problem together with its associated variable symbol $\lambda^{k,p}$, self.lambda_vars[k,self.p_k[k]].
Again, we do not numerically solve for $\lambda^{k,p}$ yet.
def DantzigWolfe_column_generation(self, tol=1e-6):
pi_j = np.array(self.master.getAttr("Pi", self.linking_j))
sigma_k = np.array(self.master.getAttr("Pi", self.convexity_k))
new_x_k = [ False for _ in range(self.K) ]
for k in range(self.K):
adjusted_psi_k = self.psi_k[k] - self.D_k[k].T @ pi_j
m = grb.Model()
m.Params.OutputFlag = 0
x = m.addMVar(self.G_k[k], lb=0.0)
m.setObjective(adjusted_psi_k @ x, grb.GRB.MAXIMIZE)
m.addConstr(self.A_k[k] @ x <= self.b_k[k])
m.optimize()
if m.Status == grb.GRB.OPTIMAL and m.ObjVal - sigma_k[k] > tol:
self.p_k[k] += 1
self.x_k_p[k].append(m.X)
new_x_k[k] = m.X
return new_x_k
DantzigWolfe.column_generation = DantzigWolfe_column_generation
We now summarize the steps of the algorithm:
Dantzig—Wolfe algorithm.
Step 0. For each $k$, compute an initial extreme point $x^{k,1} \in \arg\max_{x^k \geq 0} \big\{ (\psi^k)^\top x^k \; | \; A^k x^k \leq b^k \big\}$.
Initialize the restricted master problem with the corresponding columns $x^{k,1}$, and set $p_k=1$ as the number of columns generated up to now for subproblem $k$.
Step 1. Solve the restricted master problem \begin{align*} \max_{\lambda^{k,p} \geq 0} ~ & \sum_{k\in[K]} \sum_{p\in[p_k]} \lambda^{k,p} \big( (\psi^k)^\top x^{k,p} \big) \\ \text{s.t.} ~ & \sum_{k\in[K]} \sum_{p\in[p_k]} \lambda^{k,p} (D^k x^{k,p}) \leq b^0 \\ & \sum_{p\in[p_k]} \lambda^{k,p} \leq 1 \qquad \forall k \in [K], \end{align*} and let $(\pi,\sigma^k)$ denote an optimal dual solution.
Step 2. For each $k$, solve the column generation problem \begin{align*} \max_{x^k \geq 0} ~ & \big( (\psi^k)^\top - \pi^\top D^k \big)x^k - \sigma^k \\ \text{s.t.} ~ & A^k x^k \leq b^k. \end{align*} If all problems have value $\leq 0$, stop. Otherwise, add every improving solution $x^k$ as a new extreme point, increment the corresponding $p_k$, and return to Step 1.
The dual variables $\pi_j$ are the shadow prices of the common resources. Given these prices, in Step 2 factory $k$ chooses the feasible plan with the highest profit net of the social cost of the common resources it uses.
If this value is positive, factory $k$ has found a plan that raises total profit, so the corresponding column should be added to the master problem.
If no factory can do so, the current allocation is optimal.
Dantzig–Wolfe is thus decentralization at work: the master problem sets scarcity prices, and factories optimize locally at those prices.
We finally code the full algorithm as dantzig_wolfe_solve.
def DantzigWolfe_dantzig_wolfe_solve(self, max_iter=100, tol=1e-6, verbose=1):
self.initialize_master()
new_x_k = self.initalize_extreme_points()
it = 0
while any(new_x_k) and it <= max_iter:
it += 1
self.update_master(new_x_k)
self.master.optimize()
new_x_k = self.column_generation()
if verbose >= 1:
print(f"Iter {it}: obj = {self.master.ObjVal:.6f} (+{self.K - new_x_k.count(False)} columns)")
if verbose >= 1 and it <= max_iter:
print(f"Optimal after {it} iterations. Obj = {self.master.ObjVal:.6f}")
x = np.hstack([
sum(self.lambda_vars[k, p].X * np.array(self.x_k_p[k][p]) for p in range(self.p_k[k]+1))
for k in range(self.K)
])
return self.master.ObjVal, x
DantzigWolfe.dantzig_wolfe_solve = DantzigWolfe_dantzig_wolfe_solve
ex.dantzig_wolfe_solve()
Iter 1: obj = 1872.000000 (+1 columns) Iter 2: obj = 1962.400000 (+1 columns) Iter 3: obj = 1986.666667 (+0 columns) Optimal after 3 iterations. Obj = 1986.666667
(1986.6666666666667,
array([ 0. , 0. , 0. , 0. , 0. ,
0. , 0. , 4. , 0. , 0. ,
0. , 67. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 28.66666667,
0. , 0. , 0. , 0. , 0. ,
0. , 15.16666667, 0. , 0. , 0. ,
0. , 0. ]))