一般化固有値問題
桂田 祐史

2003 年 1 月 28 日,2004 年 12 月 30 日
この文書は以下から入手可能。

http://nalab.mind.meiji.ac.jp/~mk/labo/text/generalized-eigenvalue-problem.pdf
標準固有値問題の計算法については

http://nalab.mind.meiji.ac.jp/~mk/labo/text/eigenvalues.pdf
http://nalab.mind.meiji.ac.jp/~mk/labo/text/eigenvalues-add.pdf

序1(行列束とか簡単なことを書きたい…暇が取れるか？)A, B ∈ M (n; C) とする。λ ∈ C, x ∈ Cn \ {0} が(1)Ax = λBx,x ̸= 0をみたすとき、λ を (A, B) の一般化固有値、x を λ に属する一般化固有ベクトルとよぶことにする (あまり普通の言い方ではない？)。
(A, B) の一般化固有値全体を σ(A, B) と書こう:

σ(A, B) := {λ ∈ C; ∃y ∈ Cn \ {0} s.t. Ay = λBy}.
正直に言って、あまりテキストがないので詳しいことが分からない。Golub and Van Loan
[4] くらいか。特に B が正則でもないような場合は私にとって五里霧中。

2標準固有値問題への帰着

2.1

問題点

一般固有値問題をそれと同値な標準固有値問題に変換する方法は、以下紹介するように色々ある。

1

標準固有値問題の解法として、Householder 法のように行列の成分が直接必要なアルゴリズムを利用する場合には、行列の疎性が保存されるかどうかが重要なポイントになる。
これに対して、冪乗法や Lanczos 法を利用する場合には、行列の疎性は必ずしも必要ではない。つまり掛け算x 7→ P xを効率的に行うのに必ずしも P の疎性は必要ないことに注意しよう。

P = P1 P2 · · · Prと分解できたとして、各 i につき、Pi または Pi −1 の一方が疎であれば良いわけである 1 。

2.2

B が正則な場合に使える、ある素朴な方法

e = B −1 A とおくと、B が正則な場合、Ae = λxAxとなる。


命題 2.1 (B が正則な場合) A, B は n 次正方行列、B は正則とするとき、

e := B −1 AA
とおくとき、次の (1), (2) が成り立つ。

(1) 任意の λ ∈ C, x ∈ Cn について
Ax = λBx⇐⇒e = λx.
Ax(2) 任意の λ ∈ C についてλ が (A, B) の一般化固有値⇐⇒⇐⇒

e の固有値λがA
e = 0.
det(λI − A)

ただし I は n 次単位行列である。


そこで (A, B) の一般化固有値 λ の (代数的) 重複度をdet(λI − B −1 A) = 0の重複度によって定義する 2 。これから (A, B) の一般化固有値も重複度を考えてちょうど n個だけ存在することが分かる。
つまり、P と a から b = P a を計算するとき、P −1 が疎であるならば、連立 1 次方程式 P −1 b = a を解くことによって b を求める、という手順が採用できる。
2
微少な変形をすることによっても定義できる？もっとも B が正則でないような場合まで含めて素朴に考えるのはきっと危険だ。
12

e は対称であるとは限らず、A,
注意 2.2 (計算向きの方法ではない) A, B が対称であっても Ae は疎であるとは限らず、Ae の固有値問題に帰着することによって固有値B が疎であっても Aを計算するのはあまり勧められる方法ではない (よほど小規模な問題でない限り問題外であろう 3 )。

2.3B が正定値である場合にその Cholesky 分解を用いる方法(この項の要点は式番号をふった式にある。)
B が正定値 (実対称で固有値がすべて正) である場合、Cholesky 分解
B = UT U(2)を代入して (U は正則な上三角行列であり、特に対角成分が正であるものを選ぶことができ、その場合は一意的に決定される)
Ax = U T U x.
( )−1をかけると両辺に U T( T )−1
Ax = U x.
U(

補題 2.3 P が正則であるとき、 P行列である。

)
T −1



= (P −1 ) . 特に正則な対称行列の逆行列は対称T




−1

