MCMC 法と近似精度保証
来嶋 秀治
京都大学 数理解析研究所〒 606-8502 京都市左京区北白川追分町
kijima@kurims.kyoto-u.ac.jp
概要

MCMC 法は積分計算/数え上げのための強力な乱択近似計算法である．MCMC 法において，「マルコフ連鎖を用いたランダムサンプリング法」の側面が強調されるが，モン
テカルロ法にも創意工夫が施されている．本稿では，「設計したサンプリング法からどのようにモンテカルロ法を構築するか」について，解の近似精度保証と計算量の観点から議論する．具体的な対象を例に，自己帰着可能性を利用した再帰的な MCMC 法を紹介する．
Keywords: MCMC（Markov chain Monte Carlo; マルコフ連鎖モンテカルロ）法，数え上げ問題，#P 困難，FPRAS（全多項式時間乱択近似スキーム），自己帰着可能性（self-reducibility）．

1

はじめに

MCMC（Markov chain Monte Carlo; マルコフ連鎖モンテカルロ）法は積分計算/数え上げのための強力な乱択近似計算法である．MCMC 法の「マルコフ連鎖を用いたランダム
サンプリング法」というアイデアは面白く，その側面ばかりが強調されがちだが，もちろん，それだけでは目的の積分計算/数え上げにはたどり着かない．実際，マルコフ連鎖の収束ス
ピードの問題とも相俟って，モンテカルロ法の部分にも多くの創意工夫が施されている．
「設
計したサンプリング法を用いてどのようにモンテカルロ法を構築するか」について，本稿では自己帰着可能性（self-reducibility）に焦点を当てて議論する．解の近似精度保証と計算量という観点において，この自己帰着可能性を利用した再帰的 MCMC 法は非常な成功を収めている．いくつかの具体的な対象を取り上げ，この自己帰着可能性のアイデアを紹介するのが本稿の目的である．

1.1

モンテカルロ法

モンテカルロ法は，平均，比率の計算に高い効果のある乱択計算法（randomized algorithm）
である．モンテカルロ法の説明で最もよく用いられる例は円周率 π （= 半径 1 の円の面積）の計算法であろう．互いに独立な一様ランダムな点 X ∈ [−1, +1]2 を多数生成し，そのうち，原点を中心とする半径 1 の円内にある点の割合に正方形の面積 4 を乗じた値が π の近似値として得られる．この近似値はランダム生成する点の個数が増加すると，大数の法則に従い漸近的に真の値 π となる．また，Chebyshev の不等式などの確率不等式を利用して近似値の精度を確率的に評価することもできる [15, 40, 43]．

同様に（0 で無い面積をもつ）2 次元領域の面積の計算にも応用できる．さらに，3 次元，
4 次元領域の体積計算に対しても同様の手法を考えることもできるだろう．しかし，この素朴なモンテカルロ法の前には，計算量的な意味での「次元の呪い」が立ちはだかっている．
√
例として，半径 n/2 の n 次元ユークリッド球 Bn の体積を考えよう1 ．Bn の体積 V(Bn ) は³ √ ´n³ √ ´nn√ n
π 2 2n2n 2nn¢ < ¥n¦V(Bn ) = ¡ n
= ¥n¦
!
!
Γ 2 +122√を満たす．ただし Γ はガンマ関数である2 ．球 Bn は 1 辺の長さ n の n 次元の超立方体 Cn に内√ n
接するにも関わらず，Cn 中の一様ランダムな点が Bn 内に入る確率 pn は pn = V(Bn )/ n <(b 12 c!)−1 となる．この確率は
n = 10 で pn < 1/5! ' 0.01,
n = 20 で pn < 1/10! ' 2.7 × 10−7 ,n = 100 で pn < 1/50! ' 3.3 × 10−65 ,であり，100 次元ともなると Cn 中のランダムな点が Bn に入ることは，地球上の全原子の中から 1 原子をあてるよりも難しい3 ．すなわち，Bn は 1 辺の長さ 1 の n 次元立方体を完全に含むことから明らかなように 1 以上の体積を持つにもかかわらず，素朴なモンテカルロ法を用いると，現実的な時間では高次元球の体積は “0” になってしまうのである．

1.2

数え上げ問題/積分計算の計算量

計算量的な観点において，高次元の体積計算は “難しい” ことが知られている．1979 年にValiant は数え上げ問題の計算量クラスとして#P 完全（#P 困難）4 の概念を提唱した [51]．
#P は，クラス NP に属する問題の許容解の個数を計算する問題のクラスである．数え上げ
問題と積分計算は，Pick の定理，その一般化である Ehrhart 多項式などを介して密接な関係を持つ．自然，多くの積分計算は容易に#P 困難となり得る．高次元凸多面体の体積計算もまた#P 困難である [7]．
NP 完全などの “難しい” 問題の数え上げだけでなく，クラス P の多くの問題に対する数え上げ問題も#P 完全となることが知られている．Valiant は [51] において，正方 0-1 行列のパーマネントの計算が#P 完全であることを示している．この問題は 2 部グラフの完全マッチングの数え上げ問題と等価である．このほか 0-1 ナップサック解の数え上げ，輸送問題の整数解の数え上げなども#P 完全である．逆に多項式時間で数え上げのできる対象としては，グラフの全張木の数え上げ，平面グラフの完全マッチングの数え上げがある [17]．
1球の例は，本稿で「モンテカルロ法を用いて高次元 “球” の体積計算をしよう」という意味ではない．あくまで，球という “ふっくら” した物体（= 超立方体の全ての面に内接する凸体）に対してさえ，モンテカルロ
法は非効率的になり得るという例である．
√
2正の実数 x に対して，Γ(x + 1) = xΓ(x) が成り立ち，Γ(1) = 1，Γ( 12 ) = π である．特に正の整数 n に対して Γ(n + 1) = n! となる．
3
ちなみに全宇宙の素粒子の個数は 1080 程度と言われている．
4
本来，#P は「ナンバーピー」と読むべきであるが，近年では「シャープピー」と読まれることの方が多い．#（ナンバー）は横棒が水平，]（シャープ）は縦棒が垂直であり，両者は異なる記号である．しかし，最近では ]P という記述をする文献も見かけるようになった．

1.3

乱択計算の近似精度保証

