Level-3 BLASに基づく二重対角化 アルゴリズムとその性能評価 名古屋大学 計算理工学専攻 山本有作 2006年1月5日-7日 第3回計算数学研究会 目次 1. はじめに 2. 従来の二重対角化アルゴリズム 3. Level-3 BLAS に基づく二重対角化アルゴリズム 4. 性能評価 5. 特異ベクトルの計算手法 6. まとめと今後の課題 1. はじめに • 本研究で対象とする問題 – 実正方行列 A の下二重対角行列への変換 B = U0TAV0 • A : n×n 密行列 • B : n×n 下二重対角行列 • U0,V0 : n×n 直交行列 • 応用分野 – 行列 A の特異値分解に対する前処理 • 画像処理の分野では,n = 10,000 以上の大規模行列に対す る特異値分解の需要あり • 必要な特異ベクトルは数本程度の場合が多い • 前処理を数分程度で行えることが望ましい 特異値分解の流れ 密行列 A 計算内容 計算手法 二重対角化 U0TAV0 = B ハウスホルダー法 (U0, V0: 直交行列) 二重対角行列 B 二重対角行列の 特異値・特異ベクトル計算 Bvi =σi xi BTxi =σi yi Bの特異値 {σi }, 特異ベクトル {xi }{yi } vi = V0 yi 逆変換 ui = U 0 x i Aの特異ベクトル {ui }, {vi } QR法 分割統治法 MR3アルゴリズム I-SVDアルゴリズム 逆変換 各部分の演算量 • m 組の特異値・特異ベクトルを求める場合の演算量 – – – – 二重対角化: 特異値の計算: 特異ベクトルの計算: 逆変換: (8/3) n3 O(mn) ~ O(n2) O(mn) ~ O(m2n) O(mn2) – m ≪ n の場合は演算量の大部分を二重対角化が占める。 – ハウスホルダー法の高速化が必要 2. 従来の二重対角化アルゴリズム • ハウスホルダー法による二重対角化 ベクトル – 鏡像変換 H = I – a wwT による列の消去 • Hは対称な直交行列 • 与えられたベクトルの第1成分以外をゼロにする。 左からH を乗算 • 第 k ステップでの処理 0 0 0 0 0 右からHkR を乗算 0 左からHkL を乗算 k 非ゼロ要素 ゼロにしたい部分 影響を受ける部分 ハウスホルダー法のアルゴリズム [Step 1] k = 1 から n –1 まで以下の[Step 2] ~ [Step 8]を繰り返す。 [Step 2] A(k, k+1:n) を 0 にする鏡像変換 HkR = I – akR wkR (wkR)T の計 算 [Step 3] 行列ベクトル積: p := akR A(k:n, k:n) wkR [Step 4] 行列のrank-1更新: A(k:n, k:n) := A(k:n, k:n) – p(wkR)T [Step 5] A(k+2:n, k) を 0 にする鏡像変換 HkL = I – akL wkL (wkL)T の計算 [Step 6] 行列ベクトル積: qT := akL (wkL)T A(k+1:n, k:n) [Step 7] 行列のrank-1更新: A(k +1:n, k:n) := A(k +1:n, k:n) – wkLqT 0 0 A(k, k+1:n) A(k:n, k:n) A(k +2:n, k) A(k +1:n, k :n) k 0 0 ハウスホルダー法の特徴と問題点 • 特徴 – 全演算量のほとんどが,level-2 BLAS である行列ベクトル積と 行列のrank-1更新で占められる。 • 全演算量: 約 (8/3)n3 • 行列ベクトル積の演算量: 約 (4/3)n3 • rank-1更新の演算量: 約 (4/3)n3 • 問題点 – level-2 BLASのため,行列データの再利用性なし • キャッシュミス,SMPでのメモリ競合の影響大 – 最近のマイクロプロセッサ,SMP上での高性能が期待できない。 3. Level-3 BLAS に基づく二重対角化アルゴリズム • 基本的なアイディア – 密行列 A をまず帯幅 L の下三角帯行列 C に変換 – 次にこの帯行列を下二重対角行列 B に変換 次数 n 下三角 帯行列化 約 (8/3)n3 A 0 村田法 の拡張 O(n2L) 0 C 帯幅 L 0 0 B • 二重対角化を2段階で行うことの利点 – 下三角帯行列への変換は, level-3 BLAS のみを使って実行可能 – 下三角帯行列から二重対角行列への変換の演算量は O(n2L) で あり,前半部に比べてずっと小さい。 (参考)対称行列の三重対角化の場合 • Bischof のアルゴリズム(Bishof et al., 1993, 1994) – 対称密行列 A をまず半帯幅 L の対称帯行列 C に変換 – 次にこの帯行列を三重対角行列 T に変換 次数 n 半帯幅 L 帯行列化 約 (4/3)n3 A 0 村田法 O(n2L) 0 C 0 0 T • Bischof のアルゴリズムの性能・精度(山本, 2005) – 固有値の誤差は,LAPACKで使われる Dongarra のアルゴリズム に比べ,多くの場合2~3倍程度 – 速度は2~3倍高速 – 様々なマイクロプロセッサ上で,ピークの50~70%の性能を達成 下三角帯行列化のアルゴリズム ブロックベクトル ブロック鏡像変換によるブロック列の消去 – ブロック鏡像変換 H = I – WαWT • Hは直交行列 • 与えられたブロックベクトルを上三角 左からH 行列(正確には右上三角部分のみ を乗算 非零でそれ以外が零の行列)に変形 第 K ステップでの処理 0 0 0 0 右からHK を乗算 非ゼロ要素 R 0 ゼロにしたい部分 左からHKL を乗算 0 影響を受ける部分 下三角帯行列化のアルゴリズム(続き) [Step 1] K = 1からN /L–1まで以下の[Step 2] ~ [Step 6]を繰り返す。 [Step 2] A(K, K:N) を上三角行列に変形する鏡像変換 HKR = I – WKR aKR (WKR)T の計算 [Step 3] 行列・ブロックベクトル積: P := A(K:N, K:N) WKR aKR [Step 4] 行列のrank-L更新: A(K:N, K:N) := A(K:N, K:N) – P(WKR)T [Step 5] A(K+1:N, K) を上三角行列に変形する鏡像変換 HKL = I – WKL aKL (WKL)T の計算 [Step 6] 行列・ブロックベクトル積: QT := aKL (WKL)T A(K+1:N, K:N) [Step 7] 行列のrank-L更新: A(K+1:N, K:N) := A(K+1:N, K:N) – WkLQT • 本アルゴリズムの特徴 すべて level-3 BLAS(行列乗算) – 演算が level-3 BLAS 中心のため,キャッシュの有効利用が可能 – SMPにおけるメモリ競合の影響を低減可能 下三角帯行列の二重対角化 • 対称帯行列の三重対角行列への相似変換(村田法) – 第 k 列の三重対角化 • 第 k 列の副対角要素より下を 0 にする鏡像変換を左と右からかける。 – Bulge-chasing • 上記の結果,帯より下に非零要素がはみ出すので,相似変換を繰り 返して,はみ出した分を右下に移動し,右下隅から追い出す。 • 村田法のアイディアの二重対角化への適用(Lang, 1996) – 第 k 列の二重対角化 • 第 k 列の副対角要素より下を 0 にする鏡像変換を左からかける。 – Bulge-chasing • 上記の結果,上三角部分に非零要素がはみ出すので,左右から直 交変換を繰り返し掛けることにより,はみ出した分を右下に移動し, 右下隅から追い出す。 第1列の二重対角化と bulge-chasing 左側からの 直交変換で 更新された 要素 左側からの 直交変換で 更新された 要素 第2列の二重対角化と bulge-chasing 左側からの 直交変換で 更新された 要素 左側からの 直交変換で 更新された 要素 演算量は 8n2L SMP向けの並列化 • 下三角帯行列化 – 並列版の Level-3 BLAS を使えば,逐次版のプログラムをその まま SMP 上で並列化可能 – OpenMP などを用いて手動で並列化すれば,並列化効率を更に 向上可能 SMP向けの並列化(続き) • 二重対角化向け村田法 – 第1列の二重対角化処理と第2列の二重対角化処理の並列性 第1列のbulge-chasing における,右側からの 第3の直交変換で更新 される要素 第2列のbulge-chasing における,右側からの 第1の直交変換で更新 される要素 第1列による二重対角化は,今後 より右の要素にのみ影響を及ぼす。 第1列の計算が右下まで行くのを待たずに,第2列の計算を開始できる。 – 一般の場合の並列性 • 第1列に対する bulge-chasing の第 k ステップ • 第2列に対する bulge-chasing の第 k–2 ステップ • 第3列に対する bulge-chasing の第 k–4 ステップ ・・・ が同時に実行可能 4. 性能評価 • 評価環境 – Opteron (1.8GHz), 1~4CPU • Linux 上で PGI Fortran + GOTO BLAS を使用 • ピーク性能: 3.6GFLOPS/CPU – PowerPC G5 (2.0GHz), 1~2CPU • Mac OS X 上で IBM XL Fortran + GOTO BLAS を使用 • ピーク性能: 8GFLOPS/CPU • 評価対象 – n = 2500 ~ 12500 の乱数行列の二重対角化 – 下三角帯行列化部分は並列版 GOTO BLAS,村田法部分は OpenMP により並列化 • 評価方法 – 下三角帯行列化 + 村田法の合計時間を測定し,演算量を(8/3)n3 として MFLOPS 値を算出 – ブロックサイズ L は,各 n に対して最適な値を使用 Opteron (1.8GHz) 1CPU上での性能 Opteron 1CPU 3500 L=100 Performance (GFLOPS) 3000 ピークの 80% 2500 2000 Opteron 1CPU 1500 1000 500 0 2500 5000 7500 10000 12500 Matrix size n = 12500 のとき,本アルゴリズムはピークの80%の性能を達成 Opteron (1.8GHz) 4CPU上での性能 Performance (MFLOPS) 12000 Opteron 1CPU Opteron 4CPUs 10000 L=100 ピークの 67% 286s 8000 6000 3.3倍 4000 L=100 2000 0 2500 5000 7500 10000 12500 Matrix size n = 12500 のとき,本アルゴリズムは4CPUで3.3倍の加速率 ピークの67%の性能を達成 各部分の演算時間(Opteron, n=12500) 2500 Execution time (s) 2000 村田法 下三角帯行列化 1500 1000 500 0 L=25 1CPU L=50 1CPU L=100 1CPU L=25 4CPU L=50 4CPU L=100 4CPU 大規模問題に対しては,村田法の占める時間の割合は小さい。 村田法の加速率は,L=100 のとき4CPUで3.3倍 PowerPC G5 (2.0GHz) 1CPU上での性能 PowerPC G5 1CPU 8000 Performance (GFLOPS) 7000 PowerPC G5 1CPU 6000 L=100 5000 ピークの 57% 4000 3000 2000 1000 0 2500 5000 7500 10000 12500 Matrix size n = 12500 のとき,本アルゴリズムはピークの57%の性能を達成 PowerPC G5 (2.0GHz) 2CPU上での性能 8000 Performance (MFLOPS) L=100 PowerPC G5 1CPU Power PC G5 2CPUs 7000 6000 421s 5000 ピークの 42% L=100 4000 1.5倍 3000 2000 1000 0 2500 5000 7500 10000 12500 Matrix size n = 12500 のとき,本アルゴリズムは2CPUで1.5倍の加速率 ピークの42%の性能を達成 各部分の演算時間(PowerPC G5, n=12500) 1800 1600 村田法 下三角帯行列化 Execution time (s) 1400 1200 1000 800 600 400 200 0 L=25 1CPU L=50 1CPU L=100 1CPU L=25 2CPU L=50 2CPU L=100 2CPU 大規模問題に対しては,村田法の占める時間の割合は小さい。 村田法の加速率は,L=100 のとき2CPUで1.4倍 5. 特異ベクトルの計算手法 • 方法1: 二重対角行列の特異ベクトルを計算して2回逆変換 n L 0 0 0 0 A A の特異ベクトル {ui }{vi } C 2mn2 C の特異ベクトル {zi }{wi } B 2mn2 B の特異ベクトル {xi }{yi } 特異値 {σi } QR法 DC法 MR3 I-SVD • 長所 – 二重対角行列の特異値・特異ベクトルを求める任意の手法が適用可能 • 短所 – 逆変換の演算量が 4mn2 (ただし,m が小さければ影響は少ない) – 村田法の変換をすべて記憶するため,n2 の記憶領域が余計に必要 5. 特異ベクトルの計算手法 • 方法2: 下三角帯行列の特異ベクトルを直接計算 n L 0 0 0 A A の特異ベクトル {ui }{vi } C 2mn2 0 C の特異ベクトル {zi }{wi } 帯行列用 逆反復法 O(mnL2+ m2n) B 特異値 {σi } QR法 dqds法 mdLVs法 二分法 • 長所 – 村田法の逆変換のための演算量(2mn2)・記憶領域(n2)が不要 • 短所 – 帯行列用逆反復法(O(mnL2))と特異ベクトルの直交化(O(m2n))が必要 – 特異値σi は C の高精度な固有値でないため,ツイスト分解は使用不可 5. 特異ベクトルの計算手法 • 方法3: 村田法を使わず帯行列 C の特異値・特異ベクトルを直接計算 n L 0 帯行列版 mdLVs法 特異値 {σi } 0 A A の特異ベクトル {ui }{vi } C 2mn2 帯行列版 ツイスト分解 C の特異ベクトル {zi }{wi } • 長所 – 村田法の逆変換のための演算量(2mn2)・記憶領域(n2)が不要 – うまく行けば,直交化が不要なアルゴリズムを作れる可能性あり • 問題点 – そもそも帯行列版の mdLVs 法,ツイスト分解は構成できるか? – 演算量,演算の安定性は?
© Copyright 2024 ExpyDoc