Transcript Miyazaki-PBV200407

東京大学
The University of Tokyo
池内研究室
Computer Vision Laboratory
Numerical Recipes in C
[日本語版]
初版
宮崎大輔
2004年7月1日（木）
PBVセミナー
9th Physics Based Vision Seminar
2004/July/1
NUMERICAL RECIPES in C [日本語版]
ニューメリカルレシピ・イン・シー C言語による数値計算のレシピ
• William H. Press, Saul A.
Teukolsky William T.
Vetterling, Brian P.
Flannery
• 丹慶勝市，奥村晴彦，佐藤
俊郎，小林誠
• 1993
• 技術評論社
• ISBN4-87408-560-1
2
東京大学
The University of Tokyo
池内研究室
Computer Vision Laboratory
第16章 ２点境界値問題
Two point boundary value problems
9th Physics Based Vision Seminar
2004/July/1
２点境界値問題
(two point boundary value problem)
• 微分方程式：
• 境界条件：
を解が点x1,x2で満たす
ねらい撃ち法，射撃法
(shooting method)
緩和法
(relaxation method)
4
9th Physics Based Vision Seminar
2004/July/1
ねらい撃ち法
• 点x1の境界条件を満たす任
意のベクトルをVとする
• Vが決まると点x1での関数
の値が分かり，それを初期
値としてRunge-Kutta法な
どを使って計算していくと，
点x2での関数の値が分かる
• このとき，点x2での境界条
件が満たされているか調べ
る
• 点x2での「ずれのベクトル」をFと
する
– F=0なら境界条件を満たしている
– F≠0なら境界条件を満たしていな
い
• F=0となるようにNewtonRaphson法で解を求める
• 計算は以下の通り
– もし満たされていたら計算終
了
– もし満たされていなかったらV
を変化させ，もう一度RungeKuttaを計算しなおす
5
9th Physics Based Vision Seminar
2004/July/1
適合点へのねらい撃ち
(shooting to a fitting point)
• x1からx2へと解く代わりに，まずx1からx1とx2の間に
ある点xfへと解き，次にx2からxfへと（逆向きに）解い
ていく
6
9th Physics Based Vision Seminar
2004/July/1
緩和法(relaxation methods)
• 常微分方程式を差分方程式
(finite difference equations,
FDE)に換える
例：
• 点kでΔyを計算：
点1でΔyを計算：
• k=1,2,…,Mの点の格子（メッシュ），
各点xk，最初の境界点x1，最後の
境界点xM
• 点xkでの従属変数をykとする
• 点kで成り立つ：
点1で成り立つ：
点Mで成り立つ：
• yの初期値を与え，各反復で新し
い値をy+Δyとする
点MでΔyを計算：
• この連立１次方程式をGaussの
消去法と後退代入で解くことがで
きる
7
9th Physics Based Vision Seminar
2004/July/1
メッシュ点の自動割当て
• 解が急激に変化する所は格子点を多く取り，その他
のところは少なくしてよい→動的に格子（メッシュ）を
割り当てる(allocate)
8
9th Physics Based Vision Seminar
2004/July/1
内点での境界条件，特異点の処理
•
でD=0の場合→特異点
• 物理の問題では，D→0のときN→0となることが多い
– 正則条件(regularity condition)
• 特異点があるとき以下で解ける
– 緩和法
– 適合点へのねらい撃ち法
• 特異点の場所が未知の場合はどうするか？
• 詳しい方法は省略するが，未知の内点xsにおいて制
約条件N=0,D=0を加えれば良い
9
東京大学
The University of Tokyo
池内研究室
Computer Vision Laboratory
第17章 偏微分方程式
Partial differential equations
9th Physics Based Vision Seminar
2004/July/1
偏微分方程式
•
•
初期値問題(initial value problem),
Cauchy(コーシー)問題(Cauchy
problem)
–
双曲型(hyperbolic)
• 波動方程式(wave equation)
–
放物型(parabolic)
• 拡散方程式(diffusion equation)
•
境界条件
–
–
Dirichlet(ディリクレ)条件(Dirichlet
condition)：境界点での値
Neumann(ノイマン)条件(Neumann
condition)：境界点での勾配
境界値問題(boundary value problem)
–
楕円型(elliptic)
• Poisson(ポアソン)方程式(Poisson
equation)
• Laplace(ラプラス)方程式(Laplace’s
equation)
 2u  2u

