M.Hiroi's Home Page
http://www.geocities.jp/m_hiroi/

Python3 Programming

お気楽 NumPy プログラミング超入門

[ Home | Light | Python3 ]

●連立方程式の解法 (2)

拙作のページ 連立方程式の解法 で説明した「ガウスの消去法」は、拡大係数行列を階段行列に変形することで解を求めています。このような方法を「直接法 (Direct Method)」といいます。これに対し、適当な初期解から始めて、繰り返し計算することで真の解に収束させていく方法が考えられます。これを「反復法 (Iterative Method)」といいます。

直接法は厳密解を求めることができますが、係数行列が大きくなると解くのに時間がかかるようになります。厳密解を求めるのが難しい場合、反復法を使うと現実的な時間で近似解を求めることが可能です。今回は基本的な反復法である「ヤコビ法」と「ガウス・ザイデル法」について簡単に説明します。

●ヤコビ法

連立方程式 Ax = b の係数行列 A を対角行列 D とそれ以外の要素を持つ行列 A' に分解します。すると、方程式は次のように変形することができます。

Ax = b => (A' + D)x = b => Dx = b - A'x => x = D-1(b - A'x)

ここで、最後の式を漸化式 xi+1 = D-1(b - A'xi) と考えて、反復処理により解を求めるのが「ヤコビ法 (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 の対角成分 aii の絶対値が、他の行の成分 aij の絶対値の合計よりも大きいことを「対角優位」といいます。

|aii| > Σaij (i, j = 1, 2, ..., n, j != i)

上図の三元連立方程式でいえば、|a1| > |b1| + |c1|, |b2| > |a2| + |c2|, |c3| > |a3| + |b3| を満たすとき、ヤコビ法の解は収束することになります。

収束の判定ですが、一般的には xi+1 と xi の差分が許容誤差 ε に収まったときに収束と判定します。具体的な方法を以下に示します。

  1. Σ |xi+1 - xi| <= ε
  2. Σ |(xi+1 - xi) / xi+1| <= ε
  3. max |(xi+1 - xi) / xi+1| <= ε

1 は差分の合計値、2 は差分を解で割った値の合計値、3 は差分を解で割った値の最大値が ε 以下になったとき収束と判定します。今回は 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
3.12370788957e-13 0.3612337112426758
jacobi
2.70284279515e-08 0.01547098159790039
gauss_seidel
7.00478601555e-09 0.022022724151611328
----- 1000 -----
gauss
8.4839817463e-13 2.0641372203826904
jacobi
2.75979188328e-08 0.06420445442199707
gauss_seidel
2.03730985771e-09 0.06721687316894531
----- 1500 -----
gauss
1.64524860749e-12 5.417415380477905
jacobi
4.37321936338e-08 0.14313483238220215
gauss_seidel
3.10799718514e-09 0.12006211280822754
----- 2000 -----
gauss
2.46680139423e-12 10.766422271728516
jacobi
2.8758863701e-08 0.2580835819244385
gauss_seidel
4.11025852704e-09 0.18568634986877441

実行環境 : Windows 10, Intel Core i5-6200U 2.30GHz

変数の個数が増えるにしたがい、ガウスの消去法の実行時間は大幅に増加します。反復法でも遅くはなりますが、ガウスの消去法よりはずっと高速です。ヤコビ法とガウス・ザイデル法を比較した場合、変数の個数が少ないときはヤコビ法のほうが速くなりました。Python の実行速度は遅いので、二重の for ループになるガウス・ザイデル法よりも、NumPy で高速に処理できるヤコビ法のほうが速くなるのでしょう。興味のある方はいろいろ試してみてください。

●SOR 法

「SOR 法 (Successive Over-Relaxation)」は、ヤコビ法の漸化式に加速パラメータ w を導入することで収束を速くする方法です。具体的には、ヤコビ法で求めた値 xi+1 と xi の差分に w を掛け算し、その値を xi に足し算します。SOR 法の漸化式は次のようになります。

xi+1 = xi + w * (D-1(b - A'xi) - xi)
=> xi+1 = (1 - w) * xi + w * D-1(b - A'xi)

漸化式は 参考 URL 6 を参考にさせていただきました。また、w が 0 < w < 2 のとき SOR 法は収束するそうです。w が 1 のときはヤコビ法と同じになります。一般に、w を 1 より大きくすると収束が速くなり、1 より小さくすると遅くなるといわれています。ところが 参考 URL 3 によると、1 より小さくしたほうが収束が速くなる場合もあるそうです。これは後で試してみましょう。連立方程式によって w の最適値は変化するでしょうから、w の選択はけっこう難しい問題なのかもしれません。

●SOR 法のプログラム

それではプログラムを作りましょう。次のリストを見てください。

リスト : 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 法の漸化式をそのままプログラムしただけなので、難しいところはないと思います。

●簡単なテスト (2)

それでは実行してみましょう。次に示す方程式を 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
2.23327647078e-08 0.011053323745727539
gauss_seidel
5.27912784798e-09 0.022308349609375
----- 1000 -----
jacobi
1.54129253341e-08 0.0484774112701416
gauss_seidel
6.8278132177e-10 0.06422185897827148
----- 1500 -----
jacobi
2.44906892969e-08 0.10906100273132324
gauss_seidel
1.11494481226e-09 0.12014651298522949
----- 2000 -----
jacobi
3.19235474893e-08 0.1821589469909668
gauss_seidel
1.37436358821e-09 0.1852719783782959

ヤコビ法の場合、w を 0.9 にすると収束回数が 10 回以上減少し、その分だけ実行速度も速くなりました。ガウス・ザイデル法の場合、w を 0.95 にすると収束回数は数回減少しますが、実行速度はほとんど変わりませんでした。今回のテストでは w を 1 より増やすと、どちらの方法でも収束回数は大幅に増加します。w の設定は教科書通りにはいきませんね。収束の判定方法や許容誤差の値を変えると、異なる結果になるかもしれません。興味のある方はいろいろ試してみてください。

●参考文献・URL

  1. 奥村晴彦,『C言語による最新アルゴリズム事典』, 技術評論社, 1991
  2. 連立一次方程式(反復法), (山本昌志さん)
  3. 数値演算法 (8) 連立方程式を解く -2-, (fussy さん)
  4. ヤコビ法 - Wikipedia
  5. ガウス=ザイデル法 - Wikipedia
  6. SOR法 - Wikipedia

●プログラムリスト

#
# jacobi.py : 連立方程式の解法 (反復法)
#
#             Copyright (C) 2018 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)」といい AT と表記します。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 は正方行列を表します。

  1. trace(AT) = trace(A)
  2. |AT| = |A|, (|A| は A の行列式)
  3. (AB)T = BTAT
  4. (AT)-1 = (A-1)T

簡単な例を示しましょう。

>>> 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 とその転置行列 AT が等しい行列のことを「対称行列 (symmetric matrix)」といい、その中で対角成分 aii 以外の要素が 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 の基本的な性質を以下に示します。

  1. |A| は対角成分 aii を乗算した値になる
  2. An は対角成分 aii を n 乗したものになる
  3. A-1 は対角成分 aii を逆数にしたものになる
  4. 対角行列の固有値は対角成分 aii になる

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]])

AT と 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 となるような正方行列を「置換行列」といいます。置換行列を掛け算すると、列または行を入れ替える働きをします。

直交行列の基本的な性質を以下に示します。

  1. 直交行列の行列式は ±1
  2. 直交行列の逆行列も直交行列
  3. 対称行列は直交行列で対角化できる

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 = λx を満たす数 λ を A の「固有値」、ゼロではないベクトル x を A の「固有ベクトル」といいます。このとき、固有値は |A - λI| = 0 の根になります。この式を「固有方程式」とか「特性方程式」といいます。

Ax = λx を変形すると (A - λI)x = 0 になります。行列 (A - λI) に逆行列が存在する場合、両辺にその逆行列を掛け算すると x = 0 になってしまいます。つまり、ゼロではないベクトル x が存在するためには、逆行列が存在しないこと (|A - λ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 - λ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] = 3 * [√2 / 2, √2 / 2]
 [2, 1]]

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] = -1 * [-√2 / 2, √2 / 2]
 [2, 1]]

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 の固有ベクトル x1, x2, ..., xn を列に持つ行列 X = [x1, x2, ...., xn] を考えます。すると、X によって行列 A を以下のように対角行列 Λ に変換することができます。

