Dans les grandes lignes, la programmation dynamique est une méthode de programmation qui permet d'éviter de répéter certains calculs en stockant les résultats intermédiaires.

Les étapes pour résoudre un problème à l'aide de la programmation dynamique sont :

• décomposer le problème en sous-problèmes : des versions « plus petites » du même problème ;
• établir une relation de récurrence entre ces sous-problèmes : comment la solution d'un problème se calcule à partir des solutions de problèmes plus petits ;
• résoudre les sous-problèmes dans l'ordre (des plus petits aux plus grands) en utilisant la relation de récurrence, et en stockant le résultat de chaque sous-problème pour une potentielle utilisation future.

La programmation dynamique est utile quand deux conditions sont réunies :

• la solution optimale du problème se construit à partir des solutions optimales des sous-problèmes (on parle de sous-structure optimale) ;
• les mêmes sous-problèmes reviennent plusieurs fois dans le calcul. S'ils étaient tous différents, les stocker ne servirait à rien : ce serait simplement « diviser pour régner ».

I. Exemple "simpliste" : la suite de Fibonacci

La suite de Fibonacci est définie par :

Si on utilise une fonction récursive simple pour calculer une telle suite, on tombe vite sur un problème.

def fib(n):
    if n <= 1:
        return 1
    else:
        return fib(n - 1) + fib(n - 2)

from time import time

debut = time()
print(fib(25))
print("calculé en", round(time() - debut, 3), "secondes")

Pour comprendre pourquoi c'est si lent, dessinons l'arbre des appels récursifs pour . Les cases en orange sont des calculs déjà faits ailleurs dans l'arbre :

est calculé 2 fois, 3 fois, 5 fois... Il n'y a que 6 valeurs différentes à calculer (), mais l'arbre a 15 nœuds. Et chaque fois qu'on augmente de 1, l'arbre contient une copie de l'arbre précédent et de celui d'avant : sa taille est multipliée par environ .

Cette fonction va donc faire un nombre exponentiel d'appels récursifs (de l'ordre de , où est le nombre d'or) !

La version en programmation dynamique consiste à calculer non plus directement, mais toute la suite en partant de 0, en stockant les résultats dans un tableau (une liste). Chaque case se calcule à partir des deux cases précédentes, qui sont déjà remplies :

On remplit ce tableau de gauche à droite : , puis , puis ... Chaque valeur n'est calculée qu'une seule fois.

def fib(n):
    resultats = [1, 1]
    for i in range(2, n + 1):
        resultats.append(resultats[-1] + resultats[-2])
    return resultats[-1]

debut = time()
print(fib(25))
print("calculé en", round(time() - debut, 5), "secondes")
print(fib(200))

On peut comparer les temps de calcul des deux versions, en fonction de :

Cette fois, la complexité est en . C'est bien mieux !

II. Un vrai problème : le rendu de monnaie

a) Exemples

On prend un autre système monétaire, beaucoup moins pratique, avec des coupures de 1, 3 et 4, et on essaie de rendre 6. Cette fois, l'algorithme glouton renvoie 4, 1, 1. Ce n'est pas optimal ! On aurait pu utiliser 3 et 3.

b) Définition formelle et résolution du problème

Pourquoi cette relation ? Pour rendre avec les pièces , il n'y a que deux possibilités :

• soit on n'utilise pas la pièce : on rend avec les pièces seulement, ce qui coûte pièces au minimum ;
• soit on utilise au moins une fois la pièce : on la rend (1 pièce), et il reste à rendre, toujours avec (on a le droit de réutiliser ), ce qui coûte au minimum pièces.

On ne sait pas à l'avance laquelle des deux possibilités est la meilleure : on calcule les deux, et on garde la plus petite.

Voici le tableau complet des pour les pièces 1, 3, 4 :

Par exemple, pour la case en gras : . Pour la calculer, on a eu besoin de la case juste au-dessus () et d'une case plus à gauche sur la même ligne (). Il suffit donc de remplir le tableau dans un ordre où ces cases sont toujours calculées avant : par exemple colonne par colonne, de gauche à droite.

