「Rによるモンテカルロ法入門」読書ノート：5章モンテカルロ最適化確率的勾配法

読書ノートアーカイブ：http://yagays.github.com/blog/2012/10/20/archive-introducing-monte-carlo-methods-with-r/

内容：確率的勾配法(Stochastic Gradient Methods)で関数の最大値を求める

目的関数からの直接的なシミュレーションが難しい場合，関数の表面をローカルになぞって最大値を探索するという手法が使える．探索に用いる系列は以下のように表す事ができる．

$\theta_{j+1} = \theta_j + \varepsilon_j$

これは現在のステップに摂動項 $\varepsilon_j$ を加えて次のステップとすることを表している．つまりこれはマルコフ連鎖となるのだが，今回はそれほどマルコフ性は重要にはならない．

概念的に目的関数の分布を山の斜面と考えるならば，現在の位置からちょっと動いて…を繰り返して山を登っていくような動きをする．ただし，単純に傾斜の上の方に登っていくのではなく，どちらかに進むかはランダムなので，たまに斜面を下ったりもする．その試行を繰り返していくことで，山の頂上を目指すような方法となっている．

有限差分法

有限差分法では，関数の勾配は

$\nabla h(\theta_j) \approx \frac{h(\theta_j+\beta_j\zeta_j)-h(\theta_j-\beta_j\zeta_j)}{2\beta_j}\zeta_j = \frac{\Delta h(\theta_j,\beta_j\zeta_j)}{2\beta_j}\zeta_j$

とし，摂動項 $\varepsilon_j$ を組み込んだ系列の更新は以のようになる．

$\theta_{j+1} = \theta_j + \frac{\alpha_j}{2\beta_j}\Delta h(\theta_j,\beta_j\zeta_j)\zeta_j$

例 5.7 有限差分法による関数の最大化において，複数の勾配パラメータがどのように収束に影響するかを調べる

例5.6の続き．有限差分法を用いて以下の式を最小化することを目的とする．

$h(x,y)=(x\sin(20y)+y\sin(20x))^2\cosh(\sin(10x)x)+ (x\cos(10y)-y\sin(10x))^2\cosh(\cos(20y)y)$

今回は最大化ではなく最小化することが目的なので，この関数をマイナスにする必要がある．

h <- function(x){
  -(x[1]*sin(20*x[2]) + x[2]*sin(20*x[1]))^2 * cosh(sin(10*x[1])*x[1]) -
    (x[1]*cos(10*x[2]) - x[2]*sin(10*x[1]))^2 * cosh(cos(20*x[2])*x[2])
}

アルゴリズム

勾配の値が10^-5 になるまで以下を繰り返す
1. $\zeta_j$ を生成する
2. $\theta_{j}$ を使って勾配 $\frac{\alpha_j}{2\beta_j}\Delta h(\theta_j,\beta_j\zeta_j)\zeta_j$ を生成する
3. $\theta_{j}+\varepsilon_j$ を $\theta_{j+1}$ とする
4. $\alpha_j$ と $\beta_j$ を更新する

なお，実際のソースコードでは勾配の値が発散しないように，scaleが1以上になった場合には内部でループをして値を再計算するようにしている．

$\zeta_j$ の生成方法と有限差分法における役割

さて，ここで有限差分法のポイントとなる $\zeta_j$ に関して，実際にどうやって値を生成するのか，そしてこの値がどういう意味を持っているかについて考えてみる．

まずは $\zeta_j$ の生成方法について．本文中では「 $\zeta_j$ は単位球 $\| \zeta_j \| = 1$ 上に一様に分布しています」とあるだけで，コードの方を参照してみてもzeta / sqrt(t(zeta) %*% zeta)としか書かれていない．ではこれは一体何をしているのかというと，一言で言うと「球面上の一様分布から乱数を生成する」，つまり今回の場合は単位円の円周の上からある一点を取ってきて $\zeta_j$ の値にするということをしている．

このような球面上の一様分布に関しては「Ｒによる計算機統計学」のP.96に詳しく書かれており，この本ではD次元の球面上の一様乱数を生成する方法についてrunif.sphere()関数を定義している．今回の場合はxとyの2次元について考えているので，「Ｒによる計算機統計学」に倣ってrunif.sphere(1,2)とすれば， $\zeta_j$ と同様の値が得られる．

# Ｒによる計算機統計学 P.87より
# http://personal.bgsu.edu/~mrizzo/SCR/SCRch3.R
runif.sphere <- function(n, d) {
    # return a random sample uniformly distributed
    # on the unit sphere in R ^d
    M <- matrix(rnorm(n*d), nrow = n, ncol = d)
    L <- apply(M, MARGIN = 1,
               FUN = function(x){sqrt(sum(x*x))})
    D <- diag(1 / L)
    U <- D %*% M
    U
}

次に，有限差分法における $\zeta_j$ の役割について考えてみる．

そもそも，確率的勾配法はローカルに関数の表面をなぞるように探索して，値の大きいポイントを探す確率的なアプローチだった．概念的には，前回3Dプロットをしたような山みたいになっている関数の上をあちこち動いて頂上を探し出すということだが，あちこち動きまわるためには，歩く「方向」と「距離」を考えなければならない．今回の有限差分法における勾配は $\varepsilon_j = \frac{\alpha_j}{2\beta_j}\Delta h(\theta_j,\beta_j\zeta_j)\zeta_j$ で表されていたが，その「方向」が $\zeta_j$ ，「距離」が $\frac{\alpha_j}{2\beta_j}\Delta h(\theta_j,\beta_j\zeta_j)$ に対応している．

