Descente de gradient pour la régression linéaire

Partie 1 : Génération d’un modèle linéaire

# code reproductible pour les fois où l'on fait appel à un processus aléatoire
set.seed(123456789) 

# Paramétres du modèle.
p = 1000    #nombre de paramètres de régression
n = p * 10  #nombre d'observations.

# Simulation d'un modèle linéaire. 
# Matrice de design X.

X <- cbind(rep(1, n),matrix(data=rnorm(n * p, 0, 1), nrow=n, ncol=p)) 

# Vecteur de paramètres.
theta <- matrix(data=1:(p+1), nrow=p+1, ncol=1)

# Résidus e.
e <- matrix(data=rnorm(n, 0, 1), nrow=n, ncol=1)   
# Observations Y.
Y <-X %*% theta + e 

Elle est donnée par :

$$J(\theta) = \frac{1}{n} \vert Y - X\theta \vert{_2}^2$$

# Régression linéaire
reg <- lm(Y ~ X )
summ  <-summary(reg)

# Erreur quadratique moyenne
EQM <- mean((summ$residuals^2))
EQM

On trouve :

## [1]  0.8960279

L’estimateur des moindres carrés $\hat{\theta}$ de $\theta$ est solution du système d’équations :

$$X^TX\theta = X^TY$$

# Marquer le temps de début
t1 <- Sys.time()

# Résolution du systéme linéaire.
theta_chap <- solve(t(X) %*% X , t(X) %*% Y)

# Marquer le temps de fin
t2 <- Sys.time()

# t2-t1 correspond au temps d'exécution de la commande solve
tdiff = difftime(t2, t1)
tdiff

Le temps d’exécution vaut approximativement 9 secondes.

## Time difference of 8.958424 secs

Elle est donnée par :

$$J(\hat{\theta}) = \frac{1}{n} \vert Y - X\hat{\theta} \vert{_2}^2$$

# calcul de J_theta_chap
J_theta_chap <- norm((Y-X %*% theta_chap), type = "2")/n
J_theta_chap

## [1] 0.009449069

L’erreur d’estimation est définie par :

$$ E(\hat{\theta}) = \vert \hat{\theta} - \theta \vert_1 = \sum_{i=1}^{p+1}\vert \hat{\theta}_i - \theta \vert $$

# Erreur d'estimation
diff <- norm(theta_chap - theta, type = "1") #par défaut type= 1 (norme 1)
diff

## [1] 8.235012

Partie 2 : Descente de gradient à pas fixe

On se propose de mettre en oeuvre un algorithme de descente de gradient à pas fixe pour approcher le minimum global de la fonction $J$.

On rappelle que le gradient de $J$ en $\theta$ est :

$$\nabla J(\theta)= X^T (X\theta - Y)$$

La méthode du gradient à pas fixe est décrite par l’algorithme suivant :

• En entrée, nous avons : les observations $(Y)$, la matrice de design $(X)$, le pas de gradient $(\rho)$, le nombre maximum d'itérations $(N_{max})$ et la tolérance $(\epsilon)$;
• En sortie, nous avons : l'estimation des paramètres $(\theta)$, les valeurs successives de $J$ $(F)$;
• Initialisation aléatoire de $\theta$;
• Calcul de $F(1)=J(\theta)$;
• Pour $i$ allant de $2$ à $(N_{max})$ :
• Calcul du gradient de $(J)$: $$d = \nabla J(\theta).$$
• Si $\vert d \vert_1 < \epsilon$, Arrêt des iterations.
• Sinon, Mise à jour de $\theta$. $$\theta \leftarrow \theta - \rho d.$$
• Calcul de $F(i)=J(\theta)$.

t21 <- Sys.time()
N_max = 1000
eps = 0.01
p= 0.0001

theta <- matrix(data = rnorm(1001, 0,1),
                nrow = 1001,
                ncol=1)
F_1 <- norm((Y- X%*% theta),
            type ="2") / n
F_1

## [1] 183.5106

