270 likes | 466 Views
長方行列向け特異値分解の 浮動小数点コプロセッサによる 高速化. 2007 年 1 月 18 日 HPCS2007. 深谷 猛,山本 有作 (名古屋大学) 畝山 多加志 (京都大学) 堀 玄,梅野 健 (理化学研究所). 概要. 背景と目的 CSX600 のアーキテクチャと性能 長方行列の特異値分解アルゴリズム CSX600 向けの高速化手法 性能評価 まとめと今後の課題. 1 . 背景と目的. 長方行列の特異値分解 応用 画像処理 電子状態計算 ( Filter Diagonalization Method )
E N D
長方行列向け特異値分解の浮動小数点コプロセッサによる高速化長方行列向け特異値分解の浮動小数点コプロセッサによる高速化 2007年1月18日 HPCS2007 深谷 猛,山本 有作 (名古屋大学) 畝山 多加志 (京都大学) 堀 玄,梅野 健 (理化学研究所)
概要 • 背景と目的 • CSX600のアーキテクチャと性能 • 長方行列の特異値分解アルゴリズム • CSX600向けの高速化手法 • 性能評価 • まとめと今後の課題
1. 背景と目的 • 長方行列の特異値分解 • 応用 • 画像処理 • 電子状態計算 (Filter Diagonalization Method) • 文書検索 (Latent Semantic Indexing) • 各種統計計算 (主成分分析,最小二乗法) A:m×n 密行列 U: m×n 列直交行列 V: n×n 直交行列 S: n×n 対角行列 = × × 数万×数千以上の行列の特異値分解 高速化が必要
高い浮動小数点能力を持つ専用プロセッサの登場高い浮動小数点能力を持つ専用プロセッサの登場 • Cell • 1個の汎用CPUコア + 8個の専用コア • 256GFLOPS (単精度) • GRAPE-DR • 512個の演算コア • 512GFLOPS(単精度) • 384GFLOPS(倍精度) • ClearSpeed CSX600 • 1個の汎用コア + 96個の演算コア • 48GFLOPS (倍精度)
専用プロセッサの特徴 • 多数の浮動小数点演算器 • 極めて高いGFLOPS値 • メモリアクセス性能は相対的に低下 • バンド幅はあまり大きくなっていない • PCIバス経由でコプロセッサとして使う場合は,特に顕著 メモリネックを軽減できるアルゴリズムが性能向上の鍵
データ量: O(N 2) 演算量: O(N 3) C C A B = + × Level-3 BLASの利用 • 行列乗算 C:=C+AB • 演算量に比べ,データ量は O(1/N) • キャッシュをうまく使うことで,メモリネックを軽減可能 行列乗算を効率的に使えれば,一般のアルゴリズムも高速化可能 (LAPACKの基本的な考え方) ※行列ベクトル積( y:= y + Ax)では,データ量・演算量ともにO(N 2)
本研究の目的 • 長方行列の特異値分解アルゴリズムをCSX600を利用して高速化する • 既存のアルゴリズムをそのまま使用せず,行列乗算をなるべく効率的に使う形で CSX600 向けに最適化する • 性能評価を行い,更なる高速化に向けての課題を明らかにする
2.CSX600のアーキテクチャと性能 • CSX600 チップ • 1個の主プロセッサ • 96個の演算専用プロセッサ • 64ビット • 倍精度2演算/サイクル • 128B レジスタファイル • 6KB SRAM • 250MHz で動作 • ピーク性能: 48GFLOPS • ClearSpeed Advance ボード • 2個の CSX600 プロセッサ • 1GB DRAM • PCI-X バスにより PC 本体と接続 • ピーク性能: 96GFLOPS
今回はこれを利用 CSX600の利用環境 • Software Development Kit • コンパイラ: Cn 言語によるチップ内並列プログラミング • デバッガ • シミュレータ • CSXL ライブラリ • ClearSpeed Advance ボード用の BLAS • 行列データを PC の主メモリ上に置いてコール • PC ⇔ ボード間の転送はライブラリ内で実行 • 公称性能: DGEMM(行列乗算)で sustained 50GFLOPS • 現状では非常に制約が多い • 使用可能なルーチンは DGEMM と DGETRF のみ • DGEMM では,M,N,K ≧ 448 でないと動かない • CSFFT ライブラリ
CSXL DGEMMの性能 (1) • サイズパラメータのうち2個が450の場合 • 性能は最大でも 11GFLOPS 程度 • 他の2個が450の場合でも,ほぼ同じ結果 N K C += A B × M 性能(MFLOPS) ・M = K = 450 ・1000 ≦ N ≦ 6000 N CSX600の性能が活かせない
CSXL DGEMMの性能 (2) • サイズパラメータのうち1個のみが450の場合 • M = N が大きい場合,性能は 40GFLOPS 以上 • 他の2個が大きい場合でも,ほぼ類似の性能 N K += A B C × 性能(MFLOPS) M ・K = 450 ・1000 ≦ M = N ≦ 6000 N 性能を引き出すには,サイズパラメータのうち少なくとも2個を 大きい値にとる必要あり
3. 長方行列の特異値分解アルゴリズム m×n密行列 A QR分解 A=QR 2mn2 n×n上三角行列 R 二重対角化 (8/3)n3 二重対角行列 B 特異値分解 O(n2) ~ O(n3) Bの特異値・特異ベクトル (B と R の特異値は同じ) 逆変換 2n3 ~4n3 R の左特異ベクトル R の右特異ベクトル = Qの乗算 4mn2 A の左特異ベクトル A の右特異ベクトル
M >> Nの場合の実行時間 • M = 40000 の場合 環境 ・ Xeon 3.2GHz ・ g77 -O3 + Intel MKL 各部分のアルゴリズム(詳細は後述) ・ QR分解: ブロック化 + 再帰的QR ・ 二重対角化: Bischof + 村田法 ・ 特異値分解: 分割統治法 実行時間(秒) QR分解とQの乗算の時間が支配的
4.CSX600向けの高速化手法 高速化の方針 • QR分解とQの乗算 • 行列乗算を効率的に利用できるようアルゴリズムを変更 • CSXLで実行 • 二重対角化および逆変換 • 二重対角化はBischof のアルゴリズムと村田法を利用 • CSXLは利用しない(CPUのみで実行) • 二重対角行列の特異値分解 • 行列乗算中心のアルゴリズムである分割統治法を利用 • LAPACKのDBDSDCを利用(CSXLが使える部分は使う)
level-2 BLAS CSXLを使えない ハウスホルダー変換によるQR分解 ハウスホルダー変換による列の消去 H1 A = ( I – t1y1y1T ) A = A(1) 上三角化とQR分解 Hn・・・ H2 H1 A = A(n) A = H1 H2・・・ Hn A(n) = QR ・・・ A A(1) A(2) A(n) = R
複数のハウスホルダー変換の合成 Compact WY representation Hn・・・ H2 H1 = ( I – tn yn ynT ) ・・・ ( I – t2y2y2T )( I – t1y1y1T ) = I – Yn Tn YnT ただし, Yn = [ y1 | y2 | ・・・ | yn] (m×n行列) Tn: n×n下三角行列 I – I – I – ・・・ = × × × × × × 複数のハウスホルダー変換の作用を,行列乗算として まとめて実行可能
更新 分解 更新 分解 分解 I – Y2T2Y2T I – Y3T3Y3T I – Y1T1Y1T H3H2H1= I – Y1T1Y1T H1 H2 H3 QR分解における行列乗算の利用法 Ⅰ 単純なブロック化(LAPACKで使用) L • ブロックのQR分解はlevel-2 BLAS • 演算量: L (ブロック幅) <<n のとき 2mn2 • CSXLを使うにはL≧450 ⇒ level-2 BLASの演算量の増加
分解 更新 分解 合成 I – YL TL YLT I – YR TR YRT I – YL TL YLT I – Y’L T’L Y’LT QR分解における行列乗算の利用法 Ⅱ 再帰的QR分解 (I – YL TL YLT ) ×(I – YR TR YRT ) = I – YTYT • ほぼ全ての演算が行列乗算 • 演算量: 3mn2
更新 分解 更新 分解 分解 I – Y2T2Y2T I – Y3T3Y3T I – Y1T1Y1T I – Y1T1Y1T QR分解における行列乗算の利用法 Ⅲ 本研究で採用した方法 L • ほぼ全ての演算が行列乗算 • 演算量:Ⅰ(2mn2)とⅡ(3mn2)の中間 • ブロック幅Lを適切に選ぶことで,他の2方法より高速になる可能性がある
Q = I – YTYT n I – × × Qの乗算 Ⅰ(単純なブロック化)・Ⅲ(本研究で採用した方法)の場合 Q = ( I – Yn/L Tn/L Yn/LT ) ・・・ (I – Y1T1Y1T ) L I – I – ・・・ × × × × Ⅱ(再帰的QR分解)の場合 • 行列サイズの大きいⅡの方がCSXLの性能を引き出せる • Ⅰ・Ⅲの方が使用するメモリの量が少ない
5. 性能評価 • Xeon 3.2GHz,メモリ2GB • g77 -O3 + Intel Math Kernel Library • ClearSpeed Advance ボード 評価環境 評価例題 • [-0.5,0.5]の一様乱数を要素とする m×n 乱数行列の特異値分解 • 10000 ≦ m ≦ 40000,1000 ≦ n ≦ 3000 評価項目 • ClearSpeed ボードでの3種のQR分解法の性能比較 • 特異値分解全体の ClearSpeed ボードによる高速化効果 • 精度評価
m = 10000,n = 3000 m = 20000,n = 3000 実行時間(秒) L=500 L=500 L=500 L=500 L=1000 L=1000 ClearSpeed ボードでの3種のQR分解法の性能比較 • 方法ⅠはQR分解が極めて遅い • 方法ⅡはQの乗算が最高速 • 方法ⅢはLを適切に選ぶことで,全体として最高速になる
m = 10000 m = 40000 実行時間(秒) n = 1000 n = 1000 n = 2000 n = 2000 特異値分解全体の ClearSpeed ボードによる高速化効果 • QR分解は最大で2倍程度高速化 • Qの乗算は最大で5倍程度高速化 • 特異値分解全体としては最大で2.3倍高速化
加速率 性能(MFLOPS) m n n m,nによる加速率の変化 特異値分解全体の加速率 nによる性能変化(m=30000) • 加速率はm, nが大きくなるほど増大する • m / n の増加: 特異値分解全体に対してCSXLを使える部分の割合が増加 • n の増加 : QR分解・Qの乗算の部分におけるCSXLの性能が向上 • ただし, CSXLを使わない部分(Rの特異値分解)の演算量も増加
m : m : n : n : 1000 1000 2000 2000 3000 3000 精度評価 左特異ベクトルの直交性 残差 ∥UTU – I∥F ∥US VT – A∥F • 左特異ベクトルの直交性,残差ともCSの方が精度が悪い • 悪化の程度は一桁以下
6. まとめと今後の課題 • ClearSpeed Advance ボードと CSXL を用いて,長方行列向け特異値分解の高速化を行った • QR分解のアルゴリズムを行列乗算を効果的に使う形に改良することで, ボードによる高速化の効果を高めた • 特異値分解全体では, 40000×2000 程度の行列においてボードを使うことで2.3倍の高速化が得られた • 行列の規模が大きくなるほど高速化が期待できる 本研究のまとめ 今後の課題 • より大規模な行列で性能評価 • SDKの使用による高速化 • 他の行列計算アルゴリズムへのCSX600の適用