Programmer en R/Optimiser une fonction

Les méthodes d'optimisation sont des méthodes numériques de recherche du minimum ou d'un maximum d'une fonction.

Programmer en R

• Introduction
  • Faire ses premiers pas
  • Installer R
  • Choisir son environnement de travail
  • Installer des librairies
  • Se documenter
• Manipulation et mise en forme des données
  • Manipuler les nombres
  • Manipuler les vecteurs
  • Manipuler les chaînes de caractères
  • Manipuler les facteurs
  • Manipuler les dates et les heures
  • Manipuler un tableau de données
  • Manipuler une liste
• Importation et exportation des données
  • Importer un fichier CSV
  • Résoudre les problèmes d'encodage des caractères
  • Importer un fichier TSV
  • Importer un fichier texte à largeur fixe
  • Importer des données depuis un tableur
  • Import un fichier JSON
  • Extraire des données d'une page web
• Visualisation de données
  • Comprendre la grammaire des graphiques
  • Représenter une carte
  • Visualiser une carte choroplèthe
• Modélisation statistique
  • Estimer un modèle de régression linéaire
  • Estimer un modèle de régression circulaire
• Communication
• Programmation avancée
  • Programmer une boucle
  • Programmer une fonction
• Outils mathématiques
  • Optimiser une fonction
  • Manipuler les matrices
  • Les fonctions statistiques de base
  • Les loi de probabilités, ajustement et test

Soit une fonction ƒ à valeurs réelles, définie sur un espace vectoriel de dimension n :

y = ƒ(x₁, x₂, …, x_n)

On recherche le minimum de ƒ, c'est-à-dire un vecteur (x₁*, x₂*, …, x_n*) tel que la valeur de ƒ soit la plus petite possible.

Les méthodes de régression consistent à ajuster un modèle à des données. Elles consistent à minimiser une fonction d'écart, à l'optimiser.

Méthode d'optimisation sans contrainte

Fonction sans paramètre

Dans le cas d'une optimisation sans contrainte, R propose la commande optimise() ou optimize() (synonymes). Sa syntaxe globale est :

optimize(fct, intervalle)
optimize(fct, intervalle, maximum = TRUE)
optimize(f = fct, interval = intervalle)

où

• fct est la fonction à minimiser, voir Programmation procédurale > Fonctions définies par l'utilisateur ;
• intervalle est l'intervalle sur lequel il effectue la recherche ;
• maximum = est un booléen : « vrai » si l'on cherche la maximum, « faux » (valeur par défaut) si l'on cherche le minimum.

Si par exemple on désire connaître le minimum de la fonction

ƒ(x) = sin(x)

sur l'intervalle [-π ; π], on peut écrire

f <- function(x) sin(x)

minf <- optimize(f, c(-pi, pi))
maxf <- optimize(f, c(-pi, pi), maximum = TRUE)

print(minf$minimum)
print(minf$objective)

print(maxf$maximum)
print(maxf$objective)

ce qui nous donne

[1] -1.570798
[1] -1
[1] 1.570798
[1] 1

Donc, le minimum est en −1,570,798, qui est une approximation de -π/2 ≃ −1,570,796, et la fonction y vaut -1. Le maximum est en 1,570,798, qui est une approximation de π/2, et la fonction y vaut 1.

On peut aussi utiliser la syntaxe

minf <- optimize(f, lower = -pi, upper = pi)
maxf <- optimize(f, lower = -pi, upper = pi, maximum = TRUE)

On peut lui indiquer la tolérance visée sur le résultat, avec l'option tol =. Par exemple, pour avoir une précision de 10^-6 :

optimize(f, c(-pi, pi), tol = 1e-6)

Nous disposons également de la commande optim() ; elle ne cherche que le minimum, elle ne permet pas par défaut de définir un intervalle, et il faut lui donner un point de départ ; par contre, elle permet de travailler avec des fonctions de plusieurs variables. Sa syntaxe globale est :

optim(x0, f)

où

• x₀ est un vecteur de valeurs initiales ;
• f est la fonction à minimiser.

Si par exemple on désire connaître le minimum de la fonction

ƒ(x) = 5x² + 3x - 4

on peut écrire (en utilisant la méthode de Horner pour calculer le polynôme)

f <- function(x) (5*x + 3)*x - 4
x0 = 0
solution = optim(x0, f)

Le solveur nous donne un message d'avertissement qui est sans conséquence.

L'abscisse de la solution est donnée par

> solution$par
[1] -0.3

qui est la solution exacte, et la valeur de la fonction en ce point est

