Dynamic Programming

dynamic programming is two things wearing one name. to bellman it was a mathematical theory of multistage decision processes — sibling to linear programming in the “programming means planning” sense.¹ 𐃏 to a computer scientist it is a technique: solve a problem by combining solutions to subproblems, and never solve the same subproblem twice (Cormen, Thomas H. and Leiserson, Charles E. and Rivest, Ronald L. and Stein, Clifford, 2009). the two are the same idea at different altitudes, and this page covers both. worked implementations also live in this github repo.

bellman’s principle of optimality#

bellman’s 1957 book states the principle that makes multistage optimisation tractable:

an optimal policy has the property that whatever the initial state and initial decision are, the remaining decisions must constitute an optimal policy with regard to the state resulting from the first decision.

in other words: tails of optimal solutions are optimal solutions of tails. whenever it holds, the value of a state depends only on the state — not on how you got there — and the whole problem collapses into a functional equation. for a finite-horizon deterministic process with state \(s\), decisions \(a \in A(s)\), transition \(s’ = f_k(s, a)\) and stage reward \(r_k(s, a)\):

\begin{equation} V_k(s) = \max_{a \in A(s)} \bigl\{ r_k(s, a) + V_{k+1}(f_k(s, a)) \bigr\}, \qquad V_N(s) = r_N(s). \end{equation}

solve it backwards from the horizon and \(V_0(s_0)\) is the optimal total reward; the maximising \(a\) at each state is the optimal policy. every recurrence on this page is an instance of this equation with different decorations — and when the transition becomes random the same equation grows an expectation, which is where this page ends up.

the two ingredients#

DP applies when a problem has both of (Cormen, Thomas H. and Leiserson, Charles E. and Rivest, Ronald L. and Stein, Clifford, 2009):

• optimal substructure — an optimal solution is built from optimal solutions to subproblems. this is the principle of optimality in CS clothing. it fails when subproblems interact: the longest simple path in a general graph has no such structure (gluing two longest simple paths can revisit vertices), which is one reason it is NP-hard while shortest path is easy.
• overlapping subproblems — the naive recursion meets the same subproblem many times. the total number of distinct subproblems is small (typically polynomial), so caching turns an exponential tree into a polynomial table.

the contrast with divide-and-conquer is exactly the second ingredient. mergesort’s two halves are disjoint — no subproblem recurs, so plain recursion is already optimal and a cache would be dead weight. fibonacci’s two subcalls overlap almost entirely — the naive recursion does \(2F(n+1) - 1\) calls, \(\Theta(\varphi^n)\) with \(\varphi \approx 1.618\), to fill what is really an \((n+1)\)-entry table. 𐃏

memoisation vs tabulation#

two mechanical ways to exploit overlap:

• memoisation (top-down): keep the natural recursion, add a cache keyed by subproblem. pro — only reachable subproblems get solved; the recursion order finds dependencies for you. con — recursion overhead, stack depth, cache-hostile access patterns.
• tabulation (bottom-up): enumerate subproblems in an order where dependencies come first, fill an explicit table. pro — no recursion, tight loops, and the table’s shape exposes space optimisations. con — you must know the dependency order, and you fill every cell whether needed or not.

fibonacci both ways, instrumented:

calls = 0

def fib_naive(n):
    global calls
    calls += 1
    return n if n < 2 else fib_naive(n - 1) + fib_naive(n - 2)

def fib_memo(n, memo=None):
    global calls
    if memo is None:
        memo = {}
    if n in memo:
        return memo[n]
    calls += 1
    memo[n] = n if n < 2 else fib_memo(n - 1, memo) + fib_memo(n - 2, memo)
    return memo[n]

def fib_tab(n):
    table = [0] * (n + 1)
    if n >= 1:
        table[1] = 1
    for i in range(2, n + 1):
        table[i] = table[i - 1] + table[i - 2]
    return table[n]

def fib_rolling(n):
    a, b = 0, 1          # F(i-2), F(i-1)
    for _ in range(n):
        a, b = b, a + b
    return a

