ホロノミック勾配法勉強メモ

ホロノミック勾配法の勉強メモです。

ホロノミック勾配法のよく知られている応用として

・確率分布の正規化定数の評価

・何らかの最適化

がありますが、本記事は前者についての個人メモです。後者についてはデンソーITラボの
ホロノミック勾配法 – Tech.D-ITlab | Denso IT Laboratory researcher's blog sites
こちらにまとめがあります。

ここでいう正規化定数とは、例えば指数型分布族の形を思い出すと
$\displaystyle p(x|\theta) = \frac{\exp(\eta(\theta) \cdot T(x))}{Z(\theta)}$
のような形をしていて、この分母の
$\displaystyle Z(\theta) := \int \exp(\eta(\theta) \cdot T(x)) dx$
を指します。（ベイズ統計の正規化定数とは違うもの）

単なる積分評価なのでモンテカルロ法やらの数値積分でも評価は可能ですが、計算時間や精度の問題があるため、積分によっては別のアプローチの方が良いことがあり、ホロノミック勾配法もそのうちの1つのアプローチ方法です。

ホロノミック勾配法では大きくは以下の2ステップで $\displaystyle Z(\theta)$ を評価します：

ステップ1. $\displaystyle Z(\theta)$ に関する微分方程式系を作る

ステップ2. 上記の微分方程式系を数値解析で解く

まずステップ1.の微分方程式系を作れる（存在する）かという疑問だが、 $\displaystyle Z(\theta)$ が名前の由来になっている"ホロノミック"であれば可能であることが知られています。

ではどんなものがホロノミックであるかというと

a. 有理函数

b. exp(有理函数), log(有理函数)

c. holonomic + holonomic, holonomic * holonomic

d. $\displaystyle \int \rm{holonomic}$

などがホロノミックであるので、統計で使われているほとんどのクラスはホロノミックと言えそう。

ここで、 $\displaystyle Z(\theta)$ がホロノミックであればグレブナー基底とかで機械的に微分方程式系を作ることが可能なことが知られています。たとえば、Macaulay 2にD-加群のアルゴリズムが入っているDmodulesパッケージがあるのでそれを使えば良いらしい。ただし、次元が高かったりすると計算量が大きくてかなりきついようで、手計算で微分方程式が出せるのであればそれを使うに越したことはない。

ステップ2.については初期値を適当に決めて（例えば級数展開をして気合で計算して）、Runge-Kutta 法などの数値解析手法で解くのですが、ここの部分の初期値や数値解析手法で時間や精度が決まってくると思うので、実務上はこの辺りが肝になってきそうである。初期値を決めた後はRのhgmパッケージでも計算できます。

では具体的な例を見てみよう。
方向統計学と呼ばれる方向にも興味がある時に使う分布で、von Mises-Fisher分布というものがあるのでそれを使ってみる。ここでは3変数von Mises-Fisher分布（on $\displaystyle S^2$ ）の分子だけ知っていたとしよう。つまり、
$\displaystyle p(x|\mu, \tau) = \frac{\exp(\tau \mu^T x)}{Z(\tau)}$
として $\displaystyle Z(\tau)$ をホロノミック勾配法で数値計算して求めてみよう。
ここで $\displaystyle \mu, \tau$ はそれぞれmean directionおよびconcentration parameterと呼ばれるもので、特に $\displaystyle \tau$ は以下の図のように $\displaystyle \tau =0$ の時は全ての方向に一様で、 $\displaystyle \tau = \infty$ の時はデルタ関数となる。

f:id:lambdahat:20180611154704p:plain — von Mises-Fisher分布のconcentration parameter

（図は http://people.csail.mit.edu/jstraub/download/straub2017vonMisesFisherInference.pdf より引用）

まずステップ1.については適当に微分や積分をすると、正規化定数Zに関して以下のような微分方程式を得ることができる：
$\displaystyle \frac{d^2}{d \tau^2}Z(\tau) = Z(\tau) - \frac{2}{\tau}\frac{d}{d \tau} Z(\tau)$

ステップ2.について、まず初期値を決めるためにベッセル関数を雑に近似させたものを使った。その後、上記の微分方程式と初期値をRのhgmパッケージのhgm.Rhgm関数に入れた。（hgmパッケージの計算例と実質同じ例です）

install.packages("hgm")
library(hgm)

#微分方程式を記述
dZ.fun = function(tau, Z, fn.params=NULL){
dZ = array(0, c(1, 2))
dZ[1,1] = Z[2]
dZ[1,2] = Z[1] - 2*Z[2]/tau
return(dZ)
}

#Zの初期値を決めるためにターゲットを雑に近似
Z.asymp = function(tau){
c(2*pi*exp(tau)/tau, 2*pi*(tau-1)*exp(tau)/tau^2)
}

#tauの適当な初期値とターゲットとしているtauを指定
tau_ini = 1
tau_target = 10

#Zの初期値を計算
Z_ini = Z.asymp(tau_ini)

#ホロノミック勾配法によるZ(10)の近似値
hgm.Rhgm(tau_ini, Z_ini, tau_target, dZ.fun)[1]

上記で計算された結果を見ると

[1] 13839.66

という数値が出てくる。
これがどれだけ正しいかを見るために、以下正確な値を計算しよう。

Z.exact = function(tau){
  c( 4*pi*(sinh(tau)/tau), 4*pi*(cosh(tau)/tau - sinh(tau)/tau^2) )
}

Z.exact(tau_target)[1]

この結果を見ると

[1] 13839.64

となっているので初期値が雑でも結構正確な値となっていることが分かる。

最後にこの方法が特異モデルにも使えるかが気になるところだが、微分方程式に特異点がある場合は何らかの対処が必要なようだ。ただ代数解析とも関係があるのでそのあたりもうまくやれそうな気もするが、かなりまじめに勉強をしないと分からなそうだということだけ分かった。

参考文献：
[1] JST CREST 日比チーム (2011), グレブナー道場, 共立出版.
[2] 日比, 竹村, 原, 東口, 清, グレブナー道場著者一同 (2015), グレブナー教室, 共立出版.

あざらしの日常

あざらしが腹を叩く片手間で書いたメモ

ホロノミック勾配法勉強メモ