Transcript Miyazaki-PBV200404

Copyright: Daisuke Miyazaki
http://www.cvl.iis.u-tokyo.ac.jp/
Numerical Recipes in C
[日本語版]
初版
宮崎大輔
2004年4月2日（金）
PBVセミナー
第10章 関数の最大・最小
Minimization or maximization of
functions
7th Physics Based Vision Seminar
2004/Apr/2
黄金分割法(golden section search)
黄金比： 0.38197
a
•
•
•
•
初め：点1,3,2で囲まれてる
点4を計算，点2を捨てる
点5を計算，点1を捨てる
点6を計算，点4を捨てる
0.23606
b
0.38197
c
d
• 現在の区間:a,b,c,d
• もしf(b)<f(c)なら，新しい区
間をa,b,cにする
• もしf(b)>f(c)なら，新しい区
間をb,c,dにする
4
7th Physics Based Vision Seminar
2004/Apr/2
逆放物線補間（逆２次補間）(inverse parabolic interpolation)
Brent(ブレント)の方法(Brent’s method)
• 3点で放物線を当てはめ，極小を求めていく
5
7th Physics Based Vision Seminar
2004/Apr/2
滑降シンプレックス法
(downhill simplex method)
•
シンプレックス(simplex)：N次元
ではN+1個の点（頂点）とそれら
を結ぶ辺，面からなる多面体
– ２次元のシンプレックスは三角形
– ３次元のシンプレックスは四面体
•
下の操作（右図参照）を適切に選
ぶことを繰り返せば必ず極小に
収束する
a.
b.
c.
d.
最悪の点の反射
最悪の点の反射と膨張
最悪の点からの１次元の収縮
最良の点への全次元の収縮
6
7th Physics Based Vision Seminar
2004/Apr/2
方向集合法(direction set methods)
Powell法
• 共役な方向(conjugate directions)とは，「互いに干
渉しない」方向のこと
• 数学的な定義: u  A  v  0
を満たすu,vを共役(conjugate)という
ただし，AはHesse(ヘッセ)行列(Hessian matrix)
7
7th Physics Based Vision Seminar
2004/Apr/2
方向集合法(direction set methods)
Powell法
•
アルゴリズム
– まず，方向集合uiを基底ベクトルに初期化する
1. 出発点をP0とおく
2. i=1,…,Nについて，Pi-1を方向uiに沿った最小に移動し，
その点をPiとおく
3. i=1,…,N-1について，ui←ui+1とおく
4. uN←PN-P0とおく
5. PNを方向uNに沿った最小に移動し，その点をP0とおく
8
7th Physics Based Vision Seminar
2004/Apr/2
方向集合法(direction set methods)
Powell法
方向集合ui(i=1,2)：
9
7th Physics Based Vision Seminar
2004/Apr/2
方向集合法(direction set methods)
Powell法
•
•
さきほどのアルゴリズムをそのまま実装すると問題
が発生する
そこで
– 各段階でu1を捨てるのではなく、関数fが最も速く減少し
た方向を捨てる
– ときには、新しい方向を全く加えない
10
7th Physics Based Vision Seminar
2004/Apr/2
最急降下法(steepest descent method)
• 点P0から出発
• 点Piから下り坂の方向
Pi+1とする
• これを必要な回数だけ行う
に移動し、その点を
11
7th Physics Based Vision Seminar
2004/Apr/2
共役勾配法(conjugate gradient method)
• 各ループで，共役な方向hiを以下のように求め，それ
に沿って直線上の最小化を行う方法
– Fletcher-Reeves(フレッチャー・リーブズ)法
(Fletcher-Reeves method)
gi 1  gi 1
i 
gi  gi
– Polak-Ribiere(ポラック・リビエール)法
(Polak-Ribiere method)
( gi 1  gi )  gi 1
i 
gi  g i
12
7th Physics Based Vision Seminar
2004/Apr/2
可変計量法(variable metric methods)
準Newton(ニュートン)法(quasi-Newton methods)
• 関数fがTaylor級数の２次までで近似できるとする
• 任意の点をxiとし，極小点をxとした場合，
が成り立つ
ここでA-1はHesse行列の逆行列である
• A-1を求めたい→反復により近似を高める
というHiを構成すればよい
13
7th Physics Based Vision Seminar
2004/Apr/2
可変計量法(variable metric methods)
準Newton(ニュートン)法(quasi-Newton methods)
• 更新式は
• DFP法[Davidon-Fletcher-Powell(ダヴィドン・フレッチャー・パウエル)法
(Davidon-Fletcher-Powell algorithm)]
• BFGS法[Broyden-Fletcher-Goldfarb-Shanno(ブロイデン・フレッチャー・
ゴルトファルプ・シャノ)法(Broyden-Fletcher-Goldfarb-Shanno
algorithm)]
14
7th Physics Based Vision Seminar
2004/Apr/2
線形計画法(linear programming)
線形最適化(linear optimization)
• x1,…,xNについて以下の目的関数
(objective function)を最大化する問題
• ただし，以下の主条件を満たす
• さらに，以下の副条件を満たす(b≧0)
• （実行）可能ベクトル (feasible vector)：上
記の主条件と副条件を満たすx1,…,xN
• 最適可能ベクトル(optimal feasible
vector)：目的関数を最大化する可能ベク
トル
15
7th Physics Based Vision Seminar
2004/Apr/2
制限標準型(restricted normal form)の
シンプレックス法(simplex method)
1. まず，x2とx3を0に設定する
2. [z,x2]が2，つまり正なのでx2を増
やせばzが増える
• x3を増やすとzは減る
• x3は0なのでそれ以上減らせな
い(x3≧0の条件に反する)
•
これらの条件をシンプレックス
表(タブローtableau)で表す
3. [x1,x2]が-6，つまり負なので増や
しすぎるとx1≧0の条件に反する
• x2を増やしてもx4は増えるので
x4≧0は満たされる
4. x2は2÷6=1/3まで増やせる
• x2=1/3となる
• zは1/3×2=2/3増える
• x1は0になる
5. 次スライドにつづく
16
7th Physics Based Vision Seminar
2004/Apr/2
制限標準型(restricted normal form)の
シンプレックス法(simplex method)
6. x1とx2を入れ替えたタブ
ローを作る
2 1
11
Maximize z   x1  x3
3 3
3
1 1
1
x2   x1  x3
3 6
6
1
7
x4  9  x1  x3
2
2
7. [z,x1]と[z,x3]が負なのでzを
これ以上増やせない
•
•
•
終了
z=2/3
x2=1/3, x4=0, x1=x3=0
17
7th Physics Based Vision Seminar
2004/Apr/2
一般のシンプレックス法(simplex method)
• 人工変数（人為変数）(artificial
variables)ziを導入する
•
• 以下の補助目的関数(auxiliary objective
function)を最大化させる
スラック変数(slack
variables)yiを不等式の
片側に加えて，不等式の • zi≧0なので，z’が最大となるのは，zi=0の
制約条件を等式にする
とき
– zi=0となったらスラック変数付きの条件と
同じになる
– zi=0とならないziがあったら解なし
• zi=0となったら，元々の目的関数を最大
化させ，解を求める
18
7th Physics Based Vision Seminar
2004/Apr/2
アニーリング（焼きなまし）法
(simulated annealing)
• 焼きなまし(annealing)
–
–
–
–
高温では液体中の分子は自由に運動
ゆっくり冷やすと原子は完全に整列し，純粋な結晶になる
ゆっくり冷やすと自然は最小エネルギー状態を見つける
急冷するとエネルギーの高い多結晶質・非晶質になる
• Boltzmann(ボルツマン)の確率分布によれば，系は
坂を下るばかりでなく，ときには坂を登るが，温度が
低いほど坂を登る可能性が小さい
• Metropolis(メトロポリス)の方法
– 下り坂なら必ず下る
– 上り坂なら確率p=exp[-(エネルギーの差)/kT]で坂を登る
19
7th Physics Based Vision Seminar
2004/Apr/2
アニーリング法は巡回セールスマン問題(travelling
salesman problem)を事実上解決した
• NP完全(NP-complete):計算時間はO(eN)
• 組替え
– 経路の一部を切り離し，同じ場所に逆向きに挿入
– 経路の一部を切り離し，別の２地点間に挿入
• 目的関数は巡回の道のり
– 重みをつけることにより，川を渡るのをなるべく減らすよう
な巡回を見つけることもできる
• アニーリングスケジュール
– これはいろいろやってみるしかない
– ここでは：配置を100N回生成するか，組替えが10N回成
功するか，ごとにTを10%減らす
20
7th Physics Based Vision Seminar
a.
b.
c.
2004/Apr/2
ほぼ最短経路
川越えのペナルティを大きくした
川越のペナルティを負にした
21
第11章 固有値問題の数値計算法
Eigensystems
7th Physics Based Vision Seminar
2004/Apr/2
固有値問題
• AはN×N行列，xは固有ベクトル(eigenvalue)，λは
固有値(eigenvalue)，とすると
• 両辺にτxを加えると (A  1)x  (  )x
となる．つまり，定数に単位行列を掛けたものを行列
に加えると，行列の固有値はτだけシフトする
– シフト，は固有値・固有ベクトルを求める際に，収束を速く
するときに使う
23
7th Physics Based Vision Seminar
2004/Apr/2
定義と基本事項
•
のとき，その行列は対称(symmetric)とよぶ
•
のとき，つまり行列がそのHermite(エルミート)共役(Hermitian conjugate,
転置行列の各成分を複素共役にしたもの)と等しいとき，
Hermite(Hermitian)または自己随伴(self-adjoint)とよぶ
•
のとき，つまりAT=A-1のとき，直交(orthogonal)という
• 同様に，A†=A-1のとき，ユニタリ(unitary)という
•
のように，Hermite共役と交換する（可換である, commute)とき，正規
(normal)という
• 実数の行列では，Hermite≡対称≡正規，ユニタリ≡直交≡正規
24
7th Physics Based Vision Seminar
2004/Apr/2
行列の対角化
• 行列Aの相似変換(similarity transform)はある変換
行列Zにより以下のように変換されることをいう
行列の固有値は不変
• 固有値を求めるには
のような相似変換でAを対角行列へ変換する
すると固有ベクトルは
XRは固有ベクトルを列にして作られた行列
25
7th Physics Based Vision Seminar
2004/Apr/2
実対称行列のJacobi法(Jacobi method)
•
Jacobi回転(Jacobi rotation)
• a’pq=0となる回転角φは
tan 
•
この回転により非対角要素を0にしていく
これでAを変換:
※Aは対称
sgn( )
|  |   2 1

