M.Hiroi's Home Page

Functional Programming

お気楽 Haskell プログラミング入門

[ PrevPage | Haskell | NextPage ]

線形代数編

●複素数

Prelude> import Prelude hiding ((<>))
Prelude> :m + Numeric.LinearAlgebra
Prelude Numeric.LinearAlgebra> a = 2 :+ 1
Prelude Numeric.LinearAlgebra> b = 1 :+ (-1)
Prelude Numeric.LinearAlgebra> a + b
3.0 :+ 0.0
Prelude Numeric.LinearAlgebra> a - b
1.0 :+ 2.0
Prelude Numeric.LinearAlgebra> a * b
3.0 :+ (-1.0)
Prelude Numeric.LinearAlgebra> a / b
0.5 :+ 1.5

Prelude Numeric.LinearAlgebra> polar a
(2.23606797749979,0.4636476090008061)
Prelude Numeric.LinearAlgebra> polar b
(1.4142135623730951,-0.7853981633974483)
Prelude Numeric.LinearAlgebra> mkPolar (fst $ polar a) (snd $ polar a)
2.0 :+ 1.0
Prelude Numeric.LinearAlgebra> mkPolar (fst $ polar b) (snd $ polar b)
1.0000000000000002 :+ (-1.0)
Prelude Numeric.LinearAlgebra> mkPolar 1 (pi/4)
0.7071067811865476 :+ 0.7071067811865475
Prelude Numeric.LinearAlgebra> exp (0 :+ (pi/4))
0.7071067811865476 :+ 0.7071067811865475
Prelude Numeric.LinearAlgebra> exp (0 :+ pi) + 1
0.0 :+ 1.2246467991473532e-16

Prelude Numeric.LinearAlgebra> magnitude a
2.23606797749979
Prelude Numeric.LinearAlgebra> magnitude b
1.4142135623730951
Prelude Numeric.LinearAlgebra> phase a
0.4636476090008061
Prelude Numeric.LinearAlgebra> phase b
-0.7853981633974483
Prelude Numeric.LinearAlgebra> conjugate a
2 :+ (-1)
Prelude Numeric.LinearAlgebra> conjugate b
1 :+ 1
Prelude Numeric.LinearAlgebra> a * conjugate a
5.0 :+ 0.0
Prelude Numeric.LinearAlgebra> b * conjugate b
2.0 :+ 0.0

Prelude Numeric.LinearAlgebra> sqrt ((-1) :+ 0)
0.0 :+ 1.0
Prelude Numeric.LinearAlgebra> (sqrt ((-1) :+ 0)) ^ 2
(-1.0) :+ 0.0
Prelude Numeric.LinearAlgebra> sqrt ((-2) :+ 0)
0.0 :+ 1.4142135623730951
Prelude Numeric.LinearAlgebra> (sqrt ((-2) :+ 0)) ^ 2
(-2.0000000000000004) :+ 0.0

●複素ベクトル

Prelude Numeric.LinearAlgebra> a = fromList [1 :+ 2, 3 :+ 4] :: Vector C
Prelude Numeric.LinearAlgebra> a
[1.0 :+ 2.0,3.0 :+ 4.0]
Prelude Numeric.LinearAlgebra> b = fromList [5 :+ 6, 7 :+ 8] :: Vector C
Prelude Numeric.LinearAlgebra> b
[5.0 :+ 6.0,7.0 :+ 8.0]
Prelude Numeric.LinearAlgebra> a + b
[6.0 :+ 8.0,10.0 :+ 12.0]
Prelude Numeric.LinearAlgebra> a - b
[(-4.0) :+ (-4.0),(-4.0) :+ (-4.0)]
Prelude Numeric.LinearAlgebra> a * b
[(-7.0) :+ 16.0,(-11.0) :+ 52.0]
Prelude Numeric.LinearAlgebra> a / b
[0.2786885245901639 :+ 6.557377049180328e-2,0.4690265486725664 :+ 3.5398230088495575e-2]

