200 likes | 360 Views
Level-3 BLAS に基づく特異値分解 アルゴリズムの SMP 上での性能. 名古屋大学 計算理工学専攻 山本有作 日本応用数理学会年会 2006 年 9 月 18 日. 目次. 1 . はじめに 2 . 従来の特異値分解アルゴリズムとその問題点 3 . Level-3 BLAS に基づく特異値分解アルゴリズム 4 . 性能評価 5 . まとめと今後の課題. 1 . はじめに. 本研究で対象とする問題 実正方行列 A の特異値分解 A = U S V T A : n × n 密行列 S : n × n 対角行列
E N D
Level-3 BLASに基づく特異値分解アルゴリズムのSMP上での性能 名古屋大学 計算理工学専攻 山本有作 日本応用数理学会年会 2006年9月18日
目次 1. はじめに 2. 従来の特異値分解アルゴリズムとその問題点 3.Level-3 BLAS に基づく特異値分解アルゴリズム 4. 性能評価 5. まとめと今後の課題
1. はじめに • 本研究で対象とする問題 • 実正方行列 A の特異値分解 A = US VT • A: n×n密行列 • S: n×n対角行列 • U,V: n×n直交行列 • 応用分野 • 統計計算 (主成分分析,最小2乗法) • 信号処理 (独立成分分析など) • 画像処理 (圧縮,ノイズ除去) • 電子状態計算
本研究の目的 • 共有メモリ型並列計算機(SMP)上で高性能を達成できる特異値分解ソルバを作成し,評価 • 背景 • 問題の大規模化 • CPUのマルチコア化などによる SMP 環境の普及 デュアルコア Xeon Cell プロセッサ (1+8 コア)
2. 従来の特異値分解アルゴリズムとその問題点2. 従来の特異値分解アルゴリズムとその問題点 計算内容 計算手法 密行列 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 }
実行時間(全特異ベクトル) n = 5000,Xeon 2.8GHz(1~4PU) LAPACK での実行時間(秒) 各部分の演算量と実行時間 演算量 密行列 A (8/3) n3 二重対角化 O(n2) ~ O(n3) 二重対角行列の 特異値・特異ベクトル計算 4mn2 逆変換 (左右 m 本ずつの 特異ベクトル) ・二重対角化が実行時間の 大部分を占める ・速度向上率が低い Aの特異ベクトル {ui }, {vi }
二重対角化の性能が出ない原因 • 二重対角化の演算パターン • 左右からのハウスホルダー変換による行・列の消去 • 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
3.Level-3 BLAS に基づく特異値分解アルゴリズム • 2段階の二重対角化アルゴリズム(Bischof et al., ’93) • 密行列 A をまず帯幅 L の下三角帯行列 C に変換 • 次にこの帯行列を下二重対角行列 B に変換 • 二重対角化を2段階で行うことの利点 • 前半の変換は,level-3BLAS (行列乗算)のみを使って実行可能 キャッシュの有効利用が可能 • 後半の変換は level-2 BLAS が中心だが,演算量は O(n2L) • 前半部に比べてずっと小さい。 次数 n 下三角 帯行列化 村田法 の拡張 0 0 O(n2L) 約 (8/3)n3 0 0 C A B 帯幅 L
下三角帯行列化のアルゴリズム ブロックベクトル ブロック鏡像変換によるブロック列の消去 • ブロック鏡像変換 H = I – WαWT • Hは直交行列 • 与えられたブロックベクトルを上三角 行列(正確には右上三角部分のみ 非零でそれ以外が零の行列)に変形 第 Kステップでの処理 左からH を乗算 0 0 0 0 0 0 左からHKL を乗算 右からHKR を乗算 ゼロにしたい部分 影響を受ける部分 非ゼロ要素
下三角帯行列化のアルゴリズム(続き) [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におけるメモリ競合の影響を低減可能
本アルゴリズムの長所と短所 • 長所 • Level-3 BLAS の利用により,二重対角化の性能を向上可能 • 同様のアイディアに基づく三重対角化アルゴリズムでは,高い単体性能・並列性能を確認済み • 短所 • 特異ベクトル計算のための計算量・記憶領域が増大 • 2段階の逆変換あるいは帯行列の特異値分解が必要 • 詳しくは次のスライドで説明 • 二重対角化の高速化効果が大きければ,計算量増大を考慮しても全体としては高速化できると予想 • 特に,求める特異ベクトルが少ない場合は効果が大きいはず。
特異ベクトルの計算手法 n • 方法1: 下三角帯行列の特異ベクトルを直接計算 • 長所 • 固有ベクトルの逆変換は1段階のみ • 逆変換の演算量は 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}
SMP 上での level-3 BLAS の高速性を鑑み,方法2を採用 特異ベクトルの計算手法(続き) n • 方法2: 二重対角行列の特異ベクトルを計算して2回逆変換 • 長所 • 二重対角行列の特異値・特異ベクトルを求める任意の手法が適用可能 • 短所 • 逆変換の演算量が 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 }
アルゴリズムの全体像 分割統治法 (LAPACK DBDSDC) • 2段階の二重対角化と2段階の逆変換 • 二重対角行列の特異値分解には分割統治法を使用 • 特徴 • 演算量が 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
村田法の並列化 • パイプライン型の並列化 • 第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列の計算を開始できる。
4. 性能評価 • 評価環境 • 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ごとに最適値を使用
Xeon での実行時間 • プロセッサ数を変えたときの実行時間 • 結果 • Level-3 アルゴリズムでは PU 数に応じて実行時間が順調に減少 • 4PUでの加速率は 3~3.2 倍 • 4PU の場合は level-3 アルゴリズムが従来法より高速 n = 1200 n = 2500 n = 5000 実行時間(秒) PU数
HPC2500 での実行時間 • プロセッサ数を変えたときの実行時間 • 結果 • Level-3 アルゴリズムは従来法に比べて最大3.5倍高速 • プロセッサ数が多いとき加速率が鈍るのは,非並列化部分(ブロック鏡像変換の作成など)の影響と思われる。 n = 5000 n = 10000 n = 20000 実行時間(秒) 3.5倍 PU数
両手法の実行時間の内訳 • Xeon,n=5000の場合 • 考察 • Level-3 アルゴリズムでは,どの部分の実行時間も順調に減少 • 逆変換1(村田法の逆変換)の占める時間が大きい。 この部分について,さらに高速化が必要 • 必要な特異ベクトルの本数が少ない場合,level-3 アルゴリズムはさらに有利
両手法の実行時間の内訳 • HPC2500,n=10,000の場合 • 考察 • Level-3 アルゴリズムでは,どの部分の実行時間も順調に減少 • 従来法は,二重対角化の部分の加速が鈍い。 • ただし,32PUで6倍程度は加速 • メモリバンド幅が大きいためと思われる。