data Complex a = !a +: !a
mkPolar :: Floating a => a -> a -> Complex a
polar :: RealFloat a => Complex a -> (a, a)
ghci> import Prelude hiding ((<>)) ghci> :m + Numeric.LinearAlgebra ghci> a = 2 :+ 1 ghci> b = 1 :+ (-1) ghci> a + b 3.0 :+ 0.0 ghci> a - b 1.0 :+ 2.0 ghci> a * b 3.0 :+ (-1.0) ghci> a / b 0.5 :+ 1.5 ghci> polar a (2.23606797749979,0.4636476090008061) ghci> polar b (1.4142135623730951,-0.7853981633974483) ghci> mkPolar (fst $ polar a) (snd $ polar a) 2.0 :+ 1.0 ghci> mkPolar (fst $ polar b) (snd $ polar b) 1.0000000000000002 :+ (-1.0) ghci> mkPolar 1 (pi/4) 0.7071067811865476 :+ 0.7071067811865475 ghci> exp (0 :+ (pi/4)) 0.7071067811865476 :+ 0.7071067811865475 ghci> exp (0 :+ pi) + 1 0.0 :+ 1.2246467991473532e-16 ghci> magnitude a 2.23606797749979 ghci> magnitude b 1.4142135623730951 ghci> phase a 0.4636476090008061 ghci> phase b -0.7853981633974483 ghci> conjugate a 2 :+ (-1) ghci> conjugate b 1 :+ 1 ghci> a * conjugate a 5.0 :+ 0.0 ghci> b * conjugate b 2.0 :+ 0.0 ghci> sqrt ((-1) :+ 0) 0.0 :+ 1.0 ghci> (sqrt ((-1) :+ 0)) ^ 2 (-1.0) :+ 0.0 ghci> sqrt ((-2) :+ 0) 0.0 :+ 1.4142135623730951 ghci> (sqrt ((-2) :+ 0)) ^ 2 (-2.0000000000000004) :+ 0.0
ghci> a = fromList [1 :+ 2, 3 :+ 4] :: Vector C ghci> a [1.0 :+ 2.0,3.0 :+ 4.0] ghci> b = fromList [5 :+ 6, 7 :+ 8] :: Vector C ghci> b [5.0 :+ 6.0,7.0 :+ 8.0] ghci> a + b [6.0 :+ 8.0,10.0 :+ 12.0] ghci> a - b [(-4.0) :+ (-4.0),(-4.0) :+ (-4.0)] ghci> a * b [(-7.0) :+ 16.0,(-11.0) :+ 52.0] ghci> a / b [0.2786885245901639 :+ 6.557377049180328e-2,0.4690265486725664 :+ 3.5398230088495575e-2] ghci> c = (1 :+ (-1)) :: C ghci> a + (scalar c) [2.0 :+ 1.0,4.0 :+ 3.0] ghci> a - (scalar c) [0.0 :+ 3.0,2.0 :+ 5.0] ghci> a * (scalar c) [3.0 :+ 1.0,7.0 :+ 1.0] ghci> a / (scalar c) [(-0.5) :+ 1.5,(-0.5) :+ 3.5] ghci> a <.> b 70.0 :+ (-8.0) ghci> b <.> a 70.0 :+ 8.0 ghci> (asRow (cmap conjugate a)) <> (asColumn b) (1><1) [ 70.0 :+ (-8.0) ] ghci> (asRow a) <> (asColumn (cmap conjugate b)) (1><1) [ 70.0 :+ 8.0 ] ghci> udot a b (-18.0) :+ 68.0 ghci> (asRow a) <> (asColumn b) (1><1) [ (-18.0) :+ 68.0 ] ghci> a <.> a 30.0 :+ 0.0 ghci> b <.> b 174.0 :+ 0.0 ghci> norm_2 a 5.477225575051661 ghci> norm_2 b 13.19090595827292
ghci> a = fromLists [[1 :+ 1, 1 :+ (-1)], [2 :+ 3, 2 :+ 1]] :: Matrix C ghci> a (2><2) [ 1.0 :+ 1.0, 1.0 :+ (-1.0) , 2.0 :+ 3.0, 2.0 :+ 1.0 ] ghci> tr a (2><2) [ 1.0 :+ (-1.0), 2.0 :+ (-3.0) , 1.0 :+ 1.0, 2.0 :+ (-1.0) ] ghci> (tr $ tr a) == a True ghci> b = fromLists [[1 :+ 0, 1 :+ 1], [1 :+ (-1), 0 :+ 1]] :: Matrix C ghci> b (2><2) [ 1.0 :+ 0.0, 1.0 :+ 1.0 , 1.0 :+ (-1.0), 0.0 :+ 1.0 ] ghci> tr (a + b) (2><2) [ 2.0 :+ (-1.0), 3.0 :+ (-2.0) , 2.0 :+ (-0.0), 2.0 :+ (-2.0) ] ghci> tr a + tr b (2><2) [ 2.0 :+ (-1.0), 3.0 :+ (-2.0) , 2.0 :+ 0.0, 2.0 :+ (-2.0) ] ghci> tr $ (scalar c) * a (2><2) [ 2.0 :+ (-0.0), 5.0 :+ (-1.0) , 0.0 :+ 2.0, 3.0 :+ 1.0 ] ghci> (scalar (conjugate c)) * tr a (2><2) [ 2.0 :+ 0.0, 5.0 :+ (-1.0) , 0.0 :+ 2.0, 3.0 :+ 1.0 ] ghci> tr (a <> b) (2><2) [ 1.0 :+ 1.0, 5.0 :+ (-2.0) , 1.0 :+ (-3.0), (-2.0) :+ (-7.0) ] ghci> tr b <> tr a (2><2) [ 1.0 :+ 1.0, 5.0 :+ (-2.0) , 1.0 :+ (-3.0), (-2.0) :+ (-7.0) ] ghci> det a (-4.0) :+ 2.0000000000000004 ghci> det (tr a) (-4.0) :+ (-2.0) ghci> conjugate $ det a (-4.0) :+ (-2.0000000000000004) ghci> inv a (2><2) [ (-0.3) :+ (-0.4), 0.3 :+ (-9.999999999999999e-2) , 9.999999999999994e-2 :+ 0.8, (-9.999999999999998e-2) :+ (-0.3) ] ghci> tr $ inv a (2><2) [ (-0.3) :+ 0.4, 9.999999999999994e-2 :+ (-0.8) , 0.3 :+ 9.999999999999999e-2, (-9.999999999999998e-2) :+ 0.3 ] ghci> inv $ tr a (2><2) [ (-0.29999999999999993) :+ 0.4, 9.999999999999998e-2 :+ (-0.7999999999999999) , 0.3 :+ 9.999999999999999e-2, (-9.999999999999999e-2) :+ 0.3 ] ghci> u = fromList [1 :+ 2, 3 :+ 4] :: Vector C ghci> u [1.0 :+ 2.0,3.0 :+ 4.0] ghci> v = fromList [5 :+ 6, 7 :+ 8] :: Vector C ghci> v [5.0 :+ 6.0,7.0 :+ 8.0] ghci> (flatten (a <> (asColumn u))) <.> v 184.0 :+ (-126.0) ghci> u <.> flatten ((tr a) <> (asColumn v)) 184.0 :+ (-126.0)
今回は「固有値 (eigenvalue)」と「固有ベクトル (eigenvector)」について取り上げます。固有値と固有ベクトルは行列の性質を表す重要な指標のひとつで、行列を対角行列に変換する「対角化」の基礎になります。対角行列は単純な構造をしているので、いろいろな計算が簡単になるという利点があります。
たとえば、対角行列 \(A\) の逆行列は対角成分を逆数にしたものですし、\(A\) の n 乗は対角成分を n 乗したものになります。また、\(A\) の行列式も対角成分を掛け算するだけで求めることができます。
与えられた行列の固有値と固有ベクトルを求める問題のことを「固有値問題」といいます。Haskell (hmatrix) には固有値や固有ベクトルを求める関数があらかじめ用意されていますが、Haskell とアルゴリズムのお勉強ということで、あえてプログラムを作ってみましょう。
最初に固有値と固有ベクトルの基本を簡単に説明します。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 よりも少なくなります。
それでは実際に hmatrix で固有値と固有ベクトルを求めてみましょう。hmatrix に用意されている関数 eig, eigevalues を使います。
ghci> import Prelude hiding ((<>)) ghci> :m + Numeric.LinearAlgebra ghci> a = matrix 2 [1,2,2,1] ghci> a (2><2) [ 1.0, 2.0 , 2.0, 1.0 ] ghci> 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 ]) ghci> eigenvalues a [3.0000000000000004 :+ 0.0,(-0.9999999999999996) :+ 0.0] ghci> b = diag $ vector [1,2,3] ghci> b (3><3) [ 1.0, 0.0, 0.0 , 0.0, 2.0, 0.0 , 0.0, 0.0, 3.0 ] ghci> 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 ]) ghci> eigenvalues b [1.0 :+ 0.0,2.0 :+ 0.0,3.0 :+ 0.0] ghci> c = fromLists [[1 :+ 0, 2 :+ (-3)], [2 :+ 3, 4 :+ 0]] :: Matrix C ghci> c (2><2) [ 1.0 :+ 0.0, 2.0 :+ (-3.0) , 2.0 :+ 3.0, 4.0 :+ 0.0 ] ghci> 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 ]) ghci> 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\) がエルミート行列または実対称行列の場合、\(A^{\mathrm{T}}\) は \(A\) になるので (A + tr A) / 2 = A になります。つまり、A のデータ型が Herm に変わるだけで A の各成分の値は変化しません。
簡単な実行例を示します。
ghci> sym a Herm (2><2) [ 1.0, 2.0 , 2.0, 1.0 ] ghci> eigSH $ sym a ([3.0,-1.0],(2><2) [ 0.7071067811865475, -0.7071067811865475 , 0.7071067811865475, 0.7071067811865475 ]) ghci> eigenvaluesSH $ sym a [3.0,-1.0] ghci> sym b Herm (3><3) [ 1.0, 0.0, 0.0 , 0.0, 2.0, 0.0 , 0.0, 0.0, 3.0 ] ghci> 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 ]) ghci> eigenvaluesSH $ sym b [3.0,2.0,1.0] ghci> sym c Herm (2><2) [ 1.0 :+ 0.0, 2.0 :+ (-3.0) , 2.0 :+ 3.0, 4.0 :+ 0.0 ] ghci> 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 ]) ghci> eigenvaluesSH $ sym c [6.405124837953327,-1.4051248379533274]
n 次の正方行列 A の固有ベクトル \(x_1, x_2, \ldots, x_n\) を列に持つ行列 \(X = [x_1, x_2, \ldots, x_n]\) を考えます。すると、\(X\) によって行列 \(A\) を以下のように対角行列 Λ に変換することができます。
固有値と固有ベクトルの定義 \(Ax_i = \lambda_i x_i\) から \(AX\) は次のように変形することができます。
左側から両辺に \(X\) の逆行列を掛け算すれば \(X^{-1}AX = \varLambda\) となります。
それでは実際に ghci で試してみましょう。
ghci> a (2><2) [ 1.0, 2.0 , 2.0, 1.0 ] ghci> (xs, ys) = eigSH $ sym a ghci> xs [3.0,-1.0] ghci> ys (2><2) [ 0.7071067811865475, -0.7071067811865475 , 0.7071067811865475, 0.7071067811865475 ] ghci> (inv ys) <> a <> ys (2><2) [ 3.0, 0.0 , 0.0, -1.0 ] ghci> b = matrix 2 [4,-2,1,1] ghci> b (2><2) [ 4.0, -2.0 , 1.0, 1.0 ] ghci> (xs, ys) = eig b ghci> xs [3.0 :+ 0.0,2.0 :+ 0.0] ghci> ys (2><2) [ 0.8944271909999159 :+ 0.0, 0.7071067811865475 :+ 0.0 , 0.4472135954999579 :+ 0.0, 0.7071067811865475 :+ 0.0 ] ghci> (inv (cmap realPart ys)) <> b <> (cmap realPart ys) (2><2) [ 3.0, 0.0 , 0.0, 2.0 ]
数学 (線形代数) の世界では、正則行列 \(P\) を用いて正方行列 \(A\) を \(P^{-1}AP\) に変換することを「相似変換」といい、\(B = P^{-1}AP\) が成り立つ正方行列 \(A, B\) の関係を「相似」といいます。相似な行列 \(A\) と \(B\) の間は様々な性質が保存されます。保存される主な性質を以下に示します。
つまり、相似変換を行ってもこれらの値が変化することはありません。たとえば、正方行列 \(A\) を相似変換により対角行列に変換します。対角行列は対角成分が固有値になるので、トレースは固有値の総和になります。相似変換でトレースは保存されるので、行列 \(A\) のトレースも固有値の総和であることがわかります。
要素がすべて実数の対称行列を「実対称行列」といいます。実対称行列には次に示す性質があります。
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
このように、実対称行列の固有ベクトルは一次独立になります。
それでは実際に ghci で試してみましょう。
ghci> a (2><2) [ 1.0, 2.0 , 2.0, 1.0 ] ghci> (xs, ys) = eigSH $ sym a ghci> xs [3.0,-1.0] ghci> ys (2><2) [ 0.7071067811865475, -0.7071067811865475 , 0.7071067811865475, 0.7071067811865475 ] ghci> zs = toColumns ys ghci> (zs !! 0) <.> (zs !! 1) 0.0 ghci> (tr ys) <> a <> ys (2><2) [ 2.9999999999999996, 0.0 , 0.0, -0.9999999999999998 ] ghci> b = matrix 3 [1,4,5,4,2,6,5,6,3] ghci> b (3><3) [ 1.0, 4.0, 5.0 , 4.0, 2.0, 6.0 , 5.0, 6.0, 3.0 ] ghci> (xs, ys) = eigSH $ sym b ghci> xs [12.175971065046912,-2.5072879670936405,-3.668683097953268] ghci> ys (3><3) [ 0.49659978454619114, 0.8095854617397508, 0.31298567719355963 , 0.5773502691896256, -0.5773502691896257, 0.5773502691896255 , 0.6481167492476515, -0.10600965430705465, -0.7541264035547062 ] ghci> zs = toColumns ys ghci> (zs !! 0) <.> (zs !! 1) 1.249000902703301e-16 ghci> (zs !! 0) <.> (zs !! 2) -2.220446049250313e-16 ghci> (zs !! 1) <.> (zs !! 2) 1.942890293094024e-16 ghci> (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 ] ghci> c = matrix 4 [5,1,1,1,1,6,1,1,1,1,7,1,1,1,1,8] ghci> 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 ] ghci> (xs, ys) = eigSH $ sym c ghci> xs [9.803886359051248,6.507748705363649,5.3922752902729805,4.296089645312117] ghci> 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 ] ghci> (inv ys) <> c <> ys (4><4) [ 9.8038863590, -2.2204460492e-15, 9.9920072216e-16, -3.33066907387e-16 , 8.8817841970e-16, 6.5077487053, -7.7715611723e-16, 2.2204460492e-16 , 6.6613381477e-16, -1.55431223447e-15, 5.3922752902, -3.1918911957e-16 , -5.5511151231e-16, 3.33066907387e-16, -1.80411241501e-16, 4.2960896453 ]
行列の固有値と固有ベクトルを求める簡単な方法に「累乗法 (power method)」があります。「べき乗法」と呼ばれることもあります。参考文献『Ｃ言語による最新アルゴリズム事典』によると、『適当なベクトル x を選び、これに行列 A を何度も掛け算すると、x に含まれる A の絶対値最大の固有値に対応する固有ベクトルが最も速く成長する』とのことです。x が単位ベクトルとすると、x が収束したときの Ax の大きさが固有値の絶対値 |λ| となります。
それでは実際に試してみましょう。
ghci> iter n x = if n == 0 then [x] else let x1 = a <> x in x : iter (n - 1) (x1 / (scalar (norm_2 x1))) ghci> 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 ]] ghci> 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 商」を用いるほうが良いようです。
固有値の定義により \((Ax, x) = (\lambda x, x)\) になります。内積は結合法則 \((kx, y) = (x, ky) = k(x, y)\) が成り立つので、\(\frac{(Ax, x)}{(x, x)} = \frac{\lambda (x, x)}{(x, x)} = \lambda\) となり、固有値 \(\lambda\) を求めることができます。あとは、\(|\lambda_i - \lambda_{i-1}|\) が許容誤差 \(\epsilon\) の中に入るまで繰り返し \(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 を再帰呼び出しして処理を繰り返します。
それでは実際に試してみましょう。
ghci> power $ matrix 2 [2,1,1,2] Just (2.999999999426406,(2><1) [ 0.7071107728137574 , 0.7071027895368046 ],11) ghci> 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 番目の固有値を \(\lambda_1\) と \(\lambda_2\) とすると、\(|\frac{\lambda_2}{\lambda_1}|\) が 1 に近い値だと累乗法の収束は極端に遅くなります。ご注意くださいませ。
参考文献『Ｃ言語による最新アルゴリズム事典』によると、行列 \(A\) が対称行列であり、固有値 \(\lambda_i \ (|\lambda_1| \gt |\lambda_2| \gt \ldots \gt |\lambda_n|)\) がすべて異なる場合、固有ベクトル \(x_i\) は互いに直交し、次式が成り立つとのことです。
したがって、\(\lambda_1\) と \(x_1\) が求まったならば、行列 \((A - \lambda_1 x_1 {x_1}^{T})\) に再度累乗法を適用すると、次に大きな固有値 \(\lambda_2\) と固有ベクトル \(x_2\) を求めることができます。これを繰り返すことで、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 を連結するだけです。
それでは実際に試してみましょう。
ghci> powerAll $ matrix 2 [1,2,2,1] Just ([2.9999999988528114,-1.0000000003823961],(2><2) [ 0.7071027895368046, -0.7070948061697462 , 0.7071107728137574, 0.7071187560005527 ]) ghci> powerAll $ matrix 2 [2,1,1,2] Just ([2.999999999426406,1.0000000001911984],(2><2) [ 0.7071107728137574, 0.7070948061697185 , 0.7071027895368046, -0.70711875600058 ]) ghci> 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 ]) ghci> powerAll $ matrix 4 [5,1,1,1,1,6,1,1,1,1,7,1,1,1,1,8] Just ([9.80388635545,6.50774868672,5.39227529304,4.29608965698],(4><4) [ 0.3320050413176, 0.1359552625975, 0.2259674702248, 0.905657440661 , 0.401118359142, 0.2261797515731, 0.801715189075, -0.3810552949782 , 0.506575827611, 0.671328714829, -0.51764254730, -0.15732141501 , 0.687210044584, -0.692586921039, -0.1955444993576, -9.91535915646e-2 ])
正常に動作しているようです。興味のある方はいろいろ試してみてください。
\(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『行列の固有値問題 (PDF)』によると、逆行列を求めるよりも連立方程式 \(Ax_{i+1} = x_i\) を解くほうが良いようです。この場合、行列を 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 を返します。
それでは実行してみましょう。
ghci> powerInv $ matrix 2 [1,2,2,1] Just (-1.0000000003823961,(2><1) [ -0.7071027895368047 , 0.7071107728137576 ],11) ghci> powerAll $ matrix 2 [1,2,2,1] Just ([2.9999999988528114,-1.0000000003823961],(2><2) [ 0.7071027895368046, -0.7070948061697462 , 0.7071107728137574, 0.7071187560005527 ]) ghci> powerInv $ matrix 2 [2,1,1,2] Just (1.000000000191198,(2><1) [ 0.7071107728137573 , -0.7071027895368047 ],11) ghci> powerAll $ matrix 2 [2,1,1,2] Just ([2.999999999426406,1.0000000001911984],(2><2) [ 0.7071107728137574, 0.7070948061697185 , 0.7071027895368046, -0.70711875600058 ]) ghci> powerInv $ matrix 3 [1,4,5,4,2,6,5,6,3] Just (-2.5072879960642904,(3><1) [ -0.8096263206940646 , 0.5772748809066349 , 0.10610811803826374 ],21) ghci> 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 ]) ghci> 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) ghci> powerAll $ matrix 4 [5,1,1,1,1,6,1,1,1,1,7,1,1,1,1,8] Just ([9.80388635545,6.50774868672,5.39227529304,4.29608965698],(4><4) [ 0.3320050413176, 0.1359552625975, 0.2259674702248, 0.905657440661 , 0.401118359142, 0.2261797515731, 0.801715189075, -0.3810552949782 , 0.506575827611, 0.671328714829, -0.51764254730, -0.15732141501 , 0.687210044584, -0.692586921039, -0.1955444993576, -9.91535915646e-2 ])
正常に動作していますね。