aqq  a pp
2a pq
• 最終的に対角行列Dが求ま
る（対角要素が固有値）
Vは
Piは個々のJacobi回転行列
で，Vの列は固有ベクトル
• 左上から右下に消していく
P12,P13,…,P1n,P23,P24,…
26
7th Physics Based Vision Seminar
2004/Apr/2
Householder(ハウスホルダー)法
• 対称行列Aを３重対角にしたい
• 最初のHouseholder行列P1をA
の左にかけると
• Aの左右にかける(なお，
Householder行列は直交行列，
PT=P)と1行目,1列目が３重対角
になる
• 次のHouseholder行列P2により2
行目,2列目も３重対角になる
• これをn-2回続けるとAは３重対角
になる
• ３重対角行列の固有ベクトルを求
める方法があるので，それを適用
し，その固有ベクトルに以下の変
換を適用して，Aの固有ベクトル
が求まる
27
7th Physics Based Vision Seminar
2004/Apr/2
QR法(QR method)，QL法(QL method)
•
•
•
•
•
•
•
•
•
•
•
任意の実行列Aは，Qを直交行列，Rを上三角行列とし，
と分解できる
A’を定義:
Qは直交行列なので R  Q1A  QT A
よってA’はAを直交変換したもの：
対称行列，３重対角行列，Hessenberg行列は，QR変換しても対称行列，３重対
角行列，Hessenberg行列になる
QL法の場合(Lは下三角行列):
QL法は，以下を計算していく
s→∞でAsは下三角行列かそれに似たものに近づく．固有値は対角要素に絶対値
の昇順に並ぶ．
QL法・一般の行列：O(n3), Householder変換を使う
QL法・Hessenberg行列：O(n2)
QL法・３重対角行列：O(n), Jacobi回転を使う, 対角行列（３重対角かつ下三角）に
なるので固有ベクトルが求まる
28
7th Physics Based Vision Seminar
2004/Apr/2
Hermite行列
• 複素行列では，Hermite行列に対して，Jacobi法や，
Householder法とQL法，の複素数版で固有値と固
有ベクトルを求めることができる
• また，もし
がHermite行列なら，n×n
の複素数の固有値問題
は2n×2nの実数の問題
として固有値と固有ベクトルを求めることができる
29
7th Physics Based Vision Seminar
2004/Apr/2
Hessenberg(ヘッセンベルク)形への変換
• 今までは対称行列
• 非対称行列は，まず，Hessenberg形に変換
• 上Hessenberg(Hessenberg)行列は，
例
• Householder法かGauss消去法に似た手法を使う
30
7th Physics Based Vision Seminar
2004/Apr/2
実Hessenberg行列のQR法
• 非対称行列の場合，複素数の固有値を持ちうる
→シフトを使うとうまく避けられるし，収束が速くなる
→その方法の詳細は省略する
• As2  Qs1  Qs  As  QTs  QTs1  Pn1 P2  P1  As  P1T  P2T PnT1
QR法のステップ２回分
Q:直交
P:Householder行列
As,As+2:上Hessenberg行列
• 右下から左へ２番目の要素an,n-1が限
りなく0に近づいたら，右下の要素ann
が固有値
→行n，列nを消し，次の固有値を探す
→これを繰り返していく
31
7th Physics Based Vision Seminar
2004/Apr/2
逆反復法(inverse iteration)による固有値の改
良と固有ベクトルの計算
•
正しい固有値と固有ベクトルでは
つまり(A  1)  x  x
1.
をyについて解く．
•
•
2.
k回目の反復での固有ベクトル・固有値:bk, τk
yはbkよりも，固有値τkに対応する固有ベクトルに近づく
で固有ベクトルを更新
3. 収束してないなら1に戻る
4.
で固有値を更新
5. 収束してないなら1に戻る
• ステップ1を簡単に解ける形にAを変形する必要がある
→この変形はQL法よりも4倍程度以上効率が悪い
→固有値が高精度で求まっていて，小数の固有ベクトルだ
けを求めたい時に使う
32
第13章 データの統計的記述
Statistical description of data
7th Physics Based Vision Seminar
2004/Apr/2
メディアン，平均，分散，標準偏差，平均偏差
• メディアン（中央値）(median):データをソートして，真ん中の要素を取る
（データの個数が偶数のときは，真ん中２つの値の平均を使う）
• 平均値(mean):
• 分散(variance):
• 標準偏差(standard deviation):
• 平均偏差(average deviation), 平均絶対偏差(mean absolute
deviation):
ロバスト(robust)な(外れ値への依存が少ない，頑健な)推定量
34
7th Physics Based Vision Seminar
2004/Apr/2
歪度と尖度
• 歪度(skewness):
• 尖度(kurtosis):
分布の非対称さを表す
–
–
–
–
尖度が正：急尖(leptokurtic)，とがっている分布
尖度が負：緩尖(platykurtic)，丸い分布
その中間：中尖(mesokurtic)
正規分布は尖度が0
35
7th Physics Based Vision Seminar
2004/Apr/2
連続データの最頻値の推定
•
最頻値・モード(mode)をヒストグラム(histogram)で
求める方法もあるが，あまりよくない
1. まず全データxj(j=1,…,N)を昇順にソート
2. 「窓の大きさ」3以上の整数Jを決める
3. i=1,…N-Jで以下を計算
J
1

