08 - Dynamic Programming and Memoization Techniques

Compute the thirtieth Fibonacci number with the most obvious recursive definition and your program will make 2,692,537 function calls to return 832,040. The reason is simple and holds for many problems: recursion solves the same subproblems an enormous number of times. Dynamic programming is the technique that solves them only once and remembers the answer. This article applies it to six problems (Fibonacci, stair climbing, 0/1 knapsack, edit distance, common subsequence and increasing subsequence), counts the operations, compares every result with a brute force or an independent reference, and shows two cases in which the greedy choice, the one that seems obvious, goes wrong. You need the complexity from article 1; the comparison with the Floyd-Warshall table comes from article 7.

When a problem lends itself: two conditions

A dynamic solution makes sense if the problem has two properties. The first is optimal substructure: the best solution is obtained by combining best solutions of smaller problems (the best path from A to C through B contains the best path from A to B). The second is overlapping subproblems: the same small problems recur many times. Without the first, recombination is not correct; without the second, there is nothing to save, and plain recursion is enough (as in a binary search, where every subproblem is different).

The term and the method are associated with Richard Bellman, author of Dynamic Programming (Science, 1966; I verified its metadata, I did not read the text).

Memoization: remembering the answers

Memoization (from "memo", a reminder) leaves the recursion intact and adds a dictionary: before solving a subproblem you look whether the answer is already there. The listing compares three versions of Fibonacci. fib_naive counts its own calls in a global variable; fib_memo checks if n not in memo before computing, and counts too; fib_tabella fills from the bottom with no recursion. The last column is an independent formula: the number of calls of the naive version is 2·fib(n+1) - 1, and the comparison with the count is a check of the program, not a definition.

import functools

chiamate = 0

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

def fib_memo(n, memo=None, conta=None):
    memo = {} if memo is None else memo
    conta = [0] if conta is None else conta
    conta[0] += 1
    if n < 2:
        return n, conta[0]
    if n not in memo:  # subproblem n is solved only once
        a, _ = fib_memo(n - 1, memo, conta)
        b, _ = fib_memo(n - 2, memo, conta)
        memo[n] = a + b
    return memo[n], conta[0]

def fib_tabella(n):
    a, b = 0, 1
    for _ in range(n):
        a, b = b, a + b
    return a

print("n    fib(n)        naive calls      memo calls          2*fib(n+1)-1")
for n in (10, 20, 25, 30):
    chiamate = 0
    valore = fib_naive(n)
    v2, c2 = fib_memo(n)
    assert valore == v2 == fib_tabella(n)
    print(f"{n:<4} {valore:<13} {chiamate:<16} {c2:<19} {2 * fib_tabella(n + 1) - 1}")

# Output:
# n    fib(n)        naive calls      memo calls          2*fib(n+1)-1
# 10   55            177              19                  177
# 20   6765          21891            39                  21891
# 25   75025         242785           49                  242785
# 30   832040        2692537          59                  2692537

The naive calls grow like the Fibonacci numbers themselves, about 1.6 times per step (177, 21,891, 242,785, 2,692,537), whereas with memoization they are exactly 2n - 1 (19, 39, 49, 59). The reason is that the call tree stops branching: every value is computed once, and all the other calls that ask for it return immediately. In the standard library the same scheme is a decorator: functools.cache, which the functools documentation describes as an unbounded cache, equal to lru_cache(maxsize=None).

Memoization versus table, and the recursion limit

The two ways of writing the same idea have precise names. Top-down starts from the big problem and descends remembering; bottom-up (or tabulation) starts from the small problems and fills a table up to the big one. The difference is not aesthetic: the top-down version uses the call stack, so it has Python's recursion ceiling; the bottom-up one does not. The listing counts in how many ways you can climb a staircase of n steps with steps of 1 or 2, which has the same recurrence as Fibonacci.

import functools

@functools.cache
def scale(n):
    """In how many ways you climb n steps taking steps of 1 or 2."""
    return 1 if n <= 1 else scale(n - 1) + scale(n - 2)

def scale_tabella(n):
    a, b = 1, 1
    for _ in range(n - 1):
        a, b = b, a + b
    return b if n > 1 else 1

print("scale(40) memo:", scale(40), " table:", scale_tabella(40))
try:
    scale(5000)
