User

I. Introduction

En informatique, la programmation dynamique est une méthode algorithmique pour résoudre des problèmes d'optimisation. Le concept a été introduit au début des années 1950 par Richard Bellman.

La programmation dynamique consiste à résoudre un problème en le décomposant en sous-problèmes, puis à résoudre les sous-problèmes, des plus petits aux plus grands en stockant les résultats intermédiaires.

Pour mieux comprendre cette méthode algorithmique on partira d'un exemple simple de calcul du n^ème terme de la suite de Fibonacci, puis on cherchera à l'aide de la programmation dynamique à optimiser le calcul des différences divisées.

Enfin, on montrera comment évaluer un polynôme d'interpolation de Newton en mémorisant ses différences divisées successives.


Note importante :

Par la suite, pour faciliter leur comptage, les appels récursifs des fonctions seront représentés sur les lignes d'un triangle (partiel ou entier), ces lignes seront numérotées de 0 à n.

La position des appels situés sur une même ligne ne respectera pas forcément leur chronologie.



II. Principe de la méthode - Suite de Fibonacci

Prenons comme exemple le calcul du n^ème terme de la suite de Fibonacci.


II-A. Définition de la suite de Fibonacci

Elle est définie par la relation de récurrence :

F_i = F_i-1 + F_i-2

Avec comme premiers termes :

F₀ = 0
F₁ = 1


II-B. Fonction récursive - Méthode intuitive

Partant de la définition précédente, on peut alors écrire une fonction récursive en Python qui évalue le n^ème terme de cette suite :

Code Python :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
def fib_v1(n):    
    if n <= 1: # si 1er ou 2e terme : n=0 ou n=1
        return n # retourne n
    else: # sinon, appels récursifs basés sur la relation de récurrence
        return fib_v1(n-1) + fib_v1(n-2)

Représentons maintenant à l'aide d'un schéma les différents appels récursifs de fib(5) :

Note importante : pour mieux visualiser les choses on a regroupé les appels récursifs identiques sur chaque ligne.


On a donc besoin de réaliser 15 appels (1 appel initial et 14 appels récursifs) pour évaluer fib(5), les calculs étant effectués à la remontée des appels.

L'algorithme est en fait en temps exponentiel.

On constate qu'il y a :

• 3 appels de fib(0) ;
• 5 appels de fib(1) ;
• 3 appels pour fib(2) ;
• 2 appels pour fib(3).

On va donc chercher à mémoriser les valeurs de fib(i) pour éviter de les recalculer par la suite.


II-C. Programmation dynamique - Méthode ascendante

Dans la méthode ascendante, on va "de bas en haut" : il s'agit du sens contraire à celui utilisé par les méthodes récursives. On commence par calculer des solutions pour les sous-problèmes les plus petits, puis, de proche en proche, on calcule les solutions des problèmes en utilisant le principe d'optimalité et on mémorise les résultats dans un tableau.

On souhaite donc calculer le n^ème terme F_n de la suite de Fibonacci en mémorisant successivement les termes d'indice 0, 1, ..., n-1 dans une liste Python pour éviter de les recalculer :

f[0] = 0
f[1] = 1
f[2] = f[1] + f[0] = 1 + 0 = 1
f[3] = f[2] + f[1] = 1 + 1 = 2
...
f[n] = f[n-1] + f[n-2]


On obtient alors la fonction itérative :

Code Python :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
def fib_v2(n):
 
    if n==0: # si 1er terme : n=0
        return n # retourne n et sort
 
    # initialisation de la liste des termes avec des zéros : f = [0, 0, ..., 0]
    f=[0]*(n+1)
 
    # initialisation des 2 pemiers termes de la liste : f = [0, 1, ..., 0]
    # f[0]= 0 inutile car on a déjà : f[0] = 0
    f[1] = 1 
 
    for i in range(2, n+1): # parcours des indices : 2, ..., n
        f[i] = f[i-1] + f[i-2] # mémorisation du terme d'indice i : Fi = Fi-1 + Fi-2
 
    return f[n] # renvoie le terme d'indice n : Fn
    # return f # retourne la séquence de Fibonacci : [0, 1, 1, 2, ..., Fn]

