マルチコア・メニーコアプロセッサ向けの 行列計算アルゴリズム － 固有値・特異値計算を中心に － マルチコア・メニーコアプロセッサ向けの 行列計算アルゴリズム 　－ 固有値・特異値計算を中心に － 2009年 4月 22日 研究会「アルゴリズムによる計算科学の融合と発展」 筑波大学 計算科学研究センター 山本　有作　（名古屋大学）

※ 対称密行列の固有値分解は特異値分解と共通点が多い はじめに： 本研究で扱う問題 固有値分解 A = XLX–1 A： n×n 非対称密行列，L： 対角行列，X： 正則行列 応用分野 構造計算 量子化学 特異値分解 A = US VT A： n×n 密行列，S： 対角行列，U, V： 直交行列 統計計算，信号処理，画像処理 電子状態計算 （filter diagonalization 法） ※ 対称密行列の固有値分解は特異値分解と共通点が多い

本研究の目的 マルチコア／メニーコアプロセッサ上で高性能を達成できる固有値・特異値分解ソルバの開発 背景 問題の大規模化 マルチコア／メニーコア／GPGPUの普及 Quad Core Opteron ClearSpeed CSX600 NVIDIA GeForce8800 GTX

目次 はじめに マルチコアプロセッサでの特異値分解 GPUでの特異値分解 メニーコアプロセッサでの固有値分解 終わりに

マルチコアプロセッサでの 特異値分解

標準的な特異値分解アルゴリズム 計算内容 計算手法 密行列 A U0TAV0 = B (U0, V0: 直交行列) ハウスホルダー法 二重対角化 二重対角行列 B 密行列Aをまず三重対角行列Tに相似変換してからTの固有値・固有ベクトルを求めるのが最も一般的な計算法。 三重対角化には，後に述べるハウスホルダー法を使う場合がほとんど。 三重対角行列の固有値・固有ベクトルの計算には，色々なアルゴリズムがある。 QR法 分割統治法 MR3アルゴリズム I-SVDアルゴリズム 二重対角行列の 特異値・特異ﾍﾞｸﾄﾙ計算 Bvi =σi xi BTxi =σi yi Bの特異値 {σi }， 　 特異ベクトル {xi }{yi } vi = V0 yi ui = U0 xi 逆変換 逆変換 Aの特異ベクトル {ui }, {vi }

各部分の演算量と実行時間 演算量 実行時間(全特異ﾍﾞｸﾄﾙ) 密行列 A 二重対角化 (8/3) n3 二重対角行列の n = 2500，Quad Core Xeon ×2 LAPACK での実行時間（秒） 二重対角化 (8/3) n3 密行列Aをまず三重対角行列Tに相似変換してからTの固有値・固有ベクトルを求めるのが最も一般的な計算法。 三重対角化には，後に述べるハウスホルダー法を使う場合がほとんど。 三重対角行列の固有値・固有ベクトルの計算には，色々なアルゴリズムがある。 実行時間（秒） 二重対角行列の 特異値・特異ﾍﾞｸﾄﾙ計算 O(n2) ～ O(n3) 逆変換 4mn2 （左右 m 本ずつの 特異ベクトル） コア数 ・二重対角化が実行時間の 　大部分を占める Aの特異ベクトル {ui }, {vi } ・速度向上率が低い

二重対角化の性能が出ない原因 二重対角化の演算パターン 演算パターンに関する問題点 左右からのハウスホルダー変換による行・列の消去 A(k) := (I – aw wT ) A(k) = A(k) – aw (wTA(k)) 演算は level-2 BLAS（行列ベクトル積と rank-1 更新） 演算パターンに関する問題点 Level-2 BLAS はデータ再利用性が低い。 　　　　キャッシュの有効利用が困難であり，単体性能が出にくい。 　　　　プロセッサ間のアクセス競合により，並列性能向上も困難 行列ベクトル積 Rank-1更新 非ゼロ要素 ゼロにしたい部分 A(k) 右から の変換 左から の変換 影響を受ける部分 k