数え上げ困難な問題に対する効率的な近似計算法として，乱択化は非常に有望なアプロー
チである．乱択近似計算においては，出力される解が確率的であることから，近似精度保証にも，決定的な近似計算法には無い「確率的な評価」が必要となる．1983 年に Karp & Luby はRAS(randomized approximation scheme; 乱択近似スキーム) の概念を提案している [25, 26]．
任意の入力 ε, δ > 0 に対して，乱択計算法が，厳密解 A に対する近似解 Z を出力して，µ¶|Z − A|Pr≤ε ≥1−δAを満たす時，その乱択計算法は RAS と呼ばれる．すなわち，RAS で得られる近似解は高い確率で相対誤差が保証されている．特に，poly(入力サイズ, ε−1 , ln δ −1 ) 時間の RAS を FPRAS（fully polynomial-time RAS ; 全多項式時間乱択近似スキーム）と呼ぶ．
多くの乱択計算法の中で，#P 困難な問題に対する FPRAS の設計は（理論的な多項式性の意味で）最も成功しているもののひとつと言える5 ．そして，多くの数え上げ困難な問題に対して FPRAS を与えているのが MCMC 法である．

1.4

MCMC 法MCMC 法は，サンプリングにマルコフ連鎖を利用するモンテカルロ法であり，数値積分，シミュレーションなどに用いられる．大規模な空間を持つ対象に対して効果的な計算法であり，特にランダムサンプリング自体が困難な問題に対して強力に効果を発揮する．
MCMC 法におけるサンプリングの基本的なアイデアは，「（エルゴード的な）マルコフ連
鎖を繰り返し推移させ，漸近的に定常分布に従うサンプルを得る」という非常に単純なものである．したがって，適切なマルコフ連鎖を設計することで，所望の分布に対するサンプリ
ング法を構築することができ，各種の確率的アルゴリズムに組み込むことができる．この枠組みの簡明さゆえに，MCMC 法は，統計物理学，経済学，統計学，バイオインフォマティ
クス，オペレーションズリサーチなどの様々な分野に頻繁に現れる．
MCMC 法により，状態空間の大きい対象や複雑な分布に対して効率的なサンプリング法が設計できるようになり，さまざまな#P 困難な問題に対する効率的な乱択近似計算法の可能性が広がった．実際に，1991 年の Dyer, Frieze & Kannan による高次元凸体の体積計算に対する FPRAS（[8]），2004 年の Jerrum, Sinclair & Vigoda によるパーマネントの計算に対する FPRAS（[21]）といった重要な成果が MCMC 法によって得られている．

1.5

本稿のねらいと構成

本稿では FPRAS を目的とした MCMC 法の設計について，高次元凸体の体積計算 [7, 8,
24, 17]，パーマネントの計算 [4, 18, 21, 20] を例に紹介する．1.1 節の例のように，素朴なモ
ンテカルロ法では難しいであろう FPRAS の構築について，再帰的にモンテカルロ法を利用
5

近年，多くの乱択アルゴリズムに対して，脱乱択化（derandomization）の手法によって，同等の性能を持
つ決定的アルゴリズムが得られるようになってきた．代表的な例として，素数判定が挙げられる．多項式時間の乱択近似計算法があって，多項式時間の決定的近似計算法の知られていない問題は，それほど多く無い．

する巧妙なアルゴリズムを俯瞰する．この再帰的なアルゴリズムを支える自己帰着可能性に着目し，精度保証のための要点について考察する．モンテカルロ法，マルコフ連鎖を用いた
サンプリング法，自己帰着可能性，近似精度保証と計算量は互いに緊密に関連しており，それぞれを完全に切り離して議論することはできない．中でも，マルコフ連鎖の収束スピードの算定はこの分野における中心的な課題ではあるが，本稿ではこの話題は必要最小限にとどめ，冒頭にも述べたとおり，「数え上げ/積分計算のためのモンテカルロ法の設計」を中心に議論する．
以下，第 2 章ではマルコフ連鎖を用いたランダムサンプリング法について簡単に述べる．
第 3 章では高次元凸体の体積計算のための MCMC 法を紹介する．第 4 章では第 3 章の方法から抽出される自己帰着可能性の性質について考察し，いくつかの離散的対象に対する数え上げ計算法の例を紹介する．第 5 章ではパーマネントの計算法を紹介する．第 6 章でまとめと関連する話題について述べる．
本稿では実数全体の集合を R で表し，また整数（非負整数，正整数）全体の集合をそれぞれ Z （Z+ , Z++ ）で表す．

マルコフ連鎖を用いたランダムサンプリング法

2

本稿では離散時間のマルコフ連鎖を扱う．さらに，この章では有限の状態空間 Ω と推移
確率行列 P をもつマルコフ連鎖 M を対象とし，定常分布への収束に関する定義を行う6 ．
¡ ¢マルコフ連鎖 M が既約（irreducible）であるとは，任意の状態対 {x, y} ∈ Ω2 に対して，∃t > 0, Pr(X t = y | X 0 = x) > 0 が成り立つことを言う．ただし，各 t ≥ 0 に対する確率変数 X t はマルコフ連鎖 M に従う確率的な時系列とする．マルコフ連鎖 M が非周期的（aperiodic）であるとは ∀x ∈ Ω, gcd{t ∈ Z++ | Pr(X t = x | X 0 = x) > 0} = 1 が成り立つことを言う．ただし gcd は最大公約数を表す．既約で非周期的な有限マルコフ連鎖をエルゴー
ド的（ergodic）と呼ぶ．エルゴード的なマルコフ連鎖は唯一の定常分布を持ち，極限分布は定常分布に一致する．

2.1

定常分布の設計

次の定理は，マルコフ連鎖の定常分布を設計する上で重要な定理である．
定理 2.1 （[14] 参照．
）状態空間 Ω と推移確率行列 P を持つエルゴード的なマルコフ連鎖 Mの定常分布を π とする．関数 f : Ω → R++ が次の詳細均衡式（detailed balance equation）
µ ¶
Ω
(1)
f (x)P (x, y) = f (y)P (y, x),∀{x, y} ∈2を満たす時，関数 f に対して，正の定数 C が存在して，任意の x ∈ Ω に対して π(x) = C ·f (x)が成り立つ．ただし，P (x, y)，P (y, x) は，推移確率行列 P における x から y へ，y から x へのそれぞれの推移確率を表す．
¥
6

