Comment calculer le gradient pour un problème de moindres carrés non linéaires séparables?
Prenons le cas de la régression des moindres carrés non linéaire avec une variable dépendante $y_i$ et deux variables indépendantes $x_{i1}$ et $x_{i2}$ où la fonction non linéaire est une fonction linéaire de deux fonctions non linéaires $f_1$ et $f_2$ (pour simplifier je réduis cela à deux fonctions et fonctions avec un seul paramètre / coefficient mais cela peut être plus général)
$$y_i = \sum_{j=1,2} \alpha_j f_j(x_{ij},\beta_j) + \epsilon_i$$
Supposons que nous souhaitons utiliser ajuster cette fonction à certaines données avec une régression des moindres carrés, alors nous pourrions trouver la solution avec un algorithme par étapes qui alterne entre l'ajustement du $\alpha_j$ et $\beta_j$. Cela peut être une approche utile car la solution pour le$\alpha_j$ quand le $\beta_j$ sont fixes est facilement trouvée par la régression des moindres carrés ordinaires.
Pour effectuer l'étape d'optimisation pour le $\beta_j$nous avons besoin de connaître le gradient de la fonction de perte. Il existe des solveurs qui peuvent estimer les dérivées de manière informatique, mais les algorithmes seront plus rapides et plus précis lorsque nous pouvons fournir les dérivées nous-mêmes.
Comment décrivons-nous le dérivé $\frac{\partial L}{\partial \beta_j}$ de la fonction de perte de la somme des carrés des résidus $$L = \Vert y - \hat{y}\Vert ^2$$
quand
$$\hat y = F (F^T F)^{-1} F^T y$$
où le $F$ est la matrice des régresseurs $f(x_{ij}, \beta_{j})$
$$F = \begin{bmatrix} f(x_{{11}}, \beta_1) & f(x_{12}, \beta_2) \\ f(x_{{21}}, \beta_1) & f(x_{22}, \beta_2) \\ f(x_{{31}}, \beta_1) & f(x_{32}, \beta_2) \\ \vdots & \vdots \\ f(x_{{n1}}, \beta_1) & f(x_{n2}, \beta_2) \\ \end{bmatrix}$$
Il devrait y avoir un moyen simple d'exprimer
$$\frac{\partial L}{\partial \beta_j}$$
en terme de $\frac{\partial f(x_{ij})}{\partial \beta_j}$
Réponses
Une question connexe existe sur math.stackexchange.com Dérivée de projection par rapport à un paramètre:$D_{a}: X(a)[ X(a)^TX(a) ]^{-1}X(a)^Ty$
La réponse suggère d'utiliser la règle du produit qui conduit à:
$$\begin{align}\hat{y}^\prime =(X(X^TX)^{-1}X^Ty)^\prime&=X^\prime(X^TX)^{-1}X^Ty\\&-X(X^TX)^{-1}(X^{\prime T}X+X^TX^\prime)(X^TX)^{-1}X^Ty\\&+X(X^TX)^{-1}X^{\prime T}y\prime.\end{align}$$
Ensuite, nous calculons la dérivée de la fonction de perte comme
$$L^\prime = \left( \sum (y-\hat{y})^2 \right)^\prime = \sum -2(y-\hat{y})\hat{y}^\prime$$
Où $^\prime$ désigne le dérivé de l'un des $\beta_j$
Exemple:
Dans l'exemple ci-dessous, nous adaptons la fonction
$$y_i = \alpha_{1} e^{\beta_1 x_{1,i}} + \alpha_2 e^{\beta_2 x_{2,i}}$$
Dans ce cas $X^\prime = \frac{\partial}{\beta_j} X$ sera le même que $X$ mais avec le $i$-ème colonne multipliée par $x_i$ et les autres zéro.
Vous trouverez ci-dessous un code R qui illustre le calcul. Il s'agit d'une méthode de descente de gradient qui utilise la fonction frpour calculer la fonction de coût et la fonction grpour calculer le gradient. Dans cette fonction, grnous avons calculé les dérivées comme ci-dessus. La valeur de la fonction de coût en fonction de$\beta_1$ et $\beta_2$est illustré dans la figure ci-dessous. La ligne noire épaisse montre le chemin suivi par la méthode de descente de gradient.
set.seed(1)
# model some independent data t1 and t2
x1 <- runif(10,0,1)
x2 <- runif(10,0,0.1)+x1*0.9
t1 <- log(x1)
t2 <- log(x2)
# compute the dependent variable y according to the formula and some added noise
y <- round(1*exp(0.4*t1) - 0.5*exp(0.6*t2) + rnorm(10, 0 ,0.01),3)
###############################
# loss function
fr <- function(p) {
a <- p[1]
b <- p[2]
u1 <- exp(a*t1)
u2 <- exp(b*t2)
mod <- lm(y ~ 0 + u1 + u2)
ypred <- predict(mod)
sum((y-ypred)^2)
}
# gradient of loss function
gr <- function(p) {
a <- p[1]
b <- p[2]
u1 <- exp(a*t1) ### function f1
u2 <- exp(b*t2) ### function f2
X <- cbind(u1,u2) # matrix X
Xa <- cbind(t1*u1,0*u2) # derivative dX/da
Xb <- cbind(0*u1,t2*u2) # derivative dX/db
### predicted y
mod <- lm(y ~ 0 + u1 + u2)
ypred <- predict(mod)
### computation of the derivatives of the projection
dPa <- Xa %*% solve(t(X) %*% X) %*% t(X) %*% y -
X %*% solve(t(X) %*% X) %*% (t(Xa) %*% X + t(X) %*% Xa) %*% solve(t(X) %*% X) %*% t(X) %*% y +
X %*% solve(t(X) %*% X) %*% t(Xa) %*% y
dPb <- Xb %*% solve(t(X) %*% X) %*% t(X) %*% y -
X %*% solve(t(X) %*% X) %*% (t(Xb) %*% X + t(X) %*% Xb) %*% solve(t(X) %*% X) %*% t(X) %*% y +
X %*% solve(t(X) %*% X) %*% t(Xb) %*% y
### computation of the derivatives of the squared loss
dLa <- sum(-2*(y-ypred)*dPa)
dLb <- sum(-2*(y-ypred)*dPb)
### result
return(c(dLa,dLb))
}
# compute loss function on a grid
n=201
xc <- 0.9*seq(0,1.5,length.out=n)
yc <- 0.9*seq(0,1.5,length.out=n)
z <- matrix(rep(0,n^2),n)
for (i in 1:n) {
for(j in 1:n) {
z[i,j] <- fr(c(xc[i],yc[j]))
}
}
# levels for plotting
levels <- 10^seq(-4,1,0.5)
key <- seq(-4,1,0.5)
# colours for plotting
colours <- function(n) {hsv(c(seq(0.15,0.7,length.out=n),0),
c(seq(0.2,0.4,length.out=n),0),
c(seq(1,1,length.out=n),0.9))}
# empty plot
plot(-1000,-1000,
xlab=expression(n[1]),ylab = expression(n[2]),
xlim=range(xc),
ylim=range(yc)
)
# add contours
.filled.contour(xc,yc,z,
col=colours(length(levels)),
levels=levels)
contour(xc,yc,z,add=1, levels=levels, labels = key)
# compute path
# start value
new=c(0.9,1.1)
maxstep <- 0.001
# make lots of small steps
for (i in 1:5000) {
### safe old value
old <- new
### compute step direction by using gradient
grr <- -gr(new)
lg <- sqrt(grr[1]^2+grr[2]^2)
step <- grr/lg
### find best step size (yes this is a bit simplistic and computation intensive)
min <- fr(old)
stepsizes <- maxstep*10^seq(-2,0.001,length.out1=100)
for (j in stepsizes) {
if (fr(old+step*j)<min) {
new <- old+step*j
min <- fr(new)
}
}
### plot path
lines(c(old[1],new[1]),c(old[2],new[2]),lw=2)
}
# finish plot with title and annotation
title(expression(paste("Solving \n", sum((alpha[1]*e^{beta[1]*x[i,1]}+alpha[2]*e^{beta[2]*x[i,2]}-y[i])^2,i==1,n))))
points(0.9,1.1)
text(0.9,1.1,"start",pos=2,cex=1)
points(new[1],new[2])
text(new[1],new[2],"end",pos=4,cex=1)
Voir pour une vitrine historique de cette méthode:
«La différenciation des pseudo-inverses et des problèmes de moindres carrés non linéaires dont les variables sont séparées» par GH Golub et V. Pereyra dans SIAM Journal on Numerical Analysis Vol. 10, n ° 2 (1973), pp. 413-432