マルチコアプロセッサのメモリ階層 主メモリアクセスの場合，実効性能は実質的に転送速度で決まる データ転送速度 非常に大 レジスタ 8～128 ワード コア1 コア2 データ転送速度大 本当はキャッシュも複数階層からなるが，ここでは単純化 キャッシュ 数100Kバイト ～ 数Mバイト 帯域の 奪い合い データ転送速度小 主メモリ 数100Mバイト ～ 数Gバイト 主メモリアクセスの場合，実効性能は実質的に転送速度で決まる Byte/Flop値: 1演算を行う時間で転送できるバイト数 < 1 マルチコアでは複数のコアが帯域を奪い合うため，問題が更に深刻化 主メモリのデータ転送速度の遅さをカバーするには，キャッシュにデータがある間に演算をまとめて行う（データの再利用）ことが必要 ～

BLASの利用による高性能化 BLASとは BLASの種類 A A A C A B Basic Linear Algebra Subprograms の略 個々のマシン向けにチューニングした基本行列演算のライブラリ BLASの種類 Level-1 BLAS:　ベクトルとベクトルの演算 内積　c := xT y AXPY演算　y: = ax + y など Level-2 BLAS:　行列とベクトルの演算 行列ベクトル積　y := Ax 行列のrank-1更新　A := A + xyT Level-3 BLAS:　行列と行列の演算 行列乗算　C := AB など A = × A A = ＋ × C A B = ×

演算をできる限り level-3 BLAS で行うことが高性能化のポイント 演算量：　O(N)，　データ量：　O(N) データ再利用性：　なし 並列粒度：　O(N/p)　（N： ベクトルの次元，p： プロセッサ台数） Level-2 BLAS 演算量：　O(N2)，　データ量：　O(N2) ベクトルデータのみ再利用性あり 並列粒度：　O(N2/p) Level-3 BLAS 演算量：　O(N3)，　データ量：　O(N2) O(N)回のデータ再利用が可能 行列乗算のサイズを大きくするほど，再利用性が向上 Byte/Flop 値の小さいマシンでも性能を出せる。 並列粒度：　O(N3/p) 演算をできる限り level-3 BLAS で行うことが高性能化のポイント = A × A = A ＋ × C A B = ×

Level-3 BLAS に基づく特異値分解アルゴリズム ２段階の二重対角化アルゴリズム（Bischof et al., ’93） 密行列 A をまず帯幅 L の下三角帯行列 C に変換 次にこの帯行列を下二重対角行列 B に変換 二重対角化を２段階で行うことの利点 前半の変換は，level-3 BLAS （行列乗算）のみを使って実行可能 　　キャッシュの有効利用が可能 後半の変換は level-2 BLAS が中心だが，演算量は O(n2L) 前半部に比べてずっと小さい。 次数 n 下三角 帯行列化 村田法 (8/3)n3 8n2L A C B 帯幅 L

下三角帯行列化のアルゴリズム ブロック鏡像変換によるブロック列の消去 ブロック鏡像変換 H = I – WaWT 第 K ステップでの処理 与えられたブロックベクトルを上三角 行列（正確には右上三角部分のみ 非零でそれ以外が零の行列）に変形 ブロックベクトル 行列 左からH を乗算 第 K ステップでの処理 右からHKR を乗算 左からHKL を乗算 非ゼロ要素 ゼロにしたい部分 影響を受ける部分 ・ 消去演算は行列乗算のみで行える ・ 行列乗算のサイズ ～ L

Level-3 BLAS に基づく特異値分解アルゴリズム 特異ベクトル計算を含めたアルゴリズムの全体像 長所 O(n3) の演算量を持つ部分はすべて level-3 BLAS で実行可能 短所 逆変換の演算量が 8mn2 （従来法の2倍）。ただし level-3 化可能 n L 特異値 {σi } (8/3)n3 8n2L QR法 DC法 MR3 I-SVD A C B 4mn2 4mn2 A の特異ﾍﾞｸﾄﾙ {ui }{vi } C の特異ﾍﾞｸﾄﾙ {zi }{wi } B の特異ﾍﾞｸﾄﾙ {xi }{yi }