Note : on peut également retourner la séquence de Fibonacci f au lieu du nombre f[n].

Représentons maintenant l'enchaînement des calculs avec mémorisation des résultats intermédiaires pour fib(5) :

On constate que le calcul de Fib(n) est cette fois en temps linéaire.


Note importante :

Dans notre exemple d'optimisation on utilise la méthode dite « ascendante », mais il existe également une autre méthode dite « descendante » nécessitant d'écrire un algorithme récursif.



III. Évaluation des différences divisées


III-A. Définition

Les différences divisées d'ordre 0, 1 et 2 sont définies comme suit :



On obtient ainsi le schéma récursif suivant :




Note importante :

Elles interviennent notamment dans le polynôme d'interpolation de Newton N(x) associé à k+1 points (x₀, y₀), …,(x_k, y_k) et défini par :

N(x) = [y₀] + [y₀, y₁](x - x₀) + [y₀, y₁, y₂](x - x₀)(x - x₁) + ... + [y₀, y₁, ..., y_k](x - x₀)(x - x₁)...(x - x_k-1)


III-B. Fonction récursive - Méthode intuitive

En utilisant la relation de récurrence :

[y_i, ..., y_i+n] = ([y_i+1, ..., y_i+n] - [y_i, ..., y_i+n-1])/(x_i+n - x_i)

avec :

[y_i]= y_i


On peut alors écrire la fonction récursive :

Code Python :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
def difference_divisee_v1(x , y, i, n):
    # x, y : listes des valeurs
    # i, n : indice mini et ordre de la différence divisée [yi, yi+1, ..., yi+n]
    if n==0: # si on est à l'ordre 0
        return y[i] # on renvoie yi : [yi] = yi
    else: # sinon, appels récursifs selon la relation de récurrence [yi, yi+1, ..., yi+n] = ([yi+1, yi+1, ..., yi+n] - [yi, yi+1, ..., yi+n-1]) / (xi+n - xi)
        return (difference_divisee_v1(x , y , i+1, n-1) - difference_divisee_v1(x , y , i, n-1))/(x[i+n] - x[i])

Représentons maintenant à l'aide d'un schéma les différents appels récursifs nécessaires pour évaluer [y₀, ..., y₃] :

On a donc besoin de réaliser 2³⁺¹ - 1 = 15 appels (1 appel initial et 14 appels récursifs) pour évaluer [y₀, ..., y₃], les calculs étants effectués à la remontée des appels.

Par conséquent, l'évaluation des différences divisées à l'ordre n nécessite 2ⁿ⁺¹ - 1 appels, et "seulement" 2ⁿ - 1 calculs car les valeurs des différences divisées d'ordre 0 (ligne n) sont déjà connues.

Donc l'algorithme est en temps exponentiel en n.

On constate également qu'il y a 3 appels de [y₁] et [y₂], et 2 pour [y₁, y₂].

On va donc chercher à mémoriser les valeurs de [y_i, ..., y_i+k] pour éviter de les recalculer par la suite.


III-C. Programmation dynamique - Méthode ascendante

On reprend la relation de récurrence :

[y_i, ..., y_i+n] = ([y_i+1, ..., y_i+n] - [y_i, ..., y_i+n-1])/(x_i+n - x_i)

avec :

[y_i]= y_i


Si maintenant on part de la base, on peut alors évaluer la différence divisée en remontant le triangle ligne par ligne et en mémorisant à chaque fois les résultats intermédiaires :

On a donc besoin de réaliser 3×4/2 = 6 calculs pour évaluer [y₀, ..., y₃].

Donc le calcul des différences divisées à l'ordre n nécessite cette fois seulement n(n+1)/2 opérations.

On obtient ainsi la fonction itérative :

