95
第 10 章 固有値の計算
本章では行列の固有値と固有ベクトルの計算法を学ぶ。固有値問題は対象と
する行列が対称か非対称かで解法が大きく異なるが,本章では実対称行列(ま
たはエルミート行列)向けの解法のみを取り上げる。最初に固有値の応用と
数学的性質についての簡単な復習を行った後,すべての固有値・固有ベクト
ルを同時に求める方法であるヤコビ法について説明する。次に,絶対値最大
の固有値とそれに対応する固有ベクトルを求めるためのべき乗法を説明する。
10.1
固有値とその応用
10.1.1
定義
正方行列 A に対し ,
Ax =
x
を満たす 0 でないベクトルが存在するとき, を
x を固有ベクト ルと呼ぶ。n
ど
(10.1)
A の固有値と呼ぶ。また,
n 行列の場合,固有値は重複度を含めてちょう
n 個存在する。よく知られているように,固有値は方程式 det(A ; I ) = 0
の解として特徴付けられる。
10.1.2
応用
固有値・固有ベクトルの計算は,自然科学や工学において次のような様々
な応用がある。
時間に依存しない Schr
odinger 方程式の解の計算。量子力学において,
波動関数とエネルギー順位は固有値問題 H =
E の解として与えら
れる。これは偏微分方程式の固有値問題であるが,これをポアソン方程
式の場合と同様に差分法で離散化すると,行列の固有値問題となる。
構造力学によるビルや橋梁の固有振動数の計算。地震などの揺れが建造
物の固有振動数と近い場合,建造物の揺れは大きく増幅され( 共鳴),
倒壊しやすくなる。そこで,設計に当たって固有振動数を計算しておけ
ば,どのような揺れに対して弱いかを予測できる。また,共鳴を回避す
るように設計し直すこともできる。
96
第 10 章 固有値の計算
統計計算。主成分分析,正準相関分析などの統計解析手法では,固有値
の計算が必要となる。
情報検索。固有値と密接な関係にある特異値を用いて,文書間の類似性
を自動的に判断し ,大量の文書の中から自分の探したい情報を持つ文
書を自動的に抽出することができる。
なお,応用では,一部の固有値・固有ベクトルのみを求めればよい場合も多
い。たとえば量子力学では,基底状態を求めたい場合,もっとも小さい固有
値とそれに対応する固有ベクトルのみを求めればよい(ただし励起状態も必
要な場合は,より大きな固有値・固有ベクトルも求める必要がある)。また,
主成分分析では,大きい方から数個の固有値を求めればよい場合が多い。
10.1.3
数学的性質
以下では,特に A が実対称行列またはエルミート行列の場合を扱う。上で
挙げた Schroedinger 方程式の解の計算では行列はエルミート行列である。ま
た,固有振動数の計算,統計計算などでは,多くの場合,行列は実対称行列
となる。
A が実対称行列またはエルミート 行列ならば,その
すべての固有値は実数である。そこで,n 個の固有値に小さい順に番号を付け,
固有ベクト ルの直交性
1
とする。また,
i
2
に属する固有ベクトルを
る固有ベクトルは直交する。すなわち,
i
(10.2)
n
x
と書く。異なる固有値に属す
i
6=
j
ならば
xx
i
j
= 0 が成り立つ
x
はエルミート共役を表す)
。さらに,固有ベクトルの組 f i gn
i=1 は
n
の任意のベクトルはこれによ
正規直交系をなすように取ることができ,
(ただし
R
り展開できる。
対角化
固有ベクトルを並べてできる行列を
X = (x1 x2 : : : x
n
)
(10.3)
とし ,固有値を対角成分に並べてできる対角行列を
0
BB 01
=B
B@ ...
0
0
2
..
.
..
.
0
1
C
0 C
.. C
C
. A
0
(10.4)
n
とすると,n 組の固有値・固有ベクトルに対する式 (10.1) はまとめて
AX = X (10.5)
10.2. ヤコビ法
97
と書ける。X の列は正規直交系をなすから,X は直交行列( A がエルミート
行列の場合はユニタリ行列)で,X
ら X をかけると
X = I が成り立つ。そこで,上式に左か
X AX = (あるいは X ;1AX = )
という式が得られる。これは
(10.6)
A を相似変換により対角行列 に変換したこ
とになるので,対角化と呼ぶ。実対称行列またはエルミート行列に対しては,
対角化と,すべての固有値・固有ベクトルの組を求めることとは等価である。
行列 A に対し,正則な行列
P による相似変換を行って得られる行列 B = P ;1 AP の固有値・固有ベクト
相似変換を行った行列の固有値・固有ベクト ル
ルを考える。式 (10.1) を変形すると
x
(P ;1 AP )(P ;1 ) =
すなわち
P ;1 x
B (P ;1 x) = P ;1 x
(10.7)
(10.8)
x
となるが,これは B の固有値が ,固有ベクトルが P ;1 であることを示し
ている。すなわち,相似変換により固有値
はP
;1
x に変換される。
は不変であり,固有ベクトル
x
10.2 ヤコビ法
10.2.1
原理
基本方針
ヤコビ法は,実対称行列( またはエルミート行列)のすべての固
有値・固有ベクトルを求めるための解法である。この解法では,固有値が相
似変換によって変化しないことを利用し ,ある直交行列の列 G1 G2 G3 : : :
を用いた相似変換
A1
A2
A3
=
=
=
..
.
G1 AG1
G2 A1 G2 G3 A2 G3 (10.9)
により,行列 A が徐々に対角行列に近づくようにする。Ak が十分に対角行列
に近くなると(すなわち,非対角成分の大きさが無視できる程度になると ),
その対角成分を固有値と見なすことができる。
98
G
k
第 10 章 固有値の計算
の定め方
いま,Ak が与えられたとき,Ak+1 を求めるには次のように
する。まず,Ak の絶対値最大の要素を apq とする。Ak ! Ak+1 の変形では,
この要素を相似変換により 0 にすることを考える。そのため,Gk として次の
形の行列を用いる。
01
BB
BB
BB
BB
BB
B
G =B
BB
BB
BB
BB
BB
B@
..
.
1
cos sin 1
..
k
.
1
; sin cos 1
..
.
1
1
CC
CC
CC
CC
CC
CC
CC
CC
CC
CC
CC
CA
(10.10)
これは,第 p 行,第 q 行の 2 つの変数が作る平面内での回転行列であり,直
Givens 変換と呼ぶ。G の形から
明らかなように,Givens 変換 A +1 = G A G では,変化を被る要素は第 p
行または第 q 行に属する要素と第 p 列または第 q 列に属する要素のみである。
A の要素を a ,A +1 の要素を a0 と書くと,a0 は次のように表される。
交行列である。この Gk による相似変換を
k
k
ij
k
a0
a0
a0
a0
a0
a0
a0
a0
ij
=
pj
=
qj
=
ip
=
iq
=
pp
=
qq
=
pq
=
k
k
k
k
ij
a
a
a
a
a
a
a
ij
pj
pj
ip
ip
pp
ij
( i 6= p q かつ j 6= p q のとき)
cos ; aqj sin ( j 6= p q )
sin + aqj cos ( j 6= p q )
cos ; aiq sin ( i 6= p q )
sin + aiq cos ( i 6= p q )
cos2 + aqq sin2 ; 2apq sin cos + a cos + 2a sin cos 1
(a ; a ) sin 2 + a cos 2
2
pp
2
sin
pp
2
qq
qq
pq
pq
(10.11)
(10.12)
(10.13)
(10.14)
(10.15)
(10.16)
(10.17)
(10.18)
ここで,0 にしたい要素は a0pq であるから,
a0
pq
すなわち
=
1
(app ; aqq ) sin 2 + apq cos 2 = 0,
2
tan 2 =
となるように回転角 を選ぶ。
a
;2apq
; aqq
pp
(10.19)
(10.20)
10.2. ヤコビ法
99
この相似変形によって,本当に Ak+1 は
非対角成分が減少することの証明
A
k
に比べて対角行列に近づいたのかど うかを調べよう。そのため,対角行列
への近さの尺度として,Ak の非対角要素の 2 乗和
XX
n
r
k
n
=1 j =1
X
a2 ;
ij
を考え,rk と
(10.21)
ii
=1
i
a2
i
r +1 の大きさを比較する。いま,行列 A が対称である( 証
明せよ)ことを用いると,r は行列のトレース1 を用いて次のように書き直
k
k
k
せる。
r
k
0
X @X
n
=
i
=1
i
=
j
X
n
=
=1
1
X
a a A ; a2
n
ij
=1
a2
ii
X=12
i
i
=1
ii
=1
i
X
(A2k )ii ;
Tr(A2k ) ;
ji
a
(10.22)
ii
この式を rk+1に対して適用し ,公式 Tr(AB ) = Tr(BA) と Gk Gk =
ると,
r +1
k
=
=
=
=
X
Tr(A2k+1 ) ;
=1
a0 2
ii
X
i
Tr(Gk Ak Gk Gk Ak Gk ) ;
=1
X
i
Tr(Ak Gk Gk Ak Gk Gk ) ;
Tr(A2k ) ;
X
=1
I を用い
a0 2
ii
a0 2
ii
i
a0 2
ii
X 2=1 X 0 2
= rk +
aii ; aii
i=1
i=1
2
2
2
2
= rk + (app + aqq ) ; (a0pp + a0qq )
i
となる。
この式をさらに簡単化するため,Ak+1 =
GAG
k
k
k
p 列と第 q 列を抜き出してできる次の式を考える2 。
a0
a0
pp
qp
a0
a0
pq
qq
! =
cos
sin
; sin
cos
P
!
a
a
pp
qp
a
a
!
pq
qq
(10.23)
の第 p 行と第 q 行,第
cos
; sin
sin
cos
!
(10.24)
n n の A に対し ,Tr(A) ni=1 aii を A のト レ ースと呼ぶ 。nP n の行列 A,
n
BPに対し
) = Tr(BA) という公式が成り立つ。実際,Tr(AB ) =
(AB )ii =
i=1
;P,nTr(AB
;
P
P
P
n
n
n
n
a
b
=
b
a
=
(BA)kk = Tr(BA) である。
ik
ki
ki
ik
i=1
k=1
k=1
i=1
k=1
1
2 以下の計算は,式 (10.26) を導き出すのが目的である。行列のトレースに不慣れな読者は式
(10.16),(10.17),(10.18) を用いて直接計算しても式 (10.26) が得られるが,計算はずっと面
倒になる。
100
第 10 章 固有値の計算
k+1 = G
Ak G
k と書くと,G
k が直交行列であることを用いて,
この式を A
k
a0 2 + a0 2 + a0 2 + a0 2
pp
pq
qp
qq
Tr(A2k+1 )
=
Ak G
kG
Ak G
k)
Tr(G
k
k
=
kG
Ak G
kG
)
Tr(Ak G
k
k
=
Tr(Ak Ak )
=
a2
=
pp
+ a2pq + a2qp + a2qq
(10.25)
さらに,aqp = apq ,a0qp = a0pq = 0 を用いると次の式が得られる。
a0 2 + a0 2 = a2
pp
qq
pp
+ a2qq + 2a2pq
(10.26)
これを式 (10.23) に代入すると,最終的に
r +1 = r ; 2a2
k
k
(10.27)
pq
となる。したがって,Ak+1 の非対角成分の 2 乗和 rk+1 は
の 2 乗和 rk に比べて着実に減少し ,Ak+1 は
ていることがわかる。
A
k
A
k
の非対角成分
に比べて対角行列に近づい
ダウンロード

講義ノート