Rで多変数ニュートンラフソン法を行う

UPDATE: 2022-08-29 18:33:40

はじめに

このノートではRを使ったニュートンラフソン法についてまとめておく。

1変数ニュートンラフソン法

\(n\)回微分可能な関数\(f(x)\)において、異なる2点\(a,b\)を取る時、テイラー展開で近似できる。

\[ f(x) = f(a) + \frac{f^{ (1)} (a)}{1!}(x-a)^{1} + \frac{f^{ (2) }(a)}{2!}(x-a)^{2} + \frac{f^{ (3) }(a)}{3!}(x-a)^{3} + \dots \]

1次までの項で近似し、\(a = x_{i}, x=x_{i+1}\)とおくと、

\[ f(x_{i+1}) \approx f(x_{i}) + f^{(1)}(x_{i})(x_{i+1}-x_{i}) \]

となり、

\[ f(x_{i + 1}) = x_{i} - \frac{f(x_{i})}{f^{(1)}(x_{i})} \]

上記の漸化式が得られる。これを許容誤差や試行回数を設定し、繰り返すことで、\(f(x) = 0\)となる\(x\)を求めることができる。これがニュートンラフソン法である。

下記の関数の\(f(x)=0\)となる\(x\)を求める。

f <- function(x) x*log(x)-x 
curve(f(x), 0, 5)
abline(h = 0)
abline(v = 0)
title(main = 'f(x) = x*log(x) - x')

さきほどの計算式をRで実装していく。

x = 10
iter = 5
h = 0.01

for (i in seq_along(1:iter)) {
  df = (f(x + h) - f(x)) / h
  x = x - f(x) / df
  print(sprintf("iter%02d: x=%2.5f", i, x))
}

## [1] "iter01: x=4.34417"
## [1] "iter02: x=2.95865"
## [1] "iter03: x=2.72790"
## [1] "iter04: x=2.71832"
## [1] "iter05: x=2.71828"

組み込み関数のuniroot()の結果と、ほとんど同じ値が得られていることがわかる。

# (1~10)の範囲でf(x)=0の解を求める
uniroot(f, c(1, 10))

## $root
## [1] 2.718293
## 
## $f.root
## [1] 1.070079e-05
## 
## $iter
## [1] 8
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05

例えば下記の多項式は初期値によって結果が変わってしまう可能性がある。

f <- function(x) 4*x^3 - 10*x
curve(f(x), -2, 2)
abline(h = 0)
abline(v = -1.58114, col = 'red')
abline(v = 0, col = 'green')
abline(v =  1.58114, col = 'blue')
title(main = 'f(x) = 4x^3 - 10x')

初期値は10として始める。

x <- 10
iter <- 10
h <- 0.01

for (i in seq_along(1:iter)) {
  df = (f(x + h) - f(x)) / h
  x = x - f(x) / df
  print(sprintf("iter%02d: x=%2.5f", i, x))
}

## [1] "iter01: x=6.72599"
## [1] "iter02: x=4.57141"
## [1] "iter03: x=3.17736"
## [1] "iter04: x=2.31179"
## [1] "iter05: x=1.82838"
## [1] "iter06: x=1.62515"
## [1] "iter07: x=1.58324"
## [1] "iter08: x=1.58116"
## [1] "iter09: x=1.58114"
## [1] "iter10: x=1.58114"

初期値は-10として始める。

x <- -10
iter <- 10
h <- 0.01

for (i in seq_along(1:iter)) {
  df = (f(x + h) - f(x)) / h
  x = x - f(x) / df
  print(sprintf("iter%02d: x=%2.5f", i, x))
}

## [1] "iter01: x=-6.71938"
## [1] "iter02: x=-4.56055"
## [1] "iter03: x=-3.16408"
## [1] "iter04: x=-2.29793"
## [1] "iter05: x=-1.81654"
## [1] "iter06: x=-1.61873"
## [1] "iter07: x=-1.58208"
## [1] "iter08: x=-1.58113"
## [1] "iter09: x=-1.58114"
## [1] "iter10: x=-1.58114"

初期値は-0.5として始める。

x <- -0.5
iter <- 10
h <- 0.01

