応用数理工学特論　第10回 日立製作所の山本有作です。 「～」について発表いたします。 2008年7月3日 計算理工学専攻 山本有作

２. 共有メモリ型並列計算機での特異値分解の高速化 目次 １.　はじめに ２.　共有メモリ型並列計算機での特異値分解の高速化 ～ 正方行列の特異値分解 ～ ３.　SIMD超並列プロセッサでの特異値分解の高速化 ～ 長方行列の特異値分解 ～ 本発表では，はじめに，研究の背景を述べてから，スパースソルバの概要，並列化手法，そして本研究で工夫した点の一つであるRISCプロセッサ向けの最適化についてご説明します。最後に，並列計算機SR2201上での性能評価とまとめを述べます。

１.　はじめに

特異値分解とは 任意の実正方行列 A は，直交行列 U，対角行列 S，直交行列 VT の積に分解できる。 参考 A が実対称行列の場合，A=UDUT A が任意の実行列の場合，A=SDS–1 A が任意の実行列の場合，A=URUT A ：n×n 実行列 U ：n×n 直交行列 V ： n×n 直交行列 S ： n×n 対角行列 n = × × n n

特異値と特異ベクトル A=USV T において， 特異ベクトルの性質 V の各列 vi を右特異ベクトル U の各列 ui を左特異ベクトル S の要素 si を特異値，と呼ぶ。 特異ベクトルの性質 AV=US より， Avi = si ui UTA=S V T より， uiTA = si viT

固有値分解（対角化）との関係 A=USV T より， ATA= VSUTUSV T = VS 2V T 　 よって，V は ATA を対角化する直交行列 　 A の右特異ベクトル = ATA の固有ベクトル A の特異値 = ATA の固有値の平方根 一方， AAT= USV TVSU T = US 2U T 　よって，V は ATA を対角化する直交行列 　 A の左特異ベクトル = AAT の固有ベクトル

特異値分解を求めるアルゴリズム 二重対角行列への変換 二重対角行列の特異値分解 U0TAV0 = B (U0, V0: 直交行列) よって，U= U0U1 ，V= V0V1 とおくと（逆変換）， U0TAV0 = B　(U0, V0: 直交行列) U1TBV1 = S　(U1, V1: 直交行列) U1TU0TAV0V1 = S A = USV T

特異値分解の応用 信号処理 画像処理 電子状態計算 文書検索 各種統計計算 Filter Diagonalization Method Latent Semantic Indexing 各種統計計算 主成分分析 独立成分分析 最小二乗法

対象とする計算機

共有メモリ型並列計算機 アーキテクチャ 例１： マルチコアプロセッサ 例２： 共有メモリ型のスーパーコンピュータ 複数のプロセッサ（PU）がバスを 　　通してメモリを共有 PUはそれぞれキャッシュを持つ 例１： マルチコアプロセッサ Intel社 デュアルコアXeon 1チップに2個のプロセッサを共有メモリ型で集積 例２： 共有メモリ型のスーパーコンピュータ 富士通 PrimePower HPC2500 最大128個のPUを共有メモリで結合 キャッシュ PU0 PU1 PU2 PU3 バス メモリ デュアルコア Xeon PrimePower HPC2500

SIMD超並列プロセッサ 1チップに100個単位のプロセッサを搭載 全プロセッサが別々のデータに対して同じ命令を実行する SIMD（Single Instruction Multiple Data）型の並列計算機 線形計算に最適 汎用プロセッサの10～100倍の性能 ClearSpeed社CSX600:　96プロセッサ，48GFLOPS（倍精度） 東大GRAPE-DR:　512プロセッサ，256GFLOPS（倍精度） 演算性能は極めて高いが，データ転送速度は汎用プロセッサと同程度 データ再利用性の向上が，共有メモリ型並列計算機以上に重要

512コアプロセッサ（GRAPE-DR）

