Dynamic Programming

Dynamic Programming

You are an expert in exact combinatorial optimization, specifically in designing and implementing dynamic programming (DP) algorithms. This skill covers the full design loop — choosing a state space, writing the Bellman recursion, deciding between memoization and tabulation, recovering the optimal solution, and controlling the curse of dimensionality — plus labeling algorithms for resource-constrained shortest paths, the DP family that powers column-generation pricing. Use the framework below to assess whether the problem decomposes, design the state deliberately, implement against the patterns given, and validate against brute force on tiny instances.

Initial Assessment

Before writing any code, establish the following. Each answer changes a design decision downstream.

• Optimal substructure. Can an optimal solution be assembled from optimal solutions of subproblems? If swapping in a better sub-solution can break feasibility or optimality of the whole, DP does not apply directly and the state must be enriched until it does.
• Overlapping subproblems. Count distinct subproblems vs total recursive calls. If every call produces a fresh subproblem (no overlap), plain recursion or branch-and-bound is the right tool; DP buys nothing.
• State-space size, numerically. Multiply out the state dimensions for the target instance size before coding. n=100 items × W=10^9 capacity is 10^11 states — dead on arrival; the same knapsack with W=10^4 is trivial. This single estimate decides feasibility.
• Integer or discretizable data. Pseudo-polynomial DPs (knapsack-style) need integer resource values. If weights/durations are floats, ask what scaling factor is acceptable and what error bound the user needs.
• Value only, or solution too. Recovering the optimal solution costs either full-table memory (parent pointers) or extra recomputation (divide-and-conquer recovery). A bound inside branch-and-bound often needs only the value.
• Standalone or subroutine. A DP called once can afford O(n²) time. A pricing DP called thousands of times inside column generation, or a bound called millions of times inside branch-and-bound, must be lean and warm-startable.
• Acyclic structure. Is there a natural stage ordering (items, periods, nodes of a DAG)? Cyclic state graphs need label-correcting/label-setting treatment or a resource that strictly increases along every transition to guarantee termination.
• Memory budget. A 2D table of 10^8 float64 entries is 800 MB. Decide early whether rolling arrays (value only) suffice or whether full recovery is required.
• Exactness requirement. If the exact state space is too large, is an approximation acceptable? Profit-scaling FPTAS for knapsack (Ibarra & Kim 1975), state-space relaxation for routing, or coarser time discretization all trade accuracy for size — but change what you can claim about the answer.
• One run or a family of runs. If the user will re-solve with slightly different data (duals changing each pricing iteration, capacities varying), structure the code so the instance-independent parts (graph construction, sorting) are reused.
• Python performance ceiling. Pure-Python nested loops handle ~10^6–10^7 state transitions per second; numpy-vectorized inner loops reach 10^8–10^9. Estimate transition count and pick the implementation style accordingly.
• Validation oracle. What independent check exists? Brute-force enumeration up to n≈15–20, a MIP model of the same problem, or known optima from a benchmark library. Plan the cross-check before trusting any DP.

Algorithm Anatomy

Dynamic programming (Bellman 1957, Dynamic Programming) solves a problem by ordering decisions into stages and exploiting the principle of optimality: an optimal policy has the property that, whatever the initial state and first decision, the remaining decisions form an optimal policy for the state that results from the first decision. Formally, with states $s$, actions $A(s)$, transition $\tau(s,a)$, and stage cost $c(s,a)$, the value function satisfies the Bellman recursion

$$ V(s) \;=\; \min{a \in A(s)} \Big\{\, c(s,a) + V\big(\tau(s,a)\big) \,\Big\}, \qquad V(s) = vT(s) \ \text{ for terminal } s , $$

and the optimal solution is recovered by following the argmin actions from the initial state. Maximization flips min to max; nothing else changes.

The five design decisions

| Component | Question it answers | Failure mode if wrong | |---|---|---| | State | What is the minimal information that makes the future independent of the past? | Missing information → wrong answers; extra information → exponential blowup | | Stages / transitions | In what order are decisions made, and what states do they lead to? | Cycles in the state graph → non-terminating recursion | | Recursion | How does $V(s)$ combine stage cost and successor values? | Double-counted or dropped costs | | Base cases | What are the terminal states and their values? | Off-by-one stage errors, infeasible states valued 0 instead of $+\infty$ | | Recovery | How is the optimal solution rebuilt from the table? | Value correct, reported solution inconsistent with it |

State design discipline. The state is correct exactly when two different decision histories that reach the same state always have the same set of optimal completions. Test it on paper: construct two histories mapping to one state and ask whether any future constraint or cost can tell them apart. If yes, the distinguishing quantity belongs in the state. This is the Markov property of the decomposition, and forgetting it is the most common DP bug — for example, a knapsack state of "items considered" without "capacity used", or a shortest-path label without the time-window resource.

Memoization vs tabulation

