08 - Programmazione dinamica e memoization

Calcola il trentesimo numero di Fibonacci con la definizione ricorsiva più ovvia e il tuo programma farà 2.692.537 chiamate di funzione per restituire 832.040. Il motivo è semplice e vale per molti problemi: la ricorsione risolve gli stessi sottoproblemi un numero enorme di volte. La programmazione dinamica è la tecnica che li risolve una volta sola e ricorda la risposta. Questo articolo la applica a sei problemi (Fibonacci, scalata di gradini, zaino 0/1, distanza di edit, sottosequenza comune e sottosequenza crescente), conta le operazioni, confronta ogni risultato con una forza bruta o un riferimento indipendente, e mostra due casi in cui la scelta golosa, quella che sembra ovvia, sbaglia. Ti servono la complessità dell'articolo 1; il paragone con la tabella di Floyd-Warshall viene dall'articolo 7.

Quando un problema si presta: due condizioni

Una soluzione dinamica ha senso se il problema ha due proprietà. La prima è la sottostruttura ottima: la soluzione migliore si ottiene combinando soluzioni migliori di problemi più piccoli (il miglior percorso da A a C passando per B contiene il miglior percorso da A a B). La seconda è la sovrapposizione dei sottoproblemi: gli stessi problemi piccoli ricorrono molte volte. Senza la prima la ricombinazione non è corretta; senza la seconda non c'è niente da risparmiare, e basta la ricorsione semplice (come in una ricerca binaria, dove ogni sottoproblema è diverso).

Il termine e il metodo sono associati a Richard Bellman, autore di Dynamic Programming (Science, 1966; ne ho verificato i metadati, non letto il testo).

Memoization: ricordare le risposte

La memoization (da «memo», promemoria) lascia intatta la ricorsione e aggiunge un dizionario: prima di risolvere un sottoproblema si guarda se la risposta è già lì. Il listato mette a confronto tre versioni di Fibonacci. fib_naive conta le proprie chiamate in una variabile globale; fib_memo controlla if n not in memo prima di calcolare, e conta anch'essa; fib_tabella riempie dal basso senza ricorsione. L'ultima colonna è una formula indipendente: il numero di chiamate della versione naive è 2·fib(n+1) - 1, e il confronto con il conteggio è una verifica del programma, non una definizione.

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:  # il sottoproblema n si risolve una volta sola
        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)        chiamate naive   chiamate con memo   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)        chiamate naive   chiamate con memo   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

Le chiamate naive crescono come i numeri di Fibonacci stessi, circa 1,6 volte per ogni passo (177, 21.891, 242.785, 2.692.537), mentre con la memoizzazione sono esattamente 2n - 1 (19, 39, 49, 59). Il motivo è che l'albero delle chiamate smette di diramarsi: ogni valore viene calcolato una volta, e tutte le altre chiamate che lo richiedono tornano subito. Nella libreria standard lo stesso schema è un decoratore: functools.cache, che la documentazione di functools descrive come una cache senza limiti, uguale a lru_cache(maxsize=None).

Memoization contro tabella, e il limite di ricorsione

I due modi di scrivere la stessa idea hanno nomi precisi. Dall'alto (top-down) si parte dal problema grande e si scende ricordando; dal basso (bottom-up, o tabulazione) si parte dai problemi piccoli e si riempie una tabella fino a quello grande. La differenza non è estetica: la versione dall'alto usa lo stack delle chiamate, quindi ha il tetto di ricorsione di Python; quella dal basso no. Il listato conta in quanti modi si sale una scala di n gradini con passi da 1 o da 2, che ha la stessa ricorrenza di Fibonacci.

import functools

