Jak obliczyć gradient dla rozdzielnego nieliniowego problemu najmniejszych kwadratów?
Rozważmy przypadek nieliniowej regresji metodą najmniejszych kwadratów z jedną zmienną zależną $y_i$ i dwie zmienne niezależne $x_{i1}$ i $x_{i2}$ gdzie funkcja nieliniowa jest funkcją liniową dwóch funkcji nieliniowych $f_1$ i $f_2$ (dla uproszczenia sprowadzam to do dwóch funkcji i funkcji z tylko jednym parametrem / współczynnikiem, ale może być bardziej ogólne)
$$y_i = \sum_{j=1,2} \alpha_j f_j(x_{ij},\beta_j) + \epsilon_i$$
Powiedzmy, że chcemy użyć funkcji dopasowania tej funkcji do niektórych danych z regresją najmniejszych kwadratów, a następnie moglibyśmy znaleźć rozwiązanie za pomocą algorytmu krokowego, który zmienia się między dopasowywaniem $\alpha_j$ i $\beta_j$. Może to być przydatne podejście, ponieważ rozwiązanie dla$\alpha_j$ kiedy $\beta_j$ są stałe, można łatwo znaleźć za pomocą zwykłej regresji metodą najmniejszych kwadratów.
Aby wykonać krok optymalizacji dla $\beta_j$musimy znać gradient funkcji straty. Istnieją solwery, które potrafią oszacować pochodne obliczeniowo, ale algorytmy będą szybsze i dokładniejsze, gdy będziemy mogli sami zapewnić pochodne.
Jak opiszemy pochodną $\frac{\partial L}{\partial \beta_j}$ funkcji straty sumy kwadratów reszt $$L = \Vert y - \hat{y}\Vert ^2$$
gdy
$$\hat y = F (F^T F)^{-1} F^T y$$
gdzie $F$ jest macierzą regresorów $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}$$
Powinien istnieć jakiś prosty sposób wyrażenia
$$\frac{\partial L}{\partial \beta_j}$$
pod względem $\frac{\partial f(x_{ij})}{\partial \beta_j}$
Odpowiedzi
Podobne pytanie istnieje na math.stackexchange.com Pochodna projekcji w odniesieniu do parametru:$D_{a}: X(a)[ X(a)^TX(a) ]^{-1}X(a)^Ty$
Odpowiedź sugeruje zastosowanie reguły iloczynu, która prowadzi do:
$$\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}$$
Następnie obliczamy pochodną funkcji straty jako
$$L^\prime = \left( \sum (y-\hat{y})^2 \right)^\prime = \sum -2(y-\hat{y})\hat{y}^\prime$$
Gdzie $^\prime$ oznacza pochodną dowolnego z $\beta_j$
Przykład:
W poniższym przykładzie dopasowujemy funkcję
$$y_i = \alpha_{1} e^{\beta_1 x_{1,i}} + \alpha_2 e^{\beta_2 x_{2,i}}$$
W tym przypadku $X^\prime = \frac{\partial}{\beta_j} X$ będzie taki sam jak $X$ ale z $i$-ta kolumna pomnożona przez $x_i$ a pozostałe zero.
Poniżej znajduje się kod R, który ilustruje obliczenia. Jest to metoda zstępowania gradientu, która używa funkcji frdo obliczenia funkcji kosztu i funkcji grdo obliczenia gradientu. W tej funkcji grobliczyliśmy pochodne jak powyżej. Wartość funkcji kosztu w funkcji$\beta_1$ i $\beta_2$pokazano na poniższym rysunku. Gruba czarna linia pokazuje ścieżkę, po której następuje metoda gradientu spadku.
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)
Zobacz historyczną wizytówkę tej metody:
„The Differentiation of Pseudo-Inverses and Nonlinear Least Squares Problems Whose Variables Separate” autorstwa GH Goluba i V. Pereyry w SIAM Journal on Numerical Analysis Vol. 10, nr 2 (1973), str. 413-432