正方行列向け特異値分解の CUDAによる高速化 2009年 6月 24日 GPGPU 研究会　 筑波大学 計算科学研究センター 山本　有作　 （名古屋大学） 深谷　猛 （名古屋大学） 畝山　多加志 （京都大学） 中村　佳正　 （京都大学）

研究背景 ◆ GPGPU （General Purpose GPU） NVIDIA GeForce8800 GTX （単精度演算：345.6GFLOPS） 一般の科学技術計算へのGPUの利用 単純計算ではCPUの数十倍の性能 CPUを上回る性能向上 開発環境の整備 ・GPGPUが盛ん ・CUDAによってより身近に ・ ◆ CUDA　（Compute Unified Device Architecture） GPGPUのための統合開発環境 2006年にNVIDIA社が発表 C/C++の拡張によりGPGPUのプログラミングが可能 GPU向けにチューニング済みのライブラリ（BLAS，FFT） CUDAの利用例：行列計算，重力多体計算，分子軌道計算，流体計算など （Cf. http://www.nvidia.co.jp/object/cuda_home_jp.html）

研究目的 ◆ 正方行列の特異値分解 ◆ 研究目的 （応用分野：画像処理，情報検索，主成分分析など） A ： n×n 密行列 U ： n×n 直交行列 S ： n×n 対角行列 V ： n×n 直交行列 （応用分野：画像処理，情報検索，主成分分析など） ◆ 研究目的 CUDAによる正方行列の特異値分解計算の高速化 実装コストと汎用性の考慮 GPUに適したアルゴリズムを選択 性能評価と課題の発見 ※ GPUを使う部分の計算は単精度で行う

発表の流れ 研究背景と目的 高速化の方針 アルゴリズムの選択 実装の概要 性能評価 終わりに

高速化の方針

CUDAによる行列計算の高速化 ◆ GPU向けにプログラムの移植 ◆ CUBLASの利用 今回はできるだけCUBLASを使って高速化を行う 拡張されたC言語でコーディングして，nvccでコンパイル 自由度の高いプログラミングが可能 GPUの性能を十分に引き出すには様々なチューニングが必要 （スレッド数，条件分岐，メモリアクセスの連続性など） 今回はできるだけCUBLASを使って高速化を行う ◆ CUBLASの利用 CUDAで提供されているBLAS（Basic Linear Algebra Subprograms） 標準のC言語のプログラム中で利用可能　（gccでコンパイル可能） 限られた基本演算（行列ベクトル積，行列乗算など）のみ GPU向けにチューニング済み

CUBLASの特徴 ◆ 仕様の特徴 ◆ 性能 × ＝ ユーザー自身がデータ転送を制御 ボードのデータはプログラム終了まで保持 メインメモリ グラフィックボード (1)Set Data GPU用メモリ GPU (3)Get Data (2)SGEMM etc. PCI-Express ユーザー自身がデータ転送を制御 ボードのデータはプログラム終了まで保持 メモリ配置の工夫によるデータ転送コストの削減 ◆ 性能 × ＝ SGEMM（行列乗算） SGEMV（行列ベクトル積） （GFLOPS） （Size） GeForce8800 GTX & CUBLAS ver. 1.0 ※転送時間は含んでいない SGEMMの有効利用

特異値分解計算の特徴 A B ◆ 計算手順と特徴 (a) 二重対角化： (b) 二重対角行列の特異値分解： (c) 逆変換： U1, V1：直交行列，　B：下二重対角行列 大半がBLASルーチンにより計算可能 演算量：O(n3) A B (b) 二重対角行列の特異値分解： 様々なアルゴリズム（QR法，分割統治法，MR3，I-SVDなど） 丸め誤差の影響を受けやすい 演算量： O(n2) ～ O(n3) (c) 逆変換： 大半がBLASルーチンにより計算可能 演算量：O(n3)

Core2 Duo (1.86GHz) & Intel MKL ver. 8.1 特異値分解計算の特徴 ◆ 計算時間の内訳 (n) Core2 Duo (1.86GHz) & Intel MKL ver. 8.1 ◆ GPUを使う効果とコスト 二重対角化と逆変換 ： 少ない実装コストで大きな効果が期待できる Bの特異値分解 ： 実装コストは大きいが効果は小さい可能性が高い 大半がBLASルーチンで計算可能　⇒　CUBLASの利用が可能 計算時間全体に占める割合が大きい　⇒高速化の効果が高い 様々なアルゴリズムが存在　⇒　各々について検討が必要 複雑な演算パターン　⇒　プログラムの移植が必要　⇒　チューニングコスト大 丸め誤差の影響を受けやすい　⇒　単精度演算に適していない

