汎用非線形最小化関数 nlm の編集

//間瀬茂
COLOR(red){SIZE(30){汎用最適化関数 nlm() の使用法}}

統計学では、最尤推定や最小自乗法等、ある目的関数の最適化を行なう必要が
しばしば起こる。関数 nlm (Non-Linear Minimization)
はニュートン法に基づく非線形関数の最小化関数で、r-help 記事によればしばしば optim() 関数の諸手法よりも良い結果を与えるらしい。

#contents
~

*SIZE(20){nlm()} の使用法

**命令書式

nlm(f, p, hessian = FALSE, typsize=rep(1, length(p)), fscale=1,
         print.level = 0, ndigit=12, gradtol = 1e-6,
         stepmax = max(1000 * sqrt(sum((p/typsize)^2)), 1000),
         steptol = 1e-6, iterlim = 100, check.analyticals = TRUE, ...)

**引数:

- COLOR(red){f}　: 最小化すべき関数。もしこれが gradient 属性（f の偏微分ベクトル関数の定義本体、名前ではない）を持つか、gradient, hessian 属性（hessian 属性は f の二階偏微分行列関数の定義本体、名前ではない）の双方を持てば最適化の途中で用いられる。さもなければ gradient, hessian の値が数値微分で使われる。
- COLOR(red){p}　: パラメータの初期値。
- COLOR(red){hessian}　: もし 'TRUE' ならば最終段階での f の hessian が返される。
- COLOR(red){typsize}　: 最小値における各パラメータのサイズの推定値。
- COLOR(red){fscale}　: 最小値における 'f' のサイズの推定値。
- COLOR(red){print.level}　: 最小化の各段階での出力の水準を決める。既定値 0 は未出力、値 1 は初期値と最終値を出力。値 2 はすべての途中結果を出力。 
- COLOR(red){ndigit}　: 関数 f の有効桁数。
- COLOR(red){gradtol}　: アルゴリズムの収束判定で用いられる許容値で、スケール化された gradiennt が零に十分近いかどうかを決める。スケール化された gradient は各 p[i] の変化に対する f の値の増分を、 p[i] の変化量で割った相対変化の尺度。
- COLOR(red){stepmax}　: 許される最大のスケール化されたステップ長を与える正スカラー。stepmax は、目的関数がオーバーフローするようなステップを避けたり、アルゴリズムが関心のあるパラメータ空間を無視するのを避けたり、またアルゴリズムの発散を避けるのに使われる。stepmax はこれらのうち最初の二つが起こるのを防ぐに十分なほど小さく取るべきであるが、事前に想定される合理的なステップよりも大きくなるようにとるべきである。
- COLOR(red){steptol}　: 正のスカラーで許される最小の相対的ステップ長を与える。
- COLOR(red){iterlim}　: 正整数で、プログラム終了までに実行される最大繰返し回数を与える。
- COLOR(red){check.analyticals}　: 論理値で、もし解析的な gradient と hessian が与えられている時、初期値においてそれらが数値的な微分値と整合性があるかどうかを検査する。不正確な gradient と hessian 式が与えられた場合をチェックするのを助ける。
- COLOR(red){...}　: f に対する追加引数。

** 詳細

もし gradient や hessian が与えられても好ましくないモードや長さを与えるようなら、check.analyticals = TRUE (既定動作) が指定されていれば警告とともに無視される。
hessian は、 gradient が与えられていなかったり、チェックをパスしない限り、無視される。

原典で与えられている三種類の手法のうち、直線探索である手法 １ だけが使われる。

**返り値

次の成分を含むリスト（str 関数を使うと簡潔に一覧できる）

- COLOR(red){ minimum}　: f の最小値の推定値。
- COLOR(red){estimate}　: f の最小値が得られたパラメータ値。
- COLOR(red){gradient}　: 最小値における gradient の値。
- COLOR(red){hessian}　: f の最初値における hessian の値（要求された時だけ）。
- COLOR(red){code}　: 最適化過程の終了状態を示す整数値。
-- 1:　 gradient の相対値は零に近い。おそらく最小値が得られた。
-- 2:　 引き続く繰返しが許容値以内。おそらく最小値が得られた。
-- 3:　 最後の大局的ステップが推定値よりも低い点を見つけるのに失敗。推定値が近似的な局所的最小値であるか、steptol が小さ過ぎる。
-- 4:　 繰返し回数が上限を越えた。
-- 5:　 ５回連続して最大ステップサイズ stepmax を超過。関数が下に非有界か、ある方向で有限値に上から漸近的に近付いているか、stepmax が小さ過ぎる。
- COLOR(red){iterations}　:　実行される繰返し回数。