for f in (fib_naive, fib_memo):
    calls = 0
    print(f"{f.__name__}(30) = {f(30):>6}   calls = {calls}")
print(f"fib_tab(30)   = {fib_tab(30):>6}   table cells = 31")
print(f"fib_rolling(90) = {fib_rolling(90)}   space = 2 words")

fib_naive(30) = 832040   calls = 2692537
fib_memo(30) = 832040   calls = 31
fib_tab(30)   = 832040   table cells = 31
fib_rolling(90) = 2880067194370816120   space = 2 words

\(2{,}692{,}537 = 2F(31) - 1\) calls collapse to \(31 = n + 1\) — the count of distinct subproblems, exactly as predicted. the rolling version previews the space trick below: the recurrence only ever looks back two cells, so keep two cells. 𐃏

the subproblem dag#

every DP has a hidden graph: one node per distinct subproblem, one edge from each problem to each subproblem it consults. optimal substructure makes the values well-defined; the absence of cycles makes them computable. the two strategies are then just two graph traversals:

• memoisation = depth-first search of this dag from the goal node, caching post-order results.
• tabulation = processing nodes in a (reverse) topological order, so every edge points to an already-filled cell.

running time reads straight off the graph: \(\sum_{\text{nodes}} (\text{work per node})\), which for constant-time combination is \(O(V + E)\) of the subproblem dag (Cormen, Thomas H. and Leiserson, Charles E. and Rivest, Ronald L. and Stein, Clifford, 2009). fibonacci: \(n+1\) nodes, \(2n - 2\) edges, linear. knapsack: \((n+1)(W+1)\) nodes, out-degree \(\le 2\). single-source shortest paths in a dag is this picture with nothing hidden — relax vertices in topological order — and conversely every DP can be read as a shortest (or longest) path computation over its subproblem dag.

worked problems#

0/1 knapsack#

\(n\) items with values \(v_i\) and weights \(w_i\), capacity \(W\); take each item at most once, maximise value. as a mathematical program it is the one-constraint 0/1 integer program \(\max\, v^{\top}x\) subject to \(w^{\top}x \le W\), \(x \in \{0,1\}^n\) — the integer programming page attacks it with branch and bound driven by LP relaxation bounds; DP is the other classic exact method. subproblem: \(\mathrm{OPT}(i, w)\) = best value using only the first \(i\) items with capacity \(w\). the decision for item \(i\) is binary — skip it or take it:

\begin{equation} \mathrm{OPT}(i, w) = \begin{cases} 0 & i = 0 \\ \mathrm{OPT}(i-1,\, w) & w_i > w \\ \max\bigl\{ \mathrm{OPT}(i-1,\, w),\; v_i + \mathrm{OPT}(i-1,\, w - w_i) \bigr\} & w_i \le w \end{cases} \end{equation}

from functools import lru_cache

def knapsack(values, weights, W):
    n = len(values)
    dp = [[0] * (W + 1) for _ in range(n + 1)]
    for i in range(1, n + 1):
        for w in range(W + 1):
            dp[i][w] = dp[i - 1][w]                       # skip item i
            if weights[i - 1] <= w:                        # or take it
                dp[i][w] = max(dp[i][w], dp[i - 1][w - weights[i - 1]] + values[i - 1])
    # traceback: which items were taken?
    taken, w = [], W
    for i in range(n, 0, -1):
        if dp[i][w] != dp[i - 1][w]:
            taken.append(i - 1)
            w -= weights[i - 1]
    return dp[n][W], sorted(taken)

def knapsack_rolling(values, weights, W):
    dp = [0] * (W + 1)
    for v, wt in zip(values, weights):
        for w in range(W, wt - 1, -1):    # reverse! forward would reuse item i
            dp[w] = max(dp[w], dp[w - wt] + v)
    return dp[W]

def knapsack_memo(values, weights, W):
    @lru_cache(maxsize=None)
    def best(i, w):                        # best value using items[0..i) with capacity w
        if i == 0:
            return 0
        skip = best(i - 1, w)
        if weights[i - 1] <= w:
            return max(skip, best(i - 1, w - weights[i - 1]) + values[i - 1])
        return skip
    ans = best(len(values), W)
    return ans, best.cache_info().currsize

