Oct. 25, 2007
乱択近似(randomized approximation)

1. 離散アルゴリズム

1-1. MCMC法における近似精度保証

マルコフ連鎖モンテカルロ法(Markov chain Monte Carlo)*来嶋 秀治京都大学 数理解析研究所

1: イントロダクション
モンテカルロ法とは

3モンテカルロ法 – 円周率の計算を例に
円周率 π の値 (=3.141...)を求めよ．

π = (半径1の円の面積)

14モンテカルロ法 – 円周率の計算を例に
円周率 π の値 (=3.141...)2xxxを求めよ．

xxx

π = (半径1の円の面積)π=4x(円の中の点の数)(打った点の数)

1xxxxxxxx

5

モンテカルロ法 – 円周率の計算を例に
円周率 π の値 (=3.141...)2xxxを求めよ．

xxxπ = (半径1の円の面積)
モンテカルロ法
乱数を使用して，•比• 平均

1xxxxxxなどを近似的に計算する手法．

xx2次元領域の面積計算
モンテカルロ法
乱数を使用して，•比• 平均などを近似的に計算する手法．
(の面積) =672次元領域の面積計算
モンテカルロ法
乱数を使用して，•比xxx• 平均xなどを近似的に計算する手法．
(の面積) = (xxの面積)x(領域内の点の数)(打った点の数)xxxxxx

2次元領域の面積計算
モンテカルロ法
乱数を使用して，•比8重積分に強いという特徴がある．
(Wikipediaより)xx• 平均xなどを近似的に計算する手法．
(の面積) = (xxxの面積)x(領域内の点の数)

x(打った点の数)⇒ 高次元体の体積計算
xxxxx

9

高次元の呪い – 球を例にn次元球の体積n次元ユークリッド球の体積1外接n次元超立方体102

高次元の呪い – 球を例にn次元球の体積n次元ユークリッド球の体積上の一様ランダムな点がに入る確率1< 10-60 (n=100の時)

高次元の呪い – 球を例にn次元球の体積n次元ユークリッド球の体積

112⇒ ふっくらした高次元体• 凸•に内接する(全ての面に接する)上の一様ランダムな点がに入る確率1< 10-60 (n=100の時)⇒ (素朴なモンテカルロ法では)多項式時間での(近似)体積計算は困難！

12数え上げ，積分の計算困難性
n次元凸体の体積の計算 => #P困難
乱択近似(randomized approximation)#P: 数え上げ問題に対する計算量クラス．
#P困難 ⊂ NP困難FPRAS (Fully Polynomial-time Randomized Approximation Scheme)(全多項式時間乱択近似計算法)
A: 真の値

cf. FPTAS

Z: 近似値
MCMC法による解決
マルコフ連鎖を用いたランダム生成法 + モンテカルロ法

MCMC法における近似精度保証 – Table of Contents
1. イントロダクション
2. マルコフ連鎖を用いたランダム生成法
3. 高次元凸体の体積計算
4. 離散構造と数え上げ
4.1. ナップサック解の数え上げ
4.2. 2元分割表の数え上げ
5. パーマネントの計算
6. まとめ13

MCMC法マルコフ連鎖を用いたランダム生成法 + モンテカルロ法

2: マルコフ連鎖を用いたランダム生成法

15

マルコフ連鎖の定常分布• マルコフ連鎖 M （エルゴード的）¾ 状態空間：Ω （有限）¾ 推移規則： 推移確率行列 P によって推移する
Pij = Pr(推移 i → j ) = P( i, j)

0.2
s10.2

s20.2

0.4
0.4 0.1
0.60.3
状態遷移図

s30.6

16マルコフ連鎖の極限分布
任意の初期状態推移を繰り返す
状態を出力
状態空間

xxx・・・

xx

漸近的に定常分布「マルコフ連鎖を用いたサンプリング法」のアイデア
1. 目的の分布を定常分布にもつマルコフ連鎖を設計する。
2. 十分な回数推移させて、定常分布からサンプリングする。

状態空間

xxx・・・

xx

漸近的に定常分布17

MCMC法における近似精度保証 – Table of Contents
1. イントロダクション
2. マルコフ連鎖を用いたランダム生成法
3. 高次元凸体の体積計算
4. 離散構造と数え上げ
4.1. ナップサック解の数え上げ
4.2. 2元分割表の数え上げ
5. パーマネントの計算
6. まとめ18

3: 高次元凸体の体積計算[Dyer, Frieze, Kannan 1991]
再帰的なモンテカルロ法

問題: 高次元凸体の体積計算
K: n次元凸体
体積V(K)を計算

2021

仮定: ふっくらとした凸体
K: n次元凸体体積V(K)を計算
仮定

R

Kはふっくらとしている
1R = O(nc)K上の一様ランダム生成 – Ball walk

[Kannan, Lovasz, Simonovitz 1993]
22

Bal-walk
X: 現在の点
X’: 次の時刻の点
Step 1.
:

中の一様ランダムな点を選ぶ．
Step 2.

定常分布はK上の一様分布

xx’