* 例１ demo(nlm) より 
Helical Valley Function. Page 362 Dennis + Schnabel

nlm1 <- function () {
   theta <- function(x1,x2) atan(x2,x1)/(2*pi)
   f <- function(x) {                          # 目的関数定義
     f1 <- 10*(x[3] - 10*theta(x[1],x[2]))
     f2 <- 10*(sqrt(x[1]^2+x[2]^2)-1)
     f3 <- x[3]
     return(f1^2+f2^2+f3^2)
   }
   # x3=0 における曲面を見る
   x <- seq(-1, 2, len=50)
   y <- seq(-1, 1, len=50)
   # z は x,y のすべての組み合わせに対する関数値を順に並べたベクトル
   z <- apply(as.matrix(expand.grid(x, y)), 1, function(x) f(c(x, 0))) 
   contour(x, y, matrix(log10(z), 50, 50))     # (対数値)曲面の等高線図
   nlm.f <- nlm(f, c(-1,0,0), hessian=TRUE, print=0)
   points(rbind(nlm.f$estim[1:2]), col = "red", pch = 20) # 等高線図に最適パラメータをプロット
   nlm.f   
 }

> nlm.f <- nlm1()
 > nlm.f            # 非線形最小化の結果 
 $minimum           # 最小値
 [1] 1.238237e-14 
 $estimate          # その時のパラメータ値   
 [1]  1.000000e-00  3.071089e-09 -6.063445e-09
 $gradient          # 収束点での偏微分ベクトル gradient (ゼロに近いはず)
 [1] -3.755775e-07  3.485886e-06 -2.202374e-06
 $hessian           # 収束点での二階導関数行列 hessian (正定値符合行列のはず)
               [,1]          [,2]          [,3]
 [1,]  2.000000e+02   -0.04065842  9.774602e-07
 [2,] -4.065842e-02  506.60589960 -3.183099e+02
 [3,]  9.774602e-07 -318.30988572  2.020000e+02
 $code              # 収束判定コード (おそらく成功)
 [1] 2
 $iterations        # 終了までの反復回数 
 [1] 27
 
 > str(nlm.f)       # 結果の簡潔な要約
 List of 6
  $ minimum   : num 1.24e-14
  $ estimate  : num [1:3]  1.00e-00  3.07e-09 -6.06e-09
  $ gradient  : num [1:3] -3.76e-07  3.49e-06 -2.20e-06
  $ hessian   : num [1:3, 1:3]  2.00e+02 -4.07e-02  9.77e-07 -4.07e-02 5.07e+02 ...
  $ code      : int 2
  $ iterations: int 27

> eigen(nlm.f$hessian) 
 $values                 # hessian の固有値はすべて正（正定値符合行列！）
 [1] 707.173143 200.000000   1.432756
 $vectors
               [,1]          [,2]         [,3]
 [1,]  6.782635e-05  1.000000e-00 0.0001091525
 [2,] -8.460531e-01 -8.043624e-07 0.5330985699
 [3,]  5.330986e-01 -1.285070e-04 0.8460531324

#ref(nlm1.png, left)

* 例２ demo(nlm) より
### the Rosenbrock banana valley function