values  = [ 60, 100, 120,  80,  30]
weights = [ 10,  20,  30,  15,   5]
W = 50

best, items = knapsack(values, weights, W)
print(f"optimal value = {best}, items taken = {items}")
print(f"weights used  = {[weights[i] for i in items]}, total = {sum(weights[i] for i in items)}")
print(f"rolling-array answer = {knapsack_rolling(values, weights, W)} (O(W) space)")
memo_ans, touched = knapsack_memo(values, weights, W)
print(f"memoised answer = {memo_ans}, states touched = {touched} of {(len(values)+1)*(W+1)} in the full table")

optimal value = 270, items taken = [0, 1, 3, 4]
weights used  = [10, 20, 15, 5], total = 50
rolling-array answer = 270 (O(W) space)
memoised answer = 270, states touched = 36 of 306 in the full table

three points worth staring at:

• complexity is \(O(nW)\) time and space — pseudo-polynomial. it is polynomial in the numeric value of \(W\) but exponential in its bit length, so this does not contradict knapsack being NP-hard. scale-and-round the values and the same DP yields an FPTAS.
• memoisation earned its keep here: only 36 of 306 table states are reachable from \((n, W)\) with these weights, and top-down solved exactly those. tabulation filled all 306 (cheap ones, admittedly).
• the reverse loop in the rolling version matters: iterating \(w\) downwards means dp[w - wt] still holds the previous item’s row. iterate upwards and each item can be taken repeatedly — which is not a bug so much as a different problem: that variant correctly solves unbounded knapsack.

longest common subsequence#

given \(X = x_1 \dots x_m\) and \(Y = y_1 \dots y_n\), find the longest sequence appearing in both in order (not necessarily contiguously). subproblem: \(c(i, j)\) = LCS length of the prefixes \(x_1..x_i\) and \(y_1..y_j\) (Cormen, Thomas H. and Leiserson, Charles E. and Rivest, Ronald L. and Stein, Clifford, 2009):

\begin{equation} c(i, j) = \begin{cases} 0 & i = 0 \text{ or } j = 0 \\ c(i-1,\, j-1) + 1 & x_i = y_j \\ \max\bigl\{ c(i,\, j-1),\; c(i-1,\, j) \bigr\} & x_i \ne y_j \end{cases} \end{equation}

optimal substructure argument: if \(x_i = y_j\), some LCS ends by matching them (exchange argument), and what precedes is an LCS of both prefixes shortened by one. if not, any LCS omits \(x_i\) or omits \(y_j\), so it is an LCS of one of the two smaller problems. \(O(mn)\) time and space; the answer string is recovered by walking the table backwards — no arrows need storing, since each cell’s provenance can be re-derived by comparing it with its neighbours.

def lcs(X, Y):
    m, n = len(X), len(Y)
    dp = [[0] * (n + 1) for _ in range(m + 1)]
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            if X[i - 1] == Y[j - 1]:
                dp[i][j] = dp[i - 1][j - 1] + 1
            else:
                dp[i][j] = max(dp[i - 1][j], dp[i][j - 1])
    # traceback
    out, i, j = [], m, n
    while i and j:
        if X[i - 1] == Y[j - 1]:
            out.append(X[i - 1]); i -= 1; j -= 1
        elif dp[i - 1][j] >= dp[i][j - 1]:
            i -= 1
        else:
            j -= 1
    return dp, "".join(reversed(out))

X, Y = "ABCBDAB", "BDCABA"
dp, s = lcs(X, Y)
print(f"lcs({X!r}, {Y!r}) = {s!r}, length {len(s)}")
print("     " + "  ".join(" " + c for c in " " + Y))
for i, row in enumerate(dp):
    label = (" " + X)[i]
    print(f"  {label}  " + "  ".join(f"{v:2d}" for v in row))

