8 Régularisation des moindre carrés : Ridge, Lasso et elastic net

Exercice 1 (Questions de cours) A, B, B, B, A (pour un bon choix de \(\lambda\)) et B, A, C et D.

Exercice 2 (Projection et régression ridge)

Exercice 3 (Variance des valeurs ajustées avec une régression ridge)

Exercice 4 (Nombre effectif de paramètres de la régression ridge)

Rappelons que pour une valeur \(\kappa\) donnée, le vecteur de coefficients de la régression ridge s’écrit \[ \hat \beta_{\mathrm{ridge}}(\kappa) = (X'X + \kappa I)^{-1}X'Y. \] et donc l’ajustement par la régression ridge est \[ \hat Y_{\mathrm{ridge}}(\kappa)=X(X'X + \kappa I)^{-1}X'Y=H^*(\kappa)Y \]
Soit \(U_i\) le vecteur propre de \(A\) associé à la valeur propre \(d^2_i\). Nous avons donc par définition que \[ \begin{eqnarray*} AU_i&=&d^2_iU_i\\ AU_i+\lambda U_i&=&d^2_iU_i+\lambda U_i=(d^2_i+\lambda) U_i\\ (A+\lambda I_p)U_i&=&(d^2_i+\lambda) U_i, \end{eqnarray*} \] c’est-à-dire que \(U_i\) est aussi vecteur propre de \(A+\lambda I_p\) associé à la valeur propre \(\lambda+d^2_i\).
Nous savons que \(X=QD P'\) avec \(Q\) et \(P\) matrices orthogonales et \(D=\mathop{\mathrm{diag}}(d_1,\dotsc,d_p)\). Puisque \(Q\) est orthogonale, nous avons, par définition, \(Q'Q=I\). Nous avons donc que \(X'X=(QD P')'QD P'=PDQ'QDP'=PD^2P'\). Puisque \(P\) est orthogonale \(P'P=I_p\) et \(P^{-1}=P\). \[ \begin{eqnarray*} \mathop{\mathrm{tr}}(X(X'X+\lambda I_p)^{-1}X')&=&\mathop{\mathrm{tr}}((X'X+\lambda I_p)^{-1}X'X)\\ &=&\mathop{\mathrm{tr}}((PD^2P'+\lambda PP')^{-1}PD^2P')\\ &=&\mathop{\mathrm{tr}}((P(D+\lambda I_p )P')^{-1}PD^2P'). \end{eqnarray*} \] Nous avons donc \[ \begin{eqnarray*} \mathop{\mathrm{tr}}(X(X'X+\lambda I_p)^{-1}X')&=&\mathop{\mathrm{tr}}( (P')^{-1}(D+\lambda I_p )^{-1} P^{-1} PD^2P')\\ &=&\mathop{\mathrm{tr}}( (P')^{-1}(D+\lambda I_p )^{-1} D^2P')\\ &=&\mathop{\mathrm{tr}}( (D+\lambda I_p )^{-1} D^2). \end{eqnarray*} \] Selon la définition de \(H^*(\kappa)\), nous savons que sa trace vaut donc \[ \begin{eqnarray*} \mathop{\mathrm{tr}}( (D+\kappa I_p )^{-1} D^2). \end{eqnarray*} \] Comme \(D\) et \(I_p\) sont des matrices diagonales, leur somme et produit sont simplement leur somme et produit terme à terme des éléments de la diagonale, et donc cette trace (somme des éléments de la diagonale) vaut \[ \sum_{i=1}^{p}{\frac{d_j^2}{d_j^2+\kappa}}. \]

Exercice 5 (Estimateurs à rétrecissement - shrinkage)

Soit le modèle de régression \[ Y=X\beta+\varepsilon. \] En le pré-multipliant par \(P\), nous avons \[ Z=PY=PX\beta+P\varepsilon=DQ\beta+\eta=D\gamma+\eta. \] Puisque \(\varepsilon\sim\mathcal{N}(0,\sigma^2 I_n)\) et \(P\) fixé, nous avons que \(\eta=P\varepsilon\) suit une loi normale de moyenne \(\mathbf E(\eta)=P\mathbf E(\varepsilon)=0\) et de variance \(\mathop{\mathrm{V}}(\eta)=P\mathop{\mathrm{V}}(\varepsilon)P'=\sigma^2PP'=\sigma^2I_n\).

Par définition, \(Z\) vaut \(PY\) et nous savons que \(Y\sim\mathcal{N}(X\beta,\sigma^2 I_n)\), donc \(Z\sim\mathcal{N}(PX\beta,\sigma^2 PP')\), c’est-à-dire \(Z\sim\mathcal{N}(DQ\beta,\sigma^2 I_n)\) ou encore \(Z\sim\mathcal{N}(D\gamma,\sigma^2 I_n)\). En utilisant la valeur de \(D\) nous avons \[ \begin{eqnarray*} D\gamma&=& \begin{pmatrix} \Delta \gamma\\ 0 \end{pmatrix}. \end{eqnarray*} \] Donc \(Z_{1:p}\sim\mathcal{N}(\Delta\gamma,\sigma^2I_p)\).
Soit un estimateur de \(\beta\) linéaire en \(Y\)~: \(\hat \beta=AY\). Soit l’estimateur de \(\gamma=Q\beta\) linéaire en \(Y\)~: \(\hat\gamma=Q AY\). Pour calculer leur matrice de l’EQM, nous devons calculer leur biais et leur variance. Le biais de \(\hat \beta\) est \[ B(\hat \beta)=\mathbf E(\hat \beta)-\beta=\mathbf E(AY)-\beta=A\mathbf E(Y)-\beta=AX\beta-\beta. \] Le biais de \(\hat\gamma\) s’écrit \[ B(\hat\gamma)=\mathbf E(\hat \gamma)-\gamma=\mathbf E(Q\hat \beta)-\gamma=Q\mathbf E(\hat \beta)-\gamma=QAX\beta-\gamma. \] Comme \(\gamma=Q\beta\) et \(Q'Q=I_p\) nous avons \[ B(\hat\gamma)=QAXQ'\gamma-\gamma. \] La variance de \(\hat \beta\) s’écrit \[ \mathop{\mathrm{V}}(\hat \beta)=\mathop{\mathrm{V}}(AY)=A\mathop{\mathrm{V}}(Y)A'=\sigma^2 AA', \] et celle de \(\hat \gamma\) est \[ \mathop{\mathrm{V}}(\hat\gamma)=\mathop{\mathrm{V}}(Q\hat \beta)=Q\mathop{\mathrm{V}}(\hat \beta)Q'=\sigma^2 QAA'Q'. \] Nous en déduisons que les matrices des EQM sont respectivement \[ \begin{eqnarray*} \mathop{\mathrm{EQM}}(\hat \beta)&=&(AX\beta-\beta)(AX\beta-\beta)'+\sigma^2 AA',\\ \mathop{\mathrm{EQM}}(\hat \gamma)&=&(QAXQ'\gamma-\gamma)(QAXQ'\gamma-\gamma)' + \sigma^2 QAA'Q', \end{eqnarray*} \] et enfin les traces de ces matrices s’écrivent \[ \begin{eqnarray*} \mathop{\mathrm{tr}}(\mathop{\mathrm{EQM}}(\hat \beta))&=&(AX\beta-\beta)'(AX\beta-\beta)+\sigma^2\mathop{\mathrm{tr}}(AA'),\\ \mathop{\mathrm{tr}}(\mathop{\mathrm{EQM}}(\hat \gamma))&=&(QAXQ'\gamma-\gamma)'(QAXQ'\gamma-\gamma)+ \sigma^2\mathop{\mathrm{tr}}(AA').\\ \end{eqnarray*} \] Rappelons que \(\gamma=Q\beta\) et que \(Q'Q=I_p\), nous avons donc \[ \begin{eqnarray*} \mathop{\mathrm{tr}}(\mathop{\mathrm{EQM}}(\hat \gamma))&=&\gamma'(QAXQ'-I_p)'(QAXQ'-I_p)\gamma+ \sigma^2\mathop{\mathrm{tr}}(AA')\\ &=&\beta'(QAX - Q)'(QAX - Q)\beta+ \sigma^2\mathop{\mathrm{tr}}(AA')\\ &=&\beta'(AX-I_p)Q'Q(AX-I_p)\beta+ \sigma^2\mathop{\mathrm{tr}}(AA')\\ &=&\beta'(AX-I_p)(AX-I_p)\beta+ \sigma^2\mathop{\mathrm{tr}}(AA')=\mathop{\mathrm{tr}}(\mathop{\mathrm{EQM}}(\hat \beta)). \end{eqnarray*} \] En conclusion, que l’on s’intéresse à un estimateur linéaire de \(\beta\) ou à un estimateur linéaire de \(\gamma\), dès que l’on passe de l’un à l’autre en multipliant par \(Q\) ou \(Q'\), matrice orthogonale, la trace de l’EQM est identique, c’est-à-dire que les performances globales des 2 estimateurs sont identiques.
Nous avons le modèle de régression suivant~: \[ Z_{1:p}=\Delta\gamma+\eta_{1:p}, \] et donc, par définition de l’estimateur des MC, nous avons \[ \hat \gamma_{\mathrm{MC}}=(\Delta'\Delta)^{-1}\Delta'Z_{1:p}. \] Comme \(\Delta\) est une matrice diagonale, nous avons \[ \hat \gamma_{\mathrm{MC}}=\Delta^{-2}\Delta'Z_{1:p}=\Delta^{-1}Z_{1:p}. \] Cet estimateur est d’expression très simple et il est toujours défini de manière unique, ce qui n’est pas forcément le cas de \(\hat \beta_{\mathrm{MC}}\).

Comme \(Z_{1:p}\sim\mathcal{N}(\Delta\gamma,\sigma^2 I_p)\) nous avons que \(\hat \gamma_{\mathrm{MC}}=\Delta^{-1}Z_{1:p}\) suit une loi normale d’espérance \(\mathbf E(\Delta^{-1}Z_{1:p})=\Delta^{-1}\mathbf E(Z_{1:p})=\gamma\) et de variance \(\mathop{\mathrm{V}}(\hat \gamma_{\mathrm{MC}})=\sigma^2\Delta^{-2}\). Puisque \(\hat \gamma_{\mathrm{MC}}\) est un estimateur des MC, il est sans biais, ce qui est habituel.
L’EQM de \(\hat \gamma_{\mathrm{MC}}\), estimateur sans biais, est simplement sa variance. Pour la \(i^e\) coordonnée de \(\hat \gamma_{\mathrm{MC}}\), l’EQM est égal à l’élément \(i,i\) de la matrice de variance \(\mathop{\mathrm{V}}(\hat \gamma_{\mathrm{MC}})\), c’est-à-dire \(\sigma^2/\delta_i^2\). La trace de l’EQM est alors simplement la somme, sur toutes les coordonnées \(i\), de cet EQM obtenu.
Par définition \(\hat \gamma(c)=\mathop{\mathrm{diag}}(c_i)Z_{1:p}\) et nous savons que \(Z_{1:p}\sim\mathcal{N}(\Delta\gamma,\sigma^2 I_p).\) Nous obtenons que \(\hat \gamma(c)\) suit une loi normale d’espérance \(\mathbf E(\mathop{\mathrm{diag}}(c_i)Z_{1:p})=\mathop{\mathrm{diag}}(c_i)\Delta\gamma\) et de variance \[ \mathop{\mathrm{V}}(\hat \gamma(c))= \mathop{\mathrm{diag}}(c_i)\mathop{\mathrm{V}}(Z_{1:p})\mathop{\mathrm{diag}}(c_i)'= \sigma^2\mathop{\mathrm{diag}}(c_i^2). \] La loi de \(\hat \gamma(c)\) étant une loi normale de matrice de variance diagonale, nous en déduisons que les coordonnées de \(\hat \gamma(c)\) sont indépendantes entre elles.
Calculons l’EQM de la \(i^e\) coordonnée de \(\hat \gamma(c)\) \[ \mathop{\mathrm{EQM}}(\hat \gamma(c)_i)=\mathbf E(\hat \gamma(c)_i -\gamma)^2=\mathbf E(\hat \gamma(c)_i^2)+ \mathbf E(\gamma_i^2)-2\mathbf E(\hat \gamma(c)_i \gamma_i). \] Comme \(\gamma_i\) et que \(\mathbf E(\hat \gamma(c)_i^2)=\mathop{\mathrm{V}}(\hat \gamma(c)_i^2)+\{\mathbf E(\hat \gamma(c)_i^2)\}^2\), nous avons \[ \begin{align*} \mathop{\mathrm{EQM}}(\hat \gamma(c)_i)&=\sigma^2 c_i^2+(c_i\delta_i\gamma_i)^2+\gamma_i^2-2\gamma_i\mathbf E(\hat \gamma(c)_i)\\ &=\sigma^2 c_i^2+(c_i\delta_i\gamma_i)^2+\gamma_i^2-2\sigma^2 c_i\delta_i\gamma_i= \sigma^2c_i^2+\gamma_i^2(c_i\delta_i -1)^2. \end{align*} \]
De manière évidente si \(\gamma_i^2\) diminue, alors l’EQM de \(\hat \gamma(c)_i\) diminue aussi. Calculons la valeur de l’EQM quand \(\gamma_i^2=\frac{\sigma^2}{\delta_i^2}\frac{(1/\delta_i)+c_i}{(1/\delta_i)-c_i}\). Nous avons, grâce à la question précédente, \[ \begin{eqnarray*} \mathop{\mathrm{EQM}}(\hat \gamma(c)_i)&=&\sigma^2 c_i^2+(c_i\delta_i -1)^2\frac{\sigma^2}{\delta_i^2}\frac{(1/\delta_i)+c_i}{(1/\delta_i)-c_i}\\ &=&\sigma^2 c_i^2+\frac{\sigma^2}{\delta_i^2}(1 - c_i\delta_i)^2\frac{1+\delta_ic_i}{1-\delta_ic_i}\\ &=&\sigma^2 c_i^2+\frac{\sigma^2}{\delta_i^2}(1 - c_i\delta_i)(1+\delta_ic_i)\\ &=&\sigma^2 c_i^2+\frac{\sigma^2}{\delta_i^2}(1-\delta_i^2c_i^2)\\ &=&\sigma^2 c_i^2+\frac{\sigma^2}{\delta_i^2}-\sigma^2c_i^2=\frac{\sigma^2}{\delta_i^2}\\ &=&\mathop{\mathrm{EQM}}(\hat \gamma_{\mathrm{MC}}), \end{eqnarray*} \] d’où la conclusion demandée.
Par définition de \(\hat \gamma(c)\), nous avons \[ \begin{eqnarray*} \hat \gamma(c)&=&\mathop{\mathrm{diag}}(c_i)Z_{1:p}=\mathop{\mathrm{diag}}(\frac{\delta_i}{\delta_i^2+\kappa})Z_{1:p}\\ &=&(\Delta'\Delta + \kappa I_p)^{-1}\Delta'Z_{1:p}, \end{eqnarray*} \] puisque \(\Delta\) est diagonale. De plus nous avons \[ D = \bigl( \begin{smallmatrix} \Delta\\ 0 \end{smallmatrix}\bigr), \] ce qui entraîne que \(D'D=\Delta'\Delta\) et \(D'Z=\Delta' Z_{1:p}\). Nous obtenons donc \[ \hat \gamma(c)=(D'D+\kappa I_p)^{-1}D'Z. \] Rappelons que \(D=PXQ'\) avec \(P\) et \(Q\) matrices orthogonales, nous avons alors \[ \begin{eqnarray*} \hat \gamma(c)&=&(QX'P'PXQ' + \kappa I_p)^{-1} D'Z=(QX'XQ' + \kappa QQ')^{-1}D'Z\\ &=&(Q(X'X + \kappa I_p)Q')^{-1}D'Z=(Q')^{-1}(X'X + \kappa I_p)^{-1}(Q)^{-1}D'Z\\ &=&Q(X'X + \kappa I_p)^{-1}Q'D'Z. \end{eqnarray*} \] Comme \(Z=PY\) et \(D=PXQ'\), nous avons \[ \hat \gamma(c)=Q(X'X + \kappa I_p)^{-1}Q' QX'P' PY=Q(X'X + \kappa I_p)^{-1}XY. \] Enfin, nous savons que \(Q\hat\gamma=\hat \beta\), nous en déduisons que \(\hat\gamma=Q'\hat \beta\) et donc que dans le cas particulier où \(c_i=\frac{\delta_i}{\delta_i^2+\kappa}\) nous obtenons \[ \hat \beta=Q\hat \gamma(c)=(X'X + \kappa I_p)^{-1}XY, \] c’est-à-dire l’estimateur de la régression ridge.

Exercice 6 (Coefficient constant et régression sous contraintes)

Exercice 7 (Unicité pour la régression lasso, Giraud (2014))

Exercice 8

library(tidyverse)
signal <- read_csv("../donnees/courbe_lasso.csv")
donnees <- read_csv("../donnees/echan_lasso.csv")
ggplot(signal)+aes(x=x,y=y)+geom_line()+
  geom_point(data=donnees,aes(x=X,y=Y))

Nous cherchons à reconstruire le signal à partir de l’échantillon. Bien entendu, vu la forme du signal, un modèle linéaire de la forme \[ y_i=\beta_0+\beta_1x_i+\varepsilon_i \] n’est pas approprié. De nombreuses approches en traitement du signal proposent d’utiliser une base ou dictionnaire représentée par une collection de fonctions \(\{\psi_j(x)\}_{j=1,\dots,K}\) et de décomposer le signal dans cette base : \[ m(x)\approx \sum_{j=1}^K \beta_j\psi_j(x). \] Pour un dictionnaire donné, on peut alors considérer un modèle linéaire \[ y_i=\sum_{j=1}^K \beta_j\psi_j(x_i)+\varepsilon_i. \tag{1}\] Le problème est toujours d’estimer les paramètres \(\beta_j\) mais les variables sont maintenant définies par les éléments du dictionnaire. Il existe différents types de dictionnaire, dans cet exercice nous proposons de considérer la base de Fourier définie par \[ \psi_0(x)=1,\quad \psi_{2j-1}(x)=\cos(2j\pi x)\quad\text{et}\quad \psi_{2j}(x)=\sin(2j\pi x),\quad j=1,\dots,K. \]

mat.dict <- function(K,x){
    res <- matrix(0,nrow=length(x),ncol=2*K) |> as_tibble()
    for (j in 1:K){
      res[,2*j-1] <- cos(2*j*pi*x)
      res[,2*j] <- sin(2*j*pi*x)
    }
    return(res)
}

Il suffit d’ajuster le modèle linéaire où les variables explicatives sont données par le dictionnaire :

D25 <- mat.dict(25,donnees$X) |> mutate(Y=donnees$Y)
mod.lin <- lm(Y~.,data=D25)
S25 <- mat.dict(25,signal$x)
prev.MCO <- predict(mod.lin,newdata = S25)
signal1 <- signal |> mutate(MCO=prev.MCO) |> rename(signal=y)
signal2 <- signal1 |> pivot_longer(-x,names_to="meth",values_to="y")
ggplot(signal2)+aes(x=x,y=y)+geom_line(aes(color=meth))+
  scale_y_continuous(limits = c(-2,2))+geom_point(data=donnees,aes(x=X,y=Y))

Le signal estimé a tendance à surajuster les données. Cela vient du fait qu’on estime 51 paramètres avec seulement 60 observations.

On regarde tout d’abord le chemin de régularisation des estimateurs lasso

library(glmnet)
X.25 <- model.matrix(Y~.,data=D25)[,-1]
lasso1 <- glmnet(X.25,D25$Y,alpha=1)
plot(lasso1)

Il semble que quelques coefficients quittent la valeur 0 bien avant les autres. On effectue maintenant la validation croisée pour sélectionner le paramètre \(\lambda\).

set.seed(1234)
lasso.cv <- cv.glmnet(X.25,D25$Y,alpha=1)
plot(lasso.cv)

On calcule les prévisions et on trace le signal.

prev.lasso <- as.vector(predict(lasso.cv,newx=as.matrix(S25)))
signal1$lasso <- prev.lasso
signal2 <- signal1 |> pivot_longer(-x,names_to="meth",values_to="y")
ggplot(signal2)+aes(x=x,y=y)+geom_line(aes(color=meth))+
  scale_y_continuous(limits = c(-2,2))+geom_point(data=donnees,aes(x=X,y=Y))

L’algorithme lasso a permis de corriger le problème de sur-apprentissage. Les coefficients lasso non nuls sont les suivants

v.sel <- which(coef(lasso.cv)!=0)
v.sel

 [1]  1  2  4  5  6  8 21 28 30 38 40