有理数計算による対称行列の固有値問題における特性多項式の因子探し

全文

(1)情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2015-HPC-149 No.10 2015/6/26. 有理数計算による対称行列の固有値問題における特性多項式の因子探し寒川光 ,a). 概要：現在開発中の「有理数計算プログラミング環境」は，有理算術演算と浮動小数点演算をひとつのプログラムで使用することができる．本稿ではこの環境の例題として，対称行列の固有値の多重度解析を扱う．対称行列を有理算術演算によってフロベニウス変換して正確な特性多項式を得て，これに浮動小数点計算によって得られた近似固有値を組み合わせて当てはめることで因子多項式を探し出すプログラムを解説する．これは自前で作る数式処理システムの因数分解機能を代替するプログラムである．このプログラムに対し，有理算術演算と浮動小数点演算の特長を活かしたチューニングにより，1 桁近い高速化を実現する HPC プログラミングの新しい教材を作成したので報告する．キーワード：有理算術演算，対称固有値問題，特性多項式，ヘッセンベルグ行列，フロベニウス標準形，因数分解, 根と係数の関係公式. A searching program for divisors of characteristic polynomial in symmetric eigenvalue problem by using rational arithmetic Hikaru Samukawa ,a). Abstract: We are developing a rational arithmetic programming environment in which floatingpoint arithmetic and rational arithmetic can be used togather. In this paper, we show a numerical example which can be achieved by using both floating-point arithmetic and rational arithmetic. In a multiplicity analysis of a symmetric eigenvalue problem, an exact characteristic polynomial is obtained by transforming symmetric matrix to Frobenius form in rational arithmetic, then search its divisors by substituting approximated eigenvalues obtained in floating-point arithmetic into Vieta’s formula. This program can replace the factoring function of formula manipulation system. Roughly ten times performance gain is achieved to the original program by combining features of the rational arithmetic and floating-point arithmetic. We are aiming to learn to use the two arithmetics differently based on thier performance features in HPC programming education. Keywords: rational arithmetic, symmetric eigenvalue problem, characteristic polynomial, Hessenberg matrix, Frobenius canonical form, factorization, Vieta’s formula. 1. a). 1. はじめに早稲田大学理工学術院 Waseda University, Shinjuku-ku, Tokyo 169–8555, Japan [email protected]. ⓒ 2015 Information Processing Society of Japan. 数値計算の多くは浮動小数点演算によって行われており，計算の効率は高く，融通がきくが，信頼性. 1.

(2) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2015-HPC-149 No.10 2015/6/26. は低い．数値計算と対照的な方法が数式処理で，信. 探す．探索は指数時間アルゴリズムであり，次数が. 頼性は高いが柔軟性に欠け，計算効率も低い．我々. 高くなると計算量が爆発する．したがって，チュー. は，浮動小数点演算を用いる数値計算のプログラミ. ニングの題材には永く使用できる．ここでプログラ. ング環境に，有理算術演算を加えることで，浮動小. ミングに使用する数学の道具は，線形代数の授業で. 数点計算で生じた丸め誤差に起因する問題を解決す. 習う基本変換（elementary transformation）と，高. る目的で使用できる「有理数計算プログラミング環. 校時代からおなじみの根と係数の関係公式（Vieta’s. 境」を開発中である．有理数の分母・分子を約 30 万. formula）だけなので，大学の初年度，2 年次の学生. 桁まで動的に拡張可能な多桁の自然数を格納できる. の HPC プログラミング演習の教材としても使用でき. 配列で保持し，有理数の四則演算ごとに分母・分子. る（固体力学，量子力学などの専門性の高い力学の. を既約分数にする．演算結果は丸めないので計算誤. 知識が不要）．次章でプログラミング環境の概要について述べる．. 差は介入せず，桁溢れしないかぎり，計算は正確である．変数型 rational を C++ によるプログラミ. 3 章で対称行列の重複固有値の解析を有理数計算で. ング環境に追加するだけで，有理数型の変数と浮動. 行うアルゴリズムついて述べる．4 章で数値実験の. 小数点型の変数は共存できる*1 ．. 結果を示す．5 章でまとめる．. この環境は次の使用を目的としている [2]．. • 丸め誤差の理論を学ぶ前の学生を対象とした数. 2. 有理算術演算とプログラミング環境有理算術演算は計算機科学の黎明期からの研究テー. 学とプログラミングの教育. • 数値シミュレーションで用いられる，浮動小数. 「有理数計算プログラミング環境」でマである [3]．. 点演算による数値計算で発生した精度に関係す. は，有理数を，分母と分子を多桁数で保持し，これに. る問題を分析するツール. 符号を付加することで表現する．四則演算は次のよ. • 計算機援用証明（Computer Assisted Proof）. うに行う．. • ベンチマークの検収など，計算機システムの稼. b d bc ± ad ± = a c ac b d bd r 1 × r2 = × = a c ac b d bc r1 ÷ r2 = ÷ = a c ad r1 ± r 2 =. 働確認本稿で紹介する解析は，「有理数計算プログラミング環境」の必要機能の調査および作動確認を兼ねて，大学の学部学生で力学知識を使用しない HPC プログラミングの教材を開発することを目的とした．浮動小数点計算では解明困難な対称行列の固有値の多重度の解析を題材にする．まず，有理算術演算で特性多項式を正確に求める．浮動小数点数は数学的には有理数である．したがって，倍精度浮動小数点演算で生成された係数行列は正確に有理行列に変換できる．有理行列の全行列要素の分母の最小公倍数（Least. Common Multiple, LCM）を求めて掛けることで，有理行列は整数行列に変換できる．整数行列の特性. *1. 十進 BASIC は，数学教育を目的に，文教大学の白石教授が開発中のプログラミング言語であるが，これによって有理数計算を体験することができる [1]．十進 BASIC は 5 つの計算モードを持つ．十進 15 桁，十進 1000 桁，2 進（Intel x86 倍精度浮動小数点），2 進による複素数，有理数である．BASIC の言語規格により，数値は単一の表現に統一されるので，浮動小数点数と有理数を同時に使用することはできない．. ⓒ 2015 Information Processing Society of Japan. (3). n 桁と n 桁の整数の積は 2n 桁になるので，式 (1) の ac，bc，ad などの整数の桁数は演算のたびに倍になる．そこで有理数を構成するときに，分母，分子の最大公約数（Greatest Common Divisor, GCD）で割り既約な分数とする．開発中のプログラミング環境は，上記の a, b, c, d などの多桁の数を longint 型で扱う多桁数の階層を持つ．. z=. で，次数と固有値の数が一致する．これを利用して，因数分解された形の特性多項式の因子（divisor）を. (2). r1 , r2 は有理数，a, b, c, d は多桁数である．. 多項式は，最高次の係数が 1 か −1 の整式となるの浮動小数点演算で求めた近似固有値から逆算して，. (1). l−1 . di r i. (4). i=0. ここに di は各桁の数である．多桁の自然数 z を，基数を r = 232 として 32 ビット符号なし整数型の配列に格納し，桁数 l を整数型で保持する．零は桁数が零（l = 0）とする．配列長は 64 から始め，不足すると動的に倍，倍と拡張する可変長とした．最大長は 32,768 としている（10 進数で約 31 万桁）*2 ．多桁 *2. log10 4294967296 = 9.633 · · · である．計算可能な最大の数は，リコンパイルで変更可能である．. 2.