@functools.cache
def scale(n):
    """In quanti modi sali n gradini facendo passi da 1 o da 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), " tabella:", scale_tabella(40))
try:
    scale(5000)
except RecursionError as e:
    print("scale(5000) memo ricorsiva:", type(e).__name__ + ":", e)
r = scale_tabella(5000)
print("scale(5000) tabella: numero di", len(str(r)), "cifre")

# Output:
# scale(40) memo: 165580141  tabella: 165580141
# scale(5000) memo ricorsiva: RecursionError: maximum recursion depth exceeded
# scale(5000) tabella: numero di 1045 cifre

Con 40 gradini le due versioni coincidono (165.580.141 modi). Con 5.000 la versione memoizzata fallisce con RecursionError, perché per riempire la cache deve scendere di 5.000 livelli prima di poter risalire; la tabella restituisce un numero di 1.045 cifre senza problemi. Alzare il limite con sys.setrecursionlimit sposta il problema e rischia di far crollare il processo; la tabella è la scelta robusta quando la profondità cresce con l'ingresso. Quando la ricorsione resta bassa, invece, la memoization si scrive in due righe ed è più facile da leggere.

Zaino 0/1: scegliere quali oggetti portare

Hai uno zaino che regge un peso massimo e degli oggetti ciascuno con un peso e un valore; ogni oggetto si prende intero o si lascia. Quale insieme vale di più? Il primo istinto è prendere gli oggetti per valore al chilo decrescente. Con tre oggetti (peso 10 e valore 60, peso 20 e valore 100, peso 30 e valore 120) e capacità 50, il goloso prende i primi due (il rapporto migliore) e totalizza 160, ma il meglio è scegliere gli ultimi due e arrivare a 220. Il problema ha sottostruttura ottima: il miglior valore con i primi i oggetti e capacità c è il massimo fra non prendere l'oggetto i e prenderlo (valore più il meglio con i - 1 oggetti e capacità ridotta).

Nel listato tab[i][c] contiene proprio quel valore. La riga tab[i][c] = tab[i - 1][c] è il «non lo prendo», il ramo con tab[i - 1][c - p] + v è il «lo prendo», e il ciclo finale risale la tabella per ricostruire quali oggetti sono stati scelti: se tab[i][c] differisce da tab[i - 1][c], l'oggetto i è nella soluzione. Per verificare, forza_bruta prova tutti i sottoinsiemi su 3.000 istanze casuali con al più 9 oggetti, e si conta quante volte il goloso resta sotto l'ottimo.

import itertools
import random

def zaino(oggetti, capacita):
    """oggetti = [(peso, valore)]. tab[i][c] = miglior valore con i primi i oggetti e capacita' 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]  # non prendo l'oggetto i
            if p <= c and tab[i - 1][c - p] + v > tab[i][c]:
                tab[i][c] = tab[i - 1][c - p] + v  # lo prendo
    scelti, c = [], capacita
    for i in range(n, 0, -1):  # risale la tabella per ricostruire la scelta
        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):  # per valore/peso decrescente
    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)]  # (peso, valore), capacita' 50
print("goloso:", goloso(esempio, 50), " DP:", zaino(esempio, 50)[:2], " forza bruta:", 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 istanze: DP uguale alla forza bruta; il goloso e' sotto l'ottimo in", goloso_sbaglia)

print("n    capacita   celle    n*(capacita+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:
# goloso: 160  DP: (220, [1, 2])  forza bruta: 220
# 3000 istanze: DP uguale alla forza bruta; il goloso e' sotto l'ottimo in 379
# n    capacita   celle    n*(capacita+1)
# 20   1000       20020    20020
# 20   2000       40020    40020
# 20   4000       80020    80020

L'esempio dà 160 al goloso e 220 alla programmazione dinamica, con gli oggetti 1 e 2 (a partire da 0), e la forza bruta conferma 220. Su 3.000 istanze casuali la tabella coincide sempre con la forza bruta e il goloso resta sotto l'ottimo in 379 casi: il 12,6% delle istanze che ho generato, un numero che dipende dal modo in cui le ho generate e non è una proprietà del problema. Il conteggio dei passi è esattamente n per (capacità + 1): 20.020, 40.020, 80.020 al raddoppio della capacità.

Quest'ultima riga contiene una sottigliezza che conviene dichiarare: il costo cresce con il valore numerico della capacità, non con lo spazio che occupa per scriverla. Scrivere 4.000 richiede 12 bit, scrivere 4.000.000 ne richiede 22, ma la tabella diventa mille volte più grande. Un algoritmo con questa forma si dice pseudo-polinomiale: veloce quando i numeri in gioco sono piccoli, inutilizzabile quando sono grandi.

Distanza di edit: quante modifiche separano due parole

La distanza di edit (o di Levenshtein) fra due stringhe è il minimo numero di inserimenti, cancellazioni e sostituzioni di un carattere che trasformano la prima nella seconda. È il genere di calcolo con cui un correttore ortografico ordina i suggerimenti, e l'esempio classico è la coppia «kitten» e «sitting». L'algoritmo di Wagner e Fischer, descritto in The String-to-String Correction Problem (Journal of the ACM, 1974; metadati verificati, testo non letto), riempie una tabella in cui d[i][j] è la distanza fra i primi i caratteri della prima stringa e i primi j della seconda. Ogni cella è il minimo di tre mosse: cancellare (cella sopra più 1), inserire (cella a sinistra più 1), sostituire o accettare (diagonale, più 1 se i caratteri differiscono).

Il riferimento indipendente è un'altra visione del problema: le stringhe sono i vertici di un grafo, due stringhe sono collegate se differiscono per una modifica, e la distanza di edit è la distanza in archi che la BFS dell'articolo 6 calcola. Il listato confronta le due strade su tutte le 961 coppie di stringhe sull'alfabeto {a, b} lunghe al più 4, e conta le celle della tabella per tre dimensioni.

import itertools
from collections import deque

def edit(a, b):
    """Distanza di Levenshtein (Wagner-Fischer): d[i][j] = distanza fra a[:i] e b[:j]."""
    d = [[0] * (len(b) + 1) for _ in range(len(a) + 1)]
    for i in range(len(a) + 1):
        d[i][0] = i  # cancellare i caratteri
    for j in range(len(b) + 1):
        d[0][j] = j  # inserire i caratteri
    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])