def rendu_monnaie(pieces, S):
    #initialisation du tableau résultat : resultat[i][s] = m_(i+1)(s)
    #(la colonne s = 0 reste à 0 : on rend 0 avec 0 pièce)
    resultat = [[0]*(S + 1) for _ in pieces]

    #avec seulement la pièce de 1, il faut s pièces pour rendre s
    for s in range(S + 1):
        resultat[0][s] = s

    #remplissage du tableau par s et i croissants
    for s in range(1, S + 1):
        for i in range(1, len(pieces)):
            #cas où la dernière pièce est trop grande
            if pieces[i] > s:
                resultat[i][s] = resultat[i-1][s]
            #cas général : on prend la dernière pièce ou non
            else:
                resultat[i][s] = min(1 + resultat[i][s - pieces[i]],
                                     resultat[i-1][s])
    return resultat[-1][-1]

print(rendu_monnaie([1, 3, 4], 6))

Pour voir l'ordre de remplissage du tableau, on remplace la liste de listes par un tableau suivi (du module affichage du site), qui s'utilise de la même manière, mais retient les cases lues et écrites :

Cette fonction a une complexité (la taille du tableau).

c) Retrouver la répartition

Comment peut-on alors retrouver la répartition des pièces ?

On peut suivre en remontant le chemin qu'a pris le résultat : si 1 + resultat[i][s - pieces[i]] < resultat[i-1][s], c'est qu'on a pris la pièce i, et inversement, jusqu'à arriver à un bord du tableau.

Sur le tableau ci-dessus, en partant de la case en bas à droite, :

1. est plus petit que : on n'a pas utilisé la pièce 4. On monte d'une ligne, en .
2. est plus petit que : on a utilisé la pièce 3. On la note, et on recule de 3 colonnes, en .
3. est plus petit que : on a encore utilisé la pièce 3. On arrive en : il ne reste rien à rendre.

On a rendu 3 + 3.

Pour cela, on a besoin de tout le tableau : on modifie la fonction précédente pour qu'elle renvoie resultat plutôt que resultat[-1][-1].

def tableau_rendu_monnaie(pieces, S):
    resultat = [[0]*(S + 1) for _ in pieces]
    for s in range(S + 1):
        resultat[0][s] = s
    for s in range(1, S + 1):
        for i in range(1, len(pieces)):
            if pieces[i] > s:
                resultat[i][s] = resultat[i-1][s]
            else:
                resultat[i][s] = min(1 + resultat[i][s - pieces[i]],
                                     resultat[i-1][s])
    return resultat

def retrouver_resultat(resultat, pieces, S):
    s = S
    i = len(pieces) - 1
    a_rendre = []
    while i > 0 and s > 0:
        if pieces[i] > s or resultat[i-1][s] < 1 + resultat[i][s - pieces[i]]:
            i = i-1
        else:
            s = s - pieces[i]
            a_rendre.append(pieces[i])
    #il ne reste que des pièces de 1
    return a_rendre + [pieces[0]] * s

resultat = tableau_rendu_monnaie([1, 3, 4], 6)
print("i/s", list(range(7)))
for (i, l) in enumerate(resultat):
    print(i, "  ", l)
print(retrouver_resultat(resultat, [1, 3, 4], 6))

III. Un autre problème : la distance d'édition

Lorsqu'on veut mesurer la différence entre deux chaînes de caractères, il n'est pas toujours adapté de simplement comparer les lettres une à une : par exemple, "bonjour" et "onjour" sont des chaînes proches, mais n'ont aucune lettre en commun (à la même position).

La distance d'édition est une solution à ce problème.

Par exemple, pour passer de "bnjvr" à "bonjour", on peut faire 3 opérations : on aligne les deux mots, en laissant des trous (« – ») là où il faut ajouter une lettre.

On ne peut pas faire mieux : la distance d'édition entre "bnjvr" et "bonjour" est 3. C'est ce qu'utilisent les correcteurs orthographiques pour proposer les mots les plus proches d'un mot mal écrit.

