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.
Soit une fonction ƒ à valeurs réelles, définie sur un espace vectoriel de dimension "n" :
On recherche le minimum de ƒ, c'est-à-dire un vecteur ("x"1*, "x"2*, …, "xn"*) 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 codice_1 ou codice_2 (synonymes). Sa syntaxe globale est :
optimize(fct, intervalle)
optimize(fct, intervalle, maximum = TRUE)
optimize(f = fct, interval = intervalle)
où
Si par exemple on désire connaître le minimum de la fonction 
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 , qui est une approximation de -π/2 ≃ , et la fonction y vaut -1. Le maximum est en , 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 codice_4. Par exemple, pour avoir une précision de 10-6 :
optimize(f, c(-pi, pi), tol = 1e-6)
Nous disposons également de la commande codice_5 ; 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ù
Si par exemple on désire connaître le minimum de la fonction 
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 :
La commande permet donc de travailler avec une fonction de plusieurs variables, par exemple si l'on considère la fonction
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
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 :
Le logiciel R propose la commande codice_9, qui fait appel au solveur codice_5 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ù
Par exemple, si l'on veut rechercher le minimum de la fonction "y" = sin("x") sur l'intervalle [- ; ] :
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 ≃ - ; 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'codice_5 :
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 codice_2.
La fonction codice_2 utilise la méthode du nombre d'or associée à une interpolation parabolique.
Fonctions codice_5 et codice_9.
Méthode de Nelder-Mead.
Par défaut, la fonction codice_5 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 codice_2.
Méthode BFGS.
Le deuxième algorithme proposé est celui de BFGS, avec l'option codice_23. 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 codice_24 ; 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")
L'option codice_25 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, codice_26, 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 codice_27. 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 codice_25.
Méthode du recuit simulé.
Le quatrième algorithme proposé est celui du recuit simulé "(simulated annealing)", avec l'option codice_29. 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 codice_30. Ce solveur ne traite que des problèmes à une dimension, et est en fait similaire à la fonction codice_2. En particulier, il faut définir l'intervalle de recherche avec les options codice_32 et codice_33, 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 codice_9 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
alors R prend pour fonction de pénalité
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 codice_35 permet d'utiliser d'autres distributions, avec le paramètre codice_36.
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" (β0, β1, …, β"m") ; c'est l'expression de ces coefficients qui doit être linéaire, du type :
où x est un vecteur ("x"1, "x"2, …, "xn"). 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 :
ce modèle n'est pas linéaire en β1 — 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 codice_37 "(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 codice_38, mais il faut introduire explicitement les paramètres à affiner dans la formule. Par exemple :
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
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 codice_39, qui doit recevoir une liste nommée. par exemple, pour commencer avec β0 = β1 = 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 codice_40, on utilise un modèle affine (droite avec ordonnée à l'origine) du type
Pour une régression linéaire avec deux variables explicatives, on peut écrire codice_41 ; le signe plus n'indique pas ici l'addition, mais sépare les deux paramètres, le modèle sera donc
Attention : les variables du modèle — codice_42, codice_43, codice_44, codice_45 — 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 codice_46, par exemple codice_47.
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 codice_48, par exemple :
formula = I(x^2 + y^2) ~ x + y
indique un modèle