固有値問題ノートの補足

桂田祐史

2005 年 8 月 21 日, 2008 年 6 月 9 日, 2013 年 8 月 26 日

(大規模作業中) この文書は今のところ

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

においてある。将来的には

『行列の固有値問題』

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

にマージする予定である。なお、一般化固有値問題については別に

『一般化固有値問題』

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

を用意してある。

目次

1 いろは 2

2 概略 5

3 基本的な実直交行列とそれによる相似変換 5 3.1 Householder 変換 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 3.2 Givens 変換 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

4 QR 分解 8 4.1 QR 分解の定義 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 4.2 修正 Gram-Schmidt の直交化法 . . . . . . . . . . . . . . . . . . . . . . . . . . 10 4.3 Givens 変換による QR 分解 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 4.4 Householder 変換による QR 分解 . . . . . . . . . . . . . . . . . . . . . . . . . 13 4.5 各方法の比較 (1) 素朴にやってみよう… . . . . . . . . . . . . . . . . . . . . . 15 4.6 各方法の比較 (2) お勉強してみよう . . . . . . . . . . . . . . . . . . . . . . . . 17 5 Jacobi 法 17

5.1 Rutishauser の計算式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 5.2 MATLAB での実験 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

6 Lanczos アルゴリズム 24 6.1 はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 6.2 クリロフ部分空間 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 6.3 Lanczos 原理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 6.4 都合が「よい」状況下での Lanczos 法 . . . . . . . . . . . . . . . . . . . . . . 27 6.5 都合が「よくない」場合の Lanczos 法 . . . . . . . . . . . . . . . . . . . . . . 28 6.6 そのほか . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

7 摂動定理と分離定理 30

8 min-max principle 31

9 一般化固有値問題 31 9.1 正則な行列束に関する一般化固有値問題の定義 . . . . . . . . . . . . . . . . . . 31 9.2 Courant-Fischer の min-max principle . . . . . . . . . . . . . . . . . . . . . . . 32 9.3 2 分法の一般化固有値問題への拡張 . . . . . . . . . . . . . . . . . . . . . . . . 32

10 Shur 分解 32

11 Sylvester の慣性律, 2 次形式の対角化 33

11.1 本題に入る前に . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

11.1.1 行列についての色々な同値関係 . . . . . . . . . . . . . . . . . . . . . . 35 11.2 見掛けが少し異なる 3 つの Sylvester の慣性律 . . . . . . . . . . . . . . . . . . 37 11.3 Sylvester の慣性律の証明 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 11.4 2 次形式の対角化のアルゴリズム (後でどこかで使うもの) . . . . . . . . . . . 41

11.4.1 例に基づく説明 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 11.4.2 プログラム例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

A TODO 56

B ChageLog 56

1 いろは

A ∈ C^{n×n} に対して、

Ax = λx, x ≠ 0

を満たす x が存在するとき、 λ を A の固有値 (eigenvalue) と呼び、 x を λ に属する固有ベクトル (eigenvector) と呼ぶ。 λ が A の固有値であることと、

ker(λI − A) = {x ∈ C^{n}; Ax = λx}

が {0} でないことは同値である。

λ が A の固有値であるとき、 ker(λI − A) を固有空間 (eigenspace) と呼ぶ。

[figure]

定理 1.1 (固有値は固有多項式の根)λ が A の固有値であることと、 n 次代数方程式

det(λI − A) = 0

の根であることは同値である。

[figure]
[figure]

定義 1.2 (固有多項式、固有方程式)p(λ) := det(λI − A) を A の固有多項式、

p(λ) = 0

を固有方程式と呼ぶ。

[figure]
[figure]

系 1.3 () A ∈ C^{n×n} の固有値の個数は n 個以下である。

[figure]

A の固有多項式に重根がないときはちょうど n 個 (相異なる) 固有値が存在することになるが、重根がある場合はその重複度分だけ同じものをならべて、 λ_{1}, ⋯, λ_{n} と書くことが多い。以下、 n 次正方行列 A について、 A の固有値は λ_{1}, ⋯, λ_{n} である、と言ったときはそういう意味である。例えば、「対角行列 D = diag(d_{1}, ⋯, d_{n}) の固有値は対角成分 d_{1}, ⋯, d_{n} である。」…このことの証明は D の固有多項式についての等式

det(λI − D) = ∏_{i=1}^{n}(λ − d_{i})

による。ちなみにこの事実は三角行列に対して一般化できる。すなわち「A が上三角行列または下三角行列とするとき、 A の固有値は対角成分 a_{11}, ⋯, a_{nn} である」。

[figure]

{系 1.4 (固有値は固有多項式の根ということから) \atop (1)} 任意の非特異行列 P に対して、 A を相似変換した行列 P^{−1}AP と A の固有値は重複度も込めて一致する。

(2) A と A^{T} の固有値は重複度も込めて一致する。

(3) A の固有値の複素共役は重複度も込めて A^{∗} の固有値に一致する。

(4) det A は A の固有値 (重根の場合、重複度の分だけ並べる) の積に等しい。特に、

A が非特異 ⇔ 0 は A の固有値でない.

[figure]
[figure]

補題 1.5 A の相異なる固有値 λ_{1}, ⋯, λ_{m} に属する固有ベクトルは 1 次独立である。

[figure]

証明ヴァンデルモンドの行列式を使う。

[figure]

命題 1.6 A が相異なる n 個の固有値を持つとき、すなわち固有方程式が重根を持たないとき、以下の (1), (2) が成り立つ。

(1) A の各固有値 λ_{1}, ⋯, λ_{n} に属する固有ベクトルは全空間の基底をなす。

(2) A の各固有値に属する固有空間は 1 次元である。

[figure]

この命題の状況のように、 A の固有値 λ_{1}, ⋯, λ_{n} に属する固有ベクトル p_{1}, ⋯, p_{n} が全空間の基底となるとき、 P = (p_{1}, ⋯, p_{n}) とおくと、これは非特異で

P^{−1}AP = diag(λ_{1}, ⋯, λ_{n}).

このように適当な非特異行列 P に対して、 P^{−1}AP が対角行列 diag(λ_{1}, ⋯, λ_{n}) になるとき、 A は対角化可能という。このとき、 λ_{i} は A の固有値で、 p_{i} は λ_{i} に属する固有ベクトル、 p_{1}, ⋯, p_{n} は全空間の基底となる。

[figure]

命題 1.7 A の固有値全体を λ_{1}, ⋯, λ_{n} とするとき、以下がなりたつ。いずれの場合も対応する固有値の固有ベクトルは等しい。

(1) c ∈ C とするとき、 cA の固有値は cλ_{1}, ⋯, cλ_{n} である。

(2) c ∈ C とするとき、 cI + A の固有値は c + λ_{1}, ⋯, c + λ_{n} である。

(3) A が可逆であるとき、 A^{−1} の固有値は λ_{1}^{−1}, ⋯, λ_{n}^{−1} である。

[figure]

λ を A の固有値とするとき、\overline{λ} は A^{∗} の固有値であるが、その固有ベクトル x は

A^{∗}x = \overline{λ}x

を満たすので、

x^{∗}A = λx^{∗}

を満たす。 yA = λy, y ≠ 0 を満たす y のことを固有値 λ に属する行固有ベクトルとよぶ。

[figure]

定理 1.8 (Frobenius の定理) A の固有値を λ_{1}, ⋯, λ_{n} とするとき、任意の多項式 p(z) について、 p(A) の固有値は p(λ_{1}), ⋯, p(λ_{n}) である。

[figure]
[figure]

定義 1.9 (固有値の代数的多重度、幾何学的多重度) A ∈ C^{n×n} の固有値 λ に対して、固有値 λ の代数的多重度とは、固有多項式 p(z) = det(zI − A) における根 λ の多重度のことと定義する。また、固有値 λ の幾何学的多重度とは、 λ に属する固有空間 ker(λI − A) の次元のことと定義する。

[figure]

2 概略

対称な問題 (行列が実対称行列や Hermite 行列である場合) と、そうでない問題との間には大きな違いがある。

古典的な Jacobi 法は別にして、次の手順を採る。

1. 与えられた行列を「簡単化」する。具体的には、対称な問題の場合は三重対角化し、非対称な問題の場合には Hessenberg 形にする。直交行列による相似変換が基本である。

2. 簡単化された問題を反復法で解く。冪乗法の系統 (逆反復法, シフト法を含む)、二分法、 QR 法など。

この中継地点の発見をしたのは、 Givens であるという。戸川 [11] (1971) から引用する。

それを最初に行ったのはギブンスであった。彼は ◯ 1 のために、のちにギブンスの方法と呼ばれるようになった有名な算法を考案し、また ◯ 2 のために、現在ではスツルムの方法とか、バイセクション法と呼ばれている算法を考案した。それらは現在でも有力な手法として使われているが、その後、上記 ◯ 1 , ◯ 2 それぞれに対し、新しい方法が次つぎに発表され、どちらもきわめて能率よく計算できるようになった。現在、 ◯ 1 の方法としてはハウスホルダー法が決定版といわれており、 ◯ 2 の方法としては、小さい固有値を数個求めたいとか、大きい固有値を数個求めたいという場合はスツルム法、すべての固有値を求めたい場合は QR 法が最もよいとされている。問題によっては ◯ 1 の目的に (改良された) ランチョス法が用いられることがある。

ちなみに Wilkinson がかの “Algebraic Eigenvalue Problem” を著したのは 1965 年である。

3 基本的な実直交行列とそれによる相似変換

超平面に対する対称移動である Householder 変換と、回転である Givens 変換が代表的かつ —きほんて—

3.1 Householder 変換

R^{n} の超平面に関する対称移動は明らかに内積を保つので直交変換である。これを Householder 変換、鏡映変換、または基本直交変換と呼ぶ。

超平面を W_{u} = {x ∈ R^{n}; (x, u) = 0} (u ∈ R^{n} \ {0}) とおくと、 W_{u} に関する対称移動を表わす行列は

H = H(u) = I_{n} − 2\frac{uu^{T}}{∥u∥^{2}} .

この形からも、図形的な意味からも、 H^{−1} = H^{T} = H (実直交行列かつ実対称行列) が成り立つことが分かる。参考 N 次元の任意の直交変換は、 N 個以下のベクトルに関する鏡映変換の積に書ける (Cartan)。

[figure]

補題 3.1 (∥x∥ = ∥y∥ ≠ 0 なる 2 点を相互に写す鏡映変換)∥x∥ = ∥y∥, x ≠ y なる x, y ∈ R^{n} に対して、

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

とおくと Ux = y かつ Uy = x. [figure]

証明明らか。

計算の工夫 (1) H と与えられたベクトルとの積

H = H_{u} を記憶するのではなく、 u と α = \frac{2}{∥u∥^{2}} を記憶しておく。

x ∈ R^{n} が与えられたときに w := Hx を計算するには、

w := Hx = (I_{n} − 2\frac{uu^{T}}{∥u∥^{2}} ) x = x − \frac{2}{∥u∥^{2}} uu^{T} x = x − αuu^{T} x = x − (αu^{T} x)u

であるから、

β := α(u^{T} x), w = x − βu

と計算すればよい。加減算 2n − 1 回、乗算 2n + 1 回で十分である。

計算の工夫 (2) H による相似変換 H^{T} AH

Ã:= H^{T} AH

を計算するには、

α = \frac{2}{∥u∥^{2}} , p = A^{T} u, q = Au, v = αu, β = p^{T} v, Ã = A − vp^{T} − (q − βu)v^{T} .

特に A が実対称行列の場合は、

α = \frac{2}{∥u∥^{2}} , p = α(Au), β = \frac{α}{2}(p, u), q = p − βu, Ã = A − uq^{T} − qu^{T} .

3.2 Givens 変換

p, q ∈ {1, ⋯, n}, p ≠ q, θ ∈ R とするとき、「x_{p}x_{q} 平面における −θ の回転」を表わす行列 G = G(p, q, θ) = (g_{ij}) を

