Level-3 BLASに基づく特異値分解 アルゴリズムのSMP上での性能

Presentation on theme: "Level-3 BLASに基づく特異値分解 アルゴリズムのSMP上での性能"— Presentation transcript:

1 Level-3 BLASに基づく特異値分解 アルゴリズムのSMP上での性能
名古屋大学　計算理工学専攻 山本有作 日本応用数理学会年会 2006年9月18日 日立製作所の山本有作です。 「～」について発表いたします。

2 目次 １. はじめに ２. 従来の特異値分解アルゴリズムとその問題点 ３. Level-3 BLAS に基づく特異値分解アルゴリズム
１.　はじめに ２.　従来の特異値分解アルゴリズムとその問題点 ３.　Level-3 BLAS に基づく特異値分解アルゴリズム ４.　性能評価 ５.　まとめと今後の課題 本発表では，はじめに，研究の背景を述べてから，スパースソルバの概要，並列化手法，そして本研究で工夫した点の一つであるRISCプロセッサ向けの最適化についてご説明します。最後に，並列計算機SR2201上での性能評価とまとめを述べます。

3 １. はじめに 本研究で対象とする問題 応用分野 実正方行列 A の特異値分解 A = US VT A ： n×n 密行列
１.　はじめに 本研究で対象とする問題 実正方行列 A の特異値分解 A = US VT A ： n×n 密行列 S ： n×n 対角行列 U，V ： n×n 直交行列 応用分野 統計計算 （主成分分析，最小2乗法） 信号処理 （独立成分分析など） 画像処理 （圧縮，ノイズ除去） 電子状態計算 ここでは特に，Aが実対称またはエルミートの密行列の場合を考える。

4 本研究の目的 共有メモリ型並列計算機（SMP）上で高性能を達成できる特異値分解ソルバを作成し，評価 背景 問題の大規模化
CPUのマルチコア化などによる SMP 環境の普及 デュアルコア Xeon Cell プロセッサ （1+8 コア）

5 ２. 従来の特異値分解アルゴリズムとその問題点
２.　従来の特異値分解アルゴリズムとその問題点 密行列 A 計算内容 計算手法 U0TAV0 = B (U0, V0: 直交行列) ハウスホルダー法 二重対角化 二重対角行列 B QR法 分割統治法 MR3アルゴリズム I-SVDアルゴリズム 二重対角行列の 特異値・特異ﾍﾞｸﾄﾙ計算 Bvi =σi xi BTxi =σi yi Bの特異値 {σi }， 　 特異ベクトル {xi }{yi } 密行列Aをまず三重対角行列Tに相似変換してからTの固有値・固有ベクトルを求めるのが最も一般的な計算法。 三重対角化には，後に述べるハウスホルダー法を使う場合がほとんど。 三重対角行列の固有値・固有ベクトルの計算には，色々なアルゴリズムがある。 vi = V0 yi ui = U0 xi 逆変換 逆変換 Aの特異ベクトル {ui }, {vi }

6 各部分の演算量と実行時間 密行列 A 演算量 実行時間（全特異ﾍﾞｸﾄﾙ） (8/3) n3 二重対角化 二重対角行列の
n = 5000，Xeon 2.8GHz（1～4PU） LAPACK での実行時間（秒） (8/3) n3 二重対角化 二重対角行列の 特異値・特異ﾍﾞｸﾄﾙ計算 O(n2) ～ O(n3) 密行列Aをまず三重対角行列Tに相似変換してからTの固有値・固有ベクトルを求めるのが最も一般的な計算法。 三重対角化には，後に述べるハウスホルダー法を使う場合がほとんど。 三重対角行列の固有値・固有ベクトルの計算には，色々なアルゴリズムがある。 4mn2 逆変換 （左右 m 本ずつの 特異ベクトル） ・二重対角化が実行時間の 　大部分を占める ・速度向上率が低い Aの特異ベクトル {ui }, {vi }

7 二重対角化の性能が出ない原因 二重対角化の演算パターン 演算パターンに関する問題点 左右からのハウスホルダー変換による行・列の消去
A(k) := (I – a w wT ) A(k) 演算は level-2 BLAS（行列ベクトル積と rank-1 更新） ただしブロック化により半分は level-3 BLAS にすることが可能 演算パターンに関する問題点 Level-2 BLAS はデータ再利用性が低い。 　　　　キャッシュの有効利用が困難であり，単体性能が出にくい。 　　　　プロセッサ間のアクセス競合により，並列性能向上も困難 非ゼロ要素 ゼロにしたい部分 A(k) 右から の変換 左から の変換 影響を受ける部分 k