Prelude Numeric.LinearAlgebra> c = (1 :+ (-1)) :: C
Prelude Numeric.LinearAlgebra> a + (scalar c)
[2.0 :+ 1.0,4.0 :+ 3.0]
Prelude Numeric.LinearAlgebra> a - (scalar c)
[0.0 :+ 3.0,2.0 :+ 5.0]
Prelude Numeric.LinearAlgebra> a * (scalar c)
[3.0 :+ 1.0,7.0 :+ 1.0]
Prelude Numeric.LinearAlgebra> a / (scalar c)
[(-0.5) :+ 1.5,(-0.5) :+ 3.5]

Prelude Numeric.LinearAlgebra> a <.> b
70.0 :+ (-8.0)
Prelude Numeric.LinearAlgebra> b <.> a
70.0 :+ 8.0
Prelude Numeric.LinearAlgebra> (asRow (cmap conjugate a)) <> (asColumn b)
(1><1)
 [ 70.0 :+ (-8.0) ]
Prelude Numeric.LinearAlgebra> (asRow a) <> (asColumn (cmap conjugate b))
(1><1)
 [ 70.0 :+ 8.0 ]
Prelude Numeric.LinearAlgebra> udot a b
(-18.0) :+ 68.0
Prelude Numeric.LinearAlgebra> (asRow a) <> (asColumn b)
(1><1)
 [ (-18.0) :+ 68.0 ]

Prelude Numeric.LinearAlgebra> a <.> a
30.0 :+ 0.0
Prelude Numeric.LinearAlgebra> b <.> b
174.0 :+ 0.0
Prelude Numeric.LinearAlgebra> norm_2 a
5.477225575051661
Prelude Numeric.LinearAlgebra> norm_2 b
13.19090595827292

●複素行列

Prelude Numeric.LinearAlgebra> a = fromLists [[1 :+ 1, 1 :+ (-1)], [2 :+ 3, 2 :+ 1]] :: Matrix C
Prelude Numeric.LinearAlgebra> a
(2><2)
 [ 1.0 :+ 1.0, 1.0 :+ (-1.0)
 , 2.0 :+ 3.0, 2.0 :+ 1.0 ]
Prelude Numeric.LinearAlgebra> tr a
(2><2)
 [ 1.0 :+ (-1.0), 2.0 :+ (-3.0)
 , 1.0 :+ 1.0,    2.0 :+ (-1.0) ]
Prelude Numeric.LinearAlgebra> (tr $ tr a) == a
True
Prelude Numeric.LinearAlgebra> b = fromLists [[1 :+ 0, 1 :+ 1], [1 :+ (-1), 0 :+ 1]] :: Matrix C
Prelude Numeric.LinearAlgebra> b
(2><2)
 [ 1.0 :+ 0.0,    1.0 :+ 1.0
 , 1.0 :+ (-1.0), 0.0 :+ 1.0 ]
Prelude Numeric.LinearAlgebra> tr (a + b)
(2><2)
 [ 2.0 :+ (-1.0), 3.0 :+ (-2.0)
 , 2.0 :+ (-0.0), 2.0 :+ (-2.0) ]
Prelude Numeric.LinearAlgebra> tr a + tr b
(2><2)
 [ 2.0 :+ (-1.0), 3.0 :+ (-2.0)
 , 2.0 :+ 0.0,    2.0 :+ (-2.0) ]
Prelude Numeric.LinearAlgebra> tr $ (scalar c) * a
(2><2)
 [ 2.0 :+ (-0.0), 5.0 :+ (-1.0)
 , 0.0 :+ 2.0,    3.0 :+ 1.0 ]
Prelude Numeric.LinearAlgebra> (scalar (conjugate c)) * tr a
(2><2)
 [ 2.0 :+ 0.0, 5.0 :+ (-1.0)
 , 0.0 :+ 2.0, 3.0 :+ 1.0 ]
Prelude Numeric.LinearAlgebra> tr (a <> b)
(2><2)
 [ 1.0 :+ 1.0,    5.0 :+ (-2.0)
 , 1.0 :+ (-3.0), (-2.0) :+ (-7.0) ]
Prelude Numeric.LinearAlgebra> tr b <> tr a
(2><2)
 [ 1.0 :+ 1.0,    5.0 :+ (-2.0)
 , 1.0 :+ (-3.0), (-2.0) :+ (-7.0) ]

