行列の固有値問題

桂田祐史

2002 年 7 月 16 日, 2006 年 5 月 21 日

大幅な書き直しをしたい…が暇がない。困った。困った。この文書は

  http://nalab.mind.meiji.ac.jp/~mk/labo/text/

で読めます。最近

「固有値問題ノートの補足」

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

というのも書き出しました。

多分骨組みはもう一度 0 から考えて、それに上の二つの文書から必要なものを抜き出して貼り付けていくのでしょう。

目次

1 一般的な注意 2 1.1 問題の定式化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 線形代数の復習 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 心構え (どう立ち向かうべきか) . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2 解法についての概観 4 2.1 相似変換 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Jacobi 法 — 温故知新 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.3 巾乗法 (power method) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.4 固有ベクトルから固有値を求める方法 . . . . . . . . . . . . . . . . . . . . . . 8 2.5 逆反復法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.6 シフト法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.7 Wielandt の減次 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.7.1 対称の場合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.7.2 非対称の場合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3 三重対角化の手法 11

3.1 Householder 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.2 Lanczos (ランチョス) 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.3 実験プログラムたたき台 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

4 二分法 18 4.1 伝統的な説明 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 4.2 Sylvester の慣性律による説明 . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

5 QR 法 21 5.1 正則行列の QR 分解 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 5.2 QR 変換 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 5.3 QR 法の原理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 5.4 QR 分解のアルゴリズム . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

5.4.1 修正 Gram-Schmidt の方法 . . . . . . . . . . . . . . . . . . . . . . . . 27 5.4.2 Givens 変換を用いる方法 . . . . . . . . . . . . . . . . . . . . . . . . . 27 5.4.3 鏡映変換を用いる方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

5.5 QR 法の歴史 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 5.6 Octave での実験 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

6 特異値 28

7 固有値計算のためのパッケージ 28

A 参考書 29

B Strum の方法 29 B.1 スツルムの定理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 B.2 ユークリッドの互除法による Strum 列の生成 . . . . . . . . . . . . . . . . . . 32 B.3 3 重対角行列の固有多項式と Strum 列 . . . . . . . . . . . . . . . . . . . . . . 33 B.4 直交多項式の作る Strum 列 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 B.5 一般化された Strum 列 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

1 一般的な注意

• 次元の低い場合を除いて直接法はない (固有方程式と同値で、それは代数方程式なので、四則と巾根では解けない) 反復法になる。

• 対称な問題と非対称な問題の差は大きい

• 中継地点 (三重対角行列、 Hessenberg 形) を経由すべし

1.1 問題の定式化

いま A を N 次正方行列とする。この時

Ax = λx, x ≠ 0

を満たす λ ∈ C を A の固有値 (eigenvalue), x ∈ C^{N} \ {0} を A の固有値 λ に属する固有ベクトル (eigenvector) と呼ぶ。

固有値問題 — 固有値、固有ベクトルを求める問題 — は非線形問題であり、有限回の四則演算では解けない (N ≥ 5 のときは巾根を求める操作を用いても解けない)。これを解くには、何らかの意味での反復法が必要である。

これに類似した問題に、一般化固有値問題と特異値問題がある。ここでは名前をあげておくだけにとどめるが、将来問題に出会った時に、固有値問題の親戚であると気がつけば良い。いずれの問題も、固有値問題のアルゴリズムを修正したもので解くことが出来る。

一般化固有値問題行列 A, B が与えられた時、方程式

Ax = λBx, x ≠ 0

を満たすスカラー λ, ベクトル x を求める問題を一般化固有値問題という。これは 2 次形式の固有値問題などに関係して現れる。 B = I (単位行列) の場合は通常の固有値問題になる。応用上は B は多くの場合、正定値対称行列になる。 B が正則である場合、上の方程式に B^{−1} をかけて A′ = B^{−1}A とおくと

A′x = λx, x ≠ 0

となって、形式的には通常の固有値問題に帰着するが、このやり方は大抵の場合、得策ではない (A, B が対称であっても、 A′ が対称であるとは限らなくなるなどの理由がある)。簡単なことは次の文書に書いておいた。

『一般化固有値問題メモ』

     http://nalab.mind.meiji.ac.jp/~mk/labo/text/generalized-eigenvalue-problem.
     pdf

特異値問題一言でいうと正方行列でないような行列に対して固有値問題を一般化したものである。

1.2 線形代数の復習

線形代数学で学んだことをいくつか復習しておこう。

[figure]

定理 1.1 固有値は、固有方程式と呼ばれる代数方程式

det(λI − A) = 0

の根である。従って N 次正方行列の固有値は重複度を込めて数えると N 個ある。

[figure]

逆に任意の代数方程式に対して、それを固有方程式にもつ行列が存在するので、数学的には「固有値問題は代数方程式の問題と同値である」。

[figure]

定理 1.2 (i) Hermite 行列の固有値は実数であり、ユニタリ行列で対角化できる。

(ii) 実対称行列の固有値は実数であり、直交行列で対角化できる。

[figure]

1.3 心構え (どう立ち向かうべきか)

• 固有値を求めるのに、固有方程式を解こうとするのは得策ではない^{1}。固有値問題を解くには、以下に説明する固有値問題用の解法を採用すべきである。

• 実際的な観点からは、固有値問題の解法に使用される各種方法の原理を理解し、自分が解こうとしている問題にあった方法 or プログラムを選択できるようになることを目指すのが良い。

• 解こうとしている問題については、以下のことに留意して考えよう。

(i) 問題は対称 (=行列が実対称行列または Hermite 行列) かどうか。 (ii) すべての固有値が欲しいのかどうか。 (iii) 固有ベクトルは必要かどうか。

(ii), (iii) は要するに「ムダはやめよう」ということである。

• 対称な問題は非対称な問題と比べて解きやすい。

「対称な問題を解くのはサイエンスであるが、非対称な問題を解くのはアートである。」

という言葉があるくらいである。

2 解法についての概観

歴史的なことを述べると、以前は対称な問題は Jacobi 法と呼ばれるアルゴリズムで解かれるのが普通であった。しかし Jacobi 法は N が小さな (せいぜい数十) である場合には実用的であるが、 N の大きな問題を解くのに採用するのは得策ではない。今では行列の三重対角化や Hessenberg 形への変換を利用した解法が主流である（対称行列の場合は三重対角化、非対称行列の場合は Hessenberg 行列への変換を行なうことになる）。行列を三重対角化 (resp. Hessenberg 化) するとは、相似変換を施して、行列を三重対角行列 (resp. Hessenberg 行列) に変換することをいう。この変換の方法は色々あるが、いずれも O(N^{3}) の演算量で済む。後で述べるように相似変換で固有値は不変なので、変換後の「簡単な」行列の固有値を求めれば、元の行列の固有値が求まったことになる。三重対角行列や Hessenberg 行列に対する固有値問題は、元の行列よりも簡単に解ける、というのが基本的なアイディアである。

1そもそも次数 n が 4 以下でないと根の公式はないし、代数方程式のための数値解法で解ける問題の範囲は、固有値問題のための数値解法で解ける範囲よりも狭い。

色々な細かな技法 (部品) を、利用者の要求に応じて組み合わせることで、望ましい解法が出来上がる。以下、それらの原理を説明する。三重対角化や Hessenberg 化のためのアルゴリズムの解説は次節にまわして、ここでは実際に固有値を求める諸方法を解説する。

2.1 相似変換

固有値問題でもっとも基本的な方法は相似変換である。これは行列 A を正則行列 P によって、 P^{−1}AP に変換することを意味する。

[figure]

定理 2.1 相似変換により固有値は不変である。

[figure]

に注意しよう。 A が実対称な場合には P として実直交行列が使われる。この場合 P^{−1} = ^{t}P であり、計算が簡単になることの他に、様々な利点がある。

実対称行列の、実直交行列による相似変換は実対称行列である:

(P^{T} AP)^{T} = P^{T} A^{T} (P^{T})^{T} = P^{T} AP.

実直交行列による変換として、次の三つが重要である。

(1) Givens 変換 — 2 次元平面における回転

(2) 鏡映変換 (Householder 変換) — 超平面に関する対称移動

(3) QR 変換

2.2 Jacobi 法 — 温故知新

1960 年代まで主流であった Jacobi 法について簡単に説明しよう。これは実対称行列を 2 次元の回転変換により変形していって、対角成分以外の成分の絶対値を小さくしていくというものであり、行列のサイズ N が 10 程度の小さなものであれば現在でも実用的である。

なぜ「対角成分以外の成分の絶対値を小さくしていく」のか？これについて説明しよう。まず次の定理は簡単であるが重要である。

[figure]

定理 2.2 対角行列の固有値は対角成分である。

[figure]

一般の行列の場合も次の定理が成り立つ。おおざっぱにまとめると「行列の固有値は対角成分に近く、そのずれは非対角成分の大きさによる」。

[figure]

定理 2.3 (Gerschgorin の円板定理 (1931)) A = (a_{ij}) に対して

Δ_{i} = {z ∈ C; |z − a_{ii}| ≤ ∑_{1≤j≤N,j≠i} |a_{ij}|} (i = 1, 2, ⋯, N)

とおくと、

σ(A) ≡ A の固有値全体⊂ ⋃_{i=1}^{n} Δ_{i}.

もし Δ_{i} のうちで k 個が連結成分をなせば、その中に k 個の固有値がある。

[figure]

証明 Ax = λx, x = (x_{1}, …, x_{N}) ≠ 0 とする。今 x の絶対値最大の成分を x_{k} とする:

max_{i=1,2,...,N} |x_{i}| = |x_{k}|.

方程式 Ax = λx の第 k 成分を書くと

∑_{j=1}^{N} a_{kj}x_{j} = λx_{k}.

これから

(λ − a_{kk})x_{k} = ∑_{1≤j≤N,j≠k} a_{kj}x_{j}.

(以下 1 ≤ j ≤ N, j ≠ k を単に j ≠ k と書くことにする。 ) よって

|λ − a_{kk}| ≤ ∑_{j≠k} |a_{kj}| |\frac{x_{j}}{x_{k}} | ≤ ∑_{j≠k} |a_{kj}|,

したがって λ ∈ Δ_{k} ⊂ ⋃_{i=1}^{N} Δ_{i}.

つぎに

D = diag(a_{11}, ⋯, a_{NN}), A_{t} = (1 − t)D + tA (t ∈ [0, 1])

とおくと、 A_{t} に対する円板は

Δ_{i}(t) = {z ∈ C; |z − a_{ii}| ≤ t ∑_{j≠i} |a_{ij}|}