8 ３. Level-3 BLAS に基づく特異値分解アルゴリズム
２段階の二重対角化アルゴリズム（Bischof et al., ’93） 密行列 A をまず帯幅 L の下三角帯行列 C に変換 次にこの帯行列を下二重対角行列 B に変換 二重対角化を２段階で行うことの利点 前半の変換は，level-3 BLAS （行列乗算）のみを使って実行可能 　　キャッシュの有効利用が可能 後半の変換は level-2 BLAS が中心だが，演算量は O(n2L) 前半部に比べてずっと小さい。 次数 n 下三角 帯行列化 村田法 の拡張 約 (8/3)n3 O(n2L) A C B 帯幅 L

9 下三角帯行列化のアルゴリズム ブロック鏡像変換によるブロック列の消去 ブロック鏡像変換 H = I – WαWT Hは直交行列
ブロックベクトル ブロック鏡像変換によるブロック列の消去 ブロック鏡像変換 H = I – WαWT Hは直交行列 与えられたブロックベクトルを上三角 行列（正確には右上三角部分のみ 非零でそれ以外が零の行列）に変形 第 K ステップでの処理 左からH を乗算 左からHKL を乗算 右からHKR を乗算 非ゼロ要素 ゼロにしたい部分 影響を受ける部分

10 下三角帯行列化のアルゴリズム（続き） 本アルゴリズムの特徴
[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におけるメモリ競合の影響を低減可能

11 本アルゴリズムの長所と短所 長所 短所 Level-3 BLAS の利用により，二重対角化の性能を向上可能
同様のアイディアに基づく三重対角化アルゴリズムでは，高い単体性能・並列性能を確認済み 短所 特異ベクトル計算のための計算量・記憶領域が増大 ２段階の逆変換あるいは帯行列の特異値分解が必要 詳しくは次のスライドで説明 二重対角化の高速化効果が大きければ，計算量増大を考慮しても全体としては高速化できると予想 特に，求める特異ベクトルが少ない場合は効果が大きいはず。

12 特異ベクトルの計算手法 方法１： 下三角帯行列の特異ベクトルを直接計算 長所 短所 A C B 固有ベクトルの逆変換は１段階のみ
逆変換の演算量は 4mn2 （従来法と同じ） 短所 特異ベクトル計算のための実用的な手法は帯行列用逆反復法のみ 直交化が必要であり，演算量は O(mnL2+m2n) n L QR法 dqds法 mdLVs法 二分法 帯行列用 逆反復法 A C B 4mn2 O(mnL2+ m2n) A の特異ﾍﾞｸﾄﾙ {ui }{vi } C の特異ﾍﾞｸﾄﾙ {zi }{wi } 特異値 {σi }

13 SMP 上での level-3 BLAS の高速性を鑑み，方法２を採用
特異ベクトルの計算手法（続き） 方法２： 二重対角行列の特異ベクトルを計算して２回逆変換 長所 二重対角行列の特異値・特異ベクトルを求める任意の手法が適用可能 短所 逆変換の演算量が 8mn2 （従来法の2倍）。ただし level-3 化可能 村田法の変換をすべて記憶するため，n2 の記憶領域が余計に必要 n L 特異値 {σi } QR法 DC法 MR3 I-SVD A C B 4mn2 4mn2 A の特異ﾍﾞｸﾄﾙ {ui }{vi } C の特異ﾍﾞｸﾄﾙ {zi }{wi } B の特異ﾍﾞｸﾄﾙ {xi }{yi } SMP 上での level-3 BLAS の高速性を鑑み，方法２を採用

14 アルゴリズムの全体像 ２段階の二重対角化と２段階の逆変換 二重対角行列の特異値分解には分割統治法を使用 特徴 A C B
演算量が O(n3) となる部分はすべて level-3 BLAS で実行可能 SMP 向け並列化は，基本的に並列版 level-3 BLAS の使用により実現 村田法は OpenMP によるパイプライン型の並列化 分割統治法 （LAPACK DBDSDC） (8/3)n3 O(n2L) level-3 level-2 A C B O(n2) ～ O(n3) level-3 4mn2 4mn2 A の特異ﾍﾞｸﾄﾙ {ui }{vi } C の特異ﾍﾞｸﾄﾙ {zi }{wi } B の特異値 {σi } 特異ﾍﾞｸﾄﾙ {xi }{yi } level-3 level-3

15 村田法の並列化 パイプライン型の並列化 第1列の二重対角化処理と第2列の二重対角化処理の並列性 一般の場合の並列性
第1列に対する bulge-chasing の第 k ステップ 第2列に対する bulge-chasing の第 k–2 ステップ 第3列に対する bulge-chasing の第 k–4 ステップ　・・・ 　　が同時に実行可能 第1列のbulge-chasing における，右側からの 第3の直交変換で更新 される要素 第2列のbulge-chasing における，右側からの 第1の直交変換で更新 される要素 第1列による二重対角化は，今後 　より右の要素にのみ影響を及ぼす。 第1列の計算が右下まで行くのを待たずに，第2列の計算を開始できる。

16 ４. 性能評価 評価環境 評価対象・条件 Xeon (2.8GHz), 1～4PU
４.　性能評価 評価環境 Xeon (2.8GHz), 1～4PU Linux + Intel Fortran ver. 8.1 BLAS:　Intel Math Kernel Library LAPACK:　Intel Math Kernel Library ピーク性能：　5.6GFLOPS/CPU 富士通 PrimePower HPC2500 (2.0GHz), 1～32PU 富士通 Fortran BLAS:　富士通並列化版 BLAS LAPACK:　富士通並列化版 LAPACK ピーク性能：　8GFLOPS/CPU 評価対象・条件 Level-3 BLAS に基づくアルゴリズムと LAPACK の性能を比較 n = 1200 ～ の乱数行列の特異値分解（全特異ベクトルを計算） Level-3 アルゴリズムにとっては一番不利な条件 Level-3 アルゴリズムの L（半帯幅）は各 n ごとに最適値を使用

17 Xeon での実行時間 プロセッサ数を変えたときの実行時間 結果 Level-3 アルゴリズムでは PU 数に応じて実行時間が順調に減少
4PU の場合は level-3 アルゴリズムが従来法より高速 n = 1200 n = 2500 n = 5000 実行時間（秒） PU数

18 HPC2500 での実行時間 プロセッサ数を変えたときの実行時間 結果 Level-3 アルゴリズムは従来法に比べて最大3.5倍高速
プロセッサ数が多いとき加速率が鈍るのは，非並列化部分（ブロック鏡像変換の作成など）の影響と思われる。 n = 5000 n = 10000 n = 20000 実行時間（秒） 3.5倍 PU数

19 両手法の実行時間の内訳 Xeon，n=5000の場合 考察 Level-3 アルゴリズムでは，どの部分の実行時間も順調に減少
逆変換1（村田法の逆変換）の占める時間が大きい。 　　　　　この部分について，さらに高速化が必要 必要な特異ベクトルの本数が少ない場合，level-3 アルゴリズムはさらに有利

20 両手法の実行時間の内訳 HPC2500，n=10,000の場合 考察 Level-3 アルゴリズムでは，どの部分の実行時間も順調に減少
従来法は，二重対角化の部分の加速が鈍い。 ただし，32PUで6倍程度は加速 メモリバンド幅が大きいためと思われる。

21 ５. まとめと今後の課題 本研究のまとめ 今後の課題 SMP 向けに，level-3 BLAS に基づく特異値分解ソルバを開発した。
５.　まとめと今後の課題 本研究のまとめ SMP 向けに，level-3 BLAS に基づく特異値分解ソルバを開発した。 Xeon と HPC2500 で評価した結果，PU 数が多い場合は従来法より高い性能が得られた。 特に，求める特異ベクトルの本数が少ない場合は効果が大きい。 今後の課題 性能の改善 より効率の良い並列化 村田法の逆変換の高速化 I-SVD，MR3 の適用 より多様なマシン上での性能評価 マルチコアプロセッサ 専用チップ （Cell，Clear Speed など） 自動チューニング手法の適用（最適な L の自動決定） 応用プログラムへの組み込み