Prelude Numeric.LinearAlgebra> det a
(-4.0) :+ 2.0000000000000004
Prelude Numeric.LinearAlgebra> det (tr a)
(-4.0) :+ (-2.0)
Prelude Numeric.LinearAlgebra> conjugate $ det a
(-4.0) :+ (-2.0000000000000004)
Prelude Numeric.LinearAlgebra> inv a
(2><2)
 [ (-0.3) :+ (-0.4),            0.3 :+ (-9.999999999999999e-2)
 , 9.999999999999994e-2 :+ 0.8, (-9.999999999999998e-2) :+ (-0.3) ]
Prelude Numeric.LinearAlgebra> tr $ inv a
(2><2)
 [ (-0.3) :+ 0.4,               9.999999999999994e-2 :+ (-0.8)
 , 0.3 :+ 9.999999999999999e-2, (-9.999999999999998e-2) :+ 0.3 ]
Prelude Numeric.LinearAlgebra> inv $ tr a
(2><2)
 [ (-0.29999999999999993) :+ 0.4, 9.999999999999998e-2 :+ (-0.7999999999999999)
 , 0.3 :+ 9.999999999999999e-2,   (-9.999999999999999e-2) :+ 0.3 ]

Prelude Numeric.LinearAlgebra> u = fromList [1 :+ 2, 3 :+ 4] :: Vector C
Prelude Numeric.LinearAlgebra> u
[1.0 :+ 2.0,3.0 :+ 4.0]
Prelude Numeric.LinearAlgebra> v = fromList [5 :+ 6, 7 :+ 8] :: Vector C
Prelude Numeric.LinearAlgebra> v
[5.0 :+ 6.0,7.0 :+ 8.0]
Prelude Numeric.LinearAlgebra> (flatten (a <> (asColumn u))) <.> v
184.0 :+ (-126.0)
Prelude Numeric.LinearAlgebra> u <.> flatten ((tr a) <> (asColumn v))
184.0 :+ (-126.0)

●固有値と固有ベクトル

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

与えられた行列の固有値と固有ベクトルを求める問題のことを「固有値問題」といいます。Haskell (hmatrix) には固有値や固有ベクトルを求める関数があらかじめ用意されていますが、Haskell とアルゴリズムのお勉強ということで、あえてプログラムを作ってみましょう。

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

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

それでは実際に hmatrix で固有値と固有ベクトルを求めてみましょう。hmatrix に用意されている関数 eig, eigevalues を使います。

Prelude> import Prelude hiding ((<>))
Prelude> :m + Numeric.LinearAlgebra
Prelude Numeric.LinearAlgebra> a = matrix 2 [1,2,2,1]
Prelude Numeric.LinearAlgebra> a
(2><2)
 [ 1.0, 2.0
 , 2.0, 1.0 ]
Prelude Numeric.LinearAlgebra> eig a
([3.0000000000000004 :+ 0.0,(-0.9999999999999996) :+ 0.0],(2><2)
 [ 0.7071067811865476 :+ 0.0, (-0.7071067811865474) :+ 0.0
 , 0.7071067811865474 :+ 0.0,    0.7071067811865476 :+ 0.0 ])
Prelude Numeric.LinearAlgebra> eigenvalues a
[3.0000000000000004 :+ 0.0,(-0.9999999999999996) :+ 0.0]

Prelude Numeric.LinearAlgebra> b = diag $ vector [1,2,3]
Prelude Numeric.LinearAlgebra> b
(3><3)
 [ 1.0, 0.0, 0.0
 , 0.0, 2.0, 0.0
 , 0.0, 0.0, 3.0 ]
Prelude Numeric.LinearAlgebra> eig b
([1.0 :+ 0.0,2.0 :+ 0.0,3.0 :+ 0.0],(3><3)
 [ 1.0 :+ 0.0, 0.0 :+ 0.0, 0.0 :+ 0.0
 , 0.0 :+ 0.0, 1.0 :+ 0.0, 0.0 :+ 0.0
 , 0.0 :+ 0.0, 0.0 :+ 0.0, 1.0 :+ 0.0 ])
Prelude Numeric.LinearAlgebra> eigenvalues b
[1.0 :+ 0.0,2.0 :+ 0.0,3.0 :+ 0.0]