nlm2 <- function () {
   fR <- function(x) {  # 目的関数定義
     x1 <- x[1]; x2 <- x[2]
     100*(x2 - x1*x1)^2 + (1-x1)^2
   }
  fx <- function(x) {   # fR のベクトル化版 (z をスマートに定義するため)
    x1 <- x[,1]; x2 <- x[,2]
    100*(x2 - x1*x1)^2 + (1-x1)^2
  }
  x <- seq(-2, 2, len=100)
  y <- seq(-.5, 1.5, len=100)
  z <- fx(expand.grid(x, y))
  op <- par(mfrow = c(2,1), mar = .1 + c(3,3,0,0)) # グラフィックス画面を上下に二分割
  contour(x, y, matrix(log10(z), length(x)))       # 曲面の等高線を眺める
  nlm.f2 <- nlm(fR, c(-1.2, 1), hessian=TRUE)      # 最小化実行
  points(rbind(nlm.f2$estim[1:2]), col = "red", pch = 20)  # 等高線図に最適パラメータをプロット
  ## 拡大曲面を画く
  rect(.9, .9, 1.1, 1.1, border = "orange", lwd = 2)
  x <- y <- seq(.9, 1.1, len=100)
  z <- fx(expand.grid(x, y))
  contour(x, y, matrix(log10(z), length(x)))   # 曲面の等高線を眺める
  mtext("zoomed in")                           # 注釈を入れる
  box(col = "orange")                          # オレンジ色の枠を描
  points(rbind(nlm.f2$estim[1:2]), col = "red", pch = 20)   # 等高線図に最適パラメータをプロット
  par(op)       # グラフィックスパラメータ復帰
  nlm.f2
 }

> nlm.f2 <- nlm2()
 > nlm.f2
 $minimum           # 真値は当然ゼロ
 [1] 3.973766e-12
 $estimate          # 真のパラメータは当然 c(1,1) 
 [1] 0.999998 0.999996
 $gradient
 [1] -6.539279e-07  3.335998e-07
 $hessian
           [,1]      [,2]
 [1,]  802.2368 -400.0192
 [2,] -400.0192  200.0000
 $code  # 満足のいく終了
 [1] 1
 $iterations
 [1] 23

> str(nlm.f2)
 List of 6
  $ minimum   : num 3.97e-12
  $ estimate  : num [1:2] 1 1
  $ gradient  : num [1:2] -6.54e-07  3.34e-07
  $ hessian   : num [1:2, 1:2]  802 -400 -400  200
  $ code      : int 1
  $ iterations: int 23
 > det(nlm.f2$hessian) # hessian の行列式値（正、つまり正定値符合行列！）
 [1] 432.0022 
#ref(nlm2.png, left)

* 例３ demo(nlm) より (gradient を与える例)

nlm3 <- function () {
   fg <- function(x) {  # x は二次元ベクトル変数 (x[1],x[2])
     gr <- function(x1, x2) { # fg の偏微分(gradient)ベクトル関数
         c(-400*x1*(x2 - x1*x1)-2*(1-x1), 200*(x2 - x1*x1))
     }
     x1 <- x[1]; x2 <- x[2]
     res<- 100*(x2 - x1*x1)^2 + (1-x1)^2 # fg の本体
     attr(res, "gradient") <- gr(x1, x2) # fg の gradient 属性(関数本体)は gr であることを教える
     return(res)
   }
   nlm(fg, c(-1.2,1), hessian=TRUE) # 非線形回帰実行
 }

> nlm3.f <- nlm3() # 非線形回帰実行
 > nlm3.f           # 結果表示
 $minimum           # 最小値
 [1] 1.182096e-20
 $estimate          # その時のパラメータ値 
 [1] 1 1 
 $gradient          # その時の gradient 値 (数値誤差範囲内で零に近い！)
 [1]  2.583521e-09 -1.201128e-09
 $hessian           # その時の hessian 値
         [,1]    [,2]
 [1,]  802.24 -400.02
 [2,] -400.02  200.00
 $code              # ベストな終了判定コード
 [1] 1
 $iterations        # 収束までに24回の繰返しが行なわれた
 [1] 24
 
 > str(nlm3.f)  # 簡潔な要約
 List of 6
  $ minimum   : num 1.18e-20
  $ estimate  : num [1:2] 1 1
  $ gradient  : num [1:2]  2.58e-09 -1.20e-09
  $ hessian   : num [1:2, 1:2]  802 -400 -400  200
  $ code      : int 1
  $ iterations: int 24

