¿Cómo calcular el gradiente para un problema de mínimos cuadrados no lineales separables?

Question

Considere el caso de la regresión de mínimos cuadrados no lineal con una variable dependiente $y_i$ y dos variables independientes $x_{i1}$ y $x_{i2}$ donde la función no lineal es una función lineal de dos funciones no lineales $f_1$ y $f_2$ (para simplificar, reduzco esto a dos funciones y funciones con un solo parámetro / coeficiente, pero puede ser más general)

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

Digamos que deseamos usar ajustar esta función a algunos datos con regresión de mínimos cuadrados, entonces podríamos encontrar la solución con un algoritmo paso a paso que alterna entre ajustar el $\alpha_j$ y $\beta_j$. Este puede ser un enfoque útil porque la solución para el$\alpha_j$ cuando el $\beta_j$ son fijos se encuentra fácilmente mediante regresión de mínimos cuadrados ordinarios.

Para realizar el paso de optimización para el $\beta_j$necesitamos conocer el gradiente de la función de pérdida. Hay solucionadores que pueden estimar las derivadas computacionalmente, pero los algoritmos serán más rápidos y precisos cuando podamos proporcionar las derivadas nosotros mismos.

¿Cómo describimos la derivada? $\frac{\partial L}{\partial \beta_j}$ de la suma de la función de pérdida de residuos al cuadrado $$L = \Vert y - \hat{y}\Vert ^2$$

cuando

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

donde el $F$ es la matriz de los regresores $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}$$

Debería haber alguna forma sencilla de expresar

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

en términos de $\frac{\partial f(x_{ij})}{\partial \beta_j}$

Answer 1

Existe una pregunta relacionada en math.stackexchange.com Derivada de proyección con respecto a un parámetro:$D_{a}: X(a)[ X(a)^TX(a) ]^{-1}X(a)^Ty$

La respuesta sugiere usar la regla del producto que conduce a:

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

Luego calculamos la derivada de la función de pérdida como

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

Dónde $^\prime$ denota la derivada de cualquiera de los $\beta_j$

Ejemplo:

En el siguiente ejemplo, ajustamos la función

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

En este caso $X^\prime = \frac{\partial}{\beta_j} X$ será lo mismo que $X$ pero con el $i$-th columna multiplicada por $x_i$ y los demás cero.

A continuación se muestra un código R que ilustra el cálculo. Es un método de descenso de gradiente que usa la función frpara calcular la función de costo y la función grpara calcular el gradiente. En esta función grhemos calculado las derivadas como arriba. El valor de la función de costo en función de$\beta_1$ y $\beta_2$se muestra en la figura siguiente. La línea negra gruesa muestra el camino que sigue el método de descenso de gradiente.

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)

Vea una muestra histórica de este método:

"La diferenciación de pseudo-inversos y problemas de mínimos cuadrados no lineales cuyas variables se separan" por GH Golub y V. Pereyra en SIAM Journal on Numerical Analysis Vol. 10, núm. 2 (1973), págs. 413-432