では実際に $\zeta_j$ と勾配 $\varepsilon_j$ を図示して，「方向」と「距離」の意味を確かめてみよう．以下の図は， $\zeta_j$ の取りうる単位円と，赤い矢印 $\zeta_j$ ，黒い矢印 $\varepsilon_j$ を示した図になっている（ $j=1$ ）．この図を見ると，赤い矢印で方向が示され，黒い矢印でどれくらいの距離を進むかが見て取れる． $\zeta_j$ は乱数で生成された値なので，本書の「そのつど方向をランダムに選択する」という勾配の説明の通り，単位円の円周上からどこか一点を選択して進む方向を決めている．

といったように，この有限差分法では関数の山の傾斜を考えて進む方向を決めるのではなく，何も考えずに動いてみるといったアプローチとなっている．ただし，距離に関しては別で，関数を使って傾斜を見つつ調整しているので，その点は注意．

plot(runif.sphere(1000,2), xlim=c(-2,2), ylim=c(-2,2), cex=0.1)
arrows(0,0,grad[1],grad[2], lwd=2)             # 黒の矢印が勾配
arrows(0,0,zeta[1],zeta[2], lwd=2, col="red")  # 赤の矢印がゼータ

今回実験を行う4つのシナリオ

この有限差分法において， $\alpha_j$ と $\beta_j$ の系列の更新式は任意に設定できる．今回は4つのシナリオを用意して実験を行った．それぞれのシナリオの更新式は以下の通り．

シナリオ1： $1/\log(j+1)$ ， $1/\log(j+1)^{0.1}$
シナリオ2： $1/100\log(j+1)$ ， $1/\log(j+1)^{0.1}$
シナリオ3： $1/(j+1)$ ， $1/(j+1)^{0.5}$
シナリオ4： $1/(j+1)$ ， $1/(j+1)^{0.1}$

(※シナリオ1と2に関しては後述の「本問題の疑問点」を参考)

例 5.7のソースコード

前置きがだいぶ長くなってしまったが，本題のソースコードおよび4つのシナリオにおける確率的勾配法の探索過程を示した図は以下のようになる．

h <- function(x){
  -(x[1]*sin(20*x[2]) + x[2]*sin(20*x[1]))^2 * cosh(sin(10*x[1])*x[1]) -
    (x[1]*cos(10*x[2]) - x[2]*sin(10*x[1]))^2 * cosh(cos(20*x[2])*x[2])
}

start <- c(0.65, 0.8)
theta <- matrix(start, ncol=2)
dif <- 1
iter <- 1
alpha <- 1
beta <- 1

while(dif > 10^-5){
  zeta <- rnorm(2)
  zeta <- zeta / sqrt(t(zeta) %*% zeta)
  grad <- alpha[iter]*zeta * (h(theta[iter,]+beta[iter]*zeta) -
    h(theta[iter,]-beta[iter]*zeta)) / 2*beta[iter]
  scale <- sqrt(t(grad) %*% grad)
  while(scale > 1){
    zeta <- rnorm(2)
    zeta <- zeta / sqrt(t(zeta) %*% zeta)
    grad <- alpha[iter]*zeta * (h(theta[iter,]+beta[iter]*zeta) -
      h(theta[iter,]-beta[iter]*zeta)) / 2*beta[iter]
    scale <- sqrt(t(grad) %*% grad)
  }
  theta <- rbind(theta, theta[iter,]+grad)
  dif <- scale
  iter <- iter + 1

# scenario1
  alpha <- cbind(alpha, 1/log(iter+1))
  beta <- cbind(beta, 1/(log(iter+1))^(0.1))
}

hcur <- h(start)
hval <- h(start)
x <- seq(-1,1,le=435)
y <- seq(-1,1,le=435)
z <- matrix(0, nrow=435, ncol=435)
for (i in 1:435){
  for (j in 1:435){
    z[i,j] <- h(c(x[i], y[j]))
  }
}

image(x, y, z, col=terrain.colors(150))
lines(theta, lwd=2)
points(theta[1,1], theta[1,2], col="gold", pch=19)
title(main="scenario 1")

シナリオ2,3,4の場合は，whileループの最後の $\alpha_j$ と $\beta_j$ の更新式を以下の通りにする．

# scenario2
  alpha <- cbind(alpha, 1/(100*log(iter+1)))
  beta <- cbind(beta, 1/(log(iter+1))^(0.1))
# scenario3
  alpha <- cbind(alpha, 1/(iter+1))
  beta <- cbind(beta, 1/sqrt(iter+1))
# scenario4
  alpha <- cbind(alpha, 1/(iter+1))
  beta <- cbind(beta, 1/(iter+1)^(0.1))

上の図の解釈だが，この4つのシナリオの場合において

シナリオ1では，動きまわる距離が大きすぎて（ $\alpha_j$ が収束するのに時間がかかって）なかなかアルゴリズムが停止しない
シナリオ2では，動きまわる距離が小さすぎて（ $\alpha_j$ がすぐ収束してしまって）ローカルミニマムにはまっている
シナリオ3では，動きまわる距離はシナリオ2より改善しているが，最大値を見つけられていない
シナリオ4では，上手い具合に最大値を見つけることができている

ことが分かる．

なお，この解釈は本書（P.156）に書かれているものとは多少異なっているということを注意していただきたい．それに関しては次の疑問点のところで述べる．

本問題の疑問点

さて，この問題を素直に解いてもP.157の図5.7のようにはいかない．色々と試してみた結果，どうやら本書におけるシナリオ1と2の $\beta_j$ において，分母にかかっている0.1乗の表記がおかしい気がする．この0.1乗がlogの中の括弧にかかっていると考えて0.1*log()としてコードを組むと，値が発散しすぎてなかなか収束しなくなる．ということで，ここでは以下のように式を解釈して問題を解いた．正直なところ，これが正しいのかよくわかっていない．