X-1AX = Λ (ただし Λ は行列 A の固有値 λi を格納した対角行列)

固有値と固有ベクトルの定義 Axi = λixi から AX は次のように変形することができます。

AX = [Ax1, ... Axn]
   = [λ1x1, ..., λnxn]
   = [x1, ..., xn] @ [[λ1, 0, ...], ..., [0, ..., λn]]
   = XΛ

左側から両辺に X の逆行列を掛け算すれば X-1AX = Λ となります。

それでは実際に 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-1AP に変換することを「相似変換」といい、B = P-1AP が成り立つ正方行列 A, B の関係を「相似」といいます。相似な行列 A と B の間は様々な性質が保存されます。保存される主な性質を以下に示します。

つまり、相似変換を行ってもこれらの値が変化することはありません。たとえば、正方行列 A を相似変換により対角行列に変換します。対角行列は対角成分が固有値になるので、トレースは固有値の総和になります。相似変換でトレースは保存されるので、行列 A のトレースも固有値の総和であることがわかります。

●実対称行列

要素がすべて実数の対称行列を「実対称行列」といいます。実対称行列には次に示す性質があります。

  1. 固有値はすべて実数
  2. 二つの異なる固有値に対する固有ベクトルは直交する (内積が 0 になる)
  3. 適当な直交行列 L によって LTAL が対角行列になるように変換できる

3 は直交行列の定義 LT = L-1 から明らかです。2 は固有ベクトルが「一次独立」であることを言っています。ベクトル x1, ..., xi, ..., xn において、式 a1x1 + ... + aixi + ... + anxn = 0 が成り立つとき、係数 ai = 0 (i = 1, ..., n) 以外の解がないことを一次独立といいます。逆に、係数 ai != 0 (i = 1, ..., 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]]))

●参考文献・URL

  1. 奥村晴彦,『C言語による最新アルゴリズム事典』, 技術評論社, 1991
  2. 線形代数, (菅沼さん)
  3. 行列の固有値問題 (PDF), (桂田祐史さん)
  4. 固有値と固有ベクトル, (武内修さん)
  5. 固有値,固有ベクトルの定義と具体的な計算方法 | 高校数学の美しい物語
  6. 固有値 - Wikipedia

Copyright (C) 2018 Makoto Hiroi
All rights reserved.

[ Home | Light | Python3 ]