# 目的関数の内容
 # gradient の関数本体が属性として fg の定義の一部を構成することを注意。
 > fg          
 function(x) {  # x は二次元ベクトル変数 (x[1],x[2])
     gr <- function(x1, x2) { # fg の偏微分(gradient)ベクトル関数
         c(-400*x1*(x2 - x1*x1)-2*(1-x1), 200*(x2 - x1*x1))
     }
     x1 <- x[1]; x2 <- x[2]
     res<- 100*(x2 - x1*x1)^2 + (1-x1)^2 # fg の本体
     attr(res, "gradient") <- gr(x1, x2) 
     return(res)
  }

* 例４ (例３と同じ、gradient を数式微分で計算する例)
fdd は gradient 関数そのものではなく、gradient 属性を備えた目的関数自身であることを注意。

> nlm4 <- function () { 
     fd <- deriv(~ 100*(x2 - x1*x1)^2 + (1-x1)^2, c("x1", "x2")) #  目的関数の数式微分
     fdd <- function(x1, x2) {}  # 空の関数定義
     body(fdd) <- fd             # 関数の本体を fd とする
     nlm(function(x) fdd(x[1], x[2]), c(-1.2,1), hessian=TRUE) 
 }
 > nlm4.f <- nlm4()
 > nlm4.f
 $minimum
 [1] 1.182096e-20
 $estimate
 [1] 1 1
 $gradient
 [1]  2.583521e-09 -1.201128e-09
 $hessian
         [,1]    [,2]
 [1,]  802.24 -400.02
 [2,] -400.02  200.00
 $code
 [1] 1
 $iterations
 [1] 24

> str(nlm4.f)
 List of 6
  $ minimum   : num 1.18e-20
  $ estimate  : num [1:2] 1 1
  $ gradient  : num [1:2]  2.58e-09 -1.20e-09
  $ hessian   : num [1:2, 1:2]  802 -400 -400  200
  $ code      : int 1
  $ iterations: int 24

# fdd は gradient 関数そのものではなく、gradient 属性を備えた目的関数自身であることを注意。
 # (オブジェクト名の先頭のドットは紛らわしいオブジェクト名を避けるため)
 > fdd
 function (x1, x2)
 {
     .expr2 <- x2 - x1 * x1
     .expr5 <- 1 - x1
     .value <- 100 * .expr2^2 + .expr5^2
     .grad <- array(0, c(length(.value), 2), list(NULL, c("x1",
         "x2")))
     .grad[, "x1"] <- -(2 * .expr5 + 100 * (2 * ((x1 + x1) * .expr2)))
     .grad[, "x2"] <- 100 * (2 * .expr2)
     attr(.value, "gradient") <- .grad
     .value
 }

* 例５ demo(nlm) より (gradient と hessian をともに式で与える(失敗する)例)
hessian を陽に与えてもそれだけ賢くなるわけではない！？ 数値微分の値は正確な値とは当然微妙にことなるはずなので、こういうこともおき得る！？

これはhessianの定義が不完全(もしくはmatrixでの行と列の入れ間違い)だからでしょうか？
単なる誤植だったとすると、微分を明確に与えることによって数値微分よりも
曲面が複雑になってしまったからでしょうか？

hessian の定義は正しいですね(どのみち対称ですから)。失敗する理由は良くわかりませんが、真の解 (1,1) における hessian 行列の固有値は 1001.6 と 0.39936 と極めてアンバランス(ちょうどデスバレー風)なのが原因では？同じ変位でも方向により関数値が極端に変わります。

> fgh <- function(x)
   {
     gr <- function(x1, x2) {
         c(-400*x1*(x2 - x1*x1) - 2*(1-x1), 200*(x2 - x1*x1))
     }
     h <- function(x1, x2) {
         a11 <- 2 - 400*(x2 - x1*x1) + 800*x1*x1
         a21 <- -400*x1
         matrix(c(a11,a21,a21,200),2,2)
     }
     x1 <- x[1]; x2 <- x[2]
     res<- 100*(x2 - x1*x1)^2 + (1-x1)^2
     attr(res, "gradient") <- gr(x1, x2) # fgh に gradient 属性として gr を与える
     attr(res, "hessian") <- h(x1, x2)   # fgh に hessian 属性として h を与える
     return(res)
 }