性能評価 評価条件 計算機環境 問題サイズ n : 250, 5000, 7500 帯幅 L : 25, 50, 100 コア数 : 1, 2, 4, 8 Level-3 BLAS に基づく手法は Fortran で作成 従来法はLAPACKのルーチンを使用 計算機環境 Xeon X5355 (2.66GHz) Quad-core ×2ソケット Red Hat Enterprise WS release 4 Intel Math Kernel Library Intel Fortran Compiler 9 14

実験結果 Xeon での実行時間 n=2500 n=5000 n=7500 実行時間（秒） 結果（LAPACK(1CPU)に対する高速化率） コア数 ◆ LAPACK　■ Level-3 帯幅25　▲ Level-3 帯幅50　× Level-3 帯幅100 結果（LAPACK(1CPU)に対する高速化率） LAPACK ：8コアで最大 1.38倍 Level-3 BLAS ：8コアで最大 5.06倍 帯幅 L による性能差 L = 100 が最も性能が良い 15

実験結果 1.38倍 4.01倍 LAPACKとLevel-3 BLASの比較1 結果(コア数による高速化率) n = 2500，帯幅 L = 100 実行時間（秒） 1.38倍 4.01倍 コア数 LAPACK Level-3 BLAS 16

実験結果 1.31倍 4.36倍 LAPACKとLevel-3 BLASの比較2 実効時間の内訳 n = 7500，帯幅 L = 100 村田法の逆変換の実行時間が大きい　　　　今後改良要 n = 7500，帯幅 L = 100 1.31倍 実行時間（秒） 4.36倍 コア数 LAPACK Level-3 BLAS 17

実験結果 対称行列の固有値分解の結果 結果 特異値分解とほぼ同じ傾向が見られる。 n = 9000，帯幅 L = 100 実行時間（秒） コア数 LAPACK Level-3 BLAS 18

まとめ 特異値分解の標準的なアルゴリズムでは，二重対角化の部分で level-2 BLAS を多用するため，マルチコアプロセッサで十分な性能を得られない。 Level-3 BLAS を用いたアルゴリズムでは，演算量は増加するが，キャッシュの利用効率向上，アクセス競合の回避により，従来法より高い性能が得られる。 対称行列の固有値分解についても，level-3 BLAS を用いたアルゴリズムは効果的。

GPUでの特異値分解

研究背景 ◆ GPGPU （General Purpose GPU） NVIDIA GeForce8800 GTX （単精度演算：345.6GFLOPS） 一般の科学技術計算へのGPUの利用 単精度計算ではCPUの10倍以上の性能 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）

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 ユーザー自身がデータ転送を制御 ボードのデータはプログラム終了まで保持 メモリ配置の工夫によるデータ転送コストの削減 ◆ 性能 （GFLOPS） （Size） GeForce8800 GTX & CUBLAS ver. 1.0 ※転送時間は含んでいない × ＝ SGEMM（行列乗算，level-3 BLAS） SGEMV（行列ベクトル積，level-2 BLAS） Level-2 と 3 の性能差が CPU以上に大きい Level-3 BLAS活用の必要性

特異値分解の計算手順と特徴 A B (a) 二重対角化： (b) 二重対角行列の特異値分解： (c) 逆変換： ◆ 計算時間の内訳 (n) 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)

