計算機実験 (講義 L6)
計算機実験 (講義 L6)
藤堂眞治
2015-06-17/24
1 最適化問題
2 乱択アルゴリズムとモンテカルロ積分
1 / 36
計算機実験 (講義 L6)
最適化問題
最適化問題
目的関数 f (x) の最小値 (あるいは最大値) とその場所を求め
たい
▶
▶
連続最適化問題
離散最適化問題 ⇐ 難しい
真の (大局的な) 最小値 (最大値) を求めるのは難しい
一般的には極値を求めることしかできない
多次元では極小を囲い込むことができない
導関数を使う方法: 最急降下法、共役勾配法
使わない方法: 囲い込み法、Nelder-Mead の滑降シンプレッ
クス法、アニーリング
2 / 36
計算機実験 (講義 L6)
最適化問題
囲い込み法 (一次元の最適化)
f (a) > f (b) < f (c) を満たす 3 点の組 a < b < c の領域を狭
めていく
[a, b]、[b, c] の広い方 (例えば後者) を b から見て
0.382 : 0.618 (黄金比) に内分する点を x とする
▶
▶
f (b) > f (x) の場合: [b, c] を新しい領域にとる
f (b) < f (x) の場合: [a, x] を新しい領域にとる
もともとの b が [a, c] を 0.382 : 0.618 に内分する点だった場
合、新しい領域の幅は、どちらの場合も 0.618
最初の比率が黄金比からずれていたとしても、黄金比に収束
3 / 36
計算機実験 (講義 L6)
最適化問題
最初の囲い込み
1 点を選び、適当な ∆x を取る
左右に ∆x 動かしてみて、関数値が小さくなる方へ動く
どちらに進んでも関数値が大きくなる場合には、囲い込み
完了
小さくなった場合、その方向へ再び増えるまで ∆x を倍々に
増やしながら進む
最後の 3 点で極小値を囲い込むことができる
4 / 36
計算機実験 (講義 L6)
最適化問題
極小値をとる x の精度
実数の有効桁数を 16 桁 (ϵ ≈ 10−16 ) とする (倍精度)
真の極小 (x0 ) のまわりでテイラー展開
1
f (x) ≈ f (x0 ) + f ′′ (x0 )(x − x0 )2
2
f ′′ (x0 )/f (x0 ) が O(1) だとすると
√
|x − x0 | ∼ ϵ ∼ 10−8
以下になると、第二項の第一項に対する比が ϵ よりも小さく
なる
それ以上領域を狭めても、関数値は変化しない
5 / 36
計算機実験 (講義 L6)
最適化問題
最急降下法 (steepest-descent)
関数の微分の情報を使う
現在の点 x における勾配を計算
−∇f |i = −
∂f
∂xi
坂を下る方向にそって、一次元最適化
動いた先の勾配の方向でさらに最適化を繰り返す
関数値は単調減少 ⇒ 極小値に収束
6 / 36
計算機実験 (講義 L6)
最適化問題
308
細長い谷の場合
(a)
(Press et al 1988)
:
7 / 36
計算機実験 (講義 L6)
最適化問題
共役勾配法 (conjugate-gradient)
関数がある点のまわりで
1
f (x) ≈ c − bT x + xT Ax
2
と近似できるとする
この時、x における勾配は、連立方程式 Ax = b の「残差」の
形で書ける
−∇f = b − Ax
新しい勾配方向ではなく、それまでとは「共役な方向」に進
みたい
8 / 36
計算機実験 (講義 L6)
最適化問題
「共役な方向」とは
あるベクトル p にそった一次元の最適化が完了したとする
▶
▶
その点における p 方向の勾配は零。すなわち pT (∇f ) = 0
p 方向の勾配の値を変化させないようにしたい
次に、q にそって、x + ϵq と移動すると
δ(∇f ) = A × (ϵq) ∼ Aq
これが p に垂直であり続けるためには
pT Aq = 0
この関係が成り立つ時、p と q は「互いに共役」という
9 / 36
計算機実験 (講義 L6)
最適化問題
共役勾配法 (Conjugate-gradient)
初期条件と漸化式
p0 = r0 = b − Ax0
xn+1 = xn + αn pn
rn+1 = rn − αn Apn = b − Axn+1
pn+1 = rn+1 + βn pn
αn =
rnT rn
pT
n Apn
βn =
T r
rn+1
n+1
rnT rn
このように作ると、全ての i > j ≥ 1 について、自動的に
T
pT
i Apj = 0 ri rj = 0
10 / 36
計算機実験 (講義 L6)
最適化問題
共役勾配法 (Conjugate-gradient)
残差は互いに直交 ⇒ N 回反復すると残差は零 (完全な二次
形式の場合)
残差は負の勾配で表される ⇒ A を知らなくても ri は計算可
実際には、数値誤差により、共役性・直交性がくずれる
また、完全な二次形式ではない
しかし、最急降下法と比較すると圧倒的に速く収束
11 / 36
計算機実験 (講義 L6)
最適化問題
逆反復法による固有ベクトルの計算
f (x) の極小解は、連立一次方程式 Ax = b の解
▶
▶
連立方程式を解くのに共役勾配法を利用可
行列ベクトル積だけで計算できるので、A が疎行列の時、特
に有効
逆反復法
▶
▶
▶
近似固有値を µ とするとき、行列 (A − µI )−1 を考えると、固
有ベクトルは A と同じ、固有値は (λ − µ)−1 。
µ が十分に正確であれば、(λ − µ)−1 は絶対値最大の固有値。
行列 (A − µI )−1 を適当な初期ベクトルにかけ続けると λ に対
応する固有ベクトルに収束 (c.f. べき乗法)
実際には (A − µI )x′ = x という連立方程式を繰り返し解く
12 / 36
計算機実験 (講義 L6)
最適化問題
Nelder-Mead の滑降シンプレックス法
関数値のみ。導関数の情報を必要としない
プログラミングが簡単
収束は遅いが、安定に極小値が求まる
N + 1 個の頂点からなる N 次元の単体 (シンプレックス) を変
形しながら、極小値を探す
▶
▶
2 次元: 三角形
3 次元: 四面体
別名「アメーバ法」
13 / 36
計算機実験 (講義 L6)
最適化問題
Nelder-Mead の滑降シンプレックス法
N + 1 個の点 x0 , x1 , · · · , xN は f (x0 ) ≤ f (x1 ) ≤ · · · ≤ f (xN ) と
並べられているとする
最大値を取る点 xN を除く N 点の重心を xg とする
Nelder-Mead 法では以下のステップを繰り返す
▶
▶
▶
▶
xN の xg に関する対称な点と xN の関数値を比較し、小さい方
に移動 (反射)
x0 の関数値よりも小さくなるようであればさらに先に進む
(拡大)
xN の関数値が xN−1 のものよりもまだ大きい場合には xN を
xg に近づける (縮小)
それでも xN の関数値が小さくならない場合、x0 以外の点を
x0 に一様に近づける (収縮)
14 / 36
計算機実験 (講義 L6)
最適化問題
Nelder-Mead の滑降シンプレックス法
http://www.kniaz.net/software/RosNM.aspx
15 / 36
計算機実験 (講義 L6)
乱択アルゴリズムとモンテカルロ積分
乱択アルゴリズム (randomized algorithm)
実行中に乱数を参照してその値によって振る舞いをかえるア
ルゴリズム
ラスベガスアルゴリズム
▶
▶
乱数の出方によらず常に正しい結果を与える乱択アルゴリ
ズム
平均化効果を利用するアルゴリズム：クイックソート、最小
包含円問題
モンテカルロアルゴリズム
▶
▶
▶
▶
乱数の出方によっては誤った結果を与えることがあるアルゴ
リズム
標本を利用するアルゴリズム：最大カット問題
くじ引き型のアルゴリズム：素数性判定、関数の同一性、行
列積の検算
サンプリングアルゴリズム：マルコフ連鎖モンテカルロ
16 / 36
計算機実験 (講義 L6)
乱択アルゴリズムとモンテカルロ積分
乱択クイックソート
分割統治法に基づく再帰的なソートアルゴリズム
▶
▶
配列の中から要素を一つ選び、それより小さい要素からのみ
なる集合と大きい要素のみからなる集合の 2 つに分ける
それぞれの集合をソートし、結合する
ほぼ同じ大きさの集合に分けることができる場合の実行ス
テップ数 ∼ O(n log n)
最悪 (選んだ要素が常に最大 or 最小値) の場合のステップ数
∼ O(n2 )
分割に用いる要素を「ランダムに」選ぶ ⇒ 平均ステップ数
∼ O(n log n)
17 / 36
計算機実験 (講義 L6)
乱択アルゴリズムとモンテカルロ積分
クイックソート
18 / 36
計算機実験 (講義 L6)
乱択アルゴリズムとモンテカルロ積分
最小カット問題
「最小カット (min-cut)」= 連結グラフを二つの部分に分け
るために切らなければならない辺 (エッジ) の集合のうち、
もっとも小さいもの
“Randomized Min-Cut Algorithm”
▶
▶
▶
ランダムに辺を一つ選ぶ
辺の両端の頂点を一つにつぶす (自己ループは取り除く)
残りの頂点が二つになるまで繰り返す
上の試行を何回も繰り返して、得られたうちで最小なカット
を結果とする
頂点の数が n のグラフに関して、正解が得られる確率：少な
2
くとも n(n−1)
19 / 36
計算機実験 (講義 L6)
乱択アルゴリズムとモンテカルロ積分
最小カット問題
20 / 36
計算機実験 (講義 L6)
乱択アルゴリズムとモンテカルロ積分
モンテカルロ積分
円周率を与える公式
∫
π = lim
c→∞ 0
c
f (x) dx
2
cosh x
f (x) =
スタンダードな数値積分法: 台形公式 (一次式補間), シンプ
ソン公式 (二次式補間), etc
カットオフ c の値
▶
▶
誤差は c が大きくなると指数関数的に小さくなる
例えば c = 20 で誤差は 8.3 × 10−9 以下
21 / 36
計算機実験 (講義 L6)
乱択アルゴリズムとモンテカルロ積分
単純サンプリング
[0, c] と [0, 2] の一様分布から二次元上の点 (x, y ) を M 組
生成
f (x) の下に入った数 N をカウント
π ≃ 2c ×
平均値
4.8
3.12
3.154
誤差
1.3
0.11
0.011
2
1.5
1
y
M
100
10000
1000000
N
M
0.5
0
0
1
2
3
4
5
6
x
22 / 36
計算機実験 (講義 L6)
乱択アルゴリズムとモンテカルロ積分
統計誤差の評価
このモンテカルロ積分が実際に評価している積分
{
∫ c∫ 2
1 if y < f (x)
1
θ(x, y ) dx dy θ(x, y ) =
2c 0 0
0 otherwise
統計誤差の評価
▶
▶
▶
π
試行の成功確率 (success probability): q = 2c
一回の試行の平均値 (mean): µ = 2c × q = π
分散 (variance):
s 2 = (2c)2 q + 02 (1 − q) − µ2 = 2cπ − π 2 = 4c 2 q(1 − q)
▶
c = 20 の時:
q ≃ 0.0785 s 2 ≃ 116
23 / 36
計算機実験 (講義 L6)
乱択アルゴリズムとモンテカルロ積分
中心極限定理 (central limiting theorem)
M 回の試行のうち N 回成功する確率 (π の見積もり値が
m = 2cN/M となる確率)
N
M!
p(m = 2c ) =
q N (1 − q)M−N
M
N!(M − N)!
両辺の対数をとってスターリングの公式を使う
M
π
2c − π
log p(m) ≃
(m log + (2c − m) log
)
2c
m
2c − m
m に関して平均値 π の周りで二次まで展開
M
log p(m) ≃ − 2 (m − π)2
2s
分散 s 2 /M の正規分布
(中心極限定理)
√
統計誤差は M に反比例して減少 ⇒ 1 桁小さくするには
100 倍の計算が必要
24 / 36
計算機実験 (講義 L6)
乱択アルゴリズムとモンテカルロ積分
単純サンプリング (2)
y に関してあらかじめ積分
[0, c] の一様乱数 x を用いて
∫ c
f (x)
1 ∑
1
p(x) dx ≃
cf (xi ) p(x) =
M
c
0 p(x)
i
2
平均値
3.1
3.00
3.147
誤差
0.8
0.08
0.008
1.5
1
y
M
100
10000
1000000
0.5
0
0
1
2
3
x
4
5
6
25 / 36
計算機実験 (講義 L6)
乱択アルゴリズムとモンテカルロ積分
誤差の評価
関数 f (x)/p(x) の分散
∫ c(
∫ ∞
f (x) )2
2
2
s =
p(x) dx − π = c
f 2 (x) − π 2 = 4c − π 2
p(x)
0
0
c = 20 のとき s 2 ≃ 70.1
同じ試行回数 M の時, 誤差は
√
70.1/116 = 0.77 倍
もしくは M を 116/70.1 = 1.65 倍したのと同じ効果
積分次元は低ければ低いほど良い
26 / 36
計算機実験 (講義 L6)
乱択アルゴリズムとモンテカルロ積分
次元の呪い (curse of dimensionality)
n 次元超立方体 (1 辺の長さ 2, 体積 2n) に対する n 次元単位
球の体積の割合
q=
π n/2 /Γ( n2 + 1)
∼ (π/n)n/2
2n
n = 10 で 0.2%, n = 20 で 10−8 , n = 100 で 10−70
モンテカルロ積分で球の体積を計算しようとすると, 標準偏
差に対する平均値の割合は指数関数的に小さい
√
q
√
∼ q
q(1 − q)
次元が高くなるにつれて指数関数的に大きな M が必要と
なる
c.f. 通常の数値積分 (台形公式等) でも同様
27 / 36
計算機実験 (講義 L6)
乱択アルゴリズムとモンテカルロ積分
重点的サンプリング
(平均値が同じなら) 被積分関数の分散が小さければ小さいほ
ど良い (= 統計誤差が小さい)
サンプリングの分布 p(x) の形が f (x) に近い程良い
f (x) の値が大きい所はより頻繁にサンプリング
重点的サンプリング (importance sampling)
28 / 36
計算機実験 (講義 L6)
乱択アルゴリズムとモンテカルロ積分
重点的サンプリング
積分への寄与が大きな箇所をより重点的にサンプリング
p(x) = e −x
2
平均値
3.06
3.142
3.1412
誤差
0.06
0.006
0.0006
1.5
1
y
M
100
10000
1000000
0.5
0
0
1
2
3
x
4
5
6
29 / 36
計算機実験 (講義 L6)
乱択アルゴリズムとモンテカルロ積分
誤差のサンプル数依存性
関数 f (x)/p(x) の分散
∫ c(
f (x) )2
2
s =
p(x) dx − π 2 ≃ 2(2 + π) − π 2 = 0.414
p(x)
0
10
同じ試行回数
M の時, 誤差は
√
0.414/116 = 0.06 倍
random sampling (1)
random sampling (2)
importance sampling
1
もしくは M を 280 倍したのと同じ
error
0.1
0.01
0.001
0.0001
1e-05
102 103 104 105 106 107 108 109
M
10
30 / 36
計算機実験 (講義 L6)
乱択アルゴリズムとモンテカルロ積分
理想的な重点的サンプリング?
理想的には p(x) を f (x) に比例するように取れば良い
このとき f (x)/p(x) は定数 (分散 0) → 1 回のサンプリング
で厳密な結果が得られる???
実際には p(x) が確率密度となるように規格化条件から定数
c を決めておく必要あり
∫
∫
p(x) dx = c f (x) dx = 1
c は今欲しい答そのもの!
31 / 36
計算機実験 (講義 L6)
乱択アルゴリズムとモンテカルロ積分
乱数
自然乱数 (ハードウェア乱数)
▶
さいころ, コイン, ルーレット, 核分裂反応, 熱雑音, ショット
雑音 ...
モンテカルロシミュレーションにおける必要条件
▶
▶
▶
▶
多数の乱数が必要
ポータビリティ
生成速度
再現性
擬似乱数 (pseudo random number)
▶
▶
計算機でプログラムに従って生成
分布の一様性, 相関, 周期に注意する必要あり
32 / 36
計算機実験 (講義 L6)
乱択アルゴリズムとモンテカルロ積分
擬似乱数発生器
最も簡単な乱数発生器：線形合同法 (linear congruential
method)
xn+1 = (axn + c) mod m
例) a = 65539, c = 0, m = 2147483648 (周期 m − 1)
少しだけ異なる初期値 (x0 = 1, 2, 3) から始めた場合
1
1
0.9
0.9
0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0
20
40
60
80
100
0
99900
99920
99940
99960
99980
100000
数十ステップ進むとバラバラな振舞い ⇒ カオス的
33 / 36
計算機実験 (講義 L6)
乱択アルゴリズムとモンテカルロ積分
擬似乱数生成器における相関
合同乗算法で多次元超立方体中に「ランダムに」点を打つ
と、それらの点は全て比較的小数の等間隔に並んだ超平面の
上にのってしまう (多次元疎結晶構造)
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1
計算式に従って生成するため、必ず何らかの相関は残る
できる限り相関が少なく周期の長い、理想的な乱数の開発が
続けられている
▶
▶
現時点で最も品質が高いと考えられている乱数発生器：メル
センヌ・ツイスター
周期 219937 − 1、高速、日本製!
34 / 36
計算機実験 (講義 L6)
乱択アルゴリズムとモンテカルロ積分
乱数発生器の選び方
万能乱数発生器は存在しない
生成された乱数のもつ性質について, 様々な数学的に厳密な
証明, 多くのテスト結果がすでに存在するが, 特定のシミュ
レーションに使った場合の結果については何も保証してくれ
ない
自分で乱数発生器を「発明」してはいけない
自分で乱数発生器をプログラムしてはいけない (既存のライ
ブラリを使う)
初期化 (種の設定) を正しく行う
実際にそれらしい乱数が生成されているか, 目でみて確認
する
二種類以上の乱数発生器を使ってみて, 互いに一致する結果
が出るかどうか確認する
35 / 36
計算機実験 (講義 L6)
乱択アルゴリズムとモンテカルロ積分
実習・講義予定
実習 EX5
▶
▶
▶
特異値分解による行列の近似
最小二乗法によるフィッティング
特異な連立一次方程式
講義 L7
▶
▶
▶
マルコフ連鎖モンテカルロ法
モンテカルロ法による最適化
スパコンと計算物理
実習 EX6
▶
▶
▶
最適化問題
モンテカルロ法
レポート No.3
小グループ単位での面談スケジュール
▶
6/17 or 18 に ITC-LMS にアップロード予定 (必ず確認のこと)
36 / 36
© Copyright 2024 ExpyDoc
