正方行列向け特異値分解の
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)
ダウンロード

PPT file