1 / 20

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

Level-3 BLAS に基づく特異値分解 アルゴリズムの SMP 上での性能. 名古屋大学　計算理工学専攻 山本有作 日本応用数理学会年会 2006 年 9 月 18 日. 目次. １ . 　はじめに ２ . 　従来の特異値分解アルゴリズムとその問題点 ３ . Level-3 BLAS に基づく特異値分解アルゴリズム ４ . 　性能評価 ５ . 　まとめと今後の課題. １ . 　はじめに. 本研究で対象とする問題 実正方行列 A の特異値分解 A = U S V T A ： n × n 密行列 S ： n × n 対角行列

ira
Download Presentation

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

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Level-3 BLASに基づく特異値分解アルゴリズムのSMP上での性能 名古屋大学　計算理工学専攻 山本有作 日本応用数理学会年会 2006年9月18日

  2. 目次 １.　はじめに ２.　従来の特異値分解アルゴリズムとその問題点 ３.Level-3 BLAS に基づく特異値分解アルゴリズム ４.　性能評価 ５.　まとめと今後の課題

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

  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 } vi = V0 yi ui = U0 xi 逆変換 逆変換 Aの特異ベクトル {ui }, {vi }

  6. 実行時間（全特異ﾍﾞｸﾄﾙ） n = 5000，Xeon 2.8GHz（1～4PU） LAPACK での実行時間（秒） 各部分の演算量と実行時間 演算量 密行列 A (8/3) n3 二重対角化 O(n2) ～ O(n3) 二重対角行列の 特異値・特異ﾍﾞｸﾄﾙ計算 4mn2 逆変換 （左右 m 本ずつの 特異ベクトル） ・二重対角化が実行時間の 　大部分を占める ・速度向上率が低い Aの特異ベクトル {ui }, {vi }

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

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

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

  10. 下三角帯行列化のアルゴリズム（続き） [Step 1]K = 1からN /L–1まで以下の[Step 2] ～ [Step 6]を繰り返す。 [Step 2]A(K, K:N) を上三角行列に変形する鏡像変換 HKR= I – WKRaKR (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 – WKLaKL (WKL)Tの計算 [Step 6]　行列・ﾌﾞﾛｯｸﾍﾞｸﾄﾙ積：　QT := aKL (WKL)TA(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. 特異ベクトルの計算手法 n • 方法１： 下三角帯行列の特異ベクトルを直接計算 • 長所 • 固有ベクトルの逆変換は１段階のみ • 逆変換の演算量は 4mn2（従来法と同じ） • 短所 • 特異ベクトル計算のための実用的な手法は帯行列用逆反復法のみ • 直交化が必要であり，演算量はO(mnL2+m2n) 0 L 0 0 0 QR法 dqds法 mdLVs法 二分法 帯行列用 逆反復法 A C B 4mn2 O(mnL2+ m2n) A の特異ﾍﾞｸﾄﾙ {ui }{vi } C の特異ﾍﾞｸﾄﾙ {zi }{wi } 特異値 {σi}

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

  14. アルゴリズムの全体像 分割統治法 （LAPACK DBDSDC） • ２段階の二重対角化と２段階の逆変換 • 二重対角行列の特異値分解には分割統治法を使用 • 特徴 • 演算量が O(n3) となる部分はすべて level-3 BLAS で実行可能 • SMP 向け並列化は，基本的に並列版 level-3 BLAS の使用により実現 • 村田法は OpenMP によるパイプライン型の並列化 0 O(n2L) (8/3)n3 0 level-2 level-3 0 0 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ステップ　・・・ 　　が同時に実行可能 第2列のbulge-chasing における，右側からの 第1の直交変換で更新 される要素 第1列のbulge-chasing における，右側からの 第3の直交変換で更新 される要素 第1列による二重対角化は，今後 　より右の要素にのみ影響を及ぼす。 第1列の計算が右下まで行くのを待たずに，第2列の計算を開始できる。

  16. ４.　性能評価 • 評価環境 • 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 ～ 20000 の乱数行列の特異値分解（全特異ベクトルを計算） • Level-3 アルゴリズムにとっては一番不利な条件 • Level-3 アルゴリズムの L（半帯幅）は各 nごとに最適値を使用

  17. Xeon での実行時間 • プロセッサ数を変えたときの実行時間 • 結果 • Level-3 アルゴリズムでは PU 数に応じて実行時間が順調に減少 • 4PUでの加速率は 3～3.2 倍 • 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倍程度は加速 • メモリバンド幅が大きいためと思われる。

More Related