(3) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2015-HPC-149 No.10 2015/6/26. 数の階層（longint クラス）では，式 (1)(2)(3) の左辺. これを係数行列 A の特性方程式（characteristic equa-. の演算の本体である四則演算の関数や，GCD 関数を. tion）と呼ぶ．左辺の多項式の全係数は，行列 A の. もつ．. 要素 aij に対する乗算と加減算だけで得られるので，. 有理数階層（rational クラス）はその上に構築される．有理数は 2 つの longint 型の数と符号で保持. 有理数計算で正確に求められる．行列 A の係数が整数の場合は，特性多項式の係数も整数である．. される．式 (1)(2)(3) の有理数演算のエントリルーチ. 非線形方程式の解は一般には複素数であるが，特. ン，有理数演算を有理数型の変数を用いて z=x+y の. 性方程式 (7) は実対称行列から生成されたので解は. ように記述するためのオペレータオーバーロードの. 実数である [5]．方程式を解く段階で，数の範囲が有. 定義のほか，min，max，floor，ceil などの有理数関. 理数から代数的無理数に拡大される．. 数や，浮動小数点数を有理数に変換する関数 Rdset，. 浮動小数点計算で全固有値を求める場合は，ヤコ. 多項式の除算などのユーティリティルーチンを提供. ビ法，ギブンス法，ハウスホルダー法などが用いら. する．. れ，この場合に得られる固有値の精度は数桁から 10. その上に有理数 BLAS（Basic Linear Algebra Sub-. 数桁である．これらの方法は，面内回転や鏡映変換. programs）階層（rblas クラスとその並列版をサポー. などの無理数を扱う変換を用いるので，有理数計算. トする prblas クラス）をもち，数値線形代数（NLA:. で正確な解を求める目的では使えない．一方，ガウ. numerical linear algebra）アプリケーションの並列. ス消去法で用いる基本変換は，四則演算だけで構成. 化を容易にする [4]．本稿の例題では，ヘッセンベル. されるので，有理算術演算に適している．. グ変換を行う elmhes 関数やフロベニウス変換を行う. hesfrb 関数を含む．. 本章の 1 節で特性多項式を求めるアルゴリズムを述べる．2 節で数値例を示す．3 節で因子探しのアル. このプログラミング環境では，我々が通常数値計. ゴリズムについて述べる．. 算で使用している C++のプログラミング言語に，変数型として rational 型が追加され，rational 型同. 3.1 特性多項式を求めるアルゴリズム行列式 (6) を展開すれば特性多項式は得られるが，. 士の演算を z=x+y のように記述することができる．浮動小数点数を正確に計算するためには，倍精度数. 素朴に展開しても多重度は分からない．A をヘッセ. の場合は関数 Rdset によって rational 型の変数に. ンベルグ行列 H に変換すると，多重度が高い固有値. 変換してから有理算術演算を用いる．この機能によ. が存在する場合は，副対角項に現れる零要素からこ. り，丸め誤差の発生しない正確な四則演算が，既存. れを判定できる．零要素があればヘッセンベルグ行. の数値計算プログラミングの延長として可能になる．. 列は小行列に分けられ，その後の計算は各ヘッセン. 3. 実対称行列の固有値の多重度の解析のアルゴリズムと数値例. 3.1.1 ヘッセンベルグ変換. 固有値問題は，次の線形同次方程式の n 個の未知. ベルグ小行列に対して独立に行える．次数 n の実対称行列 A に基本変換行列 R を (n−2) 回，相似変換（similarity transformation）の形で掛ける*3 ．. 数 λ を決定する．. Ax = λx.. (5). −1 H = Rn−2 · · · R2 R1 AR1−1 R2−1 · · · Rn−2 .. (8). 本稿では，係数行列 A は実対称とする．行列 A は倍精度浮動小数点演算で作成され，有理数に変換さ. 次数 5 のヘッセンベルグ行列を示す．. れる．固有方程式 (5) は次式が成立したとき，非自 *3. 明解をもつ．. det(A − λI) = 0.. (6). H −1 A HH −1 x = λH −1 x | {z } =I. 式 (6) は次のように展開される．. (−1)n λn + pn−1 λn−1 + pn−2 λn−2 + · · · + p0 = 0 (7). ⓒ 2015 Information Processing Society of Japan. Ax = λx の両辺に H −1 を掛けると H −1 Ax = λH −1 x になるが，H −1 H = I = HH −1 なので. すなわち (H −1 AH)H −1 x = λH −1 x となる． H −1 AH = B ，H −1 x = y とおけば，By = λy の固有方程式になる．λ は同じなので，相似変換により固有値は保存され，固有ベクトルは H −1 が掛けられる．. 3.