グラフィックスプロセッサ（GPU） GPU（Graphics Processing Unit）の高速化 CPUを大きく上回るペースで演算性能が向上 グラフィックスメモリも大容量化・高速化 nVIDIA社 GeForce8800GTX ・ 240個の演算プロセッサ ・ 演算性能： 933GFLOPS（単精度） 　　　　　　　　 90GFLOPS（倍精度） ・ メモリ： 1GB GPUを汎用の数値計算に使うGPGPU（General-Purpose GPU）が注目を集める 計算機としてのアーキテクチャはSIMD超並列型

２. 共有メモリ型並列計算機での特異値分解の高速化 ２.　共有メモリ型並列計算機での特異値分解の高速化 ～ 正方行列の特異値分解 ～

正方行列の特異値分解 = ＜応用＞ A ：n×n 密行列 U ：n×n 直交行列 V ： n×n 直交行列 S ： n×n 対角行列 n × 各種統計計算 より応用の多い長方行列の特異値分解は，正方行列の特異値分解に帰着される。

従来の特異値分解アルゴリズムとその問題点 密行列 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 = 5000，Xeon 2.8GHz（1～4PU） LAPACK での実行時間（秒） (8/3) n3 二重対角化 密行列Aをまず三重対角行列Tに相似変換してからTの固有値・固有ベクトルを求めるのが最も一般的な計算法。 三重対角化には，後に述べるハウスホルダー法を使う場合がほとんど。 三重対角行列の固有値・固有ベクトルの計算には，色々なアルゴリズムがある。 二重対角行列の 特異値・特異ﾍﾞｸﾄﾙ計算 O(n2) ～ O(n3) 4mn2 逆変換 （左右 m 本ずつの 特異ベクトル） ・二重対角化が実行時間の 　大部分を占める ・速度向上率が低い Aの特異ベクトル {ui }, {vi }

特異値分解アルゴリズムの高速化 ハウスホルダー法による二重対角化 第 k ステップでの処理 鏡像変換 H = I – a wwT による列の消去 Hは対称な直交行列 与えられたベクトルの第１成分以外をゼロにする HA(k) = A(k) – au (utA(k)) 第 k ステップでの処理 ベクトル 左からH を乗算 行列ベクトル積 rank-1更新 A(k) 右からHkR を乗算 左からHkL を乗算 k 非ゼロ要素 ゼロにしたい部分 影響を受ける部分

ハウスホルダー法の特徴と問題点 特徴 問題点 全演算量のほとんどが，BLAS2 である行列ベクトル積と行列のrank-1更新で占められる。 rank-1更新の演算量： 約 (4/3)n3 問題点 BLAS2のため，行列データの再利用性なし キャッシュミス，メモリ競合の影響大 共有メモリ型並列計算機上で高性能が出ない理由

BLAS3 に基づく二重対角化アルゴリズム 基本的なアイディア 二重対角化を２段階で行うことの利点 密行列 A をまず帯幅 L の下三角帯行列 C に変換 次にこの帯行列を下二重対角行列 B に変換 二重対角化を２段階で行うことの利点 下三角帯行列への変換は， BLAS3 のみを使って実行可能 　　　　 キャッシュの有効利用が可能 下三角帯行列から二重対角行列への変換の演算量は O(n2L) 前半部に比べてずっと小さい。 次数 n 下三角 帯行列化 村田法 約 (8/3)n3 O(n2L) A C B 帯幅 L

下三角帯行列化のアルゴリズム ブロック鏡像変換によるブロック列の消去 ブロック鏡像変換 H = I – WαWT 第 K ステップでの処理 ブロックベクトル ブロック鏡像変換によるブロック列の消去 ブロック鏡像変換 H = I – WαWT Hは直交行列 与えられたブロックベクトルを上三角 行列（正確には右上三角部分のみ 非零でそれ以外が零の行列）に変形 HA(k) = A(k) – WαWTA(k) 第 K ステップでの処理 左からH を乗算 行列積 = BLAS3 A(k) 左からHKL を乗算 右からHKR を乗算 非ゼロ要素 ゼロにしたい部分 影響を受ける部分

