応用数理工学特論 線形計算と ハイパフォーマンスコンピューティング 応用数理工学特論 線形計算と ハイパフォーマンスコンピューティング 2008年5月29日 計算理工学専攻　張研究室 山本有作

目次 １. 固有値計算の概要 ２. 対称行列の固有値計算のための前処理 ３. 非対称行列の固有値計算のための前処理 １.　固有値計算の概要 ２.　対称行列の固有値計算のための前処理 ・ ハウスホルダー法による三重対角化 ・ 三重対角化処理の並列化 ・ 三重対角化処理のキャッシュ向け最適化 ３.　非対称行列の固有値計算のための前処理 ・ ハウスホルダー法によるヘッセンベルグ化

１. 固有値計算の概要 固有値計算とその応用 固有値・固有ベクトル計算の流れ 各部分の演算量 １.　固有値計算の概要 固有値計算とその応用 固有値・固有ベクトル計算の流れ 各部分の演算量 三重対角行列 T に対する固有値・固有ベクトル解法 ヘッセンベルグ行列 H に対する固有値・固有ベクトル解法

１. 固有値計算の概要 対象とする問題 固有値問題の分類 標準固有値問題 Ax = λx λ： 固有値， x ： 固有ベクトル １.　固有値計算の概要 対象とする問題 標準固有値問題　Ax = λx λ： 固有値，　x ： 固有ベクトル A が n×n 行列の場合，固有値は全部で n 個存在 固有値問題の分類 A が対称（エルミート）行列　⇔　A が非対称（非エルミート）行列 A が密行列　⇔　A が疎行列 全固有値・固有ベクトルを計算　⇔　一部のみを計算

固有値問題の応用 分野 固有値 固有ベクトル 電子状態計算 分子軌道法 エネルギーレベル 電子軌道 構造解析 固有振動数 振動モード 主成分分析 各成分の寄与度 主成分 理論流体力学 （乱流のｽｹｰﾘﾝｸﾞ指数） スケーリング指数 相関関数

固有値計算の難しさ 固有値の特徴付け 固有値問題と代数方程式 固有値問題の難しさ λが A の固有値 ⇔ det(A –λI) = 0 λは n 次代数方程式の解 固有値問題と代数方程式 n 次代数方程式が解ければ，固有値λは求まる（数学的には）。 逆に，任意の n 次代数方程式は， n×n 行列の固有値問題として定式化できる。 固有値問題の難しさ 任意の n 次代数方程式を有限回の演算で解くアルゴリズムは存在しない。 よって，任意の固有値問題を有限回の演算で解くアルゴリズムは存在しない。 ここでは特に，Aが実対称またはエルミートの密行列の場合を考える。 求解には反復法が必須

密行列の固有値計算の基本的な考え方 密行列に対して直接に反復法を適用すると，演算コストが大きい。 有限回の演算により，密行列をまず中間的な行列に相似変換 相似変換により，固有値は不変 実際には，相似変換の中でも，特に直交変換を利用。 中間的な行列とは，三重対角行列あるいはヘッセンベルグ行列 中間的な行列に対して反復法を適用し，固有値・固有ベクトルを計算 直交変換を逆に施すことにより，固有ベクトルを元の行列の固有ベクトルに変換

対称行列の固有値・固有ベクトル計算の流れ 実対称行列 A 計算内容 計算手法 三重対角化 Q*AQ = T (Q: 直交行列) ハウスホルダー法 三重対角行列 T QR法 DC法 二分法・逆反復法 MR3アルゴリズム 三重対角行列の 固有値・固有ﾍﾞｸﾄﾙ計算 Tyi =λi yi Tの固有値 {λi }， 　 固有ベクトル {yi } 密行列Aをまず三重対角行列Tに相似変換してからTの固有値・固有ベクトルを求めるのが最も一般的な計算法。 三重対角化には，後に述べるハウスホルダー法を使う場合がほとんど。 三重対角行列の固有値・固有ベクトルの計算には，色々なアルゴリズムがある。 xi = Qyi 逆変換 逆変換 Aの固有ベクトル {xi }

