Fahid Anas IS4

Rapport TP Algorithme Numérique pour l'Optimisation.

Professeur : François Lemaire.

TP1 /

Au sein de ce TP, nous devions dans un premier temps chercher à approximer une fonction de degré 3 passant par un ensemble de points expérimentaux.

Nous nous sommes appuyés sur la méthode des moindres carrées pour minimiser l'erreur (résidus) que l'on définit par la formule suivante :

Ri² = Yi² - f(xi)²

On commence par fixer des paramètres alpha, beta, gamma et mu et on affiche son graphe ainsi que le nuage des points expérimentaux sur l'intervalle [-1.2 ; 2.1].

import numpy as np from numpy import linalg as nla import matplotlib.pyplot as plt alpha, beta, gamma, mu = 0.5, -2, 1, 7 m = 4 def f(x): return alpha*(x**3) + beta*(x**2) + gamma*x + mu Tx = np.array ([-1.1, 0.17, 1.22, -.5, 2.02, 1.81]) p = Tx.shape[0] Ty_sur_la_courbe = np.array ([f(x) for x in Tx]) perturbations = 0.5*np.array ([-1.3, 2.7, -5, 0, 1.4, 6]) Ty_experimentaux = Ty_sur_la_courbe + perturbations x = np.linspace(-1.2,2.1) y = np.array([f(i) for i in x]) fig, ax = plt.subplots() ax.plot(x,y) plt.scatter(Tx, Ty_experimentaux) plt.show()

On peut alors calculer notre erreur initiale :

erreurInitiale = nla.norm(Ty_sur_la_courbe - Ty_experimentaux, 2) print(erreurInitiale)

On va chercher alors à minimiser cette erreur et avoir en sortie le vecteur x de taille m (où m est le nombre de paramètres) qui associé à cette minimisation.

Pour se faire on résout le système suivant AtAx= Atb, x est le vecteur contenant les 4 paramètres de notre fonction d'approximation et b l'ordonnées de nos points expérimentaux. On appellera x_optimal, le vecteur associé à la résolution de ce système précédemment énoncé.

A = np.empty([6, 4]) for i in range(len(Tx)): for j in range(0,4): A[i,j] = Tx[i]**(3-j) b = Ty_experimentaux AtA = np.dot(np.transpose(A), A) Atb = np.dot(np.transpose(A), b) alphaOpt, betaOpt, gammaOpt, muOpt = nla.solve(AtA, Atb) print(x)

On a alors pour x_optimal alpha = 1.58, beta = -3.14, gamma = -0.37 et gamma = 7.48.

On peut s’intéresser à présent à la nouvelle erreur quadratique entre les points expérimentaux et la nouvelle définition de f(x).

def f2(x): return alphaOpt*(x**3) + betaOpt*(x**2) + gammaOpt*x + muOpt x = np.linspace(-1.2,2.1) y = np.array([f2(i) for i in x]) Ty_sur_la_courbe = np.array ([f2(x) for x in Tx]) erreur2 = nla.norm(Ty_sur_la_courbe - Ty_experimentaux, 2) print(erreur2) ax2 = ax.twinx() ax.plot(x,y,'r') plt.show()

Remarquons, que l'on a considérablement diminuée la valeur de Ri passant de 4.24 à 3.37.

Dans un second temps, nous nous sommes attardés sur la méthode de Newton en une variable pour approximé la racine cubique de 2.

On sait que d'après la méthode de Newton nous fixons un U0 quelconque disons ici 5 et nous calculons les termes de la suite Un+1 = Un - f(Un)-f'(Un), où f est la fonction que l'on cherche à annuler. Ici on veut converger vers racine cubique de 2 on peut alors prendre f(x) = x³ - 2.

def f(x): return x**3 - 2 def fprime(x): return 3*(x**2) u = 5 for i in range(10): u = u - f(u)/fprime(u) print("U",i," ---- ",u) print(u)

