9 Moindres carrés
Les modèles linéaires sont les modèles les plus simples dans lesquels on raisonne en termes d’entrées et de sorties. Dans ces modèles, on dispose de variables \(x_i\), dites explicatives, et de variables \(y_i\), dites à expliquer, et l’on suppose qu’il existe une fonction inconnue \(f\) telle que \[y_i \approx f(x_i),\] et que l’on voudrait estimer. Les modèles linéaires consistent à supposer que \(f\) est affine. Les modèles plus complexes, commes les réseaux de neurones, placent \(f\) dans des classes plus riches. L’objectif de la méthode des moindres carrés est de trouver la meilleure approximation de \(f\) possible dans la classe des fonctions affines.
9.1 Ajustement affine en une dimension.
On suppose qu’il existe entre les données \(x_i\) et \(y_i\) une relation de la forme \(y_i \approx \alpha + \beta x_i\) où \(\alpha, \beta\) sont deux nombres réels. Ici, \(\approx\) signifie que la relation n’est pas parfaite : peut-être par exemple que les sorties sont bien égales à \(\alpha+\beta x_i\), mais que les observations \(y_i\) ont été polluées par du bruit ou des erreurs. Nous verrons cela plus tard.
Pour l’heure, nous voulons chercher les meilleurs \(\alpha, \beta\) possibles. On calcule la distance entre le nuage de points \((x_i, y_i)\) et la droite d’équation \(y = \alpha + \beta x\). Cette distance au carré est donnée par \[ L(\alpha, \beta) = \sum_{i=1}^n (y_i - \alpha - \beta x_i)^2.\]
On cherchera donc les \((\hat\alpha, \hat\beta)\) qui minimisent cette distance. La fonction \(L\) est manifestement une fonction quadratique qui tend vers \(+\infty\) lorsque \((\alpha, \beta) \to \infty\), par conséquent cette fonction possède un unique minimiseur \((\hat{\alpha}, \hat{\beta})\), et ce minimiseur est le seul point en lequel les dérivées partielles s’annulent (conditions de premier ordre) : \(\partial_\alpha L(\hat{\alpha}, \hat{\beta}) = 0\) et \(\partial_\beta L(\hat{\alpha}, \hat{\beta})=0\). Or, \[\partial_\alpha L(\alpha, \beta) = \sum_{i=1}^n (\alpha + \beta x_i - y_i)\] \[ \partial_\beta L(\alpha, \beta) = \sum_{i=1}^n x_i(\alpha + \beta x_i - y_i).\] Les conditions de premier ordre deviennent donc \(n\alpha + \beta (x_1 + \dotsc + x_n) - (y_1 + \dotsb + y_n) =0\) soit encore \(\alpha + \beta \bar{x} - \bar{y}=0\), et d’autre part \(\alpha (x_1 + \dotsb + x_n) + \beta(x_1^2 + \dotsb + x_n^2) - (x_1y_1 + \dotsb + x_ny_n) = 0\), soit \(\alpha \bar{x}+ \beta \overline{xx} - \overline{xy} = 0\), où \(\overline{xx}\) est la moyenne des carrés des \(x_i\) et \(\overline{xy}\) la moyenne des \(x_iy_i\). En résolvant ces équations, on trouve d’abord \(\alpha\) puis \(\beta\) : \[\beta = \frac{\overline{xy}-\bar{x}\bar{y}}{\overline{xx} - \bar{x}\bar{x}} ,~~\quad\alpha = \bar{y} - \hat{\beta}\bar{x}.\] Le coefficient \(\beta\) n’est rien d’autre que la covariance empirique des \(x_i\) et des \(y_i\), normalisé par la variance empirique des \(x_i\).
L’inégalité de Cauchy-Schwartz dit que \(\left|\overline{xy} - \bar{x}{\bar{y}}\right| \leqslant \tilde{\sigma}_x \tilde{\sigma}_y\), où l’on a noté \[\tilde{\sigma}_x^2 = \frac{1}{n}\sum_{i=1}^n x_i^2 - \left(\frac{1}{n}\sum_{i=1}^n x_i\right)^2\] l’estimateur naïf de la variance[^1]. L’inégalité n’est une égalité que si \(x\) et \(y\) sont effectivement colinéaires, c’est-à-dire si \(y_i = \hat{\alpha} + x_i \hat{\beta}\) pour tous les \(i\). La qualité de l’ajustement affine est donc bien mesurée par la quantité \[ R^2 = \frac{\overline{xy}-\bar{x}\bar{x}}{\tilde{\sigma}_x \tilde{\sigma}_y}.\]
9.2 Moindres carrés ordinaires
Dans le cadre général le nombre \(d\) de variables explicatives est plus grand que 1. On notera \(\bx = (x_1, \dotsc, x_d)\) un élément de \(\mathbb{R}^d\) ; les variables explicatives seront alors \(\bx_1,\dotsc, \bx_n\). Avec mes notations, ces vecteurs sont des vecteurs lignes1.
On cherchera donc des nombres \(\theta_i\) tels que \(y_i\) est aussi proche que possible de \[\theta_1 \bx_{i,1} + \dotsb + \theta_d \bx_{i,d} = \bx_i \bt, \] oùe paramètre \(\bt\), sera toujours vu comme le vecteur colonne2 des \(\theta_i\).
Remarque : où est passée la constante ? Dans l’équation ci-dessus, on a l’impression que le terme constant, qui correspondait à \(\alpha\) dans l’exemple en dimension 1, a disparu. Ce n’est pas le cas : intégrer la constante au modèle revient à considérer que la variable constante égale à 1 fait partie des variables explicatives. En pratique, cela revient à poser, par exemple, \(\bx_{i,1} = 1\) pour tout \(i\). Ainsi, la constante correspondra toujours à \(\theta_1\).
On pose \(X\) la matrice \(n\times d\) dont la \(i\)-ème ligne est \(\bx_i\) et \(Y\) le vecteur colonne des \(y_i\). La matrice dont les lignes sont composées des nombres réels \(\bx_i \bt\) n’est autre que la matrice \(X\bt\). De façon générale, pour n’importe quel \(\theta \in \mathbb{R}^d\), la distance entre le nuage de points \((X,Y)\) et la droite d’équation \(Y = X\theta\) est alors \(|Y - X\theta|\). On pourrait reproduire la méthode analytique ci-dessus pour trouver les paramètres optimaux, à savoir \[\hat{\bt} = \arg \min_{\theta} |Y - X\theta|^2. \tag{9.1}\] Cependant, une interprétation géométrique simplifie la tâche : le \(\hat{\bt}\) qui minimise Équation 9.1 est précisément celui qui garantit que \(X\hat{\bt}\) est la projection orthogonale de \(Y\) sur le sous-espace vectoriel \(\mathscr{V}_X = \{X\theta : \theta \in \mathbb{R}^d\}\).
Théorème 9.1 Si \(d\leqslant n\) et si \(X\) est de rang \(d\), alors
\[\hat{\bt} = (X^\top X)^{-1}X^\top Y. \tag{9.2}\]
Preuve. La projection orthogonale sur le sous-espace vectoriel engendré par les colonnes d’une matrice \(X\) est la matrice \(X(X^\top X)^{-1}X^\top\), comme démontré dans l’appendice Chapitre 20. Ainsi, la projection de \(Y\) sur ce sous-espace est \(X(X^\top X)^{-1}X^\top Y\), et c’est aussi (par définition de l’argmin) \(X \hat{\bt}\). Comme \(X\) est injective en vertu du théorème du rang, on en déduit le résultat.
L’expression Équation 9.2 possède de nombreuses expressions alternatives. Parmi elles, on pourra noter que \[\hat{\bt} = \bt + \left(\frac{1}{n}\sum_{i=1}^n \bx_i^\top\bx_i \right)^{-1}\frac{1}{n}\sum_{i=1}^n \bx_i^\top \varepsilon_i. \tag{9.3}\]
Remarque générale. Si, en dimension 1, on cherche à trouver le \(\theta\) qui résout l’équation \(y = x\theta\), on trouve évidemment \(\theta = y/x\), c’est-à-dire qu’on divise \(y\) par \(x\), ou encore qu’on multiplie \(y\) par l’inverse de \(x\). En dimension supérieure, quand on veut résoudre en \(\theta\) l’équation \(Y = X\theta\), c’est pareil. Le problème, c’est qu’on ne sait pas forcément inverser \(X\). La formule Équation 9.2 dit que même si \(X\) n’est pas inversible, on peut quand même “diviser par \(X\)” : c’est pour cela que la matrice \((X^\top X)^{-1} X^\top\) est appelée pseudo-inverse à gauche de \(X\) – parfois associée au nom de Moore-Penrose. Multiplier \(Y\) par \((X^\top X)^{-1}X^\top\) donne \((X^\top X)^{-1}X^\top Y = \hat{\bt}\) : ce vecteur ne vérifie par forcément \(X{\hat{\bt}} = Y\), mais parmi tous les vecteurs possibles, c’est celui qui rend \(X\bt\) le plus proche possible de \(Y\).
9.3 Résidus et \(R^2\)
Le vecteur \(\hat{Y} = X\hat{\bt}\) est appelé vecteur des prédictions. Le vecteur \(\hat{\varepsilon} = Y-\hat{Y} = Y - X\hat{\bt}\) est appelé vecteur des résidus. Si ce dernier est nul ou très petit, cela veut dire que les \(Y\) sont presque parfaitement des fonctions linéaires des \(X\).
Définition 9.1 Dans le cas d’une régression Équation 9.2 avec constante, le coefficient de détermination est défini par \[R^2 = \frac{\sum_{i=1}^n |\hat{y}_i - \bar{y}|^2}{\sum_{i=1}^n |y_i - \bar{y}|^2}. \] C’est un nombre entre 0 et 1.
Le numérateur est la variance empirique des prédictions \(\hat{y}_i\). Le dénominateur est la variance empirique des observations. Dans les des cas, il s’agit de la norme carrée d’un vecteur (\(\hat{Y}\) et \(Y\)) projeté sur l’espace des vecteurs de moyenne nulle. Comme \(\hat{Y}\) est déjà une projection de \(Y\) sur un certain sous-espace, on a forcément \(|\hat{Y} - \bar{y}| \leqslant |Y - \bar{y}|\), donc le coefficient \(R^2\) est toujours entre 0 et 1.
Plus le coefficient de détermination est proche de 1, meilleure est la régression – attention, cet indicateur possède de nombreuses limites.