> nlm5.f <- nlm(fgh, c(-1.2,1), hessian=TRUE)
 > nlm5.f
 $minimum        # 完全な失敗！
 [1] 2.829175
 $estimate
 [1] -0.6786981  0.4711891
 $gradient
 [1] -0.4911201  2.1115987
 $hessian
          [,1]     [,2]
 [1,] 366.1188 271.4593
 [2,] 271.4593 200.0000
 $code           # 失敗！
 [1] 4
 $iterations     # 最大繰返し回数まで実行されている！ 
 [1] 100

> eigen(nlm5.f$hessian)  # hessian は正定値符合行列でない！
 $values
 [1] 566.9414302  -0.8225813
 $vectors
           [,1]       [,2]
 [1,] 0.8039230 -0.5947334
 [2,] 0.5947334  0.8039230

> fgh               # 関数 fgh の中身は
 function(x)
 {
     gr <- function(x1, x2) {
         c(-400*x1*(x2 - x1*x1) - 2*(1-x1), 200*(x2 - x1*x1))
     }
     h <- function(x1, x2) {
         a11 <- 2 - 400*(x2 - x1*x1) + 800*x1*x1
         a21 <- -400*x1
         matrix(c(a11,a21,a21,200),2,2)
     }
     x1 <- x[1]; x2 <- x[2]
     res<- 100*(x2 - x1*x1)^2 + (1-x1)^2
     attr(res, "gradient") <- gr(x1, x2)
     attr(res, "hessian") <- h(x1, x2)
     return(res)
 }

# 反復回数を増やせばそれなりの結果が得られる
 > nlm(fgh, c(-1.2,1), iterlim=10000, hessian=TRUE)
 $minimum
 [1] 1.675511e-07
 $estimate
 [1] 0.9995911 0.9991805
 $gradient
 [1] -0.0000996657 -0.0003592498
 $hessian
           [,1]      [,2]
 [1,]  801.5865 -399.8564
 [2,] -399.8564  200.0000
 $code
 [1] 2
 $iterations
 [1] 5044

* 例６ 数式微分関数 deriv3 でgradient と hessian を同時に与える

> fd <- deriv3(~ 100*(x2 - x1*x1)^2 + (1-x1)^2, c("x1", "x2")) #  関数の数式微分
 > fdd <- function(x1, x2) {}  # 空の関数定義
 > body(fdd) <- fd             # 関数の本体を fd とする
 
 > fdd   # fdd の中身 (gradient と hessian がともに属性として備わっている)
 function (x1, x2)
 {
     .expr2 <- x2 - x1 * x1
     .expr5 <- 1 - x1
     .expr9 <- x1 + x1
     .value <- 100 * .expr2^2 + .expr5^2
     .grad <- array(0, c(length(.value), 2), list(NULL, c("x1",
         "x2")))
     .hessian <- array(0, c(length(.value), 2, 2), list(NULL,
         c("x1", "x2"), c("x1", "x2")))
     .grad[, "x1"] <- -(2 * .expr5 + 100 * (2 * (.expr9 * .expr2)))
     .hessian[, "x1", "x1"] <- -(100 * (2 * ((1 + 1) * .expr2 -
         .expr9 * .expr9)) - 2)
     .hessian[, "x1", "x2"] <- .hessian[, "x2", "x1"] <- -(100 *
         (2 * .expr9))
     .grad[, "x2"] <- 100 * (2 * .expr2)
     .hessian[, "x2", "x2"] <- 100 * 2
     attr(.value, "gradient") <- .grad
     attr(.value, "hessian") <- .hessian
     .value
 }

* 例７ example(nlm) より (追加引数つきの目的関数の例)

> f <- function(x, a) {  # ベクトル化関数 (変数次元は未設定、a はパラメータ)
    res <- sum((x - a)^2)
    attr(res, "gradient") <- 2 * (x - a)
    res
 }

# 目的関数は二変数関数 (x[1]-3)^2+(x[2]-5)^2 とされる。
 > nlm(f, c(10, 10), a = c(3, 5)) # 二次元引数 x、a はパラメータ
 $minimum  # 当然
 [1] 0
 $estimate # =a、当然
 [1] 3 5
 $gradient
 [1] 0 0
 $code
 [1] 1
 $iterations 
 [1] 1

タイムスタンプを変更しない