Prelude Numeric.LinearAlgebra> c = fromLists [[1 :+ 0, 2 :+ (-3)], [2 :+ 3, 4 :+ 0]] :: Matrix C
Prelude Numeric.LinearAlgebra> c
(2><2)
 [ 1.0 :+ 0.0, 2.0 :+ (-3.0)
 , 2.0 :+ 3.0, 4.0 :+ 0.0 ]
Prelude Numeric.LinearAlgebra> eig c
([(-1.4051248379533272) :+ 0.0,6.405124837953327 :+ 0.0],(2><2)
 [ 0.8318986235710117 :+ 0.0, 0.30781846803228086 :+ (-0.46172770204842134)
 , (-0.30781846803228097) :+ (-0.4617277020484214), 0.8318986235710119 :+ 0.0 ])
Prelude Numeric.LinearAlgebra> eigenvalues c
[(-1.4051248379533272) :+ 0.0,6.405124837953327 :+ 0.0]

eig の返り値はタプルで、第 1 要素が固有値を格納したベクトル、第 2 要素が固有ベクトル (列ベクトル) を格納した行列です。つまり、i 番目の固有値に対応する固有ベクトルは i 列目の列ベクトルになります。eigenvalues は固有値だけを返します。

行列 a の成分は実数で対称行列になっています。これを「実対称行列」といいます。実対称行列はあとで詳しく説明します。行列 b は対角行列です。対角行列の固有値は対角成分と同じになります。行列 c はエルミート行列です。エルミート行列の固有値も実数になります。

eig と eigenvalues は固有値を複素数で返しますが、エルミート行列や実対称行列の固有値は実数になることが知られています。この場合、関数 eigSH や eigenvaluesSH を使うと、固有値を実数で求めることができます。ただし、引数はエルミート行列や実対称行列を表すデータ型 Herm でなければなりません。

eigSH :: Field t => Herm t -> (Vector Double, Matrix t)
eigenvaluesSH :: Field t => Herm t -> Vector Double

関数 sym を使うと、Matrix を Herm に変換することができます。

sym :: Field t => Matrix t -> Herm t

sym は Matrix を次の式で Herm に変換します。

(A + tr A) / 2

行列 A がエルミート行列または実対称行列の場合、tr A は A になるので (A + tr A) / 2 = A になります。つまり、A のデータ型が Herm に変わるだけで A の各成分の値は変化しません。

簡単な実行例を示します。

Prelude Numeric.LinearAlgebra> sym a
Herm (2><2)
 [ 1.0, 2.0
 , 2.0, 1.0 ]
Prelude Numeric.LinearAlgebra> eigSH $ sym a
([3.0,-1.0],(2><2)
 [ 0.7071067811865475, -0.7071067811865475
 , 0.7071067811865475,  0.7071067811865475 ])
Prelude Numeric.LinearAlgebra> eigenvaluesSH $ sym a
[3.0,-1.0]

Prelude Numeric.LinearAlgebra> sym b
Herm (3><3)
 [ 1.0, 0.0, 0.0
 , 0.0, 2.0, 0.0
 , 0.0, 0.0, 3.0 ]
Prelude Numeric.LinearAlgebra> eigSH $ sym b
([3.0,2.0,1.0],(3><3)
 [ 0.0, -0.0, 1.0
 , 0.0,  1.0, 0.0
 , 1.0,  0.0, 0.0 ])
Prelude Numeric.LinearAlgebra> eigenvaluesSH $ sym b
[3.0,2.0,1.0]

Prelude Numeric.LinearAlgebra> sym c
Herm (2><2)
 [ 1.0 :+ 0.0, 2.0 :+ (-3.0)
 , 2.0 :+ 3.0, 4.0 :+ 0.0 ]
Prelude Numeric.LinearAlgebra> eigSH $ sym c
([6.405124837953327,-1.4051248379533274],(2><2)
 [ 0.3078184680322809 :+ (-0.4617277020484213), 0.46145432973433825 :+ (-0.6921814946015072)
 , 0.8318986235710117 :+ 0.0,                 (-0.5549276350125808) :+ 0.0 ])
Prelude Numeric.LinearAlgebra> eigenvaluesSH $ sym c
[6.405124837953327,-1.4051248379533274]

●行列の対角化

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 = Λ となります。

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

