Rでガウス・ニュートン法とレーベンバーグ・マーカート法を実装

これなら分かる最適化数学 4.3.1 に記載の通り、ガウス・ニュートン法を用いて、非線形最小二乗法の誤差関数

J(\mathbf{u})=\frac{1}{2}\sum_{\alpha=1}^N \sum_{l=1}^r F_l(\mathbf{x}_\alpha,\mathbf{u})^2

に最小値を取らせる変数$\mathbf{u}$の値を求めます。

この例では、観測した点$\mathbf{x}\alpha=(x{\alpha1},x_{\alpha2})^T , \alpha=1,\dots,N$に対して関数$x_2=u_1*(x_1-u_2)^2$を当てはめるとき、誤差関数を最小にするようなパラメータ$\mathbf{u}=(u_1,u_2)^T$を求める問題を想定し、$r=1$, $F_1(\mathbf{x},\mathbf{u})=u_1*(x_1-u_2)^2-x_2$としています。なお、$\mathbf{x}\alpha$を生成する際には、正規分布に従うノイズを加え、$x{\alpha2}=u_1*(x_{\alpha1}-u_2)^2+N(0,\sigma^2),\mathbf{u}=(2,0.5)^T$としています。

また、4.3.2に記載の通り、レーベンバーグ・マーカート法を用いて、同じ問題を解きます。更新式は

\mathbf{u}^{(K+1)}=\mathbf{u}^{(K)}-(\mathbf{H_u}^{(K)}+cD[\mathbf{H_u}])^{-1}\nabla_\mathbf{u}J^{(K)}

としていますが、if文による切り替えで、

\mathbf{u}^{(K+1)}=\mathbf{u}^{(K)}-(\mathbf{H_u}^{(K)}+c\mathbf{I})^{-1}\nabla_\mathbf{u}J^{(K)}

とすることもできます。

library(pracma)

# Taken and modified from: http://cran.r-project.org/web/packages/pracma/index.html
vectorfield <- function(fun, xlim, ylim, n = 16,
                        scale = 0.05, col = "green", ...) {
    stopifnot(is.numeric(xlim), length(xlim) == 2,
              is.numeric(ylim), length(ylim) == 2)

	xpts <- linspace(xlim[1], xlim[2], n)
	ypts <- linspace(ylim[1], ylim[2], n)

	M <- meshgrid(xpts, ypts)
	x <- M$X
	y <- M$Y

	f <- fun(x, y)
	stopifnot(nrow(f) <= 2)
	if (nrow(f) == 1) {
		px = matrix(1, nrow=n , ncol=n)
		py = fun(x, y);
	} else {
		px = matrix(f[1, ], nrow=n , ncol=n)
		py = matrix(f[2, ], nrow=n , ncol=n)
	}
	
	plot(xlim, ylim, type="n", xlab = "", ylab = ""); grid()
	quiver(x, y, px, py, scale = scale, col = col, ...)
}

# Taken and modified from: http://cran.r-project.org/web/packages/pracma/index.html
quiver <- function(x, y, u, v,
                    scale = 0.05, angle = 10, length = 0.1, ...) {
    stopifnot(is.numeric(x), is.numeric(y), is.numeric(u), is.numeric(v))
    
	arrows(x, y, x+scale*u, y+scale*v, angle=angle, length=length, ...)
}

frame()
set.seed(0)
par(mfrow=c(1, 2))
par(mar=c(2, 2, 1.5, 0.1))
par(mgp=c(1, 0.2, 0))
u1 <- seq(-1, 4, 0.01)
u2 <- seq(0, 2, 0.01)
x1 <- seq(0, 1, 0.05)
x2 <- 2 * (x1 - 0.5) ^ 2 + rnorm(length(x1), 0, 0.2)
f <- function(x1, x2, u) { u[1] * (x1 - u[2]) ^ 2 - x2 }
du1 <- function(x1, u) { (x1 - u[2]) ^ 2 }
du2 <- function(x1, u) { -2 * u[1] * (x1 - u[2]) }
u0 <- c(-0.8, 0.8)

