분리 가능한 비선형 최소 제곱 문제에 대한 기울기를 계산하는 방법은 무엇입니까?
하나의 종속 변수가있는 비선형 최소 제곱 회귀의 경우를 고려하십시오. $y_i$ 및 두 개의 독립 변수 $x_{i1}$ 과 $x_{i2}$ 여기서 비선형 함수는 두 개의 비선형 함수의 선형 함수입니다. $f_1$ 과 $f_2$ (단순화를 위해 매개 변수 / 계수가 하나만있는 두 개의 함수와 함수로 줄 였지만 더 일반적 일 수 있습니다)
$$y_i = \sum_{j=1,2} \alpha_j f_j(x_{ij},\beta_j) + \epsilon_i$$
이 함수를 최소 제곱 회귀를 사용하여 일부 데이터에 맞추고 싶다고 가정하면, 피팅을 번갈아 가며 단계적 알고리즘으로 솔루션을 찾을 수 있습니다. $\alpha_j$ 과 $\beta_j$. 이는 유용한 접근 방식이 될 수 있습니다.$\alpha_j$ 때 $\beta_j$ 고정은 일반 최소 제곱 회귀로 쉽게 찾을 수 있습니다.
에 대한 최적화 단계를 수행하려면 $\beta_j$손실 함수의 기울기를 알아야합니다. 계산적으로 미분을 추정 할 수있는 솔버가 있지만, 우리가 미분을 직접 제공 할 수 있으면 알고리즘이 더 빠르고 정확 해집니다.
도함수를 어떻게 설명합니까? $\frac{\partial L}{\partial \beta_j}$ 잔차 손실 제곱합 함수 $$L = \Vert y - \hat{y}\Vert ^2$$
언제
$$\hat y = F (F^T F)^{-1} F^T y$$
어디 $F$ 회귀 변수의 행렬입니다. $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}$$
표현하는 간단한 방법이 있어야합니다.
$$\frac{\partial L}{\partial \beta_j}$$
측면에서 $\frac{\partial f(x_{ij})}{\partial \beta_j}$
답변
관련 질문이 math.stackexchange.com 에 있습니다. 매개 변수에 대한 투영의 미분 :$D_{a}: X(a)[ X(a)^TX(a) ]^{-1}X(a)^Ty$
대답은 다음과 같은 제품 규칙을 사용하도록 제안합니다.
$$\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}$$
그런 다음 손실 함수의 미분을 다음과 같이 계산합니다.
$$L^\prime = \left( \sum (y-\hat{y})^2 \right)^\prime = \sum -2(y-\hat{y})\hat{y}^\prime$$
어디 $^\prime$ 미분을 나타냅니다. $\beta_j$
예:
아래 예에서는 함수를
$$y_i = \alpha_{1} e^{\beta_1 x_{1,i}} + \alpha_2 e^{\beta_2 x_{2,i}}$$
이 경우 $X^\prime = \frac{\partial}{\beta_j} X$ 다음과 같을 것입니다 $X$ 하지만 $i$-번째 열 곱하기 $x_i$ 나머지는 0입니다.
다음은 계산을 설명하는 R 코드입니다. fr비용 함수를 계산하는 함수와 gr기울기를 계산하는 함수 를 사용하는 경사 하강 법입니다 . 이 함수에서 gr우리는 위와 같이 미분을 계산했습니다. 함수로서의 비용 함수의 가치$\beta_1$ 과 $\beta_2$아래 그림에 나와 있습니다. 두꺼운 검은 색 선은 경사 하강 법이 뒤 따르는 경로를 나타냅니다.
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)
이 방법에 대한 역사적인 쇼케이스를 확인하십시오.
SIAM Journal on Numerical Analysis Vol. GH Golub 및 V. Pereyra의 "변수가 분리되는 유사 역 및 비선형 최소 제곱 문제의 미분" 10, No. 2 (1973), 413-432 쪽