各部分の演算量 N×N の実対称行列 A に対し，M 組の固有値・固有ベクトルを求める場合 三重対角行列の固有値・固有ベクトルの計算：　O(NM) ～ O(N3) 　　（解法および固有値の分布により大きく異なる。） 逆変換：　2N2 M 　M≪N の場合は，三重対角化の演算量が支配的 　今回はこの部分に絞って並列化・最適化の手法を述べる。

三重対角行列 T に対する固有値・固有ベクトル解法 概要 計算量 並列化 一部の固有値 ・固有ﾍﾞｸﾄﾙ のみの計算 QR法 相似変換によりTを対角行列に変換 O(N2) + O(N3) 容易 不可 DC法 （分割統治法） Tを再帰的に半分のｻｲｽﾞの行列に分割 ～ O(N3) 二分法・ 逆反復法 二分法で固有値，逆反復法で固有ﾍﾞｸﾄﾙを計算 O(NM) + O(NM2) 非自明 可 MR3 アルゴリズム 逆反復法の改良 （dqds，RRR，twisted分解などの利用） + O(NM) ・三重対角行列の固有値・固有ベクトル計算の部分をもう少し詳しく見ると，次のようなアルゴリズムがある。本報告ではこの部分は主題ではないので，簡単にまとめてある。やや誤解を招く部分もある。 ・QR法が並列化容易というのは，O(N^3)の固有ベクトル計算部分のこと。また，一部の固有ベクトルのみを計算する手法もないことはない。 ・Dhillonのアルゴリズムは７年前のDhillonの博士論文で出てきたものであり，ここ数年大変注目され，論文もたくさん出てきている。 これについては，11/28の応用数理学会研究会で簡単なレビューを行う予定。 dqds = differential qd with shift RRR = Relatively Robust Representation; T をLLt の形で表現して計算を進める手法

非対称行列の固有値・固有ベクトル計算の流れ 実非対称行列 A 計算内容 計算手法 ヘッセンベルグ化 Q*AQ = H (Q: 直交行列) ハウスホルダー法 ﾍｯｾﾝﾍﾞﾙｸﾞ行列 H QR法 DC法 ホモトピー法 複素周回積分法 ヘッセンベルグ行列の 固有値・固有ﾍﾞｸﾄﾙ計算 Hyi =λi yi Hの固有値 {λi }， 　 固有ベクトル {yi } 密行列Aをまず三重対角行列Tに相似変換してからTの固有値・固有ベクトルを求めるのが最も一般的な計算法。 三重対角化には，後に述べるハウスホルダー法を使う場合がほとんど。 三重対角行列の固有値・固有ベクトルの計算には，色々なアルゴリズムがある。 xi = Qyi 逆変換 逆変換 Aの固有ベクトル {xi }

各部分の演算量 N×N の実対称行列 A に対し，M 組の固有値・固有ベクトルを求める場合 QR法によるヘッセンベルグ行列の固有値・固有ベクトルの計算：　　10N3 程度 逆変換：　2N2 M

ﾍｯｾﾝﾍﾞﾙｸﾞ行列 H に対する固有値・固有ベクトル解法 概要 計算量 並列化 一部の固有値 ・固有ﾍﾞｸﾄﾙ のみの計算 QR法 相似変換によりHを対角行列に変換 10N3 程度 可能 不可 複素周回 積分法 ｺｰｼｰの定理を利用し 1 / det(A–λI) の極を求める O(N2M) 容易 可 ・三重対角行列の固有値・固有ベクトル計算の部分をもう少し詳しく見ると，次のようなアルゴリズムがある。本報告ではこの部分は主題ではないので，簡単にまとめてある。やや誤解を招く部分もある。 ・QR法が並列化容易というのは，O(N^3)の固有ベクトル計算部分のこと。また，一部の固有ベクトルのみを計算する手法もないことはない。 ・Dhillonのアルゴリズムは７年前のDhillonの博士論文で出てきたものであり，ここ数年大変注目され，論文もたくさん出てきている。 これについては，11/28の応用数理学会研究会で簡単なレビューを行う予定。

２. 対称行列の固有値計算のための前処理 ハウスホルダー法による三重対角化 三重対角化処理の並列化 三重対角化処理のキャッシュ向け最適化 ２.　対称行列の固有値計算のための前処理 ハウスホルダー法による三重対角化 三重対角化処理の並列化 共有メモリ向け 分散メモリ向け 三重対角化処理のキャッシュ向け最適化 はじめに，研究の背景として，有限要素法とスパースソルバの必要性についてご説明します。