証明
Q := Pとすると、P Q = I. これから Q P = (P Q) = I = I. ゆえに( T )−1T−1 T= Q = (P ) .
P( )−1のことをしばしば P −T と書く人が多い (個人的には悪趣味だと感じているが)。
なお、 P T
そこで
(3)
A′ := (U −1 )T AU −1 , y := U x
TTTT

とおくと明らかに(4)

Ax = λBx

(5)

x=0⇐⇒

⇐⇒

A′ y = λy,y = 0.

A が実対称であれば、A′ も実対称である。これは次の補題から分かる。



補題 2.4 P が n 次対称行列であるとき、任意の n 次正方行列 Q について P ′ := QT P Qは対称である。


証明(P ′ )T = (QT P Q)T = QT P T (QT )T = QT P Q = P ′ .
3

しかし、まったく役に立たないというわけではない。ずっと以前、某先生の卒研の学生が加藤敏夫先生の変分法の解説を勉強していて、微分方程式の固有値問題の数値計算に取り組んでいたとき、結局は Mathematica
で det(λI − B −1 A) = 0 を解くのが一番だった、ということがあった。「I 君の固有値問題」 (http://nalab.
mind.meiji.ac.jp/~mk/labo/text/eigen-I.pdf)を見よ。

3

補題 2.5 A, B は実対称行列で、B は正定値、B = U U T を Cholesky 分解 (すなわち Uは上三角行列) とするとき、A′ := (U −1 )T AU −1 とおくと、

det(λI − A′ ) = det(λI − B −1 A) (λ ∈ C).

証明

det(λI − A′ ) = det(λU U −1 − (U −1 )T AU −1 ) = det(λU −1 U − U −1 (U T )−1 A)= det(λI − (U T U )−1 A) = det(λI − B −1 A).

ここまでをまとめておこう。

4



定理 2.6 (対称な一般化固有値問題のコレスキー分解に基づく変換) A, B が n 次実対称行列で、B は正定値であるとき、B の Cholesky 分解 B = U T U (U は対角成分が正である上三角行列) を用いて、
A′ := (U −1 )T AUとおくと、次の (1)–(3) が成り立つ。

(1) A′ は実対称行列である (ゆえにその固有値はすべて実数であり、対応する固有ベクト
ルで Rn の正規直交基底をなすものが存在する)。
(2) A′ の固有値は (A, B) の一般化固有値に等しい。すなわち、σ(A′ ) := {λ ∈ C; ∃x ∈ Cn \ {0} s.t. A′ x = λx},σ(A, B) := {µ ∈ C; ∃y ∈ Cn \ {0} s.t. Ay = µBy}
とおくとき、σ(A) = σ(A, B). さらにdet(λI − A′ ) = det(λ − B −1 A)であるので重複度まで込めて等しい。

(3) A′ の固有値を重複度の分だけ並べて書いたものを
λ1 , λ2 , · · · , λn ,これらに対応する A′ の固有ベクトルからなる Rn の正規直交基底をv1 , v2 , · · · , vnとするとき、

ui := U −1 vi(i = 1, 2, · · · , n)

Aui = λi Bui(i = 1, 2, · · · , n),⟨ui , uj ⟩ = δij(i, j = 1, 2, · · · , n).

とおくと、

ただし ⟨·, ·⟩ は次式で定義される Rn の内積である。
証明⟨x, y⟩ := (Bx, y) = (x, By) = (U x, U y) (x, y ∈ Rn ).
(1), (2) と (3) の前半は済んでいる。(3) の後半については⟨ui , uj ⟩ = (U ui , U uj ) = (vi , vj ) = δij .

5

さて、実対称行列 A′ の標準固有値問題については、Rayleigh 商(A′ y, y)
y T A′ y= T(y, y)y yが重要な役目を果たすが、これについては次の補題が成り立つ。


補題 2.7 A, B が n 次実対称行列で、B は正定値であるとき、B の Cholesky 分解 B =U T U (U は対角成分が正である上三角行列) を用いて、A′ := (U −1 )T AUとおくと、任意の x ∈ Rn に対して

(Ax, x)(A′ y, y)=,(Bx, x)(y, y)(6)y := U x.



証明



最初に x ̸= 0 ⇔ y ̸= 0 であることを注意しておく。