帯行列の二重対角化（村田法） ・演算量は 8n2L ・並列化も可能 左側からの 直交変換で 更新された 要素 左側からの 直交変換で

本アルゴリズムの長所と短所 長所 短所 BLAS3 の利用により，二重対角化の性能を向上可能 特異ベクトル計算のための計算量・記憶領域が増大 同様のアイディアに基づく三重対角化アルゴリズムでは，高い単体性能・並列性能を確認済み 短所 特異ベクトル計算のための計算量・記憶領域が増大 ２段階の逆変換が必要 詳しくは次のスライドで説明 二重対角化の高速化効果が大きければ，計算量増大を考慮しても全体としては高速化できると予想 特に，求める特異ベクトルが少ない場合は効果が大きいはず。

特異ベクトルの計算手法 二重対角行列の特異ベクトルを計算して２回逆変換 長所 短所 A C B 二重対角行列の特異値・特異ベクトルを求める任意の手法が適用可能 短所 逆変換の演算量が 8mn2 （従来法の2倍）。ただし BLAS3 で実行可能 村田法の変換をすべて記憶するため，n2 の記憶領域が余計に必要 n L 特異値 {σi } QR法 DC法 MR3 I-SVD A C B 4mn2 4mn2 A の特異ﾍﾞｸﾄﾙ {ui }{vi } C の特異ﾍﾞｸﾄﾙ {zi }{wi } B の特異ﾍﾞｸﾄﾙ {xi }{yi }

性能評価

性能評価 評価環境 評価対象・条件 Xeon (2.8GHz), 1～4PU Linux + Intel Fortran ver. 8.1 BLAS:　Intel Math Kernel Library LAPACK:　Intel Math Kernel Library ピーク性能：　5.6GFLOPS/CPU 富士通 PrimePower HPC2500 (2.0GHz), 1～32PU 富士通 Fortran BLAS:　富士通並列化版 BLAS LAPACK:　富士通並列化版 LAPACK ピーク性能：　8GFLOPS/CPU 評価対象・条件 BLAS3 に基づくアルゴリズムと LAPACK の性能を比較 n = 1200 ～ 20000 の乱数行列の特異値分解（全特異ベクトルを計算） BLAS3 アルゴリズムにとっては一番不利な条件 BLAS3 アルゴリズムの L（半帯幅）は各 n ごとに最適値を使用

Xeon での実行時間 プロセッサ数を変えたときの実行時間 結果 BLAS3 アルゴリズムでは PU 数に応じて実行時間が順調に減少 4PU の場合は BLAS3 アルゴリズムが従来法より高速 n = 1200 n = 2500 n = 5000 実行時間（秒） PU数

HPC2500 での実行時間 プロセッサ数を変えたときの実行時間 結果 BLAS3 アルゴリズムは従来法に比べて最大3.5倍高速 プロセッサ数が多いとき加速率が鈍るのは，非並列化部分（ブロック鏡像変換の作成など）の影響と思われる。 n = 5000 n = 10000 n = 20000 実行時間（秒） 3.5倍 PU数

両手法の実行時間の内訳 Xeon，n=5000の場合 考察 BLAS3 アルゴリズムでは，どの部分の実行時間も順調に減少 逆変換1（村田法の逆変換）の占める時間が大きい。 　　　　　この部分について，さらに高速化が必要 必要な特異ベクトルの本数が少ない場合，BLAS3 アルゴリズムはさらに有利

両手法の実行時間の内訳 HPC2500，n=10,000の場合 考察 BLAS3 アルゴリズムでは，どの部分の実行時間も順調に減少 従来法は，二重対角化の部分の加速が鈍い。 ただし，32PUで6倍程度は加速 メモリバンド幅が大きいためと思われる。

まとめと今後の課題