except RecursionError as e:
    print("scale(5000) recursive memo:", type(e).__name__ + ":", e)
r = scale_tabella(5000)
print("scale(5000) table: a number of", len(str(r)), "digits")

# Output:
# scale(40) memo: 165580141  table: 165580141
# scale(5000) recursive memo: RecursionError: maximum recursion depth exceeded
# scale(5000) table: a number of 1045 digits

With 40 steps the two versions coincide (165,580,141 ways). With 5,000 the memoized version fails with RecursionError, because to fill the cache it must descend 5,000 levels before it can climb back; the table returns a number of 1,045 digits without problems. Raising the limit with sys.setrecursionlimit moves the problem and risks crashing the process; the table is the robust choice when the depth grows with the input. When the recursion stays shallow, on the other hand, memoization is written in two lines and is easier to read.

0/1 knapsack: choosing which items to carry

You have a knapsack that holds a maximum weight and items each with a weight and a value; each item is taken whole or left. Which set is worth the most? The first instinct is to take the items by decreasing value per kilo. With three items (weight 10 and value 60, weight 20 and value 100, weight 30 and value 120) and capacity 50, the greedy takes the first two (the best ratio) and totals 160, but the best is to choose the last two and reach 220. The problem has optimal substructure: the best value with the first i items and capacity c is the maximum between not taking item i and taking it (its value plus the best with i - 1 items and reduced capacity).

In the listing tab[i][c] holds exactly that value. The line tab[i][c] = tab[i - 1][c] is the "I do not take it", the branch with tab[i - 1][c - p] + v is the "I take it", and the final loop walks back up the table to rebuild which items were chosen: if tab[i][c] differs from tab[i - 1][c], item i is in the solution. To verify, forza_bruta tries all the subsets on 3,000 random instances with at most 9 items, and we count how many times the greedy stays below the optimum.

import itertools
import random

def zaino(oggetti, capacita):
    """oggetti = [(weight, value)]. tab[i][c] = best value with the first i items and capacity c."""
    n = len(oggetti)
    tab = [[0] * (capacita + 1) for _ in range(n + 1)]
    celle = 0
    for i, (p, v) in enumerate(oggetti, start=1):
        for c in range(capacita + 1):
            celle += 1
            tab[i][c] = tab[i - 1][c]  # I do not take item i
            if p <= c and tab[i - 1][c - p] + v > tab[i][c]:
                tab[i][c] = tab[i - 1][c - p] + v  # I take it
    scelti, c = [], capacita
    for i in range(n, 0, -1):  # walks back up the table to rebuild the choice
        if tab[i][c] != tab[i - 1][c]:
            scelti.append(i - 1)
            c -= oggetti[i - 1][0]
    return tab[n][capacita], sorted(scelti), celle

def forza_bruta(oggetti, capacita):
    best = 0
    for r in range(len(oggetti) + 1):
        for gruppo in itertools.combinations(oggetti, r):
            if sum(p for p, _ in gruppo) <= capacita:
                best = max(best, sum(v for _, v in gruppo))
    return best

def goloso(oggetti, capacita):  # by decreasing value/weight
    tot = 0
    for p, v in sorted(oggetti, key=lambda x: x[1] / x[0], reverse=True):
        if p <= capacita:
            capacita -= p
            tot += v
    return tot

esempio = [(10, 60), (20, 100), (30, 120)]  # (weight, value), capacity 50
print("greedy:", goloso(esempio, 50), " DP:", zaino(esempio, 50)[:2], " brute force:", forza_bruta(esempio, 50))

rng = random.Random(12)
goloso_sbaglia = 0
for prova in range(3000):
    ogg = [(rng.randint(1, 12), rng.randint(1, 30)) for _ in range(rng.randint(1, 9))]
    cap = rng.randint(1, 40)
    ref = forza_bruta(ogg, cap)
    assert zaino(ogg, cap)[0] == ref
    goloso_sbaglia += goloso(ogg, cap) < ref
print("3000 instances: DP equal to brute force; greedy is below the optimum in", goloso_sbaglia)

print("n    capacity   cells    n*(capacity+1)")
for n, cap in ((20, 1000), (20, 2000), (20, 4000)):
    ogg = [(rng.randint(1, 100), rng.randint(1, 100)) for _ in range(n)]
    print(f"{n:<4} {cap:<10} {zaino(ogg, cap)[2]:<8} {n * (cap + 1)}")

