長方行列向け特異値分解の 浮動小数点コプロセッサによる 高速化 2007年1月18日　HPCS2007 深谷 猛，山本 有作 （名古屋大学） 畝山 多加志 （京都大学） 堀 玄，梅野 健 （理化学研究所）

概要 背景と目的 CSX600のアーキテクチャと性能 長方行列の特異値分解アルゴリズム CSX600向けの高速化手法 性能評価 まとめと今後の課題

数万×数千以上の行列の特異値分解 高速化が必要 １.　背景と目的 長方行列の特異値分解 応用 画像処理 電子状態計算 （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バス経由でコプロセッサとして使う場合は，特に顕著 　 メモリネックを軽減できるアルゴリズムが性能向上の鍵

Level-3 BLASの利用 行列乗算 C:=C+AB = C A B 演算量に比べ，データ量は O(1/N) キャッシュをうまく使うことで，メモリネックを軽減可能 行列乗算を効率的に使えれば，一般のアルゴリズムも高速化可能 （LAPACKの基本的な考え方） = × + C A B データ量： O(N 2) 演算量： 　O(N 3) ※行列ベクトル積（ y := y + Ax）では，データ量・演算量ともにO(N 2)

本研究の目的 長方行列の特異値分解アルゴリズムをCSX600を利用して高速化する 性能評価を行い，更なる高速化に向けての課題を明らかにする

２. CSX600のアーキテクチャと性能 CSX600 チップ ClearSpeed Advance ボード 1個の主プロセッサ 96個の演算専用プロセッサ 64ビット 倍精度2演算/サイクル 128B レジスタファイル 6KB SRAM 250MHz で動作 ピーク性能： 48GFLOPS ClearSpeed Advance ボード 2個の CSX600 プロセッサ 1GB DRAM PCI-X バスにより PC 本体と接続 ピーク性能： 96GFLOPS

CSX600の利用環境 Software Development Kit CSXL ライブラリ CSFFT ライブラリ コンパイラ：　Cn 言語によるチップ内並列プログラミング デバッガ シミュレータ CSXL ライブラリ ClearSpeed Advance ボード用の BLAS 行列データを PC の主メモリ上に置いてコール PC ⇔ ボード間の転送はライブラリ内で実行 公称性能： DGEMM（行列乗算）で sustained 50GFLOPS 現状では非常に制約が多い 使用可能なルーチンは DGEMM と DGETRF のみ DGEMM では，M，N，K ≧ 448 でないと動かない CSFFT ライブラリ 今回はこれを利用

CSXL DGEMMの性能 (1) サイズパラメータのうち2個が450の場合 C += A B 性能は最大でも 11GFLOPS 程度 他の2個が450の場合でも，ほぼ同じ結果 N K C += A × B M 性能（MFLOPS） ・ M = K = 450 ・ 1000 ≦ N ≦ 6000 N CSX600の性能が活かせない

CSXL DGEMMの性能 (2) サイズパラメータのうち1個のみが450の場合 C += A B M = N が大きい場合，性能は 40GFLOPS 以上 他の2個が大きい場合でも，ほぼ類似の性能 N K C += A × B 性能（MFLOPS） M ・ K = 450 ・ 1000 ≦ M = N ≦ 6000 N 性能を引き出すには，サイズパラメータのうち少なくとも2個を 大きい値にとる必要あり

３. 長方行列の特異値分解アルゴリズム = QR分解 A=QR 2mn2 二重対角化 (8/3)n3 特異値分解 O(n2) ~ O(n3) ３.　長方行列の特異値分解アルゴリズム 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 の場合 QR分解とQの乗算の時間が支配的 環境 ・ Xeon 3.2GHz ・ g77 -O3 + Intel MKL 各部分のアルゴリズム（詳細は後述） ・ QR分解： ブロック化 + 再帰的QR ・ 二重対角化： Bischof + 村田法 ・ 特異値分解： 分割統治法 実行時間（秒） QR分解とQの乗算の時間が支配的

４. CSX600向けの高速化手法 高速化の方針 QR分解とQの乗算 二重対角化および逆変換 二重対角行列の特異値分解 行列乗算を効率的に利用できるようアルゴリズムを変更 CSXLで実行 二重対角化および逆変換 二重対角化はBischof のアルゴリズムと村田法を利用 CSXLは利用しない（CPUのみで実行） 二重対角行列の特異値分解 行列乗算中心のアルゴリズムである分割統治法を利用 LAPACKのDBDSDCを利用（CSXLが使える部分は使う）