0
x2 y 2
11
9th Physics Based Vision Seminar
2004/July/1
FTCS(Forward Time Centered Space, 前進時
間中心空間)表現
• 流束保存方程式(flux-conservative equation)の一
例
を解く方法を示す→一般の場合にも適用
可能
• xをΔx間隔，tをΔt間隔にとった格子を考える
• 点(tn,xj)でのuの値u(tn,xj)をujnと表す
• 差分化
– 時間：
– 空間：
• 以上から
• これは不安定でうまくいかない
12
9th Physics Based Vision Seminar
2004/July/1
Lax法
•
を使うとFTCS法は
• von Neumann(フォンノイマン)安定解析(von
Neumann stability analysis)を使うと，安定であるた
めの条件は
– Courant-Friedrichs-Lewy(クーラント・フリードリクス・レヴ
イ)の安定性基準，Courant条件(Courant condition)
– 波動方程式の速度v，時間格子サイズΔt，空間格子サイズ
Δx，で決まる
13
9th Physics Based Vision Seminar
2004/July/1
風上差分(upwind differencing)
• 安定性の条件はCourant条件と同じ
14
9th Physics Based Vision Seminar
2004/July/1
互い違い蛙跳び(staggered leapfrog)法
• 安定性の条件はCourant条件と同じ
15
9th Physics Based Vision Seminar
2004/July/1
２ステップLax-Wendroff(ラックス・ウェンドロフ)
(Two-Step Lax-Wendroff)スキーム
• Lax法＋互い違い蛙跳び法
1 n
vt n
n1 2
n
n



u

u

u

u

u
j

1
2
j

1
j
j

1
j
•
2
2x
u nj1  u nj 

vt n1 2 n1 2
u
 u j 1 2
x j 1 2

• まとめると u  u  vxt  12 u  u  2vxt u

• 安定性の基準はCourant条件
n1
j
n
j
n
j 1
n
j
n
j 1
 12 u
 u nj 
n
j

 u nj1 


vt n n 
u  u j 1 
2x j

16
9th Physics Based Vision Seminar
2004/July/1
拡散初期値問題
• 拡散方程式
• 差分方程式に換えると：
• 陽的と陰的のFTCSスキームの
平均
これはFTCSスキーム，全陽的
(fully explicit)
• 安定性基準は
•
Crank-Nicholson(クランク・ニコ
ルソン)スキーム(CrankNicholson scheme)
• 任意のΔtに対し安定
完全陰的(fully implicit), 後退時
間(backward time)
• ３重対角の連立方程式を解く
• 任意の刻み幅Δtに対し，無条件
に安定
17
9th Physics Based Vision Seminar
2004/July/1
Schrödinger(シュレーディンガー)方程式
•
• 陰的スキーム
• 無条件に安定だがユニタリ(unitary)ではない
• 物理学の問題では
と正規化されていないといけな
い
•
と書き直す
• ケーリーの形式(Cayley’s form)というのを用いると，次式を
計算すればよいことが分かる
–
–
–
–
Hはxについての差分近似で置き換える
複素数の３重対角の連立方程式を解けば良い
方法は安定でユニタリ
Crank-Nicholson法
18
9th Physics Based Vision Seminar
2004/July/1
多次元の初期値問題
• 1次元100個の格子点，2次元100×100の格子点，
多次元計算時間がかなりかかる
• 小さな格子(8×8とか)でまず試してみて，うまくいっ
たら格子サイズを大きくしたらいいだろう
• 詳細な定式化は本を見てください
19
9th Physics Based Vision Seminar
2004/July/1
緩和法
• 以下では，境界値問題として
の式を解くことを考える
• より一般の場合については本文を見てほしい
20
9th Physics Based Vision Seminar
2004/July/1
Jacobi(ヤコビ)法(Jacobi’s method)
•
• J×Jの格子の場合
– スペクトル半径(spectral radius)ρsは
• 収束率を表し，１回の反復で残差が減る割合である
• 1未満でないと収束しない
• 0と1の間の値
• 例：ρs=0.9なら反復ごとに残差が0.9倍になる
– 誤差が10-p倍に減少するのに必要な反復回数rは
• 例：J=400,p=15のときr=1200000
21
9th Physics Based Vision Seminar
2004/July/1
Gauss-Seidel(ガウス・ザイデル)法
(Gauss-Seidel method)
•
• J×Jの格子の場合
– スペクトル半径(spectral radius)ρsは
– 誤差が10-p倍に減少するのに必要な反復回数rは
• 例：J=400,p=15のときr=600000
22
9th Physics Based Vision Seminar
2004/July/1
SOR法
(successive overrelaxation, 逐次過緩和法)


