固有値 \(\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\) としたとき、\(|\frac{\lambda_2}{\lambda_1}| \gt |\frac{\lambda_2 - p}{\lambda_1 - p}|\) になるならば、行列 \(B\) に累乗法を適用したほうが収束が速くなるはずです。これを「原点移動 (シフト法)」といいます。
ただし、参考文献によっては、逆反復法のことをシフト法と呼ぶことがあるようです。まあ、どちらの方法も対角成分の値をシフトするので、シフト法というと混乱するような気がします。本稿では前者を逆反復法、後者を原点移動と呼ぶことにします。
プログラムは簡単です。次のリストを見てください。
リスト : 逆反復法と原点移動
# 逆反復法
function power_inv_shift(a::Matrix{Float64}, d::Float64; max_iter = 512, eps = 1e-8, tiny = 1e-10)
k, xs = power_inv(a - d * diagm(0 => ones(size(a, 1))), max_iter=max_iter, eps=eps, tiny=tiny)
k + d, xs
end
# 原点移動
function power_shift(a::Matrix{Float64}, d::Float64; max_iter = 512, eps = 1e-8, tiny = 1e-10)
k, xs = power(a - d * diagm(0 => ones(size(a, 1))), max_iter=max_iter, eps=eps, tiny=tiny)
k + d, xs
end
どの関数も引数 a が行列で、d がシフトする値です。Julia の場合、a - dI は a - d * diagm(0 => ones(size(a, 1))) で求めることができます。あとは、適切な関数 (power_inv, power) を呼び出して、求めた固有値に d を加算するだけです。
それでは実際に試してみましょう。最初は逆反復法です。ご参考までに収束の様子を表示しています。
julia> a = [1. 4. 5.; 4. 2. 6.; 5. 6. 3.]
3×3 Matrix{Float64}:
1.0 4.0 5.0
4.0 2.0 6.0
5.0 6.0 3.0
julia> eigen(a)
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> power_inv_shift(a, 12.1)
1 3.195039716555222 [0.488765, 0.580469, 0.651279]
2 13.16183309095262 [0.49664, 0.577333, 0.648101]
3 13.16290613028577 [0.4966, 0.57735, 0.648117]
4 13.162906158845166 [0.4966, 0.57735, 0.648117]
(12.175971065046905, [0.4966, 0.57735, 0.648117])
julia> power_inv_shift(a, -3.6)
1 -0.8108108108108013 [-0.175618, -0.658568, 0.731742]
2 -14.160434261633576 [0.321257, 0.571438, -0.755151]
3 -14.55800683665054 [-0.312466, -0.577721, 0.754058]
4 -14.559616430436828 [0.313018, 0.577327, -0.754131]
5 -14.55962279027976 [-0.312984, -0.577352, 0.754126]
6 -14.559622815406422 [0.312986, 0.57735, -0.754126]
(-3.668683097953267, [-0.312986, -0.57735, 0.754126])
julia> power_inv_shift(a, -2.5)
1 -89.99999999999244 [-0.810186, 0.576133, 0.108025]
2 -137.21166936789768 [0.80959, -0.577342, -0.106021]
3 -137.21247460733682 [-0.809585, 0.57735, 0.10601]
4 -137.21247463815925 [0.809585, -0.57735, -0.10601]
(-2.5072879670936414, [-0.809585, 0.57735, 0.10601])
julia> b = 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> b += 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> eigen(b)
Eigen{Float64,Float64,Matrix{Float64},Vector{Float64}}
eigenvalues:
4-element Vector{Float64}:
4.296089645312119
5.392275290272983
6.5077487053636425
9.80388635905124
eigenvectors:
4×4 Matrix{Float64}:
0.905684 -0.225903 0.135941 -0.332002
-0.380963 -0.801782 0.226102 -0.401113
-0.157381 0.517536 0.671404 -0.506561
-0.0991762 0.19563 -0.692542 -0.687225
julia> power_inv_shift(b, 5.3)
1 -0.22425084342131507 [-0.0858362, 0.903939, -0.387402, -0.159519]
2 9.39906651478699 [0.257113, 0.789462, -0.519894, -0.200881]
3 10.823449646815815 [0.223123, 0.80303, -0.516856, -0.195497]
4 10.837021987942329 [0.226165, 0.801677, -0.517567, -0.195672]
5 10.837136434204796 [0.225879, 0.801792, -0.51753, -0.195628]
6 10.8371373992314 [0.225905, 0.801781, -0.517536, -0.19563]
(5.392275290273574, [0.225903, 0.801782, -0.517535, -0.19563])
julia> power_inv_shift(b, 6.6)
1 -0.5642023346303511 [-0.366356, -0.189494, -0.505318, 0.757977]
2 -9.963758148788951 [0.148462, 0.230594, 0.66566, -0.694033]
3 -10.837782504602306 [-0.136476, -0.226509, -0.670924, 0.692769]
4 -10.839948840139247 [0.135969, 0.226141, 0.671373, -0.692554]
5 -10.839956238031057 [-0.135942, -0.226105, -0.671402, 0.692543]
6 -10.839956272963281 [0.135941, 0.226102, 0.671404, -0.692542]
(6.507748705363639, [-0.135941, -0.226102, -0.671404, 0.692542])
逆反復法の場合、固有値の近似が良ければ収束も速くなるようで、高い精度の固有値とその固有ベクトルを高速に求めることができます。
次はシフト法です。
julia> a
3×3 Matrix{Float64}:
1.0 4.0 5.0
4.0 2.0 6.0
5.0 6.0 3.0
julia> power(a)
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])
julia> power_shift(a, 3.)
1, 0.0, [1.0, 0.0, 0.0]
2, -2.0, [-0.298142, 0.596285, 0.745356]
3, 1.1555555555555561, [0.892008, 0.356803, 0.277514]
4, 4.491159135559923, [0.124621, 0.589555, 0.798056]
・・・省略・・・
30, 9.175970956226864, [0.49658, 0.577316, 0.648162]
31, 9.175971007573628, [0.496614, 0.577375, 0.648084]
32, 9.175971034691985, [0.49659, 0.577332, 0.648141]
33, 9.175971049014588, [0.496607, 0.577364, 0.648099]
(12.175971056579188, [0.496594, 0.577341, 0.648129])
julia> power_shift(a, -3.)
1, 0.0, [0.0, 0.0, 1.0]
2, 6.0, [0.507673, 0.609208, 0.609208]
3, 15.134020618556699, [0.495752, 0.576144, 0.649837]
4, 15.175889771554242, [0.496626, 0.577411, 0.648042]
5, 15.17597090760561, [0.496598, 0.577348, 0.64812]
6, 15.175971064741647, [0.4966, 0.57735, 0.648117]
(12.17597106504631, [0.4966, 0.57735, 0.648117])
julia> b
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(b)
1, 0.0, [0.0, 0.0, 0.0, 1.0]
2, 8.0, [0.122169, 0.122169, 0.122169, 0.977356]
3, 8.716417910447761, [0.206842, 0.220632, 0.234421, 0.923896]
4, 9.234074919186156, [0.258927, 0.288519, 0.32107, 0.864076]
・・・省略・・・
23, 9.803886246828288, [0.331985, 0.401085, 0.506479, 0.68731]
24, 9.803886309604092, [0.331991, 0.401094, 0.506507, 0.687282]
25, 9.803886337263988, [0.331995, 0.401101, 0.506525, 0.687263]
26, 9.803886349451385, [0.331997, 0.401105, 0.506537, 0.68725]
(9.803886354821369, [0.331999, 0.401108, 0.506545, 0.687242])
julia> power_shift(b, 4.)
1, 0.0, [0.0, 0.0, 0.0, 1.0]
2, 4.0, [0.229416, 0.229416, 0.229416, 0.917663]
3, 5.263157894736842, [0.297133, 0.339581, 0.382029, 0.806505]
4, 5.690090090090091, [0.318991, 0.378339, 0.452523, 0.741841]
5, 5.782180397336416, [0.326812, 0.392174, 0.483169, 0.711296]
6, 5.799820963675336, [0.32985, 0.397455, 0.496432, 0.6977]
7, 5.803127118525351, [0.331091, 0.399578, 0.502176, 0.691762]
8, 5.803744616019188, [0.331613, 0.40046, 0.504664, 0.689187]
9, 5.803859897236706, [0.331835, 0.400833, 0.505741, 0.688073]
10, 5.803881418854967, [0.33193, 0.400993, 0.506206, 0.687592]
11, 5.803885436750952, [0.331971, 0.401061, 0.506408, 0.687384]
12, 5.803886186863651, [0.331989, 0.401091, 0.506495, 0.687294]
13, 5.803886326904883, [0.331996, 0.401103, 0.506533, 0.687255]
14, 5.80388635304972, [0.331999, 0.401109, 0.506549, 0.687238]
(9.8038863579308, [0.332001, 0.401111, 0.506556, 0.687231])
julia> power_shift(b, -4.)
1, 0.0, [0.0, 0.0, 0.0, 1.0]
2, 12.0, [0.0824786, 0.0824786, 0.0824786, 0.989743]
3, 12.489795918367347, [0.150552, 0.157098, 0.163644, 0.962224]
4, 12.910835939843178, [0.202964, 0.219081, 0.236204, 0.924671]
・・・省略・・・
34, 13.803886271700657, [0.331985, 0.401084, 0.506479, 0.687311]
35, 13.8038863084367, [0.331989, 0.401091, 0.506498, 0.687291]
36, 13.803886329722985, [0.331992, 0.401097, 0.506513, 0.687275]
37, 13.80388634205714, [0.331994, 0.401101, 0.506525, 0.687263]
(9.803886349204083, [0.331996, 0.401104, 0.506533, 0.687254])
原点移動の場合、シフトする値によって収束回数が増えたり減ったりします。何かしらの選択基準があればよいのですが、よくわかりませんでした。シフト値の最適値は行列によって変わるでしょうから、明確な選択基準がなければ (M.Hiroi が知らないだけかもしれませんが)、ちょっと使いにくい方法のような気がします。興味のある方はいろいろ試してみてください。
実対称行列の固有値と固有ベクトルを求める簡単な方法に「ヤコビ法 (Jacobi method)」があります。なお、連立一次方程式を反復法で解く方法にもヤコビ法があります。拙作のページ ヤコビ法 (反復法) で説明しているので、興味のある方はお読みくださいませ。
ヤコビ法のポイントは、直交行列 \(X_1\) を使って相似変換 \(A_1 = {X_1}^{\mathrm{T}} A X_1\) を行うところです。このとき、\(A\) の非対角成分のひとつが 0 になるように直交行列 \(X_1\) を定めます。つまり、相似変換 \(A_n = {X_n}^{\mathrm{T}} A_{n-1} X_n\) を繰り返して \(A\) を対角化するわけです。対角化が完了したとき (非対角成分がすべて 0 になったとき)、\(A_n\) の対角成分が固有値になり、\(X_1 X_2 \cdots 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 * (cos2(r) - sin2(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 = atan(2 * z / (y - x)) / 2, (x != y)
r = π/ 4, (x == y)
実際の数値計算では、角度ではなくて sir(r) と cos(r) の値を求めたほうが誤差が少なくなるようです。sin(r) と cos(r) の求め方は 参考 URL 4 がわかりやすくまとまっていると思います。優れたドキュメントとプログラムを公開されている fussy さんに感謝いたします。今回はプログラムを簡単にするため、atan() で角度を求めることにしましょう。
それでは、これで本当に 0 になるのか実際に試してみましょう。
julia> givens(r) = [cos(r) sin(r); -sin(r) cos(r)]
givens (generic function with 1 method)
julia> a = [2. 1.; 1. 3.]
2×2 Matrix{Float64}:
2.0 1.0
1.0 3.0
julia> x = givens(0.5 * atan(2 / (3 - 2)))
2×2 Matrix{Float64}:
0.850651 0.525731
-0.525731 0.850651
julia> x' * a * x
2×2 Matrix{Float64}:
1.38197 -1.11022e-16
2.22045e-16 3.61803
julia> eigen(a)
Eigen{Float64,Float64,Matrix{Float64},Vector{Float64}}
eigenvalues:
2-element Vector{Float64}:
1.381966011250105
3.618033988749895
eigenvectors:
2×2 Matrix{Float64}:
-0.850651 0.525731
0.525731 0.850651
julia> b = [2. 1.; 1. 2.]
2×2 Matrix{Float64}:
2.0 1.0
1.0 2.0
julia> pi
π = 3.1415926535897...
julia> y = givens(pi / 4.)
2×2 Matrix{Float64}:
0.707107 0.707107
-0.707107 0.707107
julia> y' * b * y
2×2 Matrix{Float64}:
1.0 2.22045e-16
0.0 3.0
julia> eigen(b)
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> c = rand(2, 2)
2×2 Matrix{Float64}:
0.897777 0.648882
0.136918 0.703306
julia> c[1, 2] = c[2, 1]
0.13691784459147915
julia> c
2×2 Matrix{Float64}:
0.897777 0.136918
0.136918 0.703306
julia> z = givens(0.5 * atan(2 * c[1, 2] / (c[2, 2] - c[1, 1])))
2×2 Matrix{Float64}:
0.888543 -0.458794
0.458794 0.888543
julia> z' * c * z
2×2 Matrix{Float64}:
0.968474 5.55112e-17
5.55112e-17 0.632609
julia> eigen(c)
Eigen{Float64,Float64,Matrix{Float64},Vector{Float64}}
eigenvalues:
2-element Vector{Float64}:
0.6326094347289221
0.968473829063695
eigenvectors:
2×2 Matrix{Float64}:
0.458794 -0.888543
-0.888543 -0.458794
確かにギブンス回転を使って非対角成分を 0 にすることができました。
それではプログラムを作りましょう。ヤコビ法は Julia を使うと簡単にプログラムすることができます。次のリストを見てください。
リスト : ヤコビ法
# ヤコビ法
function jacobi(a::Matrix{Float64}, max_iter = 512, e = 1e-14)
n = size(a, 1)
vs = ones(n)
ks = diagm(0 => vs) # 固有ベクトル
xs = Matrix{Float64}(undef, n, n)
for i = 1 : max_iter
# 絶対値最大の非対角成分を探す
m, j = findmax(abs.(triu(a, 1)))
x, y = Tuple(j)
# println("$i, $x, $y, $(a[x, y]), $(a[y, x])")
if m < e
return diag(a), ks
end
# 角度を求める
d = a[y, y] - a[x, x]
if d == 0.0
r = pi / 4.0
else
r = 0.5 * atan(2.0 * a[x, y] / d)
end
# 直交行列の生成
fill!(xs, 0.0)
xs[diagind(xs)] = vs
xs[x, x] = cos(r)
xs[y, y] = cos(r)
xs[x, y] = sin(r)
xs[y, x] = - sin(r)
#
a = xs' * a * xs
ks = ks * xs
end
end
引数 a は実対称行列で、max_iter が繰り返しの最大値、e が 0 と判断する値です。変数 ks に固有ベクトルを格納します。ks は単位行列に初期化します。非対角成分の選択方法ですが、今回は絶対値最大の成分から順に消していくことにします。a は対称行列なので、対角成分を除いた上三角行列の中から最大値を findmax() で求めています。返り値は見つけた最大値とその位置です。
引数が行列の場合、findmax(), findmin() が返す添字のデータ型は CartesianIndex(x, y) になるので、Tuple() でタプルに変換して変数 x, y にセットしています。簡単な例を示しましょう。
julia> a = rand(5, 5)
5×5 Matrix{Float64}:
0.513379 0.525335 0.847615 0.0579795 0.558406
0.504764 0.74323 0.0838081 0.490538 0.235075
0.922555 0.323586 0.201859 0.597376 0.316467
0.139468 0.964848 0.261251 0.941041 0.0879425
0.944212 0.819577 0.485142 0.679319 0.116133
julia> m, i = findmax(a)
(0.9648481555934574, CartesianIndex(4, 2))
julia> x, y = Tuple(i)
(4, 2)
julia> a[i]
0.9648481555934574
julia> a[x, y]
0.9648481555934574
julia> m, i = findmin(a)
(0.057979472597648574, CartesianIndex(1, 4))
julia> x, y = Tuple(i)
(1, 4)
julia> a[i]
0.057979472597648574
julia> a[x, y]
0.057979472597648574
モジュール LinearAlgebra の関数 triu() は引数の行列を上三角行列に、tril() は下三角行列に変換します。二番目の引数に整数値を指定すると、取り出す三角行列の範囲を変更することができます。簡単な例を示しましょう。
julia> a = 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> triu(a)
5×5 Matrix{Float64}:
1.0 1.0 1.0 1.0 1.0
0.0 1.0 1.0 1.0 1.0
0.0 0.0 1.0 1.0 1.0
0.0 0.0 0.0 1.0 1.0
0.0 0.0 0.0 0.0 1.0
julia> tril(a)
5×5 Matrix{Float64}:
1.0 0.0 0.0 0.0 0.0
1.0 1.0 0.0 0.0 0.0
1.0 1.0 1.0 0.0 0.0
1.0 1.0 1.0 1.0 0.0
1.0 1.0 1.0 1.0 1.0
julia> triu(a, 1)
5×5 Matrix{Float64}:
0.0 1.0 1.0 1.0 1.0
0.0 0.0 1.0 1.0 1.0
0.0 0.0 0.0 1.0 1.0
0.0 0.0 0.0 0.0 1.0
0.0 0.0 0.0 0.0 0.0
julia> triu(a, -1)
5×5 Matrix{Float64}:
1.0 1.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0 1.0
0.0 1.0 1.0 1.0 1.0
0.0 0.0 1.0 1.0 1.0
0.0 0.0 0.0 1.0 1.0
最大値 a[x, y] の絶対値が e よりも小さくなったならば、非対角成分がすべて 0 になったと判断して、固有値 (a の対角成分) と固有ベクトル ks を返します。そうでなければ、直交行列 xs を生成して a を相似変換して、ks を ks * xs の値に更新します。
関数 diagind(A) は行列 A の対角成分の座標を返します。簡単な例を示しましょう。
julia> a = rand(5, 5)
5×5 Matrix{Float64}:
0.239522 0.663719 0.396802 0.924408 0.00656942
0.315202 0.387456 0.908763 0.280866 0.87524
0.980997 0.755987 0.492619 0.469537 0.191903
0.377995 0.506459 0.0658157 0.705611 0.123246
0.926698 0.77 0.784128 0.481811 0.418447
julia> diagind(a)
1:6:25
julia> a[diagind(a)]
5-element Vector{Float64}:
0.23952150080966295
0.3874563894608829
0.49261909906925605
0.7056106623795324
0.4184469857223021
julia> a[diagind(a)] = ones(5)
5-element Vector{Float64}:
1.0
1.0
1.0
1.0
1.0
julia> a
5×5 Matrix{Float64}:
1.0 0.663719 0.396802 0.924408 0.00656942
0.315202 1.0 0.908763 0.280866 0.87524
0.980997 0.755987 1.0 0.469537 0.191903
0.377995 0.506459 0.0658157 1.0 0.123246
0.926698 0.77 0.784128 0.481811 1.0
それでは実際に試してみましょう。収束の様子を見るため、0 にする対角成分の添字と値を表示しています。
julia> jacobi([2. 1.; 1. 3.])
1, 1, 2, 1.0, 1.0
2, 1, 2, -1.1102230246251565e-16, 2.220446049250313e-16
([1.38197, 3.61803], [0.850651 0.525731; -0.525731 0.850651])
julia> jacobi([2. 1.; 1. 2.])
1, 1, 2, 1.0, 1.0
2, 1, 2, 2.220446049250313e-16, 0.0
([1.0, 3.0], [0.707107 0.707107; -0.707107 0.707107])
julia> jacobi([1. 4. 5.; 4. 2. 6.; 5. 6. 3.])
1, 2, 3, 6.0, 6.0
2, 1, 3, 6.387849389602061, 6.387849389602061
3, 1, 2, -0.3837300281922351, -0.38373002819223523
4, 2, 3, -0.20518562849673092, -0.20518562849673116
5, 1, 3, 0.07764848736577329, 0.07764848736577315
6, 1, 2, 0.0010057157079803223, 0.0010057157079802084
7, 2, 3, 5.318604906830139e-6, 5.318604906608818e-6
8, 1, 3, 4.605674995833599e-9, 4.605674853394528e-9
9, 1, 2, -1.5459956007250589e-15, -1.6600536215271706e-15
([-2.50729, -3.66868, 12.176], [0.809585 0.312986 0.4966; -0.57735 0.57735 0.57735; -0.10601 -0.754126 0.648117])
julia> d = 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> d += 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> jacobi(d)
1, 1, 2, 1.0, 1.0
2, 2, 3, 1.3763819204711736, 1.3763819204711736
3, 3, 4, 1.6580314962279525, 1.6580314962279525
4, 1, 3, 0.4011483560399159, 0.4011483560399158
5, 2, 4, 0.27753210704204684, 0.2775321070420468
6, 1, 2, -0.2418160453957441, -0.24181604539574428
7, 2, 3, 0.2300801682095411, 0.23008016820954125
8, 3, 4, 0.06517159208983365, 0.06517159208983363
9, 1, 3, 0.05297008305598126, 0.05297008305598116
10, 1, 4, 0.0076388221073113265, 0.00763882210731137
11, 2, 4, -0.005417016203362336, -0.005417016203362362
12, 1, 2, -0.002742939677036741, -0.002742939677037079
13, 2, 3, -0.00013318242271641932, -0.00013318242271624547
14, 3, 4, 7.411922640029331e-5, 7.411922640029834e-5
15, 1, 4, 1.332060065787488e-5, 1.3320600657920408e-5
16, 1, 3, -5.867295797949511e-7, -5.867295799118133e-7
17, 2, 4, -3.1094169513768315e-8, -3.109416954510149e-8
18, 1, 2, -1.7534582710833612e-11, -1.7534921680978462e-11
19, 3, 4, -3.533813375839731e-12, -3.5338083403353147e-12
20, 2, 3, -6.992031425874347e-13, -6.990293571343703e-13
21, 1, 4, 1.0692825552099565e-19, 4.5636996101360793e-17
([4.29609, 5.39228, 9.80389, 6.50775],
[0.905684 0.225903 0.332002 -0.135941;
-0.380963 0.801782 0.401113 -0.226102;
-0.157381 -0.517536 0.506561 -0.671404;
-0.0991762 -0.19563 0.687225 0.692542])
ヤコビ法の場合、相似変換のとき以前 0 にした非対角成分が 0 でなくなることがあります。それでも、最初よりは 0 に近い値になるので、相似変換を繰り返し行っていけば、非対角成分は 0 に近づいていきます。最終的にはすべての非対角成分を 0 にすることができますが、それなりの回数が必要になります。特に、行列 (次元数 N) が大きくなると、収束するまでに必要な回数は大幅に増えてしまいます。このため、ヤコビ法が実用になるのは N が数十程度までといわれています。
#
# eig.jl : 固有値と固有ベクトル
#
# Copyright (C) 2018-2021 Makoto Hiroi
#
using LinearAlgebra
#
# 累乗法
#
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
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
#
# LU 分解による連立一次方程式の解法
#
# ピボット選択
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(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 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
# 逆反復法
function power_inv_shift(a::Matrix{Float64}, d::Float64; max_iter = 512, eps = 1e-8, tiny = 1e-10)
k, xs = power_inv(a - d * diagm(0 => ones(size(a, 1))), max_iter=max_iter, eps=eps, tiny=tiny)
k + d, xs
end
# 原点移動
function power_shift(a::Matrix{Float64}, d::Float64; max_iter = 512, eps = 1e-8, tiny = 1e-10)
k, xs = power(a - d * diagm(0 => ones(size(a, 1))), max_iter=max_iter, eps=eps, tiny=tiny)
k + d, xs
end
#
# ヤコビ法
#
function jacobi(a::Matrix{Float64}, max_iter = 512, e = 1e-14)
n = size(a, 1)
vs = ones(n)
ks = diagm(0 => vs) # 固有ベクトル
xs = Matrix{Float64}(undef, n, n)
for i = 1 : max_iter
# 絶対値最大の非対角成分を探す
m, j = findmax(abs.(triu(a, 1)))
x, y = Tuple(j)
# println("$i, $x, $y, $(a[x, y]), $(a[y, x])")
if m < e
return diag(a), ks
end
# 角度を求める
d = a[y, y] - a[x, x]
if d == 0.0
r = pi / 4.0
else
r = 0.5 * atan(2.0 * a[x, y] / d)
end
# 直交行列の生成
fill!(xs, 0.0)
xs[diagind(xs)] = vs
xs[x, x] = cos(r)
xs[y, y] = cos(r)
xs[x, y] = sin(r)
xs[y, x] = - sin(r)
#
a = xs' * a * xs
ks = ks * xs
end
end