Como calcular o gradiente para um problema de quadrados mínimos não linear separados?
Considere o caso de regressão de mínimos quadrados não linear com uma variável dependente $y_i$ e duas variáveis independentes $x_{i1}$ e $x_{i2}$ onde a função não linear é uma função linear de duas funções não lineares $f_1$ e $f_2$ (para simplificar, reduzo isso a duas funções e funções com apenas um parâmetro / coeficiente, mas pode ser mais geral)
$$y_i = \sum_{j=1,2} \alpha_j f_j(x_{ij},\beta_j) + \epsilon_i$$
Digamos que desejamos usar o ajuste desta função para alguns dados com regressão de mínimos quadrados, então poderíamos encontrar a solução com um algoritmo passo a passo que alterna entre o ajuste do $\alpha_j$ e $\beta_j$. Esta pode ser uma abordagem útil porque a solução para o$\alpha_j$ quando o $\beta_j$ são fixos é facilmente encontrado por regressão de mínimos quadrados ordinários.
Para realizar a etapa de otimização para o $\beta_j$precisamos saber o gradiente da função de perda. Existem solucionadores que podem estimar as derivadas computacionalmente, mas os algoritmos serão mais rápidos e precisos quando pudermos fornecer as derivadas nós mesmos.
Como descrevemos a derivada $\frac{\partial L}{\partial \beta_j}$ da soma dos quadrados da função de perda de resíduos $$L = \Vert y - \hat{y}\Vert ^2$$
quando
$$\hat y = F (F^T F)^{-1} F^T y$$
onde o $F$ é a matriz dos regressores $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}$$
Deve haver uma maneira simples de expressar
$$\frac{\partial L}{\partial \beta_j}$$
em termos de $\frac{\partial f(x_{ij})}{\partial \beta_j}$
Respostas
Uma pergunta relacionada existe em math.stackexchange.com Derivada da projeção em relação a um parâmetro:$D_{a}: X(a)[ X(a)^TX(a) ]^{-1}X(a)^Ty$
A resposta sugere o uso da regra do produto que leva 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}$$
Em seguida, calculamos a derivada da função de perda como
$$L^\prime = \left( \sum (y-\hat{y})^2 \right)^\prime = \sum -2(y-\hat{y})\hat{y}^\prime$$
Onde $^\prime$ denota a derivada de qualquer um dos $\beta_j$
Exemplo:
No exemplo abaixo, ajustamos a função
$$y_i = \alpha_{1} e^{\beta_1 x_{1,i}} + \alpha_2 e^{\beta_2 x_{2,i}}$$
Nesse caso $X^\prime = \frac{\partial}{\beta_j} X$ será o mesmo que $X$ mas com o $i$-ésima coluna multiplicada por $x_i$ e os outros zero.
Abaixo está algum código R que ilustra o cálculo. É um método de gradiente descendente que usa a função frpara calcular a função de custo e a função grpara calcular o gradiente. Nesta função gr, calculamos as derivadas como acima. O valor da função de custo em função de$\beta_1$ e $\beta_2$é mostrado na figura abaixo. A linha preta grossa mostra o caminho seguido pelo método de descida 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)
Consulte para obter uma demonstração histórica deste método:
"The Differentiation of Pseudo-Inverses and Nonlinear Least Squares Problems Whose Variables Separate" por GH Golub e V. Pereyra em SIAM Journal on Numerical Analysis Vol. 10, No. 2 (1973), pp. 413-432