La valeur théorique de la racine cubique de 2 étant approximativement égale à 1.260, on voit que dès le dixième terme de la suite on s'approche très fortement de cette valeur.

On veut maintenant effectuer le même traitement, cependant on s’appuiera ici sur les fonctions grad et éventuellement hessian.

import numdifftools as nd def f(x): return x**3 - 2 u = 5 for i in range(10): u = u - f(u)/nd.Gradient(f)([u]) print("U",i," ---- ",u) print(u)

A l'aide du package numdifftools on s'appuie sur la fonction Gradient qui nous permet d'être un substitut à la fonction fprime antérieurement utilisée.

TP 2 /

On cherche au sein de ce TP, à l'aide de la méthode de Newton, trouver les points stationnaires lié à F(a, b) qui est défini au sein du code ci-dessous. Nous commençons par une singulière estimation graphique... On cherche ici les intersections entres les plans de f(a, b) et g(a,

import autograd.numpy as np from autograd import grad, jacobian from mpl_toolkits.mplot3d import axes3d import matplotlib.pyplot as plt def f (a,b): return (a**3)*b - 3*(a**2)*(b-1) + b**2 - 1 def g (a,b): return (a**2)*(b**2) - 2 fig = plt.figure(figsize = (20,20)) ax = fig.add_subplot(1, 1, 1, projection='3d') ax.set_xlabel('$a$', labelpad=20) ax.set_ylabel('$b$', labelpad=20) ax.set_zlabel('$f(a,b)$ et $g(a,b)$', labelpad=20) xplot = np.arange (-8, 5, 0.1) yplot = np.arange (-8, 5, 0.1) ##### Graphe de f ##### X, Y = np.meshgrid (xplot, yplot) Z = f(X,Y) ax.plot_surface(X, Y, Z, cmap="spring_r", lw=0.5, rstride=1, cstride=1, alpha=0.5) ax.contour(X, Y, Z, 10, colors="green", linestyles="dashed") ax.contour(X, Y, Z, 1, colors="red", levels=np.array([0], dtype=np.float64), linestyles="solid") ##### Graphe de g ##### X, Y = np.meshgrid (xplot, yplot) Z = g(X,Y) ax.plot_surface(X, Y, Z, cmap="autumn_r", lw=0.5, rstride=1, cstride=1, alpha=0.5) ax.contour(X, Y, Z, 10, colors="k", linestyles="dashed") ax.contour(X, Y, Z, 1, colors="green", levels=np.array([0], dtype=np.float64), linestyles="solid")

On a alors représenté les fonctions f(a, b) et g(a, b) sur le plan 3D en colorant chacune des lignes de niveaux Y = 0 lié à chacune des fonctions et l'intersection de ces dernières correspondent aux coordonnées des points stationnaires solutions de l'équation F(a,b) = 0. A première vue, il semble y avoir 2 points qui sont solution cependant il peut en avoir plus !

On peut à présent s'appuyer sur la Jacobienne de la F(a,b) qui va nous permettre d'utiliser le binôme de Newton et trouver les coordonnées des deux zéros.

En effet on va chercher à résoudre le système suivant :

-Jac_F(a,b) * h = F(a,b), où h représente le vecteur de dimension 2 associé au décalage des coordonnées (x,y) pour se rapprocher d'un zéro potentiel.

from numpy import linalg as nla def Jac_F (a,b): Jacob = np.empty([2,2]) Jacob[0,0] = 3*(a**2)*b - 6*a*(b-1) Jacob[0,1] = a**3 - 3*a**2 + 2*b Jacob[1,0] = 2*a*(b**2) Jacob[1,1] = 2*b*(a**2) return Jacob def F(a,b): vec = np.empty([2,1]) vec[0,0] = f(a,b) vec[1,0] = g(a,b) return vec a = 2 b = 3 for i in range(300): Jacob=Jac_F(a,b) res = F(a,b) h = nla.solve(-Jacob, res) a = a + h[0][0] b = b + h[1][0] if i in range(295,300): print("X --> ",a, " Y --> ", b) print(F(a,b), "\n")