ハウスホルダー変換によるQR分解 ハウスホルダー変換による列の消去 H1 A = ( I – t1 y1 y1T ) A = A(1) level-2 BLAS CSXLを使えない 上三角化と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 – t2 y2 y2T )( I – t1 y1 y1T ) = I – Yn Tn YnT ただし， Yn = [ y1 | y2 | ・・・ | yn]　（m×n 行列） 　　　　　 Tn：　n×n 下三角行列 I – ・・・ I – = I – × × × × × × 複数のハウスホルダー変換の作用を，行列乗算として まとめて実行可能

QR分解における行列乗算の利用法 Ⅰ 単純なブロック化（LAPACKで使用） L ブロックのQR分解はlevel-2 BLAS I – Y1T1Y1T I – Y2T2Y2T I – Y3T3Y3T 分解 更新 H1 H2 H3 H3H2H1= I – Y1T1Y1T ブロックのQR分解はlevel-2 BLAS 演算量： L (ブロック幅) << n のとき 2mn2 CSXLを使うにはL≧450 ⇒ level-2 BLASの演算量の増加

QR分解における行列乗算の利用法 Ⅱ 再帰的QR分解 ほぼ全ての演算が行列乗算 演算量： 3mn2 I – YR TR YRT I – YL TL YLT 分解 更新 合成 (I – YL TL YLT ) ×(I – YR TR YRT ) = I – Y T YT I – YL TL YLT I – Y’L T’L Y’LT ほぼ全ての演算が行列乗算 演算量： 3mn2

QR分解における行列乗算の利用法 Ⅲ 本研究で採用した方法 L ほぼ全ての演算が行列乗算 演算量：Ⅰ(2mn2)とⅡ(3mn2)の中間 I – Y1T1Y1T I – Y2T2Y2T I – Y3T3Y3T 分解 更新 I – Y1 T1 Y1T ほぼ全ての演算が行列乗算 演算量：Ⅰ(2mn2)とⅡ(3mn2)の中間 ブロック幅Lを適切に選ぶことで，他の２方法より高速になる可能性がある

Q の乗算 Ⅰ（単純なブロック化）・Ⅲ（本研究で採用した方法）の場合 Q = ( I – Yn/L Tn/L Yn/LT ) ・・・ (I – Y1 T1 Y1T ) L I – × ・・・ I – × × × Ⅱ（再帰的QR分解）の場合 Q = I – Y T YT × I – n 行列サイズの大きいⅡの方がCSXLの性能を引き出せる Ⅰ・Ⅲの方が使用するメモリの量が少ない

５. 性能評価 評価環境 評価例題 評価項目 Xeon 3.2GHz，メモリ2GB ５.　性能評価 評価環境 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 ボードによる高速化効果 精度評価

ClearSpeed ボードでの3種のQR分解法の性能比較 実行時間（秒） L=500 L=1000 m = 10000，n = 3000 L=500 L=1000 m = 20000，n = 3000 方法ⅠはQR分解が極めて遅い 方法ⅡはQの乗算が最高速 方法ⅢはLを適切に選ぶことで，全体として最高速になる

特異値分解全体の ClearSpeed ボードによる高速化効果 実行時間（秒） m = 10000 n = 1000 n = 2000 m = 40000 n = 1000 n = 2000 QR分解は最大で2倍程度高速化 Qの乗算は最大で5倍程度高速化 特異値分解全体としては最大で2.3倍高速化

m，n による加速率の変化 特異値分解全体の加速率 nによる性能変化（m=30000） n 性能（MFLOPS） 加速率 m n m / n の増加：　特異値分解全体に対してCSXLを使える部分の割合が増加 n の増加　 ：　QR分解・Qの乗算の部分におけるCSXLの性能が向上 　　　　 ただし， CSXLを使わない部分（Rの特異値分解）の演算量も増加

精度評価 ∥US VT – A∥F ∥UTU – I∥F 左特異ベクトルの直交性 残差 左特異ベクトルの直交性，残差ともCSの方が精度が悪い m ： n ： 1000 2000 3000 m ： n ： 1000 2000 3000 左特異ベクトルの直交性，残差ともCSの方が精度が悪い 悪化の程度は一桁以下

６. まとめと今後の課題 本研究のまとめ 今後の課題 ６.　まとめと今後の課題 本研究のまとめ ClearSpeed Advance ボードと CSXL を用いて，長方行列向け特異値分解の高速化を行った QR分解のアルゴリズムを行列乗算を効果的に使う形に改良することで， ボードによる高速化の効果を高めた 特異値分解全体では， 40000×2000 程度の行列においてボードを使うことで2.3倍の高速化が得られた 行列の規模が大きくなるほど高速化が期待できる 今後の課題 より大規模な行列で性能評価 SDKの使用による高速化 他の行列計算アルゴリズムへのCSX600の適用