Pour trouver les sous-problèmes, on regarde ce qui se passe à la fin des deux mots. Pour transformer en , la dernière lettre a été obtenue :

• soit en supprimant la dernière lettre de , puis en transformant en ;
• soit en ajoutant la lettre à la fin, après avoir transformé en ;
• soit à partir de la dernière lettre : on transforme en , puis on modifie en (opération gratuite si ce sont les mêmes lettres).

Le cas de base signifie qu'à partir du mot vide, il faut ajouter les lettres ; , qu'il faut supprimer les lettres pour arriver au mot vide.

Dans le tableau des , chaque case se calcule à partir de trois cases voisines : celle au-dessus (suppression), celle à gauche (ajout) et celle en diagonale (modification).

On peut donc remplir le tableau ligne par ligne, de gauche à droite : les trois cases nécessaires sont toujours déjà remplies.

def distance_edition(s1, s2):
    """calcule la distance d'édition entre deux chaines de caractères s1 et s2, par la méthode de la programmation dynamique."""
    #initialisation du tableau résultat
    resultat = [[0]*(len(s2) + 1) for _ in range(len(s1) + 1)]
    for k in range(len(s2) + 1):
        resultat[0][k] = k
    for k in range(len(s1) + 1):
        resultat[k][0] = k

    for i in range(1, len(s1) + 1):
        for j in range(1, len(s2) + 1):
            #calcul des trois cas
            supp = resultat[i-1][j] + 1
            ajout = resultat[i][j-1] + 1

            if s1[i-1] != s2[j-1]:
                modif = resultat[i-1][j-1] + 1
            else:
                modif = resultat[i-1][j-1]
            resultat[i][j] = min(supp, ajout, modif)

    return resultat[-1][-1]

print(distance_edition("bnjvr", "bonjour"))

Le tableau a cases, chacune calculée en temps constant : la complexité est en .

Le remplissage du tableau, case par case (les cases orange sont celles qui servent à calculer la case rouge) :

IV. Dernier problème : l'algorithme de Floyd-Warshall

Petit retour sur nos amis les graphes ! On va utiliser la programmation dynamique pour calculer le plus court chemin entre deux sommets d'un graphe.

On va prendre des contraintes un peu différentes par rapport à Dijkstra : on cherche le plus court chemin entre chaque paire de sommets, pour tous les sommets du graphe.

On se donne un graphe orienté pondéré avec un ensemble de sommets et un ensemble d'arêtes (les arêtes peuvent avoir un poids négatif !). On suppose que n'a pas de cycle de poids négatif.

Soit la fonction qui renvoie :

• si ;
• si et sont reliés par une arête de poids ;
• sinon.

Cela correspond à la matrice d'adjacence.

L'idée : on autorise les sommets intermédiaires un par un. Au départ (), on n'a droit à aucun sommet intermédiaire : seules les arêtes directes comptent. Puis on autorise le sommet 1, puis les sommets 1 et 2, etc. Quand on autorise le sommet , un plus court chemin de à :

• soit ne passe pas par : rien ne change par rapport à l'étape précédente ;
• soit passe par (une seule fois, puisqu'il n'y a pas de cycle de poids négatif) : il est formé d'un plus court chemin de à , puis d'un plus court chemin de à , qui n'utilisent tous les deux que des sommets de .

On compare le chemin en pointillés (sans ) et le chemin en bleu (par ), et on garde le plus court.

Par exemple, sur le graphe suivant (dont les sommets sont numérotés à partir de 0) :

from math import inf
import pprint

graphe = {
    0 : [(3, -3), (2, 7)],
    1 : [(0, 11), (3, 1)],
    2 : [(1, 5), (4, -1)],
    3 : [(2, 2)],
    4 : [(0, 2), (3, 5)],
}

