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

Slide 1

Slide 1 text

櫻井‧杉浦法演習 III 発表杉本海人 2020/06/24 1 / 46

Slide 2

Slide 2 text

読んだ論文 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

Slide 3

Slide 3 text

概要与えられた n × n の複素行列 A, B に対して Ax = Bx ⇔ (A − B)x = 0 を満たすような（一般）固有値を求める問題（一般固有値問題） A, B はサイズの大きな疎行列で, 全部を求めたいのではなく, 一部だけ求まればいいというような状況 3 / 46

Slide 4

Slide 4 text

概要 B が逆行列を持てば B−1Ax = x となり, 行列 B−1A の固有値問題に帰着されるとはいえ, 数値計算上はこのようにしてもあまり意味がない • B−1A は依然としてサイズが大きい • A, B が対称‧疎など, 数値計算上良い性質を持っていたとしても, B−1A でそれが保持されているとは限らない →櫻井‧杉浦法では, 元の問題を, サイズの小さな行列の一般固有値問題に帰着させる（その解から元の一般固有値問題の解の一部が求まる） 4 / 46

Slide 5

Slide 5 text

流れ一般固有値問題 (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

Slide 6

Slide 6 text

用語の定義 • (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

Slide 7

Slide 7 text

例 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

Slide 8

Slide 8 text

ジョルダン標準形任意の 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

Slide 9

Slide 9 text

ワイエルシュトラスの標準形行列束 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

Slide 10

Slide 10 text

これを用いると 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

Slide 11

Slide 11 text

① 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

Slide 12

Slide 12 text

第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

Slide 13

Slide 13 text

第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

Slide 14

Slide 14 text

第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

Slide 15

Slide 15 text

② 閉路 Γ に対して, 周回積分 k = 1 2i ∫ Γ (z − )kf (z)dz を要素とするような行列 Hm, H< m を考えると, 一般固有値問題 (H< m − Hm )x = 0 の解が元の問題の閉路 Γ 内における解を与えるものになっていることを示す複素平面上で単一閉曲線（自身と交わらない閉曲線, ジョルダン曲線とも）Γ を考える. Γ 内にある固有値について考えていく. ※ 数値計算上は Γ はを中心とする円周で考えるが, 今はいったん任意の閉曲線で考えておく 15 / 46

Slide 16

Slide 16 text

複素解析 • 複素関数 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

Slide 17

Slide 17 text

複素解析 • 孤立特異点: 複素関数 f (z) で正則にならない点 • ローラン展開: f (z) の孤立特異点を a としたとき, f (z) = ∞ n=−∞ cn (z − a)n という形で展開できる（テイラー展開を負の次数方向にも広げたもの）. • 特に, 負の次数の部分, すなわち ∞ n=1 c−n (z − a)−n を主要部と呼ぶ. 17 / 46

Slide 18

Slide 18 text

複素解析 • 主要部について, 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

Slide 19

Slide 19 text

複素解析 • また, 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

Slide 20

Slide 20 text

複素解析 • 留数定理 • 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

Slide 21

Slide 21 text

本題に戻ると... 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

Slide 22

Slide 22 text

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

Slide 23

Slide 23 text

この 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

Slide 24

Slide 24 text

これらの 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

Slide 25

Slide 25 text

すると, 行列束 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

Slide 26

Slide 26 text

したがって, 行列束 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

Slide 27

Slide 27 text

残る問題 H< m , Hm の各要素である, 周回積分 k = 1 2i ∫ Γ (z − )kf (z)dz はどうやって求める？ → Γ はを中心とする円周で考える. すなわち, を円の半径として, j = + e2i N j (0 ≤ j ≤ N − 1) での値を使って近似する（詳細は次ページ） 27 / 46

Slide 28

Slide 28 text

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

Slide 29

Slide 29 text

一般固有値問題 (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

Slide 30

Slide 30 text

一般固有値問題 (A − B)x = 0 に対する櫻井‧杉浦法のアルゴリズム ⑤ 先ほど求めた ˆ k から ˆ H< m と ˆ Hm を考えて, 一般固有値問題 (ˆ H< m − ˆ Hm )x = 0 を解き, 解 1, ..., m を求める. ⑥ ˆ j ← + j (1 ≤ j ≤ m) 出力 ˆ 1, ..., ˆ m 30 / 46

Slide 31

Slide 31 text

並列化 A, B がサイズの大きな行列の場合, 計算に最も時間がかかるのは, f (j ) = u∗(jB − A)−1v を各 j に対して求める部分（前述アルゴリズムの②と③の部分） →この計算を並列計算することで, より高速に解を求められるようになる. 31 / 46

Slide 32

Slide 32 text

アルゴリズムの入力値の決定櫻井‧杉浦法では, 積分経路 Γ の中心 , 半径と, その経路の内部に含まれる固有値の個数 m を, 入力として適切に決めねばならない. 元論文ではこの入力値の取り方に関しては説明がされていないが, 経路内の固有値の個数を決定する手法を議論した論文を引用で紹介している. 32 / 46

Slide 33

Slide 33 text

固有ベクトルの求め方櫻井‧杉浦法のステップ 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

Slide 34

Slide 34 text

この 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

Slide 35

Slide 35 text

z = k を代入すると, (kB − A)qk = 0 となる. これは (k, qk ) が Ax = Bx を満たす (, x) のペアであることを意味する. Q が正則であることから qk ≠ 0 である点にも注意して, qk が k に対応する固有ベクトルである. 何らかの方法でこの qk を求めたい. 35 / 46

Slide 36

Slide 36 text

櫻井‧杉浦法のステップ 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

Slide 37

Slide 37 text

この 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

Slide 38

Slide 38 text

さらに, ヴァンデルモンド行列 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

Slide 39

Slide 39 text

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

Slide 40

Slide 40 text

櫻井‧杉浦法では, ˆ 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

Slide 41

Slide 41 text

ラグランジュ補間 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

Slide 42

Slide 42 text

ラグランジュ多項式 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

Slide 43

Slide 43 text

これを使えば, 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

Slide 44

Slide 44 text

今回の場合でも, ˆ 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

Slide 45

Slide 45 text

余談: 固有ベクトル同士の関係元の問題 (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

Slide 46

Slide 46 text

これをワイエルシュトラスの標準形の形に式変形すると（普通の対角化と変わらない）, (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