まとめと今後の課題 §まとめ §今後の課題 共有メモリ型並列計算機向けに，BLAS3 に基づく特異値分解アルゴリズムを紹介した Xeon と HPC2500 で評価した結果，PU 数が多い場合は従来法より高い性能が得られた 特に，求める特異ベクトルの本数が少ない場合は効果が大きい §今後の課題 性能の改善 より効率の良い並列化 村田法の逆変換の高速化 より多様なマシン上での性能評価 マルチコアプロセッサ SIMD超並列プロセッサ （Clear Speed，GRAPE-DR，GPU など）

３. SIMD超並列プロセッサによる特異値分解の高速化 ～ 長方行列の特異値分解 ～

長方行列の特異値分解 = ＜応用＞ A ：m×n 密行列 U ：m×n 列直交行列 V ： n×n 直交行列 S ： n×n 対角行列 × 10万 5千 （例） 画像処理 電子状態計算 （Filter Diagonalization Method） 文書検索 （Latent Semantic Indexing） 各種統計計算 （主成分分析，最小二乗法）

高い演算能力を持つSIMD超並列プロセッサ ClearSpeed CSX600 1個の汎用コア + 96個の演算コア 48GFLOPS （倍精度） GRAPE-DR 512個の演算コア 512GFLOPS（単精度） 256GFLOPS（倍精度） GeForce GTX280 （GPU） 240個の演算コア 933GFLOPS（単精度） 90GFLOPS（倍精度） 多数の演算装置の集積により極めて高いFLOPS値 メモリアクセス性能の相対的な低下によるメモリネック

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

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

ClearSpeed 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（行列乗算）で 50GFLOPS CSFFT ライブラリ 今回はこれを利用

CSXLのDGEMMの性能 A B × C += B A × C += m = k = 450 1000 ≦ n ≦ 6000 k = 450 1000 ≦ m = n ≦ 6000 性能（MFLOPS） n n,m A B × C += n k m k n B A × C += m 性能を引き出すにはサイズパラメータのうち少なくとも２個を大きい値にとる必要がある

CSX600向けの　　　　高速化手法

長方行列の特異値分解の手順 A R Q B QR分解 A = QR 二重対角化 R = U1 B V1T 特異値分解 m n Q R B QR分解 A = QR 二重対角化 R = U1 B V1T 特異値分解 B = U2 S V2T 逆変換 R = U’ S V T V = V1 V2 U’ = U1 U2 Qの乗算　 A = U S V T U = QU’

長方行列の特異値分解の手順 QR分解 A = QR 実行時間の大部分を占める 二重対角化 R = U1 B V1T 特異値分解 m ≫ n （例：m =100000, n =5000）の場合 計算量 QR分解 A = QR 2mn2 実行時間の大部分を占める 二重対角化 R = U1 B V1T (8/3)n3 特異値分解 B = U2 S V2T O(n2) ~ O(n3) 逆変換 R = U’ S V T V = V1 V2 U’ = U1 U2 2n3 ~ 4n3 Qの乗算　 A = U S V T U = QU’ 4mn2

高速化の方針 §CSX600を利用する部分 QR分解 A = QR Qの乗算 U = QU’ A = U S V T 行列乗算のみ高速化可能 Qの乗算　 A = U S V T U = QU’ 行列乗算中心のアルゴリズム §CSX600を利用しない（CPUのみで実行する）部分 二重対角化 R = U1 B V1T LAPACKのDGEBRD 逆変換 R = U’ S V T V = V1 V2 U’ = U1 U2 LAPACKのDORMBR 特異値分解 B = U2 S V2T I-SVD（IDBDSLV）

ハウスホルダー変換による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 – × × × × × × 複数のハウスホルダー変換の作用を，行列乗算として まとめて実行可能

ハウスホルダー変換のブロック化 = HL ・・・ H1 （I – Y T YT ） = = （I – Y T YT ） = 複数のハウスホルダー変換の合成 = （I – Y T YT ） 行列乗算としてまとめて作用 CSX６００で高速化