(Ax, x)(AU −1 y, U −1 y)((U −1 )T AU −1 y, y)(A′ y, y)===.
(Bx, x)(U x, U x)(y, y)(y, y)


系 2.8 A, B が n 次実対称行列で、B は正定値であるとき、max σ(A, B) =(Ax, x),x∈R \{0} (Bx, x)
maxnmin σ(A, B) =max{|λ|; λ ∈ σ(A, B)} =(Ax, x),x∈R \{0} (Bx, x)
minn(Ax, x).
x∈R \{0} (Bx, x)
maxn

より詳しく min-max 原理も成り立つ。すなわち固有値全部を大きい方から順に (重複度の分だけ繰り返して) 並べたものを

λ1 ≥ λ2 ≥ · · · ≥ λnとするとき、1 ≤ k ≤ n なる k に対して、(Ax, x)= 大きい方から k 番目の固有値 = λk ,dim V =k x∈V \{0} (Bx, x)(Ax, x)min max= 小さい方から k 番目の固有値 = λn−k+1 .
dim V =k x∈V \{0} (Bx, x)
max

min

ここで dim V = k は V が Rn のすべての k 次元部分空間を走ることを意味する。
注意 2.9 (数値計算する上での注意) ただし、A, B が疎行列であっても A′ は疎行列であるとは限らない (ほとんどの場合は密行列になる)。
なお、冪乗法の系統のアルゴリズム (逆反復法, シフト法) を採用する場合には、事前にシフトをすませておく工夫が有用である。後述の Lanczos 法の節を見よ。

6

2.4

split Cholesky 分解を用いる方法

LAPACK で採用されている方法である。残念ながら適当な参考文献が入手できていない。
LAPACK のドライバー・ルーチン spbstf のソース・プログラム dpbstf.f 中の注釈くらいしか参考になるものがない。
http://www.netlib.org/cgi-bin/netlibfiles.pl?filename=/lapack/double/dpbstf.f
なお、LAPACK のソースをサーチするには

http://www.cs.colorado.edu/~jessup/lapack/applets/search.html
が便利である。
BをB = ST Sの形に分解するのだが、通常の Cholesky 分解とは異なり、S は上三角行列ではなく、()U 0, U = (n + kd )/2 次上三角行列, L = (n − kd )/2 次下三角行列,
S=M Lのような行列である。ここで kd は半バンド幅である。

A′ = P −T AP −1とおくと、前小節の補題 2.4 により A′ は対称であることが分かる。
A, B が帯行列である場合、A′ も (やや帯幅は増えるものの) 帯行列になる。

3二分法(準備中 — 村田先生の本から)

4

Lanczos 法

Lanczos 法の利点
一般固有値問題は、行列の疎性を保存したまま標準固有値問題に帰着するのが困難であるわ
けだが、Lanczos 法 (実対称行列の 3 重対角化のアルゴリズムである) ではその問題を回避できる。その理由は、Lanczos 法は Lanczos 原理 (Arnoldi 原理と呼ぶべき？) に基づいている
ため、問題の行列 A の成分そのものは必要でなく、任意のベクトル x との積 Ax さえ計算できればよいということにある。

7

標準固有値問題への帰着
Ax = λBxe ∈ R に対して、必要ならばシフトを行う。すなわち任意に選んだ λee
Ax − λBx= (λ − λ)Bx.
そこでee def.
A= A − λB,def.
eµ = λ−λとおくと

e = µBx.
Ax
B の Cholesky 分解を
B = UT Uとして、上の式に代入すると

e = µU T U x.
Ax
e−1 をかけると (ただし Ae−1 = (A)e −1 )
左から µ−1 U A
e−1 U T U x.
µ−1 U x = U A
そこでdef.
e−1 U T ,A′ = U A

ν = µ−1 ,def.

def.

y = Uxとおくと、

A′ y = νyという標準固有値問題に帰着された。
これに Lanczos 法を適用して、対称三重対角行列 T の固有値問題

T z = νzに変換する。

