局所的な近似によるMCMCの高速化論文を読んだ

局所的な近似によるMCMCの高速化論文が話題になりました。
http://japan.zdnet.com/article/35073667/

論文Accelerating Asymptotically Exact MCMC for Computationally Intensive Models via Local Approximations
を読んで理解した範囲の内容を書きます。

アイデア

物理現象のモデルは(確率)微分方程式を使って書かれますが、そのパラメータが未知のことが多く、MCMCを使ってパラメータを推定することが行われています。
しかし現実的な問題ではモデルが巨大になって、MCMCの各ステップでそれを評価するのは計算量が多くなりすぎてしまいます。
そこでモデルを近似するということなのですが、サンプリングに近似した事後分布を使うというアイデアはApproximate Bayesian Computation(ABC)と共通しています。局所的に近似を行うことで高速に計算するというのがこの論文のテーマになっています。このアイデア自体は以前からあり、また尤度の式を変えてしまうのことに成るので正しいサンプリング、推定が出来るのかが疑問ですが漸近的(サンプル数無限大の極限)に元の分布を再現できることを示したのがこの論文での成果のようです(Appendix B)。
Metropolis-Hasting法では前のサンプル値に近いところから次のサンプル値をとってくるので局所的な近似で大丈夫そうなのは何となく理解できます。サンプル値を用いていることからモデルに確率変数が含まれるような場合にも対応できるそうです。*1

方法の概要

ベイズの定理に従い事前分布p,尤度L、モデルf, 観測データdに対し事後分布は
p(θ|d)∝L(θ|d,f)p(θ)
と書かれます。
観測ノイズを分布Pηに従う確率変数ηでモデル化したような場合観測データdは
d=f(θ)+η
となり、尤度L(θ|d,f)は
L(θ|d,f) =Pη(d−f(θ ))
と書かれます。
Metropolis-Hasting法を以下のような手順に置き換えます(Algorithm 1)。
pythonぽく書くと以下のようになります。

def RunChain(θ[0],S[0],L,d,p,f,l,T):
     for t in range(T):
          θ[t+1],S[t+1]=K_t(θ[t],S[t],L,d,p,f,l)
     return θ,S

def K(θ,S,L,d,p,f,l):
     θ_ ~l(θ,...)
     fp_ =approximate(f,θ_)
     fp=approximate(f,θ)
     a=min(1,L(θ_|d,fp_)p(θ_)/(L(θ|d,fp)p(θ)) )
     if(refinement is needed):
          θ_=θ__
          S=S+[(θ_,f(θ_)]
          K(θ_,S,L,d,p,f,l)
     else:
         u~Uniform(0,1)
          if(u<a):
               return (θ_,S)
          else:
               return (θ,S)

fp(論文では $\tilde{f}$ )を計算するところ(approximate)が普通のMCMCとの違いです。
より厳密で漸近的に事後分布に収束することが証明された手順はAlgorithm 4にあります(論文16 page)
refinementというのはサンプル値θが以下の近似で用いられるサンプル値の集合の中で偏った位置にあるときに取り直すことのようです。(2.3)

モデル近似の方法

ある時点でのサンプル値θの周辺の超球に含まれる以前のサンプリング値{θ^i}をもちいて近似モデルを作ります。
データ数N, モデルの次元dに対してM=(d+2)*(d+1)/2としたときのN*M行列
$\Phi = \left( \begin{array}{ccc} 1, \theta_1^1 & \cdots & \theta_d^1 & \frac{1}{2}(\theta_1^1)^2 & \cdots & \frac{1}{2}(\theta_d^1)^2 & \theta_1^1 \theta_2^1 & \cdots & \theta_{d-1}^1 \theta_d^1\\ \vdots & & & & \ddots & & & & \vdots \\ 1, \theta_1^N & \cdots & \theta_d^N & \frac{1}{2}(\theta_1^N)^2 & \cdots & \frac{1}{2}(\theta_d^N)^2 & \theta_1^N \theta_2^N & \cdots & \theta_{d-1}^N \theta_d^N\\ \end{array} \right)$

とサンプリング値に対する重み付け(2.2)の対角行列W, モデルfへθを代入した値y_i=f(θ_i) に対する最小二乗法
Φ^T W ΦZ=Φ^T W Y
を解いて得られるベクトル
$z_j^T=( a_j, b_j (H_j)_{1,1} \cdots (H_j)_{d,d}, (H_j)_{2,2} \cdots (H_j)_{d−1,d} )$
がfの二次近似
$\tilde{f}_j(\theta) =a_j + b_j\theta + \frac{1}{2}\theta^T H_j \theta$
の係数となります。θは正則化を行った後のものを使うようです。

既存手法との関係

Metropolis-Hasting法の中身を変えるような方法なので、その外側でHMC, NUTS, Riemannian manifold HMC, Metropolis-adjusted Langevin法などの方法と同時に使うことが出来ます。

使用例

論文では人工的なデータへの適用の他に双安定的(bistable)な大腸菌の遺伝子調節回路のパラメータ推定と拡散方程式の拡散係数推定が紹介されています。ナビエ・ストークス方程式のような一般の偏微分方程式にも使えるそうです。
以前このブログで微分方程式を解きながらパラメータを推定するデータ同化の方法について書きましたが、この論文ではMCMCの毎stepで微分方程式を解いて推定を行うような状況を想定しているようです。大腸菌の例では定常状態になったときの値のみを入力データとして使っています。

実装

内部はc++で書かれていてstanを取り込んで実装しているようです。pythonインターフェースがあるようです。しかしマニュアルは現時点ではほとんど何も書かれていないのでどう使うのかよくわかりません。stanのモデルを入れたりできるのでしょうか。
https://bitbucket.org/mituq/muq/src/30f954040dc87413a0e5cb52604361c7666f04b4/MUQ-Manual/UsersManual/?at=master
中の方を見るとパラレルテンパリングやRiemannian manifold HMCなどの手法も実装されているようです。
mituq / muq / source / MUQ / modules / Inference / src / MCMC — Bitbucket

https://bitbucket.org/mituq/muq/src/ab2b67cb6e8129a980bec99ce2e2e837f2c6e38a/MUQ/MUQ/Inference/MCMC/?at=master

Accelerating Bayesian inference in computationally expensive computer models using local and global approximations

*1:いっそ事後分布の推定をせず事前分布だけでサンプリングを行うという方法をとる場合もあったらしいのですが、これは流石に効率が悪そうです(Figure 1.(a) )