lcs('ABCBDAB', 'BDCABA') = 'BCBA', length 4
          B   D   C   A   B   A
      0   0   0   0   0   0   0
  A   0   0   0   0   1   1   1
  B   0   1   1   1   1   2   2
  C   0   1   1   2   2   2   2
  B   0   1   1   2   2   3   3
  D   0   1   2   2   2   3   3
  A   0   1   2   2   3   3   4
  B   0   1   2   2   3   4   4

the same table, drawn with its traceback:

edit distance#

the levenshtein distance between \(x_1 \dots x_m\) and \(y_1 \dots y_n\): the minimum number of single-character insertions, deletions and substitutions turning one into the other. subproblem \(d(i, j)\) = distance between the two prefixes; the last operation on the alignment’s final column is one of three choices:

\begin{equation} d(i, j) = \begin{cases} \max\{i, j\} & \min\{i, j\} = 0 \\ \min\bigl\{ d(i-1,\, j) + 1,\; d(i,\, j-1) + 1,\; d(i-1,\, j-1) + [x_i \ne y_j] \bigr\} & \text{otherwise} \end{cases} \end{equation}

where \([\cdot]\) is the iverson bracket (a substitution is free when the characters already agree). structurally this is LCS with a min instead of a max and a third diagonal option; both are special cases of sequence alignment with general gap and mismatch costs, which is how bioinformatics consumes this recurrence by the terabyte.

def edit_distance(a, b):
    m, n = len(a), len(b)
    dp = [[0] * (n + 1) for _ in range(m + 1)]
    for i in range(m + 1):
        dp[i][0] = i                      # delete everything
    for j in range(n + 1):
        dp[0][j] = j                      # insert everything
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            sub = 0 if a[i - 1] == b[j - 1] else 1
            dp[i][j] = min(dp[i - 1][j] + 1,          # delete a[i-1]
                           dp[i][j - 1] + 1,          # insert b[j-1]
                           dp[i - 1][j - 1] + sub)    # substitute (or match)
    return dp[m][n]

pairs = [("kitten", "sitting"), ("sunday", "saturday"), ("dynamic", "dynastic")]
for a, b in pairs:
    print(f"edit({a!r}, {b!r}) = {edit_distance(a, b)}")

edit('kitten', 'sitting') = 3
edit('sunday', 'saturday') = 3
edit('dynamic', 'dynastic') = 2

\(O(mn)\) time, and — remarkably — that is probably where it ends: any algorithm running in \(O(n^{2-\epsilon})\) time for some \(\epsilon > 0\) would refute the strong exponential time hypothesis.²

space optimisation#

the subproblem dag tells you what you may forget: once every consumer of a cell has been filled, the cell is garbage. concretely:

• look-back of depth 1: fibonacci keeps two scalars; LCS and edit distance keep two rows (or one row plus one scalar with careful in-place updates), dropping \(O(mn)\) space to \(O(\min\{m, n\})\).
• knapsack in one row: shown above — one \(O(W)\) array traversed in decreasing \(w\), because the recurrence reads strictly left of the write position in the previous row.
• the catch: rolling arrays keep the optimal value but destroy the table you would traceback through. hirschberg’s trick recovers the solution anyway — divide and conquer on the midpoint of one string, computing forward and backward rows to find where the optimal path crosses the middle — same \(O(mn)\) time, linear space.³

dp over time: value and policy iteration#

bellman’s own setting was sequential decision-making, where the modern habitat is the markov decision process: states \(s\), actions \(a\), transition probabilities \(P(s’ \mid s, a)\), rewards, discount \(\gamma \in [0, 1)\). the principle of optimality becomes the bellman optimality equation:

\begin{equation} V^{*}(s) = \max_{a} \Bigl\{ r(s, a) + \gamma \sum_{s’} P(s’ \mid s, a)\, V^{*}(s’) \Bigr\}. \end{equation}