gaussnewton <- function(u0, f, du1, du2) {
	
	u <- as.matrix(u0)
	eps <- 1e-12
	i <- 1
	repeat {
		oldu <- u
		# 勾配∇Jを求める。
		gradientj <- rowSums(f(x1, x2, u) * rbind(du1(x1, u), du2(x1, u)))
		# ヘッセ行列のガウスニュートン近似Hを求める。
		H <- matrix(rowSums((sapply(x1, function(x1) {
			c(du1(x1, u), du2(x1, u)) %*% t(c(du1(x1, u), du2(x1, u)))
		}))), nrow=length(u))
		# HΔu=-∇Jを解き、Δu=-H^-1 ∇Jを求める。
		deltau <- solve(H, -gradientj)
		u <- u + deltau
		xy <- cbind(oldu, u)
		lines(xy[1, ], xy[2, ], type="o", col=2)
		cat("i=", i, "gradientj=", gradientj, "u=" , u, "\n")
		# ||Δu||^2 が十分小さくなるまで繰り返す。
		if (sum(deltau * deltau) < eps) {
			break
		}
		i <- i + 1
	}
	u
}

levenbergmarquardt <- function(u0, f, du1, du2) {
	
	c <- 0.0001
	u <- as.matrix(u0)
	eps <- 1e-12
	i <- 0
	repeat {
		oldu <- u
		j <- 1 / 2 * sum(f(x1, x2, u)^2)
		# 勾配∇Jを求める。
		gradientj <- rowSums(f(x1, x2, u) * rbind(du1(x1, u), du2(x1, u)))
		# ヘッセ行列のガウスニュートン近似Hを求める。
		H <- matrix(rowSums((sapply(x1, function(x1) {
			c(du1(x1, u), du2(x1, u)) %*% t(c(du1(x1, u), du2(x1, u)))
		}))), nrow=length(u))
		repeat {
			if (T) {
				# (H+cD[H])Δu=-∇Jを解き、Δu=-(H+cD[H])^-1 ∇Jを求める。
				deltau <- solve(H + c * diag(diag(H)), -gradientj)
			} else {
				# (H+cI)Δu=-∇Jを解き、Δu=-(H+cI)^-1 ∇Jを求める。
				deltau <- solve(H + c * diag(length(u)), -gradientj)
			}
			newu <- u + deltau
			newj <- 1 / 2 * sum(f(x1, x2, newu)^2)
			# 解から離れてしまったときはcを大きくしてガウスニュートン法の影響を弱め、
			# Δuを求め直す。
			if (newj <= j) {
				break
			}
			cat("newj=", newj, "j=", j, "c=", c, "\n")
			c <- c * 10
		}
		# 解に近づくにつれ、cを小さくしてガウスニュートン法の影響を強める。
		u <- u + deltau
		c <- c / 10
		i <- i + 1
		xy <- cbind(oldu, u)
		lines(xy[1, ], xy[2, ], type="o", col=6)
		cat("i=", i, "gradientj=", gradientj, "u=" , u, "c=", c, "\n")
		# ||Δu||^2 が十分小さくなるまで繰り返す。
		if (sum(deltau * deltau) < eps) {
			break
		}
	}
	u
}

j <- outer(u1, u2, Vectorize(function(u1, u2) { 1 / 2 * sum(f(x1, x2, c(u1, u2))^2) }))
contour(u1, u2, j, xlab="u1", ylab="u2", main="J=Σ_x[(u1*(x1-u2)^2-x2)^2], ∇Ju")
image(u1, u2, j, col=rainbow(450)[256:1], add=T)
contour(u1, u2, j, xlab="", ylab="", levels=range(j)[1] + abs(diff(range(j))) * seq(0, 1, 0.1) ^ 4, add=T)
v <- Vectorize(function(u1, u2) {
	as.vector(rowSums(f(x1, x2, c(u1, u2)) * rbind(du1(x1, c(u1, u2)), du2(x1, c(u1, u2)))))
	})
par(new=T)
vectorfield(v, range(u1), range(u2), n=26, angle=20, length=0.05)
# l <- locator(1); u0 <- c(l$x, l$y)
ugn <- gaussnewton(u0, f, du1, du2)
ulm <- levenbergmarquardt(u0, f, du1, du2)
legend("bottomright", legend=c(
	"GaussNewton", 
	"LevenbergMarquardt"), 
	col=c(2, 6), pch=1, bg="white")

plot(x1, x2, main="x2=u1*(x1-u2)^2+N,u=(2 0.5)")
s <- function(x){ 2 * (x - 0.5) ^ 2 }
curve(s, col=1, add=T)
fitted <- function(x){ ugn[1] * (x - ugn[2]) ^ 2 }
curve(fitted, col=2, add=T)
fitted <- function(x){ ulm[1] * (x - ulm[2]) ^ 2 }
curve(fitted, col=6, add=T)