である。 A_{i} の固有値は t について連続的に変化し、 t = 0 のとき Δ_{i}(0) = {a_{ii}} にそれぞれ重複個数だけある。 t を増加させれば、いくつかの円板が拡大して合流するが、 k 個の Δ_{i}(t) が合流すれば、その中に A_{t} の k 個の固有値が存在する。ゆえに t = 1 に達したとき、 Δ_{i} の k 個が連結成分をなせば、その中に A_{1} = A の k 個の固有値が存在する。

[figure]

系 2.4 (狭義優対角行列は正則)N 次正方行列 A = (a_{ij}) について、

|a_{ii}| > ∑_{1≤j≤N,j≠i} |a_{ij}| (i = 1, 2, ⋯, N)

が成り立てば、 A は正則である。

[figure]

証明 0 ∉ ⋃_{i=1}^{N} Δ_{i} より、 0 は A の固有値ではない。ゆえに A は正則である。

なぜか今のカリキュラムの盲点となってしまって、紹介されることがないが、「代数方程式の根は方程式の係数の連続関数である」行列の固有方程式の係数は、行列の成分の連続関数であることは明らかだから、「行列の固有値は行列の成分の連続関数である」。 Gerschgorin の定理はあざやかに連続性を見せてくれる。

2.3 巾乗法 (power method)

ここでも問題は対称であると仮定する。

与えられた行列の絶対値最大の固有値を求めるための、累乗法あるいは巾乗法^{2}を説明する。行列 A の固有値 {λ_{i}; i = 1, 2, ⋯, n} は絶対値の順に番号づけられているとする:

|λ_{1}| ≥ |λ_{2}| ≥ ⋯ |λ_{n}|.

また、 {u_{i}; i = 1, 2, ⋯, n} は {λ_{i}} に対応する A の固有ベクトルからなる正規直交基底とする。ここで話を簡単にするために次の仮定をおく。

仮定: |λ_{1}| > |λ_{2}|.

この時、適当な x_{0} (x_{0} ∉ Span{u_{2}, u_{3}, ⋯, u_{n}}) を選んで