次の章で，連続空間のマルコフ連鎖が出てくるが，連続空間の場合でも本質的に同様の議論ができる．

Metropolis-Hastings MCMC 法のためのマルコフ連鎖の設計の多くは，この詳細均衡式を利用したものである．例として，空間 Ω 上で関数 f : Ω → R+ に比例する分布を設計する
ための Metropolis-Hastings 法を記そう．いま，状態空間 Ω 中に近傍が定義され，各状態の近傍サイズ（近傍の個数）は高々k とする．マルコフ連鎖の推移は次のように定義される．現在の状態 X の近傍 Y を（定数）確率 1/k で選ぶ．確率 f (Y )/(f (X) + f (Y )) で状態 Y に推移し，それ以外の場合は X にとどまる．このマルコフ連鎖の定常分布は f に比例する．従って，効率的に近傍を見つけることができ，状態 x ∈ Ω に対して f (x) が容易に計算できれば，所望の分布が漸近的に得られる．

2.2

マルコフ連鎖の収束

ここでは，マルコフ連鎖の収束具合を見積もるための道具を紹介する．まず，状態空間Ω 上の分布の間の距離を以下のように定義する．
定義 2.2 （[14] 参照．
）同一の有限状態空間 Ω 上の 2 つの確率分布 ν1 と ν2 が与えられた時，
ν1 と ν2 の間の全変動距離（total variation distance）を
(
)
X
1X
def.
dTV (ν1 , ν2 ) = max|ν1 (x) − ν2 (x)|(ν1 (x) − ν2 (x)) ≡
A⊆Ω
2
x∈A

x∈Ω

と定義する．
全変動距離の定義を用いて収束時間の尺度を以下のように定義する．
定義 2.3 （[14] 参照．
）任意の正数 ε < 1 に対して，状態空間 Ω を持つエルゴードマルコフ
連鎖 M の混交時間（mixing time）はdef.

τ (ε) = max{min{t | ∀s ≥ t, dTV (π, Pxs ) ≤ ε}}x∈Ωと定義される．ただし，π はマルコフ連鎖 M の定常分布とし，Pxs は初期状態 x ∈ Ω として
時刻 0 から時刻 s ≥ 0 まで推移させた時のマルコフ連鎖 M の確率分布とする．
すなわち，マルコフ連鎖を τ (ε) 回推移させると，定常分布との総変動距離が ε 以下の分布が得られる．特に，τ (ε) ≤ poly(入力サイズ, ε−1 ) の時，マルコフ連鎖は多項式時間で収束するという．

3

高次元凸体の体積計算
以下，Bn (x, r) は x ∈ Rn を中心とする半径 r ∈ R++ の n 次元ユークリッド球
def.

Bn (x, r) = {y ∈ Rn | ky − xk2 < r}を表すとする．次の条件を満たす n 次元の凸体 K ⊂ Rn を考えよう．

√
条件 1 凸体 K は n 次元球 Bn (0, n(n + 1)) に完全に含まれ，n 次元球 Bn (0, 1) を完全に含√む．すなわち，Bn (0, 1) ⊆ K ⊆ Bn (0, n(n + 1)) を満たす7 ．
また，凸体 K には帰属判定オラクル（membership oracle）を仮定する．すなわち，点 x ∈ Rnが与えられた時，x ∈ K の判定は単位時間で行えるとする．
以上の条件を満たす凸体 K の体積 V(K) の計算を考えよう．尚，有界で全次元的な凸体0
K に対して，適切なアフィン変換を施すことで条件 1 は得られる．また，アフィン変換に基づく体積変換は容易に計算できる8 ．
体積計算の基本的なアイデア 凸体 K に対して，r1 ∈ R++ を 1 < r1 <def.
K1 = K ∩ Bn (0, r1 ) とする．このとき，V(K) =√

n(n + 1) として，

V(K)
· V(K1 )V(K1 )(2)が成り立つ．したがって，V(K)/V(K1 ) と V(K1 ) の値がわかれば V(K) が計算できる．いま，
仮に K 中の点を一様サンプリングできたとしよう．すると，V(K)/V(K1 ) をモンテカルロ
法で見積もることができる．すなわち，K 中の一様ランダムな点を多数生成し，この時 K1に入る点の個数の割合は V(K1 )/V(K) の近似値である．したがって，V(K) の計算は V(K1 )の計算に帰着されることになる．
このとき K1 は明らかに条件 1 を満たしている．つまり，再帰的な構造を考えることが√できそうである．いま，適当な長さの実数列 r0 , . . . , rq は r0 = n(n + 1)，rq = 1 とし，各def.

i ∈ {1, . . . , q} に対して ri−1 > ri を満たすとする．また Ki = K ∩ Bn (0, ri ) (i ∈ {0, . . . , q})とする．このとき，式 (2) を再帰的に考えると，
V(K0 )· V(K1 )V(K1 )
V(K0 ) V(K1 )=·· V(K2 )V(K1 ) V(K2 )= ···V(K) ==

V(K0 ) V(K1 )V(Ki−1 )V(Kq−1 )········ V(Kq )V(K1 ) V(K2 )V(Ki )
V(Kq )(3)が得られる．もし，各 Ki−1 (i ∈ {1, . . . , q}) 内の一様ランダムサンプリングが行えるとすれ
ば，V(Ki )/V(Ki−1 ) (i ∈ {1, . . . , q}) をモンテカルロ法で見積もることができる．また，条件 1と rq = 1 を考えると，Kq = K ∩ Bn (0, rq ) = Bn (0, 1) となりnπ2¢V(Kq ) = V(Bn (0, 1)) = ¡ n
Γ 2 +1と計算できる．以上の議論から，式 (3) を用いて凸体 K の体積を近似的に計算できる．
7

したがって，もちろん K は有界で全次元的である．また，K が閉 (closed) もしくは開 (open) であることは必ずしも必要ではないが，後のマルコフ連鎖の議論の簡便のため，ここでは開とする．
8
アフィン変換についての議論は割愛する．[17] を参照されたい．

