前回は固有値と固有ベクトルの基礎知識を簡単に説明しました。今回は固有値と固有ベクトルを求めるプログラムを作ってみましょう。
行列の固有値と固有ベクトルを求める簡単な方法に「累乗法 (power method)」があります。「べき乗法」と呼ばれることもあります。参考文献『Ｃ言語による最新アルゴリズム事典』によると、『適当なベクトル \(x\) を選び、これに行列 \(A\) を何度も掛け算すると、\(x\) に含まれる \(A\) の絶対値最大の固有値に対応する固有ベクトルが最も速く成長する』とのことです。\(x\) が単位ベクトルとすると、\(x\) が収束したときの \(Ax\) の大きさが固有値の絶対値 \(|\lambda|\) となります。
それでは実際に試してみましょう。
>>> a = np.array([[1., 2.], [2., 1.]]) >>> x = np.array([1., 0]) >>> for i in range(16): ... x1 = a @ x ... k = np.linalg.norm(x1) ... print(i, x, k) ... x = x1 / k ... 0 [ 1. 0.] 2.2360679775 1 [ 0.4472136 0.89442719] 2.86356421266 2 [ 0.78086881 0.62469505] 2.98369553145 3 [ 0.6804511 0.73279349] 2.99817295964 4 [ 0.71578195 0.69832385] 2.99979680303 5 [ 0.70419091 0.71001067] 2.99997742018 6 [ 0.70807608 0.70613615] 2.9999974911 7 [ 0.70678338 0.70743003] 2.99999972123 8 [ 0.70721455 0.706999 ] 2.99999996903 9 [ 0.70707086 0.70714271] 2.99999999656 10 [ 0.70711876 0.70709481] 2.99999999962 11 [ 0.70710279 0.70711077] 2.99999999996 12 [ 0.70710811 0.70710545] 3.0 13 [ 0.70710634 0.70710722] 3.0 14 [ 0.70710693 0.70710663] 3.0 15 [ 0.70710673 0.70710683] 3.0 >>> x array([ 0.7071068 , 0.70710676]) >>> np.linalg.norm(a @ x) 2.9999999999999991
このように、\(x_{i+1} = Ax_{i}\) を繰り返すことにより、固有値と固有ベクトルの値は理論値に近づいていくことがわかります。なお、固有ベクトルは定数倍してもいいので、-1 倍した固有ベクトルが求まる場合もあります。また、実際に固有値を求めるときは、以下に示す「Rayleigh 商」を用いるほうが良いようです。
固有値の定義により \((Ax, x) = (λ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\) を掛け算すればいいわけです。
NumPy を使うと、累乗法はとても簡単にプログラムすることができます。次のリストを見てください。
リスト : 累乗法
def power(a, max_iter = 512, eps = 1e-8):
x = np.full(len(a), 0.0)
x[0] = 1.0
k = 0
for i in range(max_iter):
# print(i, k, x)
x1 = a @ x
n = np.linalg.norm(x1)
k1 = (x1 @ x) / (x @ x)
if abs(k1 - k) < eps:
return k1, x1 / n
k = k1
x = x1 / n
引数 max_iter が繰り返しの最大回数、eps が許容誤差です。変数 x が固有ベクトル、変数 k が固有値を表します。x は大きさ 1 のベクトルに初期化します。今回は単純に [1, 0, ...] としましたが、他の方法で初期化してもかまいません。
あとはループの中で x1 = a @ x を計算し、k1 = (x1 @ x) / (x @ x) で Rayleigh 商を求めます。変数 n にはベクトル x1 の大きさをセットします。abs(k1 - k) が eps 未満になったら収束したとみなして、固有値 k1 と固有ベクトル x1 / n を返します。そうでなければ、k を k1 に x を x1 / n に更新して処理を繰り返します。
それでは実際に試してみましょう。ご参考までに、k と x の値を表示しています。
>>> power(np.array([[2, 1], [1, 2]])) 0 0 [ 1. 0.] 1 2.0 [ 0.89442719 0.4472136 ] 2 2.8 [ 0.78086881 0.62469505] 3 2.9756097561 [ 0.73279349 0.6804511 ] 4 2.99726027397 [ 0.71578195 0.69832385] 5 2.99969521487 [ 0.71001067 0.70419091] 6 2.9999661304 [ 0.70807608 0.70613615] 7 2.99999623665 [ 0.70743003 0.70678338] 8 2.99999958185 [ 0.70721455 0.706999 ] 9 2.99999995354 [ 0.70714271 0.70707086] 10 2.99999999484 [ 0.70711876 0.70709481] (2.9999999994264059, array([ 0.70711077, 0.70710279])) >>> power(np.array([[1, 4, 5], [4, 2, 6], [5, 6, 3]])) 0 0 [ 1. 0. 0.] 1 1.0 [ 0.15430335 0.6172134 0.77151675] 2 10.2380952381 [ 0.56819047 0.56819047 0.59524716] 3 12.0549048316 [ 0.4795596 0.57547152 0.6624614 ] 4 12.1682902317 [ 0.50058982 0.57864665 0.64387724] 5 12.1754173848 [ 0.49562055 0.57679433 0.64936013] 6 12.1759269017 [ 0.49684795 0.57755071 0.64774786] 7 12.1759673145 [ 0.49653459 0.577283 0.64822661] 8 12.1759707355 [ 0.49661745 0.57737194 0.6480839 ] 9 12.1759710356 [ 0.49659487 0.57734345 0.64812659] 10 12.1759710624 [ 0.49660118 0.57735238 0.64811379] (12.175971064806813, array([ 0.49659938, 0.57734962, 0.64811764]))
このように、累乗法を使うと絶対値で最大の固有値と固有ベクトルを簡単に求めることができます。ただし、累乗法には収束が遅いという欠点があります。絶対値最大の固有値と 2 番目の固有値を \(\lambda_{1}\) と \(\lambda_{2}\) とすると、\(|\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}}^{\mathrm{T}})\) に再度累乗法を適用すると、次に大きな固有値 \(\lambda_{2}\) と固有ベクトル \(x_{2}\) を求めることができます。これを繰り返すことで、n 次の正方行列であれば n 個の固有値と固有ベクトルを求めることができます。
プログラムは次のようになります。
リスト : 累乗法 (2)
def power_all(a, max_iter = 512, eps = 1e-8):
n = len(a)
a = a.astype(np.float_)
xs = np.empty(n) # 固有値
yss = np.empty((n, n)) # 固有ベクトル
for i in range(n):
x, ys = power(a, max_iter, eps)
xs[i] = x
yss[:, i] = ys
zs = ys.reshape((n, 1))
a -= x * zs @ zs.T
return xs, yss
変数 xs は固有値を格納する配列、yss は固有ベクトルを格納する配列です。power() で絶対値最大の固有値 x と固有ベクトル ys を求め、xs と yss にセットします。ys は列ベクトルとして yss にセットしていることに注意してください。
そして、reshape() で ys を列ベクトルに変換して変数 zs にセットし、行列 a から x * zs @ zs.T を引き算します。これで次の固有値と固有ベクトルを求めることができます。
それでは実際に試してみましょう。
>>> power_all(np.array([[1, 2], [2, 1]]))
(array([ 3., -1.]), array([[ 0.70710279, -0.70709481],
[ 0.70711077, 0.70711876]]))
>>> power_all(np.array([[2, 1], [1, 2]]))
(array([ 3., 1.]), array([[ 0.70711077, 0.70709481],
[ 0.70710279, -0.70711876]]))
>>> power_all(np.array([[1, 4, 5], [4, 2, 6], [5, 6, 3]]))
(array([ 12.17597106, -3.66868309, -2.50728797]),
array([[ 0.49659938, -0.31301741, -0.80956631],
[ 0.57734962, -0.57732402, 0.57738552],
[ 0.64811764, 0.75413333, 0.10596392]]))
>>> 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.]])
>>> power_all(c)
(array([ 13.39054122, 9.5403944 , 8.43473667, 7.35663186, 6.27769584]),
array([[ 0.29109109, 0.09514325, 0.13338346, 0.2190289 , 0.91675417],
[ 0.33664233, 0.132605 , 0.22641654, 0.83277347, -0.35258174],
[ 0.39909575, 0.21875529, 0.74684927, -0.461795 , -0.1477534 ],
[ 0.49000445, 0.62317649, -0.57460249, -0.18066388, -0.09349382],
[ 0.63447312, -0.73291258, -0.20734644, -0.11234155, -0.0683794 ]]))
正常に動作しているようです。興味のある方はいろいろ試してみてください。
\(Ax = \lambda x\) の両辺に左側から \(A^{-1}\) を掛け算すると、\(x = \lambda A^{-1}x\) => \(A^{-1}x = (1/\lambda)x\) となります。\(A^{-1}\) に累乗法を適用すると \(1 / \lambda\) の最大値、つまり最小の固有値とその固有ベクトルを求めることができます。これを「逆累乗法」とか「逆べき乗法」といいます。
具体的には、\(x_{i+1} = A^{-1}x_{i}\) のように \(A^{-1}\) を掛け算していくことになりますが、参考 URL 2 『行列の固有値問題』によると、逆行列を求めるよりも連立方程式 \(Ax_{i+1} = x_{i}\) を解くほうが良いようです。この場合、行列を LU 分解しておいたほうが効率的です。
プログラムは次のようになります。
リスト : 逆累乗法
def power_inv(a, max_iter = 512, eps = 1e-8):
_, a1, idx = lu_pv(a)
x = np.full(len(a), 0.0)
x[0] = 1.0
k = 0
for i in range(max_iter):
x1 = lu_solver(a1, x[idx])
n = np.linalg.norm(x1)
k1 = (x1 @ x) / (x @ x)
if abs(k1 - k) < eps:
return 1 / k1, x1 / n
k = k1
x = x1 / n
最初に関数 lu_pv() で行列 a を LU 分解します。この関数は拙作のページ「連立方程式の解法」で作成したものです。ループの中で x1 = a^{-1} @ x を計算する代わりに lu_solver() を呼び出します。この関数も拙作のページで作成したものです。収束判定は累乗法のプログラムと同じです。収束したら固有値 1 / k1 と固有ベクトル x1 / n を返します。
それでは実行してみましょう。
>>> power_inv(np.array([[2, 1], [1, 2]]))
(1.0000000001911982, array([ 0.70711077, -0.70710279]))
>>> power_all(np.array([[2, 1], [1, 2]]))
(array([ 3., 1.]), array([[ 0.70711077, 0.70709481],
[ 0.70710279, -0.70711876]]))
>>> power_inv(np.array([[1, 2], [2, 1]]))
(-1.0000000003823963, array([-0.70710279, 0.70711077]))
>>> power_all(np.array([[1, 2], [2, 1]]))
(array([ 3., -1.]), array([[ 0.70710279, -0.70709481],
[ 0.70711077, 0.70711876]]))
>>> power_inv(np.array([[1, 4, 5], [4, 2, 6], [5, 6, 3]]))
(-2.5072879960642909, array([-0.80962632, 0.57727488, 0.10610812]))
>>> power_all(np.array([[1, 4, 5], [4, 2, 6], [5, 6, 3]]))
(array([ 12.17597106, -3.66868309, -2.50728797]), array([[ 0.49659938, -0.31301741, \
-0.80956631],
[ 0.57734962, -0.57732402, 0.57738552],
[ 0.64811764, 0.75413333, 0.10596392]]))
>>> c = np.ones((4, 4))
>>> c += np.diag([4,5,6,7])
>>> c
array([[ 5., 1., 1., 1.],
[ 1., 6., 1., 1.],
[ 1., 1., 7., 1.],
[ 1., 1., 1., 8.]])
>>> power_inv(c)
(4.2960898996174643, array([ 0.90578078, -0.38061757, -0.15760282, -0.09926121]))
>>> power_all(c)
(array([ 9.80388636, 6.50774869, 5.39227529, 4.29608966]),
array([[ 0.33200504, 0.13595526, 0.22596747, 0.90565744],
[ 0.40111836, 0.22617975, 0.80171519, -0.38105529],
[ 0.50657583, 0.67132871, -0.51764255, -0.15732142],
[ 0.68721004, -0.69258692, -0.1955445 , -0.09915359]]))
正常に動作していますね。
固有値 \(\lambda_{i}\) の近似値 \(\lambda'\) が分かっている場合、それよりも高い精度の固有値と固有ベクトルを逆累乗法を使って求めることができます。これを「逆反復法」といいます。原理を簡単に説明しましょう。
行列 \(A\) の固有値が \(\lambda_{1}, \ldots, \lambda_{i}, \ldots, \lambda_{n}\) とすると、行列 \(A' = A - \lambda'I\) の固有値は \(\lambda_{1} - \lambda', \ldots, \lambda_{i} - \lambda', \ldots, \lambda_{n} - \lambda'\) となります。\(\lambda'\) は \(\lambda_{i}\) に近い値なので、\(|\lambda_{i} - \lambda|\) は最小値であることが期待できます。つまり、行列 \(A'\) に逆累乗法を適用すれば、\(\lambda_{i} - \lambda’\) を求めることができるわけです。この値を k とすれば、\(\lambda_{i} = \lambda' + k\) となります。
なお、同様な原理で累乗法を高速化する手法があります。\(B = A - pI\) としたとき、\(|\lambda_{2} / \lambda_{1}| \gt |(\lambda_{2} - p) / (\lambda_{1} - p)|\) になるならば、行列 \(B\) に累乗法を適用したほうが収束が速くなるはずです。これを「原点移動 (シフト法)」といいます。ただし、参考文献によっては、逆反復法のことをシフト法と呼ぶことがあるようです。まあ、どちらの方法も対角成分の値をシフトするので、シフト法というと混乱するような気がします。本ページでは前者を逆反復法、後者を原点移動と呼ぶことにします。
プログラムは簡単です。次のリストを見てください。
リスト : 逆反復法と原点移動
# 逆反復法
def power_inv_shift(a, d, max_iter = 512, eps = 1e-8):
k, xs = power_inv(a - d * np.eye(len(a)), max_iter, eps)
return k + d, xs
# 原点移動
def power_shift(a, d, max_iter = 512, eps = 1e-8):
k, xs = power(a - d * np.eye(len(a)), max_iter, eps)
return k + d, xs
どの関数も引数 a が行列で、d がシフトする値です。NumPy の場合、a - dI は a - d * np.eys(len(a)) で求めることができます。あとは、適切な関数 (power_inv, power) を呼び出して、求めた固有値に d を加算するだけです。
それでは実際に試してみましょう。最初は逆反復法です。ご参考までに収束の様子を表示しています。
>>> a = np.array(np.array([[1, 4, 5], [4, 2, 6], [5, 6, 3]]))
>>> a
array([[1, 4, 5],
[4, 2, 6],
[5, 6, 3]])
>>> np.linalg.eigh(a)
(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]]))
>>> power_inv_shift(a, 12.1)
0 0 [ 1. 0. 0.]
1 3.19503971656 [ 0.48876527 0.58046886 0.65127906]
2 13.161833091 [ 0.49664006 0.57733326 0.64810104]
3 13.1629061303 [ 0.49659958 0.57735036 0.64811683]
4 13.1629061588 [ 0.49659979 0.57735027 0.64811675]
(12.175971065046905, array([ 0.49659978, 0.57735027, 0.64811675]))
>>> power_inv_shift(a, -3.6)
0 0 [ 1. 0. 0.]
1 -0.810810810811 [-0.17561811 -0.65856793 0.73174214]
2 -14.1604342616 [ 0.32125734 0.5714376 -0.75515084]
3 -14.5580068367 [-0.31246551 -0.57772093 0.75405824]
4 -14.5596164304 [ 0.31301836 0.57732696 -0.75413068]
5 -14.5596227903 [-0.31298362 -0.57735173 0.75412613]
6 -14.5596228154 [ 0.31298581 0.57735018 -0.75412642]
(-3.668683097953267, array([-0.31298567, -0.57735027, 0.7541264 ]))
>>> power_inv_shift(a, -2.5)
0 0 [ 1. 0. 0.]
1 -90.0 [-0.81018636 0.57613253 0.10802485]
2 -137.211669368 [ 0.80959024 -0.5773415 -0.10602089]
3 -137.212474607 [-0.80958549 0.57735022 0.10600973]
4 -137.212474638 [ 0.80958546 -0.57735027 -0.10600965]
(-2.5072879670936414, array([-0.80958546, 0.57735027, 0.10600965]))
>>> b = np.ones((4, 4))
>>> b += np.diag([4, 5, 6, 7])
>>> b
array([[ 5., 1., 1., 1.],
[ 1., 6., 1., 1.],
[ 1., 1., 7., 1.],
[ 1., 1., 1., 8.]])
>>> np.linalg.eigh(b)
(array([ 4.29608965, 5.39227529, 6.50774871, 9.80388636]),
array([[ 0.90568356, -0.22590297, 0.13594053, 0.33200196],
[-0.38096261, -0.80178162, 0.22610179, 0.40111308],
[-0.15738124, 0.51753551, 0.67140433, 0.50656131],
[-0.09917619, 0.19562996, -0.69254197, 0.68722531]]))
>>> power_inv_shift(b, 5.3)
0 0 [ 1. 0. 0. 0.]
1 -0.224250843421 [-0.08583621 0.90393881 -0.38740235 -0.15951861]
2 9.39906651479 [ 0.25711316 0.78946174 -0.51989385 -0.2008815 ]
3 10.8234496468 [ 0.22312286 0.80302951 -0.51685641 -0.19549741]
4 10.8370219879 [ 0.22616501 0.80167731 -0.51756663 -0.19567235]
5 10.8371364342 [ 0.22587936 0.801792 -0.51753032 -0.19562843]
6 10.8371373992 [ 0.22590517 0.80178073 -0.51753581 -0.19563028]
(5.3922752902735738, array([ 0.22590277, 0.80178171, -0.51753547, -0.19562994]))
>>> power_inv_shift(b, 6.5)
0 0 [ 1. 0. 0. 0.]
1 2.0 [ 0.11396058 0.22792115 0.68376346 -0.68376346]
2 128.961038961 [ 0.13603794 0.226141 0.67135609 -0.69255681]
3 129.053814725 [ 0.13594015 0.22610146 0.67140469 -0.69254181]
4 129.053816485 [ 0.13594053 0.22610179 0.67140433 -0.69254197]
(6.5077487053636487, array([ 0.13594053, 0.22610179, 0.67140433, -0.69254197]))
逆反復法の場合、固有値の近似が良ければ収束も速くなるようで、高い精度の固有値とその固有ベクトルを高速に求めることができます。
次はシフト法です。
>>> a
array([[1, 4, 5],
[4, 2, 6],
[5, 6, 3]])
>>> power(a)
0 0 [ 1. 0. 0.]
1 1.0 [ 0.15430335 0.6172134 0.77151675]
2 10.2380952381 [ 0.56819047 0.56819047 0.59524716]
3 12.0549048316 [ 0.4795596 0.57547152 0.6624614 ]
4 12.1682902317 [ 0.50058982 0.57864665 0.64387724]
5 12.1754173848 [ 0.49562055 0.57679433 0.64936013]
6 12.1759269017 [ 0.49684795 0.57755071 0.64774786]
7 12.1759673145 [ 0.49653459 0.577283 0.64822661]
8 12.1759707355 [ 0.49661745 0.57737194 0.6480839 ]
9 12.1759710356 [ 0.49659487 0.57734345 0.64812659]
10 12.1759710624 [ 0.49660118 0.57735238 0.64811379]
(12.175971064806813, array([ 0.49659938, 0.57734962, 0.64811764]))
>>> power_shift(a, 3)
0 0 [ 1. 0. 0.]
1 -2.0 [-0.2981424 0.59628479 0.74535599]
2 1.15555555556 [ 0.89200774 0.35680309 0.27751352]
3 4.49115913556 [ 0.12462143 0.58955523 0.79805647]
・・・省略・・・
30 9.17597100757 [ 0.49661378 0.57737533 0.6480837 ]
31 9.17597103469 [ 0.49658965 0.57733203 0.64814076]
32 9.17597104901 [ 0.49660713 0.57736354 0.6480993 ]
(12.175971056579188, array([ 0.49659446, 0.57734061, 0.64812943]))
>>> power_shift(a, -3)
0 0 [ 1. 0. 0.]
1 4.0 [ 0.52981294 0.52981294 0.66226618]
2 15.1228070175 [ 0.49837286 0.57706331 0.64701038]
3 15.1759039829 [ 0.49662808 0.57728693 0.64815149]
4 15.1759709733 [ 0.49660199 0.57735059 0.64811477]
5 15.1759710649 [ 0.4965998 0.57735017 0.64811682]
(12.175971065046681, array([ 0.49659979, 0.57735027, 0.64811675]))
>>> b
array([[ 5., 1., 1., 1.],
[ 1., 6., 1., 1.],
[ 1., 1., 7., 1.],
[ 1., 1., 1., 8.]])
>>> power(b)
0 0 [ 1. 0. 0. 0.]
1 5.0 [ 0.94491118 0.18898224 0.18898224 0.18898224]
2 6.5 [ 0.75537858 0.35071148 0.37768929 0.4046671 ]
・・・省略・・・
21 9.80388631692 [ 0.33201277 0.40113193 0.50661039 0.68717291]
22 9.8038863405 [ 0.33200906 0.40112538 0.50659404 0.68719058]
23 9.80388635088 [ 0.33200663 0.40112113 0.50658311 0.68720229]
(9.8038863554521107, array([ 0.33200504, 0.40111836, 0.50657583, 0.68721004]))
>>> power_shift(b, 4)
0 0 [ 1. 0. 0. 0.]
1 1.0 [ 0.5 0.5 0.5 0.5]
2 5.5 [ 0.35634832 0.4454354 0.53452248 0.62360956]
・・・省略・・・
10 5.80388620677 [ 0.33201468 0.40113442 0.50662344 0.68716091]
11 5.80388633062 [ 0.33200744 0.40112224 0.5065882 0.6871975 ]
12 5.80388635374 [ 0.33200433 0.40111702 0.50657294 0.6872133 ]
(9.8038863580604492, array([ 0.33200298, 0.40111478, 0.50656634, 0.68722012]))
>>> power_shift(b, -4)
0 0 [ 1. 0. 0. 0.]
1 9.0 [ 0.98198051 0.10910895 0.10910895 0.10910895]
2 9.78571428571 [ 0.91057245 0.22764311 0.23848326 0.24932341]
・・・省略・・・
31 13.8038863159 [ 0.33201513 0.40113615 0.50661728 0.68716422]
32 13.8038863341 [ 0.33201185 0.40113038 0.50660412 0.68717889]
33 13.8038863446 [ 0.3320094 0.40112607 0.50659403 0.68719002]
(9.8038863506769474, array([ 0.33200757, 0.40112284, 0.50658631, 0.68719848]))
原点移動の場合、シフトする値によって収束回数が増えたり減ったりします。何かしらの選択基準があればよいのですが、よくわかりませんでした。シフト値の最適値は行列によって変わるでしょうから、明確な選択基準がなければ (M.Hiroi が知らないだけかもしれませんが)、ちょっと使いにくい方法のような気がします。興味のある方はいろいろ試してみてください。
実対称行列の固有値と固有ベクトルを求める簡単な方法に「ヤコビ法 (Jacobi method)」があります。なお、連立一次方程式を反復法で解く方法にもヤコビ法があります。拙作のページ「連立方程式の解法 (2)」で説明しているので、興味のある方はお読みくださいませ。
ヤコビ法のポイントは、直交行列 \(X_{1}\) を使って相似変換 \(A_{1} = {X_{1}}^{\mathrm{T}}AX_{1}\) を行うところです。このとき、\(A\) の非対角成分のひとつが 0 になるように直交行列 \(X_{1}\) を定めます。つまり、相似変換 \(A_{n} = {X_{n}}^{\mathrm{T}}A_{n-1}X_{n}\) を繰り返して \(A\) を対角化するわけです。対角化が完了したとき (非対角成分がすべて 0 になったとき)、\(A_{n}\) の対角成分が固有値になり、\(X_{1}X_{2} \ldots X_{n}\) が固有ベクトルになります。これを式で表すと次のようになります。
直交行列 \(X\) には「ギブンス回転」を使います。ギブスン回転は N 次元の行列 \(A\) の (i, j) 平面の回転を表します。
X[i, i] = cos(x) X[j, j] = cos(x) X[i, j] = sin(x) X[j, i] = -sin(x) X の対角成分 は 1, 残りの非対角成分は 0
i j
+---------------------------------------------
| 1, 0, ...
| 0, 1, ...
| ...
i | cos(x), ..., sin(x),
|
| ...
|
j | -sin(x), ..., cos(x),
| ...
| ..., 1, 0
| ..., 0, 1
+---------------------------------------------
二次元行列の場合、ギブンス回転は回転行列と同じになります。ここでは二次元行列で説明します。非対角成分を 0 にする角度 r は次のように求めることができます。
行列 A = [[x, z], X = [[ cos(r), sin(r)],
[z, y]] [-sin(r), cos(r)]]
XTAX = [[cos(r), -sin(r)], [[x, z], [[cos(r), sin(r)],
[sin(r), cos(r)] @ [z, y]] @ [-sin(r), cos(r)]]
非対角成分 => z * (cos^{2}(r) - sin^{2}(r)) + (x - y) * sin(r)cos(r) = 0
三角関数の公式
tan(x) = sin(x) / cons(x)
cos2(x) - sin2(x) = cos(2x)
sin(x)cos(x) = sin(2x) / 2
を使って
z * cos(2r) = (y - x) * sin(2r) / 2
2 * z / (y - x) = sin(2r) / cons(2r) = tan(2r)
r = arctan(2 * z / (y - x)) / 2, (x != y)
r = π/ 4, (x == y)
実際の数値計算では、角度ではなくて sir(r) と cos(r) の値を求めたほうが誤差が少なくなるようです。sin(r) と cos(r) の求め方は参考 URL 4 『固有値問題 (1) 対称行列の固有値』がわかりやすくまとまっていると思います。優れたドキュメントとプログラムを公開されている fussy さんに感謝いたします。今回はプログラムを簡単にするため、arctan() で角度を求めることにしましょう。
それでは、これで本当に 0 になるのか実際に試してみましょう。
>>> def givens(x): return np.array([[np.cos(x), np.sin(x)], [-np.sin(x), np.cos(x)]])
...
>>> a = np.array([[2, 1],[1, 3]])
>>> a
array([[2, 1],
[1, 3]])
>>> x = givens(0.5 * np.arctan(2 / (3 - 2)))
>>> x
array([[ 0.85065081, 0.52573111],
[-0.52573111, 0.85065081]])
>>> x.T @ a @ x
array([[ 1.38196601e+00, -1.11022302e-16],
[ 2.22044605e-16, 3.61803399e+00]])
>>> np.linalg.eig(a)
(array([ 1.38196601, 3.61803399]), array([[-0.85065081, -0.52573111],
[ 0.52573111, -0.85065081]]))
>>> b = np.array([[2, 1],[1, 2]])
>>> b
array([[2, 1],
[1, 2]])
>>> y = givens(np.pi / 4)
>>> y
array([[ 0.70710678, 0.70710678],
[-0.70710678, 0.70710678]])
>>> y.T @ b @ y
array([[ 1.00000000e+00, 2.22044605e-16],
[ 0.00000000e+00, 3.00000000e+00]])
>>> np.linalg.eig(b)
(array([ 3., 1.]), array([[ 0.70710678, -0.70710678],
[ 0.70710678, 0.70710678]]))
>>> c = np.random.rand(4).reshape((2,2))
>>> c
array([[ 0.79031302, 0.44585626],
[ 0.49301966, 0.0917365 ]])
>>> c[0, 1] = c[1, 0]
>>> c
array([[ 0.79031302, 0.49301966],
[ 0.49301966, 0.0917365 ]])
>>> z = givens(0.5 * np.arctan(2 * c[0, 1] / (c[1, 1] - c[0, 0])))
>>> z
array([[ 0.88828207, -0.45929833],
[ 0.45929833, 0.88828207]])
>>> z.T @ c @ z
array([[ 1.04523554e+00, -5.55111512e-17],
[ 0.00000000e+00, -1.63186025e-01]])
>>> np.linalg.eig(c)
(array([ 1.04523554, -0.16318603]), array([[ 0.88828207, -0.45929833],
[ 0.45929833, 0.88828207]]))
確かにギブンス回転を使って非対角成分を 0 にすることができました。
それではプログラムを作りましょう。ヤコビ法は NumPy を使うと簡単にプログラムすることができます。次のリストを見てください。
リスト : ヤコビ法
def jacobi(a, max_iter = 512, e = 1e-14):
n = len(a)
ks = np.eye(n) # 固有ベクトル
for i in range(max_iter):
# 絶対値最大の非対角成分を探す
x, y = divmod(np.fabs(np.triu(a, 1)).argmax(), n)
# print(i, x, y, a[x, y], a[y, x])
if np.fabs(a[x, y]) < e:
return np.diag(a), ks
# 角度を求める
d = a[y, y] - a[x, x]
if d == 0.0:
r = np.pi / 4
else:
r = 0.5 * np.arctan(2 * a[x, y] / d)
# 直交行列の生成
xs = np.eye(n)
xs[x, x] = np.cos(r)
xs[y, y] = np.cos(r)
xs[x, y] = np.sin(r)
xs[y, x] = - np.sin(r)
#
a = xs.T @ a @ xs
ks = ks @ xs
引数 a は実対称行列で、max_iter が繰り返しの最大値、e が 0 と判断する値です。変数 ks に固有ベクトルを格納します。ks は単位行列に初期化します。非対角成分の選択方法ですが、今回は絶対値最大の成分から順に消していくことにします。
a は対称行列なので、対角成分を除いた上三角行列の中から最大値を求めています。argmax() は行列を平坦化して検索するので、返り値は 1 次元での添字になります。これを divmod() で行と列の位置に変換します。
NumPy の関数 triu() は引数の行列を上三角行列に、tril() は下三角行列に変換します。二番目の引数に整数値を指定すると、取り出す三角行列の範囲を変更することができます。簡単な例を示しましょう。
>>> a
array([[1, 1, 1, 1, 1],
[1, 1, 1, 1, 1],
[1, 1, 1, 1, 1],
[1, 1, 1, 1, 1],
[1, 1, 1, 1, 1]])
>>> np.triu(a)
array([[1, 1, 1, 1, 1],
[0, 1, 1, 1, 1],
[0, 0, 1, 1, 1],
[0, 0, 0, 1, 1],
[0, 0, 0, 0, 1]])
>>> np.tril(a)
array([[1, 0, 0, 0, 0],
[1, 1, 0, 0, 0],
[1, 1, 1, 0, 0],
[1, 1, 1, 1, 0],
[1, 1, 1, 1, 1]])
>>> np.triu(a, 1)
array([[0, 1, 1, 1, 1],
[0, 0, 1, 1, 1],
[0, 0, 0, 1, 1],
[0, 0, 0, 0, 1],
[0, 0, 0, 0, 0]])
>>> np.triu(a, -1)
array([[1, 1, 1, 1, 1],
[1, 1, 1, 1, 1],
[0, 1, 1, 1, 1],
[0, 0, 1, 1, 1],
[0, 0, 0, 1, 1]])
最大値 a[x, y] の絶対値が e よりも小さくなったならば、非対角成分がすべて 0 になったと判断して、固有値 (a の対角成分) と固有ベクトル ks を返します。そうでなければ、直交行列 xs を生成して a を相似変換して、ks を ks @ xs の値に更新します。
それでは実際に試してみましょう。収束の様子を見るため、0 にする対角成分の添字と値を表示しています。
>>> a = np.array([[2, 1], [1, 3]])
>>> a
array([[2, 1],
[1, 3]])
>>> jacobi(a)
0 0 1 1 1
1 0 1 -1.11022302463e-16 2.22044604925e-16
(array([ 1.38196601, 3.61803399]), array([[ 0.85065081, 0.52573111],
[-0.52573111, 0.85065081]]))
>>> b = np.array([[2, 1], [1, 2]])
>>> jacobi(b)
0 0 1 1 1
1 0 1 2.22044604925e-16 0.0
(array([ 1., 3.]), array([[ 0.70710678, 0.70710678],
[-0.70710678, 0.70710678]]))
>>> c = np.array(np.array([[1, 4, 5], [4, 2, 6], [5, 6, 3]]))
>>> jacobi(c)
0 1 2 6 6
1 0 2 6.3878493896 6.3878493896
2 0 1 -0.383730028192 -0.383730028192
3 1 2 -0.205185628497 -0.205185628497
4 0 2 0.0776484873658 0.0776484873658
5 0 1 0.00100571570798 0.00100571570798
6 1 2 5.31860490683e-06 5.31860490661e-06
7 0 2 4.60567499583e-09 4.60567485339e-09
8 0 1 -1.54599560073e-15 -1.66005362153e-15
(array([ -2.50728797, -3.6686831 , 12.17597107]), array([[ 0.80958546, 0.31298568, \
0.49659978],
[-0.57735027, 0.57735027, 0.57735027],
[-0.10600965, -0.7541264 , 0.64811675]]))
>>> d = np.ones((4, 4))
>>> d += np.diag([4, 5, 6, 7])
>>> d
array([[ 5., 1., 1., 1.],
[ 1., 6., 1., 1.],
[ 1., 1., 7., 1.],
[ 1., 1., 1., 8.]])
>>> jacobi(d)
0 0 1 1.0 1.0
1 1 2 1.37638192047 1.37638192047
2 2 3 1.65803149623 1.65803149623
3 0 2 0.40114835604 0.40114835604
4 1 3 0.277532107042 0.277532107042
5 0 1 -0.241816045396 -0.241816045396
6 1 2 0.23008016821 0.23008016821
7 2 3 0.0651715920898 0.0651715920898
8 0 2 0.052970083056 0.052970083056
9 0 3 0.00763882210731 0.00763882210731
10 1 3 -0.00541701620336 -0.00541701620336
11 0 1 -0.00274293967704 -0.00274293967704
12 1 2 -0.000133182422716 -0.000133182422717
13 2 3 7.41192264003e-05 7.41192264e-05
14 0 3 1.33206006579e-05 1.3320600658e-05
15 0 2 -5.86729579798e-07 -5.86729579994e-07
16 1 3 -3.1094169516e-08 -3.10941693013e-08
17 0 1 -1.75345826883e-11 -1.75351360697e-11
18 2 3 -3.53381338362e-12 -3.53411855827e-12
19 1 2 -6.992031324e-13 -6.996357785e-13
20 0 3 1.12335090335e-19 1.69864496903e-16
(array([ 4.29608965, 5.39227529, 9.80388636, 6.50774871]),
array([[ 0.90568356, 0.22590297, 0.33200196, -0.13594053],
[-0.38096261, 0.80178162, 0.40111308, -0.22610179],
[-0.15738124, -0.51753551, 0.50656131, -0.67140433],
[-0.09917619, -0.19562996, 0.68722531, 0.69254197]]))
ヤコビ法の場合、相似変換のとき以前 0 にした非対角成分が 0 でなくなることがあります。それでも、最初よりは 0 に近い値になるので、相似変換を繰り返し行っていけば、非対角成分は 0 に近づいていきます。最終的にはすべての非対角成分を 0 にすることができますが、それなりの回数が必要になります。特に、行列 (次元数 N) が大きくなると、収束するまでに必要な回数は大幅に増えてしまいます。このため、ヤコビ法が実用になるのは N が数十程度までといわれています。
#
# eig.py : 固有値と固有ベクトルを求める
#
# Copyright (C) 2018-2023 Makoto Hiroi
#
import numpy as np
#
# 累乗法
#
def power(a, max_iter = 512, eps = 1e-8):
x = np.full(len(a), 0.0)
x[0] = 1.0
k = 0
for i in range(max_iter):
# print(i, k, x)
x1 = a @ x
n = np.linalg.norm(x1)
k1 = (x1 @ x) / (x @ x)
if abs(k1 - k) < eps:
return k1, x1 / n
k = k1
x = x1 / n
def power_all(a, max_iter = 512, eps = 1e-8):
n = len(a)
a = a.astype(np.float_)
xs = np.empty(n) # 固有値
yss = np.empty((n, n)) # 固有ベクトル
for i in range(len(a)):
x, ys = power(a, max_iter, eps)
xs[i] = x
yss[:, i] = ys
zs = ys.reshape((n, 1))
a -= x * zs @ zs.T
return xs, yss
#
# LU 分解による連立方程式の解法
#
def lu_solver(xs, ys):
n = len(xs)
zs = ys.astype(np.float_)
# 前進代入
for i in range(1, n):
zs[i] -= xs[i, :i] @ zs[:i]
# 後退代入
for i in range(n - 1, -1, -1):
zs[i] -= xs[i, i+1:] @ zs[i+1:]
zs[i] /= xs[i, i]
return zs
# ピボット選択
def select_pivot(xs, idx, i):
k = np.abs(xs[i:,i]).argmax() + i
if k != i:
temp = xs[i].copy()
xs[i] = xs[k]
xs[k] = temp
idx[i], idx[k] = idx[k], idx[i]
return -1
return 1
# ピボット選択付き
def lu_pv(xs):
n = len(xs)
zs = xs.astype(np.float_)
idx = list(range(n))
det = 1
for i in range(n):
det *= select_pivot(zs, idx, i)
if zs[i, i] == 0: break
for j in range(i + 1, n):
temp = zs[j, i] / zs[i, i]
zs[j, i+1:] -= temp * zs[i, i+1:]
zs[j, i] = temp
return det * np.diag(zs).prod(), zs, idx
#
# 逆べき乗法
#
def power_inv(a, max_iter = 512, eps = 1e-8):
_, a1, idx = lu_pv(a)
x = np.full(len(a), 0.0)
x[0] = 1.0
k = 0
for i in range(max_iter):
# print(i, k, x)
x1 = lu_solver(a1, x[idx])
n = np.linalg.norm(x1)
k1 = (x1 @ x) / (x @ x)
if abs(k1 - k) < eps:
return 1 / k1, x1 / n
k = k1
x = x1 / n
#
# 逆反復法
#
def power_inv_shift(a, d, max_iter = 512, eps = 1e-8):
k, xs = power_inv(a - d * np.eye(len(a)), max_iter, eps)
return k + d, xs
#
# 原点移動
#
def power_shift(a, d, max_iter = 512, eps = 1e-8):
k, xs = power(a - d * np.eye(len(a)), max_iter, eps)
return k + d, xs
#
# ヤコビ法
#
def jacobi(a, max_iter = 512, e = 1e-14):
n = len(a)
ks = np.eye(n) # 固有ベクトル
for i in range(max_iter):
# 絶対値最大の非対角成分を探す
x, y = divmod(np.fabs(np.triu(a, 1)).argmax(), n)
print(i, x, y, a[x, y], a[y, x])
if np.fabs(a[x, y]) < e:
return np.diag(a), ks
# 角度を求める
d = a[y, y] - a[x, x]
if d == 0.0:
r = np.pi / 4
else:
r = 0.5 * np.arctan(2 * a[x, y] / d)
# 直交行列の生成
xs = np.eye(n)
xs[x, x] = np.cos(r)
xs[y, y] = np.cos(r)
xs[x, y] = np.sin(r)
xs[y, x] = - np.sin(r)
#
a = xs.T @ a @ xs
ks = ks @ xs