import time, math import numpy as np import matplotlib.pyplot as plt ############################################################################### # 02-descente_gradient.py # @title: Fondamentaux - Apprentissage par descente de gradient # @project: Mes scripts de ML # @lang: fr # @authors: Philippe Roy # @copyright: Copyright (C) 2023 Philippe Roy # @license: GNU GPL ############################################################################### ### # Commandes NumPy : # - np.array : créer un tableau à partir d'une liste de listes # - np.c_ : concatène les colonnes des tableaux # - np.ones : créer un tableau de 1 # - np.linalg.inv : inversion de matrice # - .T : transposé de matrice # - .dot : produit de matrice ### ############################################################################### # Initialisation ############################################################################### # Init du temps t_debut = time.time() # Init des plots fig = plt.figure(figsize=(15, 5)) fig.suptitle("Descente de gradient") donnees_ax = fig.add_subplot(141) # Observations : x1 et cibles : y model_ax = fig.add_subplot(142) # Modèle : theta0, theta1 couts_ax = fig.add_subplot(143) # Coûts : RMSE, MSE, ... app_ax = fig.add_subplot(144) # Taux d'apprentissage : eta i_list=[] # Itération couts_2d=[] couts_delta=[] couts_mse=[] # MSE couts_rmse=[] # RMSE eta_list=[] # Taux d'apprentissage ############################################################################### # Observations ############################################################################### # Observations d'apprentisage m = 1000 # Nombre d'observations bg = 1 # Quantité du bruit gaussien x1 = 2*np.random.rand(m, 1) # Liste des observations x1 y = 4 + 3*x1 + bg * np.random.rand(m, 1) # Liste des cibles y X = np.c_[np.ones((m, 1)), x1] # Matrice des observations, avec x0=1 donnees_ax.plot(x1, y, 'b.', label="Observations") exact_solution = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y) # Equation normale # Nouvelles observations x1_new=np.array([[0], [2]]) X_new = np.c_[np.ones((2, 1)), x1_new] # Matrice des observations, avec x0=1 ############################################################################### # Phases d'apprentissage et d'inférence ############################################################################### # Phase d'apprentissage par descente de gradient # - theta : vecteur paramètres du modèle # - gradient : gradient du coût en fonction de theta # - eta : taux d'appentissage eta = 0.01 # Taux d'appentissage (valeur par défaut : 0.1, hyperparamètre) n = 5000 # Nombre d'itérations (valeur par défaut : 1000 , hyperparamètre) # Calcul du coût (Mean Square Error (MSE)) def mse(theta): return np.sum((np.dot(X, theta) - y)**2)/m def rmse(theta): return math.sqrt(np.sum((np.dot(X, theta) - y)**2)/m) # Initialisation aléatoire theta= np.random.randn(2,1) theta0=[theta[0]] theta1=[theta[1]] # Descente du gradient for i in range(n): i_list.append(i) # Calcul du gradient du pas gradients = 2/m * X.T.dot(X.dot(theta) - y) eta_list.append(eta) theta = theta - eta * gradients theta0.append(theta[0]) theta1.append(theta[1]) # Calcul de l'erreur avec la norme du vecteur 2D (Objectif -> Theta) dans le plan (theta0, theta1) couts_2d.append(math.sqrt((theta[0]-exact_solution[0])**2+(theta[1]-exact_solution[1])**2)) # Calcul du RMSE à la main : # FIXME : étrange, je n'ai pas le même résultat qu'avec 'couts_rmse.append(rmse(theta))' delta = 0 for j in range (m): delta= delta +((theta[0] + theta[1]*x1[j])-(exact_solution[0] + exact_solution[1]*x1[j]))**2 delta = math.sqrt(delta/m) couts_delta.append(delta) # Calcul du RMSE et du MSE par les matrices couts_mse.append(mse(theta)) couts_rmse.append(rmse(theta)) # Prédiction du pas # Phase d'inférence (dernier pas) y_predict=X_new.dot(theta) donnees_ax.plot(x1_new, y_predict, 'c-', linewidth=0.5) ############################################################################### # Résultats ############################################################################### # Cartes des coûts (cc) # FIXME : la carte des coûts ne correspond au RMSE calculé # resol_cc=0.000000001 # resol_cc=0.1 # theta0_cc=np.linspace(theta[0]-resol_cc, theta[0]+resol_cc, 100).reshape(-1, 1) # theta1_cc=np.linspace(theta[1]-resol_cc, theta[1]+resol_cc, 100).reshape(-1, 1) theta0_cc=np.linspace(-5, 7, 100).reshape(-1, 1) theta1_cc=np.linspace(-5, 7, 100).reshape(-1, 1) theta0_cc_mg, theta1_cc_mg = np.meshgrid(theta0_cc, theta1_cc) r, c = theta0_cc_mg.shape rmse_cc = np.array([[0 for _ in range(c)] for _ in range(r)]) for i in range(r): for j in range(c): rmse_cc[i,j] = round(float(rmse(np.array([theta0_cc_mg[i,j], theta1_cc_mg[i,j]]))),3) # Plot des données donnees_ax.set_title("Données") donnees_ax.plot(x1_new, y_predict, 'r-', label="Prédictions") donnees_ax.set_xlabel(r'$x_1$') donnees_ax.set_ylabel(r'$y$', rotation=0) donnees_ax.legend() # Plot des paramètres du modèle model_ax.set_title("Paramètres du modèle") # model_ax.set(xlim=(theta[0]-resol_cc, theta[0]+resol_cc), ylim=(theta[1]-resol_cc, theta[1]+resol_cc)) model_ax.plot(theta0, theta1, '.', ls=':', color='r', fillstyle='none', label="Chemin", markevery=10) model_ax.plot(exact_solution[0], exact_solution[1], "o", color='k', fillstyle='full', label="Equation normale") model_ax.set_xlabel(r'$\theta_0$') model_ax.set_ylabel(r'$\theta_1 $', rotation=0) # levels = np.arange(20, 200, 0.25) # cc=model_ax.contour(theta0_mg, theta1_mg, carte_cout, levels) # model_ax.clabel(cc, inline=True, fontsize=10) # FIXME : la carte des coûts ne correspond au RMSE calculé # for i in range (0, 5000, 10): # Texte sur les points # model_ax.text(theta0[i], theta1[i], str(i)+" - "+str(round(float(rmse(np.array([theta0[i],theta1[i]]))),3)), va='center', ha='center') # cc_img = model_ax.pcolor(theta0_cc_mg, theta1_cc_mg, rmse_cc) # fig.colorbar(cc_img, ax=model_ax) model_ax.legend() # Plot du coût couts_ax.set_title("Coûts") couts_ax.plot(i_list, couts_2d, '.', ls=':', color='c', fillstyle='none', label="Coûts vecteur 2D", markevery=10) couts_ax.plot(i_list, couts_delta, '.', ls=':', color='r', fillstyle='none', label="Coûts RMSE à la main", markevery=10) couts_ax.plot(i_list, couts_mse, '.', ls=':', color='b', fillstyle='none', label="Coûts MSE", markevery=10) couts_ax.plot(i_list, couts_rmse, '.', ls=':', color='g', fillstyle='none', label="Coûts RMSE", markevery=10) couts_ax.set_xlabel(r'$i$') couts_ax.set_ylabel("Coûts") couts_ax.legend() # Plot du taux d'apprentissage app_ax.set_title("Taux d'apprentissage") app_ax.plot(i_list, eta_list, '.', ls=':', color='b', fillstyle='none', label="Taux d'apprentissage", markevery=10) app_ax.set_xlabel(r'$i$') app_ax.set_ylabel(r'$\eta$', rotation=0) # app_ax.legend() plt.show() # Performances print ("Theta th : theta0 : "+str(4)+" ; theta1 : "+str(3)) print ("Theta : theta0 : "+str(round(float(theta[0]),3))+" ; theta1 : "+str(round(float(theta[1]),3))) print ("Erreurs : theta0 : "+str(round(float(theta[0]-4),3))+" ; theta1 : "+str(round(float(theta[1]-3),3))) print ("Coûts RMSE : "+str(round(float(rmse(theta)),3))) print ("Temps : "+str(time.time()-t_debut))