p xi  xi J  
2
 N xi J  xi 
4. 最大となる1/2(xi+xi+J)を最頻値とする
36
7th Physics Based Vision Seminar
2004/Apr/2
二つの分布が同じ平均値・分散を持つかどうか
の検定
• Studentのt(Student’s t):２つの分布の平均値の差の有意性
(significance)を調べるための統計量
• 分散が同じで平均値が異なる場合
NA,NBは各標本のデータ数
• Studentの分布(Student’s distribution)A(t|ν):
• 統計量tは上の値を用い，自由度νはNA+NB-2を使う
• A(t|ν)>0.99のとき，二つの平均値は有意に異なる
37
7th Physics Based Vision Seminar
2004/Apr/2
二つの分布が同じ平均値・分散を持つかどうか
の検定
• Studentのt(Student’s t):２つの分布の平均値の差
の有意性(significance)を調べるための統計量
• 分散が異なり平均値が異なる場合
ν=
のときのStudentのt分布に従う
38
7th Physics Based Vision Seminar
2004/Apr/2
二つの分布が同じ平均値・分散を持つかどうか
の検定
• Studentのt(Student’s t):２つの分布の平均値の差
の有意性(significance)を調べるための統計量
• 対になった標本の場合
– 例：２人の就職志望者，10人の採用委員．10人のつけた
点の平均値に有意な差があるか？
•
ν=N-1
Nは各標本のデータ数．同じiどうしが対
39
7th Physics Based Vision Seminar
2004/Apr/2
二つの分布が同じ平均値・分散を持つかどうか
の検定
• F検定(F-test):二つの標本の分散が異なるかどうか
• F≫1,F≪1のとき，非常に有意
• F=両分散の比
ν1=N1-1
ν2=N2-1
有意性Q(F|ν1,ν2):
40
7th Physics Based Vision Seminar
2004/Apr/2
カイ２乗検定(chi-square test)
• 離散分布間の差異の検定
• 事象がi番のビンに入った回数（観測度数）をNi，何らかの既
知の分布を仮定して導かれる期待度数をniとする
• Niは整数，niは必ずしも整数ではない
• カイ２乗統計量
• χ2が大きいと，Niの分布とniの分布は異なる
• 有意性：カイ２乗確率関数(chi-square probability
function)Q(χ2|ν)
ν=NB-1
NBはビンの数
41
7th Physics Based Vision Seminar
2004/Apr/2
カイ２乗検定(chi-square test)
• ２組の離散分布間の差異の検定
• Ri:1組目のデータ，ビンiの度数
Si:2組目のデータ，ビンiの度数
• カイ２乗統計量
• あるバードウォッチャーが観察した鳥の数を種ごとに分類した
とき，今年と去年の分布が同じかどうか知りたいとする．この
場合，各ビンは種に対応する．
– もしバードウォッチャーが各年に見た最初の1000羽をデータとするな
ら，自由度はNB-1
– もし任意の日付を選んで各年について同じ日付をとり，その日に見た
すべての鳥をデータにするなら，自由度の数はNB
– 後者では年によって鳥の全数が変わったかどうかも検定することにな
る．これが余分の自由度の起源である．
42
7th Physics Based Vision Seminar
2004/Apr/2
Kolmogorov-Smirnov(コルモゴロフ-スミルノフ)
検定(Kolmogorov-Smirnov test)
• 連続分布間の差異の検定
• K-S統計量
データSN(x)と既知の累積分布関数P(x)との
比較
二つの異なるデータSN1(x),SN2(x)の比較
• 有意確率
– 1分布の場合N=Ne
– 2分布の場合
• なお，
これは単調増加関数で，
43
7th Physics Based Vision Seminar
2004/Apr/2
２分布の分割表による解析
• 二つの分布の連関の度合(measures of association)を知り
たい
• 連関がある→赤なら男性，男性なら赤，と言える状態
44
7th Physics Based Vision Seminar
2004/Apr/2
カイ２乗に基づく連関の測度
• Nij:第１の変数xがi番目の値をと
り，第２の変数yがj番目の値をと
る事象の個数
• N:Nijの合計つまり事象の総数
• nij:Nijの期待値




  Nij   Nij    Nij   Nij 
 j
 i
 i
  j