g_{ij} = \{ \begin{matrix}cos θ & (i = j ∈ \{p, q\}) \\ sin θ & ((i, j) = (p, q)) \\ − sin θ & ((i, j) = (q, p)) \\ 1 & (i = j ∉ \{p, q\}) \\ 0 & (その他) \end{matrix}

で定義する。これは明らかに実直交行列である。

Givens 変換の行列

[figure]

  1  % Givens.m
  2  %  G(p,q, θ) ∈ M(n;R) (p ≠ q, c=cos θ, s=sin θ)
  3  function G = Givens(n,p,q,c,s)
  4    G = eye(n,n);
  5    G(p,p) =   c; G(p,q) = s;
  6    G(q,p) = - s; G(q,q) = c;

[figure]

計算の工夫 (1) G^{T} を与えられた行列 A にかける

B := G(p, q, θ)^{T} A

は p, q 行目を除き A と変わらない。変更される行については、

b_{pj} = a_{pj} cos θ − a_{qj} sin θ (j = 1, 2, ⋯, n),
b_{qj} = a_{pj} sin θ + a_{qj} cos θ (j = 1, 2, ⋯, n).

θ を

cos θ = \frac{a_{qr}}{\sqrt{a_{pr}^{2} + a_{qr}^{2}}} , sin θ = \frac{a_{pr}}{\sqrt{a_{pr}^{2} + a_{qr}^{2}}}

を満たすように取れば、 b_{pr} = 0 となる。

計算の工夫 (2) 与えられた行列 A を G で相似変換 (G^{T} AG)

B := G^{T} AG, G := G(p, q, θ)

は p, q 行目, p, q 列目を除き A と変わらない。

b_{ij} = a_{ij} (i ≠ p, q かつ j ≠ p, q),
b_{pj} = a_{pj} cos θ − a_{qj} sin θ (j ≠ p, q),
b_{qj} = a_{pj} sin θ + a_{qj} cos θ (j ≠ p, q),
b_{ip} = a_{ip} cos θ − a_{iq} sin θ = b_{pi} (i ≠ p, q),
b_{iq} = a_{ip} sin θ + a_{iq} cos θ = b_{qi} (i ≠ p, q),
b_{pp} = a_{pp} cos^{2} θ − 2a_{pq} cos θ sin θ + a_{qq} sin^{2} θ,
b_{qq} = a_{pp} sin^{2} θ + 2a_{pq} cos θ sin θ + a_{qq} cos^{2} θ = a_{pp} + a_{qq} − b_{pp},
b_{pq} = (a_{pp} − a_{qq}) cos θ sin θ + a_{pq}(cos^{2} θ − sin^{2} θ),
b_{qp} = b_{pq}.

4 QR 分解

4.1 QR 分解の定義

R^{n} のベクトル列 a_{1}, ⋯, a_{n} が 1 次独立であるとき、 Gram-Schmidt の直交化法、つまり k = 1, 2, ⋯, n の順に

b_{k} := a_{k} − ∑_{j=1}^{k−1}(a_{k}, q_{j}), q_{k} := \frac{b_{k}}{∥b_{k}∥}

と計算することで正規直交基底 q_{1}, ⋯, q_{n} が作れる。 A =: (a_{1}, ⋯, a_{n}), Q := (q_{1}, ⋯, q_{n}), R := (r_{ij}), r_{jk} := (a_{k}, q_{j}) とおくと、

(1)  A = QR, Q は実直交行列, R は対角成分が正である上三角行列

と書ける。この形の分解は一意的である。実際

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

とすると、左から Q_{2}^{−1} = Q_{2}^{T}, 右から R_{1}^{−1} (これは対角成分が正の上三角行列である) をかけて、

Q_{2}^{T} Q_{1} = R_{2}R_{1}^{−1}.

左辺は実直交行列の積として実直交行列 (長さを変えない変換だから) で、右辺は対角成分が正の上三角行列の積として対角成分が正の上三角行列である。簡単な計算で対角成分が正の上三角実直交行列は単位行列に他ならないことが分かる。

以上をまとめて次の命題を得る。

[figure]

命題 4.1 (QR 分解) n 次正則行列 A に対して、

A = QR, Q は実直交行列, R は対角成分が正の上三角行列

を満たす Q, R が一意的に存在し、それは A の列ベクトルに対する Gram-Schmidt の直交化法により得られる。

[figure]

(1) の形の分解を QR 分解と呼ぶ^{1}。

1これを Schmidt 分解、あるいは Gram-Schmidt の分解と呼ぶことを提案する本もある (シャトラン [1], 齋藤 [9])。以下に見るように (1) の形の分解を得るアルゴリズムには色々あるので、個人的には、有名なアルゴリズムの名前に冠される “Gram-Schmidt” は適切ではないように感じている。

GramSchmidt の直交化による QR 分解

[figure]

 1  % GramSchmidt.m --- GramSchmidt()
 2

3 % 使用例

 4  % n=3; a=rand(n,n);
 5  % [q r]=GramSchmidt(a)
 6  % a-q*r
 7
 8  % A=(a1 a2 ...,an) の列ベクトルから正規直交基底
 9  % q1, q2,..., qn を並べた Q を求める。
10  function [q, r] = GramSchmidt(A)
11    [n n] = size(A);
12    %
13    q=zeros(n,n);
14    r=zeros(n,n);
15    for k=1:n
16      b=A(:,k); % 第 k 列 a_k
17      for j=1:k-1
18        r(j,k)=q(:,j)'*A(:,k); % r_{jk} = a_k と q_j の内積
19        b=b-r(j,k)*q(:,j);     % b_k = a_k-(a_k,q_j)q_j
20      end
21      r(k,k)=norm(b);  % r_{kk} = ||b_k||
22      q(:,k)=b/r(k,k); % q_k = b / r_{kk}
23    end

[figure]

 GramSchmidt.m の実行結果

[figure]

 1  >> n=4;a=rand(n,n)
 2  a =
 3      0.9827    0.1146    0.5668    0.9616
 4      0.8066    0.6649    0.8230    0.0589
 5      0.7036    0.3654    0.6739    0.3603
 6      0.4850    0.1400    0.9994    0.5485
 7  >> [q r]=GramSchmidt(a)
 8  q =
 9      0.6415   -0.6491   -0.3412   -0.2251
10      0.5266    0.7307   -0.1256   -0.4160
11      0.4593    0.1624   -0.0214    0.8730
12      0.3166   -0.1355    0.9313   -0.1185
13  r =
14      1.5319    0.6358    1.4229    0.9870
15           0    0.4517    0.2074   -0.5971
16           0         0    0.6196    0.1676
17           0         0         0    0.0086
18  >> a-q*r
19  ans =
20     1.0e-15 *
21           0         0         0    0.1110
22           0         0         0   -0.0139
23           0         0         0         0
24      0.0555         0         0    0.1110
25  >>

[figure]

注意 4.2 上の説明では R として対角成分がすべて正であるものを採用したが、応用上はあまり意味がないので、巷のソフトウェア (例えば MATLAB) の QR 分解用の関数 (手続き) の仕様では、必ずしもそうはなっていない。

4.2 修正 Gram-Schmidt の直交化法

実は前小節にあげた Gram-Schmidt の直交化法は、浮動小数点演算システムで計算する場合、丸め誤差の観点からは良くない (と言われている^{2})。代りに次にあげる修正 Gram-Schmidt の直交化法が利用される。

修正 Gram-Schmidt の直交化法による QR 分解

[figure]

 1  % ModifiedGramSchmidt.m --- ModifiedGramSchmidt()
 2

3 % 使用例

 4  % n=3; a=rand(n,n);
 5  % [q r]=ModifiedGramSchmidt(a)
 6  % a-q*r
 7
 8  % A=(a1 a2 ...,an) の列ベクトルから正規直交基底
 9  % q1, q2,..., qn を並べた Q を求める。
10  function [q,r] = ModifiedGramSchmidt(A)
11    [n n] = size(A);
12    %
13    q=zeros(n,n); r=zeros(n,n);
14    for j=1:n
15      r(j,j)=norm(A(:,j)); % 第 k 列 a_k
16      q(:,j)=A(:,j)/r(j,j);
17      for k=j+1:n
18        r(j,k)=q(:,j)'*A(:,k); % r_{jk} = a_k と q_j の内積
19        A(:,k)=A(:,k)-r(j,k)*q(:,j);
20      end
21    end

[figure]

2øA.Björk, Solving linear least squares problems by Gram-Schmidt orthogonalization, BIT, Vol.7 (1967), pp.1–21.

 ModifiedGramSchmidt.m の実行結果

[figure]

 1  >> n=4;a=rand(n,n)
 2  a =
 3      0.6273    0.6552    0.5947    0.7764
 4      0.6991    0.8376    0.5657    0.4893
 5      0.3972    0.3716    0.7165    0.1859
 6      0.4136    0.4253    0.5113    0.7006
 7  >> [q r]=ModifiedGramSchmidt(a)
 8  q =
 9      0.5700   -0.2729   -0.7192   -0.2886
10      0.6352    0.7226    0.2609   -0.0788
11      0.3609   -0.5861    0.6431   -0.3356
12      0.3759   -0.2447    0.0323    0.8932
13  r =
14      1.1005    1.1995    1.1492    1.0839
15           0    0.1046   -0.2985   -0.1386
16           0         0    0.1972   -0.2886
17           0         0         0    0.3008
18  >> a-q*r
19  ans =
20     1.0e-15 *
21           0         0   -0.1110         0
22           0         0         0         0
23           0         0         0   -0.0555
24           0         0         0         0
25  >>

[figure]

( 「良くない」のはなぜか？本当にそうなのか、どこかでチェックしたい。 )

4.3 Givens 変換による QR 分解

(説明を書く暇がないのでとりあえず算法だけでも)

Givens 変換による QR 分解

[figure]

 1  % QR_Givens.m --- QR_Givens()

2 % 注意: このプログラムは原理の説明用であり、計算の効率は非常に低い

 3

4 % 使用例

 5  % n=3; a=rand(n,n); [q r]=QR_Givens(a)
 6  % norm(a-q*r)
 7  % norm(q*q')
 8
 9  % A=(a1 a2 ...,an) の列ベクトルから正規直交基底
10  % q1, q2,..., qn を並べた Q を求める。
11  function [Q,R] = QR_Givens(A)
12    [n n] = size(A);
13    R=A; Q=eye(n,n);
14    %
15    for q=1:n-1
16      for p=q+1:n
17        d=sqrt(R(p,q)*R(p,q)+R(q,q)*R(q,q));
18        c=R(q,q)/d;
19        s=R(p,q)/d;
20        G=Givens(n,p,q,c,s); % xp, xq 平面でθの回転 (c=cos θ, s=sin θ)
21        R=G'*R;
22        Q=Q*G;
23      end
24    end

[figure]

 QR Givens.m の実行結果

[figure]

 1  >> n=4;a=rand(n,n)
 2  a =
 3      0.6273    0.6552    0.5947    0.7764
 4      0.6991    0.8376    0.5657    0.4893
 5      0.3972    0.3716    0.7165    0.1859
 6      0.4136    0.4253    0.5113    0.7006
 7  >> [q r]=QR_Givens(a)
 8  q =
 9      0.5700   -0.2729   -0.7192   -0.2886
10      0.6352    0.7226    0.2609   -0.0788
11      0.3609   -0.5861    0.6431   -0.3356
12      0.3759   -0.2447    0.0323    0.8932
13  r =
14      1.1005    1.1995    1.1492    1.0839
15     -0.0000    0.1046   -0.2985   -0.1386
16      0.0000    0.0000    0.1972   -0.2886
17      0.0000   -0.0000    0.0000    0.3008
18  >> a-q*r
19  ans =
20     1.0e-15 *
21           0    0.1110    0.1110    0.1110
22           0    0.1110    0.1110    0.0555
23     -0.0555         0         0   -0.0278
24     -0.0555    0.0555    0.1110         0
25  >>

[figure]

4.4 Householder 変換による QR 分解

算法の基本的な考え方は、 Householder 変換による実対称行列の三重対角化と同じである。

A_{0} := A, (2) A_{k+1} := U_{k}A_{k}, U_{k} := I − 2u_{k}u_{k}^{T} (k = 0, 1, …, n − 2) (3)

で {A_{k}}_{k=0}^{n−1} を定めたとき (ただし u_{k} はこれから後で定める単位ベクトル)、 A_{k} は “k 列だけ上三角”、すなわち

A_{k} = ( \begin{matrix}\begin{aligned}α_{1} & & & & ∗ \\ & . & . & . \\ 0 & & & & α_{k} \end{aligned} & B_{k} \\ O & C_{k} \end{matrix} )

という形にすることを目指す。もしこれができれば

UA = R, U := U_{n−2} ⋯ U_{1}U_{0}, R := A_{n−1}

で、 U は実直交行列、 R は上三角行列となるので、

A = QR, Q := U^{T}

という QR 分解が得られる。 u_{k} を

u_{k} = ( \begin{matrix}0 \\ v & k \end{matrix} ) , 0 ∈ R^{k}, v_{k} ∈ R^{n−k}

の形に取ると、

U_{k} = ( \begin{matrix}I_{k} & O \\ O & V_{k} \end{matrix} ) , V_{k} := I_{n−k} − 2v_{k}v_{k}^{T}

という形なるので、

A_{k+1} = U_{k}A_{k} = ( \begin{matrix}\begin{aligned}α_{1} & & & & ∗ \\ & . & . & . \\ 0 & & & & α_{k} \end{aligned} & B_{k} \\ O & V_{k}C_{k} \end{matrix} ) .

我々の目標を達成するには

V_{k}C_{k} = ( α_{k+1} ; 0 ; ⋮ ; 0 ∗ )

という形をしていればよい。それには、

v_{k} := \frac{1}{∥x_{k} − y_{k}∥}(x_{k} − y_{k}), x_{k} := C_{k} の第 1 列, y_{k} := ( α_{k+1} ; 0 ; ⋮ ; 0 ) , α_{k+1} = ± ∥x_{k}∥

と取ればよい。ここで ± となっているのは、どちらを選んでもよいが、通常丸め誤差を抑えるために (桁落

ちがおきないように)、 x_{k} の第 1 成分の符号と逆になるように取るべし、という意見がある。次のプログラムではそうしてあるが、その結果は MATLAB の組み込み関数である qr() の結果と一致した。

Householder 変換による QR 分解のサンプル・プログラム

[figure]

1  % qrhouse.m --- Householder 法による QR 分解
2  % 試すには n=5;a=rand(n,n);a=(a+a')/2;qr_house(a)
3  function [q,r] = qrhouse(a)
4    [n,n]=size(a);
5    U=eye(n,n);
6    for k=0:n-2
7      x=a(k+1:n,k+1);
8      alpha=-sign(x(1))*norm(x);
9      x(1)=x(1)-alpha;
10      v=x/norm(x);
11      V=eye(n-k,n-k)-2*v*v';
12      a(k+1:n,k+1:n)=V*a(k+1:n,k+1:n);
13      U(k+1:n,:)=V*U(k+1:n,:);
14    end
15    q=U';
16    r=a;
17  %end function

[figure]

実行結果

[figure]

 1  >> n=5;a=rand(n,n);a=(a+a')/2
 2  a =
 3      0.6273    0.7683    0.5569    0.5571    0.3849
 4      0.7683    0.3716    0.4683    0.7887    0.6153
 5      0.5569    0.4683    0.7764    0.6480    0.2756
 6      0.5571    0.7887    0.6480    0.7036    0.3125
 7      0.3849    0.6153    0.2756    0.3125    0.5668
 8  >> [q r]=qrhouse(a)
 9  q =
10     -0.4739    0.2935    0.1766    0.5104   -0.6306
11     -0.5804   -0.6970    0.3665   -0.2008    0.0519
12     -0.4206   -0.1361   -0.7970    0.3050    0.2764
13     -0.4209    0.4579   -0.1876   -0.7491   -0.1295
14     -0.2908    0.4470    0.4051    0.2121    0.7117
15  r =
16     -1.3238   -1.2876   -1.2151   -1.3813   -0.9518
17      0.0000    0.5390    0.1514   -0.0125    0.0431
18      0.0000   -0.0000   -0.3587   -0.1344    0.2448
19     -0.0000    0.0000    0.0000   -0.1372    0.0431
20      0.0000    0.0000         0   -0.0000    0.2283
21  >> [q2 r2]=qr(a)
22  q2 =
23     -0.4739    0.2935    0.1766    0.5104   -0.6306
24     -0.5804   -0.6970    0.3665   -0.2008    0.0519
25     -0.4206   -0.1361   -0.7970    0.3050    0.2764
26     -0.4209    0.4579   -0.1876   -0.7491   -0.1295
27     -0.2908    0.4470    0.4051    0.2121    0.7117
28  r2 =
29     -1.3238   -1.2876   -1.2151   -1.3813   -0.9518
30           0    0.5390    0.1514   -0.0125    0.0431
31           0         0   -0.3587   -0.1344    0.2448
32           0         0         0   -0.1372    0.0431
33           0         0         0         0    0.2283
34  >> norm(q-q2)
35  ans =
36     3.2522e-15
37  >> norm(r-r2)
38  ans =
39     1.2993e-15
40  >>

[figure]

4.5 各方法の比較 (1) 素朴にやってみよう…

文献に書いてあることを理解して要約すべきだと思うが…MATLAB 時代の良いことは試すことは簡単にできるということだ。

 compare.m

[figure]

 1  function compare(n)
 2    a=rand(n,n);
 3  %% qr()
 4    fprintf('MATLAB の qr()\n');
 5    tic
 6    [q0 r0]=qr(a);
 7    toc
 8    fprintf('||Q R-A||=%e\n\n', norm(q0*r0-a));
 9  %%
10    fprintf(' 普通の Gram-Schmidt\n');
11    tic
12    [q1 r1]=GramSchmidt(a);
13    toc
14    fprintf('MATLAB の関数の結果との差: %e %e\n', norm(q1-q0), norm(r1-r0));
15    fprintf('||Q R-A||=%e\n\n', norm(q1*r1-a));
16  %%
17    fprintf(' 修正 Gram-Schmidt\n');
18    tic
19    [q2 r2]=ModifiedGramSchmidt(a);
20    toc
21    fprintf('MATLAB の関数の結果との差: %e %e\n', norm(q2-q0), norm(r2-r0));
22    fprintf('||Q R-A||=%e\n\n', norm(q2*r2-a));
23  %%
24    fprintf('Givens 変換による QR 分解\n');
25    tic
26    [q3 r3]=QR_Givens(a);
27    toc
28    fprintf('MATLAB の関数の結果との差: %e %e\n', norm(q3-q0), norm(r3-r0));
29    fprintf('||Q R-A||=%e\n\n', norm(q3*r3-a));
30  %%
31    fprintf('Householde 変換による QR 分解\n');
32    tic
33    [q4 r4]=qrhouse(a);
34    toc
35    fprintf('MATLAB の関数の結果との差: %e %e\n', norm(q4-q0), norm(r4-r0));
36    fprintf('||Q R-A||=%e\n\n', norm(q4*r4-a));

[figure]

ここでは手抜きをして、誤差と言うよりも残差を計算している。そのうち誤差を計算できるようにしてみたいが、分解の一意性がないのでこまったな…

n  qr()  Gram-Schmidt 修正 Gram-Schmidt  Givens  Householder
10 2.65 × 10^{−15} 3.40 × 10^{−16}  3.36 × 10^{−16}  1.50 × 10^{−15} 3.19 × 10^{−15}
20 5.46 × 10^{−15} 9.26 × 10^{−16}  8.84 × 10^{−16}  3.85 × 10^{−15} 4.13 × 10^{−15}
40 1.10 × 10^{−14} 1.69 × 10^{−15}  1.73 × 10^{−15}  1.32 × 10^{−14} 1.41 × 10^{−14}
80 1.65 × 10^{−14} 3.30 × 10^{−15}  3.20 × 10^{−15}  2.63 × 10^{−14} 1.77 × 10^{−14}
160 4.10 × 10^{−14} 6.75 × 10^{−15}  6.37 × 10^{−15}  8.32 × 10^{−14} 1.18 × 10^{−13}
320 1.66 × 10^{−13} 1.34 × 10^{−14}  1.34 × 10^{−14}  ?? × 10^{−14} ?? × 10^{−13}

少なくともこの実験結果だけからは、 Gram-Schmidt の直交化法、修正 Gram-Schmidt の直交化法による結果は他の方法による結果よりも精度が良さそうである。 MATLAB の qr() による結果はそれよりも悪い。 MATLAB の qr() による結果は Householder 法による結果に近い。

Gram-Schmidt の直交化法と修正 Gram-Schmidt の直交化法による結果の差は少ない。その二つの差よりも、 Givens 変換や Householder 変換を用いた方法との差の方がずっと大きい。普通の Gram-Schmidt の直交化法ではうまく解けないような条件の悪い問題を作って、それについて実験しないといけないのか？考え込まされる結果である。

4.6 各方法の比較 (2) お勉強してみよう

さて、自分でやってよく分からないということで、少し本をひもといてみることにしたら、最初にピンポンで、 Trefethen and Bau III [12] に色々載っていました。それからぐぐってみたら、この本に書いてあることは今では世界の常識みたい^{3}。この手のことを宿題でやらされている人がたくさんいるんだ…

今忙しいのでできないけれど... I shall return.

5 Jacobi 法

誰だったか忘れたが (戸川先生?)、 Jacobi 法をモグラ叩きに例えた人がいる。これはなかなかイメージが湧いてきて良いと思う。

[figure]

{補題 5.1 (実行列の要素の平方和)A = (a_{ij}) ∈ R^{n×n} について、 \atop ∑_{i,j=1}^{n} a_{ij}^{2} = tr(A^{T} A).}

[figure]

証明

tr(A^{T} A) = ∑_{i=1}^{n} A^{T} A の (i, i) 要素
= ∑_{i=1}^{n} ∑_{j=1}^{n} A^{T} の (i, j) 要素A^{T} の (j, i) 要素
= ∑_{i=1}^{n} ∑_{j=1}^{n} a_{ji}a_{ij}
= ∑_{i,j=1}^{n} a_{ji}^{2}.

実行列 A を実直交行列 U で相似変換した B = U^{T} AU について、

tr(B^{T} B) = tr(U^{T} A^{T} U · U^{T} AU) = (U^{T} A^{T} AU)

3Trefethen and Bau III の本の実験の mgs.m, clgs.m という名前のプログラムが世の中には一杯あるみたい。

となるが、良く知られた公式 tr(AB) = tr(BA) より^{4}、

tr(B^{T} B) = (A^{T} AUU^{T}) = tr(A^{T} A)

であるから、 B の要素の平方和= A の要素の平方和,

つまり

実直交行列による相似変換で要素の平方和は不変である。

[figure]

補題 5.2 (Jacobi のもぐらたたき一発分の効果)n ∈ N, n ≥ 2, A = (a_{ij}) は n 次実対称行列、 p, q ∈ {1, 2, …, n}, p ≠ q, θ ∈ R とするとき、 U := G(n, p, q, θ), B = (b_{ij}) := U^{T} AU とおくと、次の (1), (2) が成り立つ。

(1) (Givens 変換した行列の成分と平方和の性質)

b_{ij} = a_{ij} (i ≠ p, q かつ j ≠ p, q)
b_{pj} = a_{pj} cos θ − a_{qj} sin θ (j ≠ p, q)
b_{qj} = a_{pj} sin θ + a_{qj} cos θ (j ≠ p, q)
b_{ip} = a_{ip} cos θ − a_{iq} sin θ = b_{pi} (i ≠ p, q)
b_{iq} = a_{ip} sin θ + a_{iq} cos θ = b_{qi} (i ≠ p, q)
b_{pp} = a_{pp} cos^{2} θ − 2a_{pq} cos θ sin θ + a_{qq} sin^{2} θ
b_{qq} = a_{pp} sin^{2} θ + 2a_{pq} cos θ sin θ + a_{qq} cos^{2} θ = a_{pp} + a_{qq} − b_{pp}
b_{pq} = (a_{pp} − a_{qq}) cos θ sin θ + a_{pq}(cos^{2} θ − sin^{2} θ)
b_{qp} = b_{pq}

となるが、

∑_{i,j=1}^{n} a_{ij}^{2} = ∑_{i,j=1}^{n} b_{ij}^{2} (これはもっと一般に成り立つ) ,

∑_{i≠j} a_{ij}^{2} − ∑_{i≠j} b_{ij}^{2} = 2(a_{pq}^{2} − b_{pq}^{2}).

(2) 特に

θ = \frac{1}{2} Arctan (\frac{−2a_{pq}}{a_{pp} − a_{qq}} )

と取ると、 b_{pq} = b_{qp} = 0 かつ

∑_{i≠j} a_{ii}^{2} = ∑_{i≠j} b_{ii}^{2} + 2a_{pq}^{2}.

[figure]

4これから、相似変換で trace が不変であることが導かれる。今必要なのはそれではないけれど。

証明単純な計算なので略する。

[figure]

定理 5.3 (古典 Jacobi 法の原理)A を実対称行列とする。 A = A_{0} とおき、以下次の手順で A_{1}, A_{2}, ⋯ を計算する。 A_{k} = (a_{ij}^{(k)}) = (a_{ij}) の非対角要素のうちで、絶対値が最大の要素 a_{pq} に対して (絶対値最大の要素が二つ以上ある場合、任意の一つを選ぶ)、 () で回転行列 U_{k+1} を定め、

A_{k+1} = U_{k+1}^{T} A_{k}U_{k+1}

とおく。すると、

lim_{k→∞} A_{k} = D (対角行列)

となり、 D の対角成分が A の固有値、

V =^{def.} lim_{k→∞} V_{k}, V_{k} =^{def.} U_{1}U_{2} ⋯ U_{k}

の各列が固有ベクトルを与える。

[figure]

この定理に基づき、 A_{k} が十分対角行列に近づいたときに計算を打ち切り、 A の近似固有値として a_{11}^{(k)}, ⋯, a_{nn}^{(k)}, それらに属する近似固有ベクトルとして V_{k} の各列を採用する。これを古典 Jacobi 法と呼ぶ。

有限ステップで停止したときの誤差評価として、次の定理がある。

[figure]

定理 5.4 (近似固有値を非対角要素の平方和で評価)n 次実対称行列について、古典 Jacobi 法を適用して作った行列 A_{1}, A_{2}, ⋯ があるとき、 A_{k} = (a_{ij}^{(k)}) の非対角要素の平方和を F^{(k)} とすると、次の評価が成り立つ。

(1) A の固有値を適当に並べたものを λ_{1}, ⋯, λ_{n} とするとき、

|a_{ii}^{(k)} − λ_{i}| ≤ \sqrt{F^{(k)}} (i = 1, 2, ⋯, n).

(2) N = n(n − 1)/2 とおくとき、

F^{(k+1)} ≤ (1 − \frac{1}{N} ) F^{(k)} (k = 1, 2, ⋯).

(3) ある定数 C が存在して、十分大きな k に対して、

F^{(k+N)} ≤ C(F^{(k)})^{2}. [figure]

証明 (1) は摂動定理を使うだけ。 (2) は、まず (p, q) の選び方より

(a_{pq}^{(k)})^{2} ≥ \frac{F^{(k)}}{n(n − 1)}

であり、

F^{(k)} − F^{(k+1)} = 2(a_{pq}^{(k)})^{2} ≥ \frac{2}{n(n − 1)}F^{(k)} = \frac{1}{N} F^{(k)}.

(3) については省略。絶対値最大の要素の探索はかなり手間がかかる (計算量が大きくなる) ので、色々な変種が

考えられている。ある「閾値^{しきい}」 ε_{1} ≥ ε_{2} ≥ ⋯ と選び、 |a_{pq}^{(k)}| > ε_{k} なる a_{pq}^{(k)} に対して回転を施す閾値 Jacobi 法や、

[figure]

for p:=1 to n-1 do
for q:=p+1 to n

a_{pq} を叩く [figure]

を繰り返す特別巡回 Jacobi 法などがある。

注意 5.5 (Jacobi 法の評判) Jacobi 法の評判をいくつか引用しておく。戸川 [11] (1971) から、

最近の十数年間に数値計算の “常識” がずいぶん変わった。ある時期まで “最良の方法” と考えられていたものが後退し、それに代わって新しい解法が登場した。その中でも、固有値計算の分野における世代交代は特に著しいものがあった。たとえば、ヤコビの方法は対称行列の固有値計算法の決定版といわれ、どこの計算センターでも標準的なサブルーチンとして使っていたが、最近の文献では “ヤコビの方法は、公式が比較的簡単であるという小さな利点を除けば、もはや歴史的価値しかない” と極論されるところまで転落した。

森・杉原・室田 [6] (1994) から、

実対称行列の固有値 (と固有ベクトル) を求める古典的な方法に Jacobi 法がある。 20 年以上前には決定版のようにもてはやされたが、大規模行列は苦手という理由で、現在では (Rayleigh-Ritz の技法 (→ §6.5) と組み合わせて小規模な行列の固有値を求める場合を除いて) 実用的価値はあまりないと言われるようになってしまった。しかし、行列の次数が小さいなら十分実用であるし、固有値計算法の古典としての意義は大きい。

山本 [7] (2003) から、

古典 Jacobi 法と特別巡回 Jacobi 法は最終的には 2 次収束することが知られている。しかし、 30×30 程度の行列でも、数千回の回転を要するのが普通で、次数の大きい行列に対しては、計算時間および精度の点で、次節以下に述べる Householder– Givens 法, QR 法の法がすぐれている (表 7.1 参照)。しかし、固有値・固有ベクトルがすべて同時に求まる点で、この方法は捨て難い。特に、固有ベクトルの直交性はよく保たれる。実際、近年、国外で Jacobi 法の安定性が再評価されはじめているのは興味深い。

5.1 Rutishauser の計算式

θ を atan() を用いて求めたりする必要はない。

z = \frac{a_{qq} − a_{pp}}{2a_{pq}} (= cot 2θ),
t = \frac{sign (z)}{|z| + \sqrt{1 + z^{2}}} (= tan θ),
c = \frac{1}{\sqrt{1 + t^{2}}} (= cos θ),
s = ct (= sin θ),
u = \frac{s}{1 + c} (= tan θ/2),
b_{pp} = a_{pp} − ta_{pq},
b_{qq} = a_{qq} + ta_{pq},
b_{pq} = b_{qp} = 0,
b_{pj} = b_{jp} = a_{pj} − s(a_{qj} + ua_{pj}) (j ≠ p, q),
b_{qj} = b_{jq} = a_{qj} + s(a_{pj} − ua_{qj}) (j ≠ p, q).

5.2 MATLAB での実験

まず補題 5.2 (1) にあるように、与えられた (n 次) 実対称行列 A と p, q ∈ {1, 2, …, n}, p ≠ q に対して、 B := U^{T} AU, U := G(p, q, θ) を計算する関数は次のようになる (θ を与える代りに cos θ, sin θ を与えるようにしてある)。

jacobi0.m

[figure]

1  % jacobi0.m
2  % Jacobi 法の基礎ステップ --- x_p,x_q 平面の角θの回転 (Givens 変換)
3  %  c=cos θ, s=sin θ
4  function b=jacobi0(a,p,q,c,s)
5    b=a;
6    b(p,:)=a(p,:)*c-a(q,:)*s;
7    b(q,:)=a(p,:)*s+a(q,:)*c;
8    b(:,p)=b(p,:)';
9    b(:,q)=b(q,:)';
10    b(p,p)=a(p,p)*c*c-2*a(p,q)*s*c+a(q,q)*s*s;
11    b(q,q)=a(p,p)+a(q,q)-b(p,p);
12    b(p,q)=(a(p,p)-a(q,q))*c*s+a(p,q)*(c*c-s*s);
13    b(q,p)=b(p,q);

[figure]

b_{pq} = b_{qp} = 0 とするためには、補題 5.2 (2) にあるように θ を選ばねばならないが、これについては Rutishauser の計算式を一部利用すると、次のような関数を得る。

jacobi1.m

[figure]

1  % jacobi1.m
2  % Jacobi 法の 1 ステップ --- (p,q) 成分を叩く (0 になるよう Given 変換をする)
3  function b=jacobi1(a,p,q)
4    % cos θ,sin θを得る --- Rutishauser の計算式
5    z=(a(q,q)-a(p,p))/(2*a(p,q));   % cot 2 θ
6    t=sign(z)/(abs(z)+sqrt(1+z^2)); % tan θ
7    c=1/sqrt(1+t^2);                % cos θ
8    s=c*t;                          % sin θ
9    % u=s/(1+c); % tan(θ/2)
10    %
11    b=jacobi0(a,p,q,c,s);

12 % 丸め誤差により 0 になっていないのを強制的にクリア

13    b(p,q)=0;
14    b(q,p)=0;

[figure]

特別巡回 Jacobi 法は次のようになる。

jacobi.m

[figure]

1  % jacobi.m
2  % 巡回 Jacobi 法
3  function lambda=jacobi(a,iter)
4    [n,n]=size(a);
5    for k=1:iter
6      for p=1:n-1
7        for q=p+1:n
8          a=jacobi1(a,p,q);
9        end
10      end
11      %a
12    end
13    lambda=sort(diag(a));

[figure]

次の実験プログラムは、与えられた自然数 n に対して、 n 次実対称行列 A を作り、繰り返しの数を変えながら巡回 Jacobi 法による近似固有値を求めて、それを信頼できる値 (MATLAB の eig の結果) と比較している。

 jacobiexp.m

[figure]

 1  % jacobiexp.m
 2  %
 3  function jacobiexp(n)

4 % n 次実対称行列を生成

 5    a=rand(n,n);
 6    a=(a+a')/2;

7 % チェック用に MATLAB の関数で固有値を求めておく

 8    eigen=eig(a);

9 % タイマーをリセット

10    tic
11    for i=1:100

12 % jacobi() で i 回反復したときの近似固有値の精度を求める

13      error=norm(jacobi(a,i)-eigen);

14 % 結果を表示

15      fprintf('%d %e\n', i, error);

16 % 精度が十分ならば繰り返しをやめる

17      if (error < 1e-13)
18        break;
19      end
20    end

21 % 経過時間を表示

22    toc

[figure]

実行結果

[figure]

>> jacobiexp(5)
1 1.741441e-01
2 7.220988e-02
3 3.632924e-04
4 2.167315e-13
5 2.914633e-15
elapsed_time =
    0.1253
>> jacobiexp(10)
1 2.753837e-01
2 5.652731e-02
3 8.028047e-05
4 3.017277e-09
5 8.183436e-15
elapsed_time =
    0.5129
>> jacobiexp(20)
1 7.215286e-01
2 1.455863e-01
3 1.153771e-02
4 2.517012e-05
5 7.633347e-11
6 3.977614e-14
elapsed_time =
    6.6447
>>

[figure]

6 Lanczos アルゴリズム

6.1 はじめに

与えられた実対称行列を三重対角行列に相似変換するための Lanczos アルゴリズムを理解したい。

ゆくゆくは CG 法や直交関数系の三項漸化式の話まで込めて議論したい。

6.2 クリロフ部分空間

[figure]

定義 6.1 A ∈ R^{n×n}, x ∈ R^{n}, k ∈ N に対して

K_{k}(A, x) := Span{x, Ax, A^{2}x, …, A^{k−1}x}

とおき、 A の k 次 Krylov 部分空間とよぶ。

混乱が生じないときは (A, x を固定して考えることが多いので) K_{k}(A, x) を K_{k} と略記する。

[figure]

まず定義から容易に分かることをまとめておく。

[figure]

命題 6.2 A ∈ R^{n×n}, x ∈ R^{n} として、任意の k ∈ N に対して K_{k} := K_{k}(A, x) とおくとき、次の (1)–(4) が成り立つ。

(1) 任意の k ∈ N に対して AK_{k} ⊂ K_{k+1}.

(2) 任意の k ∈ N に対して K_{k} = {f(A)x; f(x) ∈ R[x], deg f(x) ≤ k − 1}.

(3) 任意の k ∈ N に対して K_{k} は R^{n} の線形部分空間である。また K_{k} は k につき単調

増加である:

K_{1} ⊂ K_{2} ⊂ ⋯ K_{j} ⊂ K_{j+1} ⊂ ⋯ ⊂ R^{n}.

(4) 任意の k ∈ N に対して

dim K_{k+1} = \{ \begin{matrix}dim K_{k} + 1 & (A^{k}x ∉ K_{k} のとき) \\ dim K_{k} & (A^{k}x ∈ K_{k} のとき). \end{matrix}

[figure]

証明

(1)

AK_{k} = A Span{x, Ax, A^{2}x, ⋯, A^{k−1}x} = Span{Ax, A^{2}x, A^{3}x, ⋯, A^{k}x}
⊂ Span{x, Ax, A^{2}x, A^{3}x, ⋯, A^{k}x} = K_{k+1}.

K_{k} = {c_{1}x + c_{2}Ax + c_{3}A^{2}x + ⋯ + c_{k}A^{k−1}x; (c_{j}) ∈ R^{k}}
= {(c_{1}I + c_{2}A + c_{3}A^{2} + ⋯ + c_{k}A^{k−1})x; (c_{j}) ∈ R^{k}}
= {f(A)x; f(x) ∈ R[x], deg f(x) ≤ k − 1}.

(3) Span であるから R^{n} の線形部分空間であることは明らかである。また K_{k} を生成する集合 {x, Ax, …, A^{k−1}x} 自体が k について単調増加であることから、 K_{k} が k について単調増加であることも明らかである。

(4) 前項より dim K_{k} ≤ dim K_{k+1}. K_{k+1} = Span{K_{k}, A^{k}x} であるから、 A^{k}x ∈ K_{k} であれば K_{k+1} = K_{k}. A^{k}x ∉ K_{k} であれば dim K_{k+1} = dim K_{k} + 1.

[figure]

{命題 6.3 A ∈ R^{n×n}, x ∈ R^{n}, x ≠ 0 とするとき、 ∃m ∈ \{1, 2, …, n\} s.t. \atop dim K_{j} = j (1 ≤ j ≤ m), K_{j} = K_{m} (j > m).}

[figure]

証明まず K_{1} = Span{x} であるから dim K_{1} = 1. 命題 6.2 から ∃m ∈ N s.t.

K_{1} ⊊ K_{2} ⊊ ⋯ ⊊ K_{m} = K_{m+1}.

A^{m}x ∈ K_{m+1} = K_{m} = Span{x, Ax, A^{2}x, …, A^{m−1}x}. これから容易に A^{j}x ∈ K_{m} (j ≥ m) が導かれる。ゆえに ∀j ≥ m に対して K_{j} = Span{x, Ax, …, A^{j−1}x} ⊂ K_{m}. もちろん K_{j} ⊃ K_{m} であるから K_{j} = K_{m}.

6.3 Lanczos 原理

以下、まずは m = n, すなわち

(4)  K_{n} = R^{n}

の場合を考える。

[figure]

定理 6.4 A ∈ R^{n×n}, x ∈ R^{n} に対して K_{n}(A, x) = R^{n} が成り立つとするとき、 K_{0} = {0} とおき、 k = 1, 2, …, n の順に v_{k} ∈ K_{k} \ K_{k−1} を任意に取り、

w_{k} := v_{k} − ∑_{j=1}^{k−1}(v_{k}, u_{j})u_{j},
u_{k} := \frac{w_{k}}{∥w_{k}∥}

とおくと (要するに Gram-Schmidt の直交化法)、 {u_{j}}_{j=1}^{n} は以下の (i), (ii), (iii) をみたす。

(i) (u_{j}, u_{k}) = δ_{jk} (1 ≤ j, k ≤ n). ここで δ_{jk} は Kronecker のデルタである。

(ii) 任意の k ∈ {1, …, n} に対して {u_{j}}_{j=1}^{k} は K_{k}(A, x) の基底である。

(iii) U := (u_{1} u_{2} ⋯ u_{n}) とおくと U^{T} AU =: B は Hessenberg 行列で^{a}、特に A が実対称行列である場合は三重対角行列になる。すなわち {α_{i}}_{i=1}^{n}, {β_{i}}_{i=1}^{n−1} が存在して、

U^{T} AU = ( \begin{matrix}α_{1} & β_{1} & & & & & 0 \\ β_{1} & α_{2} & β_{2} \\ & ... & ... & ... \\ & & β_{i−1} & α_{i} & β_{i} \\ & & & ... & ... & ... \\ & & & & β_{n−2} & α_{n−1} & β_{n−1} \\ 0 & & & & & β_{n−1} & α_{n} \end{matrix} ) .

aB = (β_{ij}) が j > i + 1 ⟹ β_{ij} = 0 をみたすとき、 B は Hessenberg 行列という。
[figure]

証明仮定 K_{n}(A, x) = R^{n} より、 dim K_{j}(A, x) = j (j = 1, 2, …, n) が成り立つので、 K_{j}\K_{j−1} は空でなく、 Gram-Schmidt の直交化法により R^{n} の正規直交基底が得られる。また作り方から明らかに任意の k に対して u_{k} ∈ K_{k}, dim K_{k} = k であるから、 {u_{j}}_{j=1}^{k} は K_{k} の基底になっている。

さて、任意の j に対して、 u_{j} ∈ K_{j} であるので、 Au_{j} ∈ K_{j+1} であるから、 Au_{j} は u_{1}, u_{2}, …, u_{min{j+1,n}} の線形結合になる。すなわち

Au_{j} = ∑_{k=1}^{min{j+1,n}} β_{jk}u_{k}

をみたす β_{jk} (k = 1, 2, …, min{j + 1, n}) が存在する。 k > j + 1 に対して β_{jk} := 0 とおいて、行列 B := (β_{jk}) を定めると、

AU = UB.

{u_{j}} が正規直交基底であることから U は実直交行列で U^{T} = U^{−1}. これを上の等式の両辺に左からかけて U^{T} AU = B. もしも A が実対称であれば、

B^{T} = (U^{T} AU)^{T} = U^{T} A^{T} (U^{T})^{T} = U^{T} AU = B.

すなわち B も対称である。対称な Hessenberg 行列は三重対角行列に他ならない。

注意 6.5 この定理の {u_{j}} は次の意味で「符号をのぞき一意的」である。すなわち v_{j} の取り方に任意性があるようだが、 K_{j} の中で K_{j−1} に直交する方向は 1 次元しかないので、 v_{j} をどう取り直しても結果として得られるものは u_{j}, −u_{j} のいずれか一方である。

6.4 都合が「よい」状況下での Lanczos 法

A を n 次の実対称行列で、 x ∈ R^{n} に対して、「都合よく」

K_{n}(A, x) = R^{n}

が成り立つと仮定する。前節の定理から、

U^{T} AU = ( \begin{matrix}α_{1} & β_{1} & & & & & 0 \\ β_{1} & α_{2} & β_{2} \\ & ... & ... & ... \\ & & β_{i−1} & α_{i} & β_{i} \\ & & & ... & ... & ... \\ & & & & β_{n−2} & α_{n−1} & β_{n−1} \\ 0 & & & & & β_{n−1} & α_{n} \end{matrix} )

をみたす実直交行列 U = (u_{1} u_{2} ⋯ u_{n}) が存在する。これから

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}, (5)

⋮
Au_{n−1} = β_{n−2}u_{n−2} + α_{n−1}u_{n−1} + β_{n−1}u_{n},
Au_{n} = β_{n−1}u_{n−1} + α_{n}u_{n}

という条件が導かれる。また

(6)  u_{1} = ±\frac{1}{∥x∥}x

である。実は {α_{i}}, {β_{i}}, {u_{j}} はこの条件 (5), (6) だけで定めることができる。最初の式と u_{1} との内積を取ると

α_{1} = u_{1}^{T} Au_{1}.

v_{2} := Au_{1} − α_{1}u_{1} とおくと、 β_{1}u_{2} = v_{2} であるから、

|β_{1}| = ∥v_{2}∥.

ゆえに

β_{1} = ±∥v_{2}∥, u_{2} := \frac{1}{β_{1}} v_{2} = ±\frac{1}{∥v_{2}∥}v_{2} (複号同順).

(5) の第 2 式と u_{2} との内積を取ると、

α_{2} = u_{2}^{T} Au_{2}.

v_{3} := Au_{2} − β_{1}u_{1} − α_{2}u_{2} とおくと、 β_{2}u_{3} = v_{3} であるから、

|β_{2}| = ∥v_{3}∥.

ゆえに

β_{2} = ±∥v_{3}∥, u_{3} := \frac{1}{β_{2}} v_{3} = ±\frac{1}{∥v_{3}∥}v_{3} (複号同順).

以下同様にして...

6.5 都合が「よくない」場合の Lanczos 法

A を n 次実対称行列とする。前節では

K_{n}(A, x) = R^{n}

を満たす x ∈ R^{n} が得られている場合に A の三重対角化を行ったが、この条件はそれほど安易に成り立つものではない。

任意の x に対して K_{n}(A, x) ⊊ R^{n} となる場合 A の最小多項式 φ(x) = x^{r} + ⋯ の次数 r が n より小さい場合は、

A^{r} = ∑_{j=1}^{r} c_{j}A^{j−1}

をみたす {c_{j}} が存在する。すると任意の x ∈ R^{n} に対して

A^{r}x = ∑_{j=1}^{r} c_{j}A^{j−1}x ∈ K_{r}(A, x).

したがって K_{r+1}(A, x) = K_{r}(A, x). ゆえに dim K_{n}(A, x) ≤ r < n となる。従って K_{n}(A, x) = R^{n} をみたす x は存在しない。

うまく x を取ると K_{n}(A, x) = R^{n} となる場合であっても一方 A の最小多項式の次数 r が n に等しい場合は、 K_{n}(A, x) = R^{n} をみたす x が存在することがわかるが (証明？ )、実際にどうやってそのような x を見つけるか、あまり簡単ではない。極端な話、 x を A の固有ベクトルに選んでしまった場合は dim K_{n}(A, x) = 1 であり、 Lanczos 法の反復は 1 ステップで止まってしまう。もっと極端な話、 A = λI の場合、いかなる x ≠ 0 も A の固有ベクトルになってしまう。おそらく最小多項式の次数が n である場合は、 generic に K_{n}(A, x) = R^{n} となるのであろうが… 一発でうまくやることはあきらめてここでは A の最小多項式の問題には触れないことにする。最初に選んだ x ≠ 0 に対して K_{n}(A, x) ≠ R^{n} と仮定する。 m := dim K_{n}(A, x) とおくと、 K_{m}(A, x) = K_{n}(A, x) である。 K_{m}(A, x) のことを単に K_{m} と書くことにする。

明らかに K_{m} は A 不変、すなわち

AK_{m} ⊂ K_{m}

をみたすが、 K_{m} の直交補空間

K_{m}^{⊥} := {y ∈ R^{n}; ∀z ∈ K_{m} z^{T} y = 0}

も A 不変である。実際 y ∈ K_{m}^{⊥} とするとき、任意の z ∈ K_{m} に対して、

(Ay, z) = (y, Az) = 0 (K_{m} の A 不変性より Az ∈ K_{m} なので)

がなりたつので Ay ∈ K_{m}^{⊥} となり、 K_{m}^{⊥} が A 不変であることが分かる。

老婆心: y ∈ K_{n}(A, x)^{⊥} とするとき、任意の j ∈ N に対して A^{j}y ∈ K_{n}(A, x)^{⊥} であるから、 K_{n}(A, y) ⊂ K_{n}(A, x)^{⊥}.

そこで最初に選んだ x_{1} ≠ 0 に対して、 m := dim K_{n}(A, x_{1}) < n であった場合は、まず K_{n}(A, x_{1}) の正規直交基底を取り、その後で x_{2} ∈ K_{n}(A, x_{1})^{⊥}, x_{2} ≠ 0 から始めて、 K_{n}(A, x_{2}) ⊂ K_{n}(A, x_{1})^{⊥} の正規直交基底を取り、以下この操作を続けていけば、最後には

R^{n} = K_{n}(A, x_{1}) ⊕ K_{n}(A, x_{2}) ⊕ ⋯ ⊕ K_{n}(A, x_{r})

とできる。それと並行して、 K_{n}(A, x_{j}) の正規直交基底 u_{k}^{j} (k = 1, 2, …, m_{j}) が得られる。この正規直交基底を

u_{1}^{1}, u_{2}^{1}, ⋯, u_{m_{1}}^{1}, u_{1}^{2}, u_{2}^{2}, ⋯, u_{m_{2}}^{2}, ⋯, u_{1}^{r}, u_{2}^{r}, ⋯, u_{m_{r}}^{r}

と並べたものを {u_{j}}_{j=1}^{n} として、 U := (u_{1} u_{2} ⋯ u_{n}) とおくと、

U^{T} AU = B = ( \begin{matrix}B_{1} \\ & B_{2} \\ & & . & . & . \\ & & & & & B_{r} \end{matrix} ) , B_{j} = ( \begin{matrix}α_{1}^{j} & β_{1}^{j} \\ β_{1}^{j} & α_{2}^{j} & β_{2}^{j} \\ & ... & ... & ... \\ & & β_{m_{j}−1}^{j} & α_{m_{j}}^{j} \end{matrix} ) .

こうして得られた B は確かに三重対角行列である。

余談 6.1 三重対角行列

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

の固有値の計算法の一つである二分法の適用には、条件 b_{j} ≠ 0 (j = 1, 2, …, n − 1) が必要であり、 b_{j} = 0 となる要素がある場合は行列をブロック分けして処理する、ということになっているが、 Lanczos 法では b_{j} が 0 になるかどうかは三重対角化の際にチェックできること、というかチェックしなければならないわけである。素朴な疑問をぶつぶつととはいえ、ここで説明した方法を実行するには、 {u_{k}} を記憶しておく必要がある。本当に Lanczos 法が必要とされる大規模な計算でそれが実行できるものなのだろうか？どうもうまく行かなければ u_{1} を取り直して再計算と考えている人達がいるような気がするが、それは dim K_{n} = R^{n} が成り立つような問題ばかりを相手にしているからか？

6.6 そのほか

7 摂動定理と分離定理

森・杉原・室田 [6] の演習問題から

[figure]

定理 7.1 () A = (a_{ij}), B = (b_{ij}) を実対称行列とし、 C = A − B とおく。それぞれ固有値の大きい方から r 番目のものを λ(·) とするとき、次が成立する。

(1) λ_{r}(B) + λ_{n}(C) ≤ λ_{r}(A) ≤ λ_{r}(B) + λ_{1}(C).

(2) |λ_{r}(A) − λ_{r}(B)| ≤ (∑_{i=1}^{n} ∑_{j=1}^{n}(a_{ij} − b_{ij})^{2})^{1/2}.

[figure]

山本 [7] から

[figure]

定理 7.2 (摂動定理) A = (a_{ij}), B = (b_{ij}) は n 次実対称行列で、それぞれ固有値 λ_{1} ≤ ⋯ ≤ λ_{n}, μ_{1} ≤ ⋯ ≤ μ_{n} を持つものとすれば、

|λ_{i} − μ_{i}| ≤ ∥A − B∥_{F} ≤ \sqrt{∑_{i,j=1}^{n}(a_{ij} − b_{ij})^{2}} (i = 1, 2, ⋯, n).

[figure]

この定理の一般化に Hoffman-Wielandt (1953)

∑_{i=1}^{n} |λ_{i} − μ_{i}|^{2} ≤ ∥A − B∥_{E}^{2} = ∑_{i=1}^{n} ν_{i}^{2} (ν_{1}, ⋯, ν_{n} は A − B の固有値).

[figure]

定理 7.3 (分離定理) A を n 次実対称行列, B を A の第 n 行, 第 n 列を除いて得られる n − 1 次行列とする (n > 1 と仮定する)^{a}。 A, B の固有値をそれぞれ λ_{1} ≤ ⋯ ≤ λ_{n}, μ_{1} ≤ ⋯ ≤ μ_{n} とすれば、

λ_{1} ≤ μ_{1} ≤ λ_{2} ≤ μ_{2} ≤ ⋯ ≤ λ_{n−1} ≤ μ_{n−1} ≤ λ_{n}.

aMATLAB 風に書けば、 B = A(1 : n − 1, 1 : n − 1).

[figure]

8 min-max principle

これは結構覚えにくい上に、本によって色々な形があって、混乱しがちである。まとめてみた (結局は伊理先生 [2] のまとめ方に近い)。

人名は本に載っていたものだが、根拠のあるものかどうかは分からない (出典がばらばらなので)。

[figure]

定理 8.1 (min-max principle) A を n 次実対称行列、その n 個の固有値 (重複度の分

だけ並べる) を λ_{1} ≥ λ_{2} ≥ ⋯ ≥ λ_{n} として、 R_{A}(x) := \frac{x^{T} Ax}{x^{T} x} とおくとき、次の (1)-(4) が成り立つ。

(1) (Fischer-Poincaré)

min_{dim V =k} max_{x∈V \{0}} R_{A}(x) = 小さい方から k 番目の固有値= λ_{n−k+1} (1 ≤ k ≤ n).

(2)

max_{dim V =k} min_{x∈V \{0}} R_{A}(x) = 大きい方から k 番目の固有値= λ_{k} (1 ≤ k ≤ n).

(3) (Courant-Fischer)

min_{dim S=j−1} max_{x∈S^{⊥}\{0}} R_{A}(x) = 大きい方から j 番目の固有値= λ_{j} (1 ≤ j ≤ n).

(4) (Courant-Weyl)

max_{dim S=j−1} min_{x∈S^{⊥}\{0}} R_{A}(x) = 小さい方から j 番目の固有値= λ_{n−j+1} (1 ≤ j ≤ n).

[figure]

9 一般化固有値問題

別に、桂田 [13] というノートを用意してある。

9.1 正則な行列束に関する一般化固有値問題の定義

A, B ∈ C^{n×m} とするとき、

{A − λB; λ ∈ C}

を行列束と呼ぶ。行列束 {A − λB} が正則であるとは、 A と B が同じ次数の正方行列であり、かつ

∃λ ∈ C s.t. det(A − λB) ≠ 0

が成り立つことと定義する。そうでないとき、行列束は特異であると定義する。正則な行列束 {A − λB} に対して、

Ax = λBx, x ≠ 0

を満たす x が存在するような λ の値と、そのときの x を求める問題を一般化された固有値問題と呼ぶ。

sp[A, B] =^{def.} {λ ∈ C; det(A − λB) = 0}.

とおき、 sp [A, B] の元を、行列束の固有値と呼ぶ。応用上は A, B が対称で、特に B が半正定符号であることが多い。

9.2 Courant-Fischer の min-max principle

[figure]

定理 9.1 A, B を n 次 Hermite 行列、 B は正定値とする。 (A, B) の固有値の大きいものから r 番目のものを λ_{r}(A, B) とするとき、

λ_{r}(A, B) = max_{dimS=r} min_{x∈S\\{0\}} \frac{x^{∗}Ax}{x^{∗}Bx} = max_{dimS=n−r+1} max_{u∈S\\{0\}} \frac{u^{∗}Au}{u^{∗}Bu}.

[figure]

9.3 2 分法の一般化固有値問題への拡張

A を実対称, B を正値対称行列とするとき、

Ax = λBx, x ≠ 0

を満たす λ, x を求めよ、という一般化固有値問題を考える。

(1) 固有値 λ はすべて実数である。

(2) (A, B) の固有値の正のものの個数 π̂(A, B) は π(A) に等しい。

(3) 任意に与えた実数 σ に対して、 A − σB = LDL^{T} と LDL^{T} 分解すれば、 D の対角要素のうち正, 0, 負のものの個数が、 σ より大きい、等しい、小さい固有値 λ の個数を与える。

10 Shur 分解

任意の A ∈ C^{n×n} に対して、 unitary 行列 U と上三角行列 S が存在して、

U^{H}AU = S.

ここで U^{H} は U の Hermite 共役 (共役転置) \overline{U}^{T} を表わす。これを A の Schur 分解という。

(事実自体は、線形代数の教科書、例えば佐武 [10] IV§3 にも書いてあるが、 Schur 分解という名前は冠されていないことが多いようである。 )

U = (u_{1}, u_{2}, ⋯, u_{n}) の列ベクトル u_{j} を Schur ベクトルと呼ぶ。

A が Hermite 行列ならば、固有ベクトルの直交性により、 S が対角行列になり、 u_{j} は A の固有ベクトルになる。

実行列 A に対しては、 U を実直交行列 Q に制限し、 S を対角ブロックの大きさが 2 以下のブロック上三角行列^{5} に緩和した形の実 Schur 分解

Q^{T} AQ = S

が存在する (実 Schur 分解)。 S の対角ブロックの固有値が A の固有値である。

(実) Schur 分解の存在は数学的に保証されているが、それを有限回の四則と開平演算で計算することはできない

Hermite 対称でもないのに、固有ベクトルではなく Schur ベクトルを計算するのは、これが数値的に安定に計算できるからだという。

  MATLAB では schur() という関数が用意されている。 [u,t]=schur(a) として u*t*u'-a

を計算すると…

11 Sylvester の慣性律, 2 次形式の対角化

(この節は下書きです。 )

数値線形代数というよりは、線形代数の話題かもしれないが…

11.1 本題に入る前に

自分自身が最初のうち景色が見えなかったので、色々書いてみる。

Sylvester の慣性律とは、 2 次形式を平方完成したときの係数の符号が、平方完成のやり方 (それはたくさんある) によらずに決まる、という事実を定式化したものである。

平方完成のことを「2 次形式の対角化」ともいう。それは議論を行列の言葉で書くことが出来て、与えられた対称行列 (2 次形式の係数行列) A に対して、適当な正則行列 P を見つけて、 P^{T} AP を対角行列にする操作に相当しているからである。

この変換 A ↦ P^{T} AP は、相似変換 A ↦ P^{−1}AP と似ているが、一応は別物であることに注意する (前者は正則な線形変数変換によるもので、後者は基底の取り替えによるものである)。ややこしいのは、対称行列の固有値問題では、直交行列 (P^{T} = P^{−1} を満たす行列) で対角化するため、両者が重なってしまうことである (個人的にはかなり混乱したところ)。

何に使われるかについても、あらかじめ、少しは知っておいた方が良いかもしれない。大学 1,2 年次に遭遇する適用例としては、多変数関数の極値問題を微分法で扱う際に、 2 階

までの Taylor 展開をしたとき、 2 次の項が 2 次形式になっている、というのがある (大抵の微積のテキストには書いてあるが、例えば桂田 [4])。 f が n 変数関数で、内点 a で極値を取るには、 f′(a) = 0 が必要で、そのとき

f(a + h) = f(a) + h の 2 次形式+ O(∥h∥^{2}) (h → 0).

5最初にこれを聞いたとき、 Hessenberg 行列のことかと思ったが、違う!

そこで 2 次形式の符号に興味が出て来るわけである。どういう場合があって、どうやれば判定できるか。他にも Morse の Lemma などで、関数の「様子」を理解するために使われる。

数値線形代数においては、固有値を求めるための 2 分法というアルゴリズムが、 Sylvester の慣性律で理解出来る、という応用もある。

もう少し式を使って具体的な話として書いてみる。

「2 次形式の対角化」は、次の二つの同値な問題である。

(a) (平方完成) 2 次形式 A[x] = ∑_{i,j=1}^{n} a_{ij}x_{i}x_{j} が与えられたとき、

A[x] = ∑_{j=1}^{p} α_{j}y_{j}^{2} − ∑_{k=1}^{q} α_{p+k}y_{p+k}^{2} (ここで α_{j} > 0)

となる正則線形変換 x = P y を求める問題

(b) (行列の「対角化」 ) 実対称行列 A が与えられたとき、

P^{T} AP = D,
D := diag(α_{1}, ⋯, α_{p}, −α_{p+1}, ⋯, −α_{p+q}, 0, ⋯, 0) (ここで α_{j} > 0)

を満たす正則行列 P を求める問題

p, q は P によらず A だけで定まる。さらに実は p + q = rank A である。

P として一般の正則変換を許せば (つまり P ∈ GL(n; R))、 α_{j} = 1 (j = 1, ⋯, p + q) と出来る。

P として直交変換に限っても対角化は可能で (もちろん p, q も変わらない)、そのときは D の対角成分 α_{1}, ⋯, α_{p}, −α_{p+1}, ⋯, −α_{p+q}, 0, ⋯, 0 は A の固有値である。

三角行列の積の形に分解する話との関係。

• もしも δ_{k} := det A_{k} ≠ 0 (k = 1, ⋯, n) ならば、 L := P^{−T} (:= (P^{−1})^{T}) は単位下三角行

列に取れて、

A = LDL^{T}

となる (A の LDL^{T} 分解)。 L は例えば Gauss の消去法などのアルゴリズムで計算出来る。

• (さらに) もしも δ_{k} > 0 (k = 1, ⋯, n) ならば、\sqrt{D} が意味を持つので、 L\sqrt{D} を新たに

L と置き直すことで

A = LL^{T}

となる (A の Cholesky 分解)。

det A_{k} ≠ 0 (k = 1, ⋯, n) でない場合にどうするかについて述べてみよう。その後で何をしたいかによってやることが異なる、と理解すべきである。

連立 1 次方程式 Ax = b を解くために A を何らかの意味で分解するのが目的ならば、適当な置換行列 P (上の議論と文字がかぶるけど大目に見て下さい) と、単位下三角行列 L が取れて

P A = LDL^{T}

となる、というのが一つの目標となるだろう。行列の符号の判別をしたい場合は、

(i) A が正値であるためには、ピボッティングなしの Gauss の消去法で対角線から下を消去できて、対角成分がすべて正数になることが必要十分。

(ii) A が負値であるためには、ピボッティングなしの Gauss の消去法で対角線から下を消去できて、対角成分がすべて負数になることが必要十分。

(iii) A が正値でも負値でもなく、さらに行列式が 0 でないならば、 A は不定符号であるが、これは逆は成り立たない。

(iv) 一般の場合に A の符号数 (p, q) を求めるには、少し工夫が必要である。伊理 [2] または [3] を見よ。

11.1.1 行列についての色々な同値関係

1. A, B ∈ M(n; C) に対して、

A ∼_{s} B ⇔^{def.} ∃U ∈ U(n; C) s.t. U^{∗}AU = B (U^{−1}AU = B).

代表元としては Schur 標準形というのか？上三角に出来る (もとが対称であれば対角行列に出来る)。実数バージョンもある。

2. (今回は関係ない) A, B ∈ M(n; C) に対して、

A ∼_{特} B ⇔^{def.} ∃U_{1}, U_{2} ∈ U(n; C) s.t. U_{2}^{∗}AU_{1} = B.

特異値分解に対応。

3. A, B ∈ M(n; C) に対して、

A ∼_{s} B (A と B は相似) ⇔^{def.} ∃P ∈ GL(n; C) s.t. P^{−1}AP = B.

代表元としては Jordan 標準形が有名。実数の範囲で考えると、有理標準形とかいうのになる？「基底の取り替えでなるべく簡単な行列にする」 4. A と B が n 次の実対称行列とするとき、

A ∼_{c} B (A と B は合同) ⇔^{def.} ∃P ∈ GL(n; k) s.t. P^{T} AP = B.

代表元として diag(\overbrace{1, ⋯, 1}^{p}, \overbrace{−1, ⋯, −1}^{q}, 0, ⋯, 0) が取れる。伊理先生は “Sylvester 形” と呼んでいた。これは符号数 (p, q) で指定できることに注意。「2 次形式を正則線形変換で平方完成 (対角化) する」 Hermite 行列版がある。

5. k = R or C とする。 A, B ∈ M(n; k) に対して、 A ∼_{r} B ⇔^{def.} ∃P, Q ∈ GL(n; k) s.t. QAP = B.

代表元として diag(\overbrace{1, ⋯, 1}^{r}, 0, ⋯, 0) が取れる (1 の個数 r は実は rank A に等しい, 伊理先生はこの代表元を「階数標準形 (rank normal form) と呼んでおく」 )。

実対称行列や Hermite 行列に対して、直交変換や unitary 変換で相似変換する話は、 1, 3, 4 にまたがっている話になる。

2 はある意味で 1 を緩くしたものである。 3, 4 もある意味で 1 を緩くしたものである。 U ∈ U(n) を G ∈ GL(n) に変えるのだが、そ

のときに U^{∗}(= U^{−1}) を P^{T} にするか、 P^{−1} にするか。 5 はとても緩い。例えば実対称行列について、 A ∼_{c} B ⇒ A ∼_{r} B. そうか、実対称行列について、直交変換に限った 1 (あるいは 3) をして、一般の正則線形変

換で 4 をして、それから 5 をして、ということで、標準形が、

diag(λ_{1}, ⋯, λ_{n}) \{ \begin{matrix}λ_{j} > 0 & (1 ≤ j ≤ p) \\ λ_{j} < 0 & (p + 1 ≤ j ≤ p + q) \\ λ_{j} = 0 & (p + q < j ≤ n) \end{matrix}

から

diag(\overbrace{1, ⋯, 1}^{p}, \overbrace{−1, ⋯, −1}^{q}, 0, ⋯, 0)

を通じて

diag(\overbrace{1, ⋯, 1}^{rank A}, 0, ⋯, 0)

となるわけだ。もう p + q = rank A は当たり前に見えて来る。

書くとすごくゴチャゴチャしている感じがするけれど、頭の中でやると訳が分からなくなるのは仕方ないのが分かる。

絵的に書くと

Q^{T} AQ = B ⇒ P^{T} AP = B ⇒ QAP = B.

A ∼_{s} diag(λ_{1}, ⋯, λ_{n}) ∼_{c} diag(1, ⋯, 1, −1, ⋯, −1, 0, ⋯, 0) ∼_{r} diag(1, ⋯, 1, 0, ⋯, 0).

11.2 見掛けが少し異なる 3 つの Sylvester の慣性律

以下に現れる行列は (時々「実」を書き漏らすけれど) すべて実行列とする。実対称行列 A に対して、 A の固有値のうちで正であるものの個数、負であるものの個数、

0 であるものの個数をそれぞれ π(A), ν(A), ζ(A) と書くことにする。 A は n 次正方行列として、

π(A), ν(A), ζ(A) ≥ 0, π(A) + ν(A) + ζ(A) = n

である。

まず実対称行列の実直交行列による対角化に関する常識的事項 (線形代数で習うはず) を復習する (これがなくても Sylvester の慣性律の証明は出来るわけだが、使った方が私には見通しが良いので…)。

A が n 次実対称行列ならば、実直交行列 U (U^{−1} = U^{T}) と実数 λ_{1}, ⋯, λ_{n} が存在して、

U^{T} AU = diag (λ_{1}, ⋯, λ_{n}) .

ゆえに ∀x ∈ R^{n} に対して y := U^{−1}x とおくと、 x = Uy で、

(Ax, x) = (AUy, Uy) = (U^{T} AUy, y) = ∑_{j=1}^{n} λ_{j}y_{j}^{2}.

こうして 2 次形式 (Ax, x) を平方和 — 定数× (何か)^{2} の和 — に変形できる。平方完成と呼んでも良いだろうが、これを 2 次形式の対角化という。標語的に

平方完成は対角化である

と言っても良いであろう (筋の通った主張にするには色々言葉を足さなければいけないが、迷子にならないために手短に言い切っておく)。

ところで、上の対角化の議論の要点は

(i) U^{T} AU が対角行列 diag (λ_{1}, ⋯, λ_{n}) である

(ii) U は正則

の二つである、と言えるだろう。

(Ax, x) が平方和の形に書けるためには、 U が実直交行列であることは必要ない (だから以下では U でなく、 P と書いたりする)。

U が実直交行列でなければ、 λ_{j} は A の固有値とは限らないが、 λ_{j} (j = 1, ⋯, n) のうちの正であるものの個数、負であるものの個数は、 A で定まる。それを定式化したものが、次の Sylvester の慣性律と呼ばれる定理である。

[figure]

定理 11.1 (Sylvester の慣性律, (Sylvester’s law of inertia)) 任意の実対称行列 A に対して、正則行列 P があって、 P^{T} AP = diag (λ_{1}, ⋯, λ_{n}) となるとき、

n_{p} := ♯{j ∈ {1, ⋯, n}; λ_{j} > 0}, n_{n} := ♯{j ∈ {1, ⋯, n}; λ_{j} < 0}, n_{z} := ♯{j ∈ {1, ⋯, n}; λ_{j} = 0} (♯ は集合の要素数を表す)

は P に依らず、 A だけで定まる。 (実対称行列が正則行列で対角化されるとき、対角成分のうち正であるものの個数、負であるものの個数、 0 であるものの個数は、正則行列によらず、 A だけで定まる。 )

[figure]

この定理は、少し後で紹介する「良くテキストに載っている証明」からすると自然な主張であるが、やや正体が分かりづらいのでは？と思う。次の定理とセットにして覚えることを勧める。

[figure]

定理 11.2 (老婆心版 Sylvester の慣性律)任意の実対称行列 A に対して、正則行列 P があって、 P^{T} AP = diag (λ_{1}, ⋯, λ_{n}) となるとき、

n_{p} := ♯{j ∈ {1, ⋯, n}; λ_{j} > 0},
n_{n} := ♯{j ∈ {1, ⋯, n}; λ_{j} < 0},
n_{z} := ♯{j ∈ {1, ⋯, n}; λ_{j} = 0}

とおくと、

n_{p} = π(A), n_{n} = ν(A), n_{z} = ζ(A). [figure]

証明定理 11.1 を認めた上での証明を与える。 P として U^{T} AU = diag[λ_{1}, ⋯, λ_{n}] となるような実直交行列 U を取れば、 P^{T} AP = U^{T} AU = U^{−1}AU は A の相似変換であり、 λ_{1}, ⋯, λ_{n} は A の固有値である。ゆえに n_{p}, n_{n}, n_{z} はそれぞれ π(A), ν(A), ζ(A) に等しい。

杉原・室田 [8] には、次の形の Sylvester の慣性律が載っている。

[figure]

定理 11.3 (杉原・室田版 Sylvester の慣性律) A を実対称行列とするとき，任意の正則行列 P に対して、 B := P^{T} AP とおくと、

π(A) = π(B), ν(A) = ν(B), ζ(A) = ζ(B)

が成り立つ。

[figure]

結論部分が「〜に依らず、〜で定まる」でなくて、具体的な等式である点は、定理 11.2 と同じである。対角化でなくて、変換した行列 B の固有値の話にしてしまうのも、少なくとも彼らの目的 (固有値の数値計算アルゴリズムの議論をする) にとっては使いやすいようである。

ここでは、定理 11.2 と定理 11.3 が同等であること (一方を認めれば他方がすぐに導かれること) を見てみよう。

まず、定理 11.3 を認めよう。 n 次実対称行列 A, 正則行列 P に対して、 ∃λ_{1}, ⋯, λ_{n} s.t. P^{T} AP = diag[λ_{1}, ⋯, λ_{n}] となったとすると、 B := P^{T} AP の固有値は (対角成分である) λ_{1}, ⋯, λ_{n} であるから、

n_{p} = π(B) = π(A), n_{n} = ν(B) = ν(A), n_{z} = ζ(B) = ζ(A).

逆に定理 11.2 を認めよう。実対称行列 A, 正則行列 P に対して、 B := P^{T} AP とおく。 B は実対称行列であるから、適当な実直交行列 U で対角化できる: ∃λ_{1}, ⋯, λ_{n} ∈ R s.t. U^{T} BU = diag[λ_{1}, ⋯, λ_{n}]. このとき

diag (λ_{1}, ⋯, λ_{n}) = U^{T} BU = U^{T} P^{T} AP U = (P U)^{T} A(P U) = P^{′T} AP′,
P′ := P U.

P と U が正則であるから、 P′ は正則である。定理 11.2 から、 (B の固有値である) λ_{1}, ⋯, λ_{n} のうちで正であるものの個数、負であるものの個数、 0 であるものの個数は、それぞれ π(A), ν(A), ζ(A) に等しい。ゆえに π(B) = π(A), ν(B) = ν(A), ζ(B) = ζ(A).

11.3 Sylvester の慣性律の証明

それでは Sylvester の慣性律を証明しよう。前項で 3 つのバージョンの同値性は分かっているので、どれか一つだけ証明すれば十分であるが、独立に証明することを目指す。

まず定理 11.1 を証明しよう。

定理 11.1 の証明線型形式の独立性 (行列のランクというべき？ ) の議論をしている。 x = Sy, x = Rz という正則変数変換で

(Ax, x) = α_{1}y_{1}^{2} + ⋯ + α_{r}y_{r}^{2} − α_{r+1}y_{r+1}^{2} − ⋯ − α_{p}y_{p}^{2} (α_{j} > 0)
= β_{1}z_{1}^{2} + ⋯ + β_{s}z_{s}^{2} − β_{s+1}z_{s+1}^{2} − ⋯ − β_{q}z_{q}^{2} (β_{j} > 0)

と対角化できたとする。一般性を失うことなく p ≥ q として良い。目標は p = q, r = s を証明することである。

各 y_{j}, z_{k} は変数 x についての線形形式と見なせる。 y_{j} 同士、 z_{k} 同士は 1 次独立である。もしも r < s ならば、 y_{1}, ⋯, y_{p} は y_{1}, ⋯, y_{r}, z_{s+1}, ⋯, z_{q} の 1 次結合では表せない。 (も

し出来たとすると、独立な p 個の 1 次形式 y_{1}, ⋯, y_{p} が、それより少ない r + (q − s) 個の 1 次形式で表せることになって矛盾する。 r + (q − s) = (r − s) + q < q ≤ p に注意)。例えば y_{k}, r + 1 ≤ k ≤ p が y_{1}, ⋯, y_{r}, z_{s+1}, ⋯, z_{q} で表せなかったとすると、

y_{1} = ⋯ = y_{r} = z_{s+1} = ⋯ = z_{q} = 0, y_{k} = 1

は x_{1}, ⋯, x_{n} に関する連立 1 次方程式として解を持つ。しかし、

(Ax, x) = α_{1} · 0^{2} + ⋯ + α_{r} · 0^{2} − α_{r+1}y_{r+1}^{2} − ⋯ − α_{k} · 1^{2} − ⋯ α_{p}y_{p}^{2} ≤ −α_{k} < 0,
(Ax, x) = β_{1}z_{1}^{2} + ⋯ + β_{s}z_{s}^{2} − β_{s+1} · 0^{2} − ⋯ − β_{q} · 0^{2} ≥ 0

となって矛盾する。

同様に r > s ならば、 y_{1}, ⋯, y_{r} が y_{r+1}, ⋯, y_{p}, z_{1}, ⋯, z_{s} で表されない。例えば y_{k} (1 ≤ k ≤ r) が y_{r+1}, ⋯, y_{p}, z_{1}, ⋯, z_{s} で表されないとすると、

y_{r+1} = ⋯ = y_{p} = z_{1} = ⋯ = z_{s} = 0, y_{k} = 1

は x_{1}, ⋯, x_{n} に関する連立 1 次方程式として解を持つ。しかし、

(Ax, x) = α_{1} · y_{1}^{2} + ⋯ + α_{k} · 1^{2} + ⋯ α_{r} · y_{r}^{2} − α_{r+1} · 0^{2} − ⋯ − α_{p} · 0^{2} ≥ α_{k} > 0,
(Ax, x) = β_{1} · 0^{2} + ⋯ + β_{s} · 0^{2} − β_{s+1} · y_{s+1}^{2} − ⋯ − β_{q} · β_{q}^{2} ≤ 0

となって矛盾する。

ゆえに r = s. 必要ならば (−Ax, x) を考えることによって、負の係数を持つ項の個数が等しいことも示せる。

上で書いたように、実は p = q = rank A であるから、それを利用するともう少し見通しの良い証明が得られる。証明 (準備中)

この証明はやや分り辛い (少なくとも私にとって)。そこでこれとは独立に定理 11.2 の証明を与える (これは自前ででっち上げたものなので、後でもう一度チェックすること)。

定理 11.2 の証明 A を n 次実対称行列とする。

P := {V; V は R^{n} の部分空間, ∀x ∈ V \ {0} (Ax, x) > 0} ,
N := {V; V は R^{n} の部分空間, ∀x ∈ V \ {0} (Ax, x) < 0} ,
Z := {V; V は R^{n} の部分空間, ∀x ∈ V \ {0} (Ax, x) = 0}

とおく。 {0} ∈ P, N, Z であるから、 P, N , Z ≠ ∅.

N_{p} := max{dim V; V ∈ P}, N_{n} := max{dim V; V ∈ N }, N_{z} := max{dim V; V ∈ Z}

が定まる (念のため: dim{0} = 0 である)。主張: N_{p} = π(A), N_{n} = ν(A), N_{z} = ζ(A). A は実対称行列であるから、固有ベクトルからなる正規直交基底 v_{1}, ⋯, v_{n} が存在する。

Av_{j} = λ_{j}v_{j} として、

V_{p} := Span{v_{j}; λ_{j} > 0}, V_{n} := Span{v_{j}; λ_{j} < 0}, V_{z} := Span{v_{j}; λ_{j} = 0}

とおくと、 V_{p} ∈ P, V_{n} ∈ N, V_{z} ∈ Z. dim V_{p} = π(A), dim V_{n} = ν(A), dim V_{z} = ζ(A) であるから、 N_{p}, N_{n}, N_{z} の最大性によって

(♯)  N_{p} ≥ π(A), N_{n} ≥ ν(A), N_{z} ≥ ζ(A).

さて、一般に

W_{p} ∈ P, W_{n} ∈ N , W_{z} ∈ Z ⇒ W_{p} ∩ W_{n} = W_{N} ∩ W_{z} = W_{z} ∩ W_{p} = {0}

が成り立つことは容易に証明できる。 W_{p}, W_{n}, W_{z} として、特に

dim W_{p} = N_{p}, dim W_{n} = N_{n}, dim W_{z} = N_{z}

を満たすものを取って、

W_{p} ⊕ W_{n} ⊕ W_{z} ⊂ R^{n}

において、次元を調べる。 (♯) を用いると

n = π(A) + ν(A) + ζ(A) ≤ N_{p} + N_{n} + N_{z} ≤ n.

左辺と右辺が一致するので、不等式はすべて等式で

N_{p} = π(A), N_{n} = ν(A), N_{z} = ζ(A) (主張の証明終わり).

さて、正則行列 P による変換 x = P y で、

(Ax, x) = ∑_{j=1}^{n} μ_{j}y_{j}^{2}

になったとする。簡単のために順番を修正して、最初の p 個は正，次の q 個は負、残りは 0, と出来る。つまり

(Ax, x) = α_{1}y_{1}^{2} + ⋯ + α_{p}y_{p}^{2} − α_{p+1}y_{p+1}^{2} − ⋯ − α_{p+q}y_{p+q}^{2}, α_{j} > 0 (1 ≤ j ≤ p + q)

として良い。

(a) p 次元空間 {P y; y_{p+1} = ⋯ = y_{n} = 0} では (Ax, x) > 0 になるので、 p ≤ π(A).

(b) q 次元空間 {P y; y_{1} = ⋯ = y_{p} = y_{p+q+1} = ⋯ = y_{n} = 0} では (Ax, x) < 0 になるので、 q ≤ ν(A).

(c) n−(p+q) 次元空間 {Sy; y_{1} = ⋯ = y_{p+q} = 0} では (Ax, x) = 0 になるので、 n−(p+q) ≤ ζ(A).

(c) から p + q + ζ(A) ≥ n であるが、 (a) から π(A) ≥ p, (b) から ν(A) ≥ q であるから、

n = π(A) + ν(A) + ζ(A) ≥ p + q + ζ(A) ≥ n.

両端が等しいので、途中はすべて等式で p = π(A), q = ν(A).

定理 11.3 の直接的な証明が見たければ、杉原・室田 [8] を見よ (もちろん前項の議論から、間接的には証明は済んでいる)。

11.4 2 次形式の対角化のアルゴリズム (後でどこかで使うもの)

2 次形式の対角化 (平方完成)、すなわち与えられた実対称行列 A に対して、 R^{T} AR = D, D = 対角行列となる正則行列 R を求めるアルゴリズムについて考える。これが解決すれば、 2 次形式の符号数 (n_{+}, n_{−}, n_{0}) を求める 1 つのアルゴリズムが得られるわけである。

古典的な 2 次形式論では、 Lagrange の方法と言うのだそうである。でも折角だから行列の変形操作の話として書き表したい (線形代数、あるいは数値線形代数として説明したいので)。

伊理 [2], [3] に書いてあるはずのこと。 11.4.1 例に基づく説明

例 11.4 (係数行列の前進消去だけで済む場合 (det A_{k} ≠ 0 のとき))まずピボットに 0 が現れない単純な例から。

A = (\begin{matrix}1 & −2 & 1 \\ −2 & 6 & −6 \\ 1 & −6 & 12 \end{matrix})

に対して Gauss の消去法の前進消去過程で、

A = (\begin{matrix}1 & −2 & 1 \\ −2 & 6 & −6 \\ 1 & −6 & 12 \end{matrix}) → (\begin{matrix}1 & −2 & 1 \\ 0 & 2 & −4 \\ 0 & −4 & 11 \end{matrix}) → (\begin{matrix}1 & −2 & 1 \\ 0 & 2 & −4 \\ 0 & 0 & 3 \end{matrix}) .

これはつまり

L := (\begin{matrix}1 & 0 & 0 \\ −2 & 1 & 0 \\ 1 & −2 & 1 \end{matrix}) , U := (\begin{matrix}1 & −2 & 1 \\ 0 & 2 & −4 \\ 0 & 0 & 3 \end{matrix})

と置くと

L^{−1}A = U, すなわち A = LU

である、ということを意味している (A の LU 分解)。手で検算することを考えて

L^{−1} = (\begin{matrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ 3 & 2 & 1 \end{matrix})

を残しておこう (私は軟弱なので MATLAB で計算している)。

閑話休題。連立 1 次方程式を解くならば、これから後退代入の過程に入るわけだが、今の目標は 2 次形式の対角化である。行に関して行った操作を列に関しても行うことで、対角線の上側も消去されて

L^{−1}A (L^{−1})^{T} = (\begin{matrix}1 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 3 \end{matrix})

となる ((L^{−1})^{T} は単位上三角行列なので、上三角行列 U = (L^{−1}A) に右から掛けると対角線の上側しか変化しない。一方で、結果が対称行列になることが分かるので、実は対角行列である。結局、対角線の上側をゼロ・クリアするだけである。 )。

D := (\begin{matrix}1 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 3 \end{matrix})

とおくと、

L^{−1}A(L^{−1})^{T} = D あるいは A = LDL^{T}

を意味している。

理解できていれば、 A を前進消去して上三角行列 U に変換した段階で、 L と D が求まる (求まったことが分かる)。念のため LDL^{T} を計算して A になることを確認すれば完璧だろう。 R := (L^{−1})^{T} とおき、 x = Ry と変数変換すると、

A[x] = (Ax, x) = (ARy, Ry) = (R^{T} ARy, y) = (Dy, y) = D[y]

となるはず。

ここで遠い記憶が蘇る。中国からの留学生 C さんのチューターをしていたとき、 C さんが僕の知らない計算をしていたのである。

(A | I)

と A と I を並べて、掃き出し (前進消去) をしていく。

• A だけでなく I も並べて同時に変形していくこと

• 左側のブロックについては、行について行った変形を列についてもしていく (上の例題では行・列の交換が起こらないのでさぼっても良いのだけど)

の二点が要点。

(L^{−1}A | L^{−1})

となるわけで、 U(= L^{−1}A) だけでなく、 L^{−1} の具体形も得られることに注意しよう。行列 A (あるいは 2 次形式 A[x] = (Ax, x)) の符号数が知りたいだけならば必要ないが、対角化するための変数変換 x = Ry, R = (L^{−1})^{T} が具体的に知りたい場合は、便利であろう (L^{−1} があれば、その転置行列はもう分っているに等しい)。つまり中国では、こういう筆算も教えているのでしょう。似たようなのを古屋先生の本か何かで見たような気がする。いつの間にか、日本の線形代数のテキストでは落ちてしまったわけだね。

上の例でこれを実行してみると

(A | I) = (\begin{matrix}1 & −2 & 1 & 1 & 0 & 0 \\ −2 & 6 & −6 & 0 & 1 & 0 \\ 1 & −6 & 12 & 0 & 0 & 1 \end{matrix}) → (\begin{matrix}1 & −2 & 1 & 1 & 0 & 0 \\ 0 & 2 & −4 & 2 & 1 & 0 \\ 0 & −4 & 11 & −1 & 0 & 1 \end{matrix}) → (\begin{matrix}1 & −2 & 1 & 1 & 0 & 0 \\ 0 & 2 & −4 & 2 & 1 & 0 \\ 0 & 0 & 3 & 3 & 2 & 1 \end{matrix}) .

これから

D = (\begin{matrix}1 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 3 \end{matrix}) , L^{−1} = (\begin{matrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ 3 & 2 & 1 \end{matrix}) .

確かに

L^{−1}A(L^{−1})^{T} = (\begin{matrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ 3 & 2 & 1 \end{matrix}) (\begin{matrix}1 & −2 & 1 \\ −2 & 6 & −6 \\ 1 & −6 & 12 \end{matrix}) (\begin{matrix}1 & 2 & 3 \\ 0 & 1 & 2 \\ 0 & 0 & 1 \end{matrix}) = (\begin{matrix}1 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 3 \end{matrix}) = D.

めでたし、めでたし。例 11.5 (楽屋裏の紹介と例の追加)上の例は逆算して作った。 L と D を適当に決めて、 LDL^{T} を計算して A を決めたのである。この例は正値対称行列だったけれど、不定符号が欲しければ、

[figure]

  L=[1,0,0;2,1,0;3,2,1]
  d=diag([1,-2,3])
  a=L*d*L'

[figure]

として

A = (\begin{matrix}1 & 2 & 3 \\ 2 & 2 & 2 \\ 3 & 2 & 4 \end{matrix})

とか。

(A | I) = (\begin{matrix}1 & 2 & 3 & 1 & 0 & 0 \\ 2 & 2 & 2 & 0 & 1 & 0 \\ 3 & 2 & 4 & 0 & 0 & 1 \end{matrix}) → (\begin{matrix}1 & 2 & 3 & 1 & 0 & 0 \\ 0 & −2 & −4 & −2 & 1 & 0 \\ 0 & −4 & −5 & −3 & 0 & 1 \end{matrix}) → (\begin{matrix}1 & 2 & 3 & 1 & 0 & 0 \\ 0 & −2 & −4 & −2 & 1 & 0 \\ 0 & 0 & 3 & 1 & −2 & 1 \end{matrix})

であるから

D = (\begin{matrix}1 & 0 & 0 \\ 0 & −2 & 0 \\ 0 & 0 & 3 \end{matrix}) , L^{−1} = (\begin{matrix}1 & 0 & 0 \\ −2 & 1 & 0 \\ 1 & −2 & 1 \end{matrix}) .

実際

L^{−1}A (L^{−1})^{T} = (\begin{matrix}1 & 0 & 0 \\ −2 & 1 & 0 \\ 1 & −2 & 1 \end{matrix}) (\begin{matrix}1 & 2 & 3 \\ 2 & 2 & 2 \\ 3 & 2 & 4 \end{matrix}) (\begin{matrix}1 & −2 & 1 \\ 0 & 1 & −2 \\ 0 & 0 & 1 \end{matrix}) = (\begin{matrix}1 & 0 & 0 \\ 0 & −2 & 0 \\ 0 & 0 & 3 \end{matrix}) = D.

例 11.6 (ピボットの交換が必要な場合 (∀k det A_{k} ≠ 0 が成り立たないとき))

A = (\begin{matrix}0 & 1 \\ 1 & 2 \end{matrix})

を考えよう。 (1, 1) 成分が 0 なので、このままでは掃き出せない。これも連立 1 次方程式を解くだけならば、第 1 行と第 2 行を交換して、つまり

P := (\begin{matrix}0 & 1 \\ 1 & 0 \end{matrix})

として、

P A = (\begin{matrix}1 & 2 \\ 0 & 1 \end{matrix})

を LU 分解すれば良い。我々の目標は対角化なので (少しくどいかな？ )、

Ã:= P AP^{T} = P AP = (\begin{matrix}2 & 1 \\ 1 & 0 \end{matrix})

という対称行列を考える。この LU 分解は

L := (\begin{matrix}1 & 0 \\ 1/2 & 1 \end{matrix}) , L^{−1}Ã = (\begin{matrix}2 & 1 \\ 0 & −1/2 \end{matrix}) .

ゆえに Ã の対角化は (対角線の上側をゼロ・クリアした)

L^{−1}Ã(L^{−1})^{T} = (\begin{matrix}2 & 0 \\ 0 & −1/2 \end{matrix}) .

検算してみよう。

D := (\begin{matrix}2 & 0 \\ 0 & −1/2 \end{matrix}) と置くと LDL^{T} = (\begin{matrix}2 & 1 \\ 1 & 0 \end{matrix}) = Ã.

Ã = P AP^{T}, L^{−1}Ã(L^{−1})^{T} = D であるから、

L^{−1}P AP^{T} (L^{−1})^{T} = D.

ゆえに

R := P^{T} (L^{−1})^{T} = ⋯ (ちびっと計算)⋯ = (\begin{matrix}0 & 1 \\ 1 & −1/2 \end{matrix})

とおくと、

R^{T} AR = D.

ゆえに

x = Ry すなわち(\begin{matrix}x_{1} \\ x_{2} \end{matrix}) = (\begin{matrix}y_{2} \\ y_{1} − y_{2}/2 \end{matrix})

と変数変換すると対角化されるはずである。確認しよう。

A[x] = 2x_{1}x_{2} + 2x_{2}^{2}
= 2y_{2}(y_{1} − y_{2}/2) + 2(y_{1} − y_{2}/2)^{2} = 2y_{1}y_{2} − y_{2}^{2} + 2(y_{1}^{2} − y_{1}y_{2} + y_{2}^{2}/4)
= 2y_{1}^{2} − y_{2}^{2}/2
= D[y].

とりあえず、めでたし、めでたし。

中国流の筆算を考えてみよう。左のブロック (もともと A の部分) には、行に関する基本変形と列に関する基本変形の両方を行うが (列に関する変形は要するに対称になるような操作をするわけで、 0 に掃き出したところと対称な場所にあるところは 0 にする、行交換をしたら列交換もする、だけであり、計算不要とも言える)、右のブロック (もともと I の部分) には、行に関する基本変形のみを行う。つまり

(A | I) → (R_{1}AR_{1}^{T} | R_{1}I) → (R_{2}R_{1}AR_{1}^{T} R_{2}^{T} | R_{2}R_{1}I) → ⋯
→ (R_{k} ⋯ R_{2}R_{1}AR_{1}^{T} R_{2}^{T} ⋯ R_{k}^{T} | R_{k} ⋯ R_{2}R_{1}I)
= (D | R) , D = R_{k} ⋯ R_{2}R_{1}AR_{1}^{T} R_{2}^{T} ⋯ R_{k}^{T} , R := R_{k} ⋯ R_{2}R_{1}.

これは

R^{T} AR = D

を意味している。上の A に対して実行すると、これだけで良い！

[figure]

\begin{aligned}(A | I) & = ( \begin{matrix}0 & 1 & 1 & 0 \\ 1 & 2 & 0 & 1 \end{matrix} ) → ( \begin{matrix}1 & 2 & 0 & 1 \\ 0 & 1 & 1 & 0 \end{matrix} ) → ( \begin{matrix}2 & 1 & 0 & 1 \\ 1 & 0 & 1 & 0 \end{matrix} ) \\ & → ( \begin{matrix}2 & 0 & 0 & 1 \\ 0 & −1/2 & 1 & −1/2 \end{matrix} ) . \end{aligned}

[figure]

これから、

D = (\begin{matrix}2 & 0 \\ 0 & −1/2 \end{matrix}) , R^{T} = (\begin{matrix}0 & 1 \\ 1 & −1/2 \end{matrix}) , R = (\begin{matrix}0 & 1 \\ 1 & −1/2 \end{matrix}) .

こうして上で「(ちびっと計算)」と書いたところが、枠内の計算で済んでしまった。

例 11.7 (対角成分がすべて 0 — 自明でない最も簡単な場合)A = (\begin{matrix}0 & a \\ a & 0 \end{matrix}) とする。 A の対角成分はすべて 0 なので、これまでと同じようには掃き出せない。連立 1 次方程式を解く場合は、単に第 1 行と第 2 行を交換すれば良いが、対角化をする場合はそれは役に立たない。 (P = (\begin{matrix}0 & 1 \\ 1 & 0 \end{matrix}) とするとき、 P A = (\begin{matrix}a & 0 \\ 0 & a \end{matrix}) であるが、右から P^{T} をかける必要があり、そう

すると P AP^{T} = (\begin{matrix}0 & a \\ a & 0 \end{matrix}) = A と元に戻ってしまう。 )

ここでは、直交行列による対角化を思い出そう (必ず出来るはずだ)。固有値は ±a で、固有ベクトルはそれぞれ (\begin{matrix}1 \\ 1 \end{matrix}), (\begin{matrix}−1 \\ 1 \end{matrix}) である。そこで R = \frac{1}{\sqrt{2}} (\begin{matrix}1 & −1 \\ 1 & 1 \end{matrix}) とすると、 R は原点の回りの角 π/4 の回転であり、 R^{T} = R^{−1} は角 −π/4 の回転である。

R^{T} AR = (\begin{matrix}a & 0 \\ 0 & −a \end{matrix}) .

あるいは R = \frac{1}{\sqrt{2}} (\begin{matrix}1 & 1 \\ 1 & −1 \end{matrix}) とすると、 R^{T} = R^{−1} = R. これは鏡映変換で、こちらを使っても

R^{T} AR = (\begin{matrix}a & 0 \\ 0 & −a \end{matrix}) .

対角化が目的の場合、必ずしも固有値を求める必要はないので、\frac{1}{\sqrt{2}} はなくても良い。 R =

(\begin{matrix}1 & −1 \\ 1 & 1 \end{matrix}) とすれば、

R^{T} AR = (\begin{matrix}2a & 0 \\ 0 & −2a \end{matrix}) .

あるいは R = (\begin{matrix}1 & 1 \\ 1 & −1 \end{matrix}) として

R^{T} AR = (\begin{matrix}2a & 0 \\ 0 & −2a \end{matrix}) .

後者は対称行列なので、扱いやすいかもしれない。中国流筆算では、

(A | I) = (\begin{matrix}0 & a & 1 & 0 \\ a & 0 & 0 & 1 \end{matrix}) → (R^{T} AR | R^{T}) = (\begin{matrix}2a & 0 & 1 & 1 \\ 0 & −2a & 1 & −1 \end{matrix})

とするところだろう。変数変換は x = Ry, すなわち

(\begin{matrix}x_{1} \\ x_{2} \end{matrix}) = (\begin{matrix}1 & 1 \\ 1 & −1 \end{matrix}) (\begin{matrix}y_{1} \\ y_{2} \end{matrix}) = (\begin{matrix}y_{1} + y_{2} \\ y_{1} − y_{2} \end{matrix})

で、実際に

A[x] = 2ax_{1}x_{2}
= 2a(y_{1} + y_{2})(y_{1} − y_{2}) = 2a(y_{1}^{2} − y_{2}^{2})
= 2ay_{1}^{2} − 2ay_{2}^{2}

と対角化が出来る。

なお、

R = (\begin{matrix}1 & 1 \\ −1 & 1 \end{matrix})

の逆行列は

R^{−1} = \frac{1}{2} (\begin{matrix}1 & 1 \\ −1 & 1 \end{matrix}) = \frac{1}{2}R.

上の例は簡単であるが、一般の場合に基本操作として使うことが出来る。それを説明しよう。

A = ( \begin{matrix}D & O & O \\ O & A′ & B^{T} \\ O & B & C \end{matrix} ) , D = 対角行列, A′ = (\begin{matrix}0 & a \\ a & 0 \end{matrix}) , C = 対称行列

とする (つまり掃き出しの過程で対角成分が 0 のブロックが出て来た場合を考える)。

S := ( \begin{matrix}I & O & O \\ O & R & O \\ O & O & I \end{matrix} ) , R := (\begin{matrix}1 & 1 \\ 1 & −1 \end{matrix})

とすると、 R^{T} = R であるから、 S^{T} = S で、

S^{T} AS = ( \begin{matrix}D & O & O \\ O & R^{T} A′R & R^{T} B^{T} \\ O & BR & C \end{matrix} ) , R^{T} A′R = (\begin{matrix}2a & 0 \\ 0 & −2a \end{matrix}) .

ピボットが見つからないという障害物を突破できた。すなわち

(⋆)  S^{T} AS = ( \begin{matrix}D & O & O \\ O & \begin{matrix}2a & 0 \\ 0 & −2a \end{matrix} & (BR)^{T} \\ O & BR & C \end{matrix} ) .

後で使いそうなので

B = (\begin{matrix}b_{11} & b_{12} \\ b_{21} & b_{22} \\ ⋮ & ⋮ \\ b_{r1} & b_{r2} \end{matrix}) とするとき BR = (\begin{matrix}b_{11} + b_{12} & b_{11} − b_{12} \\ b_{21} + b_{22} & b_{21} − b_{22} \\ ⋮ & ⋮ \\ b_{r1} + b_{r2} & b_{r1} − b_{r2} \end{matrix})

であることを注意しておく。 (⋆) の形になれば、 2 段階の掃き出しが出来て

S^{T} AS → ( \begin{matrix}D & O & O \\ O & \begin{matrix}2a & 0 \\ 0 & −2a \end{matrix} & O \\ O & O & C′ \end{matrix} ) .

少し脇道にそれるが、後のために注意しておくと

S^{−1} = ( \begin{matrix}I & O & O \\ O & \frac{1}{2}R & O \\ O & O & I \end{matrix} )

である。

例 11.8

A = ( \begin{matrix}1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 2 & 1 & 2 \\ 0 & 2 & 0 & 3 & 1 \\ 0 & 1 & 3 & 1 & 3 \\ 0 & 2 & 1 & 3 & 2 \end{matrix} ) .

(A | I) = (\begin{matrix}1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 2 & 1 & 2 & 0 & 1 & 0 & 0 & 0 \\ 0 & 2 & 0 & 3 & 1 & 0 & 0 & 1 & 0 & 0 \\ 0 & 1 & 3 & 1 & 3 & 0 & 0 & 0 & 1 & 0 \\ 0 & 2 & 1 & 3 & 2 & 0 & 0 & 0 & 0 & 1 \end{matrix}) → (\begin{matrix}1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 4 & 0 & 4 & 3 & 0 & 1 & 1 & 0 & 0 \\ 0 & 0 & −4 & −2 & 1 & 0 & 1 & −1 & 0 & 0 \\ 0 & 4 & −2 & 1 & 3 & 0 & 0 & 0 & 1 & 0 \\ 0 & 3 & 1 & 3 & 2 & 0 & 0 & 0 & 0 & 1 \end{matrix})

→ (\begin{matrix}1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 4 & 0 & 4 & 3 & 0 & 1 & 1 & 0 & 0 \\ 0 & 0 & −4 & −2 & 1 & 0 & 1 & −1 & 0 & 0 \\ 0 & 0 & 0 & −2 & −1/2 & 0 & −3/2 & −1/2 & 1 & 0 \\ 0 & 0 & 1 & 0 & −1/4 & 0 & −3/4 & −3/4 & 0 & 1 \end{matrix})

→ (\begin{matrix}1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 4 & 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 \\ 0 & 0 & −4 & −2 & 1 & 0 & 1 & −1 & 0 & 0 \\ 0 & 0 & 0 & −2 & −1/2 & 0 & −3/2 & −1/2 & 1 & 0 \\ 0 & 0 & 1 & 0 & −1/4 & 0 & −3/4 & −3/4 & 0 & 1 \end{matrix})

→ (\begin{matrix}1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 4 & 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 \\ 0 & 0 & −4 & −2 & 1 & 0 & 1 & −1 & 0 & 0 \\ 0 & 0 & 0 & −2 & −1/2 & 0 & −3/2 & −1/2 & 1 & 0 \\ 0 & 0 & 0 & −1/2 & 0 & 0 & −3/4 & −3/4 & 0 & 1 \end{matrix})

→ (\begin{matrix}1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 4 & 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 \\ 0 & 0 & −4 & 0 & 0 & 0 & 1 & −1 & 0 & 0 \\ 0 & 0 & 0 & −2 & −1/2 & 0 & −3/2 & −1/2 & 1 & 0 \\ 0 & 0 & 0 & −1/2 & 0 & 0 & −3/4 & −3/4 & 0 & 1 \end{matrix})

→ (\begin{matrix}1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 4 & 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 \\ 0 & 0 & −4 & 0 & 0 & 0 & 1 & −1 & 0 & 0 \\ 0 & 0 & 0 & −2 & −1/2 & 0 & −3/2 & −1/2 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1/8 & 0 & −1/8 & −7/8 & −1/4 & 1 \end{matrix})

→ (\begin{matrix}1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 4 & 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 \\ 0 & 0 & −4 & 0 & 0 & 0 & 1 & −1 & 0 & 0 \\ 0 & 0 & 0 & −2 & 0 & 0 & −3/2 & −1/2 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1/8 & 0 & −1/8 & −7/8 & −1/4 & 1 \end{matrix}) .

これから

D = (\begin{matrix}1 & 0 & 0 & 0 & 0 \\ 0 & 4 & 0 & 0 & 0 \\ 0 & 0 & −4 & 0 & 0 \\ 0 & 0 & 0 & −2 & 0 \\ 0 & 0 & 0 & 0 & 1/8 \end{matrix}) , S^{T} = (\begin{matrix}1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 1 & 0 & 0 \\ 0 & 1 & −1 & 0 & 0 \\ 0 & −3/2 & −1/2 & 1 & 0 \\ 0 & −1/8 & −7/8 & −1/4 & 1 \end{matrix}) .

実際

S^{T} AS = D

が成り立つ。

より一般の場合、すなわち

A = ( \begin{matrix}D & O \\ O & A′ \end{matrix} ) , D = k − 1 次対角行列, A′ = 対称行列

で、 A′ の対角成分はすべて 0 であるが、非対角成分に 0 でないものがある場合を考える。まずは a_{ij} ≠ 0 となる (i, j) ∈ {k, ⋯, n} × {k, ⋯, n} を探す。

[figure]

    found = 0;
    for (i=k; i<=n; i++) {
      for (j=k; j<=n; j++) { // j=i+1 からで良いのかな？
       if (a[i][j] != 0) {
         found = 1;
         break;
       }
      }
    }

[figure]

a_{ij} = a_{ji} = a ≠ 0 として、 i 行と k 行、 i 列と k 列を交換、 j 行と k + 1 行、 j 列と k + 1 列を交換、すると、

A → P_{j,k+1}P_{ik}AP_{ik}P_{j,k+1} = ( \begin{matrix}D & O & O \\ O & \begin{matrix}0 & a \\ a & 0 \end{matrix} & (B′)^{T} \\ O & B′ & C′ \end{matrix} )

の形になる。ここで P_{rs} は左からかけると、 r 行と s 行を交換することになる置換行列とする。

11.4.2 プログラム例

/*
* testsign-v1.c
*/
#include <stdio.h>
#include <math.h>
#define MAXN (100)
void printmatrix(double a[MAXN+1][MAXN+1], int n)
{
int i, j;
for (i=1; i<=n; i++) {
for (j=1;j<=n; j++) {
printf("%g ", a[i][j]);
}
printf("\n");
}
}
void find_sign(double a[MAXN+1][MAXN+1], int n,
int *np, int *nn, int *nz)
{
int i,j,k,p,ll,found,verbose=1;
double c,alpha,beta,t;
for (k=1; k<n; k++) {

// 対角成分で 0 でないものを探す

found = 0;
for (p=k; p<=n; p++)
if (a[p][p] != 0) {
found = 1;
break;
}
if (found) {
if (p!=k) {
if (verbose) {

printf("ピボットは 0 だけど、対角成分で 0 でないものを見つけた。 \n");

  printmatrix(a,n);
}
printf("%d 行,%d 行を交換\n", k, p);

// p 行と k 行を交換

for (ll=k;ll<=n;ll++) {
  t=a[p][ll]; a[p][ll]=a[k][ll]; a[k][ll]=t;
}
printf("%d 列,%d 列を交換\n", k, p);

// p 列と k 列を交換

for (ll=k; ll<=n; ll++) {
  t=a[ll][p]; a[ll][p]=a[ll][k]; a[ll][k]=t;
}
if (verbose) {
  printmatrix(a,n);
}
}

// p==k のときは、もともと a[k][k]!=0 であるわけで交換の必要はない // a[k][k] をピボットとして掃き出す

    for (i=k+1; i<=n; i++) {
      c=a[i][k]/a[k][k];

a[i][k]=0; // さぼっても良い。 c を書いておいて後で使う手も

      for (j=k+1;j<=n;j++)
        a[i][j] -= c*a[k][j];
    }
  }

else { // 非対角成分で 0 でないものを探す

    found = 0;
    for (i=k;i<=n;i++) {
      for (j=k;j<=n;j++) { // 対称なんだし右上だけ: j=i+1 からで良いな
        if (a[i][j] != 0) {
          found = 1;
          break;
        }
      }
if (found) break;
    }
    if (found) {
if (verbose) {

printf("非対角成分で 0 でないものを見つけた。 \n");

printmatrix(a, n);
}
      if (!(i==k && j==k+1)) {

// i 行と k 行, i 列と k 列を交換

printf("%d 行 (列),%d 行 (列) を交換\n", k, i);
for (ll=k;ll<=n;ll++) {
  t=a[i][ll]; a[i][ll]=a[k][ll]; a[k][ll]=t;
}
for (ll=k; ll<=n; ll++) {
  t=a[ll][i]; a[ll][i]=a[ll][k]; a[ll][k]=t;
}
// j 行と k+1 行, j 列と k+1 列を交換
printf("%d 行 (列),%d 行 (列) を交換\n", k+1, j);
for (ll=k;ll<=n;ll++) {
  t=a[j][ll]; a[j][ll]=a[k+1][ll]; a[k+1][ll]=t;
}
for (ll=k; ll<=n; ll++) {
  t=a[ll][j]; a[ll][j]=a[ll][k+1]; a[ll][k+1]=t;
}
printmatrix(a, n);
}
      // a[k][k+1] == a[k+1][k] != 0 となっている
      a[k][k] = a[k][k+1];
      a[k+1][k+1] = - a[k][k+1];
      a[k][k+1] = a[k+1][k] = 0;
      for (j=k+2; j<=n; j++) {
        alpha = a[k][j]; beta = a[k+1][j];
        a[k  ][j] = a[j][k  ] = alpha + beta;
        a[k+1][j] = a[j][k+1] = alpha - beta;
      }

// a[k][k] をピボットとして掃き出す

      for (i=k+1; i<=n; i++) {
        c=a[i][k]/a[k][k];
        a[i][k]=0;
        for (j=i+1;j<=n;j++)
          a[i][j] -= c*a[k][j];
      }
      // a[k+1][k+1] をピボットとして掃き出す
      for (i=k+2; i<=n; i++) {
        c=a[i][k+1]/a[k+1][k+1];
        a[i][k+1]=0;
        for (j=i+1;j<=n;j++)
          a[i][j] -= c*a[k+1][j];
      }
    }
    else {

// もう右下は 0 である

      break;
    }
  }
}

// 対角成分のうち正、負、零のものの個数を数える

*np=0;
*nn=0;
*nz=0;
for (k=1;k<=n;k++)
  if (a[k][k] > 0)
    (*np)++;
  else if (a[k][k] < 0)
    (*nn)++;
  else
    (*nz)++;
}
int main(void)
{
int n, np, nn, nz;
double a[MAXN+1][MAXN+1];
// 例題 1
n=3;
a[1][1] =  1; a[1][2] = -2; a[1][3] =  1;
a[2][1] = -2; a[2][2] =  6; a[2][3] = -6;
a[3][1] =  1; a[3][2] = -6; a[3][3] = 12;
printf("例題 1\n");
printmatrix(a,n);
find_sign(a, n, &np, &nn, &nz);
printmatrix(a,n);
printf("np=%d, nn=%d, nz=%d\n", np, nn, nz);
// 例題 2
n=3;
a[1][1] =  1; a[1][2] = 2; a[1][3] = 3;
a[2][1] =  2; a[2][2] = 2; a[2][3] = 2;
a[3][1] =  3; a[3][2] = 2; a[3][3] = 4;

printf("\n 例題 2\n");
printmatrix(a,n);
find_sign(a, n, &np, &nn, &nz);
printmatrix(a,n);
printf("np=%d, nn=%d, nz=%d\n", np, nn, nz);
// 例題 3
n=2;
a[1][1] =  0; a[1][2] = 1;
a[2][1] =  1; a[2][2] = 2;
printf("\n 例題 3\n");
printmatrix(a,n);
find_sign(a, n, &np, &nn, &nz);
printmatrix(a,n);
printf("np=%d, nn=%d, nz=%d\n", np, nn, nz);
// 例題 4
n=2;
a[1][1] = 0; a[1][2] = 2;
a[2][1] = 2; a[2][2] = 0;
printf("\n 例題 4\n");
printmatrix(a,n);
find_sign(a, n, &np, &nn, &nz);
printmatrix(a,n);
printf("np=%d, nn=%d, nz=%d\n", np, nn, nz);
// 例題 5
n=3;
a[1][1] = 0; a[1][2] = 0; a[1][3] = 0;
a[2][1] = 0; a[2][2] = 0; a[2][3] = 2;
a[3][1] = 0; a[3][2] = 2; a[3][3] = 1;
printf("\n 例題 5\n");
printmatrix(a,n);
find_sign(a, n, &np, &nn, &nz);
printmatrix(a,n);
printf("np=%d, nn=%d, nz=%d\n", np, nn, nz);
// 例題 6
n=3;
a[1][1] = 0; a[1][2] = 0; a[1][3] = 0;
a[2][1] = 0; a[2][2] = 0; a[2][3] = 2;
a[3][1] = 0; a[3][2] = 2; a[3][3] = 0;
printf("\n 例題 6\n");
printmatrix(a,n);
find_sign(a, n, &np, &nn, &nz);
printmatrix(a,n);
printf("np=%d, nn=%d, nz=%d\n", np, nn, nz);
return 0;
}

\overline{例題 1}

1 -2 1
-2 6 -6
1 -6 12
1 -2 1
0 2 -4
0 0 3
np=3, nn=0, nz=0

例題 2

1 2 3
2 2 2
3 2 4
1 2 3
0 -2 -4
0 0 3
np=2, nn=1, nz=0

例題 3

0 1
1 2

ピボットは 0 だけど、対角成分で 0 でないものを見つけた。

0 1
1 2

1 行,2 行を交換 1 列,2 列を交換

2 1
1 0
2 1
0 -0.5
np=1, nn=1, nz=0

例題 4

0 2
2 0

非対角成分で 0 でないものを見つけた。

0 2
2 0
2 0
0 -2
np=1, nn=1, nz=0

例題 5

0 0 0
0 0 2
0 2 1

ピボットは 0 だけど、対角成分で 0 でないものを見つけた。

0 0 0
0 0 2
0 2 1

1 行,3 行を交換 1 列,3 列を交換

1 2 0
2 0 0
0 0 0
1 2 0

0 -4 0
0 0 0
np=1, nn=1, nz=1

例題 6

0 0 0
0 0 2
0 2 0

非対角成分で 0 でないものを見つけた。

0 0 0
0 0 2
0 2 0
1 行 (列),2 行 (列) を交換
2 行 (列),3 行 (列) を交換
0 2 0
2 0 0
0 0 0
2 0 0
0 -2 0
0 0 0
np=1, nn=1, nz=1

A TODO

1. LR 法

2. Gerschgorin の定理の後半

3. Courant-Fischer の mini-max

B ChageLog

思い立ってこれからは加筆修正内容を記録することにした。

• 2005 年 8 月 24 日数日来余暇に遊んできた Lanczos 法の説明を書く。

• 2005 年 8 月 25 日木村 [5] を読んでいて、 Householder 変換による QR 分解の話が載っていたので、この際書くことにした。

• 2005 年 8 月 26 日 Householder 変換の原理を確認するための MATLAB プログラムを書いて実験する。

• 2005 年 8 月 26 日試しに既に書いてあった Jacobi 法のプログラムを実験してみたら、バグで動かなかったので、デバッグする。プログラム以外にも誤植があったので訂正する。やはり実験までしないとミスは減らない…

• 2005年8月27日昨日のデバッグで分かったこと: Octave では function [q r]=myqr(a) のような書き方ができるが、 MATLAB では function [q,r]=myqr(a) のようにカンマが必要。他のプログラムも修正。 • 2005 年 8 月 27 日本来は「本に書いてあること」を理解すべきだとは思うが、 QR 分

解をする各種手法の比較実験をしてみた。中途半端な実験しかできていないせいだろうとは思うがやや意外な結果 (修正 Gram-Schmidt の直交化法を用いる理由は見えなかった)。やってみて分かったのは、何を尺度にして評価するのか、実は分かっていなかったということ。

• 2005 年 9 月 6 日 QR 分解の各方法の比較の話、ちょっと気になって Trefethen and Bau III を見たら、ありゃりゃ。やはりこの本の内容は日本でも常識になっていないといけないのだろう。翻訳出さないかなあ。忙しいので、自分で試してここに書き足すのはいつの日か。

• 2013 年 8 月 26 日 Sylvester の慣性律について少し書く (Cholesky 分解ノートからこちらにコピーした)。やっと少し分かった。

• 2013 年 9 月 4 日いよいよ一般の場合の 2 次形式の対角化、符号数の求め方、の話に手をつける。伊理先生の本を読んでプログラムを書こうかと思ったが、まずは実例攻撃でいくことにした。そのうちに古い記憶が蘇った。

参考文献

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

[2] 伊理^{いり}正夫^{まさお}. 一般線形代数. 岩波書店, 2003. 伊理正夫, 線形代数 I, II, 岩波講座応用数学, 岩波書店 (1993, 1994) の単行本化.

[3] 伊理^{いり}正夫^{まさお}. 線形代数汎論. 朝倉書店, 2009. 「一般線形代数」のリニューアル.

[4] 桂田祐史. 多変数の微分積分学 1 講義ノート (2013 年度版). http://nalab.mind.meiji.

ac.jp/~mk/lecture/tahensuu1-2013/tahensuu1-new-text.pdf, 2013.

[5] 木村英紀^{ひでのり}. 線形代数数理科学の基礎. 東京大学出版会, 2003.

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

[7] 山本哲朗^{てつろう}. 数値解析入門 [新訂版]. サイエンス社, 2003. 1976 年初版発行の定番本の待望の改訂版.

[8] 杉原正顯^{まさあき}, 室田^{むろた}一雄. 線形計算の数理. 岩波書店, 2009.

[9] 齋藤^{さいとう}正彦^{まさひこ}. 線型代数入門. 東京大学出版会, 1966.

[10] 佐武一郎. 線型代数学. 裳華房, 1958, 1974. もともと『行列と行列式』という書名であったのを、テンソル代数の章を加筆した機会に改題した。

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

[12] Lloyd Nicholas Trefethen and David Bau III. Numerical Linear Algebra. SIAM, 1997.

[13] 桂田祐史. 一般化固有値問題. http://nalab.mind.meiji.ac.jp/~mk/labo/text/
    generalized-eigen/generalized-eigenvalue-problem.pdf, 2003.