for (i in seq_along(1:iter)) {
  df = (f(x + h) - f(x)) / h
  x = x - f(x) / df
  print(sprintf("iter%02d: x=%2.5f", i, x))
}

## [1] "iter01: x=0.13743"
## [1] "iter02: x=-0.00237"
## [1] "iter03: x=0.00000"
## [1] "iter04: x=-0.00000"
## [1] "iter05: x=0.00000"
## [1] "iter06: x=-0.00000"
## [1] "iter07: x=0.00000"
## [1] "iter08: x=-0.00000"
## [1] "iter09: x=0.00000"
## [1] "iter10: x=-0.00000"

このように初期値で値が変わってしまう。polyroot()を使えば、先程の結果をすべて返してくれる。

# 4*x^3 - 10*x + 0の係数を左から並べる
polyroot(c(0, -10, 0, 4))

## [1]  0.000000+0i  1.581139+0i -1.581139+0i

関数の最大値、最小値を求めたい場合は、optimize()で求めることができ、$minimumが最小となる\(x\)で、$objectiveはその時の\(y\)の値を表す。

f <- function(x) x*log(x)-x 
minval <- optimize(f, c(0, 5), maximum = FALSE)
minval

## $minimum
## [1] 1.000002
## 
## $objective
## [1] -1

可視化しておく。

curve(f(x), 0, 5)
abline(h = minval$objective)
abline(v = minval$minimum)
title(main = 'f(x) = x*log(x) - x')

2変数ニュートンラフソン法

\(n\)回微分可能な関数\(f(x, y)\)において、テイラー展開は下記の通りとなる。

\[ f(a+h, b+k) = f(a,b) + \frac{1}{1!}f_{x}(a,b)(x-a) + \frac{1}{1!}f_{y}(a,b)(y-b) + \frac{1}{2!}\left \{ f_{xx}(a,b)(x-a)^{2} + 2f_{xy}(a,b)(x-a)(y-b) + f_{yy}(a,b)(y-b)^{2} \right \} + \dots \]

\(h\)は\(x\)軸、\(k\)は\(y\)軸の微小変化を表し、\(x=a+h \Leftrightarrow h=x-a, y=b+k \Leftrightarrow k=y-b\)として1次の項で近似する。また、\(n\)回微分可能な関数\(g(x, y)\)があると、

\[ \begin{eqnarray} f(a+h,b+k) &\approx& f(a,b) + f_{x}(a,b)h + f_{y}(a,b)k \\ g(a+h,b+k) &\approx& g(a,b) + g_{x}(a,b)h + g_{y}(a,b)k \\ \end{eqnarray} \] となる。ここで、

\[ \begin{eqnarray} \Delta f &=& f(a+h,b+k) - f(a,b) \\ \Delta g &=& g(a+h,b+k) - g(a,b) \end{eqnarray} \]

とおくと、

\[ \begin{eqnarray} \Delta f &\approx& f_{x}(a,b)h + f_{y}(a,b)k \\ \Delta g &\approx& g_{x}(a,b)h + g_{y}(a,b)k \end{eqnarray} \]

と表せる。ここで、表記を行列に直し、

\[ \begin{pmatrix} \Delta f \\ \Delta g \end{pmatrix} = \begin{pmatrix} f_{x} & f_{y} \\ g_{x} & g_{y} \end{pmatrix} \begin{pmatrix} h \\ k \end{pmatrix} \]

また、1変数のときと同様に、

\[ \begin{eqnarray} \Delta f &=& 0 - f(x_{i},y_{i}) \\ \Delta f &=& 0 - g(x_{i},y_{i}) \\ h &=& x_{i+1} - x_{i} \\ k &=& y_{i+1} - y_{i} \\ H &=& \begin{pmatrix} f_{x} & f_{y} \\ g_{x} & g_{y} \end{pmatrix} \end{eqnarray} \]

を使って、書き直す。\(H\)はヤコビアンである。

\[ \begin{pmatrix} 0 - f(x_{i},y_{i}) \\ 0 - g(x_{i},y_{i}) \end{pmatrix} = H \begin{pmatrix} x_{i+1} - x_{i} \\ y_{i+1} - y_{i} \end{pmatrix} \]

