線形不等式制約付きの最適化関数 constrOptim のバックアップ(No.1)

• バックアップ一覧
• 差分 を表示
• 現在との差分 を表示
• ソース を表示
• 線形不等式制約付きの最適化関数 constrOptim へ行く。
  • 0 (1970-01-01 (木) 18:00:00)
  • 1 (2015-03-01 (日) 01:15:59)

線形不等式制約付きの最適化関数 constrOptim()

constrOptim 関数は adaptive barrier アルゴリズム(制約領域の境界に近付くと目的関数値が増加するような付加項を適宜付け加える?)を用い、線形不等式制約(つまり無限もしくは有界な多角形内での最小値求める)下で関数を最小化する。内部的に optim() 関数を使う。矩形制約付きの最適化は otpim の method="L-BGFS-B" が使える。

• 書式 (method オプションの既定値の与え方にびっくり！)

    constrOptim(theta, f, grad, ui, ci, mu = 1e-04, control = list(),
                method = if(is.null(grad)) "Nelder-Mead" else "BFGS",
                outer.iterations = 100, outer.eps = 1e-05, ...)

• 引数
  • theta 初期値。制約式を満たす範囲内にとる必要。
  • f 最小化すべき関数
  • grad f のグラディエント関数
  • ui 制約(以下を参照)
  • ci 制約(以下を参照)
  • mu (小さい)調整用パラメータ
  • control optim() に引き渡されるパラメータ
  • method optim() に引き渡されるパラメータ
  • outer.iterations barrier アルゴリズムの繰返し回数
  • outer.eps barrier アルゴリズムの相対的収束判定用の判定値
  • ... optim に引き渡されるその他のパラメータ

• 詳細
  • 制約範囲は ui %*% theta - ci >= 0 で定義される領域。初期値はこの領域の内部にある必要がある。しかし最適値は境界上にある可能性がある。
  • 制約条件を強調するため対数的なバリヤが加えられ、それから optim が適用される。バリヤ関数は、繰返しの度に目的関数が減少するように選ばれる。制約領域内の最小値は典型的にはかなり早く見つかるが、境界上の最小値に対してはかなりの繰返しが必要になる。
  • 調整用のパラメータ mu はバリヤ項に掛けられる。その正確な値はしばしば重要ではない。mu が増大すると、変形された目的関数は本来のそれにより近くなるが、同時に制約領域の境界の近くでより滑らかでなくなる。
  • 目的関数が無限大値を取ることを許す任意の optim の手法が使える(現在のところ L-BFGS-B 以外のすべて)。method="Nelder-Mead を除くすべてのメソッドではグラディエントを与える必要がある。
  • optim と同様、既定では最小化を行なう。最大化するにはパラメータ control$fnscale を負の値にする。

• 返り値 optim と同じだが、追加の二つの情報を含む。barrier.value は最適値でのバリヤ関数の値、outer.iterations は外側の繰返し数(optim の呼出し回数)。

• 参考文献 K. Lange. Numerical Analysis for Statisticians. Springer 2001, p185ff

• 例
  • Rosenbrock の Banana 関数

• 制約無しの最小化

   > fr <- function(x) {   ## Rosenbrock Banana function
             x1 <- x[1]
             x2 <- x[2]
             100 * (x2 - x1 * x1)^2 + (1 - x1)^2
            }
   > grr <- function(x) { ## fr のグラディェント関数
              x1 <- x[1]
              x2 <- x[2]
              c(-400 * x1 * (x2 - x1 * x1) - 2 * (1 - x1),
                200 *      (x2 - x1 * x1))
            }
   > optim(c(-1.2,1), fr, grr) # 制約無しの最小化
   $par                        # 真の解は 1,1 
   [1] 1.000260 1.000506
   $value                      # 真の最小値は 0
   [1] 8.825241e-08
   $counts
   function gradient
        195       NA
   $convergence
   [1] 0
   $message
   NULL

• (無限)矩形型の制約 x1<=1, x2<=1、最適値は境界上にある

   > constrOptim(c(-1.2,0.9), fr, grr, ui=rbind(c(-1,0),c(0,-1)), ci=c(-1,-1))
   $par                     # 最小値を与えるパラメータ値(矩形境界線上)
   [1] 0.9999761 0.9999522  # ほぼ境界線上
   $value                   # 最小値
   [1] 5.708627e-10         # 不思議なことに制約無しの場合よりも零に近い！
   $counts                  # 関数評価回数
   function gradient
          6        1
   $convergence             # 完全な収束   
   [1] 0
   $message
   NULL
   $outer.iterations        # 内部での optim 関数の呼出し回数 
   [1] 14
   $barrier.value           # 収束点でのバリヤ関数の値(零に近い方が良い?)
   [1] -0.0001999198

• 領域 x1<=0.9, x1-x2>=0.1 での最適化(大局解はこの中に無い)

   > constrOptim(c(.5,0), fr, grr, ui=rbind(c(-1,0),c(1,-1)), ci=c(-0.9,0.1))
   $par
   [1] 0.8891335 0.7891335  # ほぼ境界上?
   $value
   [1] 0.01249441
   $counts
   function gradient
         18        1
   $convergence
   [1] 0
   $message
   NULL
   $outer.iterations
   [1] 5
   $barrier.value
   [1] -7.399944e-05

• 線形、2次計画法問題 適切な初期値が必要。パッケージ quadprog 中の例 example(solve.QP)。
  • グラディエントを使わない。

    > fQP <- function(b) {-sum(c(0,5,0)*b)+0.5*sum(b*b)}
    > Amat <- matrix(c(-4,-3,0,2,1,0,0,-2,1),3,3)
    > bvec <- c(-8,2,0)
    > constrOptim(c(2,-1,-1), fQP, NULL, ui=t(Amat),ci=bvec)
    $par
    [1] 0.4761374 1.0477253 2.0954507
    $value
    [1] -2.380952
    $counts
    function gradient
         252       NA
    $convergence
    [1] 0
    $message
    NULL
    $outer.iterations
    [1] 5
    $barrier.value
    [1] -0.0006243786

• グラディエントを使うと

  > gQP <- function(b) {-c(0,5,0)+b}
  > constrOptim(c(2,-1,-1), fQP, gQP, ui=t(Amat), ci=bvec)
  $par
  [1] 0.4761908 1.0476188 2.0952376
  $value
  [1] -2.380952
  $counts
  function gradient  # グラディエントを与えた分、呼出し回数が劇的に減少している
        21        1
  $convergence
  [1] 0
  $message
  NULL
  $outer.iterations
  [1] 5
  $barrier.value
  [1] -0.0006243894

• 最大化を行なう例

  > hQP <- function(b) {sum(c(0,5,0)*b)-0.5*sum(b*b)}
  > constrOptim(c(2,-1,-1), hQP, NULL, ui=t(Amat), ci=bvec,
                control=list(fnscale=-1))
  $par
  [1] 0.4763538 1.0477181 2.0954838
  $value
  [1] 2.380751
  $counts
  function gradient
       300       NA
  $convergence
  [1] 0
  $message
  NULL
  $outer.iterations
  [1] 1
  $barrier.value
  [1] -0.001142111