有理算術演算を用いた対称行列の特性多項式の因子探し

全文

(1)HPCS2015 2015/5/19. 2015年ハイパフォーマンスコンピューティングと計算科学シンポジウム High Performance Computing Symposium 2015. 有理算術演算を用いた対称行列の特性多項式の因子探し早稲田大学理工学術院：寒川光. H3 =. 我々は，浮動小数点演算を用いる数値計算の汎用的なプログラミング環境に，有理算術演算を加えたプログラミング環境を開発中である．浮動小数点演算は，倍精度を想定しているが，これに 4 倍精度や多倍長精度を含めることの有効性を検討している．. 有理算術演算の実装分母と分子を，桁数を拡張可能な自然数として配列に保持し，これに符号を付加することにする．. b d bc ± ad ± = a c ac b d bd r1 × r2 = × = a c ac b d bc r1 ÷ r2 = ÷ = a c ad. r1 ± r2 =. (1). a=. (3). di ri−1. (4). i=1. によって表わし，基数（radix）は r = 232 とする．符号は別に用意することにして，各桁の数 di を 32 ビットの符号なし整数の配列要素に格納する．ここで配列長は動的に拡張可能としておく． 64 ビットの算術演算が使えれば，式 (1)(2)(3) の四則演算は，筆算で行うのと同じ方法で実装でき，正確な計算が実現される．. 実対称行列の固有値の多重度の解析 Ax = λx.. (5). 係数行列 A は実対称を想定する．行列は倍精度浮動小数点演算で作成し，有理数に変換する．固有方程式 (5) は det(A − λI) = 0 が成立したとき，非自明解をもつ．行列式は次のように多項式に展開される．. ` ´ (−1)n λn − pn−1 λn−1 − pn−2 λn−2 − · · · − p0 = 0. (6). 左辺の多項式の全係数は，行列 A の要素 aij に対する乗算と加減算だけで得られるので，行列 A の係数が整数の場合は，特性多項式の係数も整数である．有理算術演算では A を基本変換を相似変換の形で用いてヘッセンベルグ行列に変換し，さらに有理（フロベニウス）標準形に変換する．行列が縮重して，多重度が高い固有値を含む場合には，ヘッセンベルグ行列の副対角項に零要素が現れ，ヘッセンベルグ行列は小行列に分割できる．したがって正確な特性多項式を，重根のない形の多項式（square free）として得られる．. 数値例正方形の断面をもつ四角柱で，表面の温度を与えられた場合の熱伝導問題を，断面領域で考える．境界条件処理後の内側のメッシュが 2 × 2 節点と 5 × 5 節点の場合を示す．. 100◦ C. ？. 13. 100◦ C. 9 5. 14. 10 6. 15. 11 7. 16. 12 8. 43 44 45 46 47 48 49 36. 42. 29. 35. 22. 28. 15. 21. 8 40◦ C. (a) 境界条件. 1. 2. 3. 4. (b) 内側 2 × 2 分割. 14 1. 2. 3. 4. 5. 6. 7. (c) 内側 5 × 5 分割. −1 4 0 −1. −1 0 4 −1. 0 1 −1 A −1 4. (7). 行列は物理形状の対称性のために縮重して，ヘッセンベルグ行列は次のようになる．. 0. 4 B −1 H=@. −2 4 −2. 0 −1 4 0. ⓒ 2015 Information Processing Society of Japan. 1 0 0 C −1 A 4. −2 4 −2. 0 −1 4. « ,. X3 =. „ 0. −1. 0 0 −2. 24 −22 12. « (9). が得られ，特性多項式が得らる．. p3 (λ) = 48 − 44λ + 12λ2 − λ3 = (2 − λ)(4 − λ)(6 − λ). 小さいほうの小行列 H1 から特性多項式 p1 (λ) = −λ + 4 が得られる．図の右のように 5 × 5 のメッシュに切ると，行列のサイズは n = 25 になる．ヘッセンベルグ行列は 5 つの小行列 H13 , H9 , H1 , H1 , H1 に分かれる．H9 から得られる特性多項式の係数は 6 桁に及ぶ*1 ．固有値の分布を表に示す． λ1 λ2 λ4 λ5 λ7 λ9 λ11 λ16 λ18 λ20 λ22 λ23 λ25 H13. λi λ3 λ6 λ8 λ10 λ12 λ17 λ19 λ21. λ13. λ14. λ15. H1. H1. H1. λ24 H9. approx. ei .535898384862 1.26794919243 2.0 2.26794919243 3.0 3.26794919243 4.0 4.73205080757 5.0 5.73205080767 6.0 6.73205080757 7.46410161514. exact√λ 2(2 −√ 3) 3− 3 2√ 4− 3 3√ 5− 3 4√ 3+ 3 5√ 4+ 3 6√ 5 + √3 2(2 + 3). 因子探しの方法. 表から，2 行目の 1.26794919243 と 8 行目の 4.73205080757 を加えれば 6 になることに気付く．そこで，因数分解の代替として，根と係数の関係公式を用いて，浮動小数点計算で得られた近似固有値を利用した因子探しプログラムを作成する．近似固有値を測定データと見做し，根と係数の関係公式を理論式と考えれば，因子探しは「逆問題」である．モニック 2 次方程式 “x2 + rx + s = 0” に対する根と係数の関係公式は “α + β = −r”，“αβ = s” なので，2 つの候補 ei と ej を近似固有値から選び，r = −(ei + ej ) と s = ei × ej を作成し，整数性を判定する．整数 r と s が取り出せれば，除算 pr (λ) ÷ (x2 + rx + s) を試み，剰余を求め，零ならこの 2 次多項式は特性多項式の因子であることが分かる．3 次式の場合は “x3 + rx2 + sx + t = 0” に対する公式 “α + β + γ = −r, αβ + βγ + γα = s, αβγ = −t” を使う．倍精度浮動小数点数は数学的には有理数である．係数行列が有理数行列の場合でも，行列要素の分母の最小公倍数（LCM）を掛ければ整数行列に変換できる．桁数の多い数値になるので，倍精度浮動小数点計算で得られた近似固有値では整数性の判定ができない．固有値計算を 4 倍精度，8 倍精度と上げてゆけば，桁数の多い特性多項式係数にも対応できるが，この方法では 8 倍精度， 16 倍精度，. . . と変数型を追加してゆく必要がある．そこで有理数計算で得られている正確な多項式を用いて，倍精度で得られた近似値を反復改良する*3 ．この方法で，倍精度浮動小数点計算で得られた対称行列を有理数変換して，正確な固有値を因数分解の形で求めてみる．使用する数学は，線形代数の授業で習う「基本変換」と高校時代からおなじみの「根と係数の関係公式」だけである．汎用プログラミング言語に正確な計算を実現する有理算術演算が組み込まれたら，数理科学科のプログラミング教育は，数学の授業内容にナイーブなプログラミング演習を実施でき，数学の理解が深まることが期待される*4 ．. まとめ. 熱伝導係数を 1 にすると丸め誤差は発生せず，係数行列は次のように正確に得られる．. 0 4 −1 A = @ −1 0. −1. 多重度 5（λ11 = · · · = λ15 ）が存在する*2 ．. 固有値問題は線形同次方程式の n 個の未知数 λ を決定する．. 100◦ C. „ 4. (2). r1 , r2 は有理数，a, b, c, d は整数とする．正の整数を， l X. 副対角項 k4 = 0 になるので，H は H3 と H1 に分かれる．大きいほうの行列から. (8). 熱伝導係数を 0.1 とすると，係数行列は有理数になるが，係数行列を倍精度浮動小数点演算で作成すると，丸めの影響を受ける．したがってこの行列を有理数変換し，分母の LCM 倍した整数行列の桁数は大きくなり，因数分解の形も変わる．このような計算を通して，有理算術演算をもつプログランミング環境に，多倍精度数を実装することがそれほど有効ではないとの結論に達した． p9 (λ) = −λ9 + 36λ8 − · · · + 356256λ2 − 293772λ + 102960 m × m のメッシュに切ると，多重度 m が現れる． *3 「反復改良は，高校の『数学 B』の教科書に載っている「2 分法」で *1 *2. *4. 可能である．浮動小数点演算では，丸め誤差の理論を行わずに数値計算の演習を行うことが難しい．. 85.

(2)