On obtient un zéro qui a pour coordonnées (-0.774, 1.826), bien sûr celui-ci n'est que un des deux zéros. Pour trouver le second on s'appuie sur les estimations faites auparavant qui joueront le rôle de notre U0 ici. Cependant n'ayant plus les valeurs expérimentales je suis contraint de passer à la suite prématurément.

On décide dans la suite du TP de passer par la bibliothèque autograd et utiliser la fonction jacobian et effectuer le même traitement que précédemment.

import torch import torch.autograd.functional as ag import autograd.numpy as np def F(u): a = u[0] b = u[1] vec = torch.tensor(np.empty([2,1])) vec[0,0] = f(a,b) vec[1,0] = g(a,b) return vec u = torch.tensor(np.array([-1.2, 0.3])) for i in range(300): res = F(u) Jacob = ag.jacobian(F, u) Jacob = np.concatenate((Jacob[0], Jacob[1])) h = nla.solve(-Jacob, res) u[0] = u[0] + h[0][0] u[1] = u[1] + h[1][0] if i in range(295,300): print("X --> ",u[0], " Y --> ", u[1]) print(F(u), "\n")

Pour finir, nous mettons en oeuvre la dérivation automatique qui nous permettra d'avoir les éléments constituants la Jacobienne de F.

def F(u) : a = u[0] b = u[1] # Pour f(a,b) t1 = a ** 2 t2 = b ** 2 t3 = (a*b - 3*b + 3)*t1 + t2 - 1 # Pour g(a,b) t1 = a ** 2 t2 = b ** 2 t4 = t1 * t2 - 2 # Résultat return np.array ([t3, t4]) def Jac_F(u): a = u[0] b = u[1] vec = np.empty([2,2]) # Pour f(a,b) t1 = a ** 2 t2 = b ** 2 t3 = (a*b - 3*b + 3)*t1 t5 = t3 + t2 - 1 # Pour g(a,b) t1 = a ** 2 t2 = b ** 2 t4 = t1 * t2 - 2 # Pour Jac_f(a,b) dt1_da = 2*a dt2_db = 2*b dt3_da = b*t1 + (t3/t1)*dt1_da dt3_db = (a-3)*t1 dt5_da = dt3_da + dt2_da dt5_db = dt3_db + dt2_db #Pour Jac_g(a,b) dt4_da = dt1_da + dt2_da dt4_db = dt1_db + dt2_db vec[0][0] = dt5_da vec[0][1] = dt5_db vec[1][0] = dt4_da vec[1][1] = dt4_db return vec

Voici ci-dessus l'alternative à la computation de la jacobienne via la méthode du forward mode ou dérivation automatique.

TP 3 /

Au sein de ce TP, nous devions également rechercher les points stationnaires du fonction polynomiale de degré 3. Nous nous sommes appuyé ici sur la matrice Hessienne ainsi que le nabla lié à la fonction f.

On peut commencer par une représentation graphique de notre fonction afin de pouvoir estimer graphiquement ces points stationnaires si il y en a.

import numpy as np from mpl_toolkits.mplot3d import axes3d import matplotlib.pyplot as plt def f(a,b): return a**3 + 2*(a**2) - 2*a*b + b**2 + a*(b**3) - 2*b + 5 fig = plt.figure(figsize = (20,20)) ax = plt.axes(projection='3d') ax.set_xlabel('$a$', labelpad=20) ax.set_ylabel('$b$', labelpad=20) ax.set_zlabel('$f(a,b)$', labelpad=20) xplot = np.arange (-1, 7, 0.1) yplot = np.arange (-1, 7, 0.1) X, Y = np.meshgrid (xplot, yplot) Z = f(X,Y) zmin = 3 zmax = 14 m, n = Z.shape for i in range (m) : for j in range (n) : Z[i,j] = max(min(Z[i,j], zmax), zmin) ax.plot_surface(X, Y, Z, cmap="spring_r", lw=0.5, rstride=1, cstride=1, alpha=0.5) ax.contour(X, Y, Z, 30, colors="k", linestyles="dotted") plt.show ()

