分離可能な非線形最小二乗問題の勾配を計算する方法は？
1つの従属変数を持つ非線形最小二乗回帰の場合を考えてみましょう。 $y_i$ および2つの独立変数 $x_{i1}$ そして $x_{i2}$ ここで、非線形関数は2つの非線形関数の線形関数です。 $f_1$ そして $f_2$ （簡単にするために、これを2つの関数と1つのパラメーター/係数のみを持つ関数に減らしますが、より一般的にすることもできます）
$$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$ そして他はゼロ。
以下は、計算を説明するいくつかの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)
この方法の歴史的なショーケースについては、を参照してください。
「変数が分離する疑似逆行列問題と非線形最小二乗問題の微分」、GHGolubとV.PereyraによるSIAMJournal on Numerical AnalysisVol。10、No。2（1973）、pp.413-432