2014年6月11日（水） - 青柳睦

2014年6月11日（水）
数値解析講義
前期水曜３時限
情報理学コース
青柳睦
aoyagi@cc.九州大学のドメイン
研究室
箱崎キャンパス情報基盤研究開発センター５階５０２，５００
伊都キャンパスウエスト2号館 1006室右ドア側
1
2．連立一次方程式の解法
講義内容
:
2.5
2.6
2.7
2.8
2.9
｛特異値分解法｝
共役勾配法 →演習
反復法 2.7.1 ヤコビ法 →演習
ガウス・ザイデル法
｛ＳＯＲ法（概要のみ解説）｝
3．固有値と固有ベクトルの数値計算
今日の講義
3.1 固有値問題の序論
3.2 ヤコビの方法（対角行列への帰着） →演習
3.3 ギブンスの方法（３重対角行列への帰着）→演習？
3.4 DKA法（３重対角行列の固有値の計算 )
3.5 ハウスホルダー法（３重対角行列への帰着）
4. その他の数値計算法
4.1 ・・・資料作成中・・
4.2
2
3．固有値と固有ベクトルの数値計算
3．1
固有値問題の序
固有値 λ
Ax = λ x
A : n × n 行列
x : 0でない n 次元ベクトル
det ( A − λ I ) ≠ 0 なら A − λ I の逆行列が有
0
( A − λI ) x =
x は0のみとなる（前提に矛盾）
det ( A − λ I ) =
0 : 固有方程式
( n 次方程式 )
3
det ( A − λ I ) =
0 : 固有方程式( n 次方程式)
固有値は n 次の代数方程式の解
一般に ( n ≥ 5 ) 有限回の代数処理では(厳密には)求まらない(*)
反復法による多数の繰り返し演算が必要
(実際には、要求される精度に達した時点で打ち切る)
固有値問題は、連立一次方程式同様、工学・理学・社
会学の分野で非常に重要（固有振動解析、量子力学、他・・）
* ガロア、アベル理論
4
相似変換
−1
B = P AP
（本題への準備
P : n×n
）
正則(逆行列のある)行列P
による相似変換で、
固有多項式は変わらない
固有値も変わらない
【証明】
(
)
det ( P AP − λ P P )
det ( P ( A − λ I ) P )
det ( ( A − λ I ) PP )
det ( B −
λ I ) det P −1 AP − λ I
=
=
=
=
−1
−1
−1
= det ( A − λ I )
−1
（公式）行列積の行列式
det ( XY ) = det (YX )
det ( XYZ ) = det ( X (YZ ) )
= det ( (YZ ) X )
= det (YZX ) 5
A1 ≡ A
対角形に収束す
るまで繰り返す

A2 = P −1 A1 P 
 相似変換の繰り返し
−1
A3 = P A2 P 


一旦、３重対角
形に変換する
3.2 ヤコビ法
*



*
0




*







*
0


* 

対角行列
* *


0 
* * *


* * 


   

 0
 * *


*
*


3.5 ハウスホルダー法
3重対角行列
3.4 DKA法
対角要素→固有値
3.3 ギブンス法
バイセクション法
固有値は３重対角行列に適した
別の反復法により求める
6
３．２ヤコビの方法
1




0



 p
cos φ  sin φ


P=

 


 q
− sin φ  cos φ


 
 0


1


p
PP = I
P −1 = t P 直交行列
t
相似変換の行列 P
q
Pij = δ ij
ただし、
Ppp = cos φ
Pqq = cos φ
Ppq = sin φ
Pqp = − sin φ
（毎回の）変換後の行列 B = P −1 AP の特定の成分を0にする
7
相似変換 B = P −1 AP において
1) A が対称 A = A ならば
t
−1
t
=
B t ( P=
AP)
t
−1
t
PA =
( P −1 ) P=
AP B
よって B も対称
2)
A の p, q 行目と
p, q 列目のみが変更される
以上のように回転行列による相似変換を（三重対角or
対角形になるまで）繰り返す。
GIVENS法：三重対角形に変換する（次節 3.3）。
JACOBI法：対角形に変換する（本節 3.2）。
8
まず左からP-1 をかけると、
P A=
−1
1
  a11

 
0




  a p1
cos φ  − sin φ






 

  aq1
sin φ  cos φ


0


 

1   an1

 a1 p  a1q  a1n 

 

 
 a pp  a pq  a pn 

  
 
 aqp  aqq  aqn 


   
 anp  anq  ann 
 a11

 
 a′p1

 
 aq′1

 
 an1

=
a′pk a pk cos φ − aqk sin φ
 ′
=
 aqk a pk sin φ + aqk cos φ