Graphiquement, on peut apercevoir deux points stationnaires qui annulent f(a,b).

Écrivons à présent les fonctions nabla_f et H_f qui renvoient le gradient de f(a,b) et la matrice Hessienne respectivement.

def nabla_f(a, b): Grad = np.zeros((2,1), dtype=np.float64) Grad[0][0] = 3*(a**2) + 4*a -2*b + b**3 Grad[1][0] = -2*a + 2*b + 3*a*(b**2) - 2 return Grad def H_f(a,b): Hess = np.zeros((2,2), dtype = np.float64) Hess[0][0] = 6*a + 4 Hess[0][1] = -2 + 3*(b**2) Hess[1][0] = -2 + 3*(b**2) Hess[1][1] = 2 + 6*a*b return Hess

Nous sommes désormais à même de pouvoir calculer les coordonnées des points stationnaires de f(a, b) en résolavant le système suivant :

-H_f(a, b).h = nabla_f(a, b)

Ui = np.array([0.2, 0.9]) h = np.linalg.solve(-H_f(-1.5,0.8), nabla_f(-1.5,0.8)) for i in range (0,300): Ui[0] = Ui[0] + h[0][0] Ui[1] = Ui[1] + h[1][0] h = np.linalg.solve(-H_f(Ui[0],Ui[1]), nabla_f(Ui[0],Ui[1])) print(Ui)

Il semble y avoir un point stationnaires de coordonnées (-0.2113 ; 1.5665); Vérifions ce résultat.

def f(a,b): return a**3 + 2*(a**2) - 2*a*b + b**2 + a*(b**3) - 2*b + 5 print(f(0.2255, 0.9318))

Le résultat ne semblent pas très pertinent...

Une question intéressant peut-être soulevée. Quels points parmi ces points stationnaires sont des minimas locaux. En effet il peut y avoir des points de col (ou de selles). L'intitulé suivant peut-être une belle définition.

On dit que a est un point de col si il est critique (i.e dfa = 0) et si f ne présent pas d'extremums locaux en a. Plus simplement, au point a, il existe des directions vers lesquelles on monte en altitude et d'autres vers lesquelles nous descendons.

Si a est un point critique et que en ce point la Hessienne est définit positive, alors ce points est un minimum local.

ut = np.transpose(Ui) utH = np.dot(ut, H_f(Ui[0], Ui[1])) utHu = np.dot(utH, Ui) print(utHu)

Le résultat n'ayant pas été très pertinent nous avons un résultat difficile à interpreter. L'erreur est surement à la résolution du système -H_f*h = nabla_f.

Nous allons essayer ici de reprendre le problème initial, cela-dit, nous nous appuierons sur les fonctions np.gradient et np.hessian sans oublier de définir notre fonction avec un seul paramètre u qui contiendra les deux paramètres a et b.

import torch.autograd.functional as ag import autograd.numpy as np import numdifftools as nd def f(u): a = u[0] b = u[1] return a**3 + 2*(a**2) - 2*a*b + b**2 + a*(b**3) - 2*b + 5 u = torch.tensor(np.array([0.2, 0.9])) grad = nd.Gradient(f)(u) print(grad) hess = ag.hessian(f, u) print(hess.ndim) h = np.linalg.solve(-hess, grad) for i in range(100): u = u + h grad = nd.Gradient(f)(u) hess = ag.hessian(f, u) h = np.linalg.solve(-hess, grad) print(u) print(f(u))

On obtient les mêmes coordonnées que précédemment, cependant cela ne semble pas correspondre à un point stationnaire pour la fonction f. Je suis incapable de savoir d'où cela peut provenir...

