Wie berechnet man den Gradienten für ein trennbares nichtlineares Problem der kleinsten Quadrate?

Question

Betrachten Sie den Fall einer nichtlinearen Regression der kleinsten Quadrate mit einer abhängigen Variablen $y_i$ und zwei unabhängige Variablen $x_{i1}$ und $x_{i2}$ wobei die nichtlineare Funktion eine lineare Funktion zweier nichtlinearer Funktionen ist $f_1$ und $f_2$ (Der Einfachheit halber reduziere ich dies auf zwei Funktionen und Funktionen mit nur einem Parameter / Koeffizienten, aber es kann allgemeiner sein)

$$y_i = \sum_{j=1,2} \alpha_j f_j(x_{ij},\beta_j) + \epsilon_i$$

Angenommen, wir möchten diese Funktion an einige Daten mit Regression der kleinsten Quadrate anpassen, dann könnten wir die Lösung mit einem schrittweisen Algorithmus finden, der zwischen dem Anpassen der $\alpha_j$ und $\beta_j$. Dies kann ein nützlicher Ansatz sein, da die Lösung für die$\alpha_j$ wenn die $\beta_j$ fest sind, wird leicht durch gewöhnliche Regression der kleinsten Quadrate gefunden.

Um den Optimierungsschritt für die durchzuführen $\beta_j$Wir müssen den Gradienten der Verlustfunktion kennen. Es gibt Löser, die die Ableitungen rechnerisch schätzen können, aber die Algorithmen sind schneller und genauer, wenn wir die Ableitungen selbst bereitstellen können.

Wie beschreiben wir die Ableitung? $\frac{\partial L}{\partial \beta_j}$ der Summe der quadratischen Residuenverlustfunktion $$L = \Vert y - \hat{y}\Vert ^2$$

wann

$$\hat y = F (F^T F)^{-1} F^T y$$

bei dem die $F$ ist die Matrix der Regressoren $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}$$

Es sollte eine einfache Möglichkeit geben, sich auszudrücken

$$\frac{\partial L}{\partial \beta_j}$$

bezüglich $\frac{\partial f(x_{ij})}{\partial \beta_j}$

Answer 1

Eine verwandte Frage existiert auf math.stackexchange.com. Ableitung der Projektion in Bezug auf einen Parameter:$D_{a}: X(a)[ X(a)^TX(a) ]^{-1}X(a)^Ty$

Die Antwort schlägt vor, die Produktregel zu verwenden, die zu Folgendem führt:

$$\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}$$

Dann berechnen wir die Ableitung der Verlustfunktion als

$$L^\prime = \left( \sum (y-\hat{y})^2 \right)^\prime = \sum -2(y-\hat{y})\hat{y}^\prime$$

Wo $^\prime$ bezeichnet die Ableitung zu einem der $\beta_j$

Beispiel:

Im folgenden Beispiel passen wir die Funktion an

$$y_i = \alpha_{1} e^{\beta_1 x_{1,i}} + \alpha_2 e^{\beta_2 x_{2,i}}$$

In diesem Fall $X^\prime = \frac{\partial}{\beta_j} X$ wird das gleiche sein wie $X$ aber mit dem $i$-te Spalte multipliziert mit $x_i$ und die anderen null.

Unten finden Sie einen R-Code, der die Berechnung veranschaulicht. Es ist eine Gradientenabstiegsmethode, die die Funktion frzum Berechnen der Kostenfunktion und die Funktion grzum Berechnen des Gradienten verwendet. In dieser Funktion haben grwir die Ableitungen wie oben berechnet. Der Wert der Kostenfunktion als Funktion von$\beta_1$ und $\beta_2$ist in der folgenden Abbildung dargestellt. Die dicke schwarze Linie zeigt den Pfad, dem die Gradientenabstiegsmethode folgt.

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)

Sehen Sie sich ein historisches Schaufenster dieser Methode an:

"Die Differenzierung von Pseudo-Inversen und nichtlinearen Problemen der kleinsten Quadrate, deren Variablen sich trennen" von GH Golub und V. Pereyra im SIAM Journal on Numerical Analysis Vol. 10, No. 2 (1973), S. 413-432