(4) 情報処理学会研究報告 IPSJ SIG Technical Report. ⎛. h11. ⎜ ⎜ k2 ⎜ H=⎜ ⎜ ⎜ ⎝. Vol.2015-HPC-149 No.10 2015/6/26. h12. h13. h14. h15. h22. h23. h24. h25. k3. h33. h34. h35. k4. h44. h45. k5. ⎞. p1 λ + p0 として得られる．最高次の係数は “−1” で. ⎟ ⎟ ⎟ ⎟. ⎟ ⎟ ⎠. あるが n が偶数の場合は “+1” とする．次数 n の特性多項式は式 (7) の左辺である．. (9). ヘッセンベルグ行列を経由してフロベニウス行列を求める相似変換は，1 つの行列 T にまとめられ，. h55. 対角項の下の副対角項 hs+1,s を ks+1 で表す．ヘッセンベルグ行列への変換プログラムは，数値計算法のテキスト Numerical Recipes に記載されており，有理. これを用いて相似変換の形で F = T AT −1 と記述できる．多重度が高い場合は F はブロック対角行列となる．. ⎛. 数計算でも同様のプログラムで実現できるため，ア. ⎜ ⎜ T AT −1 = ⎜ ⎜ ⎝. ルゴリズムの詳細は割愛する [6]．行列が縮重していて，複数の重複固有値が存在す. ⎞. F1 ... る場合は，副対角項に零要素が現れる（証明は省略する）．kr+1 = 0 の場合は次のように書ける．. H=. Hr. Br. 0. Hn−r. ⎟ ⎟ ⎟. ⎟ ⎠. F2 .. (11). Fm 特性多項式は各フロベニウス小行列 Fi から得られる特性多項式の積になる*4 ．. ,. p(λ) =. Hr は次数 r，Hn−r は次数 (n − r) で，どちらも上. m. pi (λ). (12). i=1. ヘッセンベルグ形である．行列式 (6) は次のように. 有理数変数を導入することにより，これまでの浮動. 計算できる．. 小数点数に対する数値計算のプログラムを融通する. det(H − λI) = det(Hr − λI) det(Hn−r − λI),. ことで正確な計算を実現することが，「有理数計算プ. つまり，kr+1 = 0 が検出されれば，後続の計算は 2 つの小行列に対して独立に進められる．特性多項式を求めるために，上ヘッセンベルグ小行列 Hr と Hn−r をフロベニウス標準形 F に変換する．相似変換操作を (n − 1) 回繰り返して，次の行列 X が得られる．. ⎛. 0. ⎜ ⎜ k2 ⎜ X=⎜ ⎜ ⎜ ⎝. x0. 0. 0. 0. 0. 0. k3. 0. 0. k4. 0. ⎟ x1 ⎟ ⎟ x2 ⎟ ⎟. ⎟ x3 ⎠. k5. x4. D. 図 1 に，内部のメッシュが 2 × 2 節点と 5 × 5 節. 度勾配と熱伝導率に比例して熱が出入りする 5 点差分スキームを用いた．未知数は内部の節点にあるので，行列の行・列番号は，境界節点の番号を詰めて振られる．解析領域を m × m に分割すると，内部は. XD の形で掛けて得られる．ここに対角行列 D. は d1 = 1, d2 = k2 , d3 = k2 k3 , d4 = k2 k3 k4 , d5 =. k2 k3 k4 k5 , ⎛ 0 ⎜ ⎜ 1 ⎜ F =⎜ ⎜ ⎜ ⎝. 正方形の断面をもつ四角柱で，表面の温度を与えられた場合の熱伝導問題を，断面領域で考える [7]．. 1 つの差分要素を示したが，上下左右の節点から温. フロベニウス行列は対角行列を相似変換 F = −1. 3.2 数値例. 点の場合を示す．メッシュ分割された (b) と (c) に. ⎞. 0. ログラミング環境」の目的である．. (m − 1) × (m − 1) にメッシュ切りされ，行列のサイズは n = (m − 1) × (m − 1) になる．3 × 3 のメッシュでは内部は 2 × 2 になるので，行列の行・列番号では節点 6 が最初の行になり，節点 7 が 2 番目，節. である．. p0. ⎞. 0. 0. 0. 0. 0. 0. 1. 0. 0. 1. 0. ⎟ p1 ⎟ ⎟ p2 ⎟ ⎟, ⎟ p3 ⎠. 1. p4. 点 10 が 3 番目，節点 11 が 4 番目になる．. pi =. n−2−i. 3.2.1 整数行列. kn−j. 熱伝導係数を 1 にすると次の行列が得られる． 0 1 4 −1 −1 0 B C B −1 4 0 −1 C C A=B (13) B C 0 4 −1 A @ −1 0 −1 −1 4. xi .. j=0. (10) 次数 5 の特性多項式が −λ5 + p4 λ4 + p3 λ3 + p2 λ2 +. ⓒ 2015 Information Processing Society of Japan. *4. 各特性多項式 pi (λ) は無平方（squarefree）である．. 4.