def floyd_warshall(G):
    #construction de la matrice d'adjacence
    mat_adj = [[inf] * len(G) for _ in G]
    for (s, voisins) in G.items():
        for (v, w) in voisins:
            mat_adj[s][v] = w
        mat_adj[s][s] = 0

    #les chemins en passant par aucun sommet sont donnés par la matrice d'adjacence
    distances = [mat_adj]
    #k : chemins en passant par les sommets 0, ..., k
    for k in range(len(G)):
        new_dist = [[None] * len(G) for _ in G]
        for i in range(len(G)):
            for j in range(len(G)):
                new_dist[i][j] = min(distances[-1][i][j], distances[-1][i][k] + distances[-1][k][j])
        distances.append(new_dist)

    return distances

distances = floyd_warshall(graphe)
pprint.pp(distances[-1])

Il y a étapes (une par valeur de ), et à chaque étape on calcule cases en temps constant : la complexité est en . Pour comparer, appliquer l'algorithme de Dijkstra depuis chacun des sommets coûterait , et Dijkstra ne fonctionne pas avec des poids négatifs.

Que se passe-t-il s'il y a un cycle de poids négatif ?

graphe2 = {
    0 : [(1, -1)],
    1 : [(0, -1)]
}
pprint.pp(floyd_warshall(graphe2))

Avec un cycle de poids négatif, on peut faire baisser le poids d'un chemin autant qu'on veut en tournant dans le cycle : les « plus courts chemins » n'existent plus. On le repère dans le résultat : une case de la diagonale devient négative ( signifie qu'on peut revenir en avec un poids négatif). C'est un moyen simple de détecter les cycles de poids négatif.

En fait, on peut même modifier le tableau sur place (les cases et ne changent pas pendant l'étape ). Voici le remplissage de ce tableau unique, étape par étape :

Reconstruction des chemins : on peut retrouver un plus court chemin en regardant, pour chaque , si le chemin passe par ou non.

def chemin(distances, x, y, k):
    """Renvoie un chemin le plus court entre x et y en utilisant les sommets jusqu'à k.
    Attention, ne termine pas si on a un cycle de poids négatif."""
    if (k == -1):
        #le chemin est direct entre x et y
        return []
    #cas où on n'a pas besoin du sommet k
    if distances[k+1][x][y] == distances[k][x][y]:
        return chemin(distances, x, y, k - 1)
    else:
        return chemin(distances, x, k, k - 1) + [k] + chemin(distances, k, y, k - 1)

print(chemin(distances, 1, 0, 4))

Pour reconstruire les chemins, on peut aussi stocker en plus un tableau des prédécesseurs (voir ci-dessous).

V. Retrouver le chemin avec le tableau des prédécesseurs

C'est utile dans le cas où le calcul de la relation de récurrence n'est pas rapide (par exemple, un maximum sur valeurs).

L'idée est de créer un tableau de même dimension que le tableau de stockage du résultat, qui, pour chaque case, indique quelle case a servi à la calculer. Par exemple, pour la distance d'édition :

def distance_edition(s1, s2):
    """calcule la distance d'édition entre deux chaines de caractère s1 et s2, par la méthode de la programmation dynamique."""
    #initialisation du tableau résultat
    resultat = [[0]*(len(s2) + 1) for _ in range(len(s1) + 1)]
    predecesseur = [["_"]*(len(s2) + 1) for _ in range(len(s1) + 1)]
    for k in range(len(s2) + 1):
        resultat[0][k] = k
    for k in range(len(s1) + 1):
        resultat[k][0] = k

    for i in range(1, len(s1) + 1):
        for j in range(1, len(s2) + 1):
            #calcul des trois cas
            supp = resultat[i-1][j] + 1
            ajout = resultat[i][j-1] + 1

            if s1[i-1] != s2[j-1]:
                modif = resultat[i-1][j-1] + 1
            else:
                modif = resultat[i-1][j-1]
            resultat[i][j] = min(supp, ajout, modif)
            if resultat[i][j] == supp:
                predecesseur[i][j] = "↑"
            elif resultat[i][j] == ajout:
                predecesseur[i][j] = "←"
            else:
                predecesseur[i][j] = "↖"

    return resultat, predecesseur

resultat, predecesseur = distance_edition("bnjvr", "bonjour")
for l in resultat:
    print(l)
print()
for l in predecesseur:
    print(l)

Pour retrouver les opérations, on part de la case en bas à droite et on suit les flèches jusqu'à la case en haut à gauche. Voici le chemin obtenu (surligné), dans le tableau des prédécesseurs :

En lisant le chemin dans l'ordre (de haut en bas), chaque flèche est une opération :

• ↖ sur deux lettres égales : on garde la lettre (b, n, j, r) ;
• ↖ sur deux lettres différentes : modification (v → o) ;
• ← : ajout d'une lettre du deuxième mot (o, puis u) ;
• ↑ : suppression d'une lettre du premier mot (il n'y en a pas ici).