あとは、これをヤコビアンの逆行列\(H^{-1}\)を使って変形すれば、2変数のニュートンラフソン方の更新式が得られる。

\[ \begin{pmatrix} x_{i+1} \\ y_{i+1} \end{pmatrix} = \begin{pmatrix} x_{i} \\ y_{i} \end{pmatrix} -H^{-1} \begin{pmatrix} f(x_{i},y_{i}) \\ g(x_{i},y_{i}) \end{pmatrix} \]

Rで実装していく。解を求める関数を可視化しておく。

fxy <- function(x,y) x^2 + y^2 - 1
x <- y <- seq(-10, 10, 0.1)
z <- outer(x,y,fxy)

library(plot3D)
persp3D(
  x = x,
  y = y,
  z = z,
  color.palette = heat.colors,
  phi = 45,
  theta = 200,
  main = "3D perspective")

gxy <- function(x,y) x^3 - y
z <- outer(x, y, gxy)
persp3D(
  x = x,
  y = y,
  z = z,
  color.palette = heat.colors,
  phi = 45,
  theta = 200,
  main = "3D perspective")

2変数ニュートンラフソン法で\(f(x,y)=0, g(x,y)=0\)を満たす\(x,y\)を求める。

# fxy <- function(x,y) x^2 + y^2 - 1
# gxy <- function(x,y) x^3 - y

iter <- 5
h <- 0.01
init <- 1
I <- t(t(rep(init,2)))

for (i in seq_along(1:iter)) {
  # 初期化
  df <- dg <- c()
  
  for (j in 1:length(I)) {
    vec = I
    vec[j] = vec[j] + h
    df <- c(df, (fxy(vec[1], vec[2]) - fxy(I[1], I[2])) / h )
    dg <- c(dg, (gxy(vec[1], vec[2]) - gxy(I[1], I[2])) / h )
  }
  H <- t(matrix(c(df, dg), ncol = length(I), nrow = length(I)))
  cat('---- Jacobian ----\n')
  print(H)
  cat('------------------\n')
  I <- I - solve(H) %*% c(fxy(I[1], I[2]), gxy(I[1], I[2]))
  print(sprintf("%d times: (x=%2.3f, y=%2.3f)", i, I[1], I[2]))
}

## ---- Jacobian ----
##        [,1]  [,2]
## [1,] 2.0100  2.01
## [2,] 3.0301 -1.00
## ------------------
## [1] "1 times: (x=0.877, y=0.626)"
## ---- Jacobian ----
##          [,1]      [,2]
## [1,] 1.763102  1.261873
## [2,] 2.331421 -1.000000
## ------------------
## [1] "2 times: (x=0.830, y=0.564)"
## ---- Jacobian ----
##          [,1]      [,2]
## [1,] 1.669524  1.138811
## [2,] 2.090507 -1.000000
## ------------------
## [1] "3 times: (x=0.826, y=0.564)"
## ---- Jacobian ----
##          [,1]      [,2]
## [1,] 1.662164  1.137205
## [2,] 2.072117 -1.000000
## ------------------
## [1] "4 times: (x=0.826, y=0.564)"
## ---- Jacobian ----
##          [,1]      [,2]
## [1,] 1.662064  1.137248
## [2,] 2.071867 -1.000000
## ------------------
## [1] "5 times: (x=0.826, y=0.564)"

1回目の更新の計算イメージは下記の通り。このような感じで繰り返されて更新される。

更新イメージ

nleqslvパッケージののnleqslv()の結果とほとんど同じ値が得られていることがわかる。

# fxy <- function(x,y) x^2 + y^2 - 1
# gxy <- function(x,y) x^3 - y
library(nleqslv)
f <- function(x){ 
  r = c(
    x[1]^2 + x[2]^2 - 1, 
    x[1]^3 - x[2]
    )
  return(r)
}
nleqslv(c(1, 1), f)$x

## [1] 0.8260314 0.5636242

正規分布のパラメタ推定

2変数ニュートンラフソン法を使えば、正規分布から生成されたデータから、パラメタを推定できる。解析的に正規分布の最尤推定量は計算できるので、わざわざ数値的に計算する必要はないが、練習のために行っておく。正規分布は、

