M.Hiroi's Home Page

Julia Language Programming

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


Copyright (C) 2019-2021 Makoto Hiroi
All rights reserved.

線形代数編

●固有値と固有ベクトル

今回は「固有値 (eigenvalue)」と「固有ベクトル (eigenvector)」について取り上げます。固有値と固有ベクトルは行列の性質を表す重要な指標のひとつで、行列を対角行列に変換する「対角化」の基礎になります。対角行列は単純な構造をしているので、いろいろな計算が簡単になるという利点があります。たとえば、対角行列 \(A\) の逆行列は対角成分を逆数にしたものですし、\(A\) の n 乗は対角成分を n 乗したものになります。また、\(A\) の行列式も対角成分を掛け算するだけで求めることができます。

与えられた行列の固有値と固有ベクトルを求める問題のことを「固有値問題」といいます。Julia には固有値や固有ベクトルを求める関数 eigen(), eigevals(), eigevecs() などがモジュール LinearAlgebra に用意されていますが、Julia とアルゴリズムのお勉強ということで、あえてプログラムを作ってみましょう。

●固有値と固有ベクトルの定義

最初に固有値と固有ベクトルの基本を簡単に説明します。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] = 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 よりも少なくなります。

それでは実際に Julia で固有値と固有ベクトルを求めてみましょう。LinearAlgebra に用意されている関数 eigen(), eigevals(), eigevecs() を使います。 

julia> using LinearAlgebra

julia> f = eigen([1 2; 2 1])
Eigen{Float64,Float64,Matrix{Float64},Vector{Float64}}
eigenvalues:
2-element Vector{Float64}:
 -1.0
  3.0
eigenvectors:
2×2 Matrix{Float64}:
 -0.707107  0.707107
  0.707107  0.707107

julia> f.values
2-element Vector{Float64}:
 -1.0
  3.0

julia> f.vectors
2×2 Matrix{Float64}:
 -0.707107  0.707107
  0.707107  0.707107

julia> val, vec = eigen([1 2; 2 1])
Eigen{Float64,Float64,Matrix{Float64},Vector{Float64}}
eigenvalues:
2-element Vector{Float64}:
 -1.0
  3.0
eigenvectors:
2×2 Matrix{Float64}:
 -0.707107  0.707107
  0.707107  0.707107

julia> val
2-element Vector{Float64}:
 -1.0
  3.0

julia> vec
2×2 Matrix{Float64}:
 -0.707107  0.707107
  0.707107  0.707107

julia> eigvals([1 2; 2 1])
2-element Vector{Float64}:
 -1.0
  3.0

julia> eigvecs([1 2; 2 1])
2×2 Matrix{Float64}:
 -0.707107  0.707107
  0.707107  0.707107

julia> a = diagm(0 => [1,2,3])
3×3 Matrix{Int64}:
 1  0  0
 0  2  0
 0  0  3

julia> eigen(a)
Eigen{Float64,Float64,Matrix{Float64},Vector{Float64}}
eigenvalues:
3-element Vector{Float64}:
 1.0
 2.0
 3.0
eigenvectors:
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

julia> eigvals(a)
3-element Vector{Float64}:
 1.0
 2.0
 3.0

julia> eigvecs(a)
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

eigen() の返り値は固有値と固有ベクトルを格納したデータ型 Eigen を返します。固有値と固有ベクトルはメンバ変数 values と vectors にセットされています。また、val, vec = eigen(...) のようにイテレータで固有値と固有ベクトルを取り出すこともできます。

eigvals() は固有値を、eigvecs() は固有ベクトルを返します。固有値を xs、固有ベクトルを ys とすると、固有値 xs[i] に対応する固有ベクトルは ys[:, i] になります。つまり、ys は固有ベクトルを列に持つ行列になります。そして、対角行列の固有値は対角成分の値となります。

●行列の対角化

n 次の正方行列 A の固有ベクトル \(x_1, x_2, \ldots, x_n\) を列に持つ行列 \(X = [x_1, x_2, \ldots, x_n]\) を考えます。すると、\(X\) によって行列 \(A\) を以下のように対角行列 Λ に変換することができます。