ハウスホルダー法による三重対角化 鏡像変換による列の消去 第 k ステップでの処理 鏡像変換 H = I – a uut Hは対称な直交行列 与えられたベクトルの第１成分以外をゼロにする。 第 k ステップでの処理 ベクトル 左からH を乗算 左からH を乗算 右からH を乗算 k 非ゼロ要素 ゼロにしたい部分 影響を受ける部分

各ステップでの処理の詳細 鏡像変換ベクトルの作成 σ= sqrt(d td) σ= sqrt(d td) u = (d1 – sgn(d1)σ, d2, … , dN – k ) α= 2 /∥u∥2 HCH = (I – a uu t)C (I – a uu t) 　 　 = C – a uu tC – aC uu t + a2 uu tC uu t 　 　 = C – up t – pu t + (1/2) a uu tpu t + (1/2) a up tuu t　 　　( p ≡ aC u ) 　　 = C – uq t – qu t 　 　 　( q ≡ p – (1/2) a (p tu) u ) C d H による両側からの乗算 k 実際には，C の対称性を利用し， HCH は下三角部分のみ計算する。 p ≡ aC u の計算でも，C の上三角部分は下三角部分で代用する。

ハウスホルダー法のアルゴリズム [Step 1] k = 1からN –2まで以下の[Step 2] ～ [Step 8]を繰り返す。 　[Step 2]　内積 σ(k) = sqrt(d (k)t d (k)) の計算 　[Step 3]　鏡像変換ベクトル作成： 　　u(k) = (d (k)1 – sgn(d (k)1)σ(k), d (k)2, … , d (k)N – k ) 　[Step 4]　規格化定数の計算：　α(k) = 2 /∥u(k)∥2 　[Step 5]　行列ベクトル積：　p(k) =α(k)C (k)u(k) 　[Step 6]　ベクトルの内積：　β(k) =α(k) u(k)tp(k) /2 　[Step 7]　ベクトルの加算：　q(k) = p(k) –β(k)u(k) 　[Step 8]　行列のrank-2更新：　C (k) = C (k) – u(k)q(k)t – q(k)u(k)t ピボット列 ピボット行 d (k)1 C (k) d (k) k

ハウスホルダー法の特徴 演算量 データの再利用性 並列粒度 行列サイズが N のときの演算量： 約 (4/3)N3 rank-2更新の演算量： 約 (2/3)N3 データの再利用性 演算はほとんどBLAS2のため，行列データの再利用性なし キャッシュマシンでは高い性能が期待できない。 　　ブロック化など，アルゴリズム上の工夫が必要 並列粒度 BLAS2のため O(N2/p) であり，比較的高い。

共有メモリ向けの並列化 基本方針 PU間でのアクセス競合の回避 行列ベクトル積 および rank-2更新 の部分について，BLAS2レベルで並列化 並列BLASを使えば，並列化は極めて容易 PU間でのアクセス競合の回避 BLAS2はデータ再利用性がないため， 　 バスを通した共有メモリへのアクセスが 　 頻発 アクセス競合による性能低下を回避 　 するには，後に述べるキャッシュ向け 　 最適化技法が有効 キャッシュ PU0 PU1 PU2 PU3 バス メモリ

分散メモリ向けの並列化 (1-1) ブロックサイクリック列分割による並列化 （Dongarra et al., 1992） 各ステップの処理 行列を幅 b の列ブロックに分割し，ブロックを周期的にプロセッサに割り当てる。 各ステップの処理 C (k) 鏡像変換 ベクトル k (1) 鏡像変換ベクトル u(k) およびα(k) の作成 ・第 k 列を担当するPUが実行 ﾌﾞﾛｰﾄﾞｷｬｽﾄ (2) u(k) およびα(k) のブロードキャスト ・第 k 列を担当するPUが他の全PUに送信 ・二進木を使った場合，通信にかかる時間は O(N log p) 　（ベクトル u(k) の長さは平均 N/2） u(k)