•
• ω:過緩和パラメータ(overrelaxation parameter)
• 0<ω<2の場合だけ収束
1 
1
H xn,y1   H xn1, y  H xn11, y  H xn, y1  H xn,y11   xy   (1  )H xn, y
4 
4
– 0<ω<1:劣緩和(underrelaxation),Gauss-Seidel法より収束が遅い
– 1<ω<2:過緩和(overrelaxation),Gauss-Seidel法より収束が速い
• J×Jの格子の場合
– 最適なωは
• 例:J=400のときω=1.984
– スペクトル半径(spectral radius)ρSORは
– 誤差が10-p倍に減少するのに必要な反復回数rは
• 例：J=400,p=15のときr=2000
• Chebyshev(チェビシェフ)加速(Chebyshev acceleration)によりωを反復
ごとに変化させることによりさらに収束を速くする
23
9th Physics Based Vision Seminar
2004/July/1
交互方向陰的方法(ADI)
(alternating-direction implicit method)
• SORの(8log10J)/J倍の演算回数で済む
– 例:J=400のとき反復回数は104
• Poisson方程式を解くためにNumerical Recipes in C[日本語版]に載って
いるプログラムを実行するのに最低限必要な知識
–
–
–
–
–
–
–
tridag関数のプログラムにバグがあるので注意
uが求めたい値(初期値は0にでもしておく)
物体領域内では，a=-1,b=2,c=-1,d=-1,e=2,f=-1,g=ρ
物体領域外では，a=0,b=1,c=0,d=0,e=1,f=0,g=0
jmaxは格子サイズjmax×jmax
epsは収束判定に使う値(例:eps=10-5)

  21  cos 

J
  21  cos

J 1 
J


• 例:J=400のときα=6.16847e-5,β=3.999938315
• 僕の場合は，smallval=1e-5,α=smallval,β=4-smallvalにした
– N=2k,Nはln(4J/π)に近い2の累乗
• 例：J=400のときln(4J/π)は6.23なのでN=2k=23=8にした
24
9th Physics Based Vision Seminar
2004/July/1
ADI法のアルゴリズムの概要
•
1.
2.
3.
4.
5.
6.
7.
以下にADI法のアルゴリ
ズムを示すが，分かりや
すく説明するために多少
間違ったことが書かれて
いるので詳しくは本文を読
むこと
ux, y  0
 x, y  0
1
u  u  x, y  x, y 
2  r x1, y x1, y
 x, y   x, y  2rux, y
1
u  u  x, y 
ux , y 
2  r x, y1 x, y1
 x, y   x, y  2rux, y
ux , y 
• ただし，rはk=1の場合，以
下を使う
1.
2.
r1 
r2 
    2 
2
    2 
2


    2   2
2
    2   2
2
– kが1じゃないときは本文を見
てください
– rはN=2k個の数列を繰り返す
– 例：k=3のとき(N=2k=23=8)
• 1回目の反復計算ではr1
を使う
• 8回目ではr8を使い，9回
目ではr1を使う
3に戻る
25
9th Physics Based Vision Seminar
2004/July/1
Multigrid Methods
• ADI法よりも高速な手法が提案されている
– multigrid method
– FMG(full multigrid) method
26
9th Physics Based Vision Seminar
2004/July/1
不完全Choleski(コレスキー)共役勾配法(incomplete
Choleski conjugate gradient method, ICCG)
• 格子が小さいならAx=bの形に書き換えてICCG法な
どで解けばよい
• 格子が大きいなら記憶容量の関係で不可能
– 例えばJ=400の場合は，格子は400×400=160000となる
が，行列Aのサイズが，160000×160000となる．double
型が8バイトのとき，行列Aのメモリサイズは，190GBとな
り，32ビット計算機のメモリ空間には入りきらない
– Aはほとんどの問題で疎行列なので工夫すればなんとか
なることもある
27
9th Physics Based Vision Seminar
2004/July/1
高速な方法
• Fourier(フーリエ)変換法(Fourier transform
method)
• 巡回還元法(cyclic reduction method)
• FACR(Fourier Analysis and Cyclic Reduction):上
の２つを合わせた方法，詳細は省略
• いずれも，境界がx軸とy軸に平行でなければならな
い→長方形
• multigrid法のほうが高速の場合もある
28
9th Physics Based Vision Seminar
2004/July/1
Fourier変換法
• 時間空間ではなく周波数空
間で解く
• 時間空間の場合
• 周波数空間の場合
– Jはx方向のサイズ，Lはy方
向のサイズ，mはxの位置，n
はyの位置，Δはサンプリング
間隔，̂ は  をFourier変換し
たもの， ûは u をFourier変換
したもの
1.  をFourier変換して ̂ を計
算する
2. û を左下の式で計算
3. û を逆Fourier変換して u を
求める
• この方法は周期的境界条
件に対して有効
•
•
Dirichlet境界条件u=0のと
きはサイン変換を使う
Neumannの境界条件
∇u=0のときはコサイン変
換を使う
29
9th Physics Based Vision Seminar
2004/July/1
巡回還元(cyclic reduction, CR)
•
2u 2u
Poisson方程式 x2  y2  (x, y)
を差分方程式で書くと u j1,l  u j1,l  u j,l 1  u j,l 1  4u j,l  2  j,l
• これを行列の形で書くと
 u j 1,0  
 u j ,0   u j 1,0   2  j ,0 

 
 


   
        