マルコフ連鎖 しかし，K 中の一様ランダム生成はどのようにして得られるのであろうか．
ここでは ball-walk と呼ばれる（連続空間）マルコフ連鎖を紹介しよう．以下，n 次元ユーク
リッド球中の一様乱数オラクルを仮定する9 ．
この ball walk は一種の Metropolis-Hastings 法で，現在の点 X ∈ K から次の時刻の状態√
0
X への推移は次のように定義される10 ．まず，次の状態の候補 Y を Bn (X, n) の中から一様ランダムに選ぶ．次に，Y ∈ K なら X 0 = Y とし，Y 6∈ K なら X 0 = X とする．このとき，明らかに X 0 ∈ K である．このマルコフ連鎖の定常分布が K 中の一様分布であることは，離散のマルコフ連鎖の詳細均衡式 (1) を注意深く連続空間に拡張することで示される．また，
このマルコフ連鎖の 1 回の推移にかかる計算量は，球中の乱数生成の時間と K の帰属判定の時間である．
精度保証のために 以上，基本的なアイデアを紹介したが，実際には様々なパラメータを設定する必要がある．まずは，Ki を定める半径の列 r0 , . . . , rq について考えよう．半径の比
ri−1 /ri をあまり大きくしすぎると，モンテカルロ法で V(Ki )/V(Ki−1 ) を見積もる時に高次元の呪いにはまり，1.1 節の例のように多項式時間では有効な値が得られない可能性がある．
一方，ri−1 /ri をあまり小さくすると，反復回数が多くなってしまう．
√
現実的な方法としては，ri = n 2 ri−1 と設定すると，この問題を解決できる．このとき，V(Ki )
1≥
V(Ki−1 )2が成り立ち，モンテカルロ法の近似精度の問題は回避される．また，再帰の回数 q もlm » ln √n(n + 1) ¼ » ln √n(n + 1) ¼√
q ≤ log2 n1 n(n + 1) == O(n ln n)
= n·
1
ln 2
ln 2 nを満たし，反復回数は n の多項式で押さえられる．
他には，マルコフ連鎖の推移回数とモンテカルロ法のサンプリング回数を議論する必要がある．マルコフ連鎖の推移回数については，総変動距離が十分小さな ε0 になる回数を決めればよい．条件 1 の下で，ball-walk が多項式時間収束することが示されている [23]．モンテ
カルロ法のサンプリング回数については，マルコフ連鎖の推移回数から得られる総変動距離ε0 を考慮したベルヌーイ試行 —すなわち Ki−1 の一様ランダムな点が Ki に入るか否か— を考えると，Chernoﬀ 限界などの確率不等式を用いて十分なサンプリング回数を算定することができる．これらの算定を注意深く組み合わせることで，FPRAS が得られる．
自己帰着可能性 本章で紹介したアルゴリズムの基本的アイデアにおいて，重要な役割を果たしているのが，自己帰着可能性である．すなわち，Ki の体積計算という問題は，条件 1 を含めて再帰的構造を持つ．条件 1 を満たすことにより，

1. ball-walk が高速に収束する，
2. メンバーシップオラクルをもつ11 ，という性質が保証される．これらの性質が再帰計算の効率性に大きく寄与している．
9

球中のランダムサンプリングについての議論も割愛する．[11] を参照されたい．
ここでは簡便のため K 内の一様ランダム生成法について述べるが，各 Ki についても同じ方法で得られる．
11x ∈ Ki の判定は，x ∈ K の判定時間と，kx − 0k2 < ri の判定時間で行える．後者の時間についても，xがランダムな点であることを考えると，（期待値として）十分速い時間で計算することができる．
10

4

自己帰着可能性と離散構造

前章では，対象が連続的な高次元凸体の体積計算に対して，自己帰着可能性を利用した再帰計算を取り扱った．これを参考に，改めて自己帰着可能性について考察してみると，集合 Ω の部分集合 ∆ が，帰着先の集合 Ω0 と（一対一）対応する時，|Ω| |Ω1 ||Ωi−1 ||Ωq−1 |········ |Ωq ||Ω1 | |Ω2 |
|Ωi |
|Ωq |
|Ωi |
|Ωq−1 ||Ω0 | |Ω1 |
·
···
···
· |Ωq |=|∆0 | |∆1 ||∆i ||∆q−1 ||Ω| =という再帰式が得られる．これを基に，各 Ωi のランダムサンプリングが行え，かつ最終的な |Ωq | が容易に計算可能であれば，元問題 |Ω| に対する近似計算が行えるという構造である．
再帰的に現れる各 Ωi が共通の構造をもっているならば，共通の手法に基づくランダムサ
ンプリング法が得られるであろう．このとき，次の点に気をつけなければならない．

1. 帰着先もまた “性質の良い” 問題であること．
2. （状態空間のサイズに関する）スケーリングが良く効くこと．
“性質の良い” とは，例えば先の例ではモンテカルロ法が多項式時間で良い近似精度の解を返す必要がある．たとえ共通の構造を持っていても，再帰で起こり得る全ての場合において取扱い易い構造でなければ，この再帰計算法はうまく働かない．一方，スケーリングに関しては，離散構造が往々にして自明な良い再帰構造をもつ．
本章では，“良い” 自己帰着可能性を持つ離散的対象の例を 2 つあげよう．

4.1

0-1 ナップサック問題の解の個数
def.