高速化の方針 ◆ GPUの使い方 ◆ 特異値分解の計算 できるだけCUBLASを使う GPU向けのプログラムの移植は最小限にする SGEMMをできるだけ利用する メインメモリとGPU用メモリ間のデータ通信コストを抑える GPU向けのプログラムの移植は最小限にする ◆ 特異値分解の計算 二重対角化と逆変換の計算にはGPUを利用する Bの特異値分解の計算はCPUのみで行う

アルゴリズムの選択 （二重対角化・逆変換）

従来法 ◆ 鏡像変換による二重対角化 鏡像変換 ： I － A B ・・・ ◆ 逆変換 に対して

CUBLASによる高速化の効果があまり期待できないのでは・・・？ 従来法の特徴 ◆ 二重対角化 計算の逐次性 鏡像変換の作成には直前のAの情報（ベクトル1本分）が必要 鏡像変換の作成のための通信コスト 鏡像変換の作成はBLASルーチンのみで行えない（CPUで行う必要がある） 鏡像変換の作用ごとに2回の通信（Aの情報の取得，鏡像変換の転送） 通信するデータ量はベクトル1本分程度 逆変換では鏡像変換をブロック化して、SGEMMが使える 鏡像変換の作用はSGEMVが中心 ：行列ベクトル積　（SGEMV） ：rank-1 更新　（SGER） CUBLASによる高速化の効果があまり期待できないのでは・・・？

Bischofの手法 A C B C ◆ 二段階の二重対角化 (a-1) 下三角帯行列化： (a-2) 村田法： U11, V11：直交行列，　C：半帯幅がLの下三角帯行列 大部分を行列乗算により計算可能 演算量： 8n3/ 3　　（ただし n ≫ L ） (a-2) 村田法： B C U12, V12：直交行列，　B：下二重対角行列 演算量： 8n2L (b) 二重対角行列の特異値分解： (c-1) 村田法の逆変換： 演算量： 4n3 (c-2) 帯行列化の逆変換： 演算量： 4n3

Bischofの手法 ◆ ブロック鏡像変換による下三角帯行列化 ブロック鏡像変換 ： I － A ・・・ C

Bischofの手法 C ◆ 村田法による帯行列の二重対角化 ・・・ サイズの小さい鏡像変換によるbulge-chasing 第1列の二重対角化 ・・・ C 第2列の二重対角化 ・・・ ・・・

二重対角化全体としてはCUBLASを効果的に使えるのでは・・・？ Bischofの特徴 ◆ 二重対角化 二重対角化の大部分は帯行列化 帯行列化の演算はSGEMMが中心 ブロック鏡像変換の作成のための通信コスト 一回の通信量の増加（ベクトル→行列） 通信回数の減少　（1/L） 村田法にCUBLASを使っても効果が乏しい 鏡像変換のサイズが小さい 演算はSGEMVが中心 二重対角化全体としてはCUBLASを効果的に使えるのでは・・・？ ◆ 逆変換 演算量が従来法の2倍に増加 逆変換のコストの増加の影響は・・・？

Bischofの手法の方が効果が期待できる 性能予測による比較 ◆ CUBLASの性能測定 （ブロック）鏡像変換の作用では，段々と行列のサイズが小さくなる 段々とサイズを変化させてSGEMM，SGEMVを実行 演算量の合計 実行時間の合計 （FLOPS） ◆ 両手法の性能予測 予測時間 (sec) n = 5120 L = 64 SGEMV : 3.54 GFLOPS SGEMM : 95.20 GFLOPS 村田法はCPUでの実際の実行時間 各ステップの演算量と演算の種類から時間を予測 Bischofの手法の方が効果が期待できる

実装の概要