\[ f(x; \mu, \sigma^2)=\frac{1}{\sqrt{2 \pi \sigma^2}}exp \left[-\frac{(x-\mu)^2}{2 \sigma^2}\right] \]

であり、対数尤度関数は、

\[ \begin{eqnarray} L &=& - \frac{n}{2}log(2 \pi \sigma^2) - \frac{1}{2 \sigma^2} \sum (x_{i}-\mu)^2 \\ \frac{ \partial L(\mu, \sigma^2) }{ \partial \mu } &=& - \frac{1}{\sigma^2} \sum (x_{i}-\mu) = 0 \\ \frac{ \partial L(\mu, \sigma^2) }{ \partial \sigma^2 } &=& - \frac{n}{2\sigma^2} + \frac{1}{2\sigma^4}\sum (x_{i}-\mu)^2 = 0 \\ \end{eqnarray} \]

となる。これを変形すれば、パラメタの最尤推定量の式が出てくるので、この式に対してニュートンラフソン法でパラメタを計算する。

\[ \begin{eqnarray} f(\mu, \sigma^2) &=& \sum x_{i} - n\mu = 0 \\ g(\mu, \sigma^2) &=& \sum (x_{i} - \mu)^2 - n \sigma^2 = 0 \end{eqnarray} \] Rで実装していく。ここでは学習率etaを利用してみた。

# mu = 5, sigma2 = 9
set.seed(1989)
x <- rnorm(100, mean = 5, sd = sqrt(9))

fmu <- function(z) sum(x) - length(x)*z[1]
gsig <- function(z) sum((x-z[1])^2) - length(x)*z[2]

iter <- 1000
h <- 0.01
I <- t(t(c(1,1)))
eta <- 0.01

for (i in seq_along(1:iter)) {
  # 初期化
  df <- dg <- c()
  
  for (j in 1:length(I)) {
    vec = I
    vec[j] = vec[j] + h
    df <- c(df, (fmu(vec) - fmu(I)) / h )
    dg <- c(dg, (gsig(vec) - gsig(I)) / h )
  }
  H <- t(matrix(c(df, dg), ncol = length(I), nrow = length(I)))

  I <- I - eta*solve(H) %*% c(fmu(I), gsig(I))
  if(i %% 50 == 0) {
    print(sprintf("%d times: (mu=%2.3f, sigma2=%2.3f)", i, I[1], I[2]))
  }
}

## [1] "50 times: (mu=2.476, sigma2=1.042)"
## [1] "100 times: (mu=3.368, sigma2=3.199)"
## [1] "150 times: (mu=3.908, sigma2=5.283)"
## [1] "200 times: (mu=4.235, sigma2=6.828)"
## [1] "250 times: (mu=4.433, sigma2=7.867)"
## [1] "300 times: (mu=4.553, sigma2=8.533)"
## [1] "350 times: (mu=4.625, sigma2=8.949)"
## [1] "400 times: (mu=4.669, sigma2=9.206)"
## [1] "450 times: (mu=4.695, sigma2=9.364)"
## [1] "500 times: (mu=4.711, sigma2=9.459)"
## [1] "550 times: (mu=4.721, sigma2=9.518)"
## [1] "600 times: (mu=4.727, sigma2=9.553)"
## [1] "650 times: (mu=4.730, sigma2=9.574)"
## [1] "700 times: (mu=4.732, sigma2=9.587)"
## [1] "750 times: (mu=4.734, sigma2=9.595)"
## [1] "800 times: (mu=4.735, sigma2=9.599)"
## [1] "850 times: (mu=4.735, sigma2=9.602)"
## [1] "900 times: (mu=4.735, sigma2=9.604)"
## [1] "950 times: (mu=4.735, sigma2=9.605)"
## [1] "1000 times: (mu=4.736, sigma2=9.606)"

実際の値とほとんど一致していることがわかる。(この精度ではだめなのかもしれないが。)

list(
  mean = mean(x),
  sigma2 = var(x)
  )

## $mean
## [1] 4.735709
## 
## $sigma2
## [1] 9.703747