# Output:
# greedy: 160  DP: (220, [1, 2])  brute force: 220
# 3000 instances: DP equal to brute force; greedy is below the optimum in 379
# n    capacity   cells    n*(capacity+1)
# 20   1000       20020    20020
# 20   2000       40020    40020
# 20   4000       80020    80020

The example gives 160 to the greedy and 220 to dynamic programming, with items 1 and 2 (counting from 0), and brute force confirms 220. On 3,000 random instances the table always coincides with brute force and the greedy stays below the optimum in 379 cases: 12.6% of the instances I generated, a number that depends on the way I generated them and is not a property of the problem. The step count is exactly n times (capacity + 1): 20,020, 40,020, 80,020 as the capacity doubles.

This last line contains a subtlety that is worth stating: the cost grows with the numeric value of the capacity, not with the space needed to write it. Writing 4,000 takes 12 bits, writing 4,000,000 takes 22, but the table becomes a thousand times larger. An algorithm of this form is called pseudo-polynomial: fast when the numbers involved are small, unusable when they are large.

Edit distance: how many changes separate two words

The edit (or Levenshtein) distance between two strings is the minimum number of insertions, deletions and substitutions of a character that transform the first into the second. It is the kind of computation with which a spell checker orders its suggestions, and the classic example is the pair "kitten" and "sitting". The Wagner and Fischer algorithm, described in The String-to-String Correction Problem (Journal of the ACM, 1974; metadata verified, text not read), fills a table in which d[i][j] is the distance between the first i characters of the first string and the first j of the second. Each cell is the minimum of three moves: delete (cell above plus 1), insert (cell to the left plus 1), substitute or accept (diagonal, plus 1 if the characters differ).

The independent reference is another view of the problem: the strings are the vertices of a graph, two strings are linked if they differ by one edit, and the edit distance is the distance in edges that the BFS of article 6 computes. The listing compares the two routes on all 961 pairs of strings over the alphabet {a, b} of length at most 4, and counts the cells of the table for three sizes.

import itertools
from collections import deque

def edit(a, b):
    """Levenshtein distance (Wagner-Fischer): d[i][j] = distance between a[:i] and b[:j]."""
    d = [[0] * (len(b) + 1) for _ in range(len(a) + 1)]
    for i in range(len(a) + 1):
        d[i][0] = i  # delete the characters
    for j in range(len(b) + 1):
        d[0][j] = j  # insert the characters
    for i in range(1, len(a) + 1):
        for j in range(1, len(b) + 1):
            costo = 0 if a[i - 1] == b[j - 1] else 1
            d[i][j] = min(d[i - 1][j] + 1, d[i][j - 1] + 1, d[i - 1][j - 1] + costo)
    return d[len(a)][len(b)], (len(a) + 1) * (len(b) + 1)

def vicini(s, alfabeto, massimo):
    for i in range(len(s) + 1):
        for c in alfabeto:
            if len(s) < massimo:
                yield s[:i] + c + s[i:]
    for i in range(len(s)):
        yield s[:i] + s[i + 1:]
        for c in alfabeto:
            yield s[:i] + c + s[i + 1:]

def distanze_bfs(sorgente, alfabeto, massimo):
    dist, coda = {sorgente: 0}, deque([sorgente])
    while coda:
        u = coda.popleft()
        for v in vicini(u, alfabeto, massimo):
            if v not in dist:
                dist[v] = dist[u] + 1
                coda.append(v)
    return dist

print("kitten -> sitting:", edit("kitten", "sitting")[0])
# independent reference: BFS on the graph of strings, one edge per edit
stringhe = ["".join(p) for r in range(5) for p in itertools.product("ab", repeat=r)]
coppie = 0
for a in stringhe:
    ref = distanze_bfs(a, "ab", 8)
    for b in stringhe:
        assert edit(a, b)[0] == ref[b]
        coppie += 1
print("pairs of strings over {a,b}, length <= 4, compared with the BFS:", coppie)
print("lengths     table cells")
for m, n in ((10, 10), (100, 100), (1000, 1000)):
    print(f"{m}x{n:<7} {edit('a' * m, 'b' * n)[1]}")