Code Python :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
def difference_divisee_v2(x , y , j, n):
    # x, y : listes des valeurs
    # j, n : indice mini et ordre de la différence divisée [yj, yj+1, ..., yj+n]
    dy = y.copy() # on copie la liste y dans dy : différences divisées à l'ordre 0
 
    for k in range(1,n+1): # on parcourt les ordres suivants des différences divisées : ordre 1, 2, ..., n       
 
        for i in range(j, j + n-k+1): # on parcourt les indices de x et dy jusqu'à j+n-k         
            h = (x[i+k] - x[i]) # on détermine la différence entre xi+k et xi pour le calcul de la différence divisée à l'ordre k
            dyi = (dy[i+1] - dy[i])/h # on évalue la différence divisée à l'ordre k : [yi, ..., yi+k])
            dy[i] = dyi # on remplace l'élément d'indice i de la liste dy par la valeur de la différence divisée dyi        
 
    return dy[j] # on renvoie la différence divisée

Note : on a voulu également optimiser l'espace mémoire en utilisant une liste à une seule dimension pour enregistrer les résultats intermédiaires.



IV. Évaluation du polynôme d'interpolation de Newton

Nous avons vu précédemment que le polynôme d'interpolation de Newton N(x) associé à k+1 points (x₀, y₀), …,(x_k, y_k) est défini par :

N(x) = [y₀] + [y₀, y₁](x - x₀) + [y₀, y₁, y₂](x - x₀)(x - x₁) + ... + [y₀, y₁, ..., y_k](x - x₀)(x - x₁)...(x - x_k-1)


Pour éviter de recalculer à chaque fois les différences divisées successives [y₀], [y₀, y₁], [y₀, y₁, y₂], ..., on doit d'abord les mémoriser en remontant le triangle ligne par ligne :

On peut ainsi écrire une fonction Python qui génère la liste des différences divisées à partir des 2 listes de valeurs x et y :

Code Python :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
def differences_divisees(x,y):    
    n=len(x)-1 # nombre de valeurs de la série - 1, il correspond au degré du polynôme, exemple pour n=2 : a2x^2 + a1x + a0, où a2, a1 et a0 représentent les coefficients
    dy = y.copy() # on copie la liste y dans dy : différences divisées à l'ordre 0
 
    diff_divisees = [dy[0]] # on initialise la liste permettant d'enregistrer les différences divisées : [y0]
 
    for k in range(1,n+1): # on parcourt les ordres suivants des différences divisées : ordre 1, 2, ..., n       
 
        for i in range(n-k+1): # on parcourt les indices de x et dy jusqu'à n-k         
            h = (x[i+k] - x[i]) # on détermine la différence entre xi+k et xi pour le calcul de la différence divisée à l'ordre k
            dyi = (dy[i+1] - dy[i])/h # on évalue la différence divisée à l'ordre k : [yi, ..., yi+k])
 
            dy[i] = dyi # on remplace l'élément d'indice i de la liste dy par la valeur de la différence divisée dyi        
 
        diff_divisees =  diff_divisees + [dy[0]] # on ajoute à la liste la différence divisée d'ordre k.
 
    return diff_divisees # on renvoie la liste des différences divisées

Évaluons maintenant le polynôme d'interpolation en se basant sur sa formule générale.

La fonction evaluer_polynome_newton prend en arguments 2 listes x et y, représentant la série de points (x_i,y_i), et renvoie la valeur du polynôme pour valeur_x :

Code Python :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
def evaluer_polynome_newton(x, y, valeur_x):
 
    diff_divisees = differences_divisees(x,y) # on génère la liste des différences divisées : [[y0], [y0, y1], [y0, y1, y2], ...]
    k=len(x)-1 # indice maxi de la liste des x
 
    valeur_polynome_newton = y[0] # on initialise la variable avec y[0]
    produit_facteurs=1 # on initialise le produit de facteurs
 
    for i in range(k): # parcours des indices : 0, 1, ..., k-1
        produit_facteurs = produit_facteurs*(valeur_x - x[i]) # on évalue le produit de facteurs : (x-x0)(x-x1)...(x-xi)
        valeur_polynome_newton = valeur_polynome_newton + diff_divisees[i+1]*produit_facteurs # on évalue progressivement le polynôme de Newton
 
    return valeur_polynome_newton # renvoie la valeur du polynôme de Newton