for (i in (2:N_max)){
  # calcul du gradient.
  d = t(X) %*% (X %*% theta - Y) 
  if (norm(d) < eps)
    {break}
  else{
    # mise a jour des parametres.
     theta <- theta - p * d
    }
  # calcul de la valeur de la fonction F.
  F_1 <- norm((Y- X %*% theta), type="2") / n
  print(F_1)
}
t22 <- Sys.time()

# pour un pas de 0.001 F(i) =  0.009465875, on n'a pas convergence : l'algo diverge

# pour un pas de 0.00001 F(i) =  0.009465875 (même résulat avec i = 491)

# pour un pas de 0.0001 F(i) =  0.009465875, i = 76 et le temps d'exécution est de 8.994529 secs


# j(theta) = 0.009465875 = F(76)

# nombre d'iterations i = 76
i

Le nombre d’itérations est égal à 76.

## [1] 76

# Temps d'exécution
temps_dexecution <- difftime(t22, t21)
temps_dexecution

Le temps d’exécution vaut 8 secondes environ.

## Time difference of 8.237663 secs

# On trouve presque la même chose j_theta_chap = 0.009465875

Partie 3 : Descente de gradient à pas optimal

La méthode du gradient à pas optimal est décrite par l’algorithme suivant :

• En entrée, nous avons : les observations $(Y)$, la matrice de design $(X)$, le pas de gradient $(\rho)$, le nombre maximum d'itérations $(N_{max})$ et la tolérance $(\epsilon)$;
• En sortie, nous avons : l'estimation des paramètres $(\theta)$, les valeurs successives de $J$ $(F)$;
• Définition de la matrice $A=X^TX$;
• Initialisation aléatoire de $\theta$;
• Calcul de $F(1)=J(\theta)$;
• Pour $i$ allant de $2$ à $(N_{max})$ :
• Calcul du gradient de $(J)$ : $$d = \nabla J(\theta).$$
• Si $\vert d \vert_1 < \epsilon$, Arrêt des iterations.
• Sinon,
  • Calcul du pas : $$\rho = \frac{\vert d \vert_2}{\langle Ad, d\rangle}.$$
  • Mise à jour de $\theta$ : $$\theta \leftarrow \theta - \rho d.$$
• Calcul de $F(i)=J(\theta)$.

t31 <- Sys.time()
N_max = 1000
eps = 0.01
p= 0.0001
A= t(X) %*% X
t32 <- Sys.time()
theta <- matrix(data = rnorm(1001, 0,1),
             nrow = 1001,
             ncol=1)

F_1 <- norm((Y- X%*% theta),
            type ="2") / n
F_1

## [1] 183.5064

for (i in (2:N_max)){
  # calcul du gradient.
  d = t(X) %*% (X %*% theta - Y) 
  if (norm(d) < eps)
    {break}
  else{
    # calcul jour du pas.
    p = norm(d, type="2")^2 / sum((A %*% d) * d)
    # mise a jour des parametres.
    theta <- theta - p * d
    }
  # calcul de la valeur de la fonction.
  F_1 <- norm((Y- X %*% theta), type="2") / n
  print(F_1)
}
t33 <- Sys.time()

F_1

## [1] 0.009449069

Le nombre d’itérations qu’il a fallu est égal à 45.

# nombre d'iteration i = 45
## [1] 45

Calcul du temps d’exécution de $J(\theta)$

# Temps d'exécution de tout le script
temps_dexecution <- difftime(t33, t31)
temps_dexecution

## Time difference of 13.55263 secs

Le temps d’exécution de $J(\theta)$ vaut 5 secondes environ.

# temps d'exécution de J_theta
tmp <- difftime(t33, t32)
tmp

## Time difference of 4.879829 secs

# On trouve presque la même chose i.e J_theta_chap = 0.009465875

Conclusion : Nous trouvons quasiment la même valeur pour $J(\hat{\theta})$ que ce soit avec la résolution du système linéaire ou avec la méthode de descente de gradient à pas fixe ou à pas optimal.

Outil utilisé : RStudio, latex

License

Copyright 2020-present Mamoudou KOUME.