分散メモリ型スパースソルバの開発と評価

密行列固有値解法の最近の発展（II）
ーマルチシフトQR法ー
名古屋大学計算理工学専攻
山本有作
2006年1月28日
LA研究会
目次
１. はじめに
２. 非対称行列向けのマルチシフトQR法
３. 対称行列向けのマルチシフトQR法
１. はじめに
• 本研究で対象とする問題
– 標準固有値問題 Ax = lx
• A ： n×n 密行列（実対称または実非対称）
• 応用分野
– 対称行列
• 分子計算（分子軌道法，FMO-MO法）
• 統計計算
– 非対称行列
• MHD，化学工学，量子化学，制御理論
– Cf. Bai and Demmel: A test matrix collection for non-Hermitian
eigenvalue problems.
固有値計算の流れ
密行列 A
三重対角化/Hessenberg化
三重対角行列 T
/ Hessenberg行列 H
計算内容
計算手法
QTAQ = T
(Q: 直交行列)
ハウスホルダー法
固有値の計算
A の固有値 {li }
固有ベクトルの計算
Tui =li ui
QR法高い信頼性
分割統治法
MR3アルゴリズム
I-SVDアルゴリズム
T or Hの固有ﾍﾞｸﾄﾙ {ui }
逆変換
Aの固有ベクトル {vi }
vi = Q ui
逆変換
各部分の実行時間（非対称行列）
• 全固有値を求める場合の演算量
– Hessenberg 化：
– QR法：
n3
(10/3)
10n3 （経験値）
Execution time (min)
45
QR
40
Hessenberg reduction
35
30
25
– 演算量（時間）の大部分をQR法が占める。
– QR法の高速化が必要
20
15
10
5
O
Origin2000 1PU （R12000，400MHz）上での実行時間
K. Braman et al: “The Multishift QR Algorithm I” より抜粋
Hessenberg 化の時間は推定値
LM
10
00
TU
B1
00
0
RD
B1
25
TO
0
LS
20
0
PD 0
E2
9
M
HD 61
32
00
B
0
各部分の実行時間（対称行列）
• 全固有値を求める場合の演算量
– 三重対角化：
– QR法：
n = 100,000の行列の固有値計算
(4/3) n3
O(n2)
5555 min.
– QR法の演算量は三重対角化に比べ小さい
– しかし，三重対角化部分は並列化が容易
– QR法も並列化できることが望ましい
8 min.
10 min.
QR法は PrimePower HPC2500 上での実測値
三重対角化の時間は推定値
（1CPUでの性能4GFLOPS，並列化効率70%と仮定）
10 min.
1CPU
三重
対角化
QR法
三重対角化部分
のみ1000CPU
QR法の特徴と高性能計算
高性能計算の条件
QR法の特徴
•計算の逐次性
（bulge-chasing 型演算）
•計算の並列性
•低いデータ再利用性
•高いデータ再利用性
（キャッシュの有効利用）
オリジナルのQR法のままでは高性能化が困難
高性能計算の条件を満たす新しいアルゴリズムが必要
マルチシフトQR法
２. 非対称行列向けのマルチシフトQR法
• QR法（Francis, 1961 & Kublanovskaya, 1961）
– 行列 A0 から始めて次のように QR 分解と相似変換を繰り返す。
• A0 = Q1R1
• A1 = R1Q1 (= Q1–1A0 Q1 )
• A1 = Q2R2
• A2 = R2Q2 (= Q2–1A1 Q2 = Q2–1Q1–1A0 Q1Q2 )
• QR法の収束
– 適当な条件の下で，Ak は（ブロック）上三角行列に収束
– A0 の固有値を絶対値の大きい順に l1, l2, … , ln とすると，対角要素 aij
は li に１次収束
– 非対角要素 aij （j < i）は収束率 rij ≡ |li| / |lj| で 0 に１次収束
ダブルシフトQR法
• シフトの導入
– Ak から固有値 li の近似値 s を引いた行列に対して QR 法を適用
• Ak – s I = Qk Rk
• Ak+1 = Rk Qk + s I (= Qk–1Ak Qk )
rij ≡ |li – s| / |lj – s| ≒ 0 より収束が加速
• ダブルシフトQR法
– 共役複素数の固有値ペアに対し，シフト s, s による2反復をまとめて行う
• (Ak – s I)(Ak – s I) = Qk Rk
• Ak+2 = Qk–1Ak Qk
– シフトは右下隅の 2×2 行列の固有値に取る
→ 局所的に2次収束
– 複素数の固有値を持つ場合でも，実数演算だけで QR 法を実行可能
Hessenberg 形と Implicit Q 定理
• Hessenberg 形の利用
– Ak が Hessenberg 形のとき，
• QR法の1ステップは O(n2) で実行可能
• Ak+1 も Hessenberg 形
• Implicit Q 定理の利用により，QR 分解を陽に行う
ことなく1ステップの計算を実行可能
0
Hessenberg 形
aij = 0 if i > j+1
A0 を Hessenberg 形に相似変換してからQR法を適用
• Implicit Q 定理
– U，V が直交行列で，G = UtAU， H = VtAV が共に既約な Hessenberg 行
列であるとする。このとき，もし U と V の第1列が等しいならば，±1 を
要素に持つ対角行列 D が存在して，V = UD，H = DGD が成り立つ
– すなわち，行列 A を既約な Hessenberg 行列に相似変換する直交行列
は，第1列目だけが与えられれば（実質的に）一意に定まる
Implicit shift QR法
• 1ステップの計算
• (Ak – s2I)(Ak – s1I) の第1列を e1の定数倍にするハウスホル
ダー変換H0 を求める
• Ak’ = H0tAk H0
• 直交行列による相似変換を繰り返すことにより， Ak’ を再び
Hessenberg 行列に変形する（bulge-chasing）
Ak
H0 AkH0
0
0
0
• この変換の性質
– ある直交行列 H による相似変換
– H の第1列は，Qk の第1列に等しい
• 上記第3のステップは第2行・第2列以降のみに影響
– H は Ak を Hessenberg 形に変換
Implicit Q 定理より，Ak+1 = H tAk H はQR法の1ステップと等価
0
Ak+1
0
演算パターンとその特徴
• Bulge-chasing の演算パターン
– 3×3のハウスホルダー変換を左右からかけることにより，bulge を1つ
右下に動かす
左から Hl を乗算
0
0にしたい要素
右から Hl を乗算
0
0
更新される要素
• 演算パターンの特徴
– 並列性は O(n)
• 並列アルゴリズム： R. Suda et al.（1999），G. Henry et al. (2002)
– QR法の1ステップで，各行列要素は3回のみ更新
データ再利用性が低く，キャッシュの有効利用が困難
マルチシフトQR法
• 基本的なアイディア（Bai & Demmel, 1989）
– m 個のシフトに対するQR法のステップを同時に行う
• (Ak – smI) ・・・ (Ak – s2I)(Ak – s1I) = Qk Rk
• Ak+m = Qk–1Ak Qk
– シフト s1, s2, …, sm は，の右下隅の m×m 行列の固有値に取る
• マルチシフトQR法の収束（Watkins & Elsner, 1991）
– 上記のようにシフトを取ったとき，マルチシフトQR法は局所的に2次収束
– 特に，Q1, Qm+1, Q2m+1, … , Qk の積を Q(k) とするとき，Q(k) の最後の m 列
の張る部分空間は，A1のある不変部分空間に2次収束する
マルチシフトQR法の演算パターン
• 1ステップの計算（implicit 型）
• (Ak – s2I) ・・・ (Ak – s2I)(Ak – s1I) の第1列を e1の定数倍にするハウスホル
ダー変換H0 を求める
• Ak’ = H0tAk H0
• 直交行列による相似変換を繰り返すことにより， Ak’ を再びHessenberg
行列に変形する（bulge-chasing）
•
演算パターンの特徴
– 並列粒度は m/2 倍増加
– QR法の1ステップで，各行列要素は m+1回更新
– データ再利用性の向上により，キャッシュの有効利用が可能
（WY representation を用いた BLAS3化）
マルチシフトQR法の実際の性能
• m の取り方に関する問題点
– m が 6 程度を超えると，収束までの反復回数が急速に増大
– 一方，BLAS3 で性能を出すには m を数十程度に取ることが必要
このマルチシフトQR法では，計算機性能を引き出せない
K. Braman et al.: “The Multishift QR Algorithm I” (2002) より引用
収束性悪化の原因は？
収束性悪化の解析
• シフトの伝播に関する理論的解析（Watkins, 1996）
– bulge を形作る (m+1)×(m+1) 行列を Bl とするとき，シフト
s1, s2, …, sm は，行列束 (Bk, N) の m 個の有限固有値
• ただし N は固有値 0 の m+1 次ジョルダン細胞
0
Bl
• 実験からの考察
– (Bl, N) の固有値は，m が大きいとき，摂動に対して極めて
敏感になる
m = 4 の場合
Bulge chasing の過程でシフトが (Bk, N) の固有値として伝わる際に，
丸め誤差による摂動でシフトが劣化し，2次収束を阻害するのでは
ないか？
多数の小さい bulge を使ったマルチシフトQR法
• 基本的なアイディア（K. Braman et al., 2002）
– 従来法と同様にしてシフト s1, s2, …, sm を計算
– シフトを2つずつ組にしてダブルシフトQR法を m/2 回実行
• この結果得られる Ak+1 は，無限精度演算では従来法と同じ
– m/2 個の bulge をできるだけ密に詰めて chasing を行うこ
とで，データ参照の局所性を向上
• 3行離れていれば，互いの干渉なしに計算が可能
大きな bulge を避けることで，丸め誤差の影響を抑えつつデータの
再利用性を高め，性能を向上
0
m = 6 の場合
収束性に関する実験結果
• シフトの数と総演算量との関係
– DHSEQR：従来のマルチシフトQR法（LAPACK）
– TTQR：新手法
K. Braman et al.: “The Multishift QR Algorithm I” (2002) より引用
新手法では，m を数十まで増やしても収束性が低下しない
Level-3 BLAS の利用
• 更新処理の分割
– m/2 個の bulge をそれぞれ r 行 chasing する際，
まず対角ブロックのみを更新
– 非対角ブロックは，更新に使ったハウスホル
ダー変換を1個の行列に蓄積し，後でまとめて
更新
非対角ブロックの更新で level-3 BLAS が使用可能
0
Bulge（3×3）
最初に更新
まとめてBLAS3で更新
• r の決め方
– 演算量の点から，r ～ 3m とするのが最適
– このとき，演算量は従来のマルチシフトQR法
の約2倍（非ゼロ構造を利用した場合）
蓄積されたハウス
ホルダー変換の非
ゼロ構造
従来のマルチシフトQR法との組み合わせ
• 基本的なアイディア（Kressner, 2005）
– 3×3の bulge を m/2 個同時に chasing する代
わりに，(q+1)×(q+1) の bulge を m/q 個同時に
chasing する
– q ～ 6 であれば，収束性は悪化しない
データ再利用性をさらに向上
0
Bulge（(q+1)×(q+1)）
最初に更新
まとめてBLAS3で更新
各手法の性能
• 比較する対象の手法
– DHSEQR：従来のマルチシフトQR法（LAPACK）
– MTTQR：多数の小さいシフトを使ったマルチシフトQR法
（ns： bulge 1個あたりのシフト数）
D. Kressner: “On the Use of Larger Bulges in the QR Algorithm”より引用
ほとんどの場合，ns = 4 or 6 とした新手法が最高速
共有メモリ向けの並列化手法
• Level-3 BLAS 内部での並列化
– 非対角ブロックの更新において，並列化 BLAS
を利用することにより，容易に並列化可能
• より効率的な並列化
– 非対角ブロックの更新において，次の対角ブロッ
クの更新に必要な部分のみを先に更新
0
Bulge（3×3）
最初に更新
並列化困難な対角ブロックの更新を，非対角ブ
ロックの更新と並列に実行可能
まとめてBLAS3で更新
次の対角ブロックの
更新で必要
シフト数の最適決定
• m の最適化
– m を大きくすると，局所性は向上するが，シフト
の計算量が増加
– 最適な m の値は計算機と問題サイズに依存
• r （1回の bulge chasing の行数）の最適化
– 演算量最小化のためには r ～ 3m が最適
– BLAS3 の実行時間は必ずしも演算量に比例し
ないため，実行時間の最適値はこれとは異なる
m，r に対する自動チューニングの必要性
n
m の最適値
1000～1999
60
2000～2499
116
2500～3999
150
4000～
180
実験的に求めた m の最適値
（Origin2000上，Braman et al.
より引用）

Download Report