> solution$value
[1] -4.45

Si l'on affiche la totalité de la solution, on obtient

> solution
$par
[1] -0.3

$value
[1] -4.45

$counts
function gradient 
      26       NA 

$convergence
[1] 0

$message
NULL

qui nous indique :

• $count : le nombre de fois (26) qu'il a fallu appeler la fonction f ; c'est globalement le nombre d'itérations qui a été nécessaire pour trouver la solution ;
• $convergence : la qualité de la convergence ; « 0 » indique que la méthode a bien convergé ;
• $message : un message éventuel du solveur (aucun).

La commande permet donc de travailler avec une fonction de plusieurs variables, par exemple si l'on considère la fonction

ƒ(x, y) = x² + y² - 10x - 8y + 41

on peut trouver son minimum avec

f <- function(x) x[1]^2 + x[2]^2 - 10*x[1] - 8*x[2] + 16
x0 <- c(0, 0)

solution <- optim(x0, f)

print(solution$par)
print(solution$value)

ce qui donne

[1] 4.999483 4.000116
[1] -25

La solution exacte est x = 5 et y = 4.

Fonction avec paramètre

Si la fonction possède des paramètres, on doit les passer en utilisant leur nom exact. Par exemple, si la fonction est

ƒ_{a, b, c}(x, y) = (x - a)² + (y - b)² + c

alors le script suivant permet de trouver le minimum :

f <- function(x, a, b, c) (x[1] - a)^2 + (x[2] - b)^2 + c
x0 <- c(0, 0)

solution <- optim(x0, f, a = 5, b = 4, c = -25)

print(solution$par)
print(solution$value)

ce qui donne

[1] 4.999483 4.000116
[1] -25

Méthode d'optimisation avec contrainte linéaire

On recherche le minimum de ƒ, mais avec des contraintes linéaires, de type :

inégalités : u_i×x_i ≥ c_i

Le logiciel R propose la commande constrOptim(), qui fait appel au solveur optim() mais en ajoutant une fonction de pénalité pour gérer les contraintes. La syntaxe de base est :

constrOptim(x0, f, NULL, u, c)
constrOptim(par = x0, fn = f, ui = u, ci = c)

où

• x₀ est le vecteur de départ, qui doit être dans la zone de faisabilité, c'est-à-dire vérifier les inégalités de contrainte ;
• f est la fonction ;
• NULL indique que l'on ne définit pas la fonction gradient ;
• u et c sont des vecteurs représentant les contraintes.

Par exemple, si l'on veut rechercher le minimum de la fonction y = sin(x) sur l'intervalle [-π ; π] :

• la fonction f est sin
• les contraintes linéaires sont
  1×x ≥ -π,
  -1×x ≥ -π
  donc on a en notation matricielle u×x =c avec
  $\mathbf {u} ={\begin{pmatrix}1\\-1\end{pmatrix}}$ et $\mathbf {c} ={\begin{pmatrix}-\pi \\-\pi \end{pmatrix}}$ ;
• on part de x₀ = 0.

x0 <- matrix(0, 1, 1)
f <- sin
u <- matrix(c(1, -1), 2, 1)
c <- matrix(c(-pi, -pi), 2, 1)
solution <- constrOptim(x0, f, NULL, u, c)

La résolution provoque des messages d'avertissement, qui sont sans conséquence ici. L'abscisse de la solution est donnée par

> solution$par
         [,1]
[1,] -1.57066

qui est une approximation de -π/2 ≃ -1,570,8 ; la valeur de la fonction en ce point est donnée par

> solution$value
     [,1]
[1,]   -1

Si l'on affiche la totalité de la solution, on obtient

$par
         [,1]
[1,] -1.57066

$value
     [,1]
[1,]   -1

$counts
function gradient 
      62       NA 

$convergence
[1] 0

$message
NULL

$outer.iterations
[1] 3

$barrier.value
              [,1]
[1,] -0.0008014317

qui nous indique, en plus des données d'optim() :

• $outer.iterations : le nombre d'itération « extérieures », c'est-à-dire ici d'appels au solveur optim() ;
• $barrier.value : c'est, lors de la dernière itération, la valeur du paramètre de pénalisation qui permet d'éviter de sortir des limites imposées.

Les différents algorithmes d'optimisation

Jusqu'ici, nous avons utilisé les méthodes d'optimisation par défaut. Selon les problèmes que l'on aborde, on peut avoir intérêt à utiliser d'autres algorithmes.

Fonction optimize()

La fonction optimize() utilise la méthode du nombre d'or associée à une interpolation parabolique.