Bischofの手法の実装 ◆ GPUを利用する部分 (a-1) 下三角帯行列化 (a-2) 村田法 (b) 二重対角行列の特異値分解 CPU CUBLAS (a-2) 村田法 GPU向けに移植 (nvcc) (b) 二重対角行列の特異値分解 CPU (c-1) 村田法の逆変換 CPU CUBLAS (c-2) 帯行列化の逆変換 CUBLAS

下三角帯行列化 CPU GPU ＝ SGEMM ＝

帯行列化の逆変換 CPU GPU SGEMM ※ V も同様にして計算

村田法の逆変換 CPU GPU SGEMM ※ V も同様にして計算 ・・・ SGEMM SGEMMの性能が出るサイズに合成 ・・・

GPUへの村田法の移植 ◆ bulge-chasingのパイプライン式並列化 第1列の二重対角化 更新範囲が重ならない ・・・ 第2列の二重対角化 ・・・ ◆ GPU上での並列計算方法　（GeForce8800 GTX） 16個のMPでbulge-chasingをパイプライン式に並列実行 MP内の8個のSPで鏡像変換の作用を並列計算（MP内の共有メモリを利用）

特異値分解計算の全体の様子 A A C Q C B H B H U’ V’ U2 V2 U2 V2 Q U V U V CPU GPU ブロック鏡像変換の作成 ブロック鏡像変換の作用 CUBLAS C 鏡像変換の作成・作用 B H B H 鏡像変換の合成 U’ V’ 合成結果の作用 CUBLAS 特異値分解 U2 V2 U2 V2 Q U V ブロック鏡像変換の作用 CUBLAS U V

性能評価

評価方法 ◆ 評価問題 ◆ 手法 ◆ 実行環境 [-0.5 , 0.5]の乱数を要素とするn×nの正方行列の特異値分解 二重対角行列の特異値分解は全てI-SVDの倍精度計算で行う 二重対角化と逆変換の計算法は以下の通り（全て単精度で行う） 従来法（LAPACKのルーチン）をCPUのみで実行 Bischofの手法（自作プログラム）をCPUのみで実行 Bischofの手法（自作プログラム）をCPUとGPUで実行 ※ Bischofの手法における半帯幅Lは32, 64, 128 ◆ 実行環境 CPU ： Intel Core2 Duo (1.86GHz), Intel MKL ver. 8.1, gcc オプション -O3 GPU ： NVIDIA GeForce8800 GTX, CUBLAS ver. 1.0, nvcc ver. 1.0

二重対角化と逆変換の評価 ◆ n=1280 実行時間 (sec) 2.9倍の高速化 CPU（2コア） CPU（1コア） & GPU

特異分解計算全体の評価 ◆ n=1280 実行時間 (sec) 1.9倍の高速化 CPU（2コア） CPU（1コア） & GPU

二重対角化と逆変換の評価 ◆ n=5120 実行時間 (sec) 7.0倍の高速化 CPU（2コア） CPU（1コア） & GPU

特異分解計算全体の評価 ◆ n=5120 実行時間 (sec) 4.2倍の高速化 CPU（2コア） CPU（1コア） & GPU

Bischofの手法のステップごとの評価 ◆ L=64の場合 Speedup (n) Speedup = CPU (2コア)での実行時間 / CPU (1コア) & GPUでの実行時間

性能予測の評価 Bischofの手法の選択は適切であったと言える 予測性能は実際の性能の上限を見積もっている 時間 (sec) n = 5120 , L = 64 村田法はCPU 予測 実際 予測性能は実際の性能の上限を見積もっている 従来法の予測性能　＜　Bischofの手法の実際の性能 Bischofの手法の選択は適切であったと言える

精度評価 n n

終わりに

まとめと今後の課題 ◆ まとめ ◆ 今後の課題 CUDAを用いた正方行列の特異値分解計算の高速化 CUBLASのSGEMMの利用を中心とした高速化 性能予測により，Bischofの手法を二重対角化・逆変換の手法として選択 メモリ配置の工夫によるデータ転送コストの削減 CPU上での事前計算による，CUBLASのSGEMMの有効利用 数値実験による性能評価（特異値分解全体で最大4倍程度の高速化） ◆ 今後の課題 CPUとGPUの仕事の分担のさらなる効率化 適切な半帯幅の決定方法 村田法の高速化 最新バージョンのCUDAや他のアクセラレータを用いた性能評価 対称行列の固有値計算への適用