A Blog Entry on Bayesian Computation by an Applied Mathematician
$$
$$
概要
MALA とは
MALA は，提案核を対象分布 \(\pi\) から Langevin 拡散の形で構成する，problem-specific な Metropolis-Hastings 法であり，ランダムウォーク MH よりも速く収束することが期待される．(Corenflos and Finke, 2024) でも state-of-the-art と呼ばれているサンプリング法である．
(Fearnhead et al., 2018) において，MALA は BPS と比較されている．モデルは AR(1) を用いており，低次元ではほとんど変わらないが，高次元では BPS の方が自己相関時間が５倍良いという結論が得られている．
MALA の他に，(Neal, 1994) などの hybrid-Monte Carlo アルゴリズムも MALA と「提案核をよく考えている」という点で関係が深いが，本論文では扱わない．
MALA のエルゴード性について
Langevin 拡散 \[
dL_t=dB_t+\frac{1}{2}\nabla\log\pi(L_t)dt
\] は，例えば１次元では，尾部確率が指数減衰するならば，指数エルゴード性を持つ．つまり，\(\pi(x)\,\propto\,e^{-\gamma\lvert x\rvert^\beta}\;(\beta\ge1)\) など．
このように，確率分布の裾の重さに影響されるという現象が普遍的である．
続いて，Langevin 拡散の離散化が，元の Langevin 拡散に収束するかどうかも別問題である．
最後に，Langevin 拡散の離散化を Metropolis 法によって修正したアルゴリズムである MALA が指数収束する条件を示す．
導入
MCMC 手法の指数エルゴード性について
(Gelfand and Smith, 1990) 以来，MH 法の利用が爆発している
There has recently been a real explosion in the use of the Hastings and Metropolis algorithms
乱歩 MH などは対象分布に依らずに実装できる汎用アルゴリズムであるが，対象分布の情報を取り入れて収束を高速にすることが考えられる．その例が MALA である．
Langevin 動力学は極めて優秀なモデリング手法であるが，シミュレーションのための離散化が入ることで，収束は遅くなるどころかそもそも収束性が失われることさえある．
(Neal, 1994) を引用しながら，hybrid Monte Carlo が Langevin アルゴリズムに似ていると触れながら，話が対称性と非対称性の話に変わる．
However, it is worth remarking that often methods induced by non-reversible methods can be shown to converge more quickly than their reversible counterparts (see, for example, (Sheu, 1991)).
しかし (Sheu, 1991) は \(P_t\) の評価を上下からしているのみで，非対称性どころか収束の議論は全くしていない．議論に飛躍があるか，引用が間違っているかが考えられる．
Langevin 拡散について
本論文が研究する MH 法は，提案を Langevin 拡散により取り入れたものである．
密度 \(\pi\) が可微分であるとき，\(\nabla\log\pi\) が考えられる．これに対して，Langevin 拡散 \(\{L_t\}\subset\mathbb{R}^k\) とは \[
dL_t=dW_t+\frac{1}{2}\nabla\log\pi(L_t)dt
\] と定義される．
\(\pi\) が十分滑らかであるとき，これはエルゴード性を持つ： \[
\|P^t_L(x,-)-\pi\|_\mathrm{TV}\to0.
\]
本論文では，この収束が指数レートで起こる場合，そして同様の積率収束も起こる場合に興味がある．
指数収束する例
１次元クラス \(\mathcal{E}(\beta,\gamma)\)
\(k=1\) とする．\(\pi\in\mathcal{E}(\beta,\gamma)\subset C^\infty(\mathbb{R})\) であるとは，可微分であり，加えてある \(x_0\in\mathbb{R}\) と \(\gamma,\beta>0\) が存在して， \[
\pi(x)\,\propto\,e^{-\gamma\lvert x\rvert^\beta},\quad\lvert x\rvert\ge x_0,
\] が成り立つことをいう．このとき， \[
\nabla\log\pi(x)=-\gamma\beta\operatorname{sgn}(x)\lvert x\rvert^{\beta-1},\quad\lvert x\rvert>x_0.
\]
2.3 節にて，指数収束するための必要十分条件は \(\beta\ge1\) であることを見る．この結論は乱歩 MH アルゴリズムと全く同様である (Mengersen and Robert, 1996, p. 定理3.5)．
多次元指数分布族 \(\mathcal{P}_m\)
(G. O. Roberts and Tweedie, 1996) で乱歩 MH 法について考えたように，多次元の統計モデルとしては指数分布族を考えることとする．特に，次の場合を考える： \[
\pi(x)\,\propto\,e^{-p(x)}\quad \lvert x\rvert\ge x_0,
\] ただし \(p\) は \(m\) 次の多項式で，次の表示を持つ： \[
p=p_m+q_{m-1}.
\] \(\pi\in\mathcal{P}_m\subset C^\infty(\mathbb{R}^k)\) とは，\(p_m\xrightarrow{\lvert x\rvert\to\infty}\infty\) を満たすことをいう．特に，\(m\ge2\) であることに注意．
Langevin 拡散 \(L_t\) の離散近似
実際の実装で \(L_t\) をそのまま用いることは出来ず，これを離散化することが必要である．
ULA (unadjusted Langevin algorithm)
The unadjusted Langevin algorithm (ULA) is a discrete-time Markov chain \(U_n\) which is the natural discretization of ordinary Langevin diffusion \(L_t\).
(Parisi, 1980), (Grenander and Miller, 1994) で考えられたものである．
\[
U_n\sim\mathrm{N}_k\left(U_{n-1}+\frac{h}{2}\nabla\log\pi(U_{n-1}),hI_k\right)
\] によって構成できる． (Besag, 1994) が指摘したように，元の \(L_t\) と似た（しかしステップサイズ \(h\) に依存する変換を受けた）不変分布を持つ．
\(\pi=\mathrm{N}_1(0,1)\) であり，ステップサイズを \(h=2\) とすると，\(\nabla\log\pi(x)\)=-x$ であるから， \[
U_{n-1}+\frac{h}{2}\nabla\log\pi(U_{n-1})=U_{n-1}-\frac{2}{2}\frac{U_{n-1}}{1}=0
\] より，\(U_n\sim\mathrm{N}_1(0,2)\) となる．
MALA (Metropolis-Adjusted Langevin Algorithm)
(Besag, 1994) に従い，修正を加える．\(U_n\) を提案として，Metropolis-Hastings 法を実行するのである．
\[
q(M_{n-1},-):=\mathrm{N}_k\left(M_{n-1}+\frac{1}{2}h\nabla\log\pi(M_{n-1}),hI_k\right)
\] を提案核とし， \[
\alpha(M_{n-1},U_n):=1\land\frac{\pi(U_n)q(U_n,M_{n-1})}{\pi(M_{n-1})q(M_{n-1},U_n)}
\] を採択確率とする．
この場合，必ずエルゴード性が成り立つ．(G. O. Roberts and Tweedie, 1996) の結論と同様，\(\ell_k\)-既約かつ周期的だからである． 特に，殆ど至る所の \(x\in\mathbb{R}^k\) について， \[
\|P^n_M(x,-)-\pi\|_\mathrm{TV}\to0.
\] また，本論文が提示する指数収束の条件の下では，任意の \(x\in\mathbb{R}^k\) で起こる．
ULA が推移的であるとき，MALA も指数エルゴード性を持たないことが判る．だが，これは尾部確率が指数よりも重い場合にしか起こらない．
MALTA (Metropolis-Adjusted Langevin Truncated Algorithm)
\[
T_n\sim\mathrm{N}_k\biggr(M_{n-1}+R(M_{n-1}),hI_k\biggl)
\] \[
R(x):=\frac{D\nabla\log\pi(x)}{2(D\lor\lvert\nabla\log\pi(x)\rvert)}.
\] を提案とし，MH 法を適用する．
MALTA は大変ロバストで，指数収束が起こりやすい．
Langevin 拡散アルゴリズムの指数収束
一般の収束結果
\[
dL_t=dW_t+\frac{1}{2}\nabla\log\pi(L_t)dt
\]
\(\nabla\log\pi\) は \(C^1(\mathbb{R}^k)\)-級で， ある \(N,a,b>0\) が存在して \[
\nabla\log\pi(x)\cdot x\le a\lvert x\rvert^2+b,\qquad\lvert x\rvert>N,
\] を満たすとする．このとき，Langevin 拡散 \(\{L_t\}\) について次が成り立つ：
- \(L\) は爆発的でなく，\(\ell_k\)-既約で，周期的で，強 Feller であり，任意のコンパクト集合は小集合である．
- \(L\) は \(\pi\)-不変であり，任意の \(x\in\mathbb{R}^k\) について \[
\|P^t_L(x,-)-\pi\|_\mathrm{TV}\to0.
\]
\(\lvert L_t\rvert\) と OU 過程を比較することで，\(\lvert L_t\rvert\) が爆発しないことが判る．
続いて，ドリフト係数が局所有界であるため，(Bhattacharya, 1978) と同様の議論により強 Feller である．加えて，\(\ell_k\)-既約であることも判る．
Theorem 2.1. If, in addition to the hypothesis (A), \(a_{ij}(-)\) and \(b_i(-)\) are bounded on \(\mathbb{R}^k\), then for each \(x\) in \(\mathbb{R}^k\) there exists a unique probability measure \(P_x\) on \((\Omega',\mathcal{M}')\) such that (i) \(P_x(X_0=x)=1\), (ii) for every bounded real \(f\) on \(\mathbb{R}^k\) having bounded continuous first and second order derivatives, the process \(f(X_t)-\int^t_0Lf(X_s)ds\) is a martingale under \(P_x\). Further, (a) \(X\) is strong Markov and strong Feller, and (b) support of \(P_x\) is \(\Omega_x':=\left\{\omega\in\Omega\mid\omega(0)=x\right\}\). (Bhattacharya, 1978)
任意の \(\ell_k\)-正集合 \(A\subset\mathbb{R}^k\) に対して， \[
\widetilde{A}:=\left\{\omega\in\Omega_x'\mid \omega(t)\in A\right\}
\] は \(P_x\)-正集合である．実際，\(A\) が \(\mathbb{R}^k\) 上の開集合を含むために \(\widetilde{A}\) は \(\Omega_x'\) 上の開集合を含むので，\(P_x\)-零であったら \(P_x\) の台が \(\Omega_x'\) 全体であることに矛盾する．これより，\(L_t\) は \(A\) に到達可能である．よって \(L_t\) は \(\ell_k\)-既約．
（強）Feller かつ \(\ell_k\)-既約ならば，任意のコンパクト集合は小集合である (Tweedie, 1994, p. 定理5.1)．
非周期性は，任意のスケルトンも \(\ell_k\)-既約であることから従う．
(Ikeda and Watanabe, 1981, p. 5.4節) により，\(\pi\) は \(L\) の不変分布であることが，生成作用素 \[
A_Lf(x)=\left(\frac{1}{2}\nabla\log\pi(x)\right)\cdot\nabla f(x)+\frac{1}{2}\mathop{}\!\mathbin\bigtriangleup f(x)
\] \[
f\in C^2(\mathbb{R}^k)
\] に関して，\((\pi|A_Lf)=0\;(f\in C^2(\mathbb{R}^k))\) を満たすために従う．\(C^2(\mathbb{R}^k)\) は，爆発しない拡散過程の生成作用素の核をなす．
\(\ell_k\)-既約な Makrov 過程であって不変確率分布を持つから，再帰的ではある．(S. Meyn and Tweedie, 2009, p. 172)定理8.0.1．これは，\(\psi\)-既約な Markov 連鎖は再帰的であるか推移的であるかの２択だからである(Tweedie, 1994, p. 定理2.3)．
また連続な過程であるから，これは Harris 再帰性を意味する．
不変確率分布を持つ非周期的な Harris 再帰的な Markov 連鎖であるから，全変動収束が従う (S. P. Meyn and Tweedie, 1993, p. 定理6.1)．
Proposition 6.1. Suppose that \(\Phi\) is positive Harris recurrent, and that some skeleton chain is irreducible. If \(C\) is petite, then there exists a constant \(T>0\) and a non-trivial measure \(\mu\) such that \(P^s(x,-)\ge\mu(-)\;(s\ge T,x\in C)\). (S. P. Meyn and Tweedie, 1993)_
\(L_t\) の指数エルゴード性
指数エルゴード性の必要条件
\(V:\mathbb{R}^k\to[1,\infty)\) に対して，\(V\)-一様エルゴード的 であるとは，任意の \(x\in\mathbb{R}^k\) に対して，ある \(R>0,\rho\in(0,1)\) が存在して \[
\|P^t_L(x,-)-\pi\|_V\le V(x)R\rho^t,\qquad t\ge0,
\] が成り立つことをいう．ただし， \[
\|A\|_V:=\sup_{\lvert f\rvert\le V}\int_{\mathbb{R}^k}f(x)A(dx).
\]
(S. Meyn and Tweedie, 2009, pp. 336–337) 第 14 章にいう \(f\)-エルゴード性に当たる概念であることに注意．\(V\)-一様エルゴード性は第 16 章で別の意味に使われている．
\(\nabla\log\pi\) は \(C^1(\mathbb{R}^k)\)-級で，ある \(N,a,b>0\) が存在して \[
\nabla\log\pi(x)\cdot x\le a\lvert x\rvert^2+b,\qquad\lvert x\rvert>N,
\] を満たすとする．このとき，Langevin 拡散 \(\{L_t\}\) について次が成り立つ：
- 任意の \(V\in C^2(\mathbb{R}^k)\) であって，\(V\ge1\) かつ，ある \(b,c>0\) と非空コンパクト集合 \(C\overset{\textrm{cpt}}{\subset}\mathbb{R}^k\) とについて \[
A_LV\le-cV+b1_C
\] を満たすものについて，\(V\)-一様エルゴード的である．
- \(L\) が指数エルゴード的ならば，ある非空コンパクト集合 \(C\overset{\textrm{cpt}}{\subset}\mathbb{R}^k\) について，ある \(\kappa>1,\delta>0\) が存在して，任意のスタート地点 \(x\in \mathbb{R}^k\) について， \[
\sup_{x\in C}\operatorname{E}[\kappa^{\tau_C^\delta}]<\infty,
\] \[
\tau^\delta_C:=\inf\left\{t\ge\delta\mid L_t\in C\right\}.
\]
非空コンパクト集合 \(C\overset{\textrm{cpt}}{\subset}\mathbb{R}^k\) は小集合であるから，これは成り立つ．(S. Meyn and Tweedie, 2009, pp. 336–337) 定理 14.0.1 が Markov 連鎖に関して与えている消息である．
任意の \(\delta>0\) について，\(\delta\)-スケルトンも幾何エルゴード的である．この \(\delta\)-スケルトンに関する \(C\) への到着時刻は \(\tau^\delta_C\) よりも必ず大きくなる．
\(L_t\) は強 Feller であるから，\(C\) は \(\delta\)-スケルトンに関しても小集合である．あとは，Markov 連鎖に関する幾何エルゴード定理より従う．
指数エルゴード性の十分条件
\(-\nabla\log\pi=O(\lvert x\rvert^{2k-1})\;(\lvert x\rvert\to\infty)\) としたとき，後述の定理 2.4 の条件である \(k<1/2\) を満たさない限り，\(k\ge1/2\) まで指数エルゴード性は成り立ち続ける．
ある \(S>0\) が存在して， \(\lvert x\rvert\ge S\) 上で \(\pi\) は有界であるとする．加えて，ある \(d\in(0,1)\) が存在して \[
\liminf_{\lvert x\rvert\to\infty}(1-d)\lvert\nabla\log\pi(x)\rvert^2+\nabla^2\log\pi(x)>0
\] を満たすならば，\(\{L_t\}\) は指数エルゴード的である．
上の定理の 1. のドリフト条件を示せば良い．\(V=\pi^{-d}\;(0<d<1)\) と取れば良い．このとき， \[
2A_LV=V\biggr(\lvert\nabla\log\pi\rvert^2(d^2-d)-d\nabla^2\log\pi\biggl).
\]
指数エルゴード性が成り立たないとき
指数エルゴード性が成り立たない場合は，返ってこない時，棄却率が高すぎる時，急に動く場合の３つが考えられる．
\(\lvert\nabla\log\pi(x)\rvert\to0\) ならば，\(\{L_t\}\) は指数エルゴード的でない．
上の定理の 2. と矛盾させれば良い．幾何エルゴード的だと仮定し，コンパクトな Kendall 集合 \(C\overset{\textrm{cpt}}{\subset}\mathbb{R}^k\) が存在するとする： \[
\sup_{x\in C}\operatorname{E}[\kappa^{\tau^\delta_C}]<\infty.
\] ある \(R>0\) が存在して， \[
\lvert\nabla\log\pi(x)\rvert\le2(\log\kappa)^{1/2}\quad(\lvert x\rvert\ge R).
\] このとき，\(S\ge\sup_{x\in C}\lvert x\rvert\lor R\) を取り，\(\lvert y\rvert\ge2S\) を満たす地点 \(y\in\mathbb{R}^k\) からスタートしたとすると，\(Z_t:=\lvert L_t\rvert\) は，\(a(x)>-(\log\kappa)^{1/2}\;(\lvert x\rvert>S)\) を満たす係数について \[
dZ_t=dW_t+a(L_t)dt
\] を満たす．\(B_t:=W_t-(\log\kappa)^{1/2}t\) とし，それぞれの過程の \(S\) への到着時刻 \(\sigma(Z),\sigma(B)\) を考えると，\(\sigma(Z)\ge\sigma(B)\;\;\text{a.s.}\) が成り立つ．このとき， \[\begin{align*}
\operatorname{P}_y[\tau_C>t]&\ge\operatorname{P}_y[\sigma(Z)>t]\\
&\ge\operatorname{P}_y[\sigma(B)>t]\\
&\ge\Phi\left(\frac{-S+t(\log\kappa)^{1/2}}{\sqrt{t}}\right)-e^{2S(\log\kappa)^{1/2}}\Phi\left(\frac{-S-t(\log\kappa)^{1/2}}{\sqrt{t}}\right).
\end{align*}\] 最後の式は Bachelier-Lévy の公式による．\(\sigma(B)\) の密度を \(f\) とすると， \[
\frac{\log f(t)}{t}\to-\frac{\log\kappa}{2}.
\] これは矛盾を起こすらしい．
指数収束する例
１次元クラス \(\mathcal{E}(\beta,\gamma)\)
多次元指数分布族 \(\mathcal{P}_m\)
6 References
Besag, J. E. (1994).
Comments on “Representations of Knowledge in Complex Systems” by U. Grenander and M. I. Miller.
Journal of the Royal Statistical Society. Series B (Methodological),
56(4), 591–592.
Fearnhead, P., Bierkens, J., Pollock, M., and Roberts, G. O. (2018).
Piecewise deterministic markov processes for continuous-time monte carlo.
Statistical Science,
33(3), 386–412.
Gelfand, A. E., and Smith, A. F. M. (1990).
Sampling-based approaches to calculating marginal densities.
Journal of the American Statistical Association,
85(410), 398–409.
Grenander, U., and Miller, M. I. (1994).
Representations of knowledge in complex systems.
Journal of the Royal Statistical Society. Series B (Methodological),
56(4), 549–603.
Mengersen, K. L., and Robert, C. P. (1996).
Testing for mixtures: A bayesian entropic approach. In
Bayesian statistics 5: Proceedings of the fifth valencia international meetings, pages 255–276.
Meyn, S. P., and Tweedie, R. L. (1993).
Stability of markovian processes II: Continuous-time processes and sampled chains.
Advances in Applied Probability,
25(3), 487–517.
Meyn, S., and Tweedie, R. L. (2009).
Markov chains and stochastic stability. Cambridge University Press.
Neal, R. M. (1994).
An improved acceptance procedure for the hybrid monte carlo algorithm.
Journal of Computational Physics,
111(1), 194–203.
Parisi, G. (1980).
A sequence of approximated solutions to the s-k model for spin glasses.
Journal of Physics A: Mathematical and General,
13(4), L115.
Roberts, Gareth O., and Rosenthal, J. S. (1998).
Optimal scaling of discrete approximations to langevin diffusions.
Journal of the Royal Statistical Society. Series B (Statistical Methodology),
60(1), 255–268.