p
=
 a1 p
 
 a′pp

′
 aqp

 anp
p
q
 a1q  a1n 


 
 a′pq  a′pn 

 
 
′  aqn
′ 
 aqq

   
 anq  ann 
q
p
q
p
p
q
q
k = 1,2,..., n
列 index
9
さらに、右からP をかけると、
B = P −1 AP
=
=
 a11

 
 a′p1

 
 aq′1

 
 an1

 a1 p
 
 a′pp

′
 aqp
 a11

 
 a′p1

 
 aq′1

 
 an1

 a1′′p  a1′′q  a1n 

 

 
 a′′pp  a′′pq  a′pn 

  
 
′′  aqq
′′  aq′n 
 aqp


   
′′  anq
′′  ann 
 anp

 anp
p
p
 a1q  a1n   1




 

0



 a′pq  a′pn  
cos φ  sin φ


 
 






′  aqn
′ 
 aqq
− sin φ  cos φ


   
 
0
 anq  ann  
1 
p
q
q
q
p
q
′′ akp cos φ − akq sin φ
=
 akp
 ′′
=
akq akp sin φ + akq cos φ
行 index
k = 1,2,..., n
k ≠ p, q
=
a′′pp a pp cos 2 φ + aqq sin 2 φ − (a pq + aqp ) cos φ sin φ

′′ a pp sin 2 φ + aqq cos 2 φ + (a pq + aqp ) cos φ sin φ
aqq
=

2
2
′′
=−
+
−
a
(
a
a
)
cos
φ
sin
φ
a
cos
φ
a
sin
φ
pp
qq
pq
qp
 pq
2
2
 aqp
′′
=
−
+
−
(
a
a
)
cos
φ
sin
φ
a
cos
φ
a
sin
φ
pp
qq
pq
qp

10
aij = a ji を用いると結局,
i, j, k = 1,2, n
i, j ≠ p, q
 bij aij 後で使う・・・・（＊）
b =
ただし、 k ≠ p, q
=
−
φ
φ
b
a
a
cos
sin
pk
kp
pk
qk

 bqk =
bkq =
a pk sin φ + aqk cos φ
k ≠ p, q

a pp + aqq a pp − aqq
b =
+
cos 2φ − a pq sin 2φ
pp

2
2

a pp + aqq a pp − aqq
b =
−
cos 2φ + a pq sin 2φ
qq

2
2

a pp − aqq

=
b
sin 2φ + a pq cos 2φ
 pq
2
bpq = 0 とするためには cot 2φ =
aqq − a pp
とすればよい．
2a pq
11
実際の計算プログラムでは、
aqq − a pp
= α とおくと
cot 2φ =
2a pq
∴ tan φ =
sgn(α )
α + α 2 +1
tan 2 φ + 2α tan φ − 1 = 0
より、
tan φ = −α ± α 2 + 1
小さい方の根は、
α > 0 ⇒ t = −α + α 2 + 1 =
≡t
−1
−α − α 2 +1
α < 0 ⇒ t = −α − α 2 + 1 =
=
1
α + α 2 +1
−1
−α + α 2 +1
このtよりBの各要素を計算する．
 bij
b
 pk
 bqk

b
 pp

b
 qq

bpq