Fonctions optim() et constrOptim()

Méthode de Nelder-Mead

Par défaut, la fonction optim() utilise la méthode de Nelder-Mead. Il consiste à prendre un simplexe — deux points sur une courbe, trois points sur une surface, … — et d'étendre ce simplexe vers les directions où la fonction est plus petite, puis à réduire ce simplexe lorsque le minimum est à l'intérieur.

Cet algorithme fonctionne lorsque la fonction n'est pas dérivable (mais elle doit être continue). Par contre, il n'est pas adapté aux problèmes à une dimension (fonction d'une seule variable), il vaut alors mieux utiliser la méthode de Brent ou bien la fonction optimize().

Méthode BFGS

Le deuxième algorithme proposé est celui de BFGS, avec l'option method = "BFGS". C'est un algorithme à direction de descente : à partir d'une point de départ, on se déplace sur le plan tangent à la surface, vers le bas. La direction de descente est l'opposé du gradient de la fonction.

Le pas de descente est fonction l'inverse de la courbure : si la courbure est importante, alors la surface s'écarte vite du plan tangent, il ne faut donc pas aller trop loin. L'algorithme calcule donc une approximation de l'inverse de cette dérivée seconde, ou plutôt, dans le cas d'une fonction de plusieurs variables, de la matrice hessienne (contenant les dérivées partielles d'ordre 2).

On peut fournir la fonction gradient, avec le troisième paramètre gr = … ; sinon, le solveur utilise une approximation par différences finies. Par exemple

f <- function(x) sin(x)
x0 = 0

solution <- optim(x0, f, method = "BFGS")

ou bien

f <- function(x) sin(x)
g <- function(x) cos(x)
x0 = 0

solution <- optim(x0, f, g, method = "BFGS")
# syntaxe alternative : optim(x0, f, gr = g, method = "BFGS")

L'option hessian = TRUE retourne également la matrice hessienne calculée à la solution. On peut aussi calculer la matrice hessienne par

optimHess(solution$par, f)

où solution$par est la solution trouvée par optimisation.

Le solveur propose une méthode BFGS modifiée, method = "L-BFGS-B", qui permet de définir des limites inférieures et supérieures aux variables.

optim(x0, f, lower = …, upper = …, method = "L-BFGS-B")

Méthode du gradient conjugué

Le troisième algorithme proposé est celui du gradient conjugué, avec l'option method = "CG". Cette méthode ne calcule pas, et donc ne stocke pas, les matrices hessiennes ; elle permet donc de traiter des problèmes plus volumineux que la méthode BFGS, mais est moins robuste.

Le solveur accepte l'option hessian = TRUE.

Méthode du recuit simulé

Le quatrième algorithme proposé est celui du recuit simulé (simulated annealing), avec l'option method = "SANN". Le but de cette méthode est d'éviter de trouver un minimum local. Pour cela, une fois qu'une solution est trouvée, le solveur repart d'un autre point pour voir s'il converge vers une autre valeur. Cette méthode est utile lorsqu'il y a de nombreux minima locaux (surface « rugueuse »).

Méthode de Brent

Le quatrième algorithme proposé est la méthode de Brent, avec l'option method = "Brent". Ce solveur ne traite que des problèmes à une dimension, et est en fait similaire à la fonction optimize(). En particulier, il faut définir l'intervalle de recherche avec les options lower = … et upper = …, mais il faut par contre toujours donner un point de départ.

optim(x0, f, lower = …, upper = …, method = "Brent")

Méthode de la barrière logarithmique

Pour la gestion des contrainte, la fonction constrOptim() utilise la méthode de la barrière logarithmique. Cette méthode consiste à ajouter à la fonction ƒ une fonction de pénalité qui augmente lorsque l'on s'approche de la frontière. Si l'inégalité de contrainte s'écrit sous la forme

h(x) ≥ 0

alors R prend pour fonction de pénalité

μ×ln(h(x))

où μ est une constante positive. Cette constante vaut par défaut 10<spu>-4, mais peut être ajustée :

constrOptim(x0, f, NULL, u, c, mu = …)

Une méthode consiste à faire des passes successives en diminuant au fur et à mesure la valeur de μ pour affiner le résultat, la fonction optimisée étant de plus en plus proche de la fonction initiale.

Modèle linéaire généralisé

La régression des moindres carrés consiste à minimiser l'estimateur de la variance ; cette démarche est pertinente si le bruit suit une distribution normale ayant le même écart type. La commande glm() permet d'utiliser d'autres distributions, avec le paramètre family.