（Ⅰ）ブロックQR分解（LAPACKで使用） I – Y1T1Y1T I – Y2T2Y2T I – Y3T3Y3T 分解 更新 H1 H2 H3 H3H2H1= I – Y1T1Y1T ・・・行列乗算 ブロックの分解（　　）は行列・ベクトル積 演算量： L (ブロック幅) << n のとき 2mn2 CSX600を使うにはL≧450 ⇒行列・ベクトル積の演算量の増加

（Ⅱ）再帰的QR分解 ・・・行列乗算 ほぼ全ての演算が行列乗算 演算量： 3mn2 Qの乗算の効率が良い I – Y2T2Y2T 更新 合成 (I – Y1T1Y1T) ×(I – Y2T2Y2T) = I – Y T YT I – Y1T1Y1T ・・・行列乗算 I – Y11T11Y11T ほぼ全ての演算が行列乗算 演算量： 3mn2 Qの乗算の効率が良い

（Ⅲ）再帰的QR分解の拡張 L ・・・行列乗算 ほぼ全ての演算が行列乗算 演算量：Ⅰ(2mn2)とⅡ(3mn2)の中間 I – Y1T1Y1T I – Y2T2Y2T I – Y3T3Y3T 分解 更新 I – Y1T1Y1T ・・・行列乗算 ほぼ全ての演算が行列乗算 演算量：Ⅰ(2mn2)とⅡ(3mn2)の中間 ブロック幅Lを適切に選ぶことで，他の２方法より高速になる可能性がある

Q の乗算 （Ⅰ）ブロックQR分解・（Ⅲ）再帰的QR分解の拡張の場合 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，メモリ8GB Intel Fortran -O3 + Intel Math Kernel Library ClearSpeed Advance ボード 評価例題 [-0.5,0.5] の一様乱数を要素とする m×n 乱数行列の特異値分解 10000 ≦ m ≦ 100000，1000 ≦ n ≦ 4000 評価項目 ClearSpeed ボードでの3種のQR分解法の性能比較 特異値分解全体の ClearSpeed ボードによる高速化効果 精度評価

Ⅲ．再帰的QR分解の拡張（数字はブロック幅） CSX600使用時の各手法の性能比較（１） m =100000 n =4000 実行時間（秒） Ⅱ．再帰的　QR分解 Ⅰ．ブロックQR分解 Ⅲ．再帰的QR分解の拡張（数字はブロック幅）

Ⅲ．再帰的QR分解の拡張（数字はブロック幅） CSX600使用時の各手法の性能比較（２） m =10000 n =4000 実行時間（秒） Ⅰ．ブロックQR分解 Ⅱ．再帰的　QR分解 Ⅲ．再帰的QR分解の拡張（数字はブロック幅）

CSX600による特異値分解全体の高速化 m = 10000 n=1000 (m:n = 10:1) 実行時間（秒） 1.3倍 1.2倍 1.8倍 3.1倍 LAPACKのDGESDD　　　　　 再帰的QR分解を利用 LAPACKのDGESDD　　　　　 再帰的QR分解を利用

高速化率＝PC単体での実行時間／PC＋CSX600での実行時間 行列サイズ による高速化率の変化（１） §再帰的QR分解を用いた場合 高速化率 m n 高速化率＝PC単体での実行時間／PC＋CSX600での実行時間

各部分の高速化率の変化 m = 50000 m = 100000 高速化率 n n

高速化率＝PC単体での実行時間／PC＋CSX600での実行時間 行列サイズ による高速化率の変化（２） §LAPACKのDGESDDを用いた場合 高速化率 m n 高速化率＝PC単体での実行時間／PC＋CSX600での実行時間

精度評価 ∥US VT – A∥F ∥UTU – I∥F 左特異ベクトルの直交性 残差 m ： n ： m ： n ： 1000 2000 3000 m ： n ： 1000 2000 3000

まとめと今後の課題

まとめと今後の課題 §まとめ CSX６００による長方行列特異値分解の高速化を行った CSX６００に適したQR分解アルゴリズムの選択を行った §今後の課題 よりCSX600に適したQR分解アルゴリズムの開発 CSX600による Rの特異値分解の高速化 他の行列計算へのCSX６００の適用