aij
bkp =
a pk cos φ − aqk sin φ
=
bkq =
a pk sin φ + aqk cos φ
=
a pp + aqq a pp − aqq
cos 2φ − a pq sin 2φ
= +
2
2
a pp + aqq a pp − aqq
cos 2φ + a pq sin 2φ
= −
2
2
a pp − aqq
sin 2φ + a pq cos 2φ
2
（次ページ参照）
i, j ≠ p, q
k ≠ p, q
k ≠ p, q
p, q の選び方
絶対値が最大の非対角
要素が0となるように選ぶ
一度0になった要素が後の変
換で非ゼロに復活してくること
もありうる。（次に「収束」の
保証について考察する）
12
計算式(参考）
b pp =
=
a pp + aqq
2
a pp + aqq
+
+
a pp − aqq
2
a pp − aqq
bqq = a pp + aqq − b pp
= aqq + a pq tan φ = aqq + ta pq
cos 2φ − a pq sin 2φ
また、
tan φ =
(1 − 2 sin 2 φ ) − a pq sin 2φ
2
2
= a pp − (a pp − aqq ) sin 2 φ − a pq sin 2φ
cot 2φ =
= a pp + 2a pq cot 2φ sin 2 φ − a pq sin 2φ
aqq − a pp
2a pq
cos 2φ sin 2 φ
= a pp + 2a pq
− a pq sin 2φ
sin 2φ
cos 2φ sin φ
= a pp + a pq
− 2a pq sin φ cos φ
cos φ
(cos 2 φ − sin 2 φ ) sin φ − 2 sin φ cos 2 φ
= a pp + a pq
cos φ
− sin φ (− cos 2 φ + sin 2 φ + 2 cos 2 φ )
= a pp + a pq
cos φ
= a pp − a pq tan φ = a pp − ta pq
より、
sgn(α )
α + α +1
2
≡t
1
より、
cos 2 φ
1
1 + tan 2 φ =
cos φ =
1+ t 2
sin φ = tan φ cos φ =
t
1+ t 2
を用いて、
b pk = bkp = a pk cos φ − aqk sin φ
bqk = bkq = a pk sin φ + aqk cos φ
k = 1,2,  n　( k ≠ p, q )
∴ b pp = a pp − ta pq
計算する。
13
ヤコビの対角化法に関する「収束について」
まず準備、
相似変換 B = P −1 AP において次が成り立つ
(a)
2
2
a
=
b
∑ ij ∑ ij
i, j
i, j
(全要素の２乗の和は変わらない)
−1 t
は直交行列
(
P
= P) だから
P
t
BB
t
−1
(P
=
AP)( P −1 AP)
t
t
PA
=
(P −1 ) P −1 AP P −1 t AAP
対角和は相似変換によって不変だから（次ページ参照）
tr At A = tr B t B
成分ごとに書くと
2
2
a
=
b
∑ ij ∑ ij
i, j
i, j
14
補足
対角和(トレース)
n
tr X = ∑ xii
i =1
XY の対角和 = YX の対角和
n
n  n
 n  n

=
=
=
tr XY=
( XY )ii ∑
 ∑ xij y ji  ∑
∑
 ∑ y ji xij  tr YX
=i 1
=i 1 =
i1

j1
 =j 1 =
相似変換によって変わらない
(
)
−1
−1
=
tr Y −1 XY tr=
Y −1 ( XY ) tr=
XY
Y
tr
=
X
YY
tr X
( )
15
(b)
2
2
bpp
+ bqq2 = a 2pp + 2a 2pq + aqq
の証明．
2
2
2
まず bpp
を示す
+ 2bpq
+ bqq2 = a 2pp + 2a 2pq + aqq
a pp + aqq a pp − aqq

cos 2φ − a pq sin 2φ
bpp = +
2
2

a pp + aqq a pp − aqq

cos 2φ + a pq sin 2φ
 bqq = −
2
2

a pp − aqq

sin 2φ + a pq cos 2φ
bpq
2

a pp + aqq
a pp − aqq
において F =
=
, G =
, H a pq とおくと
2
2
bpp =
F + G cos 2φ − H sin 2φ

F − G cos 2φ + H sin 2φ
 bqq =
b
=
 pq G sin 2φ + H cos 2φ
16
b + 2b + b = { F + ( G cos 2φ − H sin 2φ )} + ( F − ( G cos 2φ − H sin 2φ ) )
2
pp
2
pq
2
2
qq
+2 ( G sin 2φ + H cos 2φ )
2
2 F 2 + 2 ( G cos 2φ − H sin 2φ ) + 2 ( G sin 2φ + H cos 2φ )
=
2
{
2 F 2 + 2 G 2 cos 2 2φ + H 2 sin 2 2φ − 2GH cos 2φ sin 2φ
=
{
2(F
+2 G 2 sin 2 2φ + H 2 cos 2 2φ + 2GH cos 2φ sin 2φ
=
2
+ G2 + H 2
2
2
}
}
)
 a pp + aqq  2  a pp − aqq 2

2 
= 2 
 + a pq 
 +
2
2

 


2
 a 2pp + aqq

2 
= 2
+ a pq 
2


2
=a 2pp + 2a 2pq + aqq
17
2
2
2
bpp
+ 2bpq
+ bqq2 = a 2pp + 2a 2pq + aqq
において P は bpq = 0 となるように選んだから
2
2
bpp
+ bqq2 = a 2pp + 2a 2pq + aqq
（実は、次の様な式変形だけでも、
１回の相似変換で“非対角要素の二乗和”
が毎回減少することが言える・・。）
以上(a)、(b)からBの非対角要素の２乗の和は
（ｂよ)り


2
2
2
b = ∑ b −  ∑ bii + bpp + bqq 
∑
i≠ j
i, j
 i ≠ p ,q



2
2
2
2
2
= ∑ aij −  ∑ aii + a pp + 2a pq + aqq 
i, j
 i ≠ p ,q