Régression non-linéaire

Notons tout d'abord qu'une régression linéaire est un modèle dans lequel on recherche des coefficients (β₀, β₁, …, β_m) ; c'est l'expression de ces coefficients qui doit être linéaire, du type :

y = β₀ + β₁⋅ƒ₁(x) + β₂⋅ƒ₂(x) + … + β_m⋅ƒ_m(x)

où x est un vecteur (x₁, x₂, …, x_n). Les fonctions ƒ_i ne sont pas nécessairement linéaires.

Un modèle non-linéaire est un modèle dans lequel les coefficients β_i figurent eux-même de manière non-linéaire dans une formule. Par exemple, si l'on recherche un modèle gaussien :

$y=\beta _{0}\cdot \exp \left(-\beta _{1}\cdot x^{2}\right)$

ce modèle n'est pas linéaire en β₁ — notons qu'on peut le rendre linéaire en considérant le logarithme de y.

La régression non-linéaire par la méthode des moindres carrés se fait avec la commande nls() (nonlinear least squares). Par défaut, cette commande utilise la méthode de Gauss-Newton ; dans tous les cas, il s'agit d'algorithmes itératifs. La syntaxe de base est similaire à celle de lm(), mais il faut introduire explicitement les paramètres à affiner dans la formule. Par exemple :

# ********************
# Création des données

u1 <- seq(-1, 1, 0.01) # abscisse, variable explicative
u2 <- 2*exp(-3*u1^2) + rnorm(u1, 0, 0.1)  # ordonnée, variable expliquée

# ******************************
# régression des moindres carrés
donnees <- data.frame(u1, u2)
modele <- nls(u2~b0*exp(-b1*u1^2), donnees)

print("Modèle des moindres carrés : ")
print(modele)
plot(u1, u2)
lines(u1, fitted(modele), lwd = 2)

donne pour résultat :

[1] "Modèle des moindres carrés :"
Nonlinear regression model
  model:  u2 ~ b0 * exp(-b1 * u1^2) 
   data:  donnees 
   b0    b1 
1.975 2.965 
 residual sum-of-squares: 1.932

Number of iterations to convergence: 4 
Achieved convergence tolerance: 5.148e-07

La valeur de départ de l'algorithme itératif est choisie automatiquement. On peut lui indiquer des valeurs de départ avec le paramètre start, qui doit recevoir une liste nommée. par exemple, pour commencer avec β₀ = β₁ = 1 :

modele <- nls(u2~b0*exp(-b1*u1^2), data = donnees, start = list(b0 = 1, b1 = 1))

Annexe : Formules symboliques

Le modèle de régression est une formule symbolique. Par exemple, si l'on écrit y ~ x, on utilise un modèle affine (droite avec ordonnée à l'origine) du type

y = a + bx.

Pour une régression linéaire avec deux variables explicatives, on peut écrire y ~ x1 + x2 ; le signe plus n'indique pas ici l'addition, mais sépare les deux paramètres, le modèle sera donc

y = a + b₁x₁ + b₁x₂.

Attention : les variables du modèle — x, x1, x2, y — doivent être les noms des variables utilisées dans la structure de données.

Si l'on veut un modèle linéaire sans ordonnée à l'origine, on modifie la formule en ajoutant un paramètre +0, par exemple y ~ x + 0.

Dans une formule symbolique, le signe plus indique que l'on combine linéairement deux variables. Si l'on veut appliquer une addition stricte, c'est-à-dire additionner deux termes sans les assortir d'un coefficient à définir par régression, il faut « inhiber l'interpréteur » avec la fonction I(), par exemple :

formula = I(x^2 + y^2) ~ x + y

indique un modèle

x² + y² = a + b₁x + b₂y.

Formules symboliques dans R

Syntaxe	Modèle	Commentaires
Y ~ X	Y = β0 + β1X
Y ~ 0 + X	Y = β1X
Y ~ X + I(X^2)	Y = β0 + β1X + β2X2
Y ~ X1 + X2	Y = β0 + β1X1 + β2X2	modèle du premier ordre sans interaction croisée
Y ~ .	Y = β0 + β1X1 + …	le point remplace toutes les colonnes des données (data.frame) non encore mentionnées dans la formule
Y ~ X1:X2	Y = β0 + β1X1X2	modèle ne contenant que des interactions croisée du premier ordre
Y ~ X1*X2	Y = β0 + β1X1 + β2X2 + β3X1X2	modèle du premier ordre complet équivalent de Y ~ X1 + X2 + X1:X2