# Output:
# kitten -> sitting: 3
# pairs of strings over {a,b}, length <= 4, compared with the BFS: 961
# lengths     table cells
# 10x10      121
# 100x100     10201
# 1000x1000    1002001

"kitten" and "sitting" are 3 apart (two substitutions and one insertion), and all 961 pairs coincide with the graph reference. The cells are (m + 1) times (n + 1): 121, 10,201 and 1,002,001. The quadratic cost is real: to compare two texts of a thousand characters you need a million cells, and for a million characters a thousand billion, which is why programs that compare long files use variants with restrictions. In memory, moreover, you need the whole table only if you want to rebuild the edits; for the number alone, two rows at a time are enough.

Common subsequence and increasing subsequence

A subsequence is obtained by deleting characters without changing the order of the remaining ones. The longest common subsequence (LCS) of two strings is a problem close to comparing two versions of a text: what belongs to the common subsequence stays the same, the rest is difference. The recurrence has the shape of the edit table: if the last characters coincide, the LCS grows by one with respect to the diagonal; otherwise it is the maximum of the two adjacent cells.

The longest increasing subsequence (LIS) of a list of numbers can be computed in two ways. The quadratic version compares each element with all the previous ones; the fast version keeps a list in which code[k] is the smallest final value of an increasing subsequence of length k + 1, and for each element does a binary search with bisect_left. The listing compares both and the LCS with a brute force that enumerates all the subsequences (1,500 cases per type), then counts the comparisons.

import bisect
import itertools
import math
import random

def lcs(a, b):
    t = [[0] * (len(b) + 1) for _ in range(len(a) + 1)]
    for i in range(1, len(a) + 1):
        for j in range(1, len(b) + 1):
            t[i][j] = t[i - 1][j - 1] + 1 if a[i - 1] == b[j - 1] else max(t[i - 1][j], t[i][j - 1])
    return t[-1][-1]

def lcs_brute(a, b):
    sotto = lambda s: {"".join(c) for r in range(len(s) + 1) for c in itertools.combinations(s, r)}
    return max(len(x) for x in sotto(a) & sotto(b))

def lis_quadratica(v):
    best, confronti = [1] * len(v), 0
    for i in range(len(v)):
        for j in range(i):
            confronti += 1
            if v[j] < v[i]:
                best[i] = max(best[i], best[j] + 1)
    return max(best, default=0), confronti

def lis_veloce(v):
    code = []  # code[k] = smallest tail of an increasing subsequence of length k+1
    for x in v:
        k = bisect.bisect_left(code, x)  # one binary search per element
        if k == len(code):
            code.append(x)
        else:
            code[k] = x
    return len(code)

def lis_brute(v):
    return max(r for r in range(len(v) + 1)
               if any(all(c[i] < c[i + 1] for i in range(r - 1)) for c in itertools.combinations(v, r)))

rng = random.Random(21)
for prova in range(1500):
    a = "".join(rng.choice("abc") for _ in range(rng.randint(0, 8)))
    b = "".join(rng.choice("abc") for _ in range(rng.randint(0, 8)))
    assert lcs(a, b) == lcs_brute(a, b)
for prova in range(1500):
    v = [rng.randint(0, 9) for _ in range(rng.randint(0, 10))]
    assert lis_quadratica(v)[0] == lis_veloce(v) == lis_brute(v)
print("LCS and LIS: 1500 + 1500 cases equal to brute force; lcs('AGGTAB','GXTXAYB') =", lcs("AGGTAB", "GXTXAYB"))

print("n      quadratic LIS compares     n*(n-1)/2    binary searches   bound on compares     same LIS")
for n in (500, 1000, 2000):
    v = [rng.randint(0, 10**6) for _ in range(n)]
    q, conf = lis_quadratica(v)
    limite = n * math.ceil(math.log2(q + 1))  # n searches on a list at most q long
    print(f"{n:<6} {conf:<26} {n * (n - 1) // 2:<12} {n:<17} {limite:<21} {q == lis_veloce(v)} (length {q})")

# Output:
# LCS and LIS: 1500 + 1500 cases equal to brute force; lcs('AGGTAB','GXTXAYB') = 4
# n      quadratic LIS compares     n*(n-1)/2    binary searches   bound on compares     same LIS
# 500    124750                     124750       500               3000                  True (length 39)
# 1000   499500                     499500       1000              6000                  True (length 59)
# 2000   1999000                    1999000      2000              14000                 True (length 91)