\(X^{-1}AX = \varLambda\) (ただし \(\varLambda\) は行列 \(A\) の固有値 \(\lambda_i\) を格納した対角行列)

固有値と固有ベクトルの定義 \(Ax_i = \lambda_i x_i\) から \(AX\) は次のように変形することができます。

\(\begin{eqnarray} AX &=& (Ax_1, \ \ldots, \ Ax_n) \\ &=& (\lambda_i x_1, \ \ldots, \ \lambda_n x_n) \\ &=& (x_1, \ \ldots, \ x_n) \begin{pmatrix} \lambda_1 & & & & \\ & \lambda_2 & & \LARGE{0} & \\ & & \ddots & & \\ & \LARGE{0} & & \lambda_{n-1} & \\ & & & & \lambda_n \end{pmatrix} \\ &=& X\varLambda \end{eqnarray}\)

左側から両辺に \(X\) の逆行列を掛け算すれば \(X^{-1}AX = \varLambda\) となります。

それでは実際に Julia で試してみましょう。 

julia> a = [1 2; 2 1]
2×2 Matrix{Int64}:
 1  2
 2  1

julia> xs, ys = eigen(a)
Eigen{Float64,Float64,Matrix{Float64},Vector{Float64}}
eigenvalues:
2-element Vector{Float64}:
 -1.0
  3.0
eigenvectors:
2×2 Matrix{Float64}:
 -0.707107  0.707107
  0.707107  0.707107

julia> xs
2-element Vector{Float64}:
 -1.0
  3.0

julia> ys
2×2 Matrix{Float64}:
 -0.707107  0.707107
  0.707107  0.707107

julia> inv(ys) * a * ys
2×2 Matrix{Float64}:
 -1.0  0.0
  0.0  3.0

julia> b = [4 -2; 1 1]
2×2 Matrix{Int64}:
 4  -2
 1   1

julia> xs, ys = eigen(b)
Eigen{Float64,Float64,Matrix{Float64},Vector{Float64}}
values:
2-element Vector{Float64}:
 2.0
 3.0
vectors:
2×2 Matrix{Float64}:
 0.707107  0.894427
 0.707107  0.447214

julia> inv(ys) * b * ys
2×2 Matrix{Float64}:
 2.0  0.0
 0.0  3.0

数学 (線形代数) の世界では、正則行列 \(P\) を用いて正方行列 \(A\) を \(P^{-1}AP\) に変換することを「相似変換」といい、\(B = P^{-1}AP\) が成り立つ正方行列 \(A, B\) の関係を「相似」といいます。相似な行列 \(A\) と \(B\) の間は様々な性質が保存されます。保存される主な性質を以下に示します。

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

●実対称行列

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

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

3 は直交行列の定義 \(L^{\mathrm{T}} = L^{-1}\) から明らかです。2 は固有ベクトルが「一次独立」であることを言っています。