Prelude Numeric.LinearAlgebra> a
(2><2)
 [ 1.0, 2.0
 , 2.0, 1.0 ]
Prelude Numeric.LinearAlgebra> (xs, ys) = eigSH $ sym a
Prelude Numeric.LinearAlgebra> xs
[3.0,-1.0]
Prelude Numeric.LinearAlgebra> ys
(2><2)
 [ 0.7071067811865475, -0.7071067811865475
 , 0.7071067811865475,  0.7071067811865475 ]
Prelude Numeric.LinearAlgebra> (inv ys) <> a <> ys
(2><2)
 [ 3.0,  0.0
 , 0.0, -1.0 ]

Prelude Numeric.LinearAlgebra> b = matrix 2 [4,-2,1,1]
Prelude Numeric.LinearAlgebra> b
(2><2)
 [ 4.0, -2.0
 , 1.0,  1.0 ]
Prelude Numeric.LinearAlgebra> (xs, ys) = eig b
Prelude Numeric.LinearAlgebra> xs
[3.0 :+ 0.0,2.0 :+ 0.0]
Prelude Numeric.LinearAlgebra> ys
(2><2)
 [ 0.8944271909999159 :+ 0.0, 0.7071067811865475 :+ 0.0
 , 0.4472135954999579 :+ 0.0, 0.7071067811865475 :+ 0.0 ]
Prelude Numeric.LinearAlgebra> (inv (cmap realPart ys)) <> b <> (cmap realPart ys)
(2><2)
 [ 3.0, 0.0
 , 0.0, 2.0 ]

数学 (線形代数) の世界では、正則行列 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

このように、実対称行列の固有ベクトルは一次独立になります。

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

Prelude Numeric.LinearAlgebra> a
(2><2)
 [ 1.0, 2.0
 , 2.0, 1.0 ]
Prelude Numeric.LinearAlgebra> (xs, ys) = eigSH $ sym a
Prelude Numeric.LinearAlgebra> xs
[3.0,-1.0]
Prelude Numeric.LinearAlgebra> ys
(2><2)
 [ 0.7071067811865475, -0.7071067811865475
 , 0.7071067811865475,  0.7071067811865475 ]

Prelude Numeric.LinearAlgebra> zs = toColumns ys
Prelude Numeric.LinearAlgebra> (zs !! 0) <.> (zs !! 1)
0.0
Prelude Numeric.LinearAlgebra> (tr ys) <> a <> ys
(2><2)
 [ 2.9999999999999996,                 0.0
 ,                0.0, -0.9999999999999998 ]

Prelude Numeric.LinearAlgebra> b = matrix 3 [1,4,5,4,2,6,5,6,3]
Prelude Numeric.LinearAlgebra> b
(3><3)
 [ 1.0, 4.0, 5.0
 , 4.0, 2.0, 6.0
 , 5.0, 6.0, 3.0 ]
Prelude Numeric.LinearAlgebra> (xs, ys) = eigSH $ sym b
Prelude Numeric.LinearAlgebra> xs
[12.175971065046912,-2.5072879670936405,-3.668683097953268]
Prelude Numeric.LinearAlgebra> ys
(3><3)
 [ 0.49659978454619114,   0.8095854617397508, 0.31298567719355963
 ,  0.5773502691896256,  -0.5773502691896257,  0.5773502691896255
 ,  0.6481167492476515, -0.10600965430705465, -0.7541264035547062 ]

Prelude Numeric.LinearAlgebra> zs = toColumns ys
Prelude Numeric.LinearAlgebra> (zs !! 0) <.> (zs !! 1)
1.249000902703301e-16
Prelude Numeric.LinearAlgebra> (zs !! 0) <.> (zs !! 2)
-2.220446049250313e-16
Prelude Numeric.LinearAlgebra> (zs !! 1) <.> (zs !! 2)
1.942890293094024e-16
Prelude Numeric.LinearAlgebra> (inv ys) <> b <> ys
(3><3)
 [      12.175971065046905,   1.609823385706477e-15, -6.661338147750939e-16
 , -1.5543122344752192e-15,      -2.507287967093641, 1.6653345369377348e-16
 ,  1.7763568394002505e-15, -1.3877787807814457e-16,    -3.6686830979532656 ]