u
 
  2  
 u   u
1

4
1
 j 1,l 1  
 j ,l 1   j 1,l 1   2 j ,l 1 
 u j 1,l    0  0 1  4 1 0  0  u j ,l    u j 1,l      j ,l 

u
 
  2
 u   u


1

4
1
j ,l 1 
 j 1,l 1  
 j ,l 1   j 1,l 1  
   
        


 
 
  2

u
u
u


 j ,L   j 1,L  
j ,L 
 j 1,L  
これを以下のように記す u j1  T  u j  u j1  g j 2
30
9th Physics Based Vision Seminar
2004/July/1
巡回還元(cyclic reduction, CR)
•
u j 1  T  u j  u j 1  g j 2
を3つ並べると
u j 2  T  u j 1  u j  g j 12
u j 1  T  u j  u j 1  g j 2
u j  T  u j 1  u j 2  g j 12
これらを変形して足し合わせると以下のような形にな
る(T(1)はTから計算でき，gj(1)はgj-1とgjとgj+1とTから
計算できる)
u j 2  T(1)  u j  u j 2  g(j1) 2
31
9th Physics Based Vision Seminar
2004/July/1
巡回還元(cyclic reduction, CR)
• さきほど求めた式をまた3つ並べると
u j 4  T(1)  u j 2  u j  g(j1)22
u j 2  T(1)  u j  u j 2  g(j1) 2
u j  T(1)  u j 2  u j 4  g(j1) 22
これらを変形して足し合わせると以下のような形にな
る
u j 4  T(2)  u j  u j 4  g(j2) 2
32
9th Physics Based Vision Seminar
2004/July/1
巡回還元(cyclic reduction, CR)
• 格子のサイズが2の累乗だとすると，最終的に中央
の行の方程式１個だけに変形される
• なお，u0,uJは既知の境界値
• ３重対角アルゴリズムでuJ/2を解くことができる
u0  T( f )  uJ / 2  uJ  g(Jf/ )22
33
9th Physics Based Vision Seminar
2004/July/1
巡回還元(cyclic reduction, CR)
• その一つ前のレベルf-1では，以下の３つの式が出る
が，これも３重対角アルゴリズムで解ける
u0  T( f 1)  uJ / 4  uJ / 2  g(Jf/ 21)2
uJ / 4  T( f 1)  uJ / 2  u3J / 4  g(Jf/ 21)2
←解ける
uJ / 2  T( f 1)  u3J / 4  uJ  g3( Jf / 14)2
←解ける
←解く必要なし
これを最初のレベルまで解くことにより全てのuを得
る
34
東京大学
The University of Tokyo
池内研究室
Computer Vision Laboratory
終了
9th Physics Based Vision Seminar
2004/July/1
次回
• ９月上旬を予定（ずいぶん間が空いてしまうが）
• 発表者２名＋宮崎
• 宮崎はウェーブレットをやろうかと思ってる
–
–
–
–
–
７月１日ニューメリカルレシピ
９月上旬ウェーブレット前半
１０月中旬ウェーブレット後半
２月変分法と有限要素法
３月外積展開？
36
9th Physics Based Vision Seminar
2004/July/1
ウェーブレットによる信号処理と画像処理
• 中野宏毅，山本鎭男，吉田
靖夫
• 1999
• 共立出版
• ISBN4-320-08549-3
37
東京大学
The University of Tokyo
池内研究室
Computer Vision Laboratory
© Daisuke Miyazaki 2004
All rights reserved.
http://www.cvl.iis.u-tokyo.ac.jp/