前回は固有値と固有ベクトルを求める簡単な方法として「累乗法」と「ヤコビ法」を紹介しました。今回は行列の固有値を求める方法の一つである「QR 法 (QR algorithm)」と、その基礎になる「QR 分解 (QR decomposition)」について取り上げます。
QR 法の基本的な考え方はシンプルですが、これをそのままプログラムするのは効率的ではありません。このため、対称行列であれば「三重対角行列」に、それ以外の行列は「ヘッセンベルグ行列」に相似変換してから QR 法を適用するのが主流のようです。本ページでは「実対称行列」の固有値を求める QR 法の基本について簡単に説明します。
まずは最初に QR 分解を説明しましょう。QR 分解は行列 A を直交行列 Q と上三角行列 R に分解することです。
A = QR
基本的な考え方は簡単です。A の左側から適当な直交行列 Qi をいくつか掛け算したら上三角行列 R になったとしましょう。これを式で示すと次のようになります。
Qn...Q2Q1A = R
直交行列の逆行列は転置行列なので、上の式は次のように変換できます。
A = Q1TQ2T...QnTR
直交行列の掛け算は直交行列になるので、行列 A は直交行列 Q と上三角行列 R に分解することができました。
QR 分解は「ギブンス回転」、「ハウスホルダー変換 (Householder transformation)」、「グラム・シュミット分解」などの手法を使って行うことができます。ここでは実装が簡単なギブンス回転を使って説明します。ギブンス回転の場合、直交行列 Qi を掛け算するとき、対角成分よりも下の成分の一つ aij (i < j) が 0 になるように角度を決めます。次の図を見てください。
Q @ A = [[ cos(r), sin(r)], @ [[a, b], = [[..., ...]
[-sin(r), cos(r)]] [c, d]] [c*cos(r) - a*sin(r), ...]]
c * cos(r) - a * sin(r) = 0
c / a = sin(r) / cos(r) = tan(r), r = arctan(c / a) になるが、
ここで e = √(a2 + c2) とすると
sin(r) = c / e, cos(r) = a / e
正方行列 A が二次元行列 [[a, b], [c, d]] の場合、左下の要素 c を 0 にすれば、A は上三角行列になります。この場合、条件は tan(r) = c / a になりますが、三辺の直角三角形 a, c, e (= √(a*a + c* c)) を考えると、sin(r) の値は c / e に、cos(r) の値は a / e になります。
それでは実際に試してみましょう。
>>> def givens(a, c):
... d = np.sqrt(a * a + c * c)
... return np.array([[a / d, c /d], [-c / d, a / d]])
...
>>> a = np.array([[2, 1], [1, 3]])
>>> a
array([[2, 1],
[1, 3]])
>>> q = givens(2, 1)
>>> q
array([[ 0.89442719, 0.4472136 ],
[-0.4472136 , 0.89442719]])
>>> q @ a
array([[ 2.23606798, 2.23606798],
[ 0. , 2.23606798]])
確かにギブンス回転を使って正方行列を上三角行列に変換することができました。
それではギブンス回転を使って正方行列を QR 分解するプログラムを作りましょう。次のリストを見てください。
リスト : QR 分解 (ギブンス回転, 効率は悪い)
def qr_g(a):
n = len(a)
qs = np.eye(n) # 直交行列
for y in range(n):
for x in range(y + 1, n):
# print(x, y, a[x, y])
d = np.sqrt(a[x, y] ** 2 + a[y, y] ** 2)
c = a[y, y] / d
s = a[x, y] / d
# 直交行列の生成
xs = np.eye(n)
xs[x, x] = xs[y, y] = c
xs[x, y] = -s
xs[y, x] = s
#
a = xs @ a
qs = qs @ xs.T
return qs, a
直交行列は変数 qs に保持します。あとは二重の for ループの中で、直交行列 xs を作って対角線よりも下の成分を 0 にしていくだけです。
簡単な実行例を示します。ご参考までに、0 にする成分の座標と値を表示しています。
>>> a
array([[2, 1],
[1, 3]])
>>> x, y = qr_g(a)
1 0 1
>>> x
array([[ 0.89442719, -0.4472136 ],
[ 0.4472136 , 0.89442719]])
>>> y
array([[ 2.23606798, 2.23606798],
[ 0. , 2.23606798]])
>>> x @ y
array([[ 2., 1.],
[ 1., 3.]])
>>> b
array([[2, 1],
[1, 2]])
>>> x, y = qr_g(b)
1 0 1
>>> x
array([[ 0.89442719, -0.4472136 ],
[ 0.4472136 , 0.89442719]])
>>> y
array([[ 2.23606798, 1.78885438],
[ 0. , 1.34164079]])
>>> x @ y
array([[ 2., 1.],
[ 1., 2.]])
>>> c
array([[1, 4, 5],
[4, 2, 6],
[5, 6, 3]])
>>> x, y = qr_g(c)
1 0 4
2 0 5.0
2 1 1.57181049599
>>> x
array([[ 0.15430335, 0.80178373, 0.57735027],
[ 0.6172134 , -0.53452248, 0.57735027],
[ 0.77151675, 0.26726124, -0.57735027]])
>>> y
array([[ 6.48074070e+00, 6.48074070e+00, 6.78934740e+00],
[ 1.86554783e-16, 3.74165739e+00, 1.60356745e+00],
[ -4.03004391e-16, 0.00000000e+00, 4.61880215e+00]])
>>> x @ y
array([[ 1., 4., 5.],
[ 4., 2., 6.],
[ 5., 6., 3.]])
>>> d
array([[ 6., 1., 1., 1.],
[ 1., 7., 1., 1.],
[ 1., 1., 8., 1.],
[ 1., 1., 1., 9.]])
>>> x, y = qr_g(d)
1 0 1.0
2 0 1.0
3 0 1.0
2 1 0.640056896475
3 1 0.623429200159
3 2 0.453099524675
>>> x
array([[ 0.96076892, -0.19232689, -0.15176736, -0.13000043],
[ 0.16012815, 0.9729478 , -0.1264728 , -0.1083337 ],
[ 0.16012815, 0.09050677, 0.97854228, -0.09285745],
[ 0.16012815, 0.09050677, 0.05853467, 0.98119376]])
>>> y
array([[ 6.24499800e+00, 2.40192231e+00, 2.56205046e+00, 2.72217861e+00],
[ 1.04511300e-17, 6.79932123e+00, 1.59518185e+00, 1.68568863e+00],
[ 1.10271650e-16, -1.10825269e-16, 7.60863275e+00, 1.22711418e+00],
[ -7.54245832e-18, 6.61145756e-18, 0.00000000e+00, 8.49955224e+00]])
>>> x @ y
array([[ 6., 1., 1., 1.],
[ 1., 7., 1., 1.],
[ 1., 1., 8., 1.],
[ 1., 1., 1., 9.]])
>>> e
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.]])
>>> x, y = qr_g(e)
1 0 1.0
2 0 1.0
3 0 1.0
4 0 1.0
2 1 0.693103280084
3 1 0.679643682708
4 1 0.66669690304
3 2 0.517096952939
4 2 0.50351829074
4 3 0.399484471853
>>> x
array([[ 0.96152395, -0.17508462, -0.13927232, -0.1189625 , -0.10615186],
[ 0.13736056, 0.97375831, -0.11937627, -0.10196786, -0.09098731],
[ 0.13736056, 0.08394468, 0.97968691, -0.08922188, -0.0796139 ],
[ 0.13736056, 0.08394468, 0.0572978 , 0.98274831, -0.07076791],
[ 0.13736056, 0.08394468, 0.0572978 , 0.04117894, 0.98443213]])
>>> y
array([[ 7.28010989e+00, 2.47249015e+00, 2.60985072e+00, 2.74721128e+00, 2.88457184e+00],
[ -1.09771331e-16, 7.86681590e+00, 1.72206519e+00, 1.80600987e+00, 1.88995455e+00],
[ 1.09335210e-17, 6.44539686e-18, 8.67312924e+00, 1.35131411e+00, 1.40861191e+00],
[ 9.33910629e-18, 4.63219507e-18, 0.00000000e+00, 9.55850975e+00, 1.12556437e+00],
[ 8.33341176e-18, 1.10738211e-16, 0.00000000e+00, 0.00000000e+00, 1.04812325e+01]])
>>> x @ y
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.]])
正常に動作していますね。なお、0 の精度が足りない (もっと 0 に近づけたい) 場合は、上三角行列に qr_g() を繰り返し適用すれば OK です。
>>> d
array([[ 6., 1., 1., 1.],
[ 1., 7., 1., 1.],
[ 1., 1., 8., 1.],
[ 1., 1., 1., 9.]])
>>> q, r = qr_g(d)
>>> q1, r1 = qr_g(r)
>>> q1
array([[ 1.00000000e+00, -1.67352016e-18, -1.76575958e-17, 1.20775993e-18],
[ 1.67352016e-18, 1.00000000e+00, 2.25371676e-17, -1.39902245e-18],
[ 1.76575958e-17, -2.25371676e-17, 1.00000000e+00, -1.13377354e-19],
[ -1.20775993e-18, 1.39902245e-18, 1.13377354e-19, 1.00000000e+00]])
>>> r1
array([[ 6.24499800e+00, 2.40192231e+00, 2.56205046e+00, 2.72217861e+00],
[ 2.77792038e-49, 6.79932123e+00, 1.59518185e+00, 1.68568863e+00],
[ -1.23259516e-32, 0.00000000e+00, 7.60863275e+00, 1.22711418e+00],
[ 1.39748378e-51, 0.00000000e+00, 0.00000000e+00, 8.49955224e+00]])
>>> q @ q1 @ r1
array([[ 6., 1., 1., 1.],
[ 1., 7., 1., 1.],
[ 1., 1., 8., 1.],
[ 1., 1., 1., 9.]])
実用的には、ギブンス回転よりも効率が良い「ハウスホルダー変換」を用いて QR 分解を行うのが主流のようです。ハウスホルダー変換は次回以降に取り上げたいと思います。なお、NumPy には行列を QR 分解する関数 linalg.qr() が用意されているので、私たちが作る必要はありません。簡単な使用例を示します。
>>> c
array([[1, 4, 5],
[4, 2, 6],
[5, 6, 3]])
>>> q, r = np.linalg.qr(c)
>>> q
array([[-0.15430335, 0.80178373, -0.57735027],
[-0.6172134 , -0.53452248, -0.57735027],
[-0.77151675, 0.26726124, 0.57735027]])
>>> r
array([[-6.4807407 , -6.4807407 , -6.7893474 ],
[ 0. , 3.74165739, 1.60356745],
[ 0. , 0. , -4.61880215]])
>>> q @ r
array([[ 1., 4., 5.],
[ 4., 2., 6.],
[ 5., 6., 3.]])
>>> d
array([[ 6., 1., 1., 1.],
[ 1., 7., 1., 1.],
[ 1., 1., 8., 1.],
[ 1., 1., 1., 9.]])
>>> q, r = np.linalg.qr(d)
>>> q
array([[-0.96076892, 0.19232689, 0.15176736, -0.13000043],
[-0.16012815, -0.9729478 , 0.1264728 , -0.1083337 ],
[-0.16012815, -0.09050677, -0.97854228, -0.09285745],
[-0.16012815, -0.09050677, -0.05853467, 0.98119376]])
>>> r
array([[-6.244998 , -2.40192231, -2.56205046, -2.72217861],
[ 0. , -6.79932123, -1.59518185, -1.68568863],
[ 0. , 0. , -7.60863275, -1.22711418],
[ 0. , 0. , 0. , 8.49955224]])
>>> q @ r
array([[ 6., 1., 1., 1.],
[ 1., 7., 1., 1.],
[ 1., 1., 8., 1.],
[ 1., 1., 1., 9.]])
>>> e
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.]])
>>> q, r = np.linalg.qr(e)
>>> q
array([[-0.96152395, 0.17508462, 0.13927232, 0.1189625 , -0.10615186],
[-0.13736056, -0.97375831, 0.11937627, 0.10196786, -0.09098731],
[-0.13736056, -0.08394468, -0.97968691, 0.08922188, -0.0796139 ],
[-0.13736056, -0.08394468, -0.0572978 , -0.98274831, -0.07076791],
[-0.13736056, -0.08394468, -0.0572978 , -0.04117894, 0.98443213]])
>>> r
array([[ -7.28010989, -2.47249015, -2.60985072, -2.74721128, -2.88457184],
[ 0. , -7.8668159 , -1.72206519, -1.80600987, -1.88995455],
[ 0. , 0. , -8.67312924, -1.35131411, -1.40861191],
[ 0. , 0. , 0. , -9.55850975, -1.12556437],
[ 0. , 0. , 0. , 0. , 10.48123249]])
>>> q @ r
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.]])
ところで、行列 A が正則行列であれば、QR 分解を使って連立一次方程式を解くことができます。次の式を見てください。
Ax = b => (QR)x = b => Qy = b, ただし y = Rx
Q は直交行列なので、y の値は QTb で求めることができます。R は上三角行列なので、ガウスの消去法と同様に後退代入で x を求めることができます。ガウスの消去法はピボット選択が必要になる場合がありますが、QR 分解だとその必要がないところが利点です。
プログラムと実行例を示します。
リスト : 連立一次方程式の解法
def qr_solver(a, b):
q, r = qr_g(a)
x = q.T @ b
# 後退代入
for i in range(len(a) - 1, -1, -1):
x[i] -= r[i, i+1:] @ x[i+1:]
x[i] /= r[i, i]
return x
>>> qr_solver(np.array([[1, 1], [2, 4]]), np.array([100, 272])) array([ 64., 36.]) >>> qr_solver(np.array([[1, 1, 1], [2, 4, 6], [2, 0, 4]]), np.array([10, 38, 14])) array([ 3., 5., 2.]) >>> qr_solver(np.array([[0, 2, 4], [1, 1, 1], [4, 2, 6]]), np.array([14, 10, 38])) array([ 5., 3., 2.]) >>> qr_solver(np.array([[1, 1, 1, 1], [-1, 1, -1, 1], [8, 4, 2, 1], [-8, 4, -2, 1]]), np.array([-5, -7, -31, -35])) array([ -3.89492003e-16, -9.00000000e+00, 1.00000000e+00, 3.00000000e+00]) >>> qr_solver(np.array([[1, -1, 1, -1, 1], [12, -6, 2, 0, 0], [1, 1, 1, 1, 1], [12, 6, 2, 0, 0], [4, 3, 2, 1, 0]]), np.array([1, 0, 8, 0, 1])) array([ 3.12500000e-01, -1.46651700e-16, -1.87500000e+00, 3.50000000e+00, 6.06250000e+00])
次は QR 法の基本を簡単に説明します。上 (下) 三角行列の固有値は、対角行列と同様に対角成分と一致します。QR 分解は相似変換ではないので、行列 A を QR 分解したとしても、R の対角成分は A の固有値とは一致しません。そこで、次に示すような相似変換を考えます。
Ai+1 = QiTAiQi
ここで Ai を QR 分解して、Ai = QiRi とおけば、上式は次のようになります。
Ai+1 = QiTAiQi => Ai+1 = QiTQiRiQi => Ai+1 = RiQi
つまり、Ai を QR 分解して Qi と Ri を求め、それを逆に掛け算したものが Ai+1 になります。これを繰り返すと Ai+1 の対角成分より下側の成分は 0 に収束することが知られています。つまり、Ai+1 は上三角行列になるので、その対角成分が固有値になります。これは相似変換なので、Ai+1 の固有値は A の固有値と一致します。ただし、ヤコビ法と違って直交行列が固有ベクトルと一致するとは限りません。固有ベクトルを求めるときは他の方法 (たとえば逆反復法など) を使用することになります。
これをそのままプログラムすると次のようになります。
リスト : QR 法 (ナイーブ版)
def qr_eig(a, max_iter = 512, eps = 1e-8):
d = np.diag(a)
for i in range(max_iter):
# print(i, d)
q, r = qr_g(a)
a1 = r @ q
d1 = np.diag(a1)
if np.fabs((d1 - d) / d1).sum() < eps:
return d1
a, d = a1, d1
対角成分 d1 と d の差分が許容誤差 eps 内に入ったら、収束したとして対角成分 (固有値) d1 を返します。あとは、QR 法の原理をそのままプログラムしただけなので、特に難しいところはないと思います。
それでは実行してみましょう。ご参考までに、収束の様子を表示しています。
>>> a
array([[2, 1],
[1, 3]])
>>> qr_eig(a)
0 [2 3]
1 [ 3. 2.]
2 [ 3.5 1.5]
3 [ 3.6 1.4]
4 [ 3.61538462 1.38461538]
5 [ 3.61764706 1.38235294]
6 [ 3.61797753 1.38202247]
7 [ 3.61802575 1.38197425]
8 [ 3.61803279 1.38196721]
9 [ 3.61803381 1.38196619]
10 [ 3.61803396 1.38196604]
11 [ 3.61803399 1.38196601]
array([ 3.61803399, 1.38196601])
>>> b
array([[2, 1],
[1, 2]])
>>> qr_eig(b)
0 [2 2]
1 [ 2.8 1.2]
2 [ 2.97560976 1.02439024]
3 [ 2.99726027 1.00273973]
4 [ 2.99969521 1.00030479]
5 [ 2.99996613 1.00003387]
6 [ 2.99999624 1.00000376]
7 [ 2.99999958 1.00000042]
8 [ 2.99999995 1.00000005]
9 [ 2.99999999 1.00000001]
array([ 3., 1.])
>>> c
array([[1, 4, 5],
[4, 2, 6],
[5, 6, 3]])
>>> qr_eig(c)
0 [1 2 3]
1 [ 10.23809524 -1.57142857 -2.66666667]
2 [ 12.05490483 -2.5701002 -3.48480463]
3 [ 12.16829023 -2.6954938 -3.47279643]
・・・省略・・・
27 [ 12.17597107 -3.66868303 -2.50728804]
28 [ 12.17597107 -3.66868307 -2.507288 ]
29 [ 12.17597107 -3.66868308 -2.50728798]
array([ 12.17597107, -3.66868309, -2.50728797])
>>> d
array([[ 6., 1., 1., 1.],
[ 1., 7., 1., 1.],
[ 1., 1., 8., 1.],
[ 1., 1., 1., 9.]])
>>> qr_eig(d)
0 [ 6. 7. 8. 9.]
1 [ 7.23076923 6.91232561 7.51719756 8.3397076 ]
2 [ 9.00392842 6.26755685 6.94463476 7.78387997]
3 [ 10.17466019 5.81931268 6.65717678 7.34885036]
・・・省略・・・
62 [ 10.80388636 7.50774853 6.39227547 5.29608965]
63 [ 10.80388636 7.50774858 6.39227542 5.29608965]
64 [ 10.80388636 7.50774861 6.39227538 5.29608965]
array([ 10.80388636, 7.50774864, 6.39227536, 5.29608965])
>>> e
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.]])
>>> qr_eig(e)
0 [ 7. 8. 9. 10. 11.]
1 [ 8.47169811 8.11519213 8.65508891 9.4399588 10.31806205]
2 [ 10.88373702 7.45513421 8.03552924 8.87518195 9.75041758]
3 [ 12.55274895 6.86985512 7.72900506 8.57893445 9.26945642]
・・・省略・・・
73 [ 13.39054123 9.54039425 8.43473658 7.35663212 6.27769582]
74 [ 13.39054123 9.54039429 8.4347366 7.35663206 6.27769582]
75 [ 13.39054123 9.54039432 8.43473662 7.35663201 6.27769582]
array([ 13.39054123, 9.54039434, 8.43473663, 7.35663197, 6.27769582])
ナイーブな QR 法の場合、固有値は対角線上に絶対値の大きい順に並びます。結果を見ればお分かりのように、ナイーブな QR 法は収束回数が多くなるという欠点があります。これを改良する方法に「原点移動法」と「減次 (デフレーション)」があります。
QR 法の原点移動法は次式のように行います。
Ai - μiI = QiRi Ai+1 = RiQi + μiI
Ai' = Ai - μiI = QiRi とすると
Ai+1' = RiQi
= QiTAi'Qi
= QiT(Ai - μiI)Qi
= QiTAiQi - μiQiTIQi
= Ai+1 - μiI
よって Ai+1 = Ai+1' + μiI = RiQi + μiI
Ai の対角成分から移動量 μi を引き算して QR 分解を行い、Ai+1 を求めるとき Ri @ Qi の対角成分に μi を足し算します。移動量 μ は、右下隅 2 行 2 列の行列の固有値で、右下の要素に近い値を選びます。これを「ウィルキンソンの移動法」といいます。
プログラムは簡単に修正できるので、さっそく試してみたのですが、ウィルキンソンの移動法だけでは QR 法の収束回数を大幅に減らすことはできませんでした。右下隅の対角成分の収束は速くなるのですが、それ以外の対角成分の収束は速くなりません。
そこで「減次 (デフレーション)」という方法と組み合わせることにします。デフレーションは n * n の行列 A の固有値 λn をひとつ求めたら、行列の次数を一つ減らして n-1 * n-1 の行列 A' で固有値 λn-1 を求めるという方法です。
相似変換で n * n の行列 A が上三角行列に近づいていくとき、対角成分は固有値に近づいていきます。一番下の行の成分 an,0, ..., an,n-1 が 0 に収束すると、対角成分 an,n は固有値 λn に収束します。他の非対角成分が固有値 λn に影響を与えることはありません。残りの固有値 λ1, ..., λn-1 は n-1 * n-1 の行列 A' を上三角行列に相似変換すれば求めることができるはずです。
そしてここがポイントですが、最初にウィルキンソンの移動法を適用し、そのあとデフレーションしたときも、その行列での右下隅 2 * 2 行列の固有値を求めてウィルキンソンの移動法を適用します。これで右下隅の対角成分 an-1,n-1 の収束が速くなります。
一番下の行の成分 an-1,0, ..., an-1,n-2 が 0 に収束したら、対角成分 an-1,n-1 の値が固有値 λn-1 になります。あとはこれを繰り返すことで、すべての固有値を高速に求めることができます。
それではプログラムを作りましょう。最初に行列 [[a, b], [c, d]] の固有値を求め、d に近い値を返す関数 eig22() を作ります。これは二次方程式を解くだけなので簡単です。
|[[a, b], [c, d]] - λI| = 0 => (a - λ)(d - λ) - bc = 0 => λ2 - (a + d)λ + ad - bc = 0 λ= (a + d ±√((a + d)2 - 4(ad - bc))) / 2
リスト : シフト値を求める
def eig22(a, b, c, d):
e = a + d
f = np.sqrt(e * e - 4 * (a * d - b * c))
k1 = (e + f) / 2
k2 = (e - f) / 2
if np.fabs(d - k1) < np.fabs(d - k2):
return k1
else:
return k2
次はウィルキンソンの移動法とデフレーションを行う関数 qr_eig_shiftd() を作ります。
リスト : ウィルキンソンの移動法とデフレーション
def qr_eig_shiftd(a, max_iter = 1024, e = 1e-14):
ks = []
i = 0
for size in range(len(a), 1, -1):
a = a[:size, :size]
while i < max_iter:
# シフト値の計算
k = eig22(*a[-2:, -2:].flat)
dk = np.diag(np.full(size, k))
q, r = qr_g(a - dk)
a = r @ q + dk
# print(i, size, k)
# print(a)
i += 1
if np.all(np.fabs(a[-1, :-1]) < e):
ks.append(a[-1, -1])
break
else:
raise Exception("Repeated Over!!", i)
ks.append(a[0, 0])
return np.array(ks)
引数 a が固有値を求める行列、max_iter が繰り返しの最大値、e が 0 とみなす値です。変数 ks は固有値を格納するリストです。最初の for ループで、行列 a の次元数 (size) を一つずつ減らしていきます。NumPy の場合、a の大きさは a[:size, :size] で簡単に変更することができます。次の while ループで相似変換を行います。まず、シフト値を関数 eig22() で求めて変数 k にセットし、変数 dk に対角成分が k の対角行列をセットします。これは k * np.eye(size) としてもかまいません。それから、gr_g() で行列 a - dk を QR 分解して、r @ q + dk で相似変換します。
次の if 分で、最終行の対角成分以外の要素の絶対値が e 未満なったら収束と判定し、固有値 a[-1, -1] を ks に追加してから break で while ループを脱出します。なお、a[-1, -1] の左隣の要素 a[-1, -2] だけで収束を判定する方法もあるようです。最初の for ループを抜けたとき、行列 a は 2 行 2 列になっています。a[1, 1] の固有値は ks に追加していますが、a[0, 0] の固有値は追加していません。append() で ks に a[0, 0] を追加してから return で ks を返します。
プログラムはこれで完成です。簡単な実行例を示します。ご参考までに、収束の様子を表示しています。
>>> a
array([[2, 1],
[1, 3]])
>>> qr_eig_shiftd(a)
0 2 3.61803398875
[[ 1.38196601e+00 0.00000000e+00]
[ 5.83678785e-17 3.61803399e+00]]
array([ 3.61803399, 1.38196601])
>>> b
array([[2, 1],
[1, 2]])
>>> qr_eig_shiftd(b)
0 2 1.0
[[ 3. 0.]
[ 0. 1.]]
array([ 1., 3.])
>>> c
array([[1, 4, 5],
[4, 2, 6],
[5, 6, 3]])
>>> qr_eig_shiftd(c)
0 3 8.5207972894
[[ -2.18468389 0.36261048 2.69196612]
[ 0.36261048 -2.23784521 4.17335334]
[ 2.69196612 4.17335334 10.4225291 ]]
1 3 11.6744358241
[[ -2.71116708 -0.44141708 0.10137603]
[ -0.44141708 -3.46268926 0.1521935 ]
[ 0.10137603 0.1521935 12.17385634]]
2 3 12.1753375269
[[ -2.73873218e+00 -4.63930768e-01 4.31264245e-06]
[ -4.63930768e-01 -3.43723889e+00 5.92369306e-06]
[ 4.31264245e-06 5.92369306e-06 1.21759711e+01]]
3 3 12.175971065
[[ -2.76824116e+00 -4.84740342e-01 -6.07480171e-16]
[ -4.84740342e-01 -3.40772990e+00 1.87499213e-17]
[ 3.32874030e-19 4.17057524e-19 1.21759711e+01]]
4 2 -3.66868309795
[[ -2.50728797e+00 4.99600361e-16]
[ 3.02600320e-16 -3.66868310e+00]]
array([ 12.17597107, -3.6686831 , -2.50728797])
>>> d
array([[ 6., 1., 1., 1.],
[ 1., 7., 1., 1.],
[ 1., 1., 8., 1.],
[ 1., 1., 1., 9.]])
>>> qr_eig_shiftd(d)
0 4 9.61803398875
[[ 5.39663758 0.32266086 0.32231652 -0.37256625]
[ 0.32266086 6.55308941 0.55249916 -0.63863477]
[ 0.32231652 0.55249916 8.08193777 -1.25061382]
[-0.37256625 -0.63863477 -1.25061382 9.96833524]]
1 4 10.591552409
[[ 5.32345151 0.17018919 0.06893487 0.01745754]
[ 0.17018919 6.39147123 0.15856483 0.04015604]
[ 0.06893487 0.15856483 7.48633156 0.12316191]
[ 0.01745754 0.04015604 0.12316191 10.7987457 ]]
2 4 10.8033187827
[[ 5.31261672e+00 1.33156028e-01 4.12195538e-02 -1.80853875e-06]
[ 1.33156028e-01 6.39022484e+00 1.20797339e-01 -5.30007357e-06]
[ 4.12195538e-02 1.20797339e-01 7.49327208e+00 -2.16426813e-05]
[ -1.80853875e-06 -5.30007357e-06 -2.16426813e-05 1.08038864e+01]]
3 4 10.803886359
[[ 5.30634471e+00 1.05318869e-01 2.48048781e-02 1.47853836e-17]
[ 1.05318869e-01 6.39008321e+00 9.18730579e-02 -4.24549007e-16]
[ 2.48048781e-02 9.18730579e-02 7.49968572e+00 3.91064947e-16]
[ 3.07495988e-18 1.13882502e-17 6.19396968e-17 1.08038864e+01]]
4 3 7.50724119488
[[ 5.29858409e+00 5.22316899e-02 5.73353360e-06]
[ 5.22316899e-02 6.38978085e+00 4.27866991e-05]
[ 5.73353360e-06 4.27866991e-05 7.50774870e+00]]
5 3 7.50774870534
[[ 5.29672526e+00 2.63883732e-02 -3.96894228e-17]
[ 2.63883732e-02 6.39163968e+00 1.04189796e-15]
[ 7.02984464e-17 1.04333854e-15 7.50774871e+00]]
6 2 6.39227529027
[[ 5.29608965e+00 4.19803081e-16]
[ 1.38274344e-16 6.39227529e+00]]
array([ 10.80388636, 7.50774871, 6.39227529, 5.29608965])
>>> e
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.]])
>>> qr_eig_shiftd(e)
0 5 11.6180339887
[[ 6.40992926 0.36217332 0.36378603 0.43255257 0.39145839]
[ 0.36217332 7.54973524 0.55218313 0.65656241 0.59418642]
[ 0.36378603 0.55218313 8.8970441 1.06661253 0.96528017]
[ 0.43255257 0.65656241 1.06661253 10.71437352 1.55150134]
[ 0.39145839 0.59418642 0.96528017 1.55150134 11.42891788]]
1 5 12.6637511626
[[ 6.31313998 0.18946853 0.10252768 0.05328719 0.06236709]
[ 0.18946853 7.3679007 0.19908322 0.10347045 0.12110136]
[ 0.10252768 0.19908322 8.44432735 0.23093232 0.27028217]
[ 0.05328719 0.10347045 0.23093232 9.68149993 0.79762451]
[ 0.06236709 0.12110136 0.27028217 0.79762451 13.19313204]]
2 5 13.365811482
[[ 6.30050796e+00 1.53776217e-01 6.65952049e-02 2.10645413e-02 2.24358222e-04]
[ 1.53776217e-01 7.36024741e+00 1.56010797e-01 4.93473347e-02 5.25597978e-04]
[ 6.65952049e-02 1.56010797e-01 8.42724096e+00 1.35139381e-01 1.43936823e-03]
[ 2.10645413e-02 4.93473347e-02 1.35139381e-01 9.52147103e+00 5.55418280e-03]
[ 2.24358222e-04 5.25597978e-04 1.43936823e-03 5.55418280e-03 1.33905326e+01]]
3 5 13.3905406131
[[ 6.29325740e+00 1.27774525e-01 4.58363269e-02 1.14741331e-02 1.96130417e-11]
[ 1.27774525e-01 7.35817710e+00 1.28488231e-01 3.21642497e-02 5.49800693e-11]
[ 4.58363269e-02 1.28488231e-01 8.42903296e+00 1.07399122e-01 1.83582374e-10]
[ 1.14741331e-02 3.21642497e-02 1.07399122e-01 9.52899131e+00 9.04227484e-10]
[ 1.96132353e-11 5.49797526e-11 1.83581995e-10 9.04227877e-10 1.33905412e+01]]
4 5 13.390541233
[[ 6.28850597e+00 1.06946506e-01 3.17050451e-02 6.23606749e-03 2.01056861e-16]
[ 1.06946506e-01 7.35706076e+00 1.05853177e-01 2.08202685e-02 -3.07058311e-16]
[ 3.17050451e-02 1.05853177e-01 8.43036851e+00 8.46492111e-02 -3.92739717e-16]
[ 6.23606749e-03 2.08202685e-02 8.46492111e-02 9.53352353e+00 3.82980461e-16]
[ 5.07485169e-27 1.57618168e-26 6.88738274e-26 3.16578380e-25 1.33905412e+01]]
5 4 9.53998118262
[[ 6.28228110e+00 7.00964507e-02 1.06139732e-02 -7.94557973e-07]
[ 7.00964507e-02 7.35481196e+00 5.37254679e-02 -4.02186798e-06]
[ 1.06139732e-02 5.37254679e-02 8.43197129e+00 -3.23372054e-05]
[ -7.94557973e-07 -4.02186798e-06 -3.23372054e-05 9.54039442e+00]]
6 4 9.54039442567
[[ 6.27972377e+00 4.67164901e-02 3.60448231e-03 -6.65142151e-17]
[ 4.67164901e-02 7.35531082e+00 2.74145505e-02 2.76966221e-16]
[ 3.60448231e-03 2.74145505e-02 8.43402975e+00 3.93199470e-16]
[ 3.61289724e-18 2.74785518e-17 4.35043023e-16 9.54039443e+00]]
7 3 8.43472601226
[[ 6.27820046e+00 2.33284107e-02 1.78220228e-08]
[ 2.33284107e-02 7.35612722e+00 2.72067721e-07]
[ 1.78220231e-08 2.72067721e-07 8.43473667e+00]]
8 3 8.4347366665
[[ 6.27782193e+00 1.16638031e-02 2.77970468e-16]
[ 1.16638031e-02 7.35650575e+00 -3.43316800e-16]
[ 2.52505622e-23 7.98560634e-22 8.43473667e+00]]
9 2 7.35663185484
[[ 6.27769582e+00 -4.16333634e-17]
[ 2.44075076e-17 7.35663185e+00]]
array([ 13.39054123, 9.54039443, 8.43473667, 7.35663185, 6.27769582])
結果を見ればお分かりのように、ウィルキンソンの移動法とデフレーションを適用することで、QR 法の収束回数を大幅に減少させることができました。ここまで効果が高いとは M.Hiroi も大変驚きました。なお、実用的には行列 A をそのまま QR 分解するのではなく「三重対角行列」に相似変換してから QR 法を適用したほうが良いようです。これは次回以降に取り上げたいと思います。
#
# qr.py : QR 分解と QR 法
#
# Copyright (C) 2018-2023 Makoto Hiroi
#
import numpy as np
#
# ギブンス回転による QR 分解
#
def qr_g(a):
n = len(a)
qs = np.eye(n) # 直交行列
for y in range(n):
for x in range(y + 1, n):
# print(x, y, a[x, y])
d = np.sqrt(a[x, y] ** 2 + a[y, y] ** 2)
c = a[y, y] / d
s = a[x, y] / d
xs = np.eye(n)
xs[x, x] = xs[y, y] = c
xs[x, y] = -s
xs[y, x] = s
a = xs @ a
#print(a)
qs = qs @ xs.T
return qs, a
#
# 連立一次方程式の解法
#
def qr_solver(a, b):
n = len(a)
q, r = qr_g(a)
x = q.T @ b
# 後退代入
for i in range(n - 1, -1, -1):
x[i] -= r[i, i+1:] @ x[i+1:]
x[i] /= r[i, i]
return x
#
# QR 法 (単純な方法, 収束回数が多い)
#
def qr_eig(a, max_iter = 512, eps = 1e-8):
d = np.diag(a)
for i in range(max_iter):
print(i, d)
q, r = qr_g(a)
a1 = r @ q
d1 = np.diag(a1)
if np.fabs((d1 - d) / d1).sum() < eps:
return d1
a, d = a1, d1
#
# QR 法 (ウィルキンソンの移動法とデフレーション)
#
# 行列 [[a, b], [c, d]] の固有値を求める
def eig22(a, b, c, d):
e = a + d
f = np.sqrt(e * e - 4 * (a * d - b * c))
k1 = (e + f) / 2
k2 = (e - f) / 2
if np.fabs(d - k1) < np.fabs(d - k2):
return k1
else:
return k2
# ウィルキンソンの移動法とデフレーション
def qr_eig_shiftd(a, max_iter = 1024, e = 1e-14):
ks = []
i = 0
for size in range(len(a), 1, -1):
a = a[:size, :size]
while i < max_iter:
# シフト値の計算
k = eig22(*a[-2:, -2:].flat)
dk = np.diag(np.full(size, k))
q, r = qr_g(a - dk)
a = r @ q + dk
# print(i, size, k)
#print(a)
i += 1
if np.all(np.fabs(a[-1, :-1]) < e):
ks.append(a[-1, -1])
break
else:
raise Exception("Repeated Over!!", i)
ks.append(a[0, 0])
return np.array(ks)