Prelude Numeric.LinearAlgebra> c = matrix 4 [5,1,1,1,1,6,1,1,1,1,7,1,1,1,1,8]
Prelude Numeric.LinearAlgebra> c
(4><4)
 [ 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 ]
Prelude Numeric.LinearAlgebra> (xs, ys) = eigSH $ sym c
Prelude Numeric.LinearAlgebra> xs
[9.803886359051248,6.507748705363649,5.3922752902729805,4.296089645312117]
Prelude Numeric.LinearAlgebra> ys
(4><4)
 [ 0.33200196406021854, 0.13594052686577587,   0.2259029659858124,  -0.9056835644828857
 ,  0.4011130835259563, 0.22610178940387854,   0.8017816195421144,  0.38096260921130765
 ,  0.5065613134846454,  0.6714043318139947,  -0.5175355099230781,  0.15738124053003655
 ,  0.6872253093165059, -0.6925419678295478, -0.19562995806342287, 9.917618936878829e-2 ]
Prelude Numeric.LinearAlgebra> (inv ys) <> c <> ys
(4><4)
 [      9.803886359051251,  -2.220446049250313e-15,   9.992007221626409e-16, -3.3306690738754696e-16
 ,  8.881784197001252e-16,       6.507748705363648,  -7.771561172376096e-16,   2.220446049250313e-16
 ,  6.661338147750939e-16, -1.5543122344752192e-15,       5.392275290272984,  -3.191891195797325e-16
 , -5.551115123125783e-16,  3.3306690738754696e-16, -1.8041124150158794e-16,       4.296089645312119 ]

●累乗法

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

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

Prelude Numeric.LinearAlgebra> iter n x = if n == 0 then [x] else let x1 = a <> x in 
x : iter (n - 1) (x1 / (scalar (norm_2 x1)))
Prelude Numeric.LinearAlgebra> iter 10 (col [1, 0])
[(2><1)
 [ 1.0
 , 0.0 ],(2><1)
 [ 0.4472135954999579
 , 0.8944271909999159 ],(2><1)
 [ 0.7808688094430303
 , 0.6246950475544242 ],(2><1)
 [ 0.6804510993672778
 , 0.7327934916262993 ],(2><1)
 [ 0.7157819483772108
 , 0.6983238520753277 ],(2><1)
 [ 0.7041909139949838
 , 0.7100106736147771 ],(2><1)
 [ 0.708076083151601
 , 0.706136148677213 ],(2><1)
 [ 0.7067833845845355
 , 0.7074300299501206 ],(2><1)
 [  0.707214547210912
 , 0.7069989987356877 ],(2><1)
 [ 0.7070708555277229
 , 0.7071427050202062 ],(2><1)
 [  0.707118756000577
 , 0.7070948061697219 ]]
Prelude Numeric.LinearAlgebra> map (\x -> norm_2 (a <> x)) (iter 10 (col [1, 0]))
[2.23606797749979,2.863564212655271,2.9836955314492535,2.998172959635653,2.999796803034175,2.9999774201803375,
 2.9999974911017744,2.9999997212331677,2.9999999690259034,2.999999996558434,2.999999999617604]

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

Rayleigh 商 = (Ax, x) / (x, x)   # (x, y) は x と y の内積を表す

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

●累乗法のプログラム

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

リスト : 累乗法

power :: Matrix R -> Maybe (R, Matrix R, Int)
power xs =
  iter 0 z 0
  where
    z = col (1 : replicate (rows xs - 1) 0)
    iter i x k
      | i >= 512  = Nothing
      | otherwise = let x0 = flatten x
                        x1 = xs <> x
                        n  = norm_2 x1
                        k1 = ((flatten x1) <.> x0) / (x0 <.> x0)
                    in
                      if abs (k1 - k) < 1e-8
                      then Just (k1, x1 / scalar n, i + 1)
                      else iter (i + 1) (x1 / scalar n) k1

関数 power は実対称行列 xs を受け取り、固有値、固有ベクトル (列ベクトル)、繰り返しの回数を返します。実際の処理は局所関数 iter で行います。引数 i が繰り返しの回数、引数 x が固有ベクトル、引数 k が固有値を表します。x は大きさ 1 のベクトルに初期化します。今回は単純に先頭要素を 1 に、残りの要素を 0 に初期化しています。