| Aspect | Memoization (top-down) | Tabulation (bottom-up) | |---|---|---| | States touched | Only those reachable from the initial state | Every state in a fixed enumeration order | | Implementation | Recursion + cache (dict, functools.lru_cache) | Nested loops over the table | | Python performance | Function-call overhead; recursion limit (default 1000 frames) | Fast loops; inner dimension vectorizable with numpy | | Memory control | Hard — the cache keeps everything | Easy — rolling arrays keep only the rows the recursion reads | | Best when | Reachable states are sparse or transitions irregular | The state space is dense and rectangular |

Rule of thumb in Python: prototype top-down to get the recursion right, then convert to bottom-up tabulation with a numpy-vectorized inner dimension for production use. The worked knapsack example below shows the same recursion in both styles.

Complexity landscape

| DP | State | Time | Memory | Note | |---|---|---|---|---| | 0-1 knapsack | (items considered, capacity used) | O(nW) | O(W) value-only | Pseudo-polynomial: W is a magnitude, not an input length | | Wagner-Whitin lot sizing | (first uncovered period) | O(T²), improvable to O(T log T) | O(T) | Wagelmans, van Hoesel & Kolen (1992) | | Held-Karp TSP | (visited subset, last city) | O(n²·2ⁿ) | O(n·2ⁿ) | Exact to n ≈ 20–23 in Python/numpy | | SPPRC labeling | (node, resource vector) | Output-sensitive | Dominance-dependent | Column-generation pricing workhorse (Irnich & Desaulniers 2005) | | DAG shortest path | node | O(V+E) | O(V) | Every finite acyclic DP is a shortest path on its state graph |

Curse of dimensionality. State spaces grow multiplicatively in the number of state dimensions: a vector of d resources each with R levels gives Rᵈ states per node. Remedies, in order of preference: (1) drop a dimension by proving the recursion never needs it; (2) relax the state space (project to fewer dimensions, accept a bound instead of the exact value — see Advanced Techniques); (3) coarsen discretization with a quantified error bound; (4) abandon exactness for approximate DP (Powell 2011, Approximate Dynamic Programming) or a metaheuristic.

When to use DP — and when not

• Use DP standalone when the state-space estimate is comfortably within memory/time (≤ ~10⁸ transitions in numpy), data is integer or cleanly discretizable, and a proof of optimality is wanted. For knapsack with moderate W, uncapacitated lot sizing, sequence alignment, and small-subset problems, DP beats MIP solvers outright.
• Use DP as a subroutine when an outer method needs it: labeling DPs as pricing oracles in column generation, DP over relaxed states as bounding functions in branch-and-bound, DP decoders inside metaheuristics (optimal split of a giant tour, optimal timing for a fixed sequence).
• Prefer a MIP when constraints couple decisions in ways that explode the state (many global side constraints), when data is fractional and scaling is unacceptable, or when the model will keep gaining new constraint types — a MIP absorbs a new constraint in one line; a DP may need a new state dimension.
• Prefer branch-and-bound when subproblems barely overlap but good bounds and dominance rules exist; B&B with a DP-based bound is often the best of both (see branch-and-bound).

Generic DP Engine

The skeleton every finite DP follows, independent of problem:

SOLVE-DP(initial state s0):
  define VALUE(s):
    if s in memo:            return memo[s]
    if s is terminal:        memo[s] = terminal_value(s); return it
    best = +inf  (minimization)
    for each action a in A(s):
      v = c(s, a) + VALUE(tau(s, a))
      if v  tuple[float, list[Action]]:
    """Memoized Bellman recursion over a finite acyclic state graph.

    Returns (optimal value, optimal action sequence from `initial`).
    Requires: every state reaches a terminal state; the state graph is acyclic;
    recursion depth (longest state chain) stays below Python's recursion limit.
    """
    memo: dict[State, tuple[float, Any]] = {}

    def value(s: State) -> float:
        if s in memo:
            return memo[s][0]
        if is_terminal(s):
            memo[s] = (terminal_value(s), None)
            return memo[s][0]
        best_v = math.inf if sense == "min" else -math.inf
        best_a: Any = None
        for a in actions(s):
            v = cost(s, a) + value(transition(s, a))
            if (sense == "min" and v  best_v):
                best_v, best_a = v, a
        memo[s] = (best_v, best_a)
        return best_v

    opt = value(initial)
    plan: list[Action] = []
    s = initial
    while not is_terminal(s):
        a = memo[s][1]
        plan.append(a)
        s = transition(s, a)
    return opt, plan

# Tiny demo: shortest path on a 4-node DAG, arcs[u] = [(successor, cost), ...]
arcs: dict[int, list[tuple[int, float]]] = {
    0: [(1, 2.0), (2, 5.0)],
    1: [(3, 4.0)],
    2: [(3, 2.0)],
}
val, plan = solve_dp(
    initial=0,
    is_terminal=lambda s: s == 3,
    terminal_value=lambda s: 0.0,
    actions=lambda s: [v for v, _ in arcs.get(s, [])],
    transition=lambda s, a: a,
    cost=lambda s, a: dict(arcs[s])[a],
)
print(val, plan)
# Expected: 6.0 [1, 3]  -- path 0 -> 1 -> 3 costs 2 + 4 = 6, beating 0 -> 2 -> 3 = 7