\{ \begin{matrix}y_{k+1} & := & Ax_{k} \\ x_{k+1} & := & y_{k+1}/∥y_{k+1}∥ \end{matrix}

によりベクトル列 {x_{k}; k = 0, 1, ⋯ } を定めると、

lim_{k→∞} x_{k} = ±u_{1}

となる。実際には十分大きな番号 k を取れば、 x_{k} = ±u_{1} とみなしてよい。

大雑把な説明まず適当な {c_{i}} を選ぶことにより

x_{0} = ∑_{i=1}^{n} c_{i}u_{i}

2 「羃乗法」あるいは「冪乗法」と書くのが正しい？

と展開されることに注意する。ここで c_{1} ≠ 0 である。

A^{k}x_{0} = ∑_{i=1}^{n} c_{i}A^{k}u_{i} = ∑_{i=1}^{n} c_{i}λ_{i}^{k}u_{i} = c_{1}λ_{1}^{k} \{u_{1} + ∑_{i=2}^{n} (\frac{λ_{i}}{λ_{1}} )^{k} c_{i}u_{i}\} .

ここで x_{k} = A^{k}x_{0}/∥A^{k}x_{0}∥ であることに注意すると、 x_{k} → ±u_{1} (k → ∞) となることが分かる。

2.4 固有ベクトルから固有値を求める方法

誤差がなければ、 Ax = λx の適当な (x_{i} ≠ 0 となる i に対する) 成分に対応する方程式

∑_{j=1}^{n} a_{ij}x_{j} = λx_{i}

の両辺を x_{i} で割れば λ が求まるが、 x が近似的な固有ベクトルでしかない場合には、以下に解説する Rayleigh^{3} 商を用いる方法の方が良い。

[figure]

定義 2.5 (Rayleigh 商) N 次正方行列 A と、 x ∈ C^{N} \ {0} に対して

\frac{(Ax, x)}{(x, x)}

を A の x に対する Rayleigh 商という。

[figure]
[figure]

補題 2.6 N 次正方行列 A と x ∈ C^{N} \ {0} が与えられたとき、 ∥Ax − λx∥ を最小にする λ は A の x に対する Rayleigh 商である。

[figure]

証明

∥Ax − λx∥^{2} = (Ax, Ax) − λ(Ax, x) − \overline{λ}(x, Ax) + |λ|^{2}∥x∥^{2}
= ∥x∥^{2} [λ − \frac{(x, Ax)}{∥x∥^{2}}] [\overline{λ} − \frac{\overline{(x, Ax)}}{∥x∥^{2}}] + ∥Ax∥^{2} − \frac{|(x, Ax)|^{2}}{∥x∥^{2}}
= ∥x∥^{2} |λ − \frac{(x, Ax)}{∥x∥^{2}} |^{2} + ∥Ax∥^{2} − \frac{|(x, Ax)|^{2}}{∥x∥^{2}} .

x を A の近似固有ベクトルとする時、 Rayleigh 商

λ′ = \frac{(Ax, x)}{(x, x)}

は x に対応する固有値の良い近似になる。粗く言って

固有値の誤差= O(固有ベクトルの誤差^{2})

3Lord Rayleigh — アルゴンの発見などの業績のある物理学者。

が成り立つ。詳しくは次の命題を見よ。

[figure]

命題 2.7 A が正規行列 (AA^{∗} = A^{∗}A) で、その固有値を λ_{1}, ⋯, λ_{N} とし、それに対応する固有ベクトルからなる正規直交基底を x_{1}, x_{2}, ⋯, x_{N} とする。いま単位ベクトル v の Rayleigh 商 μ = v^{∗}Av が

∥Av − μv∥ = ε, |λ_{j} − μ| ≥ δ > ε (j = 2, …, N)

を満たせば、次の意味で μ は λ_{1} にごく近い:

|λ_{1} − μ| ≤ \frac{ε^{2}}{δ(1 − ε^{2}/δ^{2})}.

さらに v は次の意味で x_{1} に近い: |c| = 1 である定数 c に対して

∥cv − x_{1}∥^{2} ≤ (ε/δ)^{2} + (ε/δ)^{4}. [figure]

証明省略。一松 [2] あるいは [3] を見よ。

2.5 逆反復法

絶対値が最小の固有値 (上の記号で λ_{n}) を求める方法である。 A^{−1} の固有値は λ_{1}^{−1}, ⋯, λ_{n}^{−1} であり、絶対値最大のものは λ_{n}^{−1} となるから、 A^{−1} に巾乗法を適用することにより、 λ_{n}^{−1} が求まり、その逆数を取れば λ_{n} が得られる。

計算上の注意: 反復の各段階で

y_{k+1} := A^{−1}x_{k}

という計算が必要になるが、これは

Ay_{k+1} = x_{k}

という y_{k+1} に関する連立一次方程式を解くことにより実行する (いつものことであるが、 A^{−1} を計算するのは馬鹿馬鹿しい)。

2.6 シフト法

既に得られている近似固有値の精度を改良したい、あるいは、ある指定した値に最も近い固有値を求める方法である。

行列 A の固有値 λ_{i} に対する近似固有値 λ′_{i} が分かっているとしよう。この時、

A′ = A − λ′_{i}I

の固有値は

λ_{1} − λ′_{i}, λ_{2} − λ′_{i}, ⋯, λ_{n} − λ′_{i}

となる。 λ′_{i} は λ_{i} の近似値ということで、絶対値が最小なのは λ_{i} − λ′_{i} であると期待できる。よって、 A′ に対して逆反復法を適用すれば、この値 (それを Δλ とおこう) が高精度に計算できる。こうして

λ_{i} := λ′_{i} + Δλ

により λ_{i} が求まる。

2.7 Wielandt の減次

行列 A ∈ C^{n×n} の固有値 λ_{i} (i = 1, 2, ⋯, n) について、

|λ_{1}| > |λ_{2}| > |λ_{3}| ≥ ⋯ ≥ |λ_{n}|

となっていると仮定する。

2.7.1 対称の場合

最初に簡単のため、 A が対称 A^{∗} = A の場合を考える。冪乗法により、 λ_{1} と対応する固有ベクトル z_{1} を求めたとする。 λ_{i} に属する固有ベクトルを z_{i} として (i = 2, 3, ⋯, n)、 z_{1}, ⋯, z_{n} が正規直交系 (z_{i}^{∗}z_{j} = δ_{ij}) であるとする (ただし i ≥ 2 について、具体的に求めるわけではない)。このとき、

B = A − λ_{1}z_{1}z_{1}^{∗}

とおくと、これも対称で、

Bz_{1} = Az_{1} − λ_{1}z_{1}z_{1}^{∗}z_{1} = λ_{1}z_{1} − λ_{1}z_{1} · 1 = 0,
Bz_{i} = Az_{i} − λ_{1}z_{1}z_{1}^{∗}z_{i} = λ_{i}z_{i} − λ_{1}z_{1} · 0 = λ_{i}z_{i} (i = 2, 3, ⋯, n).

ゆえに B について冪乗法を適用すると、 λ_{2} が得られるはずである。 B と任意に与えられたベクトル v との積は

Bv = Av − λ_{1}z_{1}(z_{1}^{∗}v)

と計算できるので、 B の成分を求めておく必要はないことに注意しよう。

2.7.2 非対称の場合

A が対角化可能としておく (これは余計な仮定か？ )。つまり、 λ_{i} に属する固有ベクトル z_{i} が存在し、 z_{1}, ⋯, z_{n} が基底となると仮定する。

よく知られているように、\overline{λ}_{1} は、 A の Hermite 共役 A^{∗} の固有値になるが、その固有ベクトル w_{1} をシフト法で求めておく (もちろん A が「対称」ならば λ_{1} = λ_{1}^{∗} で、 w_{1} = z_{1} と取れるので、対称の場合にやった操作の一般化になっている)。

必要ならば定数倍することによって、 w_{1}^{∗}z_{i} = δ_{i1} となるようにできる。実際、 A^{∗}w_{1} = \overline{λ}_{1}w_{1} に注意して、

w_{1}^{∗}z_{i} = (z_{i}, w_{1}) = (z_{i}, \overline{λ}_{1}^{−1}A^{∗}w_{1}) = λ_{1}^{−1}(Az_{i}, w_{1}) = λ_{1}^{−1}λ_{i}(z_{i}, w_{1})

であるから、

(1 − \frac{λ_{i}}{λ_{1}} ) w_{1}^{∗}z_{i} = 0.

ゆえに λ_{i} ≠ λ_{1} ならば

w_{1}^{∗}z_{i} = 0.

一方 w_{1}^{∗}z_{1} ≠ 0 である。実際、 w_{1} を w_{1} = c_{1}z_{1} + ⋯ + c_{n}z_{n} と展開すると (ここで {z_{i}} が基底であることを使っている)、直交性 w_{1}^{∗}z_{i} = 0 (i = 2, 3, ⋯, n) によって、 0 ≠ w_{1}^{∗}w_{1} = c_{1}w_{1}^{∗}z_{1} となるから。

そこで B = A − λ_{1}z_{1}w_{1}^{∗} とおくと、

Bz_{1} = Az_{1} − λ_{1}z_{1}w_{1}^{∗}z_{1} = λ_{1}z_{1} − λ_{1}z_{1} = 0

また

Bz_{i} = Az_{i} − λ_{1}z_{1}w_{1}^{∗}z_{i} = λ_{i}z_{i} − λ_{1}z_{1}0 = λ_{i}z_{i}.

ゆえに B について冪乗法を適用すると、 λ_{2} が得られる。

疑問点 A が対角化可能でない場合はどうなるだろう。 w_{1}^{∗}z_{1} ≠ 0 を保証するためだけにこの仮定をおくのは大げさのように思う。 A の固有ベクトルが n 個存在しない場合には、 B の固有ベクトルとして何が出て来るかすぐには分からないということはある。

3 三重対角化の手法

実対称行列に対する三重対角化の手法を学ぶ。

(1) Givens 法

(2) Householder 法

(3) Lanczos 法

与えられた実対称行列 A = (a_{ij}) に対し、適当な正則行列 P を求めて

P^{−1}AP = 三重対角

とする。特に P として直交行列を取る。

歴史的には、 Jacobi 法の変形として Givens 法が最初に現われたが、その後現われた Householder 法は演算回数が約半分となるなど利点が多く、普及している。大型疎行列に対しては Lanczos 法が有力な方法とされている。

3.1 Householder 法

u を R^{n} の単位ベクトルとする。 u に直交する超平面

W_{u} ≡ {x ∈ R^{n}; (x, u) = 0}

に関する対称移動を表す変換を U とすると、

U = I − 2uu^{T}

である。 (\overrightarrow{OP} = x なる点 P から W_{u} に下ろした垂線の足を Q とすると、\overrightarrow{QP} = (x, u)u. それゆえ、 Ux = x − 2Q⃗P = x − 2(x, u)u = (I − 2uu^{T})x. ゆえに U = I − 2uu^{T}.)

[figure]

定義 3.1 (鏡映変換、 Householder 行列)R^{n} の単位ベクトル u に対し、行列 U = I − 2uu^{T} で定まる線形変換をベクトル u に対する鏡映変換、 U を Householder 行列 (または基本直交行列) と呼ぶ。

[figure]
[figure]

補題 3.2 (鏡映変換の性質)U を Householder 変換とするとき、

(i) U は直交変換である。

(ii) U は対称変換である。

[figure]

証明

(i) UU^{T} = I を計算で確かめてもよいし、 U が長さを変えないことから明らかである。

(ii) これは U = I − 2uu^{T} であることから分かる。

参考 N 次元の任意の直交変換は、 N 個以下のベクトルに関する鏡映変換の積に書ける (Cartan)。

[figure]

補題 3.3 ∥x∥ = ∥y∥, x ≠ y なる x, y ∈ R^{n} に対して、

u = \frac{x − y}{∥x − y∥}, U = I − 2uu^{T}

とおくと Ux = y. [figure]

証明明らか。 Householder 行列による三重対角化

A_{0} = A,

A_{1} = U_{0}A_{0}U_{0} = [ \begin{matrix}∗ & ∗ & 0 & ⋯ & 0 \\ ∗ & ∗ & ∗ & ⋯ & ∗ \\ 0 & ∗ & ∗ & ⋯ & ∗ \\ ⋮ & ⋮ & ⋮ & ... & ⋮ \\ 0 & ∗ & ∗ & ⋯ & ∗ \end{matrix} ] , U_{0} = I − 2u_{0}u_{0}^{T} ,

A_{2} = U_{1}A_{1}U_{1} = [ \begin{matrix}∗ & ∗ & 0 & 0 & ⋯ & 0 \\ ∗ & ∗ & ∗ & 0 & ⋯ & 0 \\ 0 & ∗ & ∗ & ∗ & ⋯ & ∗ \\ 0 & 0 & ∗ & ∗ & ⋯ & ∗ \\ ⋮ & ⋮ & ⋮ & ⋮ & ... & ⋮ \\ 0 & 0 & ∗ & ∗ & ⋯ & ∗ \end{matrix} ] , U_{1} = I − 2u_{1}u_{1}^{T} ,

· · · · · · · · ·

A_{N−2} = U_{N−3}A_{N−3}U_{N−3} = [ \begin{matrix}∗ & ∗ & 0 & 0 & ⋯ & 0 & 0 & 0 & 0 \\ ∗ & ∗ & ∗ & 0 & ⋯ & 0 & 0 & 0 & 0 \\ 0 & ∗ & ∗ & ∗ & ... & 0 & 0 & 0 & 0 \\ 0 & 0 & ∗ & ∗ & ... & 0 & 0 & 0 & 0 \\ ⋮ & ⋮ & ... & ... & ... & ... & ... & ⋮ & ⋮ \\ 0 & 0 & 0 & 0 & ... & ∗ & ∗ & 0 & 0 \\ 0 & 0 & 0 & 0 & ... & ∗ & ∗ & ∗ & 0 \\ 0 & 0 & 0 & 0 & ⋯ & 0 & ∗ & ∗ & ∗ \\ 0 & 0 & 0 & 0 & ⋯ & 0 & 0 & ∗ & ∗ \end{matrix} ] , U_{N−3} = I − 2u_{N−3}u_{N−3}^{T}

のように N − 2 回の Householder 変換で三重対角化する。この 1 ステップ (A_{k−1} から A_{k} を作るところ) を発見的に説明しよう。面倒なので A_{k−1} を単に A と書き、

A = [ \begin{matrix}C & \frac{0}{b^{T}} \\ _{0 b} & B \end{matrix} ]

とブロック分けしておく。 u = u_{k} の最初の k 成分を 0 とおく、すなわち

u = ( 0 ; ⋮ ; 0 ; ∗ ; ⋮ ; ∗ ) = ( 0 ; ⋮ ; 0 ; v ) , v ∈ R^{N−k}

とすると

U = I_{N} − 2uu^{T} = [ \begin{matrix}I_{k} & O \\ O & Q \end{matrix} ] , Q = I_{N−k} − 2vv^{T}

となるので、

新A = UAU

とすると対応するブロックは

新C = C,
新B = QBQ,
新b = Qb

目標は

新b = Qb = ( s ; 0 ; ⋮ ; 0 )

の形にすることである。補題 3.3 を考えると s = ±∥b∥ とすればよい。

s の符号= −b の第 1 成分 (−b_{1}) の符号

となるように選ぶと桁落ちが起こりにくい。まとめると

(♡_{k}) \{ \begin{matrix}s & = & −sign (b_{1}) · ∥b∥, \\ 新 b & = & ( s ; 0 ; ⋮ ; 0 ) , \\ v & = & \frac{b − 新 b}{∥b − 新 b∥}, \\ Q & = & I_{N−k} − 2vv^{T} \end{matrix}

全体としては

[figure]

for k := 1 to N − 2 do begin

第 k 列の第 k + 1 成分以降を b とする。 (♡_{k}) により s, 新 b, v, Q を作る。新 B = QBQ を計算する。

end

[figure]

色々な工夫が可能である。まとめると

\{ \begin{matrix}s & = & −sign (b_{1}) × ∥b∥ \\ w & = & b − 新b = 第 1 成分は旧 b の第 1 成分 −s, その他の成分はそのまま \\ ∥w∥^{2} & = & 2\{s^{2} − (b の第一成分) × s\} \\ p & = & \frac{Bw}{(∥w∥^{2}/2)} \\ α & = & w^{T} p/∥w∥^{2} \\ q & = & p − αw \\ 新B & = & B − wq^{T} − qw^{T} (対称なので半分の計算で OK) \end{matrix}

注: この文書の過去の版で、 ∥w∥^{2} の計算式の符号を間違えて

∥w∥^{2} = 2{s^{2} + (b の第 1 成分) × s} (この式は間違っている！ )

のようにしていました。ごめんなさい。

3.2 Lanczos (ランチョス) 法

(この話題については、『固有値問題ノートの補足』 http://nalab.mind.meiji.ac.jp/~mk/
labo/text/eigenvalues-add.pdf に説明を書いた。そのうちこちらとマージする。 )

実対称行列 A が直交行列 P で三重対角化されたとする:

B = P^{T} AP = ( \begin{matrix}α_{1} & β_{1} & 0 & & 0 \\ β_{1} & α_{2} & β_{2} \\ & ... & ... & ... \\ & & β_{N−2} & α_{N−1} & β_{N−1} \\ 0 & & 0 & β_{N−1} & α_{N} \end{matrix} ) .

これから AP = P B であるが、 P = (u_{1}u_{2} ⋯ u_{N}) とおくと、

A(u_{1}u_{2} ⋯ u_{N}) = (u_{1}u_{2} ⋯ u_{N}) ( \begin{matrix}α_{1} & β_{1} & 0 & & 0 \\ β_{1} & α_{2} & β_{2} \\ & ... & ... & ... \\ & & β_{N−2} & α_{N−1} & β_{N−1} \\ 0 & & 0 & β_{N−1} & α_{N} \end{matrix} )

となるから

Au_{1} = α_{1}u_{1} + β_{1}u_{2}
Au_{2} = β_{1}u_{1} + α_{2}u_{2} + β_{2}u_{3}
⋯
Au_{i} = β_{i−1}u_{i−1} + α_{i}u_{i} + β_{i}u_{i+1}
⋯
Au_{N} = β_{N−1}u_{N−1} + α_{N}u_{N}

第 k 行と u_{k} の内積を作ると

(1)  α_{k} = u_{k}^{T} Au_{k}.

また

(2)  v_{k+1} := \{ \begin{matrix}Au_{1} − α_{1}u_{1} & (k = 0) \\ Au_{k} − (β_{k}u_{k−1} + α_{k}u_{k}) & (1 ≤ k ≤ N − 1) \end{matrix}

とおくと、 v_{k+1} = β_{k}u_{k+1}, ∥u_{k+1}∥ = 1 であるから、

(3)  β_{k} = ±∥v_{k+1}∥
(4)  u_{k+1} = \frac{v_{k+1}}{β_{k}} .

これから、最初に適当に ∥u_{1}∥ = 1 となる u_{1} を選び、以下 β_{k} ≠ 0 である限り、計算を進めることができる。

3.3 実験プログラムたたき台

実験用行列生成のヒント

[figure]

% ランダム Hessenberg 行列
function ret = rand_h(n)
  ret = triu(rand(n,n)) + diag(rand(n-1,1),-1);

% ランダム三重対角行列

function ret = rand_t(n)
  ret=diag(rand(n,1),0)+diag(rand(n-1,1),1)+diag(rand(n-1,1),-1);

% ランダム実対称三重対角行列

function ret = rand_st(n)
  u=rand(n-1,1);
  ret=diag(rand(n,1),0)+diag(u,1)+diag(u,-1);

% ランダム実対称行列

function ret = rand_s(n)
  a=rand(n,n);
  ret = (a+a')/2;

% QR 変換 (回数指定)

function QR = qr_iteration(a,n)
  QR = a;
  for i=1:n
    [q r]=qr(QR);
    QR=r*q;
  end

% 下三角部分のノルム

function n = norm_l(a)
  n = norm(a-triu(a));

[figure]

beki.m

[figure]

 1  %

2 % beki.m --- 冪乗法の系統の算法の実験

 3  %
 4
 5  format long
 6
 7  a=[1,1,1;1,2,2;1,2,3]
 8  eigen=eig(a)
 9
10  maxiter=20; r=zeros(maxiter,1);
11
12  disp("冪乗法")
13  x=ones(3,1);
14  for k=1:maxiter
15    y=a*x; x=y/norm(y); r(k)=x'*(a*x);
16  end
17  subplot(1,4,1)
18  semilogy(1:maxiter,abs(r-eigen(3)))

19 % 絶対値最大の固有値に属する固有ベクトルを保存しておく

20  u=x;
21
22  disp("逆反復法")
23  x=ones(3,1);
24  for k=1:maxiter
25    y=a\x; x=y/norm(y); r(k)=x'*(a*x);
26  end
27  subplot(1,4,2)
28  semilogy(1:maxiter,abs(r-eigen(1)))
29
30  disp("シフト法")
31  kinji=0.6;
32  ap=a-kinji*eye(3,3);
33  x=ones(3,1);
34  for k=1:maxiter
35    y=ap\x; x=y/norm(y); r(k)=x'*(ap*x)+kinji;
36  end
37  subplot(1,4,3)
38  semilogy(1:maxiter,abs(r-eigen(2)))
39
40  x=ones(3,1);
41  x=x-(u'*x)*u;
42  for k=1:maxiter
43    y=a*x; x=y/norm(y); x=x-(u'*x)*u; r(k)=x'*(a*x);
44  end
45  subplot(1,4,4)
46  semilogy(1:maxiter,abs(r-eigen(2)))

[figure]

householder.m

[figure]

1  % Householder
2  function a=householder(A)
3    [n,n] = size(A);
4    a = A;
5    for k=1:n-1
6      b = a(k+1:n,k);
7      B = a(k+1:n,k+1:n);
8      s = - sign(b(1)) * norm(b);
9      w = b;
10      w(1) = w(1) - s;
11      normw2 = 2 * s * (s - b(1));
12      p = (B * w) / (normw2/2);
13      alpha = (w'*p) / normw2;
14      q = p - alpha * w;
15      a(k+1:n,k+1:n) = B - w * q' - q * w';
16      a(k+2:n,k) = zeros(n-k-1,1);
17      a(k,k+2:n) = zeros(n-k-1,1)';
18      a(k+1,k) = s;
19      a(k,k+1) = s;
20    end

[figure]

なお、授業のレポートである福澤 [7] も参考になるかも。

4 二分法

二分法 (bisection method, Strum method) を解説する。

• これは実対称行列専用の方法である (中間値の定理を用いるので、固有多項式が実係数で、根がすべて実数でないといけない)。

• 前節で説明した方法を用いて、与えられた実対称行列を三重対角化しておくことで、二分法を効率的に実行することが出来る。 4.1 では最初から三重対角行列に対する説明になっている。

• 固有値は固有多項式の根であるが、 Strum 列の理論によって、ある区間内の固有値の個数を計算することが出来る。このことと、いわゆる二分探索 (binary search) の方法を組み合わせて得られるのが、固有値計算手法としての二分法である。

4.1 伝統的な説明

T = ( \begin{matrix}a_{1} & b_{1} & 0 & & 0 \\ b_{1} & a_{2} & b_{2} \\ & ... & ... & ... \\ & & b_{N−2} & a_{N−1} & b_{N−1} \\ 0 & & 0 & b_{N−1} & a_{N} \end{matrix} )

を実対称三重対角行列とする。 b_{k} ≠ 0 (k = 1, 2, ⋯, N − 1) と仮定する。もしある k に対して b_{k} = 0 ならば

T = ( \begin{matrix}T′ & O \\ O & T′′ \end{matrix} )

となり、 T′, T′′ の固有値を求める問題に帰着できるから、一般性は失わない。 p_{k}(λ) を λI − T の k 次首座小行列式とする (k = 0, 1, ⋯, N)。すなわち

p_{k}(λ) := \{ \begin{matrix}det(λI_{k} − T_{k}) & (k = 1, 2, ⋯, N) \\ 1 & (k = 0). \end{matrix}

ただし、

I_{k} = k 次の単位行列, T_{k} = T の k 次首座小行列= ( \begin{matrix}a_{1} & b_{1} & 0 & & 0 \\ b_{1} & a_{2} & b_{2} \\ & ... & ... & ... \\ & & b_{k−2} & a_{k−1} & b_{k−1} \\ 0 & & 0 & b_{k−1} & a_{k} \end{matrix} ) .

すぐ分かることは、

[figure]

\begin{matrix}補題 4.1 (p_{k}(λ) の計算は簡単) \\ \{ & \begin{aligned}p_{0}(λ) & = & 1 \\ p_{1}(λ) & = & λ − a_{1} \\ p_{k+1}(λ) & = & (λ − a_{k+1})p_{k}(λ) − b_{k}^{2}p_{k−1}(λ) (k = 1, 2, ⋯, N − 1) \\ p_{N}(λ) & = & det(λI − T). \end{aligned} \end{matrix}

[figure]
[figure]

補題 4.2 {p_{k}(λ)}_{k=0}^{N} は Strum 列である。すなわち

(i) p_{0}(λ) は定符号。 (ii) p_{k}(λ) の根は有限個 (k = 1, 2, …, N). (iii) p_{k}(λ), p_{k+1}(λ) は共通根を持たない。 (iv) p_{k}(λ_{0}) = 0 ⟹ p_{k−1}(λ_{0})p_{k+1}(λ_{0}) < 0 (k = 1, 2, ⋯, N − 1).

(v) p_{N}(λ_{0}) = 0 ⟹ p′_{N}(λ_{0})p_{N−1}(λ_{0}) > 0. [figure]

証明 Strum 列については別に一つの章をもうけて説明するので、そこで述べることにする。

符号の変化数 {p_{i}(λ)}_{i=0}^{N} を Strum 列とするとき、 p_{N}(x) ≠ 0 なる x に対して

N(x) ≡^{def.} “{p_{0}(x), p_{1}(x), ⋯, p_{N}(x)}” の符号の変化数

とおく。例えば

\begin{matrix}p_{0}(x) & p_{1}(x) & p_{2}(x) & p_{3}(x) & p_{4}(x) & p_{5}(x) & p_{6}(x) & p_{7}(x) & p_{8}(x) & p_{9}(x) & p_{10}(x) \\ + & − & − & − & 0 & + & + & + & + & − & − \end{matrix}

では N(x) = 3.

この符号の変化数は (特別な注意をせずに) 数値的に安定して計算できる。つまり絶対値が非常に小さくて、符号の判別がつきにくい場合も、「符号の変化数」そのものは疑いがなく計算できる。例えば

\begin{matrix}p_{k−1}(x) & p_{k}(x) & p_{k+1}(x) \\ + & 絶対値小 & − \end{matrix} または\begin{matrix}p_{k−1}(x) & p_{k}(x) & p_{k+1}(x) \\ − & 絶対値小 & + \end{matrix}

において p_{k}(x) が正であっても負であっても 0 であっても符号の変化数の計算にとっては影響がない。注意すべきは Strum 列の条件 (iv) から

\begin{matrix}p_{k−1}(x) & p_{k}(x) & p_{k+1}(x) \\ + & 絶対値小 & + \end{matrix} または\begin{matrix}p_{k−1}(x) & p_{k}(x) & p_{k+1}(x) \\ − & 絶対値小 & − \end{matrix}

のような場合 (もしこうなったら符号の変化数の計算がむつかしい) が起こり得ないことである。

Strum の定理によって、 p_{N}(a)p_{N}(b) ≠ 0 なる [a, b] において、 [a, b] 内の零点の個数は N(a)− N(b) であることが分かる。

4.2 Sylvester の慣性律による説明

二分法を Sylvester の慣性律を用いて説明するのは、森・杉原・室田 [5] による。

線形代数でおなじみの Sylvester の慣性律「行列の符号数は座標変換で変化しない」を復習しよう。正則行列 M に対して

π(M) := M の固有値のうち正のものの個数 ,
ζ(M) := M の固有値のうち 0 のものの個数,
ν(M) := M の固有値のうち負のものの個数

とおく。

[figure]

定理 4.3 (Sylvester の慣性律による説明)A を実対称行列、 S を正則行列で B = S^{T} AS とするとき

π(A) = π(B), ζ(A) = ζ(B), ν(A) = ν(B). [figure]

任意に選んだ σ ∈ R に対して、 A − σI を LDU 分解する、すなわち

A − σI = LDL^{T} , L は下三角行列, D は対角行列.

このとき定理から

π(A − σI) = π(D), ζ(A − σI) = ζ(D), ν(A − σI) = ν(D)

であるが、 π(D), ζ(D), ν(D) はすぐに分かる。従って、行列 A の固有値について、任意に与えられた実数 σ よりも大きいもの、小さいものの個数が勘定できることになる。

A が三重対角行列であれば、 LDU 分解は O(n) の計算量で済ませられることを注意しておく。

5 QR 法

この節の内容は主に一松 [2] による。

5.1 正則行列の QR 分解

A ∈ GL(n; C) の列ベクトルを a_{1}, ⋯, a_{n} とする:

(5)  A = (a_{1} a_{1} ⋯ a_{n}).

この a_{1}, ⋯, a_{n} に Gram-Schmidt の直交化法を施して q_{1}, ⋯, q_{n} を作る。つまり

v_{k} := a_{k} − ∑_{i=1}^{k−1}⟨a_{k}, q_{i}⟩q_{i}, q_{k} =^{def.} \frac{v_{k}}{∥v_{k}∥} (k = 1, 2, ⋯, n)

として計算するわけだが^{4}、中間変数 r_{ik} (1 ≤ i ≤ k ≤ n) を導入して、次のように書き換えておこう。

(6)  v_{k} := a_{k} − ∑_{i=1}^{k−1} r_{ik}q_{i}, r_{ik} := ⟨a_{k}, q_{i}⟩,
(7)  q_{k} := \frac{v_{k}}{r_{kk}} , r_{kk} := ∥v_{k}∥ (k = 1, 2, ⋯, n).

さて

Q := (q_{1} q_{2} ⋯ q_{n}), R := ( \begin{matrix}r_{11} & r_{12} & ⋯ & r_{1n} \\ & r_{22} & ⋯ & r_{2n} \\ & & ... & ⋮ \\ 0 & & & r_{nn} \end{matrix} )

とおくと、 Q は unitary 行列で、

A = QR.

ここまで A は複素行列として計算してきたが、実行列である場合は、 Q, R も実行列 (したがって Q は実直交行列) である。

4これは線型代数の教科書の、 Gram-Schmidt の直交化法の説明に載っている式であるが、数値計算に用いるには問題も多いとのこと。

[figure]

定義 5.1 (QR 分解) 正則行列 A ∈ GL(n; C) に対して、

A = QR (Q は unitary 行列, R は対角線分がすべて正の上三角行列)

となる Q, R を見い出すことを A を QR 分解すると言う。

[figure]

上の議論から任意の正則行列は QR 分解可能であることが分かったが、実はこれは一意的である。

[figure]

{命題 5.2 (QR 分解の一意性)任意の A ∈ GL(n; C) に対して、 QR 分解はただ一つしか \atop ない。}

[figure]

証明 A = QR は

a_{k} = ∑_{i=1}^{k} r_{ik}q_{i} (k = 1, 2, ⋯, n)

と書ける。 k = 1 についての条件

a_{1} = r_{11}q_{1}

から

r_{11} = |r_{11}| = \frac{∥a_{1}∥}{∥q_{1}∥} = ∥a_{1}∥,

q_{1} = \frac{a_{1}}{r_{11}} .

次に k = 2 についての条件

a_{2} = r_{12}q_{1} + r_{22}q_{2}

から、まず q_{1} との内積を取って、

r_{12} = ⟨a_{1}, q_{1}⟩.

これから

v_{2} := a_{2} − r_{12}q_{1}

は計算できて、

v_{2} = r_{22}q_{2}

であるから、

r_{22} = |r_{22}| = \frac{∥v_{2}∥}{∥q_{2}∥} = ∥v_{2}∥, q_{2} = \frac{v_{2}}{r_{22}} .

以下 k = 3, ⋯, n と順に同様に計算できる。別証明 A の二つの QR 分解

A = Q_{1}R_{1}, A = Q_{2}R_{2}

があれば、 Q_{1}R_{1} = Q_{2}R_{2} から

(8)  Q_{2}^{∗}Q_{1} = R_{2}R_{1}^{−1}.

この等式の左辺は unitary 行列である。実際

(Q_{2}^{∗}Q_{1})(Q_{2}^{∗}Q_{1})^{∗} = Q_{2}^{∗}Q_{1}Q_{1}^{∗}Q_{2}^{∗∗} = Q_{2}^{∗}Q_{2} = I.

また (8) の右辺は上三角行列で、対角成分はすべて正である。

unitary かつ上三角かつ対角成分がすべて正という行列は単位行列に限られることは容易に証明できる。ゆえに

Q_{2}^{∗}Q_{1} = R_{2}R_{1}^{−1} = I.

これから

Q_{1} = Q_{2}, R_{1} = R_{2}.

注意 5.3 (QR 分解を利用した連立 1 次方程式の解法)A = QR という QR 分解があるとき、

Ax = b

は

QRx = b

であるから

x = R^{−1}(Q^{∗}b).

上三角行列の逆行列をかける計算は、乗算 n^{2}/2 回程度の計算量で計算できるから、 x は乗算 3n^{2}/2 回程度の計算量で計算できることが分かる。 QR 分解は O(n^{3}) 程度の計算量で求まることは分かるから、 Gauss の消去法に基づく LU 分解とオーダーだけは匹敵する計算法である。 (決して有利な方法ではないが。 )

5.2 QR 変換

A = QR

という QR 分解が得られたとき、順番を変えて掛け算した行列 A_{1} を作る。

A_{1} = RQ

これは R = Q^{∗}A より

A_{1} = Q^{∗}AQ

となるので、実は A を Q で相似変換したものである。これを QR 変換と言う。

固有値を求めるために QR 変換を用いる場合、まず Hessenberg 行列に相似変換した後に使われることが多い。

• Hessenberg 行列の QR 変換は Hessenberg 行列である。

• 実対称行列の QR 変換は実対称行列である。 (∵ 対称行列の直交行列による相似変換は実対称行列であるから。 )

• Hermite 行列の QR 変換は Hermite 行列である。

• 三重対角行列の QR 変換は三重対角行列とは限らないが、実対称三重対角行列の QR 変換は実対称三重対角行列である (実対称三重対角 == 実対称かつ Hessenberg だから)。

5.3 QR 法の原理

行列 A が与えられたとき、

A_{0} = A = Q_{0}R_{0}, A_{1} = R_{0}Q_{0},
A_{1} = Q_{1}R_{1}, A_{2} = R_{1}Q_{1},
(9)  A_{2} = Q_{2}R_{2}, A_{3} = R_{2}Q_{2},
⋯
A_{k} = Q_{k}R_{k}, A_{k+1} = R_{k}Q_{k},
⋯

と QR 変換を繰り返すと、 A_{m} の対角線の下側の成分はすべて 0 に収束する。

注意 5.4 (たまにある誤解を正す)このことを「A_{k} は上三角行列に収束する」と言うことがあるが、対角線の上にある成分が特定の値に収束するわけではない。言い換えると、特定の上三角行列 U があって、

lim_{k→∞} A_{k} = U

となるわけではない。

一般性の追求はほどほどにして、次の定理を証明しよう。

[figure]

定理 5.5 (QR 法の原理) A ∈ GL(n; C) で、その固有値 λ_{j} (j = 1, 2, ⋯, n) はすべて相異なり、

|λ_{1}| > |λ_{2}| > ⋯ > |λ_{n}|

と仮定する。このとき k → ∞ とすると、 A_{k} の対角線より下のすべての成分は 0 に収束し、 A_{k} の対角成分は λ_{1}, ⋯, λ_{n} に収束する。

[figure]

注意 5.6 実対称行列については、上の定理の仮定はまあまあ一般的であるが、非対称な行列については、虚根を持たないことが必要になるので (虚根を持てば、複素共役となるので、絶対値は等しくなってしまう)、あまり一般的な仮定とは言えない。

証明

ステップ 1 (A^{k+1} の QR 分解) まず

主張 1

[figure]

\begin{matrix}P_{k} := Q_{0}Q_{1} ⋯ Q_{k−1}Q_{k} とおくと \\ (10) & A_{k} = P_{k−1}^{∗}AP_{k−1}. \end{matrix}

[figure]

実際、

A_{1} = R_{0}Q_{0} = (Q_{0}^{∗}A)Q_{0} = Q_{0}^{∗}AQ_{0},
A_{2} = R_{1}Q_{1} = (Q_{1}^{∗}A_{1})Q_{1} = Q_{1}^{∗}(Q_{0}^{∗}AQ_{0})Q_{1} = Q_{1}^{∗}Q_{0}^{∗}AQ_{0}Q_{1},
A_{2} = R_{2}Q_{2} = (Q_{2}^{∗}A_{2})Q_{2} = Q_{2}^{∗}(Q_{1}^{∗}Q_{0}^{∗}AQ_{0}Q_{1})Q_{2} = Q_{2}^{∗}Q_{1}^{∗}Q_{0}^{∗}AQ_{0}Q_{1}Q_{2},
⋯
A_{k} = R_{k}Q_{k} = (Q_{k}^{∗}A_{k−1})Q_{k} = Q_{k}^{∗}(Q_{k−1}^{∗} ⋯ Q_{1}^{∗}Q_{0}^{∗}AQ_{0}Q_{1} ⋯ Q_{k−1})Q_{k},
⋯

となる (帰納法で証明すべきかもしれないが)。

主張 2

[figure]

U_{k} := R_{k}R_{k−1} ⋯ R_{1}R_{0} とおくと A^{k+1} = P_{k}U_{k} であり、これは A^{k+1} の QR 分解である。 [figure]

まず

P_{k}U_{k} = (Q_{0}Q_{1} ⋯ Q_{k−1}Q_{k})(R_{k}R_{k−1} ⋯ R_{1}R_{0}) = (Q_{0}Q_{1} ⋯ Q_{k−1})(Q_{k}R_{k})(R_{k−1} ⋯ R_{1}R_{0})
= P_{k−1}A_{k}U_{k−1}

で、主張 1 より P_{k−1}A_{k} = AP_{k−1} だから、

P_{k}U_{k} = AP_{k−1}U_{k−1}.

この式を再帰的に用いて、

P_{k}U_{k} = AP_{k−1}U_{k−1} = A(AP_{k−1}U_{k−2}) ⋯ = A^{k}P_{0}U_{0} = A^{k}Q_{0}R_{0} = A^{k}A = A^{k+1}.

一方、 P_{k} は unitary 行列の積であるから unitary 行列、 U_{k} は対角線分が正である上三角行列の積であるから対角成分が正である上三角行列である。ゆえに A^{k+1} = P_{k}U_{k} は A^{k+1} の QR 分解である。

ステップ 2 (A_{k} の挙動) A は固有値がすべて相異なるから対角化が可能である。すなわち ∃M ∈ GL(n; C) s.t.

M^{−1}AM = Λ, Λ = diag(λ_{1}, ⋯, λ_{n}).

M を M = QR と QR 分解、 M^{−1} を M^{−1} = LU (ただし L の対角成分はすべて 1) と LU 分解して、

(11)  A^{k+1} = MΛ^{k+1}M^{−1} = (QR)Λ^{k+1}(LU) = QR(Λ^{k+1}LΛ^{−(k+1)})Λ^{k+1}U.

ここで Λ^{k+1}LΛ^{−(k+1)} は下三角行列で、対角成分はすべて 1 に等しく、 i > j のとき (i, j) 成分は

ℓ_{ij} (\frac{λ_{i}}{λ_{j}} )^{k+1}

に等しい (ただし ℓ_{ij} は L の (i, j) 成分)。そこで

Λ^{k+1}LΛ^{−(k+1)} = I + E_{k}

とおくと、 lim_{k→∞} E_{k} = O. 上式を (11) 代入して、

(12)  A^{k+1} = QR(Λ^{k+1}LΛ^{−(k+1)})Λ^{l+1}U
= QR(I + E_{k})(R^{−1}R)(Λ^{l+1}U) = QR(I + E_{k})R^{−1}(RΛ^{l+1}U)
= Q(I + RE_{k}R^{−1})(RΛ^{l+1}U)

この I + RE_{k}R^{−1} を QR 分解する:

(13)  I + RE_{k}R^{−1} = S_{k}T_{k}.

I + RE_{k}R^{−1} → I (k → ∞) で、 QR 分解の連続性は明らかだから、

lim_{k→∞} S_{k} = I, lim_{k→∞} T_{k} = I.

(13) を (12) に代入して、

(14)  A^{k+1} = Q(S_{k}T_{k})(RΛ^{k+1}U) = (QS_{k})(T_{k}RΛ^{k+1}U).

この右辺は unitary 行列と、上三角行列の積の形になっている。

ステップ 3 (QR 分解との関係) U = (u_{ij}) として、

D_{1} = diag (\frac{\overline{λ_{1}}}{|λ_{1}|}, \frac{\overline{λ_{2}}}{|λ_{2}|}, ⋯, \frac{\overline{λ_{n}}}{|λ_{n}|}) , D_{2} = diag (\frac{\overline{u_{11}}}{|u_{11}|}, \frac{\overline{u_{22}}}{|u_{22}|}, ⋯, \frac{\overline{u_{nn}}}{|u_{nn}|})

とおく。これを使って (14) を

A^{k+1} = (QA_{k}D_{2}^{∗}D_{1}^{k∗})[(D_{1}^{k}D_{2})(T_{k}R)(D_{1}^{k}D_{2})^{∗}(D_{1}Λ)^{k}(D_{2}U)]

のように書き換えると、これは A^{k+1} の QR 分解になる。すでに A^{k+1} = P_{k}U_{k} という QR 分解が得られていたので、 QR 分解の一意性から、

(15)  P_{k} = QS_{k}D_{2}^{∗}D_{1}^{k∗}.

ステップ 4 (10) と (15) から

A_{k+1} = P_{k}^{∗}AP_{k} = D_{1}^{k}D_{2}S_{k}^{∗}(Q^{∗}AQ)S_{k}D_{2}^{∗}D_{1}^{k∗}.
A = MΛM^{−1} = (QR)Λ(R^{−1}Q^{∗})

より

Q^{∗}AQ = RΛR^{−1}.

この行列は (右辺に注目することで) 上三角行列で、かつ対角線上に λ_{1}, ⋯, λ_{n} がこの順にならんでいる。 k → ∞ のとき、 S_{k} → I だから

A_{k+1} → D_{1}^{k}(D_{2}RΛR^{−1}D_{2}^{∗})D_{1}^{∗k}

この右辺は上三角行列で、対角成分は λ_{1}, λ_{2}, ⋯, λ_{n} である。この定理から k → ∞ とするときの収束の速さは

(\frac{λ_{i}}{λ_{j}} )^{k}

という因子に関わることが分かる。

5.4 QR 分解のアルゴリズム

QR 分解の存在証明として、線型代数の本でおなじみの Gram-Schimidt の直交化法を用いたが、 QR 分解を実際に計算するには、以下のような工夫をする。

5.4.1 修正 Gram-Schmidt の方法

5.4.2 Givens 変換を用いる方法

5.4.3 鏡映変換を用いる方法

5.5 QR 法の歴史

1961 年に Francis が複素行列の固有値の計算法として導入したもの (一松 [2])。

5.6 Octave での実験

MATLAB や Octave には qr() という QR 分解をする関数があって、 QR 法の実験が簡単にできる。

以下、実験をする際の参考に、小プログラム集を。

# qr_t.m --- QR 変換
function QR = qr_t(a)
  [q r]=qr(a);
  QR=r*q;
endfunction
# qr_iteration.m --- QR 変換の反復
function QR = qr_iteration(a,n)
  QR = a;
  for i=1:n
    [q r]=qr(QR);
    QR=r*q;
  end
endfunction

# norm_l.m --- 下三角部分 (対角線含まず) のノルム

function n = norm_l(a)
  n = norm(a-triu(a));
endfunction
# rand_h.m --- テスト用 Hessenberg 行列
function ret = rand_h(n)
  ret = triu(rand(n,n)) + diag(rand(n-1,1),-1);
endfunction
# rand_t.m --- テスト用三重対角行列
function ret = rand_t(n)
  ret=diag(rand(n,1),0)+diag(rand(n-1,1),1)+diag(rand(n-1,1),-1);
endfunction
# rand_s.m --- テスト用対称行列
function ret = rand_s(n)
  a=rand(n,n);
  ret = (a+a')/2;
endfunction
# rand_st.m --- テスト用対称三重対角行列
function ret = rand_st(n)
  u=rand(n-1,1);
  ret=diag(rand(n,1),0)+diag(u,1)+diag(u,-1);
endfunction
# rand_hermite.m --- テスト用 Hermite 行列
function ret = rand_hermite(n)
  a=rand(n,n)+i*rand(n,n);
  ret = (a+a')/2;
endfunction

6 特異値

いつか勉強して書きたい。

日本語の解説はかつては、一松 [2], 一松 [3] くらいしか見当たらなかったが、今では色々ある。

7 固有値計算のためのパッケージ

(準備中)

A 参考書

森・杉原・室田 [5], 一松 [2], [3] 以外には、名取 [9], 戸川 [6], シャトラン [1] など。

B Strum の方法

Strum の方法とは、実軸上の区間 [a, b] における多項式 f(x) の実根の個数に関する有名な Strum の定理をもとにした多項式の根の計算法であり、数値計算法のあちこちで顔を出す。ここで簡単にまとめておく。なお、以下の記述は森 [4], 一松 [2] から適当に抜き書きしたものであるが、実は高木 [8] の第 3 章の方が詳しいかもしれない (暇があったら書き直そう)^{5}。

B.1 スツルムの定理

[figure]

定義 B.1 (Strum 列) [a, b] を R の区間とするとき、次の 4 つの条件を満足する実係数多項式列

f_{0}(x) = f(x), f_{1}(x), f_{2}(x), ⋯, f_{ℓ}(x)

は区間 [a, b] において Strum 列をなすという。

(1) ∀x ∈ [a, b] に対して、隣り合う二つの多項式 f_{k}(x), f_{k+1}(x) は同時に 0 にはならない。

(2) x_{0} ∈ [a, b], k ∈ {1, 2, ⋯, ℓ − 1} において f_{k}(x_{0}) = 0 となれば、 f_{k−1}(x_{0})f_{k+1}(x_{0}) < 0.

(3) 列の最後の多項式 f_{ℓ}(x) は [a, b] において一定の符号を持つ。

(4) x_{0} ∈ [a, b] において f(x_{0}) = 0 ならば f′(x_{0})f_{1}(x_{0}) > 0. [figure]

5何となく、すごく高尚なことが書いてあるのではないかと、中を見もしないで敬遠していたのだが…損をした気分でいます。

[figure]

定義 B.2 (符号変化の回数)実係数多項式の列 f(x) = f_{0}(x), f_{1}(x), ⋯, f_{ℓ}(x) が区間 [a, b] で Strum 列をなすとき、 f(x) の根でない x ∈ [a, b] に対して、関数値の列

f_{0}(x), f_{1}(x), ⋯, f_{ℓ}(x)

を左から右に見ていったときの符号の変化の回数 N(x) が定義できる。

• x において、すべての k につき f_{k}(x) ≠ 0 が成り立つ場合の N(x) の定義は明らか。

• ある k について f_{k}(x) = 0 となった場合は、仮定 (Strum 列の条件 (3) と x が f(x) の根でないという仮定) から k = 0, ℓ はありえず 1 ≤ k ≤ ℓ − 1 で、その場合も Strum 列の条件 (2) によって、 f_{k−1}(x)f_{k+1}(x) < 0 となることから、向う三軒両隣では

\begin{matrix}f_{k−1}(x) & f_{k}(x) & f_{k+1}(x) \\ + & 0 & − \end{matrix} または\begin{matrix}f_{k−1}(x) & f_{k}(x) & f_{k+1}(x) \\ − & 0 & + \end{matrix}

のいずれかのパターンしかありえない。そこで、この 3 項で 1 回符号変化したと数えることにする。

[figure]

例 B.3 +, +, −, 0, +, +, − では 3 回と数える。

[figure]

定理 B.4 (Strum) 実係数多項式の列 f(x) = f_{0}(x), f_{1}(x), ⋯, f_{ℓ}(x) は区間 [a, b] で Strum 列をなし、 f(a)f(b) ≠ 0 であるとする。このとき、 f(x) ≠ 0 なる x ∈ [a, b] に対して、関数値の列

f_{0}(x), f_{1}(x), ⋯, f_{ℓ}(x)

を左から右に見ていったときの符号の変化の回数を N(x) とすると、方程式 f_{0}(x) = 0 の区間 [a, b] における解の個数は N(a) − N(b) である。

[figure]

この定理によって、区間 [a, b] 内に存在する f(x) = 0 の解を数を知ることができるが、二分探索の技法と組み合わせることで、解が存在する区間を好きなだけ小さくすることが出来、その区間内の適当な点を解の近似値として採用するという近似解法ができる。これを Strum の方法と呼ぶ。証明 f(x) = 0 の解の個数は有限個である。そのうち [a, b] に含まれるものを大きさの順に並べて

x_{1} < x_{2} < ⋯ < x_{n}

とする。 n + 1 個の区間

I_{0} = [a, x_{1}), I_{1} = (x_{1}, x_{2}), I_{2} = (x_{2}, x_{3}), I_{n−1} = (x_{n−1}, x_{n}), I_{n} = (x_{n}, b]

の和集合において (f(x) は 0 とならないので) N(x) が定義できる。以下では (1) 各区間 I_{j} において N(x) は定数である:

∃{n_{j}}_{j=0}^{n} s.t. N(x) = n_{j} (x ∈ I_{j}, j = 0, 1, ⋯, n)

(2) n_{j+1} − n_{j} = −1 (j = 0, 1, ⋯, n − 1)

であることを証明する (この二つから容易に N(a) − N(b) = n が導かれる)。

(1) の証明 I_{j} 内のある部分区間 J において、すべての多項式 f_{k}(x) ≠ 0 とすると、 f_{k}(x) の符号は J で一定であるので (中間値の定理から、符号が変化するには 0 にならないといけない)、 N(x) の値は変化しないことが分かる。

そこで、ある多項式 f_{k}(x) の解 x_{∗} ∈ I_{j} を通過した場合を考える。 Strum 列の条件から k ≠ ℓ, I_{j} の取り方から k ≠ 0 であるから、 1 ≤ k ≤ ℓ − 1 である。 Strum 列の条件 (2) から f_{k−1}(x_{∗})f_{k+1}(x_{∗}) < 0 であるが、連続性から x_{∗} の十分小さな近傍 V を取れば、そこで f_{k−1}(x) と f_{k+1}(x) は定符号となり、 V 上で f_{k−1}(x)f_{k+1}(x) < 0 がなりたつ。それゆえ、列

f_{0}(x), f_{1}(x), ⋯, f_{k−1}(x), f_{k}(x), f_{k+1}(x), ⋯, f_{ℓ}(x)

の部分列

f_{k−1}(x), f_{k}(x), f_{k+1}(x)

の符号については次の二つのいずれか一方だけしか起こらない。

(a) \begin{matrix}f_{k−1}(x) & f_{k}(x) & f_{k+1}(x) \\ V 上 − & & V 上 + \end{matrix} (b) \begin{matrix}f_{k−1}(x) & f_{k}(x) & f_{k+1}(x) \\ V 上 + & & V 上 − \end{matrix}

それゆえ、 (f_{k}(x) の符号が何であっても) 任意の x ∈ V について、部分列 f_{k−1}(x), f_{k}(x), f_{k+1}(x) における符号の変化は 1 として N(x) の値に算入される。これから x_{∗} の十分小さな近傍で N(x) の値に変化はないことが分かる。ゆえに N(x) は I_{j} 上で定数である。

(2) の証明 f(x) = 0 の解 x_{j} において考える。 Strum 列の条件 (4) は f′(x_{j}) と f_{1}(x_{j}) が同符号であることを示している。例えば f′(x_{j}) > 0 の場合、十分小さな ε > 0 を取ると、

• f(x) < 0 (x ∈ [x_{j} − ε, x_{j})), f(x_{j}) = 0, f(x) > 0 (x ∈ (x_{j}, x_{j} + ε]).

• f_{1}(x) > 0 (x ∈ (x_{j} − ε, x_{j} + ε)).

これから

n_{j+1} − n_{j} = N(x_{j} + 0) − N(x_{j} − 0) = −1

であることが分かる。 f′(x_{j}) < 0 の場合も同様の議論で n_{j+1} − n_{j} = −1 であることが分かる。

B.2 ユークリッドの互除法による Strum 列の生成

多項式 f(x) ∈ R[x] が与えられたとき、 f_{0}(x) = f(x) と f_{1}(x) = f′(x) から Euclid の互除法を行い、関数列 f_{0}(x), f_{1}(x), ⋯, f_{ℓ}(x) を作る:

(16)  f_{0}(x) := f(x), f_{1}(x) := f′(x),
f(x) = q_{1}(x)f_{1}(x) − f_{2}(x), deg f_{2}(x) < deg f_{1}(x),
f_{1}(x) = q_{2}(x)f_{2}(x) − f_{3}(x), deg f_{3}(x) < deg f_{2}(x),
f_{2}(x) = q_{3}(x)f_{3}(x) − f_{4}(x), deg f_{4}(x) < deg f_{3}(x),
⋮
(17)  f_{k−1}(x) = q_{k−1}(x)f_{k}(x) − f_{k+1}(x), deg f_{k+1}(x) < deg f_{k}(x),
⋮
f_{ℓ−2}(x) = q_{ℓ−1}(x)f_{ℓ−1}(x) − f_{ℓ}(x), deg f_{ℓ}(x) < deg f_{ℓ−1}(x),
f_{ℓ−1}(x) = q_{ℓ}(x)f_{ℓ}(x).

(普通の互除法と異なり、 f_{k−1}(x) を f_{k}(x) で割ったときの普通の剰余の (−1) 倍を f_{k+1}(x) とすることに注意しよう。そうする理由は以下の (18) を成立させるためである。 )

よく知られているように f_{ℓ}(x) は f(x) と f′(x) の最大公約多項式であるから、 f(x) ∈ R[x] が重根を持たない場合、 f_{ℓ}(x) ≡ 定数(≠ 0) となることに注意しよう。以下この場合に f(x) = 0 の解 (根) を求めることを考える。 f(x) が重根を持つ場合は f(x) の代わりに g(x) = f(x)/f_{ℓ}(x) を考えることで同様の議論ができる。

[figure]

定理 B.5 f(x) ∈ R[x] が重根を持たないとき、 f(a)f(b) ≠ 0 を満たす任意の区間 [a, b] において、 f(x) と f′(x) から互除法で作った剰余列 (詳しくは (16), (17) 参照)

f(x) = f_{0}(x), f′(x) = f_{1}(x), f_{2}(x), ⋯, f_{ℓ}(x)

は Strum 列をなす。

[figure]

証明まず f_{ℓ}(x) ≡ 定数であるから、 Strum 列の条件 (3) は満たされている。次にある x_{0} ∈ [a, b], ある k ∈ {0, 1, ⋯, ℓ − 1} に対して

f_{k}(x_{0}) = f_{k+1}(x_{0}) = 0

となったとすると、式 (17) から f_{k+2}(x_{0}) 以降の f_{j}(x_{0}) もすべて 0 になり、特に f_{ℓ}(x_{0}) = 0. これは f_{ℓ}(x) が定数関数 (≠ 0) であることに矛盾する。ゆえに Strum 列の条件 (1) が満たされる。次にある点 x_{0}, ある k ∈ {1, 2, ⋯, ℓ − 1} に対して f_{k}(x_{0}) = 0 となったとすると、式 (17) から

(18)  f_{k−1}(x_{0}) = −f_{k+1}(x_{0}).

ゆえに Strum 列の条件 (2) も満たされる (条件 (3) から上式の値は 0 にならないことに注意)。最後に x_{0} が f(x) の根であるとき、 f(x) が重根を持たないという仮定から f′(x_{0}) ≠ 0 で、

f′(x_{0})f_{1}(x_{0}) = f′(x_{0})^{2} > 0

となり、条件 (4) も満たされる。

B.3 3 重対角行列の固有多項式と Strum 列

T = ( \begin{matrix}a_{1} & b_{1} & 0 & & 0 \\ b_{1} & a_{2} & b_{2} \\ & ... & ... & ... \\ & & b_{N−2} & a_{N−1} & b_{N−1} \\ 0 & & 0 & b_{N−1} & a_{N} \end{matrix} )

を実対称三重対角行列とする。

注意 B.6 (対称性の仮定について)実は以下の議論で T の対称性はあまり本質的でなく、対角線をはさんで対称の位置にある成分が同符号、すなわち T の第 (i, j) 成分を t_{ij} と書くとき、

t_{k,k−1}t_{k−1,k} > 0 k = 2, ⋯, N

が成り立つことが本質的であるという指摘が一松先生の本にあった。なるほど。行列を Hessenberg 形にする算法を施すとき、特に与えられた行列が実対称であった場合には、結果が実対称三重対角になってしまうということで、 (以下 Strum の方法を用いて固有値を求める段階になると) 実対称性はある意味では必要性はあまりないのだが、最初に実対称性がないと三重対角化できないわけで、まあ必ずついてくる「おまけ」ということでしょう。

以下 b_{k} ≠ 0 (k = 1, 2, ⋯, N − 1) と仮定する。もしある k に対して b_{k} = 0 ならば

T = ( \begin{matrix}T′ & O \\ O & T′′ \end{matrix} )

とブロック分けでき、 T′, T′′ の固有値を求める問題に帰着できるから、一般性は失われない。 p_{k}(λ) を λI − T の k 次首座小行列式とする (k = 0, 1, ⋯, N)。すなわち

(19)  p_{k}(λ) := \{ \begin{matrix}det(λI_{k} − T_{k}) & (k = 1, 2, ⋯, N) \\ 1 & (k = 0). \end{matrix}

ただし、

I_{k} = k 次の単位行列, T_{k} = ( \begin{matrix}a_{1} & b_{1} & 0 & & 0 \\ b_{2} & a_{2} & b_{1} \\ & ... & ... & ... \\ & & b_{k−2} & a_{k−1} & b_{k−1} \\ 0 & & 0 & b_{k−1} & a_{k} \end{matrix} ) .

すぐ分かる命題を二つ。

[figure]

\begin{matrix}補題 B.7 (漸化式) (19) で定義された \{p_{j}(λ)\}_{j=0}^{N} について、以下が成り立つ。 \\ \{ & \begin{aligned}p_{0}(λ) & = & 1, \\ p_{1}(λ) & = & λ − a_{1}, \\ p_{k+1}(λ) & = & (λ − a_{k+1})p_{k}(λ) − b_{k}^{2}p_{k−1}(λ) (k = 1, 2, ⋯, N − 1), \\ p_{N}(λ) & = & det(λI − T). \end{aligned} \end{matrix}

[figure]

証明行列式の展開定理を用いる。

[figure]

補題 B.8 (Strum 列であること)(19) で定義された {p_{j}(λ)}_{j=0}^{N} を逆順に並べた多項式列

p_{N}(λ), p_{N−1}(λ), ⋯, p_{1}(λ), p_{0}(λ)

は任意の閉区間 [a, b] において Strum 列である。すなわち

(1) 隣り合う二つの多項式 p_{k}(λ), p_{k+1}(λ) は共通根を持たない。

(2) ある λ_{0} ∈ R において p_{k}(λ_{0}) = 0 ならば p_{k−1}(λ_{0})p_{k+1}(λ_{0}) < 0 (k = 1, 2, ⋯, N − 1).

(3) 列の最後の多項式 p_{0}(λ) は R において定符号である。

(4) ある λ_{0} ∈ R において p_{N}(λ_{0}) = 0 ならば p′_{N}(λ_{0})p_{N−1}(λ_{0}) > 0.

[figure]

証明

(1) もしも p_{k}(λ_{0}) = p_{k+1}(λ_{0}) = 0 とすると、漸化式から b_{k}^{2}p_{k−1}(λ_{0}) = 0. b_{k} = 0 と仮定した

から p_{k−1}(λ_{0}) = 0. これを繰り返すと

0 = p_{k+1}(λ_{0}) = p_{k}(λ_{0}) = p_{k−1}(λ_{0}) = ⋯ = p_{2}(λ_{0}) = p_{1}(λ_{0}) = p_{0}(λ_{0}).

これから

p_{0}(λ_{0}) = 0

これは p_{0}(λ) ≡ 1 に矛盾する。

(2) p_{k}(λ_{0}) = 0 を漸化式に代入すると p_{k+1}(λ_{0}) = −b_{k}^{2}p_{k−1}(λ_{0}). 前項より左辺 ≠ 0. これから p_{k+1}(λ_{0}), p_{k−1}(λ_{0}) は異符号である。

(3) p_{0}(λ) ≡ 1 であるから明らか。

(4) 漸化式

(20)  p_{k}(λ) = (λ − a_{k})p_{k−1}(λ) − b_{k−1}^{2}p_{k−2}(λ)

を微分すると、

(21)  p′_{k}(λ) = p_{k−1}(λ) + (λ − a_{k})p′_{k−1}(λ) − b_{k−1}^{2}p′_{k−2}(λ).

(21) × p_{k−1}(λ) − (20) × p′_{k−1}(λ) より

p′_{k}(λ)p_{k−1}(λ) − p_{k}(λ)p′_{k−1}(λ) = b_{k−1}^{2} (p′_{k−1}(λ)p_{k−2}(λ) − p_{k−1}(λ)p′_{k−2}(λ)) + p_{k−1}(λ)^{2}

という漸化式が得られる。ここで

q_{k}(λ) := p′_{k}(λ)p′_{k−1}(λ) − p_{k}(λ)p_{k−1}(λ)

とおくと、漸化式は

q_{k}(λ) = p_{k−1}(λ)^{2} + b_{k−1}^{2}q_{k−1}(λ) (k = 2, 3, ⋯, N).

となる。ところで

q_{1}(λ) = p′_{1}(λ)p_{0}(λ) − p_{1}(λ)p′_{0}(λ) = p′_{1}(λ) = 1 · 1 − (λ − α_{1}) · 0 = 1 > 0

であるから、以下帰納的に

q_{k}(λ) > 0 (k = 2, 3, ⋯, N)

が示せる。特に

q_{N}(λ) = p′_{N}(λ)p_{N−1}(λ) − p_{N}(λ)p′_{N−1}(λ) > 0

であるが、 p_{N}(λ) = 0 であるから

p′_{N}(λ)p_{N−1}(λ) > 0.

符号の変化数の計算に関する注意 p_{N}(a) ≠ 0 なる a に対して

N(a) := “{p_{0}(a), p_{1}(a), ⋯, p_{N}(a)}” の符号の変化数

とおく (逆順であっても符号の変化数は同じになる)。例えば

\begin{matrix}p_{0}(a) & p_{1}(a) & p_{2}(a) & p_{3}(a) & p_{4}(a) & p_{5}(a) & p_{6}(a) & p_{7}(a) & p_{8}(a) & p_{9}(a) & p_{10}(a) \\ + & − & − & − & 0 & + & + & + & + & − & − \end{matrix}

では N(a) = 3.

この符号の変化数は (特別な注意をせずに) 数値的に安定して計算できる。つまり絶対値が非常に小さくて、符号の判別がつきにくい場合も、「符号の変化数」そのものは疑いがなく計算できる。例えば

\begin{matrix}p_{k−1}(a) & p_{k}(a) & p_{k+1}(a) \\ + & 絶対値小 & − \end{matrix} または\begin{matrix}p_{k−1}(a) & p_{k}(a) & p_{k+1}(a) \\ − & 絶対値小 & + \end{matrix}

において p_{k}(a) が正であっても負であっても 0 であっても符号の変化数の計算にとっては影響がない。注意すべきは Strum 列の条件 (iv) から

\begin{matrix}p_{k−1}(a) & p_{k}(a) & p_{k+1}(a) \\ + & 絶対値小 & + \end{matrix} または\begin{matrix}p_{k−1}(a) & p_{k}(a) & p_{k+1}(a) \\ − & 絶対値小 & − \end{matrix}

のような場合 (もしこうなったら符号の変化数の計算がむつかしい) が起こり得ないことである。

B.4 直交多項式の作る Strum 列

(ここは単なる覚え書き。後で肉付けするかもしれない。 ) w: [a, b] → R は連続で、有限個の点で 0 になる他は正で、条件

sup_{k∈N} ∫_{a}^{b} x^{k}w(x) dx < ∞

を満たすような関数とする。このとき [a, b] 上の実数値連続関数全体の集合に

(u, v)_{w} := ∫_{a}^{b} u(x)v(x)w(x) dx

で定義される内積を導入して、内積空間としたものを H_{w}(a, b) とする。また (·, ·)_{w} に付随するノルムを ∥ · ∥_{w} と書く:

∥u∥_{w} = \sqrt{(u, u)_{w}}.

自然数 n に対して、関数列 {1, x, ⋯, x^{n}} から Gram-Schimidt の直交化法によって得られる直交多項式系を {p_{0}(x), p_{1}(x), ⋯, p_{n}(x)} とする。

[figure]

命題 B.9 H_{w}(a, b) において関数列 {1, x, ⋯, x^{n}} から Gram-Schimidt の直交化法によって得られる直交多項式系を {p_{0}(x), p_{1}(x), ⋯, p_{n}(x)} とし、 p_{n}(x) の最高次の係数を μ_{n} とすると、

inf {∥q∥_{w}; q(x) ∈ R[x], deg q(x) = n, q(x) の最高次係数 = 1} = ∥p_{n}/μ_{n}∥_{w}. [figure]

[figure]

命題 B.10 H_{w}(a, b) において関数列 {1, x, ⋯, x^{n}} から Gram-Schimidt の直交化法によって得られる直交多項式系を {p_{0}(x), p_{1}(x), ⋯, p_{n}(x)} とすると、 p_{n}(x) の根はすべて単根で、区間 [a, b] の内部にある。

[figure]
[figure]

命題 B.11 H_{w}(a, b) において関数列 {1, x, ⋯, x^{n}} から Gram-Schimidt の直交化法によって得られる直交多項式系を {p_{0}(x), p_{1}(x), ⋯, p_{n}(x)} とすると、次のような 3 項漸化式が存在する。

p_{k}(x) = (α_{k}x + β_{k})p_{k−1}(x) − γ_{k}p_{k−2}(x) (k = 1, 2, ⋯).

ただし p_{−1}(x) ≡ 0 とし、

α_{k} = \frac{μ_{k}}{μ_{k−1}} , β_{k} = −\frac{α_{k}(xp_{k−1}, p_{k−1})_{w}}{λ_{k−1}} ,
γ_{k} = \frac{α_{k}(xp_{k−1}, p_{k−2})_{w}}{λ_{k−2}} = \frac{μ_{k}μ_{k−2}λ_{k−1}}{μ_{k−1}^{2}λ_{k−2}} .

ただし、 λ_{k} = (p_{k}, p_{k})_{w}, μ_{k} = p_{k}(x) の最高次係数。

[figure]

[figure]

命題 B.12 H_{w}(a, b) において関数列 {1, x, ⋯, x^{n}} から Gram-Schimidt の直交化法によって得られる直交多項式系を {p_{0}(x), p_{1}(x), ⋯, p_{n}(x)} とすると (最高次の係数が正であるようにしておくと)、

p_{n}(x), p_{n−1}(x), ⋯, p_{1}(x), p_{0}(x)

は区間 [a, b] において Strum 列をなす。

[figure]

B.5 一般化された Strum 列

一松先生の本にあった Strum 列の定義の条件は、森先生の本よりも少し緩い条件であった (つまり一般化してあると言える)。明示していないけれど p_{k}(λ) の連続性は仮定されている、と思う。

次の性質を持つ関数列 {p_{k}(λ)}_{k=0}^{n} を Strum 系と言う。

(1) 各 p_{k}(λ) の解は有限個; p_{k}(λ_{0}) = p_{k+1}(λ_{0}) = 0 はありえない。

(2) p_{k}(λ_{0}) = 0 のとき p_{k−1}(λ_{0})p_{k+1}(λ_{0}) < 0.

(3) p_{0}(λ) は定符号。

もちろん対応する Strum の定理の結論は少し複雑になる。

[figure]

定理 B.13 (一般化された Strum の定理){p_{k}(λ)}_{k=0}^{n} を Strum 系とする。 p_{n}(a)p_{n}(b) ≠ 0 (a < b) のとき、

N(a) − N(b) = ∑_{p_{n}(λ_{0})=0} χ(λ_{0}).

ただし χ(λ_{0}) は、 λ を λ_{0}− から λ_{0}+ に変化させるときの p_{n}(λ)/p_{n−1}(λ) の符号変化を表す量で

χ(λ_{0}) := \{ \begin{matrix}1 & (− から + に変化) \\ 0 & (符号の変化なし) \\ −1 & (+ から − に変化). \end{matrix}
[figure]

証明略

[figure]

系 B.14 N(a) ≠ N(B) ならば [a, b] 内に p_{n}(λ) = 0 の解がある。

[figure]

[figure]

系 B.15 (実対称三重対角行列の固有値問題)n 次実対称三重対角行列

T = ( \begin{matrix}a_{1} & b_{1} & 0 & & 0 \\ b_{1} & a_{2} & b_{2} \\ & ... & ... & ... \\ & & b_{n−2} & a_{n−1} & b_{n−1} \\ 0 & & 0 & b_{n−1} & a_{n} \end{matrix} ) , b_{k} ≠ 0 (k = 1, 2, ⋯, n − 1)

から p_{k}(λ) = det(λI_{k} − T_{k}) (T_{k} は T の第 k 次首座行列) として作った Strum 列 p_{n}(λ), p_{n−1}(λ), ⋯, p_{1}(λ), p_{0}(λ) ≡ 1 について、 n 個の固有値はすべて実単根であり、各固有値 λ_{0} において、つねに χ(λ_{0}) = 1 である。ゆえに p_{n}(a)p_{n}(b) ≠ 0 となる a < b に対して

N(a) − N(b) = [a, b] 内の p_{N}(λ) = 0 の根の個数= [a, b] 内の T の固有値の個数.

そして p_{k}(λ) の隣接した根の中間に、必ず p_{k−1}(λ) = 0 の根がある (k = 1, 2, ⋯, n)。

[figure]

証明略。

参考文献

[1] F. Chatelin. Valeurs propres de matrices. Masson, Paris, 1988. (邦訳) F. シャトラン著, 伊理正夫, 伊理由美訳, 行列の固有値問題, シュプリンガー・フェアラーク東京 (1993).

[2] 一松^{ひとつまつ}信^{しん}. 数値解析. 朝倉書店, 1982/10.

[3] 一松信. 代数学入門第一〜三課. 近代科学社, 1992, 1992, 1994.

[4] 森正武^{まさたけ}. 数値解析第 2 版. 共立出版, 2002/2/25. 第 1 版は 1973 年に出版された。

[5] 森正武, 杉原正顯^{まさあき}, 室田^{むろた}一雄^{かずお}. 線形計算. 岩波講座応用数学. 岩波書店, 1994.

[6] 戸川隼人^{はやと}. マトリクスの数値計算. オーム社, 1971.

[7] 福澤誠人. 応用数理実験課題レポート実対称行列の三重対角化. http://nalab.mind.

   meiji.ac.jp/~mk/labo/text/members/20030904fukuzawa-report.pdf, 2003.

[8] 高木貞治. 代数学講義改訂新版. 共立出版, 1965.

[9] 名取亮. 線形計算. 朝倉書店, 1993.