引数 i が 512 を以上になった場合、収束しなかったので Nothing を返します。あとは繰り返しの中で x1 = xs <> x を計算して Rayleigh 商を求めます。変数 n には列ベクトル x1 の大きさをセットします。abs(k1 - k) が 1e-8 未満になったら収束したとみなして、固有値 k1 と固有ベクトル x1 / n と回数 i を返します。そうでなければ、iter を再帰呼び出しして処理を繰り返します。

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

*Main> power $ matrix 2 [2,1,1,2]
Just (2.999999999426406,(2><1)
 [ 0.7071107728137574
 , 0.7071027895368046 ],11)
*Main> power $ matrix 3 [1,4,5,4,2,6,5,6,3]
Just (12.175971064806813,(3><1)
 [ 0.49659938075279125
 ,  0.5773496196304903
 ,  0.6481176372761912 ],11)

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

●複数の固有値を求める

参考文献 1 によると、行列 A が対称行列であり、固有値 λi (|λ1| > |λ2| > ... > |λn|) がすべて異なる場合、固有ベクトル xi は互いに直交し、次式が成り立つとのことです。

A = λ1x1x1T + λ2x2x2T + ... + λnxnxnT

したがって、λ1 と x1 が求まったならば、行列 (A - λ1x1x1T) に再度累乗法を適用すると、次に大きな固有値 λ2 と固有ベクトル x2 を求めることができます。これを繰り返すことで、n 次の正方行列であれば n 個の固有値と固有ベクトルを求めることができます。

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

リスト : 累乗法 (2)

powerAll :: Matrix R -> Maybe ([R], Matrix R)
powerAll xs =
  iter 1 xs
  where
    iter i xs =
      power xs >>= \(e, v, _) ->
        if i == rows xs
        then
          return ([e], v)
        else
          iter (i + 1) (xs - (scalar e) * v <> (tr v)) >>= \(es, vs) -> return (e:es, v ||| vs)

関数 powerAll は実対称行列 xs の固有値と固有ベクトルを求めます。固有値はリストに格納して返します。実際の処理は局所関数 iter で行います。引数 i が求めた固有値と固有ベクトルの個数、引数 xs が固有値と固有ベクトルを求める行列です。

まず power を呼び出して絶対値最大の固有値 e と固有ベクトル v を求めす。i が rows xs と等しい場合、固有値と固有ベクトルを全部求めたので、return で ([e], v) を返します。そうでなければ iter を再帰呼び出します。このとき、(scalar e) * v <> (tr v) を計算し、xs から引き算します。これで次の固有値と固有ベクトルを求めることができます。あとは iter の返り値 es, vs の先頭に e と v を連結するだけです。

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

*Main> powerAll $ matrix 2 [1,2,2,1]
Just ([2.9999999988528114,-1.0000000003823961],(2><2)
 [ 0.7071027895368046, -0.7070948061697462
 , 0.7071107728137574,  0.7071187560005527 ])

*Main> powerAll $ matrix 2 [2,1,1,2]
Just ([2.999999999426406,1.0000000001911984],(2><2)
 [ 0.7071107728137574, 0.7070948061697185
 , 0.7071027895368046,  -0.70711875600058 ])

*Main> powerAll $ matrix 3 [1,4,5,4,2,6,5,6,3]
Just ([12.175971064806813,-3.668683093726788,-2.507287970031582],(3><3)
 [ 0.49659938075279125, -0.3130174080074894, -0.8095663080185014
 ,  0.5773496196304903, -0.5773240193997118,  0.5773855217274656
 ,  0.6481176372761912,   0.754133329662888, 0.10596391942918502 ])

*Main> powerAll $ matrix 4 [5,1,1,1,1,6,1,1,1,1,7,1,1,1,1,8]
Just ([9.803886355452113,6.507748686722966,5.392275293049927,4.296089656982493],(4><4)
 [ 0.33200504131765435, 0.13595526259756915,  0.22596747022482527,    0.9056574406615683
 ,  0.4011183591424394, 0.22617975157313117,   0.8017151890754934,  -0.38105529497820484
 ,  0.5065758276117462,  0.6713287148299257,   -0.517642547301511,    -0.157321415017894
 ,  0.6872100445846333, -0.6925869210396215, -0.19554449935766896, -9.915359156469715e-2 ])

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