V. Module complet de test

Il contient l'ensemble des fonctions décrites dans ce billet et le code permettant de comparer les versions récursives et optimisées :

Code Python :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
import time
 
def fib_v1(n):
 
    if n<=1: # si 1er ou 2e terme : n=0 ou n=1
        return n # retourne n
    else: # sinon, appels récursifs basés sur la relation de récurrence
        return fib_v1(n-1) + fib_v1(n-2)
 
 
def fib_v2(n):
 
    if n==0: # si 1er terme : n=0 
        return n # retourne n et sort
 
    # initialisation de la liste des termes avec des zéros : f = [0, 0, ..., 0]
    f=[0]*(n+1)
 
    # initialisation des 2 pemiers termes de la liste : f = [0, 1, ..., 0]
    # f[0]= 0 inutile car on a déjà : f[0] = 0
    f[1] = 1 
 
    for i in range(2, n+1): # parcours des indices : 2, ..., n
        f[i] = f[i-1] + f[i-2] # mémorisation du terme d'indice i : Fi = Fi-1 + Fi-2
 
    return f[n] # renvoie le terme d'indice n : Fn
    # return f # retourne la séquence de Fibonacci : [0, 1, 1, 2, ..., Fn]
 
 
def difference_divisee_v1(x , y, i, n):
    # x, y : listes des valeurs
    # i, n : indice mini et ordre de la différence divisée [yi, yi+1, ..., yi+n]
    if n==0: # si on est à l'ordre 0
        return y[i] # on renvoie yi : [yi] = yi
    else: # sinon, appels récursifs selon la relation de récurrence [yi, yi+1, ..., yi+n] = ([yi+1, yi+1, ..., yi+n] - [yi, yi+1, ..., yi+n-1]) / (xi+n - xi)
        return (difference_divisee_v1(x , y , i+1, n-1) - difference_divisee_v1(x , y , i, n-1))/(x[i+n] - x[i])
 
def difference_divisee_v2(x , y , j, n):
    # x, y : listes des valeurs
    # j, n : indice mini et ordre de la différence divisée [yj, yj+1, ..., yj+n]
    dy = y.copy() # on copie la liste y dans dy : différences divisées à l'ordre 0
 
    for k in range(1,n+1): # on parcourt les ordres suivants des différences divisées : ordre 1, 2, ..., n       
 
        for i in range(j, j + n-k+1): # on parcourt les indices de x et dy jusqu'à j+n-k         
            h = (x[i+k] - x[i]) # on détermine la différence entre xi+k et xi pour le calcul de la différence divisée à l'ordre k
            dyi = (dy[i+1] - dy[i])/h # on évalue la différence divisée à l'ordre k : [yi, ..., yi+k])
            dy[i] = dyi # on remplace l'élément d'indice i de la liste dy par la valeur de la différence divisée dyi        
 
    return dy[j] # on renvoie la différence divisée
 
 
def differences_divisees(x,y):    
    n=len(x)-1 # nombre de valeurs de la série - 1, il correspond au degré du polynôme, exemple pour n=2 : a2x^2 + a1x + a0, où a2, a1 et a0 représentent les coefficients
    dy = y.copy() # on copie la liste y dans dy : différences divisées à l'ordre 0
 
    diff_divisees = [dy[0]] # on initialise la liste permettant d'enregistrer les différences divisées : [y0]
 
    for k in range(1,n+1): # on parcourt les ordres suivants des différences divisées : ordre 1, 2, ..., n       
 
        for i in range(n-k+1): # on parcourt les indices de x et dy jusqu'à n-k         
            h = (x[i+k] - x[i]) # on détermine la différence entre xi+k et xi pour le calcul de la différence divisée à l'ordre k
            dyi = (dy[i+1] - dy[i])/h # on évalue la différence divisée à l'ordre k : [yi, ..., yi+k])
 
            dy[i] = dyi # on remplace l'élément d'indice i de la liste dy par la valeur de la différence divisée dyi        
 
        diff_divisees =  diff_divisees + [dy[0]] # on ajoute à la liste la différence divisée d'ordre k.
 
    return diff_divisees # on renvoie la liste des différences divisées
 