23

K上の一様ランダム生成
x(ふっくらとした凸体に対して)
Ball-walkを用いて，

x一様ランダムサンプリング可能

xxxx
どうやって，体積を計算するのか？

xxxx
再帰的アルゴリズム + モンテカルロ法

xxxx

24

再帰のアイデア
K: n次元凸体
体積V(K)を計算

R25

再帰のアイデア
K: n次元凸体
体積V(K)を計算

Rr1

26

再帰のアイデア
K: n次元凸体
体積V(K)を計算

Rr1
モンテカルロ法で計算

27

再帰のアイデア
K: n次元凸体
体積V(K)を計算

xxxxxxxxx

モンテカルロ法で計算

xxxxx

28

再帰のアイデア
K: n次元凸体
体積V(K)を計算

xxxxxxxxx

モンテカルロ法で計算
K1: n次元凸体
体積V(K1)を計算

xxxxx

再帰のアイデア
K: n次元凸体
体積V(K)を計算

モンテカルロ法で計算
K1: n次元凸体
体積V(K1)を計算

29

30

再帰のアイデア
K: n次元凸体
体積V(K)を計算

モンテカルロ法で計算
K1: n次元凸体
体積V(K1)を計算

K1はふっくらとしている

再帰のアイデア
K1: n次元凸体
体積V(K1)を計算

31

再帰のアイデア
K1: n次元凸体
体積V(K1)を計算

32

33

再帰のアイデア
K1: n次元凸体
体積V(K1)を計算

xxxxxxモンテカルロ法で計算

xxxx

34

再帰のアイデア
K1: n次元凸体体積V(K1)を計算

xxxxxxモンテカルロ法で計算
先程の再帰式に代入すると，

xxxx

35

再帰式

(ただし，rq = 1)が求まる！

36

再帰の回数
再帰の回数 => 少なくしたい
比=> 大きくしたい

比(再帰の回数)37

理論的な近似精度保証に向けて
Ball-walk=> 定常分布は一様．
十分な回数の推移は何回？
xx’

理論的な近似精度保証に向けて[Kannan, Lovasz, Simonovitz 1993]
38

Ball-walkx=> 定常分布は一様
xx

十分な回数の推移は何回？
xモンテカルロ法• 真の比 ρ• 近似比 ρ’ (分布の誤差に起因)|ρ’ – ρ| < εmixing time:xx多項式回の推移でx“比” の ε 近似が可能
xxxx
誤差が最大となる部分集合のとり方に対して，誤差 ε が保証される推移の回数．

xxx

39

3章のまとめ
効率的な再帰的モンテカルロ法を設計するには？
1. 帰着先の割合が大きい
2. 再帰の回数が小さい
3. 帰着先のランダムサンプリングが容易

離散構造の場合は？

MCMC法における近似精度保証 – Table of Contents
1. イントロダクション
2. マルコフ連鎖を用いたランダム生成法
3. 高次元凸体の体積計算
4. 離散構造と数え上げ
4.1. ナップサック解の数え上げ
4.2. 2元分割表の数え上げ
5. パーマネントの計算
6. まとめ40

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

4.1: 0-1ナップサック解の数え上げ

0-1ナップサックの許容解集合
Given
正数 a1,a2,…,an, b.

42

0-1ナップサックの許容解集合
Given
正数 a1,a2,…,an, b.
目的|Ωn|の計算#P完全43

[Morris Sinclair 2004]
44

0-1ナップサックの許容解集合上のマルコフ連鎖
Given
正数 a1,a2,…,an, b.
マルコフ連鎖X ∈ Ωn : 現在の状態
X’ : 次の時刻の状態

Step 1.
添え字 i ∈{1,…,n}を一様ランダムに選ぶ．

Step 2.
Xi の値を変更してもナップサック解なら，確率1/2でその解に更新．
それ以外の場合は X’ = Xとする．

定常分布は Ωn 上の一様分布

45

一様ランダムサンプリング

xxxxxxxxx

46

再帰のアイデア

xn = 1xn = 0[注] このとき，xxn = 1
xn = 0xxxxxxxx

47

再帰のアイデア

関係式

xn = 1xn = 0x
モンテカルロ法で計算可能
xn = 1xxxn = 0xxxxxx再帰式

48

自己帰着可能性 (self-reducibility)
1.
2.
3.

が，ある程度大きい．
再帰式に基づく再帰的モンテカルロ法

上の効率的なサンプリング法が鍵49

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

4.2: 2元分割表の数え上げ (問題紹介)– 輸送多面体中の整数格子点の数え上げ．

51

2元分割表
2元分割表
9 m 行 n 列の表 (mxn行列)
9 各セルには非負整数が入る
9 （与えられた）周辺和を満たす

1218
5 4 3 7 5 6 30

cf. Hitchcock 輸送問題の整数解• 供給地 2 / 需要地 6

52

2元分割表
2元分割表
9 m 行 n 列の表 (mxn行列)
9 各セルには非負整数が入る
9 （与えられた）周辺和を満たす

1218
5 4 3 7 5 6 30