nij 

N
 Nij
i
j
• Cramer(クラメール)のV
(Cramer’s V)
I,Jは行数，列数
Vは0以上1以下，連関がないとき
0，完全な連関のとき1
• 連関係数(contingency
coefficient)C
• カイ２乗統計量
0から1の値，1にはならない，上
限の値がI,Jに依存するので，異
なる大きさの表を比較できない
45
7th Physics Based Vision Seminar
2004/Apr/2
エントロピー(entropy)に基づく連関の測度
• 不確定性係数(uncertainty
coefficient)
• xのエントロピー
• yのエントロピー
• xとyの結合エントロピー
• xが与えられたときのyのエントロピー
• xとyに連関がない
→xとyが独立
→H(x,y)=H(x)+H(y)
→U(x,y)=0
• xとyに連関がある
→xとyが従属
→H(x,y)=H(x)=H(y)
→U(x,y)=1
• yが与えられたときのxのエントロピー
• ただし
Nij
pij 
N
pi. 
N
j
N
ij
p. j 
N
ij
i
N
46
7th Physics Based Vision Seminar
2004/Apr/2
線形相関係数(linear correlation coefficient), 積率相関係数(productmoment correlation coefficient), Pearson(ピアソン)のr(Pearson’s r)
•
•
線形相関係数r
rは-1～1, r=1完全な正の相関(正の傾き
の直線), r=-1完全な負の相関(負の傾き
の直線), r=0無相関(uncorrelated)
有意確率その１(N>500)
•
有意確率その２(Nが小さい)
Studentのt分布
ν=N-2
•
有意確率その３(N≧10)
平均値
rtrueは母集団の（真の）相関係数
– 測定値のrと仮説値rtrueが異なるか調べる
ための有意確率
この値が小さいと相関が強い
値が小さいほど異なる
– 二つの相関係数r1,r2の差の有意確率
47
7th Physics Based Vision Seminar
2004/Apr/2
ノンパラメトリック相関(nonparametric
correlation), 順位相関(rank correlation)
• xiをソートして順位をつける
– 同じ順位(タイ,tie)になったときは，平均順位(midrank)を使
う
• 例：[100,85,85,85,12,12,7]→[1,3,3,3,5.5,5.5,7]
• yiも同じ
• 相関を調べる
48
7th Physics Based Vision Seminar
2004/Apr/2
Spearman(スピアマン)の順位相関係数(Spearman’s
rank-order correlation coefficient)
• xiの順位をRi, yiの順位をSi
• Spearmanの順位相関係数
• 有意性
ν=N-2
のStudentのt分布
• 順位差の２乗和
• 有意性
1


D

D


Var(
D
)

erfc

2




平均
分散
49
7th Physics Based Vision Seminar
2004/Apr/2
Kendall(ケンドール)のτ(Kendall’s tau)
• xi,yiには，値が大きい・小さい・等しいの３種類の状態を当てはめる(３種類
である必要はない)
• N個のデータ点(xi,yi)の対の組み合わせは1/2N(N-1)通り
• 対で比較したとき
–
–
–
–
–
同順(concordant):xi大・yj大，xi小・yj小
逆順(discordant):xi大・yj小，xi小・yj大
余分なy:xi等・yj大，xi等・yj小
余分なx:xi大・yj等，xi小・yj等
なし:xi等・yj等
• Kedallのτはこれらの回数を数えて
-1～1，1順位が完全に一致，-1完全に逆順
• 有意性
1




Var( ) 
erfc

2




50