拙作のページ「連立方程式の解法」で説明したガウスの消去法は、拡大係数行列を階段行列に変形することで解を求めています。このような方法を「直接法 (Direct Method)」といいます。これに対し、適当な初期解から始めて、繰り返し計算することで真の解に収束させていく方法が考えられます。これを「反復法 (Iterative Method)」といいます。
直接法は厳密解を求めることができますが、係数行列が大きくなると解くのに時間がかかるようになります。厳密解を求めるのが難しい場合、反復法を使うと現実的な時間で近似解を求めることが可能です。今回は基本的な反復法である「ヤコビ法」と「ガウス・ザイデル法」について簡単に説明します。
連立方程式 \(Ax = b\) の係数行列 \(A\) を対角行列 \(D\) とそれ以外の要素を持つ行列 \(A'\) に分解します。すると、方程式は次のように変形することができます。
ここで、最後の式を漸化式 \(x_{i+1} = D^{-1}(b - A'x_i)\) と考えて、反復処理により解を求めるのが「ヤコビ法 (Jacobi Method)」です。たとえば、三元連立方程式の場合、漸化式は下図のようになります。
a1 * x + a2 * y + a3 * z = d1
b1 * x + b2 * y + b3 * z = d2
c1 * x + c2 * y + c3 * z = d3
A' = [[ 0, a2, a3], D = [[a1, 0, 0], D-1 = [[1/a1, 0, 0],
[b1, 0, b3], [ 0, b2, 0], [ 0, 1/b2, 0],
{c1, c2, 0]] [ 0, 0, c3]] [ 0, 0, 1/c3]]
xi+1 = (d1 - a2 * yi - a3 * zi) / a1
yi+1 = (d2 - b1 * xi - b3 * zi) / b2
zi+1 = (d3 - c1 * xi - c2 * yi) / c3
ヤコビ法で解が収束した場合、それが元の連立方程式の解となります。なお、行列 \(A\) によっては解が収束しない (発散する) 場合があります。収束の十分条件ですが、『係数行列 A が対角優位な行列である場合に収束する』ことが知られています。行列 \(A\) の対角成分 \(a_{ii}\) の絶対値が、他の行の成分 \(a_{ij}\) の絶対値の合計よりも大きいことを「対角優位」といいます。
上図の三元連立方程式でいえば、|a1| > |b1| + |c1|, |b2| > |a2| + |c2|, |c3| > |a3| + |b3| を満たすとき、ヤコビ法の解は収束することになります。
収束の判定ですが、一般的には \(x_{i+1}\) と \(x_i\) の差分が許容誤差 \(\epsilon\) に収まったときに収束と判定します。具体的な方法を以下に示します。
1 は差分の合計値、2 は差分を解で割った値の合計値、3 は差分を解で割った値の最大値が \(\epsilon\) 以下になったとき収束と判定します。今回は 2 の方法を使うことにします。
それではプログラムを作りましょう。NumPy を使うとヤコビ法はとても簡単にプログラムすることができます。
リスト : ヤコビ法
def jacobi(a, b, max_iter=512, eps=1e-6):
d = np.diag(a)
a1 = a - np.diag(d)
x = np.full(len(b), 0.)
for i in range(max_iter):
x1 = (b - a1 @ x) / d
if np.fabs((x1 - x) / x1).sum() < eps:
return x1
x = x1
引数 a が係数行列、b が右辺値を格納した配列、キーワード引数 max_iter は繰り返しの最大回数、eps は許容誤差を表します。最初に、a の対角成分を diag() で取り出して変数 d にセットします。そして、diag(d) で対角行列を生成し、a から引き算すれば、対角成分を 0 にした行列 a1 を生成することができます。変数 x には解を格納する配列をセットします。x は 0 で初期化します。
あとは for ループで漸化式 x = D-1(b - A'x) を繰り返し計算します。NumPy の場合、行列 A, B の割り算 A / B は要素同士の割り算になるので、対角行列 d の逆行列を求める必要はありません。x1 = (b - a1 @ x) / d で漸化式を計算することができます。収束条件のチェックも簡単です。np.fabs((x1 - x) / x1).sum() で合計値を求め、それが eps よりも小さくなったならば解は収束したので x1 を返します。for ループが終了した場合、解は収束していないので None を返します。
それでは実行してみましょう。次に示す方程式をヤコビ法で解きます。なお、変数 x1 の途中経過を表示するようにプログラムを修正しています。
9x + y + 2z = 9 x + 9y + z = 18 2x + y + 9z = -5
0 [ 0. 0. 0.] 1 [ 1. 2. -0.55555556] 2 [ 0.90123457 1.95061728 -1. ] 3 [ 1.00548697 2.01097394 -0.97256516] 4 [ 0.99268404 1.99634202 -1.00243865] 5 [ 1.00094836 2.00108385 -0.99796779] 6 [ 0.99942797 1.99966882 -1.00033118] 7 [ 1.00011039 2.00010036 -0.99983609] 8 [ 0.99995242 1.99996952 -1.00003568] 9 [ 1.00001132 2.00000925 -0.99998604] 10 [ 0.99999587 1.99999719 -1.00000354] 11 [ 1.0000011 2.00000085 -0.99999877] 12 [ 0.99999963 1.99999974 -1.00000034] 13 [ 1.0000001 2.00000008 -0.99999989] [ 0.99999997 1.99999998 -1.00000003]
厳密解は [1, -2, -1] ですが、ヤコビ法だと 13 回で収束しています。ガウス・ザイデル法を使うと、これよりも少ない階数で解を求めることができます。
ヤコビ法で漸化式を計算するとき、右辺式は xi, yi zi の値を使いますが、1 行目から順番に計算していくと、2 行目の yi+1 を計算するときには xi+1 の値がすでに求まっています。3 行目の zi+1 を計算するときは、xi+1 と yi+1 の値が求まっています。これらの値を使って漸化式を計算する方法を「ガウス・ザイデル法 (Gauss-Siedel Method)」といいます。
たとえば、三元連立方程式の場合、ガウス・ザイデル法の漸化式は下図のようになります。
a1 * x + a2 * y + a3 * z = d1
b1 * x + b2 * y + b3 * z = d2
c1 * x + c2 * y + c3 * z = d3
A' = [[ 0, a2, a3], D = [[a1, 0, 0], D-1 = [[1/a1, 0, 0],
[b1, 0, b3], [ 0, b2, 0], [ 0, 1/b2, 0],
{c1, c2, 0]] [ 0, 0, c3]] [ 0, 0, 1/c3]]
xi+1 = (d1 - a2 * yi - a3 * zi) / a1
yi+1 = (d2 - b1 * xi+1 - b3 * zi) / b2
zi+1 = (d3 - c1 * xi+1 - c2 * yi+1) / c3
ガウス・ザイデル法の解の収束条件や収束の判定方法はヤコビ法と同じです。
それではプログラムを作りましょう。一般的なプログラミング言語の場合、ヤコビ法よりもガウス・ザイデル法のほうが簡単にプログラムできるのですが、NumPy の場合は逆にちょっとだけ複雑になります。
リスト : ガウス・ザイデル法
def gauss_seidel(a, b, max_iter=512, eps=1e-6):
d = np.diag(a)
a1 = a - np.diag(d)
x = np.full(len(b), 0.)
for i in range(max_iter):
err = 0.
for j in range(len(x)):
s = x[j]
x[j] = (b[j] - a1[j] @ x) / d[j]
err += abs((x[j] - s) / x[j])
if err < eps: return x
ポイントは 2 番目の for ループで配列 x の値を逐次的に更新していくところです。これで a1[j] @ x を計算するときに計算済みの値を用いることができます。変数 err には誤差の合計値を格納します。2 番目の for ループが終了したとき、err が eps よりも小さい場合、解は収束したので x を返します。
それでは実際に試してみましょう。
9x + y + 2z = 9 x + 9y + z = 18 2x + y + 9z = -5
0 [ 0. 0. 0.] 1 [ 1. 1.88888889 -0.98765432] 2 [ 1.00960219 1.99756135 -1.00186286] 3 [ 1.00068493 2.00013088 -1.00016675] 4 [ 1.00002251 2.00001603 -1.00000678] 5 [ 0.99999973 2.00000078 -1.00000003] [ 0.99999992 2.00000001 -0.99999998]
ヤコビ法だと収束するまで 13 回かかっていたのが、ガウス・ザイデル法を使うと 5 回ですみました。
それでは簡単なテストとして、ガウスの消去法、ヤコビ法、ガウス・ザイデル法の実行時間を比較してみましょう。次のリストを見てください。
リスト : 簡単なテスト
# テストデータの作成
def test(n):
a = np.random.rand(n * n).reshape((n, n))
for i in range(n):
a[i, i] = (a[:, i].sum()) * 2
x = np.random.rand(n)
return a, x, a @ x
for n in [500, 1000, 1500, 2000]:
a, x, b = test(n)
print("-----", n, "-----")
for f, n in [(gauss, "gauss"), (jacobi, "jacobi"), (gauss_seidel, "gauss_seidel")]:
print(n)
s = time.time()
x1 = f(a, b)
e = time.time()
print(np.fabs(x1 - x).sum(), e - s)
乱数でテストデータを生成します。これを関数 test() で行っています。連立方程式 Ax = b の係数行列 A と解 x を乱数で生成して、右辺値の b を計算します。係数行列は対角優位になるよう修正しています。test() は A, x, b を返すので、A と b から x を求め、その誤差と実行時間を表示します。
実行結果は次のようになりました。
----- 500 ----- gauss 2.474088642290262e-13 0.43669724464416504 jacobi 5.408217758177958e-08 0.006071567535400391 gauss_seidel 6.633193446771435e-09 0.05102801322937012 ----- 1000 ----- gauss 7.061761657205234e-13 1.9147076606750488 jacobi 2.9102833972263244e-08 0.031321048736572266 gauss_seidel 2.1125195586975745e-09 0.11282467842102051 ----- 1500 ----- gauss 1.3118729193239975e-12 4.963993549346924 jacobi 4.275096058495992e-08 0.06757640838623047 gauss_seidel 3.088161337110387e-09 0.14022326469421387 ----- 2000 ----- gauss 1.9380514174830582e-12 10.164307355880737 jacobi 2.87893372736105e-08 0.14386820793151855 gauss_seidel 4.147552099860322e-09 0.1823873519897461
変数の個数が増えるにしたがい、ガウスの消去法の実行時間は大幅に増加します。反復法でも遅くはなりますが、ガウスの消去法よりはずっと高速です。ヤコビ法とガウス・ザイデル法を比較した場合、ヤコビ法のほうが少し速くなりました。Python の実行速度は遅いので、二重の for ループになるガウス・ザイデル法よりも、NumPy で高速に処理できるヤコビ法のほうが速くなるのでしょう。興味のある方はいろいろ試してみてください。
「SOR 法 (Successive Over-Relaxation)」は、ヤコビ法の漸化式に加速パラメータ \(w\) を導入することで収束を速くする方法です。具体的には、ヤコビ法で求めた値 \(x_{i+1}\) と \(x_{i}\) の差分に \(w\) を掛け算し、その値を \(x_{i}\) に足し算します。SOR 法の漸化式は次のようになります。
漸化式は参考 URL 6 『SOR法 - Wikipedia』を参考にさせていただきました。また、\(w\) が \(0 \lt w \lt 2\) のとき SOR 法は収束するそうです。\(w\) が 1 のときはヤコビ法と同じになります。一般に、\(w\) を 1 より大きくすると収束が速くなり、1 より小さくすると遅くなるといわれています。ところが、参考 URL 3 『数値演算法 (8) 連立方程式を解く -2-』によると、1 より小さくしたほうが収束が速くなる場合もあるそうです。これは後で試してみましょう。連立方程式によって w の最適値は変化するでしょうから、w の選択はけっこう難しい問題なのかもしれません。
それではプログラムを作りましょう。次のリストを見てください。
リスト : SOR 法
# ヤコビ法
def jacobi_sor(a, b, w, max_iter=512, eps=1e-6):
d = np.diag(a)
a1 = a - np.diag(d)
x = np.full(len(b), 0.)
for i in range(max_iter):
x1 = (1 - w) * x + w * (b - a1 @ x) / d
if np.fabs((x1 - x) / x1).sum() < eps:
return x1
x = x1
# ガウス・ザイデル法
def gauss_seidel_sor(a, b, w, max_iter=512, eps=1e-6):
d = np.diag(a)
a1 = a - np.diag(d)
x = np.full(len(b), 0.)
for i in range(max_iter):
err = 0.
for j in range(len(x)):
s = x[j]
x[j] = (1 - w) * s + w * (b[j] - a1[j] @ x) / d[j]
err += abs((x[j] - s) / x[j])
if err < eps: return x
関数 jacobi_sor() と gauss_seidel_sor() の引数 w が加速パラメータを表します。どちらの関数も SOR 法の漸化式をそのままプログラムしただけなので、難しいところはないと思います。
それでは実行してみましょう。次に示す方程式を SOR 法で解いたときの収束回数を示します。
9x + y + 2z = 9 x + 9y + z = 18 2x + y + 9z = -5
w | 0.9 0.925 0.95 0.975 1.0 1.025 1.05 1.075 1.1 -------+-------------------------------------------------------- Jacobi | 12 12 11 12 13 14 15 16 18 Gauss | 8 8 7 6 5 5 6 7 8
ヤコビ法の場合、w が 0.95 のとき、収束回数が一番少なくなりました。ガウス・ザイデル法の場合、このような簡単な方程式では SOR 法の効果はほとんどないようです。
次は test() で作成した方程式を解いてみましょう。テストプログラムと実行結果を示します。
リスト : 簡単なテスト
for n in [500, 1000, 1500, 2000]:
a, x, b = test(n)
print("-----", n, "-----")
for f, n, w in [(jacobi_sor, "jacobi", 0.9), (gauss_seidel_sor, "gauss_seidel", 0.95)]:
print(n)
s = time.time()
x1 = f(a, b, w)
e = time.time()
print(np.fabs(x1 - x).sum(), e - s)
----- 500 ----- jacobi 8.137453068557355e-09 0.007251262664794922 gauss_seidel 5.574765369627645e-09 0.05086326599121094 ----- 1000 ----- jacobi 1.6065243091129694e-08 0.0250089168548584 gauss_seidel 1.069149200031913e-08 0.11444592475891113 ----- 1500 ----- jacobi 2.3646492177862265e-08 0.08397293090820312 gauss_seidel 1.0884738805812602e-09 0.14935922622680664 ----- 2000 ----- jacobi 1.1291291129194615e-08 0.13947463035583496 gauss_seidel 1.3888327852942372e-09 0.18320822715759277
ヤコビ法の場合、\(w\) を 0.9 にすると収束回数が 10 回以上減少し、その分だけ実行速度も速くなりました。ガウス・ザイデル法の場合、\(w\) を 0.95 にすると収束回数は数回減少しますが、実行速度はほとんど変わりませんでした。今回のテストでは \(w\) を 1 より増やすと、どちらの方法でも収束回数は大幅に増加します。\(w\) の設定は教科書通りにはいきませんね。収束の判定方法や許容誤差の値を変えると、異なる結果になるかもしれません。興味のある方はいろいろ試してみてください。
#
# jacobi.py : 連立方程式の解法 (反復法)
#
# Copyright (C) 2018-2023 Makoto Hiroi
#
import numpy as np
import time
# ガウスの消去法
def gauss(xs, ys):
# 拡大係数行列の生成
n = len(xs)
zs = np.c_[xs.astype(np.float_), ys.astype(np.float_)]
# 前進消去
for i in range(n - 1):
for j in range(i + 1, n):
temp = zs[j, i] / zs[i, i]
zs[j, i:] -= temp * zs[i, i:]
# 後退代入
for i in range(n - 1, -1, -1):
zs[i, n] -= zs[i, i+1:n] @ zs[i+1:, n]
zs[i, n] /= zs[i, i]
return zs[:, n]
# ヤコビ法
def jacobi(a, b, max_iter=512, eps=1e-6):
d = np.diag(a)
a1 = a - np.diag(d)
x = np.full(len(b), 0.)
for i in range(max_iter):
x1 = (b - a1 @ x) / d
if np.fabs((x1 - x) / x1).sum() < eps:
return x1
x = x1
# ガウス・ザイデル法
def gauss_seidel(a, b, max_iter=512, eps=1e-6):
d = np.diag(a)
a1 = a - np.diag(d)
x = np.full(len(b), 0.)
for i in range(max_iter):
err = 0.
for j in range(len(x)):
s = x[j]
x[j] = (b[j] - a1[j] @ x) / d[j]
err += abs((x[j] - s) / x[j])
if err < eps: return x
# テストデータの作成
def test(n):
a = np.random.rand(n * n).reshape((n, n))
for i in range(n):
a[i, i] = (a[:, i].sum()) * 2
x = np.random.rand(n)
return a, x, a @ x
#
# SOR 法
#
# ヤコビ法
def jacobi_sor(a, b, w, max_iter=512, eps=1e-6):
d = np.diag(a)
a1 = a - np.diag(d)
x = np.full(len(b), 0.)
for i in range(max_iter):
x1 = (1 - w) * x + w * (b - a1 @ x) / d
if np.fabs((x1 - x) / x1).sum() < eps:
return x1
x = x1
# ガウス・ザイデル法
def gauss_seidel_sor(a, b, w, max_iter=512, eps=1e-6):
d = np.diag(a)
a1 = a - np.diag(d)
x = np.full(len(b), 0.)
for i in range(max_iter):
err = 0.
for j in range(len(x)):
s = x[j]
x[j] = (1 - w) * s + w * (b[j] - a1[j] @ x) / d[j]
err += abs((x[j] - s) / x[j])
if err < eps: return x
今回は「固有値 (eigenvalue)」と「固有ベクトル (eigenvector)」について取り上げます。固有値と固有ベクトルは行列の性質を表す重要な指標のひとつで、行列を対角行列に変換する「対角化」の基礎になります。
対角行列は単純な構造をしているので、いろいろな計算が簡単になるという利点があります。たとえば、対角行列 \(A\) の逆行列は対角成分を逆数にしたものですし、\(A\) の n 乗は対角成分を n 乗したものになります。\(A\) の行列式も対角成分を掛け算するだけで求めることができます。
与えられた行列の固有値と固有ベクトルを求める問題のことを「固有値問題」といいます。NumPy には固有値と固有ベクトルを求める関数 linalg.eig() が用意されていますが、Python (NumPy) とアルゴリズムのお勉強ということで、あえてプログラムを作ってみましょう。
まず最初に、特殊な行列とその性質について簡単に説明します。行列 \(A\) の行と列を入れ替えたものを「転置行列 (transposed matrix)」といい \(A^{\mathrm{T}}\) と表記します。NumPy では A.T で転置行列を求めることができます。簡単な例を示しましょう。
>>> a = np.array([[1, 2], [3, 4]])
>>> a.T
array([[1, 3],
[2, 4]])
>>> b = np.arange(9).reshape((3, 3))
>>> b
array([[0, 1, 2],
[3, 4, 5],
[6, 7, 8]])
>>> b.T
array([[0, 3, 6],
[1, 4, 7],
[2, 5, 8]])
>>> c = np.arange(6).reshape((3, 2))
>>> c
array([[0, 1],
[2, 3],
[4, 5]])
>>> c.T
array([[0, 2, 4],
[1, 3, 5]])
転置行列の基本的な性質を以下に示します。A は正方行列を表します。
簡単な例を示しましょう。
>>> a = np.array([[1, 2], [3, 4]])
>>> a
array([[1, 2],
[3, 4]])
>>> np.trace(a)
5
>>> np.trace(a.T)
5
>>> np.linalg.det(a)
-2.0
>>> np.linalg.det(a.T)
-2.0
>>> b = np.array([[5, 6], [7, 8]])
>>> b
array([[5, 6],
[7, 8]])
>>> (a @ b).T
array([[19, 43],
[22, 50]])
>>> b.T @ a.T
array([[19, 43],
[22, 50]])
>>> np.linalg.inv(a.T)
array([[-2. , 1.5],
[ 1. , -0.5]])
>>> np.linalg.inv(a).T
array([[-2. , 1.5],
[ 1. , -0.5]])
正方行列 \(A\) とその転置行列 \(A^{\mathrm{T}}\) が等しい行列のことを「対称行列 (symmetric matrix)」といい、その中で対角成分 \(a_{ii}\) 以外の要素が 0 の行列のことを「対角行列 (diagonal matrix)」といいます。単位行列は対角行列の一種です。簡単な例を示しましょう。
>>> a = np.array([[1, 3], [3, 2]])
>>> a
array([[1, 3],
[3, 2]])
>>> a.T
array([[1, 3],
[3, 2]])
>>> b = np.array([[1, 4, 5], [4, 2, 6], [5, 6, 3]])
>>> b
array([[1, 4, 5],
[4, 2, 6],
[5, 6, 3]])
>>> b.T
array([[1, 4, 5],
[4, 2, 6],
[5, 6, 3]])
>>> c = np.array([[1, 0, 0], [0, 2, 0], [0, 0, 3]])
>>> c
array([[1, 0, 0],
[0, 2, 0],
[0, 0, 3]])
>>> c.T
array([[1, 0, 0],
[0, 2, 0],
[0, 0, 3]])
対角行列 A の基本的な性質を以下に示します。
4 はあとで説明します。簡単な例を示しましょう。
>>> c
array([[1, 0, 0],
[0, 2, 0],
[0, 0, 3]])
>>> np.linalg.det(c)
6.0
>>> np.diag(c).prod()
6
>>> c @ c
array([[1, 0, 0],
[0, 4, 0],
[0, 0, 9]])
>>> c @ c @ c
array([[ 1, 0, 0],
[ 0, 8, 0],
[ 0, 0, 27]])
>>> c ** 2
array([[1, 0, 0],
[0, 4, 0],
[0, 0, 9]])
>>> c ** 3
array([[ 1, 0, 0],
[ 0, 8, 0],
[ 0, 0, 27]])
>>> np.linalg.inv(c)
array([[ 1. , 0. , 0. ],
[ 0. , 0.5 , 0. ],
[ 0. , 0. , 0.33333333]])
>>> np.diag(1 / np.diag(c))
array([[ 1. , 0. , 0. ],
[ 0. , 0.5 , 0. ],
[ 0. , 0. , 0.33333333]])
\(A^{\mathrm{T}}\) と \(A^{-1}\) が等しい行列のことを「直交行列 (orthogonal matrix)」といいます。簡単な例を示しましょう。
>>> d = np.array([[0, 1], [1, 0]])
>>> d
array([[0, 1],
[1, 0]])
>>> d.T
array([[0, 1],
[1, 0]])
>>> np.linalg.inv(d)
array([[ 0., 1.],
[ 1., 0.]])
>>> e = np.array([[0, 1, 0], [1, 0, 0], [0, 0, 1]])
>>> e
array([[0, 1, 0],
[1, 0, 0],
[0, 0, 1]])
>>> e.T
array([[0, 1, 0],
[1, 0, 0],
[0, 0, 1]])
>>> np.linalg.inv(e)
array([[ 0., 1., 0.],
[ 1., 0., 0.],
[ 0., 0., 1.]])
d と e のように、各行各列にちょうど一つだけ 1 の要素を持ち、それ以外は全て 0 となるような正方行列を「置換行列」といいます。置換行列を掛け算すると、列または行を入れ替える働きをします。
直交行列の基本的な性質を以下に示します。
3 はあとで説明します。簡単な例を示しましょう。
回転行列 R は直交行列
R = [[cos(x), -sin(x)],
[sin(x), cos(x)]]
転置行列
RT = [[cos(x), sin(x)]
[-sin(x), cos(x)]]
行列式 |R| = cos(x) * cos(x) + sin(x) * sin(x) = 1
逆行列
R-1 = [[cos(x), sin(x)],
[-sin(x), cos(x)]]
逆行列の転置 (R に戻るので直交行列)
(R-1)T = [[cos(x), -sin(x)],
[sin(x), cos(x)]]
最初に固有値と固有ベクトルの基本を簡単に説明します。n 次の正方行列 \(A\) に対して、\(Ax = \lambda x\) を満たす数 \(\lambda\) を \(A\) の「固有値」、ゼロではないベクトル \(x\) を \(A\) の「固有ベクトル」といいます。このとき、固有値は \(|A - \lambda I| = 0\) の根になります。この式を「固有方程式」とか「特性方程式」といいます。
\(Ax = \lambda x\) を変形すると \((A - \lambda I)x = 0\) になります。行列 \((A - \lambda I)\) に逆行列が存在する場合、両辺にその逆行列を掛け算すると \(x = 0\) になってしまいます。つまり、ゼロではないベクトル \(x\) が存在するためには、逆行列が存在しないこと (\(|A - \lambda I| = 0\)) が条件になるわけです。
それでは、具体的に固有値と固有ベクトルを求めてみましょう。たとえば、行列 [[1, 2], [2, 1]] の固有値は以下のようになります。
A = [[1, 2], A - λI = [[1 - λ, 2],
[2, 1]] [2, 1 - λ]]
|A - λI| = (1 - λ) * (1 - λ) - 4
= λ2 - 2 λ - 3
= (λ- 3)(λ+ 1) = 0
固有値 λ= 3, -1
固有値が求まったら、連立方程式 \((A - \lambda I)x = 0\) を解いて固有ベクトル x を求めます。
(1) λ= 3 の場合 [[1 - 3, 2 ], * [a, b] = 0 [2, 1 - 3]] -2a + 2b = 0, 2a - 2b = 0 => a = b なので x = k * [1, 1] (k は任意の数, ただし k != 0) 単位ベクトル (大きさ 1 のベクトル) で表すと [1, 1] / √(12 +12) = [√2 / 2, √2 / 2] 実際に計算すると [[1, 2], * [√2 / 2, √2 / 2] = [√2 / 2 + 2*√2 / 2, 2*√2 / 2 + √2 / 2] [2, 1]] = 3 * [√2 / 2, √2 / 2] Ax = λx を満たす
(2) λ= -1 の場合 [[1 + 1, 2 ], * [a, b] = 0 [2, 1 + 1]] 2a + 2b = 0, 2a + 2b = 0 => a = -b なので x = k * [-1, 1] (k は任意の数, ただし k != 0) 単位ベクトルで表すと [-1, 1] / √(12 +12) = [-√2 / 2, √2 / 2] 実際に計算すると [[1, 2], * [-√2 / 2, √2 / 2] = [-√2 / 2 + 2*√2 / 2, -2*√2 / 2 + √2 / 2] [2, 1]] = -1 * [-√2 / 2, √2 / 2] Ax = λx を満たす
一般に、n 次の正方行列 \(A\) は n 個の固有値とそれに対応する固有ベクトルを持ちます。なお、固有方程式が重根を持つ場合、固有値の個数は n よりも少なくなります。
それでは実際に NumPy で固有値と固有ベクトルを求めてみましょう。関数 linalg.eig() を使います。
>>> import numpy as np
>>> np.linalg.eig(np.array([[1, 2], [2, 1]]))
(array([ 3., -1.]), array([[ 0.70710678, -0.70710678],
[ 0.70710678, 0.70710678]]))
>>> np.linalg.eig(np.array([[4, -2], [1, 1]]))
(array([ 3., 2.]), array([[ 0.89442719, 0.70710678],
[ 0.4472136 , 0.70710678]]))
>>> np.linalg.eig(np.diag([1, 2, 3]))
(array([ 1., 2., 3.]), array([[ 1., 0., 0.],
[ 0., 1., 0.],
[ 0., 0., 1.]]))
>>> np.linalg.eig(np.diag([4, 3, 2, 1]))
(array([ 4., 3., 2., 1.]), array([[ 1., 0., 0., 0.],
[ 0., 1., 0., 0.],
[ 0., 0., 1., 0.],
[ 0., 0., 0., 1.]]))
eig() の返り値を (xs, ys) とすると、固有値 xs[i] に対応する固有ベクトルは ys[:, i] になります。つまり、ys は固有ベクトルを列に持つ行列になります。そして、3, 4 番目の例のように、対角行列の固有値は対角成分の値となります。
n 次の正方行列 \(A\) の固有ベクトル \(x_{1}, x_{2}, ..., x_{n}\) を列に持つ行列 \(X = [x_{1}, x_{2}, ...., x_{n}]\) を考えます。すると、\(X\) によって行列 \(A\) を以下のように対角行列 \(\varLambda\) に変換することができます。
固有値と固有ベクトルの定義 \(Ax_i = \lambda_i x_i\) から \(AX\) は次のように変形することができます。
左側から両辺に \(X\) の逆行列を掛け算すれば \(X^{-1}AX = \varLambda\) となります。
それでは実際に NumPy で試してみましょう。
>>> a = np.array([[1, 2], [2, 1]])
>>> xs, ys = np.linalg.eig(a)
>>> xs
array([ 3., -1.])
>>> ys
array([[ 0.70710678, -0.70710678],
[ 0.70710678, 0.70710678]])
>>> np.linalg.inv(ys) @ a @ ys
array([[ 3.00000000e+00, 4.44089210e-16],
[ 6.66133815e-16, -1.00000000e+00]])
>>> b = np.array([[4, -2], [1, 1]])
>>> xs, ys = np.linalg.eig(b)
>>> xs
array([ 3., 2.])
>>> ys
array([[ 0.89442719, 0.70710678],
[ 0.4472136 , 0.70710678]])
>>> np.linalg.inv(ys) @ b @ ys
array([[ 3., 0.],
[ 0., 2.]])
数学 (線形代数) の世界では、正則行列 \(P\) を用いて正方行列 \(A\) を \(P^{-1}AP\) に変換することを「相似変換」といい、\(B = P^{-1}AP\) が成り立つ正方行列 \(A, B\) の関係を「相似」といいます。相似な行列 \(A\) と \(B\) の間は様々な性質が保存されます。保存される主な性質を以下に示します。
つまり、相似変換を行ってもこれらの値が変化することはありません。たとえば、正方行列 A を相似変換により対角行列に変換します。対角行列は対角成分が固有値になるので、トレースは固有値の総和になります。相似変換でトレースは保存されるので、行列 \(A\) のトレースも固有値の総和であることがわかります。
要素がすべて実数の対称行列を「実対称行列」といいます。実対称行列には次に示す性質があります。
3 は直交行列の定義 \(L^{\mathrm{T}} = L^{-1}\) から明らかです。2 は固有ベクトルが「一次独立」であることを言っています。
ベクトル \(x_{1}, ..., x_{i}, ..., x_{n}\) において、式 \(a_{1}x_{1} + \cdots + a_{i}x_{i} + \cdots + a_{n}x_{n} = 0\) が成り立つとき、係数 \(a_{i} = 0 \ (i = 1, \ldots, n)\) 以外の解がないことを一次独立といいます。逆に、係数 \(a_{i} \ne 0 \ (i = 1, \ldots, n)\) であっても式が 0 になることを「一次従属」といいます。
簡単な例を示しましょう。
[[1, 0, 0]
[0, 2, 0], の固有値 [1, 2, 3], 固有ベクトル [1, 0, 0], [0, 1, 0], [0, 0, 1]
[0, 0, 3]]
[1, [0, [0, [0,
a1 * 0, + a2 * 1, + a3 * 0 = 0,
0] 0] 1] 0]
a1 * 1 + a2 * 0 + a3 * 0 = 0 => a1 = 0
a1 * 0 + a2 * 1 + a3 * 0 = 0 => a2 = 0
a1 * 0 + a2 * 0 + a3 * 1 = 0 => a3 = 0
このように、実対称行列の固有ベクトルは一次独立になります。
それでは実際に NumPy で試してみましょう。
>>> a
array([[1, 2],
[2, 1]])
>>> xs, ys = np.linalg.eig(a)
>>> xs
array([ 3., -1.])
>>> ys
array([[ 0.70710678, -0.70710678],
[ 0.70710678, 0.70710678]])
>>> ys[:, 0] @ ys[:, 1]
0.0
>>> ys.T @ a @ ys
array([[ 3.00000000e+00, 4.44089210e-16],
[ 6.10622664e-16, -1.00000000e+00]])
>>> b
array([[1, 4, 5],
[4, 2, 6],
[5, 6, 3]])
>>> xs, ys = np.linalg.eig(b)
>>> xs
array([ 12.17597107, -2.50728797, -3.6686831 ])
>>> ys
array([[ 0.49659978, 0.80958546, -0.31298568],
[ 0.57735027, -0.57735027, -0.57735027],
[ 0.64811675, -0.10600965, 0.7541264 ]])
>>> ys[:, 0] @ ys[:, 1]
-9.7144514654701197e-17
>>> ys[:, 0] @ ys[:, 2]
-5.5511151231257827e-17
>>> ys[:, 1] @ ys[:, 2]
-4.8572257327350599e-16
>>> ys.T @ b @ ys
array([[ 1.21759711e+01, -2.99760217e-15, 0.00000000e+00],
[ -2.02615702e-15, -2.50728797e+00, 4.44089210e-16],
[ -8.88178420e-16, 1.11022302e-15, -3.66868310e+00]])
なお、NumPy の linalg.eig() は一般の行列の固有値と固有ベクトルを求める関数ですが、対称行列であれば関数 linalg.eigh() を使うこともできます。詳細はリファレンスマニュアル Linear algebra をお読みください。
>>> a
array([[1, 2],
[2, 1]])
>>> np.linalg.eigh(a)
(array([-1., 3.]), array([[-0.70710678, 0.70710678],
[ 0.70710678, 0.70710678]]))
>>> b
array([[1, 4, 5],
[4, 2, 6],
[5, 6, 3]])
>>> np.linalg.eigh(b)
(array([ -3.6686831 , -2.50728797, 12.17597107]), array([[-0.31298568, 0.80958546,
-0.49659978],
[-0.57735027, -0.57735027, -0.57735027],
[ 0.7541264 , -0.10600965, -0.64811675]]))
>>> c = np.ones((5, 5))
>>> c += np.diag([6, 7, 8, 9, 10])
>>> c
array([[ 7., 1., 1., 1., 1.],
[ 1., 8., 1., 1., 1.],
[ 1., 1., 9., 1., 1.],
[ 1., 1., 1., 10., 1.],
[ 1., 1., 1., 1., 11.]])
>>> np.linalg.eigh(c)
(array([ 6.27769582, 7.35663185, 8.43473667, 9.54039443, 13.39054123]),
array([[ 0.91678475, -0.21893977, 0.13337217, -0.09513681, 0.29108754],
[-0.35246548, -0.83284951, 0.22633151, -0.13258644, 0.33663729],
[-0.14781784, 0.46166517, 0.74694901, -0.21865948, 0.39908692],
[-0.09351905, 0.18073897, -0.57446875, -0.6232889 , 0.4899839 ],
[-0.06839508, 0.11236447, -0.20745782, 0.73284978, 0.63449885]]))