ML1/ML1_regression.Rmd

---
title: "Régression"
output: html_notebook
---

\renewcommand{\vec}[1]{\mathbf{#1}}

# Concepts de base

$\vec{x}^{(1)} \dots \vec{x}^{(P)}$ sont des vecteurs de $\mathbb{R}^N$ associés aux valeurs, aussi appelées étiquettes, $y^{(1)} \dots y^{(P)}$ de $\mathbb{R}$. On cherche une fonction $f(\vec{x}) : \mathbb{R}^N \rightarrow \mathbb{R}$ qui modélise la relation entre les observations $\vec{x}$ et les étiquettes $y$.

La fonction $f$ peut avoir une forme paramétrique, comme par exemple :
$$
f(\vec{x}) = a_0 + a_1x_1 + a_2x_2 + \dots + a_Nx_N 
$$
Si $P=N+1$, les paramètres $a_0, a_1, \dots a_N$ sont solutions d'un système linéaire :
$$
\begin{cases}
y^{(1)} &= a_0 + a_1 x_1^{(1)} + a_2 x_2^{(1)} + \dots + a_N x_N^{(1)} \\
y^{(2)} &= a_0 + a_1 x_1^{(2)} + a_2 x_2^{(2)} + \dots + a_N x_N^{(2)} \\
\dots &= \dots \\
y^{(P)} &= a_0 + a_1 x_1^{(P)} + a_2 x_2^{(P)} + \dots + a_N x_N^{(P)} \\
\end{cases}
$$
Ce système s'écrit également sous forme matricielle :
$$
\left( \begin{array}{cccc}
1 & x^{(1)}_1 & \dots & x^{(1)}_N \\
1 & x^{(2)}_1 & \dots & x^{(2)}_N \\
\dots & \dots & \dots & \dots \\
1 & x^{(P)}_1 & \dots & x^{(P)}_N
\end{array} \right)
\left( \begin{array}{c}
a_0 \\ a_1 \\ \dots \\ a_N 
\end{array} \right)
=
\left( \begin{array}{c}
y^{(1)} \\ y^{(2)} \\ \dots \\ y^{(P)} 
\end{array} \right)
$$
Chaque ligne $i$ de la matrice du terme de gauche de l'égalité ci-dessus est le vecteur ligne $\vec{x}^{(i)T}$ avec l'addition d'un premier terme constant qui correspond au paramètre $a_0$. En nommant cette matrice $\vec{X}^T$, le système linéaire ci-dessus s'écrit :
$$
\vec{X}^T \vec{a} = \vec{y}
$$
On considère le cas particulier où $x$ est un scalaire et $f$ est un polynome de degré $N$ :
$$
f(x) = a_0 + a_1x + a_2x^2 + \dots + a_Nx^N
$$
Avec $P=N+1$ observations et les étiquettes associées $\left( x^{(k)}, y^{(k)}  \right)$, les coefficients de ce polynome sont solution d'un système linéaire :

$$
\left( \begin{array}{ccccc}
1 & x^{(1)} & (x^{(1)})^2 & \dots & (x^{(1)})^N \\
1 & x^{(2)} & (x^{(2)})^2 & \dots & (x^{(2)})^N \\
\dots & \dots & \dots & \dots \\
1 & x^{(P)} & (x^{(P)})^2 & \dots & (x^{(P)})^N
\end{array} \right)
\left( \begin{array}{c}
a_0 \\ a_1 \\ \dots \\ a_N 
\end{array} \right)
=
\left( \begin{array}{c}
y^{(1)} \\ y^{(2)} \\ \dots \\ y^{(P)} 
\end{array} \right)
$$

La matrice du terme de gauche de l'égalité ci-dessus est traditionnellement appelée "matrice de Vandermonde".

# Exemple sur un jeu de données synthétique

On considère un exemple de fonction non-linéaire.
```{r}
f <- function(x) {exp(x) * cos(2*pi*sin(pi*x))}
```
On utilise cette fonction pour générer un jeu de données avec l'ajout d'un bruit gaussien.
```{r}
gendat <- function(N, sd) {
  # N: nombre d'observations à générer
  # sd: écart type d'un bruit gaussien de moyenne nulle
  X = runif(N)
  Y = f(X) + rnorm(N, mean=0, sd=sd)
  list(X = X, Y = Y)
}
```
On se donne une fonction pour afficher simultanément un jeu de données et la fonction utilisée pour le générer.
```{r}
plt <- function(data, f) {
  xs = seq(0,1,length.out=100)
  plot(xs, f(xs), type="l")
  points(data$X, data$Y)
}
```
On affiche un exemple de jeu de données.
```{r}
set.seed(1123)
data = gendat(5,0.2)
plt(data,f)
```
La fonction ci-dessous résout le système linéaire correpsondant à la matrice de Vandermonde et retourne les coefficients du polynôme résultat.
```{r}
polyreg1 <- function(data) {
  vandermonde = outer(data$X, 0:(length(data$X)-1), "^")
  solve(vandermonde, data$Y)
}
```
La fonction ci-dessous évalue un polynôme en un point `x` étant donné ses coefficients `coef`.
```{r}
polyeval <- function(x,coef) {
  powers = 0:(length(coef)-1)
  f = function(y) { sum(coef * y^powers) }
  sapply(x,f)
}
```
La fonction ci-dessous ajoute au graphe courant le tracé en pointillés d'un polynôme défini par ses coefficients.
```{r}
pltpoly <- function(coef) {
  xs = seq(0,1,length.out=100)
  lines(xs, polyeval(xs,coef), lty="dotted")
}
```
On affiche la fonction génératrice, le jeu de donnée et le polynôme qui passe par chaque point du jeu de données.
```{r}
coef = polyreg1(data)
plt(data,f)
pltpoly(coef)
```
Ce polynôme qui passe exactement par chaque observation a peu de chance d'offrir de bonnes capacités prédictives. Vérifier par exemple que, sur notre exemple synthétique, pour cinq points  générés à partir de la fonction $f$ et avec l'ajout d'un bruit gaussien (par exemple d'écart type $0.2$), le polynôme découvert, de degré quatre, peut être très éloigné de la fonction génératrice. C'est un exemple du phénomène de sur-apprentissage. Pour limiter ce problème, on cherche à découvrir un polynôme de degré plus faible, qui ne passera donc pas exactement par toutes les observations, mais qui aura probablement une meilleure capacité à prédire les étiquettes de nouvelles observations.

# Généralisation à un espace de fonctions

On considère un espace vectoriel composé de fonctions. Une base de cet espace est un ensemble de fonctions ($f_1, f_2, \dots f_N$) tel que toute fonction de l'espace s'exprime comme combinaison linéaire des fonctions de base.
$$
f(\vec{x}) = a_1 f_1(\vec{x}) + a_2 f_2(\vec{x}) + \dots + a_N f_N(\vec{x})
$$
Pour un jeu de données $\left\{ \vec{x^{(k)}},y^{(k)}\right\}_{k=1}^{N}$ de taille $N$, les coefficients $a_i$ sont solutions d'un système linéaire.
$$
\left( \begin{array}{ccccc}
f_1(\vec{x^{(1)}}) & f_2(\vec{x^{(1)}}) & \dots & f_N(\vec{x^{(1)}}) \\
f_1(\vec{x^{(2)}}) & f_2(\vec{x^{(2)}}) & \dots & f_N(\vec{x^{(2)}}) \\
\dots & \dots & \dots & \dots \\
f_1(\vec{x^{(N)}}) & f_2(\vec{x^{(N)}}) & \dots & f_N(\vec{x^{(N)}})
\end{array} \right)
\left( \begin{array}{c}
a_1 \\ a_2 \\ \dots \\ a_N 
\end{array} \right)
=
\left( \begin{array}{c}
y^{(1)} \\ y^{(2)} \\ \dots \\ y^{(N)} 
\end{array} \right)
$$
On note ce système linéaire $\vec{Ax} = \vec{b}$.

# Approche des moindres carrés

Le système linéaire $\vec{Ax} = \vec{b}$ avec $\vec{A} \in \mathbb{R}^{M \times N}$ n'a pas de solution quand le nombre d'observations dépasse le nombre de fonctions de base (c'est-à-dire, $M>N$). Une approche est alors de chercher une approximation $\vec{Ax} \approx \vec{b}$ qui minimise la somme des carrés des erreurs : $\|\vec{A}\vec{x}-\vec{b}\|^2_2$.

$$
\begin{align*}
 & \|\vec{A}\vec{x}-\vec{b}\|^2_2 \\
= \{ & \|\vec{x}\|_2 = \sqrt{\vec{x}\cdot\vec{x}} \} \\
 & \left(\vec{A}\vec{x}-\vec{b}\right) \cdot \left(\vec{A}\vec{x}-\vec{b}\right) \\
= \{ & \text{Par définition du produit scalaire euclidien} \} \\
 & \left(\vec{A}\vec{x}-\vec{b}\right)^T \left(\vec{A}\vec{x}-\vec{b}\right) \\
= \{ & \text{propriété de la transposition} \} \\
 & \left(\vec{x}^T\vec{A}^T - \vec{b}^T \right) \left(\vec{A}\vec{x}-\vec{b}\right) \\
= \{ & \text{multiplication} \} \\
 & \vec{x}^T\vec{A}^T\vec{A}\vec{x} - \vec{x}^T\vec{A}^T\vec{b} - \vec{b}^T\vec{A}\vec{x} + \vec{b}^T\vec{b} \\
= \{ & \vec{b}^T\vec{A}\vec{x} \text{ étant une valeur scalaire, } \vec{b}^T\vec{A}\vec{x} = \left(\vec{b}^T\vec{A}\vec{x}\right)^T = \vec{x}^T\vec{A}^T\vec{b} \} \\
 & \vec{x}^T\vec{A}^T\vec{A}\vec{x} - 2\vec{x}^T\vec{A}^T\vec{b} + \vec{b}^T\vec{b}
\end{align*}
$$
Cette dernière expression quadratique en $\vec{x}$ correspond à une surface convexe. Donc, le minimum de cette expression peut être calculé en annulant sa dérivée (penser à une courbe $y = a+bx+cx^2$ dont l'unique extremum est atteint lorsque la pente est nulle).
$$
\begin{align*}
 & \vec{0} = 2\vec{A}^T\vec{A}\vec{x} - 2\vec{A}^T\vec{b} \\
=& \\
 & \vec{A}^T\vec{A}\vec{x} = \vec{A}^T\vec{b}
\end{align*}
$$
Ainsi, quand $M>N$, la solution approximée $\vec{x}$, telle que $\vec{Ax} \approx \vec{b}$ en minimisant la somme des carrés des erreurs, est la solution du système linéaire suivant où $\vec{A}^T\vec{A}$ est appelée la matrice de Gram.
$$
\vec{A}^T\vec{A} \vec{x} = \vec{A}^T\vec{b}
$$

## Approche des moindres carrés appliquée à la régression polynomiale

Pour un polynôme de degré $N-1$, les fonctions de bases mentionnées ci-dessus sont : $f_1(x)=1$, $f_2(x)=x$, $f_3(x)=x^2$,..., $f_N(x)=x^{N-1}$. Elles permettent de définir la matrice $\vec{A}$ et la matrice de Gram $\vec{A}^T\vec{A}$.
La fonction ci-dessous résout le système linéaire correspondant à la matrice de Gram pour un polynôme de degré fixé. Elle retourne les coefficients du polynôme résultat.
```{r}
polyreg2 <- function(data, degre) {
  A = outer(data$X, 0:degre, "^")
  gram = t(A) %*% A
  solve(gram, as.vector(t(A) %*% data$Y))
}
```
Sur notre exemple synthétique, on affiche la fonction génératrice, le jeu de donnée et le "meilleur" polynôme de degré 3.
```{r}
coef = polyreg2(data,3)
plt(data,f)
pltpoly(coef)
```
Ce polynôme de degré trois modélise mieux la fonction génératrice inconnue que celui de degré quatre qui ne commettait aucune erreur sur les données observées.

# Régularisation de Tikhonov

Avec moins d'observations que de fonctions de base ($M<N$), le système $\vec{A}\vec{x}=\vec{b}$ ne possède pas de solution unique. Même quand $M \geq N$, le système linéaire peut posséder une solution approchée préférable à la solution optimale. C'est en particulier vrai quand plusieurs observations sont très proches. Par exemple, soit le système linéaire suivant :
$$
\left( \begin{array}{cc}
1 & 1 \\
1 & 1.00001 \\
\end{array} \right)
\left( \begin{array}{c}
x_1 \\ x_2
\end{array} \right)
=
\left( \begin{array}{c}
1 \\ 0.99 
\end{array} \right)
$$
Sa solution est $\vec{x}^T = (1001,-1000)$. Cependant, la solution approchée $\vec{x}^T = (0.5,0.5)$  est nettement préférable. En effet, la solution optimale a peu de chance de bien s'adapter à de nouvelles observations (par exemple, l'observation $(1,2)$ serait projetée sur l'étiquette -999).

Ainsi, lorsqu'il faut choisir entre plusieurs solutions, on préférera celles dont les coefficients (ou paramètres) ont de faibles valeurs. Par exemple, on pourra vouloir minimiser $|x_1|+|x_2|+\dots$ (aussi noté $\|\vec{x}\|_1$, la "norme 1") ou encore $x_1^2+x_2^2+\dots$ (aussi noté $\|\vec{x}\|_2^2$, le carré de la "norme 2"). Dans ce dernier cas, on peut chercher à résoudre un nouveau problème de minimisation :
$$
\begin{align*}
\min_{\vec{x}} \|\vec{A}\vec{x}-\vec{b}\|^2_2 + \alpha \|\vec{x}|^2_2 \\
avec \quad 0 \leq \alpha
\end{align*}
$$
Il s'agit à nouveau d'un problème de minimisation quadratique en $\vec{x}$ dont le minimum peut être découvert par annulation de la dérivée.