def evaluer_polynome_newton_v1(x, y, valeur_x):
    k=len(x)-1 # indice maxi de la liste des x
 
    valeur_polynome_newton = y[0] # on initialise la variable avec y[0]
    produit_facteurs =1 # on initialise le produit de facteurs
 
    for i in range(k): # parcours des indices : 0, 1, ..., k-1
        produit_facteurs = produit_facteurs*(valeur_x - x[i]) # on évalue le produit de facteurs : (x-x0)(x-x1)...(x-xi)
        valeur_polynome_newton = valeur_polynome_newton + difference_divisee_v1(x,y,0,i+1)*produit_facteurs # on évalue progressivement le polynôme de Newton
 
    return valeur_polynome_newton # renvoie la valeur du polynôme de Newton
 
def evaluer_polynome_newton_v2(x, y, valeur_x):
 
    diff_divisees = differences_divisees(x,y) # on génère la liste des différences divisées : [[y0], [y0, y1], [y0, y1, y2], ...]
    k=len(x)-1 # indice maxi de la liste des x
 
    valeur_polynome_newton =y[0] # on initialise la variable avec y[0]
    produit_facteurs=1 # on initialise le produit de facteurs
 
    for i in range(k): # parcours des indices : 0, 1, ..., k-1
        produit_facteurs = produit_facteurs*(valeur_x - x[i]) # on évalue le produit de facteurs : (x-x0)(x-x1)...(x-xi)
        valeur_polynome_newton = valeur_polynome_newton + diff_divisees[i+1]*produit_facteurs # on évalue progressivement le polynôme de Newton
 
    return valeur_polynome_newton # renvoie la valeur du polynôme de Newton
 
n=40 
 
print("Comparaison de la durée d'exécution des 2 fonctions fib_v1 et fib_v2") 
print("")
 
print("Calcul de fib(" + str(n) + ") -  version n°1 :")
 
start_time = time.time()
fib_n = fib_v1(n)
print("fibo_v1(" + str(n) + ") = " + str(fib_n))
print("Durée d'exécution : %s secondes" % (time.time() - start_time))
 
print("")
print("Calcul de fib(" + str(n) + ") -  version n°2 :")
 
start_time = time.time()
fib_n = fib_v2(n)
print("fibo_v2(" + str(n) + ") = " + str(fib_n))
print("Durée d'exécution : %s secondes" % (time.time() - start_time))

Le code affiche :

Comparaison de la durée d'exécution des 2 fonctions fib_v1 et fib_v2

Calcul de fib(40) - version n°1 :
fibo_v1(40) = 102334155
Durée d'exécution : 108.11277770996094 secondes

Calcul de fib(40) - version n°2 :
fibo_v2(40) = 102334155
Durée d'exécution : 0.01198887825012207 secondes



VI. Conclusion

En partant d'un exemple simple consistant à évaluer le n^ème terme de la suite de Fibonacci, nous avons mieux compris le principe de la programmation dynamique.

Ceci nous a ensuite permis d'optimiser le calcul des différences divisées et celui du polynôme d'interpolation de Newton.

Enfin, nous avons pu vérifier sur nos fonctions que cette méthode algorithmique apporte un net gain de temps d'exécution.


Sources :

https://fr.wikipedia.org/wiki/Programmation_dynamique
https://en.wikipedia.org/wiki/Dynamic_programming
https://fr.wikipedia.org/wiki/Suite_de_Fibonacci
https://fr.wikipedia.org/wiki/Diff%C..._divis%C3%A9es
https://fr.wikipedia.org/wiki/Triangle_de_Pascal