(5) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2015-HPC-149 No.10 2015/6/26. 行列は物理形状の対称性のために縮重して，ヘッセンベルグ行列は次のようになる． 0 4 −2 0 0 B B −1 4 −1 0 H=B B −2 4 −1 @ 0 4. λi λ1. 1. λ2. C C C C A. (14). きいほうの行列から. λ12. H3 = ⎝ −1. −2. 0. 4. −1. −2. 4. ⎞ ⎠,. λ7. λ8. λ9 λ10. 4. λ6. λ5. ⎛. 0. X3 = ⎝ −1. 0 0 −2. exact λ 0. 0.62771. √ 17) √ 0.5(7 − 33). 0.99999. 1. 2.99999. 3. 0.43844. λ4. k4 = 0 になるので，H は H3 と H1 に分かれる．大. ⎛. λ3. ei 2.37e-16. λ11. 4.56155 6.39228. 0.5(5 −. √ 17) √ 0.5(7 + 33) 0.5(5 +. ⎞ H7 H2 H1 H1 H1 24 多重度は 4 で，零固有値をもつが，浮動小数点 −22 ⎠ 計算では零固有値は正確に零としては得られず 12. e1 = 2.37 · · · × 10−16 となる．. が得られ，特性多項式が得られ，因数分解できる．. 3.3 因子探しの方法 p3 (λ) = 48 − 44λ + 12λ2 − λ3 = (2 − λ)(4 − λ)(6 − λ). 小さいほうの小行列 H1 から特性多項式 p1 (λ) =. −λ + 4 が得られる．特性多項式が定義式 (6) から λ4 − 16λ3 + 92λ2 − 224λ + 192 と得られても，それが (2 − λ)(4 − λ)2 (6 − λ) と 2 乗の項を含む形の因数分解を得るまでは重根の存在をいえない．. 整数行列の特性多項式の係数は整数で得られるが，固有値は代数的無理数になる．ここでは因数分解は，高校数学の因数分解と同様に「整式を整数の範囲で因数分解」することを考える．図 1 の (c) のように 5 × 5 のメッシュに切ると，行列のサイズは n = 25 になる．ヘッセンベルグ行列は. 5 つの小行列 H13 , H9 , H1 , H1 , H1 に分かれる．H9. 3.2.2 有理数行列熱伝導係数を 1 から 0.1 にする．ただし 0.1 は倍精度浮動小数点数とするので，2 進数の丸め誤差の影. から得られる特性多項式の係数は 6 桁に及ぶ．固有値の分布を次表に示す．. 響を受ける．倍精度浮動小数点数として生成された. λi. 行列 A を有理数変換する．“0.4” は有理数では次の. λ1. 数値に変換される．. λ2. 3, 602, 879, 701, 896, 397 0.4 . 9, 007, 199, 254, 740, 992. ei .53589838. λ3. 1.2679491. λ5. λ6. 2.2679491. λ7. λ8. 3.0. 2.0. λ4. (15). 3.2679491. λ9. λ10. 4 + 2 = (numerator) × 10 の関係がある．2 進小数で ˙ 1˙ になるので，倍精度は，10 分の 1 が循環小数 0.0001. λ11. λ12. λ16. λ17. λ18. λ19. 5.0. 浮動小数点演算では丸めの影響は避けられない．行列. λ20. λ21. 5.7320508. 要素の分母の最小公倍数は，36,028,797,018,963,968. λ22. と 17 桁の数になる．したがって最小公倍数（LCM）. λ23. この分母と分子には，丸めのために (denominator) ×. 倍した行列の要素の桁数も 17 桁と大きくなるが，ヘッセンベルグ変換すると，熱伝導係数が 1 の場合と同様に，k4 = 0 が現れる．. 3.2.3 零固有値をもつ整数行列熱伝導行列の境界条件を指定しないとグラフ・ラプラシアン行列（整数行列）が生成される．図 1 の (d) に分割数 m = 4 の場合を示した．節点数は. (m − 2)2 + 4(m − 2) などになる．固有値の分布は次表のようになる．. ⓒ 2015 Information Processing Society of Japan. λ13. λ14. λ15. 4.7320508. 6.0 6.7320508. λ24. 7.4641016. λ25 H13. 4.0. H9. H1. H1. exact λ √ 2(2 − 3) √ 3− 3 2 √. 4−. 3. 3 √ 5− 3 4 √ 3+ 3 5 √. 4+. 3. 6 √ 5+ 3 √ 2(2 + 3). H1. 多重度 5（λ11 = · · · = λ15 ）が存在する．この場合は. H9 から得られる多項式は H13 から得られる多項式を割り切るので，H13 の因子は 4 次多項式を分解して得られる．11 個の固有値が有理数で，14 個の固有値が無理数である．表の第 2 列の浮動小数点計算で得られた近似値を眺めると，例えば 2 行目の 1.2679491 と 8 行目の. 4.7320508 を加えれば 5.9999999 になることに気付. 5.

(6) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2015-HPC-149 No.10 2015/6/26. く．そこで，因数分解の代替として，浮動小数点計. 有理行列の場合は，行列要素の分母の最小公倍数. 算で得られた倍精度の近似固有値を利用した因子探. （LCM）を掛けて整数行列に変換し，桁数の多い特性. しプログラムを作成する演習問題が考えられる．. 多項式に対して，大きな固有値を用いて，因子探し. 特性多項式が因数分解されて 2 次の因子（divisor）を. を行う*6 ．この場合の整数性判定は，倍精度で得ら. 持つとする．モニック 2 次方程式 “x2 +rx+s = 0” に. れた近似値を，有理算術演算による 2 分法によって. 対する根と係数の関係公式は “α+β = −r”，“αβ = s”. 精度改良する方法をとる*7 ．有理算術演算であるか. なので，2 つの候補 ei と ej を近似固有値から選び，. ら精度改良は桁数の許す範囲で実現できる．精度改. r = −(ei + ej ) と s = ei × ej を作成し，. 良を多倍精度演算で行うのがよいか，有理算術演算で行うのがよいかを考える上で興味深い題材である．. Floor(|r|) − |r| or Ceil(|r|) − |r|. 4. 因子探しプログラム. and Floor(|s|) − |s| or Ceil(|s|) − |s|. 熱伝導係数が 1 で，分割数が 12 以下の整数行列. が許容誤差の範囲内にあるかどうかテストする*5 ．こ. の場合は，固有値が 0.1 から 10 の範囲に入り，整数性の判定を倍精度浮動小数点演算で可能なので，演. のテストを整数性判定と呼ぶことにする．範囲内で. 習問題としても難易度は高くない．この場合の処理. あれば r と s を rational 型の変数に整数として取. の概要は次のようになる．. り出し，除算 pr (λ) ÷ (x2 + rx + s) を有理算術演算. 浮動小数点計算で行列 A の固有値をヤコビ法で n. で試み，剰余が零ならこの 2 次多項式は特性多項式. 個求める．有理算術演算で A → H 変換を行う*8 ．. の因子であることが分かる．. 以下，各ヘッセンベルグ小行列 Hk に対して. 公式（Vieta’s formula）は，次数 n = 5 のモニック多項式 p0 x5 + p1 x4 + p2 x3 + p3 x2 + p4 x + p5 の場合は，根を α1 から α5 とすると次のようになる．. • 小行列 Hk をフロベニウス変換して特性多項式 p(λ) を生成する（次数 nsize）． • 既知の因子多項式がある場合はこれらで特性多項式 p(λ) を割り，割り切れれば商多項式で p(λ). p0 = −1. を置き換える（nsize を小さくする）．. p1 = α1 + α2 + α3 + α4 + α5. • 近似固有値の p(λ) に対する包含（inclusion）を成. p2 = −(α1 α2 + α1 α3 + α1 α4 + · · · + α4 α5 ). 立させる*9 ．p(λ) が 1 次式なら，p(a) · p(b) < 0，. p3 = α1 α2 α3 + α1 α2 α4 + α1 α2 α5 + · · · + α3 α4 α5. n 次式の場合，p(ai ) · p(bi ) < 0 の区間 (ai , bi ) を. p4 = −(α1 α2 α3 α4 + α1 α2 α3 α5 + · · · + α2 α3 α4 α5 ). n 個見つける． • 近似固有値 ei の整数性を調べ，整数らしけれ. p5 = α1 α2 α3 α4 α5. ば ei を有理数型の整数 e として抽出し，1 次式. n 次多項式の r 次の係数は次のようになる． pr = (−1). (r+1). n Cr . (−λ + e) で p(λ) を割り，割り切れれば因子として保存し，商多項式で p(λ) を置き換える. αi1 αi2 · · · αir. (16). • 次数 d = 2, 3, . . . について以下の 2 次以上の因子探しを行う*10 ．. i1,i2,··· ,ir. 倍精度浮動小数点数の有効桁数は 16 桁程度である．. *6. p1 = α1 + α2 で ei = αi ± εi ただし εi ≥ 0 とすると，倍精度浮動小数点計算で得られる近似値を使って，倍精度演算で調べられる係数の桁数も 16 桁を超えられない．この理由を Floor 関数の場合で記すと，. p1 ≤ ei + ej < p1 + 1 の範囲でなくては整数性は判定できないので，桁数の大きな係数に対しては，相 εi も小さくしなくてはならないからである．対誤差 ei したがって整数性判定には，求める係数の桁数に応じて近似固有値の精度を上げる必要がある． *5. 最高次の係数 p0 = 1 の多項式をモニック多項式という．. ⓒ 2015 Information Processing Society of Japan. *7 *8 *9 *10. 10 分割では，H51 , H41 と 8 つの H1 の積に分かれるが， H41 から得られる次数 41 の特性多項式の定数項は 660 桁になる． 2 分法は高校の『数学 B』にある．ヘッセンベルグ変換ルーチン elmhes や，フロベニウス変換ルーチン hesfrb は rblas クラスに用意されている．解を 1 つだけ含む区間 (a, b) を見つけることを「包含」と呼ぶ．この処理は while(d<=nsize) ループで書かれ，内側に for ループによる「組合せループ」がある．「組合せループ」は n 個の固有値のうち包含が成立した ncand 個から d 個を選ぶ ncand Cd 個の「組合せ」の反復を行うが，因子多項式が見つかると nsize が次数 d だけ小さくなる．したがって，内側のループを終える場合と，継続する候補の組合わせを探す場合に分岐させて，すでにチェック済みの組合せは再チェックしないようにすると高速化できる．こ. 6.

(7) 情報処理学会研究報告 IPSJ SIG Technical Report. – d 個の近似固有値から，根と係数の関係公式の d 1 次項 ej を求めて整数性を先行判定する． j. – 先行判定に通過すれば，式 (16) によって他の係数を求め，それらの整数性を調べる．. – 全係数が整数と見なせれば，係数を有理数型の. Vol.2015-HPC-149 No.10 2015/6/26. 式以下に分解される．. 10 分割を例に計算の進行を紹介する．ヘッセンベルグ行列 H は，H51 , H41 と 8 つの H1 に分割される．H41 の特性多項式は H1 の特性多項式で割り切れて 40 次式になり，8 つの 5 次式で割り切れる．40 C5 = 658, 008 であるが，割り切れる組合せは，. 整数として抽出し，それらによって整多項式. 66,903 番目に現れる．nsize が減り. v(λ) を作り p(λ) を割る．p(λ) = v(λ)q(λ) な. で，割り切れる組合せは，36,980 番目に現れる．. ら割り切れて，v(λ) を因子多項式とし保存し，. 30 C5. 商多項式 q(λ) で p(λ) を置き換える．. 番目に現れる．このように比較的少ない反復回数で. 以上が整数行列に対する処理であるが，熱伝導係. 35 C5. = 324, 632. = 142, 506 で，割り切れる組合せは，18,032. 因子が見つかる．. 数 0.1 の浮動小数点行列を有理数変換した行列の場. H51 の特性多項式は H41 の特性多項式で割り切れ. 合は，全分母の LCM 倍した桁数の大きい整数行列. て 10 次式になり，2 つの 5 次式で割り切れる．小行. を扱う．固有値も LCM 倍されるので，近似固有値. 列から得られる特性多項式の次数は大きく，近似固. の精度を改良する処理が加わる．改良された固有値. 有値の精度も 80 桁以上が要求されるが，5 次の因子. の有効桁数は多いので，整数性判定にも有理算術演. 多項式で割り切れるので計算時間は短い．. 算を用いる必要がある．. 4.1.2 グラフ・ラプラシアン行列表 1（下）に m を 4 から 12 までの結果を示すが，. 4.1 計算結果. m = 7 を例外として，多重度 4 は最後の 3 つの小行. 計算は，Linux ワークステーション（Intel の Xeon. 列と Ha に現れる．また Ha から得られる特性多項. 2.60 GHz）の 1 スレッドで行った（乗算に FFT は. 式の定数項は零である．熱伝導行列と比較すると，次. 使用しないで測定した）．. 数の高い因子のために計算時間が長くなる．. 4.1.1 熱伝導行列. 10 分割の場合は， H52 ，H23 ，H18 と 3 つの H1. 熱伝導係数が 1 の整数行列と 0.1 を 10 進小数と. に分割され，H18 の特性多項式が 2 つの 9 次の因子. IEEE 754 標準の倍精度浮動小数点数で与えた場合の. に分解される．H23 は分解されない．H52 の特性多. 有理行列の因子探しを計算した．結果は 3 者とも同. 項式は既知の因子のうち 1 次式，2 つの 9 次式，23. じ形に分解された．行列がヘッセンベルグ小行列の. 次式で割り切れて 10 次式になるが，1 次式（零固有. 積に分かれる状態を表 1 に示した．表では 4 − λ の 1. 値）と 4 次，5 次式に分解される．したがって計算時. 次式の因子が複数現れる場合を “3*1” のように表し. 間は 23 次式まで次数を 1 次ずつ上げて行う整数性. た．すべてのケースで λ = 4 の固有値が分割数に一. 判定と，23 次式の係数の算出に大部分を費やしてい. 致する多重度を持つ．内部分割が 5 × 5 以下は，整. る*11 ．近似固有値の精度は，23 次の特性多項式の係. 数の範囲での因数分解が 2 次以下の因子に分解され. 数の最大が約 12 桁程度であるため，固有値（零から. る．6 × 6 の場合は 3 次式を含み，7 × 7 の場合は 4. 7.77 の間）の精度改良も 30 桁程度で間にあう．しか. 次式を含む．因数分解が得られると，どの近似固有. し，因数分解されない特性多項式があるため，整数. 値がどの因子に繋がるかが分かるので，特性多項式. 性判定の組合せループの反復回数が爆発して，計算. の構成を考えることができる（固有値を数値として. 時間が長くなる．. 何 10 桁求めてみても，これは分からない）．表 1 の最下段から，9 分割以下，および 11 分割では 4 次多. 4.2 プログラムの高速化. 項式以下に分解されるので，べき根による公式を正. 有理算術演算と浮動小数点演算を使用して，数値. 確な解と考えれば，正確な固有値を得ることができ. 式に素朴にプログラミングしたオリジナルのプログ. る．10 分割は 5 次多項式以下に，12 分割は 6 次多項. ラムに対して，簡単なチューニングを施す．HPC プログラミングの教材としては，計算環境の下位の階. の分岐に，解答例には goto 文を用いた．HPC プログラミングでは，コンパイル後のアセンブラ・コードを読む機会もあるので，分岐命令を読めるようになることは，HPC では重要である．「計算機ハードウェアの理解できる唯一の言語は機械語」である [8]．. ⓒ 2015 Information Processing Society of Japan. *11. n X. n Cr. 回なので，23 次の全組合せは 8,388,584 回の反. r=2. 復を行う．指数時間アルゴリズムである．. 7.

(8) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2015-HPC-149 No.10 2015/6/26. 層（longint クラス，rational クラス，rblas クラス）. 表 1 分割数 m と小行列サイズと因子多項式の最高次数 d 熱伝導行列（上）とグラフ・ラプラシアン行列（下）. は与えられたものと考えて改良を加えずに，整数性判定，根と係数の関係公式のプログラミング，近似固有値の精度改良を対象とする．有理算術演算は，浮動小数点演算と比較するとチューニングマージンが大きいので，逐次処理の計算時間 10 分の 1 を目標としてチューニングを行うことで，有理算術演算の高速化技術と背景の数学を学ぶ．. 4.2.1 整数性判定有理数型で得られた候補固有値から 1 次の係数を. m n Ha Hb Hc Hd He d. 3 9 5 3 1. 4 16 9 3 3 1. 5 25 13 9 3*1. 6 36 19 13 4*1. 2. 2. m n Ha Hb Hc Hd He Hf d. 4 12 7 2 3*1. 5 21 12 6 3*1. 6 32 18 7 4 3*1. 2. 4. 7. 2. 3 7 43 25 15 1 1 2 1 10. 7 49 25 19 5*1. 8 64 33 25 6*1. 9 81 41 33 7*1. 10 100 51 41 8*1. 4. 3. 4. 5. 11 121 55 37 15 7 7*1 4. 8 60 33 14 10 3*1. 9 77 42 32 3*1. 10 96 52 23 18 3*1. 11 117 63 42 9 3*1. 12 140 75 34 28 3*1. 14. 18. 23. 28. 34. 12 144 73 61 10*1 6. 求めて先行判定する部分で，演算に使用する桁数を少なくする．具体的には，固有値の小数点以下だけを抽出して倍精度浮動小数点数に変換して，浮動小. 例として，近似固有値を使用した対称行列の特性多. 数点演算によって判定する．. 項式の因数分解を代替する因子探しのプログラムを. 4.2.2 根と係数の関係. 紹介した．倍精度浮動小数点計算で得られた近似値. 公式 (16) によって次数ごとに求めると乗算回数が. を，必要な精度に改善して，整数性をチェックする. 多くなるので，図 2 の順序に改めた．アルゴリズム. 単純アルゴリズムであるが，プログラミング環境が. は再帰呼び出し（recursive call）を使う*12 ．. 正確な計算をサポートすると，ちょっとした工夫で. 4.2.3 計測結果. 数式処理システムの代替も可能なケースがあること. グラフ・ラプラシアン行列で分割数が 10 と 11 の場合を示す（単位：秒）． m チューニングなし. 10 11 10 分割の. 「有理数計算プログラミング環境」の目的のひとチューニングあり. 1,152 23 C23. 94,379 と 11 分割の. 192 28 C28. 10,599 の違いは指数時. 間アルゴリズムであるため 30 倍以上あり計算時間に大きな差が現れる．チューニングの効果は 10 倍に届かなかった．浮動小数点演算との大きな違いは，桁数を削減する効果にある．2 項目のチューニングを行ったが，精度改良した固有値の分母と分子の桁数を，2 分法そのままではなく，最良近似分数に変換すると高速化に効果があると考えられる．これを今後の課題としたい．このようにいろいろな例題を通して必要機能を洗い出し，プログラミング環境が装備すべき関数を充実させている．. 5. まとめと今後の課題本論文では，既存の数値計算プログラミングの環境に有理算術演算を含めることで実現可能な計算の *12. が確認された．つに「プログラミング演習を通して行う数学教育のツール」をあげられる．プログラミング環境は目的に特化した製品によって専用化されたシステムを使用する傾向が強くなってきており，近年，学生の汎用プログラミング言語離れが顕著である．数値計算は数式処理システムや MATLAB や SciLab のようなインタプリタ形式のスクリプト言語が利用される．これではなかなか HPC プログラミングに届かない．汎用プログラミング言語を授業で教えても，日常的なツールとして使う学生はいない（学生は使わないので忘れてしまう）．一方，言語そのものが高度になったので，即戦力になるようなプログラミング技術を授業で教えるには時間が足りない．本稿で紹介した有理算術演算を四則演算のレパートリーに加えたコンパイラ型のプログラミング言語なら，数学知識の習得にプログラミングを入れられるので，数値計算プログラミングの魅力を復活させられるかもしれない．現時点で保有している数学知識を，プログラミング技術で花開かせる経験を通して，HPC プログラミング好きになる学生を増やせればと考えた．. 「有理数計算プログラミング環境」では，rational 型の変数は動的に配列要素を拡張する．プロシージャフレームに拡張されるローカル変数があると，再帰関数が機能しないので，稼働を確認する目的も含めてプログラミングした．プロシージャフレーム（またはアクティベーション・レコード）については教科書を参照されたい [8]．. ⓒ 2015 Information Processing Society of Japan. 謝辞. 本研究の一部は科学研究費補助金・基盤（C). 課題番号 25330145「有理数計算ライブラリの並列化と誤差診断ツールの開発」から支援を頂いた．記して謝意を表す．. 8.