5 4 3 0 0 0 120 0 0 7 5 6 185 4 3 7 5 6 30状態A

0 0 0 1 5 6 12
5 4 3 6 0 0 18
5 4 3 7 5 6 30状態C

4 3 1 3 1 0 12
1 1 2 4 4 6 18
5 4 3 7 5 6 30状態B与えられた周辺和を満たす2元分割表の個数を求める問題⇒ #P完全 (行数を2に固定しても)

53

2元分割表
2元分割表
9 m 行 n 列の表 (mxn行列)
9 各セルには非負整数が入る
9 （与えられた）周辺和を満たす

1218
5 4 3 7 5 6 30

• 2xn分割表のFPRAS [Dyer, Greenhill 2000]• m x n分割表の再帰構造 [K, Matsui 2002]• m*x n分割表の多項式時間サンプリング法 [Cryan et al 2002]

未解決
m x n分割表の多項式時間サンプリング法

MCMC法における近似精度保証 – Table of Contents
1. イントロダクション
2. マルコフ連鎖を用いたランダム生成法
3. 高次元凸体の体積計算
4. 離散構造と数え上げ
4.1. ナップサック解の数え上げ
4.2. 2元分割表の数え上げ
5. パーマネントの計算
6. まとめ54

5: パーマネントの計算– 2部グラフの完全マッチングの数え上げ56

パーマネント
A: 正方行列 (nxn)(Sn: n次の対称群)cf.

各項は2部グラフの完全マッチングに対応

57

2部グラフの完全マッチング
A: 0-1正方行列 (nxn)#P完全 [Valiant 1979](Sn: n次の対称群)
abcd

E

aE

FbF

GcGHdH

パーマネント

完全マッチングの個数

58

完全マッチング上のマルコフ連鎖の設計は難しい
Ω := 完全マッチング全体の集合
X ∈ Ω: 現在の状態

次の時刻の状態？？？

aEbFcGdH59

状態空間: n マッチング+ n-1マッチング [Broder 1986]
Ξ := n マッチング + n-1マッチングaE

状態空間: Ξb

F

推移規則:

cGd

H

枝 e を一様ランダムに選ぶ．
• X が n マッチングの場合• X が n-1 マッチングの場合
定常分布: Ξ 上の 一様分布

60

状態空間: n マッチング+ n-1マッチング [Broder 1986]
Ξ := n マッチング + n-1マッチングaE

状態空間: ΞbF

推移規則:

cGd

H

枝 e を一様ランダムに選ぶ．
• X が n マッチングの場合• X が n-1 マッチングの場合
定常分布: Ξ 上の 一様分布
rapidly mixing [Jerrum, Sinclair 1993]
棄却サンプリング？

61

棄却サンプリングは非効率的
k-hexagon

・・・

n-1マッチングの個数 ≅ 2k ≅ 2n/6
完全マッチングの個数 = 1
棄却サンプリングでは，完全マッチングを得るのに指数時間かかる

FPRAS [Jerrum, Sinclair, Vigoda 2001]“アニーリング”を利用したアルゴリズム

62

63

FPRAS [Jerrum, Sinclair, Vigoda 2001]
問題点: (一般には) n-1マッチングの比率が大きい．
aE

マッチングに “重み” をつけて，bF重みに比例する確率分布を考える．

c

G

マッチング重み関数 w : “hole” -> R+d

H

アイデア 1.

hole: n-1マッチングで，マッチングされていない頂点対=
X が完全マッチング

Xがn-1マッチング,[u,v]がhole

重み関数 w の計算が困難．=> “アニーリング” を利用

FPRAS [Jerrum, Sinclair, Vigoda 2001]
アイデア 2.
完全2部グラフ + 枝活性度=> 活性度を “アニーリング” する．

64aEbF

cGdH• 状態空間: nマッチング+n-1マッチング• マッチング重み関数: 現在の枝活性度に関して 枝活性度の積=X が完全マッチング≅
Xがn-1マッチング,[u,v]がhole• 定常分布: (枝活性度の積) x (マッチング重み) に比例
枝活性度を “アニーリング” しながら，重み関数 w を算定．
<= “平均” の近似計算 by モンテカルロ法

6: Concluding Remarks

66

まとめ
MCMC法と近似精度保証
key word: 自己帰着可能性 (self-reducibility)
1. 帰着先の割合が大きい
2. 再帰の回数が小さい
3. 帰着先のランダムサンプリングが容易(理論的な厳密さを追求しなければ)アイデアはそれほど難しくない．
ぜひ，MCMC法を使ってみましょう！

67

今後の課題
未解決問題 [多項式時間収束のマルコフ連鎖]
9 2元分割表 (m x n)
9 マトロイドの基
9 分配束構造9 グラフ中の森

離散の難しさ (連続との比較)• 非等方性• スケーリングの難しさ

9 コーダルグラフサンドイッチ
関連する話題 統計学: bridge sampling, path sampling 統計物理: レプリカ交換法, 拡張アンサンブル法
 最適化: annealing

The endThank you for the attention.