●逆累乗法

Ax = λx の両辺に左側から A-1 を掛け算すると、x = λA-1x => A-1x = (1/λ)x となります。A-1 に累乗法を適用すると 1 / λ の最大値、つまり最小の固有値とその固有ベクトルを求めることができます。これを「逆累乗法」とか「逆べき乗法」といいます。具体的には、xi+1 = A-1xi のように A-1 を掛け算していくことになりますが、参考 URL 2 によると、逆行列を求めるよりも連立方程式 Axi+1 = xi を解くほうが良いようです。この場合、行列を LU 分解しておいたほうが効率的です。

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

リスト : 逆累乗法

powerInv :: Matrix R -> Maybe (R, Matrix R, Int)
powerInv xs =
  iter 0 z 0
  where
    tri = makeLU' xs
    z = col (1 : replicate (rows xs - 1) 0)
    iter i x k
      | i >= 512  = Nothing
      | otherwise = let x0 = flatten x
                        x1 = col $ luSolver' tri x0
                        n  = norm_2 x1
                        k1 = ((flatten x1) <.> x0) / (x0 <.> x0)
                    in
                      if abs (k1 - k) < 1e-8
                      then Just (1 / k1, x1 / scalar n, i + 1)
                      else iter (i + 1) (x1 / scalar n) k1

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

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

*Main> powerInv $ matrix 2 [1,2,2,1]
Just (-1.0000000003823961,(2><1)
 [ -0.7071027895368047
 ,  0.7071107728137576 ],11)
*Main> powerAll $ matrix 2 [1,2,2,1]
Just ([2.9999999988528114,-1.0000000003823961],(2><2)
 [ 0.7071027895368046, -0.7070948061697462
 , 0.7071107728137574,  0.7071187560005527 ])

*Main> powerInv $ matrix 2 [2,1,1,2]
Just (1.000000000191198,(2><1)
 [  0.7071107728137573
 , -0.7071027895368047 ],11)
*Main> powerAll $ matrix 2 [2,1,1,2]
Just ([2.999999999426406,1.0000000001911984],(2><2)
 [ 0.7071107728137574, 0.7070948061697185
 , 0.7071027895368046,  -0.70711875600058 ])

*Main> powerInv $ matrix 3 [1,4,5,4,2,6,5,6,3]
Just (-2.5072879960642904,(3><1)
 [ -0.8096263206940646
 ,  0.5772748809066349
 , 0.10610811803826374 ],21)
*Main> powerAll $ matrix 3 [1,4,5,4,2,6,5,6,3]
Just ([12.175971064806813,-3.668683093726788,-2.507287970031582],(3><3)
 [ 0.49659938075279125, -0.3130174080074894, -0.8095663080185014
 ,  0.5773496196304903, -0.5773240193997118,  0.5773855217274656
 ,  0.6481176372761912,   0.754133329662888, 0.10596391942918502 ])

*Main> powerInv $ matrix 4 [5,1,1,1,1,6,1,1,1,1,7,1,1,1,1,8]
Just (4.296089899617464,(4><1)
 [    0.9057807814119976
 ,  -0.38061757467977975
 ,  -0.15760282339290624
 , -9.926121059173881e-2 ],28)
*Main> powerAll $ matrix 4 [5,1,1,1,1,6,1,1,1,1,7,1,1,1,1,8]
Just ([9.803886355452113,6.507748686722966,5.392275293049927,4.296089656982493],(4><4)
 [ 0.33200504131765435, 0.13595526259756915,  0.22596747022482527,    0.9056574406615683
 ,  0.4011183591424394, 0.22617975157313117,   0.8017151890754934,  -0.38105529497820484
 ,  0.5065758276117462,  0.6713287148299257,   -0.517642547301511,    -0.157321415017894
 ,  0.6872100445846333, -0.6925869210396215, -0.19554449935766896, -9.915359156469715e-2 ])

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

●参考文献・URL

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

Copyright (C) 2021 Makoto Hiroi
All rights reserved.

[ PrevPage | Haskell | NextPage ]