the right-hand side, viewed as an operator \(T\) on value functions, is a \(\gamma\)-contraction in the sup norm, so \(V^{*}\) exists, is unique, and value iteration \(V_{k+1} = T V_k\) converges geometrically: \(\lVert V_{k+1} - V^{*} \rVert_{\infty} \le \gamma \lVert V_k - V^{*} \rVert_{\infty}\). policy iteration alternates exact policy evaluation (a linear solve) with greedy improvement, and terminates in finitely many steps because each iteration strictly improves over a finite policy set. this is DP with the table indexed by state rather than prefix — and the tabulation order replaced by iterate-to-fixpoint, since the “dag” now has cycles that the contraction tames.

import numpy as np

# a 4-state chain: states 0,1,2,3.  state 3 is the goal (absorbing, reward 10 on entry).
# actions: 0 = left, 1 = right.  moves succeed with prob 0.8, stay put with prob 0.2.
# every non-terminal step costs -1.  discount gamma = 0.9.
nS, nA, gamma = 4, 2, 0.9
P = np.zeros((nS, nA, nS))   # P[s, a, s']
R = np.zeros((nS, nA))
for s in range(3):
    for a, d in ((0, -1), (1, +1)):
        t = min(max(s + d, 0), 3)
        P[s, a, t] += 0.8
        P[s, a, s] += 0.2
        R[s, a] = -1 + 0.8 * (10 if t == 3 else 0)
P[3, :, 3] = 1.0             # absorbing, zero reward forever

# value iteration: V <- max_a [ R(s,a) + gamma * sum_s' P(s,a,s') V(s') ]
V = np.zeros(nS)
for k in range(1, 200):
    Q = R + gamma * P @ V              # Q[s,a]
    V_new = Q.max(axis=1)
    if np.abs(V_new - V).max() < 1e-10:
        break
    V = V_new
policy = Q.argmax(axis=1)
print(f"converged after {k} sweeps")
print("V* =", np.round(V, 4))
print("pi* =", ["left" if a == 0 else "right" for a in policy[:3]], "(state 3 terminal)")

# policy iteration: solve the linear system for each fixed policy, then greedify
pi = np.zeros(nS, dtype=int)           # start: always left
for it in range(1, 50):
    Ppi = P[np.arange(nS), pi]         # nS x nS
    Rpi = R[np.arange(nS), pi]
    Vpi = np.linalg.solve(np.eye(nS) - gamma * Ppi, Rpi)   # exact policy evaluation
    pi_new = (R + gamma * P @ Vpi).argmax(axis=1)
    if (pi_new == pi).all():
        break
    pi = pi_new
print(f"policy iteration: {it} iterations, V = {np.round(Vpi, 4)}")

converged after 21 sweeps
V* = [4.2911 6.276  8.5366 0.    ]
pi* = ['right', 'right', 'right'] (state 3 terminal)
policy iteration: 4 iterations, V = [4.2911 6.276  8.5366 0.    ]

both routes agree, and policy iteration gets there in 4 iterations against value iteration’s 21 sweeps — typical, since each of its steps does much more work. everything above assumed the model \(P\) is known; strike that assumption and you must sample the bellman backup from experience instead of computing its expectation, which is precisely q-learning — stochastic-approximation DP.

the same expectation-inside-a-max structure is what stochastic programming calls a multistage problem with recourse: decide now, watch randomness resolve, recover optimally, repeat. stochastic DP solves it by backward induction on states; stochastic programming unrolls the randomness into a scenario tree and hands one enormous deterministic program to a solver. same principle of optimality, opposite trade-off — DP scales in horizon but dies in state dimension (bellman also coined “the curse of dimensionality” for exactly this), while scenario methods tolerate rich state but explode with stages.

see also#

• linear programming — the other “programming” that means planning; simplex vs bellman is 1947 vs 1953
• integer programming — knapsack by branch and bound; where DP’s pseudo-polynomial trick meets NP-hardness
• stochastic programming — recourse problems and scenario trees, the optimisation view of stochastic DP
• classical algorithms — bellman-ford, floyd-warshall and friends: shortest-path DP in the wild
• q-learning — the bellman optimality backup, sampled instead of computed

References#

Cormen, Thomas H. and Leiserson, Charles E. and Rivest, Ronald L. and Stein, Clifford (2009). Introduction to Algorithms, MIT Press.