行列データは基本的にGPUに置き，必要な部分のみCPUに転送 高速化の方針 ◆ GPUの使い方 できるだけCUBLASを使う 行列乗算をできるだけ利用する 　Bischof のアルゴリズムの利用 メインメモリとGPU用メモリ間のデータ通信コストを抑える 行列データは基本的にGPUに置き，必要な部分のみCPUに転送 GPU向けのプログラムの移植は最小限にする ◆ 特異値分解の計算 二重対角化と逆変換の計算にはGPUを利用する Bの特異値分解の計算はCPUのみで行う 直交変換が中心のため，丸め誤差に強い 単精度計算でも安定性の問題がない 非線形方程式の求解あるいは漸化式の計算 安定性の面から倍精度計算が望ましい

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 サイズの小さい level-2 BLAS のコールが多数発生 (b) 二重対角行列の特異値分解： (c-1) 村田法の逆変換： 行列乗算により実行可能 演算量： 4n3 (c-2) 帯行列化の逆変換： 行列乗算により実行可能 演算量： 4n3

下三角帯行列化の実装法 CPU GPU ＝ 行列乗算 ＝ ブロック鏡像変換の作成 （BLASのみでは実行不可）

GPUでの村田法の実装法 ◆ 実装方針 ◆ 消去演算のパイプライン式並列化 ◆ GPU上での並列計算方法 （GeForce8800 GTX） サイズの小さい BLAS2 のコールが多数発生　　　CUBLASでなくCUDAで実装 ◆ 消去演算のパイプライン式並列化 第1列の二重対角化 更新範囲が重ならない ・・・ 第2列の二重対角化 ・・・ ◆ GPU上での並列計算方法　（GeForce8800 GTX） 16個のMPで消去演算をパイプライン式に並列実行 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での実行時間

精度評価 n n

まとめと今後の課題 本研究のまとめ 今後の課題 CUBLAS/CUDAを用いて正方行列の特異値分解をGPU上で実装 Bischof のアルゴリズムに基づき，演算量の大部分を行列乗算の形で実行することで，高速な CUBLAS の SGEMM ルーチンを活用 行列を GPU メモリに置くことで，データ転送量を削減 特異値分解全体でデュアルコアCPUに比べ4倍程度の高速化を実現 今後の課題 CPU と GPU の仕事の分担のさらなる効率化 新たな性能ネック部分の高速化 村田法（帯行列から二重対角行列への変換） 二重対角行列の特異値分解 最新の GPU や他のアクセラレータを用いた性能評価

メニーコアプロセッサでの 固有値分解

非対称行列の固有値計算 高速化対象 ヘッセンベルグQR法 対角成分 = 固有値 有限回演算 反復計算 密行列 ヘッセンベルグ行列 非対称行列を直交変換によりヘッセンベルグ行列に変換 ヘッセンベルグ行列の固有値を QR 法により求める 演算量の大部分を占める QR 法の高速化が必要 高速化対象 対角成分 = 固有値 有限回演算 反復計算 密行列 ヘッセンベルグ行列 右上三角行列 ハウスホルダー法 演算量 : （10/3）n3 QR 法 演算量 : 10 n3 （経験値）

QR 法の高速化 在来研究 本研究 マルチシフト QR 法 （K. Braman et al., 2002） 複数のシフトを導入 演算の大部分が行列乗算（Level-3 BLAS が利用可能） キャッシュの有効利用 並列性の向上 本研究 演算の大部分を占める行列乗算を高速化 浮動小数点アクセラレータの利用 ClearSpeed社 CSX600 1 チップに96個の演算コアを集積 最適化された行列乗算ライブラリ

マルチシフトQR法をそのまま使った場合 乱数行列 （n = 6000） マルチシフト QR 法により すべての固有値を計算 計算機環境 Xeon 3.2 GHz, Memory 8 GB ClearSpeed advance board CSX600 × 2 実行時間（秒） CSX600 （無） / （有） シフト数（ブロックサイズ）増加 問題点 実行時間の大部分を占める行列乗算（O(n3)）は，CSX600 により高速化 ただし，CSX600 の性能を引き出すには，シフト数を大きく取る必要あり Byte/Flop が 0.011（通常のCPUの1/10以下）と極端に小さい そのため，性能を出すには行列乗算のサイズが最低 1500 程度必要 その場合，行列乗算以外の部分（O(n2)）の演算量が増加してしまう

