Numerical Recipes in C [日本語版] 初版 宮崎大輔 2003年12月4日（木） PBVセミナー 第６章 特殊関数 Special functions 6th Physics Based Vision Seminar 2003/Dec/4 ガンマ関数(gamma function) • ガンマ関数 • zが整数なら • 漸化式 3 6th Physics Based Vision Seminar 2003/Dec/4 ベータ関数(beta function) • ベータ関数 • ガンマ関数との関係 4 6th Physics Based Vision Seminar 2003/Dec/4 不完全ガンマ関数 (incomplete gamma function) • 不完全ガンマ関数 • 極限値 • P(a,x)の補関数 も不完全ガンマ関数という • 極限値 5 6th Physics Based Vision Seminar 2003/Dec/4 誤差関数 • 誤差関数 • その補関数 • 極限値と対称性 • 不完全ガンマ関数との関係 6 6th Physics Based Vision Seminar 2003/Dec/4 累積Poisson(ポアソン)確率関数 (cumulative Poisson probability function) • 累積Poisson確率関数Px(<k) (x:正, 整数k≧1) • 期待平均値をxとしたとき、Poissonのランダム事象 の起きる回数が0以上k-1以下の確率 • 極限値 • 不完全ガンマ関数との関係 7 6th Physics Based Vision Seminar 2003/Dec/4 カイ２乗確率関数 • P(χ2|ν)は、あるモデルについて、観測されるカイ２乗 の値がχ2より小さくなる確率 • これを1から引いたものがQ(χ2|ν) • 整数νは自由度の個数 • 極限値 • 不完全ガンマ関数との関係 8 6th Physics Based Vision Seminar 2003/Dec/4 不完全ベータ関数(incomplete beta function) • 極限値 • 対称性 9 6th Physics Based Vision Seminar 2003/Dec/4 Studentの分布(Student’s distribution) • Studentの分布 • 応用例：二つの分布が同じ平均を持つかどうか • 極限値 • 不完全ベータ関数との関係 10 6th Physics Based Vision Seminar 2003/Dec/4 F分布の確率関数 • 応用例：二つの標本が同じ分散を持つかどうか • 極限値 • 不完全ベータ関数との関係 11 6th Physics Based Vision Seminar 2003/Dec/4 累積２項確率(cumulative binomial probability) • ある事象が1回の試行で起きる確率をpとする • そのとき、n回の試行でk回以上それが起きる確率P を累積2項確率という • 不完全ベータ関数との関係 12 6th Physics Based Vision Seminar 2003/Dec/4 Bessel(ベッセル)関数(Bessel function) • Bessel関数 • Bessel関数 13 6th Physics Based Vision Seminar 2003/Dec/4 変形Bessel関数(modified Bessel function) • 変形Bessel関数 14 6th Physics Based Vision Seminar 2003/Dec/4 球面調和関数(spherical harmonics) • 球面調和関数Ylm(θ,φ) (-l≦m≦l) • Legendre(ルジャンドル)陪多項式(associated Legendre polynomials) Plm(x)との関係 • mが負でもこの関係式を使えばOK • 最初のいくつかのLegendre陪多項式と、対応する正規化された球面調和 関数 15 6th Physics Based Vision Seminar 2003/Dec/4 楕円積分(elliptic integral) • 以下のような積分を楕円積分という A( x) B( x) S ( x) C ( x) D( x) S ( x) dx A( x) B( x) S ( x) dx ただし、A,B,C,Dは任意の多項式、Sは３次か４次の 多項式 16 6th Physics Based Vision Seminar 2003/Dec/4 一般第２種楕円積分 (general elliptic integral of the 2nd kind) el 2( y, kc , a, b) (a bx 2 )dx y 0 (1 x 2 ) (1 x 2 )(1 kc2 x 2 ) • y≧0でa,b,kcは任意の実数 • 変数変換 x tan k 2 1 kc2 2 により tan y a (b 1) sin el 2( y, kc , a, b) 0 1 1 k sin 2 2 d 17 6th Physics Based Vision Seminar 2003/Dec/4 完全楕円積分(complete elliptic integral) • 一般第２種楕円積分 el 2( y, kc , a, b) (a bx 2 )dx y 0 (1 x 2 ) (1 x 2 )(1 kc2 x 2 ) でy→∞つまりtan-1y→π/2のときを完全楕円積分とい う • より一般的に、一般完全楕円積分は cel(kc , p, a, b) a bx2 0 (1 px ) (1 x )(1 k x ) 2 0 2 2 2 c 2 dx (a cos2 b sin 2 )d (cos2 p sin 2 ) cos2 kc2 sin 2 18 6th Physics Based Vision Seminar 2003/Dec/4 第１種Legendre楕円積分 (Legendre elliptic integral of the 1st kind) • 第１種Legendre楕円積分 d F ( , k ) F ( \ ) 0 dx x 0 1 k 2 sin 2 (1 x )(1 k x ) 2 2 c 2 dy y 0 (1 y 2 )(1 k 2 y 2 ) el 2( x, kc ,1,1) • ここで、 y sin , x tan, kc2 1 k 2 , kc cos , k sin 19 6th Physics Based Vision Seminar 2003/Dec/4 第１種完全楕円積分 (complete elliptic integral of the 1st kind) K (k ) F ( , k ) cel (kc ,1,1,1) 2 20 6th Physics Based Vision Seminar 2003/Dec/4 第２種Legendre楕円積分 (Legendre elliptic integral of the 2nd kind) • 第２種Legendre楕円積分 E ( , k ) E ( \ ) 0 x 0 1 k 2 sin 2 d 1 kc2 x 2 dx (1 x 2 ) 1 x 2 y 0 1 k 2 y 2 dy 1 y2 el 2( x, kc ,1, kc2 ) • ここで、 y sin , x tan, kc2 1 k 2 , kc cos , k sin 21 6th Physics Based Vision Seminar 2003/Dec/4 第２種完全楕円積分 (complete elliptic integral of the 2nd kind) • 第２種完全楕円積分 E (k ) E ( , k ) cel (kc ,1,1, kc2 ) 2 • 以下を用いることもある B ( , k ) 0 D( , k ) 0 cos 2 1 k 2 sin 2 sin 2 1 k sin 2 2 d el 2( x, kc ,1,0) d el 2( x, kc ,0,1) 22 6th Physics Based Vision Seminar 2003/Dec/4 第３種楕円積分(elliptic integral of the 3rd kind) • 第３種楕円積分 d 0 (1 n sin 2 ) 1 k 2 sin 2 ( , n, k ) dy y 0 (1 ny2 ) (1 y 2 )(1 k 2 y 2 ) x 1 x2 0 (1 px2 ) (1 x 2 )(1 kc2 x 2 ) dx el3( x, kc , p ) • ここで、 y sin , x tan, kc2 1 k 2 , kc cos , k sin p≡n+1 23 6th Physics Based Vision Seminar 2003/Dec/4 第３種完全楕円積分 (complete elliptic integral of the 3rd kind) • 第３種完全楕円積分 ( , n, k ) cel (kc , p,1,1) 2 • p≡n+1 24 6th Physics Based Vision Seminar 2003/Dec/4 Jacobi(ヤコビ)の楕円関数 (Jacobian elliptic function) • 第１種Legendre楕円積分（ここで、y=sinφ) d F ( , k ) F ( \ ) 0 • この楕円関数 関数 dx x 0 1 k 2 sin 2 (1 x )(1 k x ) 2 2 c 2 dy y 0 (1 y 2 )(1 k 2 y 2 ) el 2( x, kc ,1,1) を考える代わりにその逆 を考える。あるいは同じことであるが とも書ける • 関数cnとdnは と定義 • 以上、sn(u,kc), cn(u,kc), dn(u, kc)をJacobiの楕円関数という 25 第７章 乱数 Random numbers 6th Physics Based Vision Seminar 2003/Dec/4 線形合同法(linear congruential generators) • 整数の列I1,I2,I3,…は0以上m-1以下 • 漸化式 により生成 • mは法(modulus)、a,cは正の整数で乗数(multiplier)、 増分(increment) • この漸化式はいつかはもとに戻り、その周期はm以 下 • m,a,cをうまく選べば最大の周期mが得られる • 最大周期の場合0以上m-1以下のすべての整数が どこかで生じる • 初期値(種,seed)I0は何を選んでも、そこから数列が 出発するというだけで結局は同じこと 27 6th Physics Based Vision Seminar 2003/Dec/4 「拙速」(quick and dirty)法 • プログラムに の一文を入れるだ け 28 6th Physics Based Vision Seminar 2003/Dec/4 さらに速い乱数生成法 • 整数の長さが32ビットでm=232にすればmod mが不 要になり、単に以下を計算するだけで良い • a=1664525 • c=10013904223 29 6th Physics Based Vision Seminar 2003/Dec/4 乗算合同法 • 以下のa,mの値を使った単純な乗算合同法は、そこ そこいい乱数を生成し、他の乱数生成法を評価する 際の良い基準である 30 6th Physics Based Vision Seminar 2003/Dec/4 切り混ぜ 1. iyの値をもとに表のどれかを選ぶ 2. 表の選んだ箱の中に入っている値を出力（乱数）として吐く。 さらにiyをその箱の中の値にする 3. その箱の中には新しく計算された乱数値を入れる 31 6th Physics Based Vision Seminar 2003/Dec/4 2通りの生成法の組み合わせ • ２つの乱数値を足した新しい乱数列は周期が長くな る（＝２つの周期の最小公倍数） • m1=2147483563, a1=40014, 周期m1-1 • m2=2147483399, a2=40692, 周期m2-1 • これらの組み合わせ：周期≒2.3×1018 32 6th Physics Based Vision Seminar 2003/Dec/4 減算法(subtractive method) • • • • • • ma[1…55]にあらかじめ値を入れておく 以下、mjは必ず範囲内に収まるようにする 乱数1: ma[1]-ma[32], ma[1]に乱数1を入れる 乱数2: ma[2]-ma[33], ma[2]に乱数2を入れる 乱数3: ma[3]-ma[34], ma[3]に乱数3を入れる ・・・ 33 6th Physics Based Vision Seminar 2003/Dec/4 データ暗号化規格 (Data Encryption Standard, DES) • 入力64ビット • ビット混ぜ16回 • 暗号化関数 (cipher function) g: 32ビット→32ビット 非線形 34 6th Physics Based Vision Seminar 2003/Dec/4 データ暗号化に基づく乱数列 • 右が乱数生成のためのg • 入力は右32ビット語 • C1, C2: 定数, ただし32ビットのうち0が 16個あり1が16個ある数 • 半語を反転(reverse half-words): 左右 16ビットを入れ替える • gの出力値と左32ビット語を32ビット XORした値を新しい右32ビット語とする • 古い右32ビット語は新しい左32ビット語 に入れる • この左右32ビット語の入れ替えを4回計 算した結果を乱数値として出力とする 35 6th Physics Based Vision Seminar 2003/Dec/4 相対的な実行時間と推奨 • ran0=乗算合同法, ran1=切り混ぜ, ran2=2通りの生 成法の組み合わせ, ran3=減算法, ranqd1=拙速法, ranqd2=232線形合同法, ran4=データ暗号化乱数列 • お薦めran1, 長周期が必要ならran2, ran4は良いけ ど遅い 36 6th Physics Based Vision Seminar 2003/Dec/4 GF(2) • 位相qのGalois体をGF(q)と表す • 簡単に言うと、GF(2)は、2進数で、加減乗除ができ る集合 • 体(field)では逆元-xが存在しなければならない つまり (-x)+x=x+(-x)=0が成り立つ ＋ 0 1 ?+0=0を満たす?は0、よって-0=0 0 0 1 ?+1=0を満たす?は1、よって-1=1 1 1 0 以上から-x=x × 0 1 • また、0+0=0、1+1=0からx+x=0 0 0 0 1 0 1 37 6th Physics Based Vision Seminar 2003/Dec/4 既約多項式 • それ以上因数分解できない多項式 • x2+1は実数の世界だと既約(=分解不可能) 0,1の世界だと(x+1)(x+1)=x2+1より既約でない(=分 解可能) • GF(2)上の多項式x2+x+1は既約 38 6th Physics Based Vision Seminar 2003/Dec/4 原始多項式(primitive polynomial) • • • • GF(2)でx2+x+1は既約多項式 この根をαとする：α2+α+1=0 よってα2=-1-α=1+α α0=1 α1=α α2=1+α α3=α+α2=1=α0 • α3で1に戻っているので、GF(2)にαを付加した体は、0,1,α,α2 の元だけで体をなす • このようにして得られる体をGF(22)と表す • このように、GF(p)にm次多項式の根を追加して得られる Galois体GF(pm)の元が0,α0,α1,α2,…,αpm-2で表すことができ るとき、その既約多項式を原始多項式という 39 6th Physics Based Vision Seminar 2003/Dec/4 x4+x+1は既約多項式かつ原始多項式 • 既約多項式x4+x+1の根をαとす る • α4+α+1=0 つまりα4=-α-1=α+1 • 右から分かる通り、α15で1に戻っ ている • よって、GF(24)は0とα0～α14の元 だけでなす体である • 24-2=14である • よってこの多項式は原始多項式 である • α0=1 α1=α α2=α2 α3=α3 α4=α+1 α5=α2+α α6=α3+α2 α7=α3+α+1 α8=α2+1 α9=α3+α α10=α2+α+1 α11=α3+α2+α α12=α3+α2+α+1 α13=α3+α2+1 α14=α3+1 α15=1 40 6th Physics Based Vision Seminar 2003/Dec/4 x4+x3+x2+x+1は既約多項式だが原始多項式で はない • 既約多項式x4+x3+x2+x+1 の根をαとする • α4+α3+α2+α+1=0 つまりα4=α3+α2+α+1 • 右から分かる通り、α5で1に 戻っている • これによって作られた Galois体GF(24)は0とα0～ α4の元だけでなす体である • 24-2≠4である • よってこの多項式は原始多 項式でない • α0=1 α1=α α2=α2 α3=α3 α4=α3+α2+α+1 α5=1 41 6th Physics Based Vision Seminar 2003/Dec/4 2を法とする原始多項式 • 表で(18,5,2,1,0)とあるの は のことである 42 6th Physics Based Vision Seminar 2003/Dec/4 ランダムなビットの生成 • • • • 原始多項式を使った乱数列は最大周期2n-1 原始多項式の例： ビットをa1,a2,…,anとし、新しいビットa0を求めたい a0を生成したあと、すべてのビット一つずらすので、 古いanは失われ、前回のa0はa1になる • 計算法の例： 43 6th Physics Based Vision Seminar 2003/Dec/4 モンテカルロ積分 44 6th Physics Based Vision Seminar 2003/Dec/4 変換法(transformation method) • 一様分布で、xとx+dxの間の数を生成する確率 • x→yに変換 つまり • この dx p( y ) を解くと x F ( y) p( y)の不定積分 dy よって 一様乱数 （入力） 変換された乱数 （出力） 45 6th Physics Based Vision Seminar 2003/Dec/4 正規（ガウス）乱数 • 変換法は２次元以上でもOK • xの同時確率密度p(x1,x2,…)dx1dx2… x→yに変換 xとyは同じ個数 yの同時確率密度： |∂( )/∂( )|はxのyについてのヤコビ行列式（ヤコビア ン） • 正規分布 の乱数を生成するBox-Muller(ボックス・マラー)法 (Box-Muller method)というものがある 46 6th Physics Based Vision Seminar 2003/Dec/4 棄却法(rejection method) • p(x)の下の領域にランダムな点を一様に選ぶことができれば、そのランダ ムな点のx座標が目的の乱数 • 0≦p(x)≦Aのとき、図の長方形内部に一様にランダムな点を生成し、p(x) より下なら採択、p(x)より上なら棄却する方法もあるが棄却する点があま りに多く、非効率 • p(x)に近い関数f(x)≧p(x)があり、f(x)の不定積分が計算できるなら – まずf(x)に関して変換法で乱数x0を生成 – 次に0～f(x0)の乱数を生成し、その値がp(x0)より下なら採択、p(x0)より上なら 棄却する 第１の乱数 （入力） 第２の乱数 （入力） 47 6th Physics Based Vision Seminar 2003/Dec/4 ガンマ分布 • 整数次a>0のガンマ分布は、平均1のPoisson(ポアソン)確率 過程でa回目の事象が起きるまでの待ち時間 • ガンマ分布の乱数は確率pa(x)dxでxとx+dxの間に入る • 以下のLorentz(ローレンツ)分布(Lorentzian distribution)を 使って棄却法でガンマ分布の乱数を生成すると効率的 – 不定積分の逆関数が単なるタンジェントなので計算しやすい 48 6th Physics Based Vision Seminar 2003/Dec/4 Poisson乱数 • 与えられた時間間隔xに、平均１のPoisson確率事象がm回起きる確率 • m≧0は整数なので、Poisson分布の密度関数px(m)dmは、ほとんど至る ところ0で、m≧0が整数のところだけ∞ • この関数を階段状の関数にしてしまう [m]はmを超えない最大の整数 • この関数に対してLorentz型の関数を 使った棄却法で乱数を生成する 49 6th Physics Based Vision Seminar 2003/Dec/4 2項乱数 • 事象が確率qで起きるとき、n回の試行でm回起きる 確率は次の2項分布に従う • 2項分布に従うmのとりうる値は0以上n以下の整数 である 50 第８章 ソーティング Sorting 6th Physics Based Vision Seminar 2003/Dec/4 ソーティング • N個の要素を大きさの順にソート（並べ替え）する • N>1000のとき、Nlog2Nのオーダーの演算回数で済 む以下の2手法がおすすめ – J.W.J. Williamsによるヒープソート(Heapsort) – C.A.R. Hoareによるクイックソート(Quicksort) • N<50のとき、単純挿入法(straight insertion sort)が 良いが、N2のオーダーの演算回数である • 50<N<1000なら、Shell(シェル)ソート(Shell’s method)が良く、N1.5のオーダーの演算回数である 52 6th Physics Based Vision Seminar 2003/Dec/4 単純挿入法(straight insertion sort) • (3 5 2 1 4)を(1 2 3 4 5)に並び替えたいとき – – – – – 3を()に入れて(3) 5を(3)に入れて(3 5) 2を(3 5)に入れて(2 3 5) 1を(2 3 5)に入れて(1 2 3 5) 4を(1 2 3 5)に入れて(1 2 3 4 5) 53 6th Physics Based Vision Seminar 2003/Dec/4 Shell(シェル)ソート(Shell’s method) • 16個の数n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11, n12, n13, n14, n15, n16をソートする場合を考えよう – (n1,n9), (n2,n10), (n3,n11), (n4,n12), (n5,n13), (n6,n14), (n7,n15), (n8,n16)それぞれを単純挿入法でソート – (n1,n5,n9,n13), (n2,n6,n10,n14), (n3,n7,n11,n15), (n4,n8,n12,n16)それぞれを単純挿入法でソート – (n1,n3,n5,n7,n9,n11,n13,n15), (n2,n4,n6,n8,n10,n12,n14,n16)を 単純挿入法でソート – (n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12,n13,n14,n15,n16)を単 純挿入法でソート 54 6th Physics Based Vision Seminar 2003/Dec/4 ヒープソート(Heapsort) • a1≧a2とa1≧a3のような関係を持つ木をヒープという • ヒープの生成 – 最初は会社のトップとして雇われる – 部下より無能なら部下の下に格下げされる – どんどん格下げされていき、上下関係の正しい所でストップ • ヒープソート – – – – – 会社のトップを退職させる ２人の部下のうち大きい方をトップに格上げする こうしてできた空席にまたその部下の大きい方を昇格させる これをヒープの下まで繰り返す 順番に退職した社員を並べたものがヒープソートの結果である 55 6th Physics Based Vision Seminar 2003/Dec/4 索引づけと順位づけ • 下図のように、左端の(14,8,32,7,3,15)を並び替え、右端の (3,7,8,14,15,32)にしたい • 配列の要素の移動に伴うオーバーヘッドを減らすため、ポイ ンタのみをソートしてやったのが左から２列目 • 右から２列目は順位表 56 6th Physics Based Vision Seminar 2003/Dec/4 クイックソート(Quicksort) 1. 2. 3. • 適当な数（ピボット）を選択する ピボットより小さい数を前方、大きい数を後方に移動させる 二分割された各々のデータを、それぞれクイックソートする (8 4 3 7 6 5 2 1)の例(下線はピボット) – – – – – – (8 4 3 7 6 5 2 1) (1 4 3 2 6 5 7) (8) (1 2) (3 4 6 5 7) (8) (1 2) (3 4 5 6) (7) (8) (1 2) (3 4 5) (6) (7) (8) (1 2) (3 4) (5) (6) (7) (8) 57 6th Physics Based Vision Seminar 2003/Dec/4 同値類(equivalence class)の決定 • 「要素3と要素7は同じクラスに属する」などのように、 それぞれの要素にクラス番号を割り振りたい そのまま クラスA クラスB 統合 クラスA クラスA クラスA 58 第９章 非線形方程式と非線形連立方 程式の解法 Root finding and nonlinear sets of equations 6th Physics Based Vision Seminar 2003/Dec/4 方程式 • 方程式 を解く • 陰関数定理(implicit function theorem) – N個の未知数があればN個の方程式を同時に満たす可能 性がある • 多次元では を満たすN次元の解ベクトルxを求める • fはN次元のベクトル値の関数で、その成分が、同時 に満たされるべき個々の方程式 60 6th Physics Based Vision Seminar 2003/Dec/4 根の探索で出会う種々の状況 • 中間値の定理(intermediate value theorem) – f(a)とf(b)が異符号なら区間(a,b)に根がある a. 関数が異符号である2点a,bの間に根x1が ある b. 関数の符号は必ずしも変化しないし、根 があるとも限らない c. 多数の根を持つ病的な関数 d. 関数はa,bで異符号だが、根ではなく特異 点 61 6th Physics Based Vision Seminar 2003/Dec/4 二分法(bisection method) • アルゴリズムの動作の様子： • 誤差は半分に減る • 二分法のように、繰り返しごとに誤差が定数倍になる 手法は１次収束(linear convergence)するという • 1次よりも高いベキで収束するとき超１次収束 (superlinear convergence)するという 62 6th Physics Based Vision Seminar 2003/Dec/4 割線法, 挟み撃ち法 • 割線法(secant method) – 収束の次数は1.618 – 破線は、根を囲むか否かによら ず、最新の２点を選ぶ • 挟み撃ち法(false position method, regula falsa) – 収束の次数は低いが超１次収束 のことが多い – 破線は必ず根を囲む２点を結ぶ – 割線法より収束は遅いが確実 63 6th Physics Based Vision Seminar 2003/Dec/4 van Wijngaarden(ヴァン・ワインハールデン)–Dekker (デッカー)–Brent(ブレント)法(Brent’s algorithm) • 3点を通る補間式は • y=0のときのxを求める • 挟み撃ち法や割線法が２点を直線で結ぶが、３点を通る逆２次 関数(xがyの２次関数)を使う方法を逆２次補間(inverse quadratic interpolation)という • 逆２次補間と二分法を組み合わせた手法をBrentの方法という – 逆２次補間が失敗するようなケースでは二分法に切り替える – 探索範囲を超える場合も適切な処理を行う – 二分法の確実性＋高次の方法の速さ 64 6th Physics Based Vision Seminar 2003/Dec/4 Newton(ニュートン)-Raphson(ラフソン)法 (Newton-Raphson method) • Taylor(テイラー)展開(Taylor series expansion) • ２次以上を省略し、f(x+δ)=0とすると つまり、xからδだけ動かせばfは0に近づく=根に近づく: Newton法 • 誤差は2次収束 65 6th Physics Based Vision Seminar 2003/Dec/4 多項式の根 • 根の近くでの直線、２次曲 線、３次曲線の振る舞い • 拡大してみると、３次曲線は 一つの根、２次曲線は二つ の根を持つことがわかる 66 6th Physics Based Vision Seminar 2003/Dec/4 多項式の減次 • 最初に多項式P(x)の根rを一つ見つける • 次に、因数分解してQ(x)を求める P(x)=(x-r)Q(x) • 多項式Q(x)はP(x)より次数が一つ減るので根を推定 しやすくなる 67 6th Physics Based Vision Seminar 2003/Dec/4 各種の方法 • 多項式の根を計算する手法に以下のような手法が ある（説明は省略） – – – – Muller(マラー)法(Muller’s method) Laguerre(ラゲール)法(Laguerre’s method) Jenkins-Traub method Lehmer-Shur algorithm 68 6th Physics Based Vision Seminar 2003/Dec/4 根の改良技法 • 大体の根が求まったら、それを改良してさらに精度 良く根を求めたい • 例えば、減次を使う場合など、減次をしていくにつれ、 徐々に誤差が累積していくので、求まった根は精度 が低い可能性があるからである • このように根を改良する方法として以下の２つの手 法がある（説明は省略） – Newton-Raphson法 – Bairstow法(Bairstow’s method) 69 6th Physics Based Vision Seminar 2003/Dec/4 非線形連立方程式のNewton-Raphson法 • • • • • 非線形連立方程式を解きたい 例として２次元の場合：f(x,y)=0, g(x,y)=0 fとgのゼロ等高線の交点が求めるべき解 ３次元以上の場合、ほぼ不可能 ただし、根の近くを初期値とすれば多次元版のNewtonRaphson法が適用できる 70 6th Physics Based Vision Seminar 2003/Dec/4 多次元方程式の解法と、多次元関数の最小化 • 共役勾配法など、多変数関数の極小を求める手法が数多く あるので、それを適用したくなる人もいるだろう • 問題が発生する例：下図の点Mでは極小なのに根ではない 71 (c) Daisuke Miyazaki 2003 All rights reserved. http://www.cvl.iis.u-tokyo.ac.jp/
© Copyright 2024 ExpyDoc
