1. はじめに

関数の最適値を与える解や方程式の解を求めるためには，最急降下法のような微分を必要とする方法が必ずしも必要とは限りません．今回は，解の存在する区間を狭めていくことによって，解を求める最適化法についてまとめました．

2. 囲い込みによる求解法

2.1 二分法(Bisection search)

これは，$f(x) = 0$となるような$x$を求める方法です．この方法を適用するには，

区間$[a,b]$において連続
$f(a)f(b)<0$(端点において関数の符号が異なる)
区間$[a,b]$にただ一つの解が存在する．

という条件が成立している必要があります．2番目の条件は中間値の定理の条件にも出てくるものと同等です．二分法では，以下のようなアルゴリズムで区間を狭めていきます．

$x' \leftarrow \frac{a + b}{2}$
もし，$f(x') < 0$なら$a \leftarrow x'$，$f(x')>0$なら，$b \leftarrow x'$と更新する．
$|b-a|$が十分小さくなる，もしくは，事前に決めた絞り込み回数に到達したら終了する．

これを，雑に図解すると，

こんな感じ．反復のたびに，区間幅が$1/2$倍されていっていることがわかります．

2.2 三分割法(Ternary search)

先ほどの二分法は，区間の中心により区間を二分割したものでした．三分割法では，名前の通り，区間を三分割することで最適値(今回は最小値)の存在する区間を絞り込みます．三分割法を適用するための条件は，

区間$[a,b]$で$f$は連続
区間$[a,b]$で関数$f(x)$は単峰である．
区間$[a,b]$に解が存在する．

です．三分割法は，以下のようなアルゴリズムで実行されます．

$a^1\leftarrow a,\ b^1\leftarrow b$と初期化する
$i$回目の探索における区間$[a^i,b^i]$を次のように計算される$x_1^i,x_2^i$によって三分割する

\begin{align*}
\left\{\begin{array}{l}
x_1^i \leftarrow a^i + \frac{1}{3}(b^i-a^i)\\
x_2^i \leftarrow a^i + \frac{2}{3}(b^i-a^i)
\end{array}
\right.
\end{align*}

3.$i+1$回目の探索区間を，$f(x_1^i) \le f(x_2^i)$なら，$[a^{i+1},b^{i+1}]\leftarrow[a^i, x_2^i]$，$f(x_1^i) > f(x_2^i)$なら，$[a^{i+1},b^{i+1}]\leftarrow[x_1^i, b^i]$と更新する．
4.$|b^{i+1}-a^{i+1}|<\epsilon$もしくは，事前に決めた反復回数に到達したら終了する．そうでない場合は，2,3のステップを繰り返す．

これも雑に図解すると，

こんな感じ．反復回数ごとに，区間幅の$1/3$が切り捨てられていることがわかります．

2.3 黄金分割法(Golden section search)

三分割法では，区間の$1/3$を切り捨てた後，三分割するための内分点を再び計算していました．黄金分割法では，一つ前のステップの内分点の片方を再利用します．ざっくりとした考え方を図解すると，

こんな感じ．$x_1^1$が次のステップの$x_2^2$として再利用されています．黄金分割法では，区間を$\gamma:1-\gamma\ (0<\gamma<1/2)$の比率で内分します．この$\gamma$は，

\left\{
\begin{array}{l}
x_1^1-a^1 = \gamma(b^1-a^1)\\
x_1^1-a^1 = (1-\gamma)(x_2^1-a^1) = (1-\gamma)^2(b^1-a^1)
\end{array}\right.

という方程式$\gamma = (1-\gamma)^2$を解くことでもとまります．解いた結果が
$$\gamma = \frac{3-\sqrt{5}}{2}$$
です．黄金分割と呼ばれる由縁は，内分の比率
$$\frac{1-\gamma}{\gamma}$$が黄金比となるからです(おそらく，多分．考案者がジョジョ五部が好きとかではないと思う)．
黄金分割法は以下のようなアルゴリズムで実行されます．

$i$回目の区間$[a^i,b^i]$を内分する点$x_1^i,x_2^i$を$\gamma$を用いて，$x^i_1 = a^i + \gamma(b^i-a^i),\ x^i_2 = a^i + (1-\gamma)(b^i-a^i)$として求める．
$i+1$回目の探索区間を，もし，$f(x_1)^i\le f(x_2^i)$なら，

\left\{\begin{array}{l}
[a^{i+1},b^{i+1}]\leftarrow[a^{i},x_2^i]\\
x_2^{i+1}\leftarrow x_1^i\\
x_1^{i+1}\leftarrow a^{i+1} + \gamma(b^{i+1}-a^{i+1})
\end{array}\right.

と更新し，もし，$f(x_1)^i > f(x_2^i)$なら，

\left\{\begin{array}{l}
[a^{i+1},b^{i+1}]\leftarrow[x_1^{i},b^{i}]\\
x_1^{i+1}\leftarrow x_2^i\\
x_2^{i+1}\leftarrow a^{i+1} + (1-\gamma)(b^{i+1}-a^{i+1})
\end{array}\right.

と更新する．
3.収束するまで，2を反復する．

3. サンプルコード

今回は，Rによる実装例です．

r, bisection.R

bisection.method.eq0 = function(f, interval, niter = 100){
  upper = interval[2]
  lower = interval[1]
  for(i in 1:niter){
    x.mean = (upper + lower)/2
    if(abs(f(x.mean) - 0)<1e-6){
      return(x.mean)
    }
    if(sign(f(lower)) == sign(f(x.mean))){
      lower = x.mean
    }
    else{
      upper = x.mean
    }
    if(abs(upper - lower) < 1e-6){
      return(list(x.mean = x.mean, cnt = i))
    }
  }
  return(list(x.mean = x.mean, cnt = i))
}
f = function(x)return(x^3 -27)
interval1 = c(-10, 10)
bisection.method.eq0(f, interval1,1000000)
f.2 = function(x)return(5*x + 14)
bisection.method.eq0(f.2, interval1)

r, ternary.R

ternary.search = function(f, interval, niter = 100, tol = 1e-6){
  upper = interval[2]
  lower = interval[1]
  for(i in 1:niter){
    a.i = lower + (1/3)*(upper - lower)
    b.i = lower + (2/3)*(upper - lower)
    if(f(a.i) <= f(b.i)){
      upper = b.i
    }
    else{
      lower = a.i
    }
    if(abs(upper - lower) < tol){
      return(list(upper = upper, lower = lower, cnt = i))
    }
  }
  return(list(upper = upper, lower = lower, cnt = i))
}

f = function(x)return(x*(x-2)*(x+2))
interval2 = c(0,10)
ternary.search(f, interval2)

r, goldensection.R

golden.section.search = function(f, interval, niter = 100, tol = 1e-6){
  lower = interval[1]
  upper = interval[2]
  gam = (3-sqrt(5))/2
  a.i = lower + gam * (upper - lower)
  b.i = lower + (1 - gam) * (upper - lower)
  cnt = 0
  while(abs(upper - lower) > tol){
    if(f(a.i) <= f(b.i)){
      upper = b.i
      b.i = a.i
      a.i = lower + gam * (upper - lower)
    }
    else{
      lower = a.i
      a.i = b.i
      b.i = lower + (1-gam) * (upper - lower)
    }
    cnt = cnt + 1
    if(cnt == niter){
      return(list(lower = lower, upper = upper, cnt = cnt))
    }
  }
  return(list(lower = lower, upper = upper, cnt = cnt))
}

f = function(x)return(x*(x-2)*(x+2))
interval2 = c(0,10)
golden.section.search(f, interval2)

囲い込みによる非線形最適化入門(私が)