分散メモリ向けの並列化 (1-2) 各ステップの処理（続き） (3) C (k) の下三角部分と u(k) との積 ・各PUは自分の行列要素を使って部分和ベクトルを計算 ・結果の集計／ブロードキャストの通信時間はO(N log p) (4) C (k) の上三角部分と u(k) との積 ・各PUは自分の行列要素を使って結果ﾍﾞｸﾄﾙの一部を計算 ・全対全通信により全PUに結果ﾍﾞｸﾄﾙの全要素を分配する 　通信時間はO(N log p) (3)と(4)の結果の集計，β(k) の計算，q(k) の計算 rank-2更新 C (k) = C (k) – u(k)q(k)t – q(k)u(k)t ・各PUで独立に計算 × ＝ 部分和の集計 ／ﾌﾞﾛｰﾄﾞｷｬｽﾄ × ＝ 全対全通信により 結果ﾍﾞｸﾄﾙを分配 以上より，ブロックサイクリック分割の通信時間は O(N log p)

分散メモリ向けの並列化 (2-1) 双方向サイクリック 分割による並列化 （Choi et al., 1995） 各ステップの処理 行列を b×b のブロックに分割し，ブロックに対して大きさ √p× √ p のプロセッサテンプレートを周期的に割り当てる。 各ステップの処理 C (k) 鏡像変換 ベクトル k (1) 鏡像変換ベクトル u(k) およびα(k) の作成 ・第 k 列を担当するPU群が（通信しながら）実行 ﾌﾞﾛｰﾄﾞｷｬｽﾄ u(k) (2) u(k) およびα(k) のブロードキャスト ・第 k 列担当のPUが，まず横方向のPUにﾌﾞﾛｰﾄﾞｷｬｽﾄ ・次に，対角ブロックを担当するPUが縦方向のPUに 　ﾌﾞﾛｰﾄﾞｷｬｽﾄ ・通信にかかる時間は O(N log p /√p) ﾌﾞﾛｰﾄﾞｷｬｽﾄ

分散メモリ向けの並列化 (2-2) 各ステップの処理（続き） (3) C (k) の下三角部分と u(k) との積 ・各PUは自分の行列要素を使って部分和ベクトルを計算 ・横方向に集計し，対角PUに集める（O(N log p /√p) ）。 (4) C (k) の上三角部分と u(k) との積 ・縦方向に集計し，対角PUに集める（O(N log p /√p) ）。 (3)と(4)の結果の集計とブロードキャスト ・対角PUにおいて，(3)と(4)の結果を集計 ・結果の部分ベクトルを，横方向のPUにブロードキャスト（O(N log p /√p) ） ・結果の部分ベクトルを，縦方向のPUにブロードキャスト（O(N log p /√p) ） × ＝ 部分和を横方向に集計し， 結果を対角PUに集める。 × ＝ 部分和を縦方向に集計し， 結果を対角PUに集める。

分散メモリ向けの並列化 (2-3) 各ステップの処理（続き） 両並列化方式の比較 (6) β(k) の計算，q(k) の計算 各PUの計算時間 各ステップの計算時間は両方式ともO(N2 / p) 　（p に反比例して減少） 通信時間 ブロックサイクリック列分割は O(N log p) 　（pとともに増加） ブロックScattered Square 分割は O(N log p /√p) 　（pとともに減少） (6) β(k) の計算，q(k) の計算 (7) rank-2更新 C (k) = C (k) – u(k)q(k)t – q(k)u(k)t ・各PUで独立に計算 PU台数が増えるほど， 双方向サイクリック分割が有利

サイクリック列分割と双方向サイクリック分割の性能比較 性能予測モデルによる推定性能 10 5 双方向分割 10 4 ｻｲｸﾘｯｸ列分割 Mflops 10 3 10 2 10 10 1 10 2 10 3 10 4 Num. of Proc. 超並列では双方向サイクリック分割が優位

SMPクラスタ上での性能評価 16-nodes 4-nodes 1-node 1832[sec] 902[sec] 447[sec] 双方向サイクリック 分割による並列化 SR8000/F1（SMPクラスタ）の１～16ノードで評価。キャッシュ向け最適化あり。 N=2000 ～ 32000 のエルミート行列での性能 16-nodes 4-nodes 1-node 1832[sec] 902[sec] 447[sec]

