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

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

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

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

-COLOR(magenta){参考文献} K. Lange. Numerical Analysis for Statisticians. Springer 2001, p185ff

-COLOR(magenta){例}
--COLOR(magenta){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

---(無限)矩形型の制約 COLOR(red){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

---領域 COLOR(red){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

--COLOR(magenta){線形、2次計画法問題} 適切な初期値が必要。パッケージ COLOR(magenta){quadprog} 中の例 COLOR(red){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

トップ   編集 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS