スライド 0

正方行列向け特異値分解の
CUDAによる高速化
2009年 6月 24日
GPGPU 研究会
筑波大学計算科学研究センター
山本
深谷
畝山
中村
有作
猛
多加志
佳正
（名古屋大学）
（名古屋大学）
（京都大学）
（京都大学）
研究背景
◆ GPGPU （General Purpose GPU）
一般の科学技術計算へのGPUの利用
• 単純計算ではCPUの数十倍の性能
• CPUを上回る性能向上
• 開発環境の整備
NVIDIA GeForce8800 GTX
（単精度演算：345.6GFLOPS）
◆ 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）
（1/35）
研究目的
◆ 正方行列の特異値分解
•
•
•
•
A ： n×n
U ： n×n
S ： n×n
V ： n×n
密行列
直交行列
対角行列
直交行列
（応用分野：画像処理，情報検索，主成分分析など）
◆ 研究目的
CUDAによる正方行列の特異値分解計算の高速化
• 実装コストと汎用性の考慮
• GPUに適したアルゴリズムを選択
• 性能評価と課題の発見
※ GPUを使う部分の計算は単精度で行う
（2/35）
発表の流れ
1. 研究背景と目的
2. 高速化の方針
3. アルゴリズムの選択
4. 実装の概要
5. 性能評価
6. 終わりに
（3/35）
高速化の方針
CUDAによる行列計算の高速化
◆ GPU向けにプログラムの移植
• 拡張されたC言語でコーディングして，nvccでコンパイル
• 自由度の高いプログラミングが可能
• GPUの性能を十分に引き出すには様々なチューニングが必要
（スレッド数，条件分岐，メモリアクセスの連続性など）
◆ CUBLASの利用
• CUDAで提供されているBLAS（Basic Linear Algebra Subprograms）
• 標準のC言語のプログラム中で利用可能（gccでコンパイル可能）
• 限られた基本演算（行列ベクトル積，行列乗算など）のみ
• GPU向けにチューニング済み
今回はできるだけCUBLASを使って高速化を行う
（5/35）
CUBLASの特徴
◆ 仕様の特徴
メインメモリ
• ユーザー自身がデータ転送を制御
PCI-Express
• ボードのデータはプログラム終了まで保持
(1)Set Data
(3)Get Data
GPU用メモリ
(2)SGEMM etc.
メモリ配置の工夫によるデータ転送コストの削減
GPU
グラフィックボード
◆ 性能
（GFLOPS）
※転送時間は含んでいない
140
120
100
80
60
40
20
0
SGEMV（行列ベクトル積）
× ＝
SGEMM（行列乗算）
×
＝
SGEMMの有効利用
0
1000
2000
3000
4000
5000 （Size）
GeForce8800 GTX & CUBLAS ver. 1.0
（6/35）
特異値分解計算の特徴
◆ 計算手順と特徴
(a) 二重対角化：
• U1, V1：直交行列， B：下二重対角行列
• 大半がBLASルーチンにより計算可能
A
B
• 演算量：O(n3)
(b) 二重対角行列の特異値分解：
• 様々なアルゴリズム（QR法，分割統治法，MR3，I-SVDなど）
• 丸め誤差の影響を受けやすい
• 演算量： O(n2) ～ O(n3)
(c) 逆変換：
• 大半がBLASルーチンにより計算可能
• 演算量：O(n3)
（7/35）
特異値分解計算の特徴
◆ 計算時間の内訳
100%
逆変換
Bの特異値分解（I-SVD）
二重対角化
80%
60%
40%
20%
Core2 Duo (1.86GHz) & Intel MKL ver. 8.1
0%
1280
2560
3840
5120
(n)
◆ GPUを使う効果とコスト
• 二重対角化と逆変換：少ない実装コストで大きな効果が期待できる
 大半がBLASルーチンで計算可能 ⇒ CUBLASの利用が可能
 計算時間全体に占める割合が大きい ⇒高速化の効果が高い
• Bの特異値分解：実装コストは大きいが効果は小さい可能性が高い
 様々なアルゴリズムが存在 ⇒ 各々について検討が必要
 複雑な演算パターン ⇒ プログラムの移植が必要 ⇒ チューニングコスト大
 丸め誤差の影響を受けやすい ⇒ 単精度演算に適していない
（8/35）
高速化の方針
◆ GPUの使い方
• できるだけCUBLASを使う
 SGEMMをできるだけ利用する
 メインメモリとGPU用メモリ間のデータ通信コストを抑える
• GPU向けのプログラムの移植は最小限にする
◆ 特異値分解の計算
• 二重対角化と逆変換の計算にはGPUを利用する
• Bの特異値分解の計算はCPUのみで行う
（9/35）
アルゴリズムの選択
（二重対角化・逆変換）
従来法
◆ 鏡像変換による二重対角化
鏡像変換：
I－
A
・・・
B
◆ 逆変換
に対して
（11/35）
従来法の特徴
◆ 二重対角化
• 計算の逐次性
 鏡像変換の作成には直前のAの情報（ベクトル1本分）が必要
• 鏡像変換の作成のための通信コスト
 鏡像変換の作成はBLASルーチンのみで行えない（CPUで行う必要がある）
 鏡像変換の作用ごとに2回の通信（Aの情報の取得，鏡像変換の転送）
 通信するデータ量はベクトル1本分程度
• 鏡像変換の作用はSGEMVが中心
：行列ベクトル積（SGEMV）
：rank-1 更新（SGER）
CUBLASによる高速化の効果があまり期待できないのでは・・・？
（12/35）
Bischofの手法
◆ 二段階の二重対角化
(a-1) 下三角帯行列化：
• U11, V11：直交行列， C：半帯幅がLの下三角帯行列
• 大部分を行列乗算により計算可能
• 演算量： 8n3/ 3
A
C
C
B
（ただし n ≫ L ）
(a-2) 村田法：
• U12, V12：直交行列， B：下二重対角行列
• 演算量： 8n2L
(b) 二重対角行列の特異値分解：
(c-1) 村田法の逆変換：
• 演算量： 4n3
(c-2) 帯行列化の逆変換：
• 演算量： 4n3
（13/35）
Bischofの手法
◆ ブロック鏡像変換による下三角帯行列化
ブロック鏡像変換：
A
I－
・・・
C
（14/35）
Bischofの手法
◆ 村田法による帯行列の二重対角化
サイズの小さい鏡像変換によるbulge-chasing
• 第1列の二重対角化
C
・・・
• 第2列の二重対角化
・・・
・
・
・
（15/35）
Bischofの特徴
◆ 二重対角化
• 二重対角化の大部分は帯行列化
 帯行列化の演算はSGEMMが中心
• ブロック鏡像変換の作成のための通信コスト
 一回の通信量の増加（ベクトル→行列）
 通信回数の減少（1/L）
• 村田法にCUBLASを使っても効果が乏しい
 鏡像変換のサイズが小さい
 演算はSGEMVが中心
二重対角化全体としてはCUBLASを効果的に使えるのでは・・・？
◆ 逆変換
• 演算量が従来法の2倍に増加
逆変換のコストの増加の影響は・・・？
（16/35）
性能予測による比較
◆ CUBLASの性能測定
（ブロック）鏡像変換の作用では，段々と行列のサイズが小さくなる
演算量の合計
段々とサイズを変化させてSGEMM，SGEMVを実行
実行時間の合計
（FLOPS）
◆ 両手法の性能予測
各ステップの演算量と演算の種類から時間を予測
• n = 5120
予測時間 (sec)
70
逆変換（帯行列化）
• L = 64
60
• SGEMV : 3.54 GFLOPS
50
• SGEMM : 95.20 GFLOPS
40
逆変換
30
村田法
20
帯行列化
• 村田法はCPUでの実際の実行時間
Bischofの手法の方が効果が期待できる
逆変換（村田法）
10
二重対角化
0
従来法
Bischof
（17/35）
実装の概要
Bischofの手法の実装
◆ GPUを利用する部分
(a-1) 下三角帯行列化
CPU
(a-2) 村田法
GPU向けに移植 (nvcc)
(b) 二重対角行列の特異値分解
CPU
(c-1) 村田法の逆変換
CPU
(c-2) 帯行列化の逆変換
CUBLAS
CUBLAS
CUBLAS
（19/35）
下三角帯行列化
CPU
GPU
＝
＝
SGEMM
（20/35）
帯行列化の逆変換
CPU
GPU
SGEMM
※ V も同様にして計算
（21/35）
村田法の逆変換
CPU
GPU
SGEMM
SGEMM
・
・
・
SGEMMの性能が出るサイズに合成
・
・
・
SGEMM
※ V も同様にして計算
（22/35）
GPUへの村田法の移植
◆ bulge-chasingのパイプライン式並列化
• 第1列の二重対角化
・・・
• 第2列の二重対角化
更新範囲が重ならない
・・・
◆ GPU上での並列計算方法（GeForce8800 GTX）
• 16個のMPでbulge-chasingをパイプライン式に並列実行
• MP内の8個のSPで鏡像変換の作用を並列計算（MP内の共有メモリを利用）
（23/35）
特異値分解計算の全体の様子
CPU
GPU
A
A
ブロック鏡像変換の作成
ブロック鏡像変換の作用
C
C
CUBLAS
鏡像変換の作成・作用
B
H
B
H
特異値分解
U2
V2
U2
鏡像変換の合成
Q
V2
合成結果の作用
Q
U’
V’
ブロック鏡像変換の作用
U
V
U
CUBLAS
CUBLAS
V
（24/35）
性能評価
評価方法
◆ 評価問題
• [-0.5 , 0.5]の乱数を要素とするn×nの正方行列の特異値分解
• n = 1280, 2560, 3840, 5120
◆ 手法
• 二重対角行列の特異値分解は全て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
（26/35）
二重対角化と逆変換の評価
◆ n=1280
実行時間 (sec)
10.00
逆変換（帯行列化）
逆変換（村田法）
逆変換
村田法
9.00
8.00
7.00
6.00
帯行列化
二重対角化
5.00
4.00
3.00
2.9倍の高速化
2.00
1.00
0.00
従来法
Bischof(32) Bischof(64) Bischof(128) Bischof(32) Bischof(64) Bischof(128)
CPU（2コア）
CPU（1コア） & GPU
（27/35）
特異分解計算全体の評価
◆ n=1280
実行時間 (sec)
10.00
逆変換（帯行列化）
逆変換（村田法）
逆変換
I-SVD
村田法
帯行列化
二重対角化
9.00
8.00
7.00
6.00
5.00
1.9倍の高速化
4.00
3.00
2.00
1.00
0.00
従来法
Bischof(32) Bischof(64) Bischof(128) Bischof(32) Bischof(64) Bischof(128)
CPU（2コア）
CPU（1コア） & GPU
（28/35）
二重対角化と逆変換の評価
◆ n=5120
実行時間 (sec)
400
逆変換（帯行列化）
逆変換（村田法）
逆変換
村田法
350
300
250
帯行列化
二重対角化
200
150
100
7.0倍の高速化
50
0
従来法
Bischof(32)
Bischof(64) Bischof(128) Bischof(32)
CPU（2コア）
Bischof(64) Bischof(128)
CPU（1コア） & GPU
（29/35）
特異分解計算全体の評価
◆ n=5120
実行時間 (sec)
400
逆変換（帯行列化）
逆変換（村田法）
逆変換
I-SVD
村田法
帯行列化
二重対角化
350
300
250
200
150
4.2倍の高速化
100
50
0
従来法
Bischof(32)
Bischof(64) Bischof(128) Bischof(32)
CPU（2コア）
Bischof(64) Bischof(128)
CPU（1コア） & GPU
（30/35）
Bischofの手法のステップごとの評価
◆ L=64の場合
Speedup
25.0
帯行列化
村田法
逆変換（村田法）
逆変換（帯行列化）
20.0
15.0
10.0
5.0
0.0
1280
2560
3840
5120
(n)
Speedup = CPU (2コア)での実行時間 / CPU (1コア) & GPUでの実行時間
（31/35）
性能予測の評価
時間 (sec)
n = 5120 , L = 64
70
60
逆変換（帯行列化）
逆変換（村田法）
逆変換
村田法
帯行列化
二重対角化
村田法はCPU
50
40
30
20
10
0
従来法
Bischof
予測
Bischof
Bischof
実際
• 予測性能は実際の性能の上限を見積もっている
• 従来法の予測性能＜ Bischofの手法の実際の性能
Bischofの手法の選択は適切であったと言える
（32/35）
精度評価
1.0E+00
CPU_LAPACK (double)
CPU_LAPACK (single)
CPU_Bischof (single)
CPU & GPU_Bischof (single)
1.0E-01
1.0E-02
1.0E-03
1.0E+00
1.0E-01
1.0E-02
1.0E-03
1.0E-04
1.0E-04
1.0E-05
1.0E-05
1.0E-06
1.0E-06
1.0E-07
1.0E-07
1.0E-08
1.0E-08
1.0E-09
1.0E-09
1.0E-10
1.0E-10
1.0E-11
1.0E-11
1.0E-12
1280
3840
2560
n
5120
1280
3840
2560
5120
n
（33/35）
終わりに
まとめと今後の課題
◆ まとめ
• CUDAを用いた正方行列の特異値分解計算の高速化
• CUBLASのSGEMMの利用を中心とした高速化
• 性能予測により，Bischofの手法を二重対角化・逆変換の手法として選択
• メモリ配置の工夫によるデータ転送コストの削減
• CPU上での事前計算による，CUBLASのSGEMMの有効利用
• 数値実験による性能評価（特異値分解全体で最大4倍程度の高速化）
◆ 今後の課題
• CPUとGPUの仕事の分担のさらなる効率化
• 適切な半帯幅の決定方法
• 村田法の高速化
• 最新バージョンのCUDAや他のアクセラレータを用いた性能評価
• 対称行列の固有値計算への適用
（35/35）

Download Report