T を計算で求めるのが楽であること
A, B が帯行列であっても、A′ が疎性を持つことは期待できないが、任意に与えられた x ∈ Rnに対して、e−1 U T x
A′ x = U A
e−1 の掛け算 A
e−1 zを計算するのは簡単である。帯三角行列 U T , U の掛け算は容易であるし、Aは対称帯行列係数の連立 1 次方程式
e =zAw
e を Cholesky 分解しておけば非常に効率的に計算可能である。したを解くことであるから、Aがって Lanczos 法に現れる計算は楽である。
8

T についての固有値問題の解法
e = µBx の固有値の逆数でまず T の固有値は A′ の固有値に等しく、一般化固有値問題 Axあることに注意しよう。
T は三重対角であるから、例えば二分法 (bisection method) を使うこともできる。
e = µBx の絶対値最小のあるいは冪乗法を行うと、T の絶対値最大の固有値、すなわち Axe に最も近い固有値が求められるこ固有値の逆数が求まる。ゆえに元の問題 Ax = λBx の、λe = 0 とすると、絶対値最小の固有値が求まることになる。とになる。シフト λ(この辺り矢川・青山 [2] を参考にしたが、絶対値最小の固有値にこだわっているのは、扱っている問題の性質によるのだと思う。この点は整理する必要がある。)

精度保証付き数値計算

5

偏微分方程式の精度保証付き数値計算で、一般化固有値問題の固有値の絶対値の最大値 (これってスペクトル半径と言って良い？) の精度の良い上界が必要になる。そのため一般化固有値問題についても研究が行われている。当の中尾充宏先生のグループが色々試していて、その成果は渡部・山本・中尾 [3] に詳しく説明されている。
また、http://www.cc.kyushu-u.ac.jp/RD/watanabe/RESERCH/paper_list.html には院生の修士論文が色々ある。
桂田研の学生向けに読書ノートがhttp://nalab.mind.meiji.ac.jp/~mk/labo/members/reading-wyn1999.pdf
においてある。

A雑題

A.1

ノルム

以下紹介する命題の証明は、桂田 [6] で読める。
Rn は自然に R 上のベクトル空間であるが、任意の p ∈ [1, ∞) に対して∥x∥p :=( n∑)1/p|xj |p(x = (x1 , · · · , xn )T ∈ Rn ),j=1

また p = ∞ に対して∥x∥p = ∥x∥∞ := max |xj | (x = (x1 , · · · , xn )T ∈ Rn )1≤j≤nとおくと (本当は n にも依存するわけだが、混乱は起らないと考えられるので、単に ∥ · ∥p と書くことにする)、任意の 1 ≤ p ≤ ∞ に対して ∥ · ∥p は Rn のノルムとなる。以下 ∥ · ∥p をRn の p ノルムと呼ぶことにする。

9

実数を成分とする m 行 n 列の行列全体の集合を Rm×n と書く。m = n ならば Rm×n は多元環、m ̸= n であっても Rm×n はベクトル空間の構造を持つ。
A ∈ Rm×n は線型写像 Rn ∋ x 7→ Ax ∈ Rm を引き起こすが、Rn の任意のノルム ∥ · ∥, Rmの任意のノルム ∥ · ∥′ に対して、∥A∥′′ :=∑ ∥Ax∥′

=∥x∥

x̸=0

∑∥Ax∥′ =∥x∥=1∑∥Ax∥′(A ∈ Rm×n )

∥x∥≤1

とおくと、

∥A∥′′ < ∞ (A ∈ Rm×n )となり、∥ · ∥′′ は Rm×n のノルムとなる。これを ∥ · ∥, ∥ · ∥′ から誘導される作用素ノルムと呼ぶ。
任意の A ∈ Rn×n に対して、A のスペクトル半径 (A の固有値の絶対値の最大値) を r(A)と書く。


補題 A.1 (作用素ノルム≧スペクトル半径) ∥ · ∥ を Rn×n の任意の作用素ノルムとする
とき、∥A∥ ≥ r(A) (A ∈ Rn×n ).

R に p ノルムを与えたとき、Rn∥A∥ :=n×nの作用素ノルム∥Ax∥px∈Rn \{0} ∥x∥psup(A ∈ Rn×n )を Rn×n の p ノルムと呼ぶことにして、やはり ∥ · ∥p という記号で表わす。
命題 A.2 (行列の p ノルム) 任意の n ∈ N, A = (aij ) ∈ R∥A∥1 = maxn
∑