The LCS of "AGGTAB" and "GXTXAYB" is 4 and all 3,000 cases agree with brute force. The quadratic version performs exactly n(n - 1)/2 comparisons (124,750, 499,500, 1,999,000 for 500, 1,000 and 2,000 elements); the fast one does n binary searches on a list at most as long as the answer (39, 59, 91), that is at most n·⌈lg(L + 1)⌉ comparisons: 3,000, 6,000 and 14,000. The two versions give the same length in all three cases. The gain is two orders of magnitude at 2,000 elements, and grows with n.

When the greedy choice is enough and when it is not

The knapsack has already shown a greedy that goes wrong. The most instructive case is that of change-making: giving change with the minimum number of coins. With coins of 1, 3 and 4 and an amount of 6 the greedy (the largest coin that fits) gives 4 + 1 + 1, three coins, whereas 3 + 3 uses two. The listing compares greedy and dynamic programming on all the amounts from 1 to 100 with this coin system, and then on all the amounts from 1 to 1,000 cents with the euro coins.

def goloso(monete, importo):
    n = 0
    for m in sorted(monete, reverse=True):
        n += importo // m
        importo %= m
    return n

def dp(monete, massimo):
    """minimum number of coins for each amount 0..massimo (INF if impossible)."""
    INF = float("inf")
    best = [0] + [INF] * massimo
    for x in range(1, massimo + 1):
        for m in monete:
            if m <= x and best[x - m] + 1 < best[x]:
                best[x] = best[x - m] + 1
    return best

print("coins {1,3,4}, amount 6: greedy", goloso([1, 3, 4], 6), " DP", dp([1, 3, 4], 6)[6])
tab = dp([1, 3, 4], 100)
print("amounts 1..100 with {1,3,4} where greedy is worse than DP:", sum(goloso([1, 3, 4], x) > tab[x] for x in range(1, 101)))
euro = [1, 2, 5, 10, 20, 50, 100, 200]
tab = dp(euro, 1000)
print("euro cents, amounts 1..1000 where greedy is worse than DP:", sum(goloso(euro, x) > tab[x] for x in range(1, 1001)))

# Output:
# coins {1,3,4}, amount 6: greedy 3  DP 2
# amounts 1..100 with {1,3,4} where greedy is worse than DP: 24
# euro cents, amounts 1..1000 where greedy is worse than DP: 0

With the coins {1, 3, 4} the greedy is worse than the DP on 24 amounts out of 100. With the euro coins (1, 2, 5, 10, 20, 50, 100, 200 cents) it is so on zero amounts out of 1,000: I verified this interval, I did not prove that it always holds; the {1, 3, 4} system shows that for an arbitrary system it is not guaranteed. The rule that follows is to state what was checked: if you cannot prove that the greedy is optimal for your problem, the table costs little and gives you the right answer.

What you can do now, and where the series arrives

You can recognize a dynamic problem from the two conditions, write the recurrence, implement it with memoization or with a table choosing according to the depth, and count the cells to state the cost. You can verify a result with a brute force or with a reference built in a different way, and you can build the counterexample that unmasks the greedy choice. With this article the series on algorithms closes: the sequence started from measurement, built the structures and closed with the technique that reuses results.

• Article 1: Big-O notation and complexity analysis, the language for stating costs.
• Article 2: sorting algorithms, bubble, merge, quick and heap.
• Article 3: linear and binary search, and hash tables.
• Article 4: lists, stacks, queues and deques.
• Article 5: binary trees, BST and AVL.
• Article 6: graphs, representation and BFS and DFS traversals.
• Article 7: shortest paths, Dijkstra, Bellman-Ford and Floyd-Warshall.
• Article 8: this page, dynamic programming and memoization.

The list follows the reading order and that order makes sense: the measurement of complexity (1) serves to judge everything else; sorting and search (2-3) and the structures (4-5) provide the bricks; graph traversals (6) and shortest paths (7) use them for network problems; dynamic programming (8) brings back the idea of the Floyd-Warshall table in problems without a graph. If you restart from a single point, restart from the first: the O notation is presupposed by every count in the other seven.

The full list of the eight articles is on the series page.