三重対角化処理のキャッシュ向け最適化 従来のハウスホルダー法の問題点 Dongarra のアルゴリズム Bischof のアルゴリズム 性能・精度評価 はじめに，研究の背景として，有限要素法とスパースソルバの必要性についてご説明します。

従来のハウスホルダー法の問題点 演算量 データの再利用性 行列サイズが N のときの演算量： 約 (4/3)N3 rank-2更新の演算量： 約 (2/3)N3 データの再利用性 演算はほとんどBLAS2のため，行列データの再利用性なし キャッシュマシンでは高い性能が期待できない。 　　ブロック化など，アルゴリズム上の工夫が必要

Dongarra のアルゴリズム = – – 原理（Dongarra et al., 1992） 問題点 複数本のピボット列・行を溜めておき，後でまとめて更新 全演算量（(4/3)n3）の半分を行列乗算（BLAS3）にでき，SMPで有効 LAPACK，ScaLAPACK 等で採用され，広く使われる。 問題点 演算の残り半分は，依然としてBLAS2のまま残る。 = – × – × 複数段分の更新（BLAS3使用）

Bischof のアルゴリズム 原理（Bishof et al., 1993, 1994） 三重対角化を２段階で行うことの利点 密行列 A をまず半帯幅Lの帯行列 B に変換し，次にこの帯行列を三重対角行列 T に変換する。 三重対角化を２段階で行うことの利点 帯行列への変換は，BLAS3のみを使って実行可能。 帯行列から三重対角行列への変換の演算量はO(N2L)であり， 　　前半部に比べてずっと小さい。 次数 N 半帯幅 L 帯行列化 村田法 約 (4/3)N3 O(N2L) A B T

Bischof のアルゴリズム ブロック鏡像変換によるブロック列の消去 ブロック鏡像変換 H = I – UαUt Hは直交行列 ブロックベクトル ブロック鏡像変換によるブロック列の消去 ブロック鏡像変換 H = I – UαUt Hは直交行列 与えられたブロックベクトルの 第１ブロックを上三角行列にし， それ以外をゼロにする。 第kステップでの処理 左からH を乗算 右からH を乗算 左からH を乗算 非ゼロ要素 ゼロにしたい部分 影響を受ける部分

Bischof のアルゴリズム すべてBLAS3（行列乗算） [Step 1]　K = 1からN /L–1まで以下の[Step 2] ～ [Step 6]を繰り返す。 [Step 2]　 D(K) を第１ブロックが上三角行列で第２ブロック以下が ゼロ行列であるブロックベクトルに変換するブロック鏡像変換 I – U(K)α(K)U(K)tを求める。 [Step 3]　行列・ブロックベクトル積：　P(K) = C (K)U(K)α(K) [Step 4]　ブロックベクトルの内積：　β(K) = α(K)tU(K)tP(K) / 2 [Step 5]　ブロックベクトルの加算：　Q(K) = P(K) –U(K)β(K) [Step 6]　行列のrank-2L更新：　C (K) = C (K) – U(K)Q(K)t – Q(K)U(K)t すべてBLAS3（行列乗算）

Bischof のアルゴリズムの特徴と問題点 データの再利用性 演算のうち，O(N3) の部分はすべて BLAS3 化 キャッシュサイズに応じて L を選ぶことで，データ再利用性を最大化することが可能 演算量 後半部（村田法）の演算量は L に比例して増加 全体としての実行時間が最小になるよう，L を適切に選ぶ必要がある。 逆変換も２段階となるため，逆変換の演算量が 4N2M に倍増

性能評価 Dongarra のアルゴリズム（LAPACK版）と Bischof のアルゴリズムの性能を，以下の条件で比較。 計算機環境 Xeon (2.8GHz) 富士通 PrimePower HPC2500 プロセッサ台数 1 ～ 32 （計算機により異なる） 問題サイズ n = 1500, 3000, 6000, 12000, 24000 全固有値・固有ベクトルを計算する場合と，固有値のみ計算の場合 アルゴリズム中のパラメータ 両アルゴリズムとも，実験により得られた最適値を採用 並列化 並列版 BLAS を利用