マルチシフトQR法の演算パターン k バルジ追跡型の演算 シフト数 m が増えることの効果 直交変換により，（m / 2）個のバルジ（出っ張り）を対角線に沿って k 行移動させる 最初に 対角ブロック内でバルジを移動 それに伴い，対角ブロック内部を更新 更新に用いた直交変換を1 個の行列に累積 最後に 非対角ブロックを まとめて更新 累積した行列 と 非対角ブロックの行列乗算 シフト数 m が増えることの効果 行列乗算のサイズが増加 対角ブロック内部の更新の演算量増大 k ボトルネック BLAS 3 バルジ（3×3） 逐次的に更新 まとめて更新 （BLAS 3 利用） 行列乗算部の性能向上 他の部分の実行時間増加

Level-3 BLASの利用可能部分の増大 アルゴリズムの改良（1） 再帰型アルゴリズムへの変換 k k / q 逐次的に更新 まとめて更新（BLAS 3 利用） （例：再帰段数 d = 1） (m / 2) / q 個の bulge を k / q 行 chasing Level-3 BLASの利用可能部分の増大

アルゴリズムの改良（2） 収束に伴う行列の分離 分離した小行列の固有値計算 更新処理の分割 ブロック化 右下　　　　　　 小行列が分離 （　　　　　　） 分離した小行列の固有値計算 ダブルシフト QR 法を適用 シフト数 m 増加により，小行列の次元増加 更新処理の分割 対角ブロックの更新 （右上三角化するまで） 更新に用いたハウスホルダー変換を 1 個の行列に累積 非対角ブロックの更新 （最後に 1 度だけ） 累積した行列 と 非対角ブロックの行列乗算 ボトルネック ブロック化 逐次的に更新 まとめて更新 （BLAS 3 利用） ・ Level-3 BLASの利用可能部分の増大 ・ 演算量削減

数値実験 テスト問題 解法 計算機環境 [0, 1] の成分を持つ 乱数行列 すべての固有値と固有ベクトルを計算 従来のマルチシフトQR法 [0, 1] の成分を持つ　　　　　乱数行列 ハウスホルダー法によりヘッセンベルグ化 すべての固有値と固有ベクトルを計算 解法 従来のマルチシフトQR法 改良型のマルチシフトQR法 計算機環境 Xeon 3.2 GHz , Memory 8 GB Fortran 77 倍精度演算 ClearSpeed advance board CSX600 × 2 Level-3 BLAS（行列乗算） CSXL Library (CSX600) Intel Math Kernel Library (CPU)

性能評価：アルゴリズムの効果 従来法，提案法ともに CSX600 を利用 実行時間（秒） 従来法に比べ，提案法は 約 1.4 倍高速 提案法 （n = 3000, m = 160, q = 4） （n = 6000, m = 200, q = 5） 従来法に比べ，提案法は 約 1.4 倍高速 対角ブロックの更新：約 1.5 倍高速化（level-3 BLAS化） 小行列の固有値計算：10 倍以上高速化（演算量削減）

性能評価：CSX600 の効果 実行時間（秒） CSX600 （+ 提案法） により n = 12000 n = 6000 （有） （有） m = 100 q = 5 m = 200 q = 5 m = 100 q = 5 m = 240 q = 6 CSX600 （+ 提案法） により n = 6000 のとき，約 3.5 倍高速化 n = 12000 のとき，約 3.8 倍高速化

まとめと今後の課題 本研究のまとめ 今後の課題 メニーコア型の浮動小数点アクセラレータ CSX600 を用いて非対称固有値計算のためのマルチシフトQR法を実装 アクセラレータの性能を引き出すには，ブロックサイズを大きく取る必要あり しかしその結果，level-3 BLAS 以外の(O(n2)の部分の)計算時間が増大 再帰型アルゴリズムへの変換などの手法を用いることで，O(n2)の部分についても level-3 BLAS 化を行い，全体を高速化 今後の課題 CPU とアクセラレータの仕事の分担のさらなる効率化 他のアクセラレータを用いた性能評価