# riferimento indipendente: BFS sul grafo delle stringhe, un arco per ogni modifica
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("coppie di stringhe su {a,b}, lunghezza <= 4, confrontate con la BFS:", coppie)
print("lunghezze   celle della tabella")
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
# coppie di stringhe su {a,b}, lunghezza <= 4, confrontate con la BFS: 961
# lunghezze   celle della tabella
# 10x10      121
# 100x100     10201
# 1000x1000    1002001

«kitten» e «sitting» distano 3 (due sostituzioni e un inserimento), e tutte le 961 coppie coincidono con il riferimento a grafo. Le celle sono (m + 1) per (n + 1): 121, 10.201 e 1.002.001. Il costo quadratico è reale: per confrontare due testi di mille caratteri servono un milione di celle, e per un milione di caratteri mille miliardi, ed è per questo che i programmi che confrontano file lunghi usano varianti con restrizioni. In memoria, poi, ti serve l'intera tabella solo se vuoi ricostruire le modifiche; per il solo numero bastano due righe alla volta.

Sottosequenza comune e sottosequenza crescente

Una sottosequenza si ottiene cancellando caratteri senza cambiare l'ordine di quelli rimasti. La più lunga sottosequenza comune (LCS) di due stringhe è un problema vicino al confronto fra due versioni di un testo: ciò che appartiene alla sottosequenza comune resta uguale, il resto è differenza. La ricorrenza ha la forma della tabella di edit: se gli ultimi caratteri coincidono, la LCS cresce di uno rispetto alla diagonale; altrimenti è il massimo fra le due celle adiacenti.

La più lunga sottosequenza crescente (LIS) di una lista di numeri si può calcolare in due modi. La versione quadratica confronta ogni elemento con tutti i precedenti; la versione veloce tiene una lista in cui code[k] è il valore finale più piccolo di una sottosequenza crescente lunga k + 1, e per ogni elemento fa una ricerca binaria con bisect_left. Il listato confronta entrambe e la LCS con una forza bruta che enumera tutte le sottosequenze (1.500 casi per tipo), poi conta i confronti.

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] = coda minima di una sottosequenza crescente lunga k+1
    for x in v:
        k = bisect.bisect_left(code, x)  # una ricerca binaria per elemento
        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 e LIS: 1500 + 1500 casi uguali alla forza bruta; lcs('AGGTAB','GXTXAYB') =", lcs("AGGTAB", "GXTXAYB"))