ベクトル a ∈ Rn++ と実数 b ∈ R++ が与えられた時，集合 Ωn = {x ∈ {0, 1}n | a · x ≤ b}の要素数を計算する問題がナップサック問題の解の数え上げ問題である．この問題は#P 完
全であることが知られている．
def.
− def.
添え字 n に注目して，∆+n = {x ∈ Ω | xn = 1}，∆n = {x ∈ Ω | xn = 0} とする．もし00x ∈ ∆+
n ならば，n 番目のアイテムを取り除いたナップサック解 x = (x1 , . . . , xn−1 , 0) は x ∈
−
+
−
+
−
∆−
n である．つまり |∆n | ≥ |∆n | が成り立つ．∆n ∪∆n = Ωn であるから，|∆n |/|Ωn | ≥ 1/2 が成
def.
n−1
り立つ．∆−| a0 ·x ≤ b}n は n−1 アイテムのナップサック問題の解集合 Ωn−1 = {x ∈ {0, 1}n−1と対応する．ただし a0 ∈ R++
は a0 = (a1 , . . . , an−1 ) とする．
したがって，Ωk (k ∈ {2, . . . , n}) のランダム生成ができるならば，|Ωn | |Ωn−1 |
|Ωi |
|Ω2 |
·
···
···
· |Ω1 ||Ωn−1 | |Ωn−2 ||Ωi−1 |
|Ω1 |
|Ωi |
|Ω2 |
|Ωn | |Ωn−1 |· − · · · − · · · − · |Ω1 |=−|∆n | |∆n−1 |
|∆i |
|∆2 |

|Ωn | =という再帰式に基づいて，再帰的なモンテカルロ法を構成できる．|Ω1 | は高々2 通りの解の実行可能性を確認すれば良いので，容易に計算できる．
0-1 ナップサック問題の解に対する素朴なマルコフ連鎖が設計でき，素朴なマルコフ連鎖
は多項式時間で収束することが示されている [42]．

4.2

2 元分割表の個数（輸送問題の整数解の個数）n2 つのベクトル r ∈ Zm
+ と s ∈ Z+ が与えられ，krk1 = ksk1 を満たすとする．2 元 m × n
分割表の数え上げ問題は，輸送多面体中の整数格子点集合¯ P)(¯ n x = r (i ∈ {1, . . . , m}),def.
¯
i
j=1 i,jΩ(r, s) = X = (xi,j ) ∈ Zm×n
¯ Pm+
¯ i=1 xi,j = sj (j ∈ {1, . . . , n})の要素数の計算であり，#P 完全であることが知られている [10]．
いま，添え字 j ∗ ∈ {1, . . . , n} は sj > 0 を満たすものとする．集合 ∆ij ∗ ⊂ Ω(r, s) を
def.

∆ij ∗ = {X ∈ Ω(r, s) | xi,n ≥ sn /m}(i ∈ {1, . . . , m})Smと定義すると， i=1 ∆ij ∗ = Ω(r, s) である．従って，少なくともひとつ i∗ ∈ {1, . . . , m} が存∗在して，|∆ij ∗ |/|Ω(r, s)| ≥ 1/m が成り立つ．このとき，r 0 ，s0 を
r0

def.

=

(r1 , . . . , ri∗ −1 , ri∗ − dsj ∗ /me, ri∗ +1 . . . , rm ),

s0

def.

(s1 , . . . , sj ∗ −1 , sj ∗ − dsj ∗ /me, sj ∗ +1 . . . , sn )=∗とすると，∆ij ∗ と Ω(r 0 , s0 ) には一対一の対応がある．したがって，再帰的なモンテカルロ法が構築できる [28]．
この問題は，モンテカルロ法を行うのに，集合分割的ではなく，Ω(r, s) を被覆する集合族を考えているのが特徴的である12 ．またナップサックの時と異なり，モンテカルロ法で得∗られる比は |∆ij ∗ |/|Ω(r, s)| ≥ 1/m しか保証されない．しかし，この場合でも多項式回のサンプリングで十分良い近似値が得られる．現在，2 元分割表に対しては m（または n）が定数の場合だけ多項式時間のランダムサンプリング法が知られている [9, 29, 6]．一般の m，n に関する多項式時間のランダム生成法は未解決である．

パーマネントの計算

5

本章では，素朴には効率的な自己帰着可能性が得られない離散構造の例として，パーマネントの計算法の例を紹介しよう．
パーマネントの計算は古くから知られる問題である13 ．n × n の 0-1 行列のパーマネントの計算は Valiant が示した最初の#P 完全問題である [51]．この問題は頂点数 n + n の 2 部グ
ラフ G = (U, V ; E) の完全マッチングの数え上げ問題と等価である．

5.1

問題

いま，n × n 行列 A = (ai,j ) に対して，パーマネント per(A) は

Qn
def. P
i=1 ai,σ(i)σ∈Snper(A) =0

すなわち，一般に ∆ij ∗ ∩ ∆ij ∗ 6= ∅．
13
平面グラフの完全マッチングの数え上げに対しては，Pfaﬃan 向き付けを利用した多項式時間の決定的ア
ルゴリズムが存在する．
12

と定義される．ただし Sn は n 次の対称群を表す．行列 A の行列式

Qn
def. P
i=1 sgn(σ)ai,σ(i)σ∈Sndet(A) =と比較すると，足し合わせの各項の値は符号を除いて等しい．特に行列 A の各成分が 0 また
は 1 のとき，枝集合 E = {(i, j) ∈ U × V | aij = 1} を持つ 2 部グラフ G = (U, V ; E) の完全マッチングの個数と per(A) は一致する．以下，2 部グラフ G = (U, V ; E) の完全マッチングの集合を Ω で表し，|Ω| の乱択近似計算法について議論する．

5.2

素朴な自己帰着枝 e ∈ E を含む完全マッチングの集合を ∆e ⊂ Ω で表す．枝 e∗ = (u, v) は全ての e に対S
して，|∆e∗ | ≥ |∆e | を満たすとする．このとき， e∈E ∆e = Ω より，|∆e∗ |/|Ω| ≥ 1/m を満たす．グラフ G0 を G から頂点 u, v を取り除いたグラフとすると，G0 の完全マッチング全体の集合 Ω0 は ∆e∗ と一対一対応をもつ．従って，2 部グラフの完全マッチングのランダムサンプ
リングが効率的に行えると |Ω| の計算が行える．

5.3

マルコフ連鎖の設計

そこで，2 部グラフの完全マッチングの一様ランダム生成法を考えよう．マルコフ連鎖を用いたランダム生成法を設計するのだが，完全マッチング M ∈ Ω の “近傍” の数え上げが難しいので，詳細均衡式を利用して定常分布を設計することは容易ではない．Broder は n − 1
マッチングを状態空間に加えたマルコフ連鎖を提案した [4]．このマルコフ連鎖を紹介する．
頂点対 (u, v) ∈ U × V に対して，Ω[u, v] を u と v にはマッチングの枝が接続していない
n − 1 マッチング全体の集合とする．このとき，必ずしも (u, v) ∈ E では無いことに注意がS
def.
必要である．表記の簡便のため，Ξ = Ω ∪ (u,v)∈U ×V Ω[u, v] とする．状態空間 Ξ をもつマ
ルコフ連鎖の現在の状態 X ∈ Ξ から次の状態 X 0 への推移は次のように定義される．まず，
e = (u, v) ∈ E を一様ランダムに選ぶ．

1. もし X ∈ Ω かつ e ∈ X なら X 0 = X − e とする．
2. もし X ∈ Ω[u, v] なら，X 0 = X + e とする．
3. もし X ∈ Ω[u, v 0 ] (v 0 6= v) なら，枝 e0 = (u0 , v) ∈ X に対して，X 0 = X + e − e0 とする．
4. もし X ∈ Ω[u0 , v] (u0 6= u) なら，枝 e0 = (u, v 0 ) ∈ X に対して，X 0 = X + e − e0 とする．
5. それ以外の場合 X 0 = X とする．
明らかに X 0 ∈ Ξ である．このマルコフ連鎖の定常分布は Ξ 上の一様分布となる．したがって，マルコフ連鎖を用いて，Ξ 上の一様ランダムな状態 X を生成し，X ∈ Ω なら出力するという方法で，完全マッチングの一様ランダムサンプリングが行える．

Jerrum & Sinclair によって，Broder のマルコフ連鎖は多項式時間で収束することが示されている [19]．しかし，実は上述の方法では完全マッチングのランダムサンプリングはうまn
くいかない．というのも，|Ω|/|Ξ| = O(2− 3 ) となる例が存在し，高次元の呪いに捕らわれてしまうからである．そこで，Broder のマルコフ連鎖を改変して，効率的なサンプリングを目指そうというのが今後の方針である．

5.4

パーマネントの計算に対する FPRAS の設計

基本的なアイデアとマッチング重み関数 2 部グラフ G = (U, V ; E) に対して，マッチング
重み関数 w : U × V → R+ が定義されているとする．このとき，Ξ 上の分布 µw を，1(M ∈ Ω)
Cdef.
µw (M ) =w(u, v)(M ∈ Ω[u, v])
Cdef. Pと定義する．ただし，C は正規化定数で C =M ∈Ξ µw (M ) である．マッチング重み関数 w を適切に設定して，定常分布が µw となるように Broder のマルコフ連鎖を Metropolis-Hastings
法で改変しようというのがアイデアである．
マッチング重み関数 w(u, v) が 1 に近い時には，Broder のマルコフ連鎖に近い挙動をし，多項式時間の収束が見込まれる反面，高次元の呪いの問題が出現する．一方，w(u, v) = 0 とすると，完全マッチングの出現確率は大きくなるが，マルコフ連鎖は既約でなくなってしま
い，収束が遅くなってしまう．
def.
マッチング重み関数 w = w∗ が w∗ (u, v) = |Ω|/|Ω[u, v]| と定義されているとすると，任P意の (u, v) ∈ U × V に対して，w∗ (u, v) · |Ω[u, v]| = |Ω| が成り立つことから， M ∈Ω µ(M ) =1/(n2 + 1) が成り立つ．すなわち，Ξ 上の µ に従うランダムな点が Ω に入る確率は 1/(n2 + 1)となり，素朴な手法で問題となった高次元の呪いを回避することができる．また，このときマルコフ連鎖も多項式時間で収束することを Jerrum, Sinclair & Vigoda は示している [21]．
したがって，w∗ を計算できれば 2 部グラフの完全マッチングのランダムサンプリングが効率的に行える．しかし，マッチング重み関数値 w∗ (u, v) の計算は容易ではない．そこで，w∗ (u, v) を近似計算しようと言うのが，Jerrum, Sinclair & Vigoda のアイデアである [21]．
枝活性度の導入と自己帰着性 完全 2 部グラフ Kn,n に対して，枝活性度関数 λ : U × V → R+が定義されているとする．完全 2 部グラフ Kn,n の完全マッチング全体の集合を Ω0 で表し，
頂点対 (u, v) ∈ U × V がマッチされていない n − 1 マッチングの集合を Ω0 [u, v] で表す．まS
def.
た，Ξ0 = Ω0 ∪ (u,v)∈U ×V Ω0 [u, v] とする．枝活性度関数の定義を拡張して，M ∈ Ξ0 に対しdef. Qても λ(M ) = e∈M λ(e) と定義する．
2 部グラフ G = (U, V ; E) の枝集合 E に対して，たとえば，λ = λ∗ を λ∗ (e) = 1 (e ∈ E)，P
λ∗ (e) = 1/n! (e 6∈ E) と定義すると， M ∈Ω0 λ(M ) は 2 部グラフ G の完全マッチングの個
数 |Ω| の近似解となる．また，λ = λ0 を λ0 ((u, v)) = 1 ((u, v) ∈ U × V ) と定義すると，
P0∗
M ∈Ω0 λ0 (M ) = |Ω | = n! と，容易に計算できる．そこで，λ を λ0 から λ へ少しづつ変化させていく．
枝活性度関数 λ をもつ完全 2 部グラフ Kn,n に対して，マッチング重み関数 w がP
0 λ(M )w(u, v) ' P M ∈Ω((u, v) ∈ U × V )(4)
M ∈Ω0 [u,v] λ(M )def.

)w(M )
を定常分布にもつ Broder 型のマルコを満たしている時，Ξ0 上の分布 µλ,w (M ) = λ(MC(λ,w)
フ連鎖が多項式時間で収束することを Jerrum, Sinclair & Vigoda は示している [21]．ただし，
M ∈ Ω[u, v] に対して便宜的に w(M ) = w(u, v) と表し，M ∈ Ω に対しては w(M ) = 1 として

def. P02いる．また，C(λ, w) =M ∈Ξ0 λ(M )w(M ) は正規化定数である．この時，µ(Ω ) は 1/O(n )となる．したがって，λ∗ に対して，(4) を満たすようなマッチング重み関数 w∗ が計算でき，
C(λ∗ , w∗ ) を計算できれば，|Ω| = C(λ∗ , w∗ ) · µλ∗ ,w∗ (Ω0 ) という形で，2 部グラフの完全マッ
チングの計算が可能である．
1
詳細は省略するが，枝活性度関数の列 λ0 , λ1 , · · · , λq を λi+1 = e− 2 λi を満たすように設定すると，λi と (4) を満たすようなマッチング重み関数 wi から，wi+1 および C(λi , wi )/C(λi+1 , wi+1 )をモンテカルロ法で推定することができ，

C(λ∗ , w∗ ) = C(λq , wq ) =C(λq , wq )C(λq−1 , wq−1 )C(λ1 , w1 )····· C(λ0 , w0 )C(λq−1 , wq−1 ) C(λq−2 , wq−2 )C(λ0 , w0 )の形で，所望の結果が得られる．

6

まとめ

本稿では，自己帰着可能性をキーワードに，数え上げ困難な問題に対する再帰的なモン
テカルロ法を利用した近似精度保証つきアルゴリズムの設計法を紹介した．
マルコフ連鎖の収束具合の算定 理論的な精度保証と計算量の評価を行うためには，マルコフ連鎖の収束具合，モンテカルロ法のサンプリング回数，スケーリングの効率性など，相互に緊密に関連する話題に注意して，確率不等式をはじめとする道具を駆使した細心の議論が必要となる．中でもモンテカルロ法のための効率的なランダム生成法を支える「マルコフ連鎖の収束スピードの算定」は，この研究における中心的な話題と言える．マルコフ連鎖の収束スピードを算定するために，コンダクタンス法 [47, 48]，カップリング法 [1, 5]，log-Sobolev
不等式 [41] など様々な算定手法が開発されている [12]．
一方，マルコフ連鎖のシミュレーション方法を工夫することで，定常分布に厳密に実現するという，完璧サンプリング法 (perfect sampler) が Propp & Wilson によって提案されている [45, 14, 34]．この奇抜な手法に対して “うまく” マルコフ連鎖が設計することができれば，非常に実用的で理論的保証を持つサンプリング法が得られる [45, 29, 30, 31, 38]．
関連する話題と様々な分野の関連 自己帰着可能性を介して，（近似）数え上げと（近似）ランダム生成法は密接に関連することが知られている [22]．自己帰着可能性は，アルゴリズム
設計では古くから良く利用される基本的な性質である．今後，計算幾何，最適化，列挙算法などで培われた技法を利用し，または還元することで，重要な未解決問題の解決の糸口が得られることが期待される．
効率的な MCMC 法の設計についての議論は，アルゴリズムの分野だけではなく，統計物
理学，統計学をはじめとする様々な分野で，理論，実用の両面から研究されている．統計学の分野では，本稿で紹介した Jerrum. Sinclair & Vigoda のパーマネントの計算法のアイデ
アとも近いブリッジサンプリング，パスサンプリングが活発に議論されている [39, 13, 44]．
また，統計物理の分野では，アニーリング法とも関係のあるレプリカ交換法，拡張アンサン
ブル法といった手法が研究されている [50]．

MCMC 法を使おう！ MCMC 法を用いた多項式時間乱択近似計算法については，本稿であげた例の他にも，Ising モデルの分配関数の計算 [19, 45, 16]，s-t 経路ネットワーク信頼性の計算 [27]，暗号の安全性の解析手法 [35]，待ち行列ネットワークの指標計算 [32, 33]，などが提案されている．また，高次元凸体のランダムサンプリング法を利用して凸計画問題を解く乱択アルゴリズムが提案されている [2]．
MCMC 法は実用的手法としても優れており，生態情報学，人工知能，機械学習，データ
マイニングなど様々な応用諸分野で用いられている．他のヒューリスティクス同様，最悪ケースを考えるとマルコフ連鎖の収束に時間がかかり，近似精度の理論的保証が困難な問題でも，実用的なインスタンスに対しては，収束が速く，精度の高い近似解が得られる場合もある．
実験的にマルコフ連鎖の収束を算定する手法なども提案されており，ヒューリスティクスとしてパラメータをチューニングすれば，実用的に十分使える道具である．また，MCMC 法の利点として並列化が容易という点もあげられる．計算が困難な問題に対する解決策の一つとして，ぜひ MCMC 法を試していただきたい．もし性能が良ければ，それから理論的な精度保証に挑戦しても遅くはない．

参考文献[1] D. Aldous, Random walks on ﬁnite groups and rapidly mixing Markov chains, in Séminariede Probabilitiés XVII 1981/1982, vol. 986 of Springer-Verlag Lecture Notes in Mathematics,
D. Dold and B. Eckmann, (ed.), Springer, 1983, 243–297.
[2] D. Bertsimas and S. Vempala, Solving convex programs by random walks, Journal of theACM, 51 (2004), 540–556.
[3] I. Bezakova, D. Stefankovic, E. Vigoda, and V. V. Vazirani, Accelerating simulated annealing
algorithm for the permanent and combinatorial counting problems, Proceedings of the 17th
ACM-SIAM Symposium on Discrete Algorithms (SODA 2006), 900–907.
[4] A. Broder, How hard is it to marry at random? (On the approximation of the permanent),
Proceedings of the 18th Annual ACM Symposium on Theory of Computing (STOC 1986),
50–58.
[5] R. Bubley and M. Dyer, Path coupling: a technique for proving rapid mixing in Markov
chains, Proceedings of the 38th Annual Symposium on Foundations of Computer Science(FOCS 1997), 223–231.
[6] M. Cryan, M. Dyer, L. A. Goldberg, M. Jerrum, and R. Martin, Rapidly mixing Markov
chains for sampling contingency tables with constant number of rows, Proceedings of the 43rd
Annual Symposium on Foundations of Computer Science (FOCS 2002), 711–720.
[7] M. Dyer and A. Frieze, On the complexity of computing the volume of a polyhedron, SIAM
Journal on Computing, 17 (1988), 967–975.
[8] M. Dyer, A. Frieze, and R. Kannan, A random polynomial time algorithm for approximating
the volume of convex bodies, Journal of ACM, 38 (1991), 1–17.
[9] M. Dyer and C. Greenhill, Polynomial-time counting and sampling of two-rowed contingency
tables, Theoretical Computer Sciences, 246 (2000), 265–278.
[10] M. Dyer, R. Kannan, and J. Mount, Sampling contingency tables, Random Structures andAlgorithms, 10 (1997), 487–506.
[11] 伏見正則，乱数, 東京大学出版会，1989．

[12] V. Guruswami, Rapidly mixing Markov chains: a comparison of techniques (a survey), 2000.
(available from http://www.cs.washington.edu/homes/venkat/pubs/papers.html)[13] A. Gelman and X. Meng, Simulating normalizing constants: from importance sampling tobridge sampling to path sampling, Statistical Science, 13 (1998), 163–185.
[14] O. Häggström, Finite Markov Chains and Algorithmic Application, Cambridge University
Press, 2002.
[15] W. Höﬀding, Probability inequalities for sums of bounded random variables, Journal of American Statistical Association, 58 (1963), 13–30.
[16] 伊庭幸人, 統計学者・数理工学者のための統計物理入門―格子スピン模型とマルコフ連鎖モンテ
カルロ法を中心にして― (改訂版), Research Memorandum, 635, Institute of Statistical Mathematics, 1997.
[17] M. Jerrum, Counting, Sampling and Integrating: Algorithms and Complexity, ETH Zürich,
Birkhauser, Basel, 2003.
[18] M. Jerrum and A. Sinclair, Approximating the permanent, SIAM Journal on Computing, 18(1989), 1149–1178.
[19] M. Jerrum and A. Sinclair, Polynomial time approximation algorithms for the Ising model,
SIAM Journal on Computing, 22 (1993), 1087–1116.
[20] M. Jerrum and A. Sinclair, The Markov chain Monte Carlo method: an approach to approximate counting and integration, in Approximation Algorithm for NP-hard Problems, D.
Hochbaum, (ed.), PWS, 1996, 482–520.
[21] M. Jerrum, A. Sinclair, and E. Vigoda, A polynomial-time approximation algorithm for thepermanent of a matrix with non-negative entries, Journal of ACM, 51 (2004), 671–697.
[22] M. Jerrum, L. Valiant, and V. Vazirani, Random generation of combinatorial structures froma uniform distribution, Theoretical Computer Science, 43 (1986), 169–188.
[23] R. Kannan, L. Lovász, and M. Simonovitz, Random walks in a convex body and an improved
volume algorithm, Random Structures and Algorithms, 4 (1993), 359–412.
[24] R. Kannan, L. Lovász, and M. Simonovitz, Random walks and an O∗ (n5 ) volume algorithm
for convex bodies, Random Structures and Algorithms, 11 (1997), 1–50.
[25] R. Karp and M. Luby, Monte Carlo algorithms for enumeration and reliability problems,Proceedings of 24th IEEE Annual Symposium on Foundations of Computer Science (FOCS
1983), 56–64.
[26] R. Karp, M. Luby, and N. Madras, Monte Carlo approximation algorithms for enumeration
problems, Journal of Algorithms, 10 (1989), 429–448.
[27] D. Karger, A randomized fully polynomial time approximation scheme for the all-terminal
network reliability problem, SIAM Journal on Computing, 29 (1999), 492–514.
[28] S. Kijima and T. Matsui, Approximate counting scheme for m × n contingency tables, IEICE
Transactions on Information and Systems, E87-D (2004), 308–314.
[29] S. Kijima and T. Matsui, Polynomial time perfect sampling algorithm for two-rowed contingency tables, Random Structures and Algorithms, 29 (2006), 243–256.
[30] S. Kijima and T. Matsui, Rapidly mixing chain and perfect sampler for logarithmic separable
concave distributions on simplex, DMTCS Proceedings Series, AD (2005), 371–382.
[31] S. Kijima and T. Matsui, Approximate/perfect samplers for closed Jackson networks, Proceedings of the 2005 Winter Simulation Conference (WSC 2005), 862–868.

[32] S. Kijima and T. Matsui, Polynomial-time randomized approximation and perfect sampler
for closed Jackson networks with single servers, Mathematical Engineering Technical Reports,
2005-12, University of Tokyo, 2005.
[33] S. Kijima and T. Matsui, Randomized approximation scheme and perfect sampler for closed
Jackson networks with multiple servers, Mathematical Engineering Technical Reports, 200634, University of Tokyo, 2006.
[34] 来嶋 秀治, 松井 知己, 完璧にサンプリングしよう !, オペレーションズ・リサーチ, 50 (2005),
169–174, 264–269, 329–334.
[35] J. Kim, R. Montenegro, and P. Tetali, A near optimal bound for Pollard’s rho to solve discretelog, Proceedings of the 48th Annual Symposium on Foundations of Computer Science (FOCS
2007), to appear.
[36] L. Lovász and S. Vempala, Hit-and-run is fast and fun, Microsoft Research Technical Report,
2003-05, Microsoft Research, 2003.
[37] L. Lovász and S. Vempala, Simulated annealing in convex bodies and an O∗ (n4 ) volume
algorithms, Microsoft Research Technical Report, 2003-31, Microsoft Research, 2003.
[38] T. Matsui, M. Motoki, N. Kamatani, and S. Kijima, Polynomial time approximate/perfect
samplers for discretized Dirichlet distribution, Mathematical Engineering Technical Reports,
2006-09, University of Tokyo, 2006.
[39] X. Meng and W. Wong, Simulated ratios of normalizing constants via a simple identity: atheoretical exploration, Statistica Sinica, 6 (1996), 831–860.
[40] M. Mitzenmacher and E. Upfal, Probability and Computing, Randomized Algorithms andProbabilistic Analysis, Cambridge University Press, 2005.
[41] R. Montenegro and P. Tetali, Mathematical aspects of mixing times in Markov chains, Now
Publishers, 2006.
[42] B. Morris and A. Sinclair, Random walks on truncated cubes and sampling 0-1 knapsack
solutions, SIAM Journal on Computing, 34 (2004), 195–226.
[43] R. Motowani and P. Raghavan, Randomized Algorithms, Cambridge University Press, 1995.
[44] 大森裕浩, マルコフ連鎖モンテカルロ法の最近の展開, 日本統計学会誌, 31 (2001), 305–344.
[45] J. Propp and D. Wilson, Exact sampling with coupled Markov chains and applications tostatistical mechanics, Random Structures and Algorithms, 9 (1996), 223–252.
[46] J. Propp and D. Wilson, How to get a perfectly random sample from a generic Markov chain
and generate a random spanning tree of a directed graph, Journal of Algorithms, 27 (1998),
170–217.
[47] A. Sinclair, Improved bounds for mixing rates of Markov chains and multicommodity ﬂow,
Combinatorics, Probability and Computing, 1 (1992), 351–370.
[48] A. Sinclair, Algorithms for Random Generation and Counting: A Markov Chain Approach,
Birkhäuser, Boston, 1993.
[49] A. Sinclair and M. Jerrum, Approximate counting, uniform generation and rapidly mixing
Markov chains, Information and Computation, 82 (1989), 93–133.
[50] 田中和之 編著, 確率的情報処理と統計力学, 臨時別冊・数理科学 SGC ライブラリ 50, サイエン
ス社, 2006.
[51] L. G. Valiant, The complexity of computing the permanent, Theoretical Computer Science, 8(1979), 189–201.
[52] V. V. Vazirani, Approximation Algorithms, Springer, 2003.