Xeon (2.8GHz) での性能 固有値のみを計算する場合は，Bischof のアルゴリズムが有利 全固有値・固有ベクトルを計算する場合 実行時間 (sec.) n = 1500 n = 3000 n = 6000 プロセッサ数 全固有値・固有ベクトルを計算する場合 実行時間 (sec.) n = 1500 n = 3000 n = 6000 プロセッサ数 全固有値のみを計算する場合 固有値のみを計算する場合は，Bischof のアルゴリズムが有利

富士通 PrimePower HPC2500 での性能 n = 6000 n = 12000 n = 24000 実行時間 (sec.) プロセッサ数 全固有値・固有ベクトルを計算する場合 n = 6000 n = 12000 n = 24000 実行時間 (sec.) プロセッサ数 全固有値のみを計算する場合 プロセッサ数が4以上ならば，全固有値・固有ベクトルを計算する場合でも Bischof のアルゴリズムが有利

プロセッサ数と必要な固有ベクトル数に応じた 最適アルゴリズム Xeon (2.8GHz) 富士通 PrimePower HPC2500 n = 1500 n = 3000 n = 6000 必要な固有ベクトルの割合 Dongarra Bischof プロセッサ数 n = 6000 n = 12000 n = 24000 必要な 固有 ベクトル の割合 Dongarra Bischof プロセッサ 数

精度評価 テスト例題 評価マシン 評価方法 フランク行列 aij = min(i, j) ２次元ラプラス方程式の５点差分により得られる行列 Harwell-Boeing Sparse Matrix Collection の行列 乱数行列 （要素が [0, 1] の一様乱数である対称行列） 評価マシン Opteron (1.6GHz) 評価方法 各解法に対し，計算した固有値λi’の相対誤差の最大値 　　max 1≦i≦N |λi’ – λi| / |λi| を評価 比較対象の固有値λiとしては，解析解あるいはオリジナルのハウスホルダー法により求めた固有値を使用

Maximum relative error フランク行列に対する計算精度 Matrix size Maximum relative error Bischof 及び Wu の方法の誤差は，Dongarra の方法と同程度

Maximum relative error ２次元５点差分の行列に対する計算精度 Matrix size Maximum relative error Bischof / Wu の方法の誤差は Dongarra の方法に比べやや大きい。 しかし，増大の程度は１桁以内に収まっている。

Harwell-Boeing Collection の行列に対する計算精度 Matrix name Maximum relative error Bischof / Wu の方法の誤差は Dongarra の方法に比べやや大きい。 しかし，増大の程度は２～３倍以内に収まっている。

乱数行列に対する計算精度 Bischof / Wu の方法の誤差は Dongarra の方法に比べやや大きい。 しかし，増大の程度は１桁以内に収まっている。

キャッシュ向け最適化に関するまとめ キャッシュ向け最適化を行った三重対角化手法として Dongarra とBischof のアルゴリズムを取り上げ，計算機環境と問題サイズを変えて性能比較を行った。 プロセッサ数が増えるに連れ，また，必要な固有ベクトルの割合が少なくなるに連れ，Bischof のアルゴリズムが有利となる傾向が見られた。 Bischof のアルゴリズムによる固有値の相対誤差は，Dongarraのアルゴリズムに比べてやや大きい。しかし，増大の程度は多くの場合２～３倍以内，最大でも１桁程度である。 はじめに，研究の背景として，有限要素法とスパースソルバの必要性についてご説明します。 より詳細なデータについては， 　　http://www.na.cse.nagoya-u.ac.jp/~yamamoto/work.html　を参照

グループ発表について 概要 スケジュール 最後の２回の授業を使って行う。 ２～３人のグループで１つの発表を行う。 講義で勉強したことをもとに，自分たちで数値実験を行った結果を発表。あるいは論文を読んでまとめて発表。 スケジュール グループ決定：　6/19 まで 発表テーマ決定：　6/26 まで 発表：　7/10，7/17 （各グループ15～20分程度）

発表テーマの例 自分の取り組んでいるシミュレーションのプログラムを，高性能化あるいは並列化した結果の報告 講義で学んだ行列計算アルゴリズムの性能評価 分散メモリ型計算機上でのLU分解の並列化 スパースソルバの性能評価 固有値計算アルゴリズムの性能評価，など 高性能計算に関する論文を読み，まとめて報告 その他