On trouve « b, ajout o, n, j, v → o, ajout u, r » : encore 3 opérations, mais pas les mêmes que dans l'exemple du III. Il peut y avoir plusieurs suites d'opérations optimales ; le tableau des prédécesseurs en donne une.

VI. Utiliser des dictionnaires pour la mémoïsation

Jusqu'ici, on a toujours rempli un tableau des plus petits sous-problèmes vers les plus grands : c'est la méthode ascendante (bottom-up). Il faut alors réfléchir à l'ordre de remplissage.

On peut aussi partir de la fonction récursive « naïve », qui découpe le problème du plus grand vers les plus petits (méthode descendante, top-down), et lui ajouter une mémoire : avant de calculer, on regarde si le résultat a déjà été calculé ; après avoir calculé, on le note. C'est la mémoïsation. Un dictionnaire est pratique pour cette mémoire : les clés sont les arguments, les valeurs sont les résultats.

Sur Fibonacci, l'arbre des appels devient beaucoup plus petit. Les appels en pointillés trouvent leur résultat dans le dictionnaire et s'arrêtent immédiatement :

Le nombre d'appels devient linéaire en , comme avec le tableau.

Voici le rendu de monnaie, avec un dictionnaire dont les clés sont les couples (i, S) :

def rendu_monnaie(pieces, S, i, dico):
    """Nombre minimal de pièces parmi pieces[0], ..., pieces[i] pour rendre S.
    dico contient les résultats déjà calculés. Premier appel : rendu_monnaie(pieces, S, len(pieces) - 1, {})"""
    if i == 0:
        return S
    if (i, S) in dico:
        return dico[(i, S)]
    #cas où la dernière pièce est trop grande
    if pieces[i] > S:
        dico[(i, S)] = rendu_monnaie(pieces, S, i - 1, dico)
    #cas général : on prend la dernière pièce ou non
    else:
        dico[(i, S)] = min(1 + rendu_monnaie(pieces, S - pieces[i], i, dico),
                           rendu_monnaie(pieces, S, i - 1, dico))
    return dico[(i, S)]

print(rendu_monnaie([1, 3, 4], 6, 2, {}))

Les deux méthodes ont la même complexité dans le pire cas. Laquelle choisir ?

Récapitulatif : la méthode

Exercices

def nb_chemins(n, m):
    c = [[1] * (m + 1) for _ in range(n + 1)]
    for i in range(1, n + 1):
        for j in range(1, m + 1):
            c[i][j] = c[i-1][j] + c[i][j-1]
    return c[n][m]

nb_chemins(3, 3)   # 20

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

lcs("programmation", "dynamique")   # 3 ("ami")

def sac_a_dos(poids, valeurs, P):
    n = len(poids)
    v = [[0] * (P + 1) for _ in range(n + 1)]
    for i in range(1, n + 1):
        for p in range(P + 1):
            v[i][p] = v[i-1][p]
            if poids[i-1] <= p:
                v[i][p] = max(v[i][p], v[i-1][p - poids[i-1]] + valeurs[i-1])
    return v[n][P]

sac_a_dos([2, 3, 4, 5], [3, 4, 5, 6], 5)   # 7 : les objets de poids 2 et 3