1 / 62

応用数理工学特論　第 10 回

応用数理工学特論　第 10 回. 2008 年 7 月 3 日 計算理工学専攻 山本有作. 目次. １ . 　はじめに ２ . 　共有メモリ型並列計算機での特異値分解の高速化 ～ 正方行列の特異値分解 ～ ３ . SIMD 超並列プロセッサでの特異値分解の高速化 ～ 長方行列の特異値分解 ～. １ . 　はじめに. n. n. 特異値分解とは. 任意の実正方行列 A は，直交行列 U ，対角行列 S ，直交行列 V T の積に分解できる。 参考 A が実対称行列の場合， A=UDU T A が任意の実行列の場合， A=SDS – 1

reece
Télécharger la présentation

応用数理工学特論　第 10 回

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. 応用数理工学特論　第10回 2008年7月3日 計算理工学専攻 山本有作

  2. 目次 １.　はじめに ２.　共有メモリ型並列計算機での特異値分解の高速化 ～ 正方行列の特異値分解 ～ ３.SIMD超並列プロセッサでの特異値分解の高速化 ～ 長方行列の特異値分解 ～

  3. １.　はじめに

  4. n n 特異値分解とは • 任意の実正方行列 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

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

  6. 固有値分解（対角化）との関係 • 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の固有ベクトル

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

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

  9. 対象とする計算機

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

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

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

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

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

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

  16. 従来の特異値分解アルゴリズムとその問題点 計算内容 計算手法 密行列 A U0TAV0 = B (U0, V0: 直交行列) ハウスホルダー法 二重対角化 二重対角行列 B QR法 分割統治法 MR3アルゴリズム I-SVDアルゴリズム Bvi =σi xi BTxi =σi yi 二重対角行列の 特異値・特異ﾍﾞｸﾄﾙ計算 Bの特異値 {σi}， 　 特異ベクトル {xi }{yi } vi = V0 yi ui = U0 xi 逆変換 逆変換 Aの特異ベクトル {ui }, {vi }

  17. 実行時間（全特異ﾍﾞｸﾄﾙ） n = 5000，Xeon 2.8GHz（1～4PU） LAPACK での実行時間（秒） 各部分の演算量と実行時間 演算量 密行列 A (8/3) n3 二重対角化 O(n2) ～ O(n3) 二重対角行列の 特異値・特異ﾍﾞｸﾄﾙ計算 4mn2 逆変換 （左右 m 本ずつの 特異ベクトル） ・二重対角化が実行時間の 　大部分を占める ・速度向上率が低い Aの特異ベクトル {ui }, {vi }

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

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

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

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

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

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

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

  25. 性能評価

  26. 性能評価 • 評価環境 • 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ごとに最適値を使用

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

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

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

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

  31. まとめと今後の課題

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

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

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

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

  36. 本章の目的 • 長方行列の特異値分解アルゴリズムをCSX600を利用して高速化する • 既存のアルゴリズムをそのまま使用せず，行列積をなるべく効率的に使う形で CSX600 向けに最適化する • 性能評価を行い，更なる高速化に向けての課題を明らかにする

  37. ClearSpeed CSX600 について

  38. CSX600のアーキテクチャと性能 • CSX600 チップ • 1個の主プロセッサ • 96個の演算専用プロセッサ • 64ビット • 倍精度2演算/サイクル • 128B レジスタファイル • 6KB SRAM • 250MHz で動作 • ピーク性能： 48GFLOPS • ClearSpeed Advance ボード • 2個の CSX600 プロセッサ • 1GB DRAM • PCI-X バスにより PC 本体と接続 • ピーク性能： 96GFLOPS

  39. 今回はこれを利用 CSX600の利用環境 • Software Development Kit • コンパイラ：　Cn 言語によるチップ内並列プログラミング • デバッガ • シミュレータ • CSXL ライブラリ • ClearSpeed Advance ボード用の BLAS • 行列データを PC の主メモリ上に置いてコール • PC ⇔ ボード間の転送はライブラリ内で実行 • 公称性能： DGEMM（行列乗算）で 50GFLOPS • CSFFT ライブラリ

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

  41. CSX600向けの　　　　高速化手法

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

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

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

  45. level-2 BLAS CSXLを使えない ハウスホルダー変換によるQR分解 ハウスホルダー変換による列の消去 H1 A = ( I – t1y1y1T ) A = A(1) 上三角化とQR分解 Hn・・・ H2 H1 A = A(n) A = H1 H2・・・ Hn A(n) = QR ・・・ A A(1) A(2) A(n) = R

  46. 複数のハウスホルダー変換の合成 Compact WY representation Hn・・・ H2 H1 = ( I – tn yn ynT ) ・・・ ( I – t2y2y2T )( I – t1y1y1T ) = I – Yn Tn YnT ただし， Yn = [ y1 | y2 | ・・・ | yn]　（m×n行列） Tn：　n×n下三角行列 I – I – I – ・・・ = × × × × × × 複数のハウスホルダー変換の作用を，行列乗算として まとめて実行可能

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

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

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

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

More Related