dans cette article, nous allons étudier et implémenter la Régression Linéaire en Python avec Numpy. A vos ceintures
Bonjour à toutes et à tous,
Dans ce nouvel article, nous allons étudier et implémenter un nouvel algorithme de base, à savoir la régression linéaire. Pour ce premier article nous nous limiterons à la régression sur une seule variable, nous verrons la régression sur plusieurs variables dans un prochain article.
Tout d'abord, "remercier" Andrew Ng pour sa série de vidéos sur la régression linéaire, même si en vrai il s'en fiche certainement. Au mieux il doit être content d'inspirer d'autres personnes à étudier le machine learning et se lancer à faire eux-même leur cours et tutos. Allez donc jeter un coup d'oeil, il y a bien plus de matériel et mieux expliqué que je l'ai fais, sur la régression linéaire.
Sommaire¶
Introduction et Notation¶
En quoi consiste la régression linéaire ? Soit un ensemble de points suivants
Il s'agit de trouver la droite d'équation $y = ax + b$
ou dans l'exemple précédent, $poids = a*taille + b$, telle que cette droite soit la mieux ajustée aux données.
Pourquoi faire cela ? Mais quelle question, parce-que :p !! Non plus sérieusement:
- pour effectuer de la prédiction: prenons l'exemple précédent, à partir de ces mesures de taille et de poids, on est en mesure d'estimer le poids d'un individu dont on ne connaît que la taille, et inversement, la taille d'un individu dont on ne connaît que le poids.
- pour étudier des liaisons entre variables: plus de détails et d'exemples dans un prochain article.
La régression linéaire est un algorithme d'apprentissage automatique ou "machine learning", et, dans l'ensemble des algorithmes de machine learning, c'est, comme son nom l'indique, un algorithme de régression, par opposition aux algorithmes de classification, et c'est de plus un algorithme supervisé, par opposition aux algorithmes "non-supervisés"
La sortie de la classification est un ensemble prédéfini, souvent restreint de valeurs:
Il n'y a pas forcément d'ordre entre les différentes valeurs de la sortie, et celle-ci est même souvent binaire.
apprentissage supervisé vs apprentissage non-supervisé¶
Pour ne pas trop encombrer cet article, nous verrons cela plus en détails dans un autre article. Sinon pour aller au plus simple, dans l'apprentissage supervisé, pendant l'apprentissage, pour chaque exemple, on a la classe d'appartenance de l'exemple.
Dans l'apprentissage non-supervisé, on ne dispose pas de la classe d'appartenance de l'exemple.
Termes et Notation¶
- x: variables d'entrée
- y: variable de sortie
- (x, y): un exemple d'apprentissage
- $(x^{(i)}, y^{(i)})$: le i-ème exemple d'apprentissage
- ensemble d'apprentissage: l'ensemble des exemples d'apprentissage
- m: nombre d'exemple d'apprentissage.
Ainsi pour l'exemple suivant:
import matplotlib.pyplot as plt
x=range(5)
y=[1.5,2.,2.5,3.,3.5]
plt.ylim(0,5)
plt.xlim(0,5)
plt.plot(x,y,'bx')
plt.show()
- $x = [0, 1, 2, 3, 4]$
- $y = [1.5, 2.0, 2.5, 3.0, 3.5]$
- un exemple d'apprentissage est par exemple $(2, 2.5)$
- $(x^{(0)},y^{(0)}) = (0, 1.5)$, $(x^{(1)},y^{(1)}) = (1, 2.0)$, etc.
- l'ensemble d'apprentissage est : $[(0, 1.5), (1, 2.0), (2, 2.5), (3, 3.0), (4, 3.5)]$
- $m=5$, il y a 5 exemples d'apprentissage
Reprenons la définition des termes et de la notation
- modèle ou "hypothèse": $h_{\theta}(x) = \theta_0 + \theta_1 x$
- paramètres du modèles: vecteur $\theta = [\theta_0, \theta_1]$
Fonction de coût¶
En introduction, nous avons défini la régression linéaire comme étant un algorithme permettant de déterminer la droite la mieux ajustée aux données ... la mieux ajustée, ça veut dire quoi ? Comment le définit-on ?
Reprenons le tout premier exemple :
Quelle est la droite la mieux ajustée aux données ? Est-ce la droite D1 ? la droite D2 ? la droite D3 ? Une autre droite peut-être ?
Comment peut-on mesurer l'ajustement d'une droite aux données ?
C'est là qu'entre la fonction de coût, on va établir une fonction qui va mesurer le coût d'une droite par rapport aux données. Intuitivement le coût est l'inverse de l'ajusement, minimiser l'un revient à maximiser l'autre.
Etablissons cette fonction de coût petit à petit et essayons de le faire intuitivement. Comment pourrait-on mesurer le coût ou l'ajustement d'une droite par rapport aux données ?
Nous pourrions mesurer les écarts sur la variable de sortie entre la droite et les données. Par exemple:
Le problème est que dans cet exemple, $e_1 = -e_5$, $e_2 = -e_4$ et $e_3 = 0$, ainsi, la somme des écarts = 0, alors que la droite n'est absolument pas ajustée aux données.
Pour remédier à cela, on met les écarts au carré.
Nous obtenons au final la fonction de coût suivante : $$ J(\theta) = \frac{1}{2m} \sum_{i=1}^m \left( h_{\theta}(x^{(i)}) - y^{(i)} \right)^2 $$
$h_{\theta}(x^{(i)})$ est la valeur obtenu pour l'exemple $x^{(i)}$ par la droite d'équation $h_{\theta}(x) = \theta_0 + \theta_1 x^{(i)}$
On multiplie par $\frac{1}{2m}$, qui est une constante, afin de simplifier des calculs ultérieurs.
Minimisation de la fonction de coût¶
On arrive à une nouvelle formulation du but de la régression linéaire: $$ \underset{\theta = [\theta_0, \theta_1]}{\operatorname{min}} J(\theta) = \frac{1}{2m} \sum_{i=1}^m \left( h_{\theta}(x^{(i)}) - y^{(i)} \right)^2 $$
où l'on rappelle que :
- $h_{\theta}(x) = \theta_0 + \theta_1 x^{(i)}$ est l'hypothèse ou le modèle utilisé
- $\theta = [\theta_0, \theta_1]$ est le vecteur des paramètres de l'hypothèse ou modèle
- m est le nombre d'exemples d'apprentissage
- $y^{(i)}$ est la "vraie valeur" associée à l'exemple $x^{(i)}$
Descente de Gradient¶
Pour minimiser la fonction de coût définie juste au-dessus, nous allons utiliser l'algorithme de la descente de gradient. Pour l'illustrer, nous allons voler ... pardon, emprunter cette illustration récupéré après une recherche Google ici
!
dans cette illustration, "loss" désigne la fonction de coût, elle est en effet souvent nommée "loss" en anglais. Et w désigne $\theta$
$\alpha$ désigne le taux d'apprentissage, typiquement $\alpha \in [0,1]$, et très souvent, $\alpha = 0.05$
Lé dérivée nous indique de combien augmente la fonction de coût (loss) quand on augmente le paramètre $\theta$ (w) d'une quantité infinitésimale (Cette quantité est négative quand la fonction diminue).
Puisque nous cherchons à minimiser le coût, et non à le maximiser, à chaque étape d'apprentissage, nous allons mettre à jour $\theta$ avec une fraction de la dérivée du coût.
Appliqué à notre cas, on obtient:
$$ \theta_{t+1} = \theta_{t} - \alpha \frac {\partial J(\theta)} {\partial \theta_{t}} $$Calcul du gradient¶
On rappelle l'expression de la fonction de coût:
$$ J(\theta) = \frac{1}{2m} \sum_{i=1}^m \left( h_{\theta}(x^{(i)}) - y^{(i)} \right)^2 $$calculons maintenat son gradient $\frac {\nabla J(\theta)} {\nabla\theta}$:
$$ \frac {\nabla J(\theta)} {\nabla\theta} = \begin{bmatrix} \frac {\partial J(\theta)} {\partial\theta_0} \\ \frac {\partial J(\theta)} {\partial\theta_1} \end{bmatrix} = \begin{bmatrix} \frac{1}{m} \sum_{i=1}^m \left( h_{\theta}(x^{(i)}) - y^{(i)} \right) \\ \frac{1}{m} \sum_{i=1}^m \left( h_{\theta}(x^{(i)}) - y^{(i)} \right)x^{(i)} \end{bmatrix}\\ $$Le gradient du coût est le vecteur dont chaque composante est la dérivée en fonction de la composante.
Voici un lien math.stackexchange pour un refresher rapide sur le gradient, la dérivée, et la différence entre les 2. Bon, vous l'auriez trouvé tout seul avec une recherche Google mais bon :p .
Voici les détails du calcul de $\frac {\partial J(\theta)} {\partial\theta_0}$, on va y aller étape par étape :
$$ \begin{align} \frac {\partial J(\theta)} {\partial\theta_0} & = \frac {\partial } {\partial\theta_0} \left ( \frac{1}{2m} \sum_{i=1}^m \left( h_{\theta}(x^{(i)}) - y^{(i)} \right)^2 \right) & \text{par définition de }J(\theta) \\ & = \frac{1}{2m} \frac {\partial } {\partial\theta_0} \left (\sum_{i=1}^m \left( h_{\theta}(x^{(i)}) - y^{(i)} \right)^2 \right) & \text{on extrait la constante} \\ & = \frac{1}{2m} \sum_{i=1}^m \frac {\partial } {\partial\theta_0} \left (\left( h_{\theta}(x^{(i)}) - y^{(i)} \right)^2 \right) & \text{car dérivée de la somme = somme des dérivées} \\ & = \frac{1}{2m} \sum_{i=1}^m \frac {\partial } {\partial\theta_0} \left (\left( \theta_0 + \theta_1 x^{(i)} - y^{(i)} \right)^2 \right) & \text{par définition de }h_{\theta} \\ & = \frac{1}{2m} \sum_{i=1}^m \frac {\partial } {\partial\theta_0} \left ( f(g(\theta_0)) \right) & \text{où }g(\theta_0)=\theta_0 + \theta_1 x^{(i)} - y^{(i)} \text{ et } f(u) = u² \\ & = \frac{1}{2m} \sum_{i=1}^m f'(g(\theta_0))g'(\theta_0) & \text{par dérivation de fonctions composées }\\ & = \frac{1}{2m} \sum_{i=1}^m 2 g(\theta_0) g'(\theta_0) & \text{car }f'(u) = 2u \\ & = \frac{1}{2m} \sum_{i=1}^m 2 \left( \theta_0 + \theta_1 x^{(i)} - y^{(i)} \right) g'(\theta_0) & \text{par définition de }g(\theta_0)\\ & = \frac{1}{2m} \sum_{i=1}^m 2 \left( \theta_0 + \theta_1 x^{(i)} - y^{(i)} \right) & \text{car }\theta_0 = cte \rightarrow g'(\theta_0)=1 \\ & = 2 \frac{1}{2m} \sum_{i=1}^m \left( \theta_0 + \theta_1 x^{(i)} - y^{(i)} \right) & \text{en sortant la constante multiplicative 2 } \\ & = \frac{1}{m} \sum_{i=1}^m \left( h_{\theta}(x^{(i)}) - y^{(i)} \right) & 2 \frac{1}{2m}= \frac{1}{m} \text{ et par définition de la fonction }h_{\theta}(x) \end{align} $$De la même façon, les détails du calcul de $\frac {\partial J(\theta)} {\partial\theta_1}$ :
$$ \begin{align} \frac {\partial J(\theta)} {\partial\theta_1} & = \frac {\partial } {\partial\theta_1} \left ( \frac{1}{2m} \sum_{i=1}^m \left( h_{\theta}(x^{(i)}) - y^{(i)} \right)^2 \right) & \text{par définition de }J(\theta) \\ & = \frac{1}{2m} \frac {\partial } {\partial\theta_1} \left (\sum_{i=1}^m \left( h_{\theta}(x^{(i)}) - y^{(i)} \right)^2 \right) & \text{on extrait la constante} \\ & = \frac{1}{2m} \sum_{i=1}^m \frac {\partial } {\partial\theta_1} \left (\left( h_{\theta}(x^{(i)}) - y^{(i)} \right)^2 \right) & \text{car dérivée de la somme = somme des dérivées} \\ & = \frac{1}{2m} \sum_{i=1}^m \frac {\partial } {\partial\theta_1} \left (\left( \theta_0 + \theta_1 x^{(i)} - y^{(i)} \right)^2 \right) & \text{par définition de }h_{\theta} \\ & = \frac{1}{2m} \sum_{i=1}^m \frac {\partial } {\partial\theta_1} \left ( f(g(\theta_1)) \right) & \text{où }g(\theta_1)=\theta_1 + \theta_1 x^{(i)} - y^{(i)} \text{ et } f(u) = u² \\ & = \frac{1}{2m} \sum_{i=1}^m f'(g(\theta_1))g'(\theta_1) & \text{par dérivation de fonctions composées }\\ & = \frac{1}{2m} \sum_{i=1}^m 2 g(\theta_1) g'(\theta_1) & \text{car }f'(u) = 2u \\ & = \frac{1}{2m} \sum_{i=1}^m 2 \left( \theta_0 + \theta_1 x^{(i)} - y^{(i)} \right) g'(\theta_1) & \text{par définition de }g(\theta_1)\\ & = \frac{1}{2m} \sum_{i=1}^m 2 \left( \theta_0 + \theta_1 x^{(i)} - y^{(i)} \right)x^{(i)} & \text{car }\theta_1 = cte \rightarrow g'(\theta_1)=x^{(i)} \\ & = 2 \frac{1}{2m} \sum_{i=1}^m \left( \theta_0 + \theta_1 x^{(i)} - y^{(i)} \right)x^{(i)} & \text{en sortant la constante multiplicative 2 } \\ & = \frac{1}{m} \sum_{i=1}^m \left( h_{\theta}(x^{(i)}) - y^{(i)} \right)x^{(i)} & 2 \frac{1}{2m}= \frac{1}{m} \text{ et par définition de la fonction }h_{\theta}(x^{(i)}) \end{align} $$Récapitulons:
l'algorithme d'apprentissage de la régression linéaire est, en pseudo-code :
début
initialiser $theta
répéter jusqu'à convergence:
$theta <- $theta - $alpha * $gradient
finNous avons calculé le gradient, pour rappel, voici sa formule :
$$ \frac {\nabla J(\theta)} {\nabla\theta} = \begin{bmatrix} \frac {\partial J(\theta)} {\partial\theta_0} \\ \frac {\partial J(\theta)} {\partial\theta_1} \end{bmatrix} = \begin{bmatrix} \frac{1}{m} \sum_{i=1}^m \left( h_{\theta}(x^{(i)}) - y^{(i)} \right) \\ \frac{1}{m} \sum_{i=1}^m \left( h_{\theta}(x^{(i)}) - y^{(i)} \right)x^{(i)} \end{bmatrix}\\ $$Nous disposons désormais de tous les éléments pour l'implémenter. Nous allons nous efforcer de vectoriser l'implémentation de celui-ci
Implémentation en Python¶
Dans un premier temps, nous allons utiliser des données jouet très très simples, nous les avons utilisé d'ailleurs précédemment.
import matplotlib.pyplot as plt
x=range(5)
y=[1.5,2.,2.5,3.,3.5]
plt.ylim(0,5)
plt.xlim(0,5)
plt.plot(x,y,'bx')
plt.show()
Les points sont parfaitement alignés, nous allons chercher à ajuster une droite sur ces points en implémentant l'algorithme de la régression linéaire. C'est parti !
class RegressionLineaire:
def __init__(self, x, y, theta, alpha=.05, nb_steps=1000):
assert len(x) == len(y)
self.theta=np.array(theta)
bias_terms=np.ones(len(x))
self.x=np.stack((bias_terms,x),axis=1)
self.y=y
self.alpha=alpha
self.nb_steps=nb_steps
self.m=len(y)
def compute_outputs(self):
output = np.dot(self.x,self.theta)
self.last_output=output
return output
def compute_cost(self):
return .5*np.sum(np.subtract(self.last_output, self.y)**2)
def compute_gradient(self):
delta = np.subtract(self.last_output, self.y)[:,np.newaxis]*self.x
gradient = np.sum(delta,axis=0)/self.m
self.gradient=gradient
return gradient
def maj_theta(self):
self.theta -= self.alpha*self.gradient
def fit(self):
hist=[]
for i in range(self.nb_steps):
self.compute_outputs()
self.compute_gradient()
self.maj_theta()
hist.append(self.compute_cost())
return hist
y=np.array([1.5,2.,2.5,3.,3.5])
x=np.array(range(5))
nb_steps=1000
model = RegressionLineaire(x=x,y=y, theta=[3.1,-2.7], nb_steps=nb_steps)
y2=model.theta[1]*x+model.theta[0]
print("La droite à l'état initial")
plt.plot(x, y,'bx')
plt.plot(x, y2, color='r', linestyle='-', linewidth=2)
plt.show()
hist=model.fit()
print("La droite après la régression linéaire")
y2=model.theta[1]*x+model.theta[0]
plt.plot(x, y,'bx')
plt.plot(x, y2, color='r', linestyle='-', linewidth=2)
plt.show()
print("l'historique de la fonction de coût")
plt.plot(range(nb_steps),hist)
plt.show()
Essayons désormais avec un autre jeu d'exemples:
X=np.linspace(0,5,100)
Y=1.5*X+2
noise=np.random.normal(0,1,100)
Y2=Y+noise
plt.scatter(X,Y2,label='y=1.5x+c+N(\mu=0,\sigma=1)')
plt.legend()
plt.title('Noise data points')
plt.show()
nb_steps=1000
model = RegressionLineaire(x=X,y=Y2, theta=[3.1,-2.7], nb_steps=nb_steps)
y2=model.theta[1]*X+model.theta[0]
print("La droite à l'état initial")
plt.plot(X, Y2,'bx')
plt.plot(X, y2, color='r', linestyle='-', linewidth=2)
plt.show()
hist=model.fit()
print("La droite après la régression linéaire")
y2=model.theta[1]*X+model.theta[0]
plt.plot(X, Y2,'bx')
plt.plot(X, y2, color='r', linestyle='-', linewidth=2)
plt.show()
print("l'historique de la fonction de coût")
plt.plot(range(nb_steps),hist)
plt.show()
Notre algorithme a l'air de bien marcher. Après 1000 étapes, la droite est presque presque parfaitement ajustée aux points.
Avec les données légèrement bruitées, le résultat semble correct aussi. Et l'algorithme est relativement rapide (contrairement à une précédente implémentation des k plus proches voisins haha)
Explication du code¶
def __init__(self, x, y, theta, alpha=.05, nb_steps=1000):
assert len(x) == len(y)
self.theta=theta
bias_terms=np.ones(len(x))
self.x=np.stack((bias_terms,x),axis=1)
self.y=y
self.alpha=alpha
self.nb_steps=nb_steps
self.m=len(y)
On définit alpha et le nombre d'itérations en paramètres par défaut.
On voulait initialement définir aussi theta avec une initialisation par défaut, mais le problème est qu'après chaque instantiation de RegessionLrineaire, theta n'était pas réinitialisé.
Expliquons maintenant:
bias_terms=np.ones(len(x))
self.x=np.stack((bias_terms,x),axis=1)
on ajoute un terme de biais aux exemples x afin de calculer $h_{\theta}(x)$ directement avec un produit scalaire.
passons maintenant à
def compute_outputs(self):
output = np.dot(self.x,self.theta)
self.last_output=output
return output
On a vectorisé l'opération tel que décrit dans le schéma juste au-dessus. On s'éloigne très légèrement de la formule $h_{\theta}(x)=\theta^T x$, mais à l'usage, on constate une très nette amélioration du temps de calcul en n'effectuant pas de transposée (devient très visible quand on définit le nombre d'itérations à 10.000 ou plus.
Et dans un sens, on colle au maths, car $\theta^T x$ est bien le produit interne, ou "dot product" en anglais de entre $\theta$ et $x$ (en espérant que je commets par d'erreur au niveau des maths)
def compute_cost(self):
return .5*np.sum(np.subtract(self.last_output, self.y)**2)
rien à dire, np.subtract effectue un soustraction terme à terme, np.sum effectue la somme des éléments du vecteur résultat. On colle parfairement à la formule mathématique.
def compute_gradient(self):
delta = np.subtract(self.last_output, self.y)[np.newaxis].T*self.x
gradient = np.sum(delta,axis=0)/self.m
self.gradient=gradient
return gradient
On rappelle la formule du gradient:
$$ \frac {\nabla J(\theta)} {\nabla\theta} = \begin{bmatrix} \frac {\partial J(\theta)} {\partial\theta_0} \\ \frac {\partial J(\theta)} {\partial\theta_1} \end{bmatrix} = \begin{bmatrix} \frac{1}{m} \sum_{i=1}^m \left( h_{\theta}(x^{(i)}) - y^{(i)} \right) \\ \frac{1}{m} \sum_{i=1}^m \left( h_{\theta}(x^{(i)}) - y^{(i)} \right)x^{(i)} \end{bmatrix}\\ $$Voici schématiquement comment on l'implémente :
en (1), on cherche à multiplier $(h_{\theta}(x^{(i)}) - y^{(i)})$ terme à terme avec les exemples d'entrée auxquels on a concaténé le terme de biais.
Pour cela, on augmente le rang du vecteur $ \begin{pmatrix} h_{\theta}(x^{(1)}) - y^{(1)} \\ h_{\theta}(x^{(2)}) - y^{(2)} \\ h_{\theta}(x^{(3)}) - y^{(3)} \\ h_{\theta}(x^{(4)}) - y^{(4)} \\ h_{\theta}(x^{(5)}) - y^{(5)} \end{pmatrix}$, de 1 à 2, on fait en sorte que chaque élément du vecteur soit une ligne d'une matrice en 2 dimensions de rang 2? de manière à ce que chaque élément puisse être multiplié terme par terme à chaque ligne de la matrice $ \begin{bmatrix} 1 & x^{(1)} \\ 1 & x^{(2)} \\ 1 & x^{(3)} \\ 1 & x^{(4)} \\ 1 & x^{(5)} \end{bmatrix}$
Illustration avec du code Python à la volée :
a = np.array(range(5))
print(a)
print(a.shape)
Ici, a un est vecteur de rang 1
a2 = a[:,np.newaxis]
print(a2)
print(a2.shape)
Et ici, on a transformé le vecteur a, de rang 1, en un matrice de rang 2, où chaque élément est une ligne à lui tout seul.
on a fait [:,np.newaxis] et pas [np.newaxis], car sinon on aurait obtenu :
a3 = a[np.newaxis]
print(a3)
print(a3.shape)
on n'aurait eu qu'une matrice de rang 2, mais avec une seule ligne contenant tous les éléments du vecteur de départ. Et la multiplication terme à terme n'aurait pas réussi:
b = np.array(range(10)).reshape((5,2))
print(b)
print(b.shape)
print(a)
a[:,np.newaxis]*b
a[np.newaxis]*b
prochaine méthode :
def maj_theta(self):
self.theta -= self.alpha*self.gradient
On applique purement et simplement l'algo : $$ \theta_{t+1} = \theta_{t} - \alpha \frac {\partial J(\theta)} {\partial \theta_{t}} $$
Ou, sous forme vectorisée
$$ \theta_{t+1} = \theta_{t} - \alpha \frac {\nabla J(\theta)} {\nabla \theta_{t}} $$Où, encore une fois : $$ \frac {\nabla J(\theta)} {\nabla\theta} = \begin{bmatrix} \frac {\partial J(\theta)} {\partial\theta_0} \\ \frac {\partial J(\theta)} {\partial\theta_1} \end{bmatrix} = \begin{bmatrix} \frac{1}{m} \sum_{i=1}^m \left( h_{\theta}(x^{(i)}) - y^{(i)} \right) \\ \frac{1}{m} \sum_{i=1}^m \left( h_{\theta}(x^{(i)}) - y^{(i)} \right)x^{(i)} \end{bmatrix}\\ $$
que l'on a calculé et implémenté juste au-dessus.
Enfin:
def fit(self):
hist=[]
for i in range(self.nb_steps):
self.compute_outputs()
self.compute_gradient()
self.maj_theta()
hist.append(self.compute_cost())
return hist
On calcule le gradient et on enregistre le coût, pour analyse et debug. On répète l'opération jusqu'à la terminaison. Ici notre critère de terminaison est simplement répéter l'opération n fois.
Voila, c'est tout pour cet article. Dans une prochaine mise à jour, ou un prochain article, il y aura l'implémentation Tensorflow, pour coller à l'article précédent "prise en main de Tensorflow".
Le prochain algorithme étudié et implémenté sera très certainement la régression logistique.
A la prochaine