1≤j≤n

∥A∥2 =√|aij | ,∥A∥∞ = max√√

i=1r(AAT ) =1≤i≤nr(AT A) ≤n×nn∑

に対して、|aij | ,j=1

∥A∥1 ∥A∥∞ .

特に A が対称ならば∥A∥1 = ∥A∥∞ ,∥A∥2 = r(A).

二つ目の等式は、より一般に A が正規行列 (AAT = AT A) ならばなりたつ。
A.2

Rayleigh 商と min-max 原理

桂田 [7] を見よ。
前項に引き続き行列 A のスペクトル半径を r(A) と書くことにする。
A ∈ Rn×n , x ∈ Rn , x ̸= 0 とするとき、RA (x) :=(Ax, x)
xT Ax= T(x, x)x x
10を A の x に対する Rayleigh 商と呼ぶ。(RA (x) という記号はここだけのものである。)• x が A の固有値 λ に属する固有ベクトルならば、RA (x) = λ である。
• A と x が与えられたとき、∥Ax − λx∥ を最小にする λ は RA (x) である。
• x が A の固有値 λ に属する固有値ベクトルに近いとき、RA (x) は λ の良い近似値となる (ここには書かないが具体的な評価式が得られる)。


命題 A.3 (min-max 原理) A ∈ Rn×n が対称であるとき、A の最大固有値、最小固有値を λmax , λmin とすると、}{(Ax, x)n; x ∈ R \ {0} = [λmin , λmax ] ,(x, x)(Ax, x)= λmax ,x∈Rn \{0} (x, x)
sup(Ax, x)= λmin ,x∈R \{0} (x, x)
infn(Ax, x)= r(A).
x∈Rn \{0} (x, x)
supより詳しく、A の固有値を (重複度の分だけ繰り返して並べて) λ1 ≥ λ2 ≥ · · · ≥ λn とす
ると、minmax RA (x) = 小さい方から k 番目の固有値 = λn−k+1maxmin RA (x) = 大きい方から k 番目の固有値 = λk

dim V =k x∈V \{0}dim V =k x∈V \{0}(1 ≤ k ≤ n),(1 ≤ k ≤ n).

ただし “dim V = k” は「Rn のすべての k 次元線型部分空間 V について」という意味である。
同じことを次のように表現してある本も多い。

(1 ≤ j ≤ n),minmax RA (x) = 大きい方から j 番目の固有値 = λjmaxmin RA (x) = 小さい方から j 番目の固有値 = λn−j+1

dim S=j−1 x∈S ⊥ \{0}dim S=j−1 x∈S ⊥ \{0}(1 ≤ j ≤ n).

(厳密には命題の後半を min-max 原理という。証明については例えば？？？を見よ。)

A.3

正値対称行列, Cholesky 分解

Cholesky 分解については桂田 [5] を見よ。
A = (aij ) ∈ M n×n が上三角行列であるとは、対角線より下にあるすべての成分が 0, つまりi>j=⇒を満たすことをいう。

11aij = 0

命題 A.4 (正値対称行列の Cholesky 分解の存在) A ∈ Rn×n が正値対称行列ならば

A = UT U(7)を満たす上三角行列 U が存在する。特に U の対角成分がすべて 0 であるものに限ると Uは一意的である。


与えられた A に対して、(7) を満たす U を求めることを「A を Cholesky 分解する」という。


命題 A.5 (正値対称行列の平方根) A ∈ Rn×n

が正値対称行列ならば

A = B2を満たす正値対称行列 B が存在する。
証明



A が実対称であるから、ある実直交行列 P を取るとP T AP = diag(λ1 , λ2 , · · · , λn ) (λi は A の固有値)となることはよく知られている。

A = P diag(λ1 , λ2 , · · · , λn )P Tであるが、B := P diag(√λ1 ,√λ2 , · · · ,√)λn P Tとおくと、B は明らかに実対称で、P T P = I に注意すると(√ √(√ √√ )√ )2T
B = P diag
λ1 , λ2 , · · · , λn P P diag
λ1 , λ2 , · · · , λn P T(√ 2 √ 2√ 2)= P diag
λ1 , λ2 , · · · , λn P T = P diag(λ1 , λ2 , · · · , λn )P T = A.

A.4misc

命題 A.6 (帯行列の積は帯行列) A, B は n 次正方行列で、それぞれ半バンド幅が mA ,mB であるとする。つまり|i − j| > mA=⇒|i − j| > mBaij = 0,=⇒bij = 0.

このとき AB の半バンド幅は mA mB 以下である。すなわち

|i − j| > mA mB=⇒(AB)ij = 0.


あれ mA + mB かな？証明を書け。

12

B

文献について一般化固有値問題のアルゴリズムについて、参考になったのは、村田 [1], 矢川・青山 [2],LAPACK USERS’ GUIDE だったりする。古い本が冴えないのが多い中、村田先生の本が役に立つのはさすがだなあ、と思う。
一般化固有値の絶対値の上界評価については、とにかく渡部・山本・中尾充宏 [3] が便利である。
ベクトル、行列のノルムについては、桂田 [6] が便利である。行列の (標準) 固有値問題のアルゴリズムについては、桂田 [7] とその続編 (将来マージ予定) の桂田 [8] を見よ。正値対称
行列の Cholesky 分解については、桂田 [5] を見よ。

参考文献
[1] 村田 健郎, 線形代数と線形計算法序説, サイエンス社 (1986).
げんき[2] 矢川 元基・青山裕司, 有限要素固有値解析, 森北出版 (2001).
[3] 渡部善隆・山本野人・中尾充宏, 一般化固有値問題の精度保証付き計算とその応用, 応用
数理学会論文誌, VOL.9 NO.3 SEPTEMBER (1999), pp.137–150.
な お 、http://www.cc.kyushu-u.ac.jp/scp/system/library/PROGRAM_LIBRARIES/Nakao/vpgep/vpgep.pdf も参照せよ。
[4] G. H. Golub and C. F. Van Loan, Matrix Computations, third edition, The Johns Hopkins
University Press (1996).
[5] 桂田 祐史, Cholesky 分解ノート, http://nalab.mind.meiji.ac.jp/~mk/labo/text/
cholesky.pdf においてある。
[6] 桂田 祐史, 連立 1 次方程式 III 「自分が使うための線形代数」, http://nalab.mind.
meiji.ac.jp/~mk/labo/text/linear-eq-3.pdf においてある。
[7] 桂田 祐史, 行列の固有値問題, http://nalab.mind.meiji.ac.jp/~mk/labo/text/eigenvalues.pdf においてある。
[8] 桂田 祐史, 固有値問題ノートの補足, http://nalab.mind.meiji.ac.jp/~mk/labo/text/eigenvalues-add.pdf においてある。
[9] 桂田 祐史, 一般化固有値問題, http://nalab.mind.meiji.ac.jp/~mk/labo/text/generalized-eigenvalue-problem.pdf においてある。
この文書である。
[10] 桂田 祐史, 『渡部・山本・中尾の一般化固有値問題の論文 (1999) を読む』, http://nalab.
mind.meiji.ac.jp/~mk/labo/members/reading-wyn1999.pdf
渡部善隆・山本野人・中尾充宏 [3] の読書ノート。
13[11] N. Yamamoto and M. T. Nakao, Numerical verifications of solutions for elliptic equations
in non-convex polygonal domains, Numerische Mathematik, Vol.65 (1993), pp. 503–521.
[12] T. Yamamoto, Error bounds for computed eigenvalues and eigenvectors, Numerische
Mathematik, Vol.34 (1980), pp. 189–199.
[13] T. Yamamoto, Error bounds for computed eigenvalues and eigenvectors II, Numerische
Mathematik, Vol.40 (1982), pp. 201–206.
[14] S. M. Rump, Verification Methods for Dense and Sparse Systems of equations, in Jürgen
Herzberger (Editor), Topics in Validated Computations (Proceedings of the IMACSGAMM International Workshop on Validated Computation, Oldenburg, Germany, August 30th–September 3rd, 1993), Elsevier Science, North Holland, Amsterdam, 1994,
pp. 63–135.

14