櫻井・杉浦法（Sakurai-Sugiura method）

演習III 論文紹介
内容に問題や不備がある場合は、お手数ですが hellorusk1998 [at] gmail.com までご連絡お願いいたします。

Transcript

1. 櫻井‧杉浦法 演習 III 発表 杉本 海人 2020/06/24 1 / 46
2. 読んだ論文 Sakurai, Tetsuya, and Hiroshi Sugiura. "A projection method for generalized eigenvalue problems using numerical integration." Journal of computational and applied mathematics 159.1 (2003): 119-128. 2 / 46
3. 概要 与えられた n × n の複素行列 A, B に対して Ax = Bx ⇔ (A − B)x = 0 を満たすような（一般）固有値 を求める問題（一般固有値問題） A, B はサイズの大きな疎行列で, 全部を求めたいのではなく, 一部だ け求まればいいというような状況 3 / 46
4. 概要 B が逆行列を持てば B−1Ax = x となり, 行列 B−1A の固有値問題に帰着される とはいえ, 数値計算上はこのようにしてもあまり意味がない • B−1A は依然としてサイズが大きい • A, B が対称‧疎など, 数値計算上良い性質を持っていたとしても, B−1A でそれが保持されているとは限らない →櫻井‧杉浦法では, 元の問題を, サイズの小さな行列の一般固有値問 題に帰着させる（その解から元の一般固有値問題の解の一部が求 まる） 4 / 46
5. 流れ 一般固有値問題 (A − B)x = 0 の解法 1 f (z) = u∗(zB − A)−1v (u, v ∈ Cn, u, v ≠ 0) という関数を考えると, f (z) = d j=1 j z − j + g(z) という形に変形できることを示す 2 閉路 Γ に対して, 周回積分 k = 1 2i ∫ Γ (z − )kf (z)dz を要素とす るような行列 Hm, H< m を考えると, 一般固有値問題 (H< m − Hm )x = 0 の解が 元の問題の閉路 Γ 内における解を与える ものになっていることを示す 5 / 46
6. 用語の定義 • (A − B)x = 0: 一般固有値問題 • A − zB: 行列束 • 行列束の固有値を以下のように定義する • p(z) = det(A − zB) の根 これは一般固有値問題 (A − B)x = 0 の解を与えるものになっている • 無限大 ∞ (重複度: n − deg(p(z))) 便宜上の解 • A や B が逆行列を持つ場合は 通常の固有値問題に帰着でき, 通常 の固有値問題では deg(p(z)) = n なので, 便宜上の 無限大 ∞ の解 は存在しない 6 / 46
7. 例 A = 1 0 0 0 1 0 0 0 0 , B = 2 0 0 0 0 0 0 0 1 の場合 p(z) = det(A − zB) = (2z − 1)(−z) なので, 行列束 A − zB の固有値は 0, 1 2 , ∞ となる 7 / 46
8. ジョルダン標準形 任意の n × n の行列 A に対して, その相異なる固有値を 1, ..., k とす ると, P−1AP = J となるような P, J がある. ただし, J は 1, ..., k のジョルダンブロック を並べたもの. ジョルダンブロックとは J = 1 0 · · · 0 0 1 0 0 0 0 . . . ... 1 0 0 0 という形をしたもの. 特に固有値が全て異なる場合, J は対角行列 8 / 46
9. ワイエルシュトラスの標準形 行列束 A − zB のジョルダン標準形のようなもの. n × n の 行列束 A − zB の相異なる固有値を 1, 2, ..., k, ∞ とし, p() の次数を d（すなわち ∞ の重複度を n − d）とすると, P(zB − A)Q = zId − Jd O O zJn−d − In−d となるような正則の複素行列 P, Q が存在する. ただし, Jd, Jn−d はそれぞれジョルダン標準形（特に k = d, すなわち固 有値に重複度が無い場合 Jd は対角行列. 今回はこのケースのみで考 える）. Jn−d の対角成分は 0. 9 / 46
10. これを用いると Q−1(zB − A)−1P−1 = zId − Jd O O zJn−d − In−d −1 (zB − A)−1 = Q zId − Jd O O zJn−d − In−d −1 P P = P1 P2 (P1 ∈ Cd×n, P2 ∈ C(n−d)×n), Q = Q1 Q2 (Q1 ∈ Cn×d, Q2 ∈ Cn×(n−d)) とおく. 10 / 46
11. ① f (z) = u∗(zB − A)−1v (u, v ∈ Cn, u, v ≠ 0) という関数を考えると, f (z) = d j=1 j z − j + g(z) という形に変形できることを示す f (z) = u∗(zB − A)−1v = u∗Q zId − Jd O O zJn−d − In−d −1 Pv = u∗ Q1 Q2 zId − Jd O O zJn−d − In−d −1 P1 P2 v = u∗Q1 (zId − Jd )−1P1v − u∗Q2 (In−d − zJn−d )−1P2v 11 / 46
12. 第1項について Id は単位行列, また Jd は固有値 1, ...d が対角要素の対角行列なので, (zId − Jd )−1 は (z − i )−1 = 1 z − i が対角要素の対角行列. P∗ = (p1...pn ), Q = (q1...qn ) とおくと u∗Q1 (zId − Jd )−1P1v = u∗q1 · · · u∗qd (zId − Jd )−1 p∗ 1 v . . . p∗ d v = d j=1 u∗qjp∗ j v z − i = d j=1 j z − i 12 / 46
13. 第2項について Jn−d は対角成分が 0 のジョルダン標準形であった. ここで, サイズ k × k の D = 0 1 0 · · · 0 0 0 1 0 0 0 0 0 . . . ... 1 0 0 0 0 に対し, (I − zD)−1 = I + zD + (zD)2 + ... + (zD)k−1 が成立する. この行列は対角成分が 0 のジョルダン標準形のジョルダンブロックで ある. 13 / 46
14. 第2項について したがって, (In−d − zJn−d ) のジョルダンブロックの最大サイズを K と すると, 第 2 項 −u∗Q2 (In−d − zJn−d )−1P2v は次数 K − 1 の z の多項式に なる（これを g(z) とおく）. 以上より, ① f (z) = u∗(zB − A)−1v (u, v ∈ Cn, u, v ≠ 0) という関数を考えると, f (z) = d j=1 j z − j + g(z) という形に変形できることを示す が完了. 14 / 46
15. ② 閉路 Γ に対して, 周回積分 k = 1 2i ∫ Γ (z − )kf (z)dz を要素とするよ うな行列 Hm, H< m を考えると, 一般固有値問題 (H< m − Hm )x = 0 の解が 元の問題の閉路 Γ 内における解を与えるものになっていることを示す 複素平面上で単一閉曲線（自身と交わらない閉曲線, ジョルダン曲線 とも）Γ を考える. Γ 内にある固有値について考えていく. ※ 数値計算上は Γ は を中心とする円周で考えるが, 今はいったん任意の閉曲線で考えておく 15 / 46
16. 複素解析 • 複素関数 f (z) が D で正則であるとは, 全ての z ∈ D において, f が 微分可能かつ f′ が連続であること. • 例) f (z) = 1 z は z = 0 を除いて C 全体で正則 • 複素積分: ∫ C f (z)dz は複素平面上の曲線上 C で f (z) の線積分を 考えたもの. • コーシーの積分定理: 単連結領域（ドーナツみたいな形をしていない領域）D において 正則な関数 f (z) と D 内の区分的に滑らかな閉曲線 C に対して, ∫ C f (z)dz = 0 16 / 46
17. 複素解析 • 孤立特異点: 複素関数 f (z) で正則にならない点 • ローラン展開: f (z) の孤立特異点を a としたとき, f (z) = ∞ n=−∞ cn (z − a)n という形で展開できる（テイラー展開を負の次数方向にも広げた もの）. • 特に, 負の次数の部分, すなわち ∞ n=1 c−n (z − a)−n を主要部と呼ぶ. 17 / 46
18. 複素解析 • 主要部について, S = {n ∈ N|c−n ≠ 0} を考えると, • S が空でない有限集合のとき, a は極と呼ばれる. また, n の最大 値を位数と呼ぶ. • f (z) = g(z) z − a (g(z) は z = a で正則) のような形でかけるとき, a は f (z) の 1 位の極になっている • （ S が空集合のとき a は除去可能な特異点と呼ばれる, S が無限 集合のとき a は真性特異点と呼ばれる） 18 / 46
19. 複素解析 • また, 1 z − a の係数 c−1 を, 孤立特異点 a における f (z) の留数と 呼ぶ. • a が f (z) の 1 位の極の場合は, わざわざローラン展開して係数を 考えなくても, lim z→a {(z − a)f (z)} で留数が求まる. • f (z) = g(z) z − a の形のとき, f (z) の a における留数は g(a) となる. 19 / 46
20. 複素解析 • 留数定理 • C を区分的になめらかな閉曲線とし, C で囲まれた領域を D と する. D 内に有限個の点 a1, ..., am があり, f (z) は C および D を含むあ る開集合で, a1, ..., am を除き, 正則であるとする. このとき ∫ C f (z)dz = 2i m k=1 Res[f (z), ak ] ただし Res[f (z), ak ] は f (z) の ak における留数 20 / 46
21. 本題に戻ると... f (z) = d j=1 j z − j + g(z) に対し, k = 1 2i ∫ Γ (z − )kf (z)dz を考えると, k = 1 2i d j=1 ∫ Γ j z − j (z − )kdz + ∫ Γ g(z)(z − )kdz 第 2 項について, 被積分関数が領域全体で正則なので, 0 になる（コー シーの積分定理）. 第 1 項について, 領域の内側にある固有値 1, ...m についてのみ正則 でない. 21 / 46
22. j z − j (z − )k の z = j における留数は j (j − )k なので, 留数定理 より k = 1 2i · 2i m j=1 j (j − )k + 0 = m j=1 j (j − )k （ただし 1, ...m は 閉曲線 Γ 内の固有値） 22 / 46
23. この k = 1 2i ∫ Γ (z − )kf (z)dz = m j=1 j (j − )k に対して, 次のような 2 種類の m × m の行列を考える. Hm = 0 1 2 · · · m−1 1 2 3 m 2 3 4 m+1 . . . ... m−1 m m+1 2m−2 , H< m = 1 2 3 · · · m 2 3 4 m+1 3 4 5 m+2 . . . ... m m+1 m+2 2m−1 23 / 46
24. これらの 2 つの行列は, ヴァンデルモンド行列 Vm = 1 1 1 1 (1 − ) (2 − ) (3 − ) (m − ) (1 − )2 (2 − )2 (3 − )2 (m − )2 . . . ... (1 − )m−1 (2 − )m−1 (3 − )m−1 (m − )m−1 を使って, 次のように書き換えられる. HM = VmDmVT m H< M = VmDm ΛVT m ただし, Dm = diag(1, ..., m ), Λ = diag(1 − , ..., m − ) 24 / 46
25. すると, 行列束 H< m − zHm について, H< m − zHm = VmDm ΛVT m − zVmDmVT m = VmDm (Λ − zI)VT m と変形できる. 固有値 1, ...m は全て異なるという仮定から, ヴァン デルモンド行列 Vm は正則. また, j = u∗qjp∗ j v ≠ 0 より Dm も正則. 結局, Λ − zI は H< m − zHm のワイエルシュトラスの標準形になっている. 25 / 46
26. したがって, 行列束 H< m − zHm の相異なる固有値は 1 − , ..., m − に なっている. 一般固有値問題 (H< m − Hm )x = 0 の解 1 − , ..., m − を求めて, それ ぞれに を足せば, 元の問題の Γ 内の解になっていると言える. ② 閉路 Γ に対して, 周回積分 k = 1 2i ∫ Γ (z − )kf (z)dz を要素とするよ うな行列 Hm, H< m を考えると, 一般固有値問題 (H< m − Hm )x = 0 の解が 元の問題の閉路 Γ 内における解を与えるものになっていることを示す が完了. 26 / 46
27. 残る問題 H< m , Hm の各要素である, 周回積分 k = 1 2i ∫ Γ (z − )kf (z)dz はどう やって求める？ → Γ は を中心とする円周で考える. すなわち, を円の半径として, j = + e2i N j (0 ≤ j ≤ N − 1) での値 を使って近似する（詳細は次ページ） 27 / 46
28. k = 1 2i ∫ Γ (z − )kf (z)dz は z = + ei で z から に座標変換すると, 1 2i ∫ Γ (( + ei) − )kf ( + ei)ieid = 1 2 ∫ Γ (( + ei) − )k+1f ( + ei)d これに台形則を適用すると, 1 2 N−1 i=0 (( + e2i N j) − )k+1f ( + e2i N j) 2 N = 1 N N−1 i=0 (j − )k+1f (j ) 28 / 46
29. 一般固有値問題 (A − B)x = 0 に対する櫻井‧杉浦法のアルゴリズム 入力: u, v ∈ Cn, N, m, , ① j ← + e2i N j (0 ≤ j ≤ N − 1) ② yj ← (jB − A)−1v (0 ≤ j ≤ N − 1) ③ fj ← u∗yj (0 ≤ j ≤ N − 1) ④ ˆ k = 1 N N−1 i=0 (j − )k+1fj により ˆ k を 0 ≤ k ≤ 2m − 1 の範囲で計算 （次ページに続く） 29 / 46
30. 一般固有値問題 (A − B)x = 0 に対する櫻井‧杉浦法のアルゴリズム ⑤ 先ほど求めた ˆ k から ˆ H< m と ˆ Hm を考えて, 一般固有値問題 (ˆ H< m − ˆ Hm )x = 0 を解き, 解 1, ..., m を求める. ⑥ ˆ j ← + j (1 ≤ j ≤ m) 出力 ˆ 1, ..., ˆ m 30 / 46
31. 並列化 A, B がサイズの大きな行列の場合, 計算に最も時間がかかるのは, f (j ) = u∗(jB − A)−1v を各 j に対して求める部分 （前述アルゴリズムの②と③の部分） →この計算を並列計算することで, より高速に解を求められるように なる. 31 / 46
32. アルゴリズムの入力値の決定 櫻井‧杉浦法では, 積分経路 Γ の中心 , 半径 と, その経路の内部に 含まれる固有値の個数 m を, 入力として適切に決めねばならない. 元論文ではこの入力値の取り方に関しては説明がされていないが, 経路内の固有値の個数を決定する手法を議論した論文を引用で紹介し ている. 32 / 46
33. 固有ベクトルの求め方 櫻井‧杉浦法のステップ 1 では, 以下のようなワイエルシュトラスの 標準形を考えた（P, Q は正則）. P(zB − A)Q = zId − Jd O O zJn−d − In−d そして, P∗ = (p1...pn ), Q = (q1...qn ) とおいて, f (z) = u∗(zB − A)−1v = d j=1 u∗qjp∗ j v z − i + g(z) = d j=1 j z − i + g(z) という式変形を考えた（ただし j = u∗qjp∗ j v とおいた）. 33 / 46
34. この Q の k (1 ≤ k ≤ d) 番目の列ベクトル qk が, k 番目の固有値 k の 固有ベクトルに対応している. （略証） P(zB − A)(q1...qn ) = zId − Jd O O zJn−d − In−d である. したがって, P(zB − A)qk は この行列の k 番目の列ベクトルを 取ることに相当するから, P(zB − A)qk = 0 · · · 0 z − k 0 · · · 0 T (zB − A)qk = P−1 0 · · · 0 z − k 0 · · · 0 T となる. 34 / 46
35. z = k を代入すると, (kB − A)qk = 0 となる. これは (k, qk ) が Ax = Bx を満たす (, x) のペアであること を意味する. Q が正則であることから qk ≠ 0 である点にも注意して, qk が k に対 応する固有ベクトルである. 何らかの方法でこの qk を求めたい. 35 / 46
36. 櫻井‧杉浦法のステップ 2 では, f (z) = u∗(zB − A)−1v に対して, k = 1 2i ∫ Γ (z − )kf (z)dz という周回積分を考えた. ここで, sk = 1 2i ∫ Γ (z − )k (zB − A)−1v dz というベクトルを考える. k と sk は k = u∗sk という関係がある. また, Γ を円周上の路とすれば, sk は k と同様, 分点を取って近似を 計算することが可能である. つまり ˆ sk = 1 N N−1 i=0 (j − )k+1(jB − A)−1v 36 / 46
37. この sk を活用して qk を計算したい. 先週の内容で述べた通り, 留数定理により, Γ の内部にある固有値を 1, ..., m として k = 1 2i ∫ Γ (z − )kf (z)dz = m j=1 j (j − )k = m j=1 u∗qjp∗ j v(j − )k と書けるのであった. k = u∗sk より, sk = m j=1 qjp∗ j v(j − )k となる. ここで k = p∗ k v とおくと, sk = m j=1 jqj (j − )k (1) 37 / 46
38. さらに, ヴァンデルモンド行列 Vm = 1 1 1 1 (1 − ) (2 − ) (3 − ) (m − ) (1 − )2 (2 − )2 (3 − )2 (m − )2 . . . ... (1 − )m−1 (2 − )m−1 (3 − )m−1 (m − )m−1 を考えると, (s0...sm−1 ) = (1q1...mqm )VT m となる. 38 / 46
39. k = p∗ k v は定数なので, 固有ベクトルを求める上では無視してよい. → 先ほどの式の k の部分を無視した, (s0...sm−1 ) = (q1...qm )VT m で考える. すなわち, 求めたい固有ベクトルの近似 ˆ q1, ..., ˆ qm は (ˆ q1...ˆ qm ) = (ˆ s0...ˆ sm−1 ) ˆ V−T m で計算できる（−T は逆行列の転置 = 転置の逆行列）. 残る ˆ V−T m はどうするか？ 39 / 46
40. 櫻井‧杉浦法では, ˆ k から ˆ H< m と ˆ Hm を考えて, 一般固有値問題 (ˆ H< m − ˆ Hm )x = 0 を解き, 解 1, ..., m を求めると, これらが 1 − , ..., m − の近似になっているので, ˆ j ← + j (1 ≤ j ≤ m) を出力するのであった. したがって, ヴァンデルモンド行列 Vm の近似 ˆ Vm は ˆ Vm = 1 1 1 1 1 2 3 m 2 1 2 2 2 3 2 m . . . ... m−1 1 m−1 2 m−1 3 m−1 m である. 40 / 46
41. ラグランジュ補間 n 個の点 {xi } での値 yi = p(xi ) が与えられているとき, n − 1 次多項式 p(x) を求めよ. 例: {(xi, yi )} = {(1, 1), (2, 4), (3, 9)} → p(x) = x2 直接 xk (0 ≤ k ≤ n − 1) の係数をそれぞれ求めるのではなく, 多項式 pi(x) の重ね合わせ p(x) = n i=1 yi i (x) で計算する. i (x) は, i (xj ) が i = j の場合だけ 1, あとは 0 をとってほしい. 41 / 46
42. ラグランジュ多項式 i (x) = j≠i (x − xj ) j≠i (xi − xj ) 例えば {(xi, yi )} = {(1, 1), (2, 4), (3, 9)} の場合, 1 (x) = x − 2 1 − 2 · x − 3 1 − 3 , 2 (x) = x − 1 2 − 1 · x − 3 2 − 3 , 3 (x) = x − 1 3 − 1 · x − 2 3 − 2 で, p(x) = 1 · 1 (x) + 42 (x) + 93 (x) = x2 42 / 46
43. これを使えば, 1 x1 x2 1 xn−1 1 1 x2 x2 2 xn−1 2 1 x3 x2 3 xn−1 3 . . . ... 1 xn x2 n xn−1 n a1 a2 a3 . . . an = b1 b2 b3 . . . bn という式で {xi } と {bi } から {ai } を求めたい場合, 直接左辺の行列の 逆行列を計算せずとも, ラグランジュ補間で多項式を計算し, その係 数から {ai } を導くことができる. これにより, bi = 1, bj = 0 (j ≠ i) を 1 ≤ i ≤ n で考えて, それぞれの ai (1 ≤ i ≤ n) を計算すれば, 左辺の行列の逆行列が求まる（直接逆行 列を計算するとべき乗の計算量が重いので, こちらの方が有利）. 43 / 46
44. 今回の場合でも, ˆ VT m = 1 1 2 1 m−1 1 1 2 2 2 m−1 2 1 3 2 3 m−1 3 . . . ... 1 m 2 m m−1 m の逆行列 ˆ V−T m を同様に求めることができるので, (ˆ q1...ˆ qm ) = (ˆ s0...ˆ sm−1 ) ˆ V−T m が完了. 44 / 46
45. 余談: 固有ベクトル同士の関係 元の問題 (A − B)x = 0 に対し, ワイエルシュトラスの標準形 P(zB − A)Q = zId − Jd O O zJn−d − In−d を考えると, Q の列ベクトル qk が行列束 A − zB の固有ベクトルとな ることを先ほど述べた. 一方, 先週の内容より H< m − zHm = VmDm ΛVT m − zVmDmVT m = VmDm (Λ − zI)VT m ただし, Dm = diag(1, ..., m ), Λ = diag(1 − , ..., m − ) 45 / 46
46. これをワイエルシュトラスの標準形の形に式変形すると（普通の対角 化と変わらない）, (VmDm )−1(zHm − H< m )V−T m = zI − Λ であり, 前と同様に考えると V−T m の列ベクトルが行列束 H< m − zHm の 固有ベクトルであることが言える. つまり, 先ほどの (q0...qm−1 ) = (s1...sm )V−T m という式は, 2 つの一般固 有値問題の固有ベクトルを (s1...sm ) という行列で変換しているとい う意味だった. 46 / 46