= ∑ aij2 − 2a 2pq
2
ij
2
ij
（
＊）
（よ
a)り
となる
（
＊）
i≠ j
bij = aij
i, j ≠ p, q
18
2
=
b
∑ ij
i≠ j
2
2
a
−
2
a
∑ ij pq
i≠ j
において p, q は a 2pq = max aij2 となるように選んだから
i≠ j
全ての i ≠ j について
aij2 ≤ a 2pq
和をとると
∑a
i≠ j
2
ij
∑
2
の項数は
n
− n だから
i≠ j
≤ n(n − 1)a
2
2
2
2
a
≤
a
∑
ij
pq
n(n − 1) i ≠ j
2
pq
よって

2 
2
=
−
≤
−
2
1
b
a
a
a

 ∑ ij
∑
∑
i≠ j
i≠ j
 n(n − 1)  i ≠ j
2
ij
2
ij
2
pq
19

2 
2
≤
−
1
b
a

 ∑ ij
∑
i≠ j
 n(n − 1)  i ≠ j
2
ij
k − 1 回相似変換を行った後の行列 Ak の要素を aij( k ) と書くと
∑(
i≠ j
(k )
ij
a
)
2

2 
≤ 1 −

 n(n − 1) 
k
2
a
∑ ij
i≠ j
n ≥ 3 に対しては 0 < 1 − 2 n(n − 1) < 1 だから
k → ∞ のとき
∑(a )
i≠ j
(k )
ij
2
→0
よって・・(上のPによる) 相似変換の繰り返しで
非対角成分は全体として0に収束する
ただし理論的に「繰り返し数の厳密な上限」を
導出することは出来ない。
（あるepsまでの近似回転数は導出できる）
⇔
もし示せたら、n次方程式が
有限回演算で解けたことに
なってしまう！（矛盾）
20
「固有値と固有ベクトル」
k 回目の相似変換の行列を Pk と書く
−1
=
Ak +1 P=
k 1, 2,
k Ak Pk ,
T = lim P1 P2  Pk
k →∞
(実際の数値計算では、必要な精度が
得られたらkを有限回で打ち切る。)
を作ると、以上の結果から
Λ =T −1 AT
ここで Λ は対角行列
 λ1



λ2
0


Λ =


0




λ
n

λi : A の固有値
21
Λ =T −1 AT
から
AT = T Λ
T の第 i列からなる列ベクトルを ti とすると
T = (t1 , t2 , , tn )
Atk = λk tk
ti が A の固有値 λi に属する固有ベクトルとなる
【以下は補足】
相似変換によって行列Aに関する次の値は不変。
① 行列式det(A)
④ 固有値 λ
② トレースｔr(A)
③階数 rank(A)
（一次独立な固有ベクトルの本数）
②④から、固有値の和は tr(A)
①④から、固有値の積は det(A)
④から、 A-1 , A2, A3 ・・など Aの関数が簡単に求まる。
22
【実験１】 10～100次元のFrank 行列＊ Aについて、Jacobiの対角化法
を用い、Iterationの回数を調べた。
次元数
10
20
30
40
50
60
70
80
90
100
EPS=10-6
Iteration
122
465
991
1561
2275
2944
3662
4315
5062
5718
Frank行列の固有値 λk
λ−k1 = 2{1 − cos(
2k − 1
π )}, 1≤ k ≤ n
2n + 1
Frank行列の固有ベクトル vk
2k − 1 π
(2k − 1)(2 j − 1) π
t
vk = {cos(
2n + 1 2
 , cos(
), 2n + 1
2
), , cos(
(2k − 1)(2n − 1) π
)}
2n + 1
2
vkのj番目の要素
（参考）Frank行列の逆行列
 1 −1 0
− 1 2 − 1

−1
A =  0 −1 2




 0 0




−1
0
0 


− 1
2 
23
【脚注*】
Frank行列は厳密な固有値、固有ベクトルの解析解が求められる
【実験２】 Frank行列Aについて、Jacobiの対角化法を用い、
収束の様子を調べた。
n − 1 n − 2  1
 n
 n − 1 n − 1 n − 2  1


A =  n − 2 n − 2 n − 2  






1


 1
1
1 1

Frank行列
Color
40次元の場合の実行例。 CまたはJavaに似た言語でグラフィッ
ク処理が簡単な「Processing2.0.x」という処理系を用いて可視化し
た（例として、計算の途中の様子を図示した。）。
[ 参考 ] http://processing.org/
|aij|の値
100< |aij|
10 < |aij| < 100
1.0 < |aij| < 10.0
0.1 < |aij| < 1.0
10-2 < |aij| < 0.1
10-3 < |aij| < 10-2
10-4 < |aij| < 10-3
|aij| < 10-4
24

Download Report