ベクトル \(x_1, \ldots, x_i, \ldots, x_n\) において、式 \(a_1 x_1 + \ldots + a_i x_i + \ldots + 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

このように、実対称行列の固有ベクトルは一次独立になります。それでは実際に Julia で試してみましょう。 

julia> a = [1 2; 2 1]
2×2 Matrix{Int64}:
 1  2
 2  1

julia> xs, ys = eigen(a)
Eigen{Float64,Float64,Matrix{Float64},Vector{Float64}}
eigenvalues:
2-element Vector{Float64}:
 -1.0
  3.0
eigenvectors:
2×2 Matrix{Float64}:
 -0.707107  0.707107
  0.707107  0.707107

julia> dot(ys[:, 1], ys[:, 2])
0.0

julia> ys' * a * ys
2×2 Matrix{Float64}:
 -1.0  0.0
  0.0  3.0

julia> b = [1 4 5; 4 2 6; 5 6 3]
3×3 Matrix{Int64}:
 1  4  5
 4  2  6
 5  6  3

julia> xs, ys = eigen(b)
Eigen{Float64,Float64,Matrix{Float64},Vector{Float64}}
eigenvalues:
3-element Vector{Float64}:
 -3.668683097953268
 -2.5072879670936397
 12.175971065046879
eigenvectors:
3×3 Matrix{Float64}:
 -0.312986   0.809585  0.4966
 -0.57735   -0.57735   0.57735
  0.754126  -0.10601   0.648117

julia> dot(ys[:, 1], ys[:, 2])
-9.71445146547012e-17

julia> dot(ys[:, 1], ys[:, 3])
1.1102230246251565e-16

julia> dot(ys[:, 2], ys[:, 3])
-1.3877787807814457e-17

julia> ys' * b * ys
3×3 Matrix{Float64}:
 -3.66868       4.44089e-16  -4.44089e-16
 -2.77556e-16  -2.50729       7.49401e-16
  8.88178e-16   7.77156e-16  12.176

julia> c = ones(5, 5)
5×5 Matrix{Float64}:
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0

julia> c += diagm(0 => [6,7,8,9,10])
5×5 Matrix{Float64}:
 7.0  1.0  1.0   1.0   1.0
 1.0  8.0  1.0   1.0   1.0
 1.0  1.0  9.0   1.0   1.0
 1.0  1.0  1.0  10.0   1.0
 1.0  1.0  1.0   1.0  11.0

julia> xs, ys = eigen(c)
Eigen{Float64,Float64,Matrix{Float64},Vector{Float64}}
eigenvalues:
5-element Vector{Float64}:
  6.277695819922925
  7.356631854844213
  8.43473666649578
  9.540394425688122
 13.390541233048946
eigenvectors:
5×5 Matrix{Float64}:
 -0.916785    0.21894   -0.133372  -0.0951368  0.291088
  0.352465    0.83285   -0.226332  -0.132586   0.336637
  0.147818   -0.461665  -0.746949  -0.218659   0.399087
  0.093519   -0.180739   0.574469  -0.623289   0.489984
  0.0683951  -0.112364   0.207458   0.73285    0.634499

julia> ys' * c * ys
5×5 Matrix{Float64}:
  6.2777       -1.27259e-15   2.63549e-15   8.90413e-17  -4.41384e-16
 -1.24455e-15   7.35663      -1.96836e-15  -2.699e-16    -6.86134e-16
  2.91423e-15  -1.65208e-15   8.43474      -7.86832e-16   7.34279e-16
  6.19199e-17  -3.10929e-16  -9.03661e-16   9.54039       2.2058e-15
 -6.33799e-16  -2.42131e-16   3.86327e-16   7.08199e-16  13.3905

●累乗法

行列の固有値と固有ベクトルを求める簡単な方法に「累乗法 (power method)」があります。「べき乗法」と呼ばれることもあります。参考文献『Ｃ言語による最新アルゴリズム事典』によると、『適当なベクトル \(x\) を選び、これに行列 \(A\) を何度も掛け算すると、\(x\) に含まれる \(A\) の絶対値最大の固有値に対応する固有ベクトルが最も速く成長する』とのことです。\(x\) が単位ベクトルとすると、\(x\) が収束したときの \(Ax\) の大きさが固有値の絶対値 \(|\lambda|\) となります。

それでは実際に試してみましょう。 

julia> a
2×2 Matrix{Int64}:
 1  2
 2  1

julia> x = [1, 0]
2-element Vector{Int64}:
 1
 0

julia> for i = 1 : 16
       x1 = a * x
       k = norm(x1)
       println("$i, $x, $k")
       global x = x1 / k
       end
1, [1, 0], 2.23606797749979
2, [0.447214, 0.894427], 2.863564212655271
3, [0.780869, 0.624695], 2.9836955314492535
4, [0.680451, 0.732793], 2.998172959635653
5, [0.715782, 0.698324], 2.999796803034175
6, [0.704191, 0.710011], 2.9999774201803375
7, [0.708076, 0.706136], 2.999997491101774
8, [0.706783, 0.70743], 2.9999997212331677
9, [0.707215, 0.706999], 2.999999969025903
10, [0.707071, 0.707143], 2.999999996558434
11, [0.707119, 0.707095], 2.999999999617604
12, [0.707103, 0.707111], 2.9999999999575113
13, [0.707108, 0.707105], 2.999999999995279
14, [0.707106, 0.707107], 2.9999999999994755
15, [0.707107, 0.707107], 2.999999999999942
16, [0.707107, 0.707107], 2.9999999999999933

julia> x
2-element Vector{Float64}:
 0.7071067976130431
 0.7071067647600516

julia> norm(a * x)
2.999999999999999

このように\(x_{i+1} = Ax_i\) を繰り返すことにより、固有値と固有ベクトルの値は理論値に近づいていくことがわかります。なお、固有ベクトルは定数倍してもいいので、-1 倍した固有ベクトルが求まる場合もあります。また、実際に固有値を求めるときは、以下に示す「Rayleigh 商」を用いるほうが良いようです。

Rayleigh 商 \(= \dfrac{(Ax, x)}{(x, x)}\)

\((x, y)\) は \(x\) と \(y\) の内積を表す

固有値の定義により \((Ax, x) = (\lambda x, x)\) になります。内積は結合法則 \((kx, y) = (x, ky) = k(x, y)\) が成り立つので、\(\frac{(Ax, x)}{(x, x)} = \lambda \frac{(x, x)}{(x, x)} = \lambda\) となり、固有値 \(\lambda\) を求めることができます。あとは、\(|\lambda_i - \lambda_{i-1}|\) が許容誤差 \(\epsilon\) の中に入るまで繰り返し \(A\) を掛け算すればいいわけです。

●累乗法のプログラム

Julia を使うと、累乗法はとても簡単にプログラムすることができます。次のリストを見てください。 

リスト : 累乗法

function power(a::Matrix{Float64}; max_iter = 512, eps = 1e-8, tiny = 1e-10)
    x = zeros(size(a, 1))
    _, i = findmax(abs.(diag(a)))
    x[i] = 1.0
    k = 0.0
    for i = 1 : max_iter
        # println("$i, $k, $x")
        x1 = a * x
        n = norm(x1)
        if n < tiny
            error("power($a)\nベクトルの大きさが $n になりました")
        end
        k1 = dot(x1, x) / dot(x, x)
        if abs(k1 - k) < eps
            return k1, x1 / n
        end
        k = k1
        x = x1 / n
    end
end

引数 max_iter が繰り返しの最大回数、eps が許容誤差です。変数 x が固有ベクトル、変数 k が固有値を表します。x は大きさ 1 のベクトルに初期化します。今回は対角成分の一番大きな要素の位置 i を求め、x[i] を 1.0 にしています。単純に x[1] = 1.0 としてもかまいません。

あとはループの中で x1 = a * x を計算して Rayleigh 商を求めます。変数 n にはベクトル x1 の大きさをセットします。n が引数 tiny より小さくなったら、0 とみなしてエラーを送出することにします。abs(k1 - k) が eps 未満になったら収束したとみなして、固有値 k1 と固有ベクトル x1 / n を返します。そうでなければ、k を k1 に x を x1 / n に更新して処理を繰り返します。

それでは実際に試してみましょう。ご参考までに、k と x の値を表示しています。 

julia> power([2. 1.; 1. 2.])
1, 0.0, [1.0, 0.0]
2, 2.0, [0.894427, 0.447214]
3, 2.8000000000000003, [0.780869, 0.624695]
4, 2.9756097560975614, [0.732793, 0.680451]
5, 2.997260273972603, [0.715782, 0.698324]
6, 2.9996952148735145, [0.710011, 0.704191]
7, 2.999966130397968, [0.708076, 0.706136]
8, 2.9999962366542348, [0.70743, 0.706783]
9, 2.999999581849771, [0.707215, 0.706999]
10, 2.999999953538855, [0.707143, 0.707071]
11, 2.9999999948376512, [0.707119, 0.707095]
(2.999999999426406, [0.707111, 0.707103])

julia> power([1. 4. 5.; 4. 2. 6; 5. 6. 3.])
1, 0.0, [0.0, 0.0, 1.0]
2, 3.0, [0.597614, 0.717137, 0.358569]
3, 10.428571428571429, [0.455378, 0.517475, 0.724465]
4, 12.000428449014567, [0.50746, 0.5946, 0.623646]
5, 12.159906346157214, [0.493337, 0.571957, 0.65535]
6, 12.174512735519642, [0.497551, 0.578981, 0.645929]
7, 12.17583872124317, [0.496317, 0.576855, 0.648774]
8, 12.175959052754244, [0.496684, 0.5775, 0.647919]
9, 12.17596997462294, [0.496575, 0.577305, 0.648176]
10, 12.175970966057655, [0.496607, 0.577364, 0.648099]
11, 12.175971056060387, [0.496598, 0.577346, 0.648122]
(12.175971064231074, [0.4966, 0.577352, 0.648115])

このように、累乗法を使うと絶対値で最大の固有値と固有ベクトルを簡単に求めることができます。ただし、累乗法には収束が遅いという欠点があります。絶対値最大の固有値と 2 番目の固有値を λ1 と λ2 とすると、|λ2 / λ1 | が 1 に近い値だと累乗法の収束は極端に遅くなります。ご注意くださいませ。

●複数の固有値を求める

参考文献『Ｃ言語による最新アルゴリズム事典』によると、行列 \(A\) が対称行列であり、固有値 \(\lambda_i \ (|\lambda_1| \gt |\lambda_2| \gt \ldots \gt |\lambda_n|)\) がすべて異なる場合、固有ベクトル \(x_i\) は互いに直交し、次式が成り立つとのことです。

\( A = \lambda_1 x_1 {x_1}^{T} + \lambda_2 x_2 {x_2}^{T} + \ldots + \lambda_n x_n {x_n}^{T} \)

したがって、\(\lambda_1\) と \(x_1\) が求まったならば、行列 \((A - \lambda_1 x_1 {x_1}^{T})\) に再度累乗法を適用すると、次に大きな固有値 \(\lambda_2\) と固有ベクトル \(x_2\) を求めることができます。これを繰り返すことで、n 次の正方行列であれば n 個の固有値と固有ベクトルを求めることができます。

プログラムは次のようになります。 

リスト : 累乗法 (2)

function power_all(a::Matrix{Float64}; max_iter = 512, eps = 1e-8, tiny = 1e-10)
    n = size(a, 1)
    xs = zeros(n)       # 固有値
    yss = zeros(n, n)   # 固有ベクトル
    for i = 1 : n
        x, ys = power(a, max_iter=max_iter, eps=eps, tiny=tiny)
        xs[i] = x
        yss[:, i] = ys
        zs = reshape(ys, n, 1)
        a -= x * zs * zs'
    end
    xs, yss
end

変数 xs は固有値を格納する配列、yss は固有ベクトルを格納する配列です。power() で絶対値最大の固有値 x と固有ベクトル ys を求め、xs と yss にセットします。ys は列ベクトルとして yss にセットしていることに注意してください。そして、reshape() で ys を列ベクトルに変換して変数 zs にセットし、行列 a から x * zs * zsT を引き算します。これで次の固有値と固有ベクトルを求めることができます。

それでは実際に試してみましょう。 

julia> power_all([1. 2.; 2. 1.])
([3.0, -1.0], [0.707103 0.707095; 0.707111 -0.707119])

julia> power_all([2. 1.; 1. 2.])
([3.0, 1.0], [0.707111 -0.707095; 0.707103 0.707119])


julia> power_all([1. 4. 5.; 4. 2. 6; 5. 6. 3])
([12.176, -3.66868, -2.50729],
 [0.4966 0.313028 -0.809564; 0.577352 0.577327 0.57739; 0.648115 -0.754127 0.105957])

julia> c = ones(5, 5)
5×5 Matrix{Float64}:
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0

julia> c += diagm(0 => [6,7,8,9,10])
5×5 Matrix{Float64}:
 7.0  1.0  1.0   1.0   1.0
 1.0  8.0  1.0   1.0   1.0
 1.0  1.0  9.0   1.0   1.0
 1.0  1.0  1.0  10.0   1.0
 1.0  1.0  1.0   1.0  11.0

julia> xs, ys = power_all(c)
([13.3905, 9.54039, 8.43474, 7.35663, 6.2777], 
[0.291085 0.0951708 … 0.218843 0.916818; 0.336633 0.132637 … 0.832933 -0.352339; … ;
 0.489964 0.623219 … -0.180824 -0.0935466; 0.634522 -0.732855 … -0.112391 -0.0684122])

julia> xs
5-element Vector{Float64}:
 13.390541225519467
  9.540394395454161
  8.434736667230538
  7.356631856639677
  6.277695841268191

julia> ys
5×5 Matrix{Float64}:
 0.291085   0.0951708   0.133323   0.218843   0.916818
 0.336633   0.132637    0.226185   0.832933  -0.352339
 0.39908    0.218796    0.746978  -0.46152   -0.147888
 0.489964   0.623219   -0.574554  -0.180824  -0.0935466
 0.634522  -0.732855   -0.20731   -0.112391  -0.0684122

正常に動作しているようです。興味のある方はいろいろ試してみてください。

●逆累乗法

\(Ax = \lambda x\) の両辺に左側から \(A^{-1}\) を掛け算すると、\(x = \lambda A^{-1}x\) => \(A^{-1}x = \frac{1}{\lambda}x\) となります。\(A^{-1}\) に累乗法を適用すると \(\frac{1}{\lambda}\) の最大値、つまり最小の固有値とその固有ベクトルを求めることができます。これを「逆累乗法」とか「逆べき乗法」といいます。

具体的には、\(x_{i+1} = A^{-1} x_i\) のように \(A^{-1}\) を掛け算していくことになりますが、参考 URL 2 『行列の固有値問題』によると、逆行列を求めるよりも連立方程式 \(Ax_{i+1} = x_i\) を解くほうが良いようです。この場合、行列を LU 分解しておいたほうが効率的です。

プログラムは次のようになります。 

リスト : 逆累乗法

function power_inv(a::Matrix{Float64}; max_iter = 512, eps = 1e-8, tiny = 1e-10)
    a1, idx = mylu_pv(a)
    x = zeros(size(a, 1))
    x[1] = 1.0
    k = 0.0
    for i = 1 : max_iter
        x1 = mylu_solver(a1, x[idx])
        n = norm(x1)
        if n < tiny
            error("power_inv($a)\nベクトルの大きさが $n になりました")
        end
        k1 = dot(x1, x) / dot(x, x)
        if abs(k1 - k) < eps
            return 1 / k1, x1 / n
        end
        k = k1
        x = x1 / n
        # println("$i $k $x")
    end
end

最初に関数 mylu_pv() で行列 a を LU 分解します。この関数は拙作のページ「ピボット選択を使った LU 分解」で作成したものです。ループの中で x1 = a-1 * x を計算する代わりに mylu_solver() を呼び出します。この関数も拙作のページで作成したものです。収束判定は累乗法のプログラムと同じです。収束したら固有値 1 / k1 と固有ベクトル x1 / n を返します。

それでは実行してみましょう。 

julia> power_inv([1. 2.;2. 1.])
(-1.0000000003823963, [-0.707103, 0.707111])

julia> power_all([1. 2.;2. 1.])
([3.0, -1.0], [0.707103 0.707095; 0.707111 -0.707119])

julia> power_inv([2. 1.;1. 2.])
(1.0000000001911982, [0.707111, -0.707103])

julia> power_all([2. 1.;1. 2.])
([3.0, 1.0], [0.707111 -0.707095; 0.707103 0.707119])

julia> power_inv([1. 4. 5.; 4. 2. 6.; 5. 6. 3])
(-2.50728799606429, [-0.809626, 0.577275, 0.106108])

julia> power_all([1. 4. 5.; 4. 2. 6.; 5. 6. 3])
([12.176, -3.66868, -2.50729], 
 [0.4966 0.313028 -0.809564; 0.577352 0.577327 0.57739; 0.648115 -0.754127 0.105957])

julia> c = ones(4, 4)
4×4 Matrix{Float64}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0

julia> c += diagm(0 => [4,5,6,7])
4×4 Matrix{Float64}:
 5.0  1.0  1.0  1.0
 1.0  6.0  1.0  1.0
 1.0  1.0  7.0  1.0
 1.0  1.0  1.0  8.0

julia> power_inv(c)
(4.296089899617464, [0.905781, -0.380618, -0.157603, -0.0992612])

julia> power_all(c)
([9.80389, 6.50775, 5.39228, 4.29609], 
 [0.331999 0.135927 0.225845 0.905707; 
  0.401108 0.226026 0.801844 -0.380878; 
  0.506545 0.67148 -0.517431 -0.157436; 
  0.687242 -0.692495 -0.195715 -0.0991969])

正常に動作していますね。

●参考文献・URL

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

●プログラムリスト

リスト : 累乗法と逆累乗法

using LinearAlgebra

# LU 分解
function mylu(xs::Matrix{Float64})
    n = size(xs, 1)
    zs = copy(xs)
    for i = 1 : n
        for j = i + 1 : n
            temp = zs[j, i] / zs[i, i]
            for k = i + 1 : n
                zs[j, k] -= temp * zs[i, k]
            end
            zs[j, i] = temp
        end
    end
    zs
end

# LU 分解による連立方程式の解法
function mylu_solver(xs::Matrix{Float64}, ys::Vector{Float64})
    n = size(xs, 1)
    zs = copy(ys)
    # 前進代入
    for i = 2 : n, j = 1 : i - 1
        zs[i] -= xs[i, j] * zs[j]
    end
    # 後退代入
    for i = n : -1 : 1
       for j = i + 1 : n
           zs[i] -= xs[i, j] * zs[j]
       end
       zs[i] /= xs[i, i]
    end
    zs
end

# ピボット選択
function select_pivot(xs::Matrix{Float64}, idx::Vector{Int}, i::Int)
    _, k = findmax(abs.(xs[i:end, i]))
    k += i - 1
    if k != i
        temp = xs[i, :]
        xs[i, :] = xs[k, :]
        xs[k, :] = temp
        idx[i], idx[k] = idx[k], idx[i]
    end
end

function mylu_pv(xs::Matrix{Float64})
    n = size(xs, 1)
    zs = copy(xs)
    idx = collect(1 : n)
    for i = 1 : n
        select_pivot(zs, idx, i)
        for j = i + 1 : n
            temp = zs[j, i] / zs[i, i]
            for k = i + 1 : n
                zs[j, k] -= temp * zs[i, k]
            end
            zs[j, i] = temp
        end
    end
    zs, idx
end

function mylu_solver_pv(xs::Matrix{Float64}, idx::Vector{Int}, ys::Vector{Float64})
    mylu_solver(xs, ys[idx])
end

# 累乗法
function power(a::Matrix{Float64}; max_iter = 512, eps = 1e-8, tiny = 1e-10)
    x = zeros(size(a, 1))
    _, i = findmax(abs.(diag(a)))
    x[i] = 1.0
    k = 0.0
    for i = 1 : max_iter
        # println("$i, $k, $x")
        x1 = a * x
        n = norm(x1)
        if n < tiny
            error("power($a)\nベクトルの大きさが $n になりました")
        end
        k1 = dot(x1, x) / dot(x, x)
        if abs(k1 - k) < eps
            return k1, x1 / n
        end
        k = k1
        x = x1 / n
    end
end

# 累乗法 (2)
function power_all(a::Matrix{Float64}; max_iter = 512, eps = 1e-8, tiny = 1e-10)
    n = size(a, 1)
    xs = zeros(n)       # 固有値
    yss = zeros(n, n)   # 固有ベクトル
    for i = 1 : n
        x, ys = power(a, max_iter=max_iter, eps=eps, tiny=tiny)
        xs[i] = x
        yss[:, i] = ys
        zs = reshape(ys, n, 1)
        a -= x * zs * zs'
    end
    xs, yss
end

# 逆累乗法
function power_inv(a::Matrix{Float64}; max_iter = 512, eps = 1e-8, tiny = 1e-10)
    a1, idx = mylu_pv(a)
    x = zeros(size(a, 1))
    x[1] = 1.0
    k = 0.0
    for i = 1 : max_iter
        x1 = mylu_solver(a1, x[idx])
        n = norm(x1)
        if n < tiny
            error("power_inv($a)\nベクトルの大きさが $n になりました")
        end
        k1 = dot(x1, x) / dot(x, x)
        if abs(k1 - k) < eps
            return 1 / k1, x1 / n
        end
        k = k1
        x = x1 / n
        # println("$i $k $x")
    end
end
初版 2019 年 1 月 12 日
改訂 2021 年 12 月 5 日