La partie sur la dérivation backward mode ne sera malheureusement pas traité au sein de ce rapport.

TP 4 /

Au sein de ce dernier TP nous allons en premier lieu chercher les points stationnaires de f(a, b) sous la contrainte suivante :

a² + b² <= 1/2

La fonction f est la même que celle du TP3 il y a donc peu d’intérêt à refaire la visualisation graphique.

On nous informe que f(a, b) a cinq points stationnaires au total, et l'un d'entre eux a pour coordonnées (0.2255, 0.9312) fortement proche de celui trouvé via méthode du Gradient ou encore via la matrice Hessienne.

Vérifions si la contrainte est respectée.

u = [0.2255, 0.9312] print((u[0]**2 + u[1]**2) <= 0.5)

La contrainte n'est pas respectée. Si on veut transformer cette contrainte en une contrainte d'égalité on pose alors la contrainte suivante :

a² + b² - 1/2 = 0

Dans la section de code ci-dessous, nous affichons quelques vecteurs Gradients (en rouge pour la contrainte et en bleu pour l'objectif). Nous avons également représenté d'un trait noir plein la contrainte d'égalité susmentionnée.

import numpy as np from mpl_toolkits.mplot3d import axes3d import matplotlib.pyplot as plt def f (a,b): return a**3 + 2*a**2 - 2*a*b + b**2 + a*b**3 - 2*b + 5 def nabla_f (a,b) : return np.array ([b**3 + 3*a**2 + 4*a - 2*b, 3*a*b**2 - 2*a + 2*b - 2], dtype=np.float64) fig = plt.figure(figsize = (20,20)) ax = plt.axes(projection='3d') ax.set_xlabel('$a$', labelpad=20) ax.set_ylabel('$b$', labelpad=20) ax.set_zlabel('$f(a,b)$', labelpad=20) xplot = np.arange (-.8, 1.2, 0.05) yplot = np.arange (-.8, 1.2, 0.05) # Tracé du graphe de f X, Y = np.meshgrid (xplot, yplot) Z = f(X,Y) m, n = Z.shape for i in range (m) : for j in range (n) : Z[i,j] = max(min(Z[i,j], 12), 3) ax.plot_surface(X, Y, Z, cmap="spring_r", lw=0.5, rstride=1, cstride=1, alpha=0.5) ax.contour(X, Y, Z, 30, colors="k", linestyles="dotted") # Le minimim local de f sans tenir compte de la contrainte zero = np.array ([0.2255014396, 0.9318083312]) ax.scatter ([zero[0]], [zero[1]], [f(zero[0],zero[1])], color='black') # Tracé de la contrainte a**2 + b**2 = r**2 (avec r**2 = 1/2) r = 1/np.sqrt(2) n = 100 angles = np.linspace (-np.pi, np.pi, n) cxplot = np.array ([r * np.cos(theta) for theta in angles]) cyplot = np.array ([r * np.sin(theta) for theta in angles]) czplot = np.array ([f (cxplot[i], cyplot[i]) for i in range (0,n)]) ax.plot (cxplot, cyplot, czplot, color='black') def c(a,b): return a**2 + b**2 - 1/2 def nabla_c (a, b) : return np.array ([2*a, 2*b]) for angle in np.linspace(0,9,9): a, b = r*np.cos(angle), r*np.sin(angle) grad_f = nabla_f (a, b) grad_f = (.25/np.linalg.norm(grad_f,2)) * grad_f ax.quiver (a, b, f(a,b), grad_f[0], grad_f[1], 0, color='blue') ax.text (a, b, f(a,b), 'un gradient') grad_c = nabla_c (a, b) grad_c = (.25/np.linalg.norm(grad_c,2)) * grad_c ax.quiver (a, b, f(a,b), grad_c[0], grad_c[1], 0, color='red') ax.text (a, b, f(a,b), 'un gradient') plt.show ()