(9) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2015-HPC-149 No.10 2015/6/26. 100◦ C. 100◦ C. 100◦ C. ？. 40◦ C. 13. 14. 15. 16. 9. 10. 11. 12. 5. 6. 7. 8. 1. 2. 3. 4. (b) 解析領域 3 × 3 分割. (a) 境界条件. 43 44 45 46 47 48 49 36 42 29 35 22 28 15 21 8 14 1 2 3 4 5 6 7 (c) 解析領域 6 × 6 分割. 11. 12. 7. 8. 9. 10. 3. 4. 5. 6. 1. 2. (d) 4 × 4 分割. 図 1 熱伝導行列とグラフ・ラプラシアン行列. p2 = −(α1 α2 + α1 α3 + α1 α4 + α1 α5 + α2 α3 + α2 α4 + α2 α5 + α3 α4 + α3 α5 + α4 α5 ) | {z } | {z } | {z } | {z } | {z } | {z } | {z } | {z } | {z } | {z } 1. 9. 13. 15. 16. 20. 22. 23. 25. 26. z }| { z }| { z }| { z }| { z }| { z }| { z }| { z }| { z }| { z }| { p3 = α1 α2 α3 + α1 α2 α4 + α1 α2 α5 + α1 α3 α4 + α1 α3 α5 + α1 α4 α5 + α2 α3 α4 + α2 α3 α5 + α2 α4 α5 + α3 α4 α5 | {z } | {z } | {z } | {z } | {z } | {z } | {z } | {z } | {z } | {z } 2. 6. 8. 10. 12. 14. 17. z }| { z }| { z }| { z }| { z }| { p4 = −(α1 α2 α3 α4 + α1 α2 α3 α5 + α1 α2 α4 α5 + α1 α3 α4 α5 + α2 α3 α4 α5 ) | {z } | {z } | {z } | {z } | {z } 3. z }| { p5 = α1 α2 α3 α4 α5 | {z }. 5. 7. 11. 19. 21. 24. 18. 4. 図 2 根と係数の関係公式の計算順序. ⓒ 2015 Information Processing Society of Japan. 9.

(10) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2015-HPC-149 No.10 2015/6/26. 参考文献 [1] [2]. [3]. [4] [5] [6]. [7] [8]. http://hp.vector.co.jp/authors/VA008683/. 寒川光, 有理数計算による実対称行列の正確な 3 重対角化による固有値の高精度計算, HPCS2013 論文集, pp. 11–22, 2013. D. Knuth., The Art of Computer Programming, Volume 2, Third Edition, Addison-Wesley, 1998, 有澤誠，和田英一監訳：The Art of Computer Programming, Third Edition ，株式会社アスキー， 2004. 寒川光, 有理数線形代数計算における有理数 blas の提案, HPCS2014 論文集, pp. 57–64, 2014. J. H. Wilkinson, Algebraic Eigenvalue Problem, Oxford University Press, 1965. W. H. Press, Teukolsky S. A., Vetterling W. T., and B. P. Flannery, Numerical Recipes in Fortran, second edition, Cambridge University Press, 1992. 寒川光, 藤野清次, 長嶋利夫, 高橋大介, HPC プログラミング—ITText シリーズ, オーム社, 2009. Patterson, D. and Hennessy, J. Computer Organization & Design: The Hardware/Software Interface, fourth edition, Morgan Kaufmann Publishers, Inc., 2011, 成田光彰訳: コンピュータの構成と設計–ハードウェアとソフトウェアのインターフェース–第 4 版,（上）（下），日経 BP 社, 2011.. ⓒ 2015 Information Processing Society of Japan. 10.

(11)