Rで黄金分割探索

ある化石的な最適化アルゴリズムを書くための準備*1なのですが、黄金分割探索をRで実装してみましょう。計算するだけならばoptimize関数を関数を用いればよいのですが、建築や美術など見た目以外でも黄金比が実用的に使えるのを実感できます。

黄金分割探索は、微分法を使わずに一変数の最適化問題を解くことができるもので、引数の有限探索区間を3分割し、内側の2箇所の区間境界における目的関数の値に応じて、左右どちらかの区間1つを探索区間から除く作業を、収束条件を満たすまで繰り返すアルゴリズムで、探索区間の分割を黄金比にするところに名前の由来があります。数値演算の教科書にはよく載っていますし、Wikipediaあたりには実装に必要な量の説明がされているので、図つきの詳しい説明はそちらを見てください。

#
# 黄金分割探索（目的関数を極大化する引数を0から1までの範囲で探す）
#
gss <- function(vf, min=0, max=1, maxit=30, steptol=1e-6, maximize=FALSE){

# a:=p2-p0, b:=p1-p2, c:=p3-p2と置いて、b/a, c/aが黄金比になるように、p0 < p2 < p4 < p1の初期値を定める。
  gr <- (1 + sqrt(5))/2 # 黄金比
  p0 <- min
  p1 <- max
  l <- (p1 - p0)
  p2 <- p0 + l*1/(1 + gr)
  p3 <- p2 + l*1/(1 + gr)/(1 + gr)

# 初期値における目的関数の値を計算
  v2 <- vf(p2)
  v3 <- vf(p3)

  it <- 0 # 繰り返し回数を初期化

# 終了条件は目的関数の引数が変化しなくなったときと、繰り返し回数が上限に達したとき
  while(steptol < p1 - p0 && it < maxit){

    # 目的関数の大小関係を示す差分を計算
    # 最大化のときは符号を逆にする
    v_diff <- (v2 - v3)*(1 - 2*maximize)

    # 探索区間を縮小し、p1,p2,p3,p4,v2,v3の値を更新する
    if(0 > v_diff){
      p1 <- p3
      p3 <- p2
      v3 <- v2
      l <- (p1 - p0)
      p2 <- p0 + l*1/(1 + gr)
      v2 <- vf(p2)
    } else {
      p0 <- p2
      p2 <- p3
      v2 <- v3
      l <- (p1 - p0)
      p3 <- p2 + l*1/(1 + gr)/(1 + gr)
      v3 <- vf(p3)
    }

    it <- it + 1
  }

# 収束時はp1=p2=p3=p4, 収束しないことは考えない
  p1
}

試しに二次関数の最小化/最大化問題を解いてみましょう。

gss(function(x){
  (x - 0.23)^2
})

gss(function(x){
  -1*(x - 0.61)^2
}, min=0.5, max=1, maxit=30, maximize=TRUE)

そこそこループするので速度的には微妙な気もしますが、0.2300004と0.6100005が計算できます。

*1:準ニュートン法と呼ばれる多変数の最適化アルゴリズム群において、一変数の最適化問題を求めたいところが出てきます。目的関数が与えられたときのことを考えましょう。もっとも基本的なものはニュートン・ラフソン法になります。真の係数を、推定量を（は反復回数で指数ではない）とすると、 になります。はステップ幅です。このヘッセ行列の部分があるニュートン・ラフソン法ではでよいのですが、二階部分の計算があると上手く計算されないことが多々あり、一階微分しかないフィッシャー情報行列の符号を反転させたもので代用したり、一回計算したら簡素な更新を続けて使いまわしたりします。代用する場合、移動量が十分に調整されないわけで、ステップ幅を毎回調整しないと上手く計算できないようです。 ステップ幅自体は一定の条件を満たせば何でもよいのですが、説明用にあまり雑なコードを書いているとあれこれ言われてしまいます。ステップ幅の計算のために、黄金分割探索のソースコードを準備します。