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

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

例
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

View Slide

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

View Slide

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

View Slide

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

View Slide

①
f (z) = d
j=1
j
z − j
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

View Slide

第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

View Slide

第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

View Slide

第2項について
したがって, (In−d
− zJn−d
) のジョルダンブロックの最大サイズを K と
すると, 第 2 項 −u∗Q2
(In−d
− zJn−d
)−1P2v は次数 K − 1 の z の多項式に
なる（これを g(z) とおく）.
以上より,
①
f (z) = d
j=1
j
z − j
が完了.
14 / 46

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

したがって, 行列束 H<
m
− zHm
の相異なる固有値は 1
− , ..., m
− に
なっている.
一般固有値問題 (H<
m
− Hm
)x = 0 の解 1
− , ..., m
− を求めて, それ
ぞれにを足せば,
元の問題の Γ 内の解になっていると言える.
②
=
1
2i
∫
Γ
(z − )kf (z)dz を要素とするよ
うな行列 Hm, H<
m
を考えると, 一般固有値問題 (H<
m
− Hm
)x = 0 の解が
元の問題の閉路 Γ 内における解を与えるものになっていることを示す
が完了.
26 / 46

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

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

Kaito Sugimoto

More Decks by Kaito Sugimoto

Other Decks in Science

Featured

Transcript