Every finite DP is a shortest (or longest) path on its state graph; if the engine above gives a different answer than your specialized tabulation on small instances, the tabulation has a bug, not the recursion.

Worked Example 1: 0-1 Knapsack

The decomposition. Items $1,\dots,n$ with profits $pi$ and integer weights $wi$, capacity $W$. Define $V(i,w)$ as the best profit achievable using a subset of the first $i$ items with total weight at most $w$ — the state $(i, w)$ records exactly what the remaining items need to know about the past:

$$ V(i, w) \;=\; \max\Big\{\, \underbrace{V(i-1,\, w)}{\text{skip item } i},\;\; \underbrace{V(i-1,\, w - wi) + pi}{\text{take item } i,\ \text{only if } w_i \le w} \,\Big\}, \qquad V(0, w) = 0 . $$

The state is two-dimensional and dense, so bottom-up tabulation with a numpy-vectorized capacity axis is the production implementation. Time O(nW), pseudo-polynomial — see Kellerer, Pferschy & Pisinger (2004), Knapsack Problems, for the full algorithmic landscape.

import numpy as np

def knapsack_dp(profits: np.ndarray, weights: np.ndarray, capacity: int) -> tuple[int, list[int]]:
    """0-1 knapsack by bottom-up tabulation; full table kept for solution recovery.

    profits, weights: integer arrays of length n. Time O(nW), memory O(nW).
    Returns (optimal profit, sorted list of chosen item indices).
    """
    n = len(profits)
    V = np.zeros((n + 1, capacity + 1), dtype=np.int64)
    for i in range(1, n + 1):
        w_i, p_i = int(weights[i - 1]), int(profits[i - 1])
        V[i] = V[i - 1]  # default: skip item i (row assignment copies values)
        if w_i  int:
    """Top-down memoized 0-1 knapsack (value only).

    Explores only reachable (i, w) states; recursion depth equals len(profits).
    """
    sys.setrecursionlimit(max(10_000, len(profits) + 100))

    @lru_cache(maxsize=None)
    def best(i: int, w: int) -> int:
        if i == 0:
            return 0
        skip = best(i - 1, w)
        if weights[i - 1]  int:
    """0-1 knapsack optimal value in O(W) memory with one rolling array.

    Correctness note: `dp[: capacity - w + 1] + p` is materialized before the
    assignment to dp[w:], so every item is used at most once.
    """
    dp = np.zeros(capacity + 1, dtype=np.int64)
    for p, w in zip(profits.tolist(), weights.tolist()):
        if w  tuple[float, list[tuple[int, int]]]:
    """Uncapacitated lot sizing by the Wagner-Whitin forward recursion.

    demand[t]: demand of period t; setup[s]: setup cost if producing in s;
    holding[t]: unit cost of carrying inventory from period t to t+1.
    Returns (optimal cost, production plan) where the plan lists
    (production period, last period covered), all 0-indexed.
    Time O(T^2), memory O(T^2) for the batch-cost matrix.
    """
    T = len(demand)
    c = np.full((T, T), np.inf)  # c[s, t] = setup in s covering demand[s..t]
    for s in range(T):
        c[s, s] = float(setup[s])
        cum_h = 0.0
        for t in range(s + 1, T):
            cum_h += float(holding[t - 1])  # sum of holding[s..t-1]
            c[s, t] = c[s, t - 1] + cum_h * float(demand[t])
    F = np.full(T + 1, np.inf)
    F[0] = 0.0
    pred = np.zeros(T + 1, dtype=np.int64)
    for t in range(1, T + 1):
        cands = F[:t] + c[np.arange(t), t - 1]  # vectorized over setup period s
        s = int(np.argmin(cands))
        F[t], pred[t] = float(cands[s]), s
    plan: list[tuple[int, int]] = []
    t = T
    while t > 0:
        s = int(pred[t])
        plan.append((s, t - 1))
        t = s
    plan.reverse()
    return float(F[T]), plan

demand = np.array([20, 50, 10, 50])
setup = np.array([100, 100, 100, 100])
holding = np.array([1.0, 1.0, 1.0, 1.0])
cost, plan = wagner_whitin(demand, setup, holding)
print(cost, plan)
# Expected: 270.0 [(0, 2), (3, 3)] -- produce 80 units in period 0 (covers periods
# 0-2: setup 100 + holding 1*50 + 2*10 = 170) and 50 units in period 3 (setup 100)

Two practical notes. First, validate against the ZIO property: in the recovered plan, every production period must coincide with zero entering inventory — an independent feasibility check that catches indexing bugs. Second, the planning-horizon theorem (Wagner & Whitin 1958) says that if the optimal $F(t)$ chooses a setup in period $s$, no later period's optimal plan ever needs a setup before $s$; this is what the O(T log T) algorithms exploit, and it also enables rolling-horizon re-solving without recomputing the past. Once capacities enter, the DP state would need the inventory level and the problem becomes NP-hard — switch to the MIP and valid-inequality machinery in lot-sizing.

Labeling Algorithms for Constrained Shortest Paths

The shortest path problem with resource constraints (SPPRC) asks for a minimum-cost source-sink path where each arc also consumes resources

…

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.