print("n      confronti LIS quadratica   n*(n-1)/2    ricerche binarie  limite sui confronti  stessa 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 ricerche su una lista lunga al piu' q
    print(f"{n:<6} {conf:<26} {n * (n - 1) // 2:<12} {n:<17} {limite:<21} {q == lis_veloce(v)} (lunghezza {q})")

# Output:
# LCS e LIS: 1500 + 1500 casi uguali alla forza bruta; lcs('AGGTAB','GXTXAYB') = 4
# n      confronti LIS quadratica   n*(n-1)/2    ricerche binarie  limite sui confronti  stessa LIS
# 500    124750                     124750       500               3000                  True (lunghezza 39)
# 1000   499500                     499500       1000              6000                  True (lunghezza 59)
# 2000   1999000                    1999000      2000              14000                 True (lunghezza 91)

La LCS di «AGGTAB» e «GXTXAYB» è 4 e tutti i 3.000 casi concordano con la forza bruta. La versione quadratica esegue esattamente n(n - 1)/2 confronti (124.750, 499.500, 1.999.000 per 500, 1.000 e 2.000 elementi); la veloce fa n ricerche binarie su una lista lunga al più quanto la risposta (39, 59, 91), cioè al più n·⌈lg(L + 1)⌉ confronti: 3.000, 6.000 e 14.000. Le due versioni danno la stessa lunghezza in tutti e tre i casi. Il guadagno è di due ordini di grandezza a 2.000 elementi, e cresce con n.

Quando la scelta golosa basta e quando no

Lo zaino ha già mostrato un goloso che sbaglia. Il caso più istruttivo è quello dei resti: dare il resto con il minimo di monete. Con monete da 1, 3 e 4 e un importo di 6 il goloso (la moneta più grande che entra) dà 4 + 1 + 1, tre monete, mentre 3 + 3 ne usa due. Il listato confronta goloso e programmazione dinamica su tutti gli importi da 1 a 100 con questo sistema di monete, e poi su tutti gli importi da 1 a 1.000 centesimi con le monete in euro.

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

def dp(monete, massimo):
    """minimo di monete per ogni importo 0..massimo (INF se impossibile)."""
    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("monete {1,3,4}, importo 6: goloso", goloso([1, 3, 4], 6), " DP", dp([1, 3, 4], 6)[6])
tab = dp([1, 3, 4], 100)
print("importi 1..100 con {1,3,4} dove il goloso e' peggio del 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("centesimi di euro, importi 1..1000 dove il goloso e' peggio del DP:", sum(goloso(euro, x) > tab[x] for x in range(1, 1001)))

# Output:
# monete {1,3,4}, importo 6: goloso 3  DP 2
# importi 1..100 con {1,3,4} dove il goloso e' peggio del DP: 24
# centesimi di euro, importi 1..1000 dove il goloso e' peggio del DP: 0

Con le monete {1, 3, 4} il goloso è peggiore del DP in 24 importi su 100. Con le monete in euro (1, 2, 5, 10, 20, 50, 100, 200 centesimi) lo è in zero importi su 1.000: ho verificato questo intervallo, non dimostrato che valga sempre; il sistema {1, 3, 4} mostra che per un sistema qualunque non è garantito. La regola che ne segue è dichiarare cosa si è controllato: se non puoi dimostrare che il goloso è ottimo per il tuo problema, la tabella costa poco e ti dà la risposta giusta.

Che cosa sai fare adesso, e dove arriva la serie

Sai riconoscere un problema dinamico dalle due condizioni, scrivere la ricorrenza, implementarla con memoization o con tabella scegliendo in base alla profondità, e contare le celle per dichiarare il costo. Sai verificare un risultato con una forza bruta o con un riferimento costruito in modo diverso, e sai costruire il controesempio che smaschera la scelta golosa. Con questo articolo la serie sugli algoritmi si chiude: la sequenza è partita dalla misura, ha costruito le strutture e ha chiuso con la tecnica che riusa i risultati.

• Articolo 1: Big-O notation e analisi della complessità, il linguaggio per dichiarare i costi.
• Articolo 2: algoritmi di ordinamento, bubble, merge, quick e heap.
• Articolo 3: ricerca lineare, binaria e tabelle hash.
• Articolo 4: liste, stack, queue e deque.
• Articolo 5: alberi binari, BST e AVL.
• Articolo 6: grafi, rappresentazione e visite BFS e DFS.
• Articolo 7: cammini minimi, Dijkstra, Bellman-Ford e Floyd-Warshall.
• Articolo 8: questa pagina, programmazione dinamica e memoization.

L'elenco segue l'ordine di lettura e quell'ordine ha un senso: la misura della complessità (1) serve a giudicare tutto il resto; ordinamento e ricerca (2-3) e le strutture (4-5) forniscono i mattoni; le visite sui grafi (6) e i cammini minimi (7) li usano per problemi di rete; la programmazione dinamica (8) ripropone l'idea della tabella di Floyd-Warshall in problemi senza grafo. Se riparti da un punto solo, riparti dal primo: la notazione O è presupposta da ogni conteggio degli altri sette.

L'elenco completo degli otto articoli è nella pagina della serie.