固有値 \(\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\) に累乗法を適用したほうが収束が速くなるはずです。これを「原点移動 (シフト法)」といいます。
ただし、参考文献によっては、逆反復法のことをシフト法と呼ぶことがあるようです。まあ、どちらの方法も対角成分の値をシフトするので、シフト法というと混乱するような気がします。本稿では前者を逆反復法、後者を原点移動と呼ぶことにします。
プログラムは簡単です。次のリストを見てください。
リスト : 逆反復法と原点移動 -- 逆反復法 powerInvShift :: Matrix R -> R -> Maybe (R, Matrix R, Int) powerInvShift xs d = powerInv(xs - (scalar d) * ident (rows xs)) >>= \(k, ys, i) -> return (k + d, ys, i) -- 原点移動法 powerShift :: Matrix R -> R -> Maybe (R, Matrix R, Int) powerShift xs d = power(xs - (scalar d) * ident (rows xs)) >>= \(k, ys, i) -> return (k + d, ys, i)
どの関数も引数 xs が行列で、d がシフトする値です。hmatrix の場合、xs - dI は xs - (scalar d) * ident (rows xs) で求めることができます。あとは、適切な関数 (powerInv, power) を呼び出して、求めた固有値に d を加算するだけです。
それでは実際に試してみましょう。最初は逆反復法です。
ghci> a = matrix 3 [1,4,5,4,2,6,5,6,3] ghci> a (3><3) [ 1.0, 4.0, 5.0 , 4.0, 2.0, 6.0 , 5.0, 6.0, 3.0 ] ghci> eigenvaluesSH $ sym a [12.175971065046909,-2.5072879670936405,-3.668683097953265] ghci> powerInv a Just (-2.5072879960642904,(3><1) [ -0.8096263206940646 , 0.5772748809066349 , 0.10610811803826374 ],21) ghci> powerInvShift a (-2.5) Just (-2.5072879670936414,(3><1) [ -0.8095854617408919 , 0.5773502691875207 , 0.10600965430980433 ],5) ghci> powerInvShift a (-3.6) Just (-3.6686830979532674,(3><1) [ -0.3129856690762894 , -0.5773502749784004 , 0.7541264024918053 ],7) ghci> powerInvShift a 12.1 Just (12.175971065046905,(3><1) [ 0.49659978454065673 , 0.577350269192263 , 0.6481167492495429 ],5) ghci> b = matrix 4 [5,1,1,1,1,6,1,1,1,1,7,1,1,1,1,8] ghci> b (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> eigenvaluesSH $ sym b [9.803886359051248,6.507748705363647,5.392275290272981,4.296089645312118] ghci> powerInv b Just (4.296089899617464,(4><1) [ 0.9057807814119976 , -0.38061757467977975 , -0.15760282339290624 , -9.926121059173881e-2 ],28) ghci> powerInvShift b 4.2 Just (4.296089645312182,(4><1) [ 0.9056835800330293 , -0.38096255422803676 , -0.15738127537430335 , -9.917620327559983e-2 ],6) ghci> powerInvShift b 5.3 Just (5.392275290273574,(4><1) [ 0.2259027659665049 , 0.8017817062691626 , -0.5175354688082241 , -0.19562994235608297 ],7) ghci> powerInvShift b 6.5 Just (6.507748705363649,(4><1) [ 0.13594052685630392 , 0.22610178938299505 , 0.6714043318290509 , -0.6925419678236284 ],5) ghci> powerInvShift b 9.8 Just (9.80388635905125,(4><1) [ 0.3320019640610339 , 0.40111308352620756 , 0.5065613134848598 , 0.6872253093158076 ],4)
逆反復法の場合、固有値の近似が良ければ収束も速くなるようで、高い精度の固有値とその固有ベクトルを高速に求めることができます。
次はシフト法です。
ghci> a (3><3) [ 1.0, 4.0, 5.0 , 4.0, 2.0, 6.0 , 5.0, 6.0, 3.0 ] ghci> power a Just (12.175971064806813,(3><1) [ 0.49659938075279125 , 0.5773496196304903 , 0.6481176372761912 ],11) ghci> powerShift a 3 Just (12.175971056579188,(3><1) [ 0.49659446247373806 , 0.5773406147155016 , 0.6481294272291191 ],33) ghci> powerShift a (-3) Just (12.175971065046681,(3><1) [ 0.49659978753546913 , 0.5773502707500971 , 0.648116745567121 ],6) ghci> b (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> power b Just (9.803886355452113,(4><1) [ 0.33200504131765435 , 0.4011183591424394 , 0.5065758276117462 , 0.6872100445846333 ],24) ghci> powerShift b 4 Just (9.80388635806045,(4><1) [ 0.3320029837555454 , 0.40111478205411727 , 0.5065663400662365 , 0.6872201201267257 ],13) ghci> powerShift b (-4) Just (9.803886350676947,(4><1) [ 0.33200756559204925 , 0.4011228418464535 , 0.5065863062581886 , 0.6871984840279777 ],34)
原点移動の場合、シフトする値によって収束回数が増えたり減ったりします。何かしらの選択基準があればよいのですが、よくわかりませんでした。シフト値の最適値は行列によって変わるでしょうから、明確な選択基準がなければ (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『固有値問題 (1) 対称行列の固有値』がわかりやすくまとまっていると思います。優れたドキュメントとプログラムを公開されている fussy さんに感謝いたします。今回はプログラムを簡単にするため、atan() で角度を求めることにしましょう。
それでは、これで本当に 0 になるのか実際に試してみましょう。
ghci> import Prelude hiding ((<>)) ghci> :m + Numeric.LinearAlgebra ghci> givens r = (2><2) [cos r, sin r, -(sin r), cos r] ghci> a = matrix 2 [2,1,1,3] ghci> a (2><2) [ 2.0, 1.0 , 1.0, 3.0 ] ghci> x = givens(0.5 * atan(2 / (3 - 2))) ghci> x (2><2) [ 0.85065080835204, 0.5257311121191336 , -0.5257311121191336, 0.85065080835204 ] ghci> (tr x) <> a <> x (2><2) [ 1.381966011250105, 2.220446049250313e-16 , -1.1102230246251565e-16, 3.6180339887498953 ] ghci> eigenvalues a [1.381966011250105 :+ 0.0,3.618033988749895 :+ 0.0] ghci> b = matrix 2 [2,1,1,2] ghci> b (2><2) [ 2.0, 1.0 , 1.0, 2.0 ] ghci> y = givens(pi / 4) ghci> y (2><2) [ 0.7071067811865476, 0.7071067811865475 , -0.7071067811865475, 0.7071067811865476 ] ghci> (tr y) <> b <> y (2><2) [ 1.0, 0.0 , 2.220446049250313e-16, 3.0 ] ghci> eigenvaluesSH $ sym b [3.0,1.0] ghci> c <- rand 2 2 ghci> c (2><2) [ 0.58682165089381, 0.9242200036180299 , 0.23022081015176177, 0.9368664384525579 ] ghci> c1 = sym c ghci> c1 Herm (2><2) [ 0.58682165089381, 0.5772204068848958 , 0.5772204068848958, 0.9368664384525579 ] ghci> c2 = unSym c1 ghci> z = givens(0.5 * atan(2 * (c2 ! 0 ! 1) / ((c2 ! 1 ! 1) - (c2 ! 0 ! 0)))) ghci> z (2><2) [ 0.8031718519624857, 0.5957474097427121 , -0.5957474097427121, 0.8031718519624857 ] ghci> (tr z) <> c2 <> z (2><2) [ 0.1586722312993933, -1.1102230246251565e-16 , -1.3877787807814457e-17, 1.365015858046975 ] ghci> eigenvaluesSH c1 [1.3650158580469747,0.15867223129939326]
確かにギブンス回転を使って非対角成分を 0 にすることができました。
それではプログラムを作りましょう。ヤコビ法は hmatrix を使うと簡単にプログラムすることができます。次のリストを見てください。
リスト : ギブンス回転
givens :: Matrix R -> Int -> Int -> Matrix R
givens xs x y =
accum (ident n) const [((x,x), cos r), ((y,y), cos r), ((x,y), sin r), ((y,x), - sin r)]
where
n = rows xs
d = (xs ! y ! y) - (xs ! x ! x)
r = if d == 0 then pi / 4 else 0.5 * atan(2.0 * (xs ! x ! y) / d)
関数 givens はギブンス回転用の行列を生成します。x, y が座標を表します。関数 accum は連想リストの値とベクトル (行列) の成分を関数に渡して評価し、その結果で置き換えたベクトル (行列) を返します。
accum :: Container c e => c e -> (e -> e -> e) -> [(IndexOf c, e)] -> c e
簡単な例を示しましょう。
ghci> b = matrix 3 [1..9] ghci> b (3><3) [ 1.0, 2.0, 3.0 , 4.0, 5.0, 6.0 , 7.0, 8.0, 9.0 ] ghci> accum b (+) [((0,0),10), ((2,1),20), ((1,2), 30)] (3><3) [ 11.0, 2.0, 3.0 , 4.0, 5.0, 36.0 , 7.0, 28.0, 9.0 ] ghci> accum b (*) [((0,0),10), ((2,1),20), ((1,2), 30)] (3><3) [ 10.0, 2.0, 3.0 , 4.0, 5.0, 180.0 , 7.0, 160.0, 9.0 ] ghci> accum b const [((0,0),10), ((2,1),20), ((1,2), 30)] (3><3) [ 10.0, 2.0, 3.0 , 4.0, 5.0, 30.0 , 7.0, 20.0, 9.0 ]
accum に const を渡すと、行列の成分は連想リストで指定した値に置き換わることに注意してください。あとは ident n で単位行列を生成し、(x,x), (y,y), (x,y), (y,x) の成分を cos r, cos r, sin r, -sin r に置き換えれば、ギブンス回転用の行列を生成することができます。
次はヤコビ法で固有値と固有ベクトルを求める関数 jacobi を作ります。
リスト : ヤコビ法
-- 非対角線の成分で絶対値最大の位置を求める
searchMaxIdx :: Matrix R -> (Int, Int)
searchMaxIdx xs = maxIndex $ cmap abs (xs - diag (takeDiag xs))
-- ヤコビ法
jacobi :: Matrix R -> Maybe (Vector R, Matrix R, Int)
jacobi xs =
iter 0 xs (ident (rows xs))
where
iter i xs ks
| i >= 512 = Nothing
| otherwise = let
(x, y) = searchMaxIdx xs
zs = givens xs x y
in
if abs (xs ! x ! y) < 1e-14
then
Just (takeDiag xs, ks, i + 1)
else
iter (i + 1) ((tr zs) <> xs <> zs) (ks <> zs)
引数 xs は実対称行列で、実際の処理は局所関数 iter で行います。引数 i が繰り返しの回数、xs が実対称行列、ks が求める固有ベクトルになります。ks の初期値は単位行列です。i が 512 以上になったら収束しなかったとして Nothing を返します。
非対角成分の選択方法ですが、今回は絶対値最大の成分から順に消していくことにします。関数 searchMaxIdx は行列 xs の対角成分以外で、絶対値最大の成分の位置を探します。xs から対角成分を取り出して、それを xs から引き算することで対角成分を 0 にします。そのあと、cmap で各成分に abs を適用して絶対値を求め、その中から関数 maxIndex で最大値の位置を求めます。
maxIndex :: Container c e => c e -> IndexOf c minIndex :: Container c e => c e -> IndexOf c maxElement :: Container c e => c e -> e minElement :: Container c e => c e -> e
IndexOf は添字を表すデータ型で、ベクトルであれば Int、行列であれば (Int, Int) になります。maxIndex, minIndex は最大値 (最小値) の添字を返します。maxElement, minElement は最大値 (最小値) を返します。簡単な例を示しましょう。
ghci> a = matrix 3 [1..9] ghci> a (3><3) [ 1.0, 2.0, 3.0 , 4.0, 5.0, 6.0 , 7.0, 8.0, 9.0 ] ghci> minIndex a (0,0) ghci> minElement a 1.0 ghci> maxIndex a (2,2) ghci> maxElement a 9.0
searchMaxIdx で絶対値最大の位置 (x, y) を求め、関数 givens でギブンス回転行列 (zs) を求めます。最大値が 1e-14 よりも小さくなったならば、非対角成分がすべて 0 になったと判断して、固有値 (xs の対角成分) と固有ベクトル ks と繰り返し回数 i を返します。そうでなければ、zs を使って xs を相似変換し、ks を ks <> zs の値に更新して処理を続けます。
それでは実際に試してみましょう。
ghci> jacobi $ matrix 2 [2,1,1,3] Just ([1.381966011250105,3.6180339887498953],(2><2) [ 0.85065080835204, 0.5257311121191336 , -0.5257311121191336, 0.85065080835204 ],2) ghci> jacobi $ matrix 2 [2,1,1,2] Just ([1.0,3.0],(2><2) [ 0.7071067811865476, 0.7071067811865475 , -0.7071067811865475, 0.7071067811865476 ],2) ghci> jacobi $ matrix 3 [1,4,5,4,2,6,5,6,3] Just ([-2.5072879670936397,-3.6686830979532643,12.175971065046907],(3><3) [ 0.8095854617397512, 0.3129856771935582, 0.49659978454619125 , -0.5773502691896246, 0.5773502691896267, 0.5773502691896257 , -0.106009654307056, -0.7541264035547062, 0.6481167492476514 ],9) ghci> jacobi $ matrix 4 [5,1,1,1,1,6,1,1,1,1,7,1,1,1,1,8] Just ([4.29608964531,5.39227529027,9.8038863590,6.50774870536],(4><4) [ 0.905683564482, 0.2259029659858, 0.332001964060, -0.1359405268657 , -0.3809626092113, 0.801781619542, 0.4011130835259, -0.226101789403 , -0.157381240530, -0.517535509923, 0.506561313484, -0.671404331813 , -9.91761893687e-2, -0.1956299580634, 0.687225309316, 0.692541967829 ],21)
ヤコビ法の場合、相似変換のとき以前 0 にした非対角成分が 0 でなくなることがあります。それでも、最初よりは 0 に近い値になるので、相似変換を繰り返し行っていけば、非対角成分は 0 に近づいていきます。最終的にはすべての非対角成分を 0 にすることができますが、それなりの回数が必要になります。特に、行列 (次元数 N) が大きくなると、収束するまでに必要な回数は大幅に増えてしまいます。このため、ヤコビ法が実用になるのは N が数十程度までといわれています。
--
-- eig.hs : 固有値と固有ベクトル
--
-- Copyright (C) 2021 Makoto Hiroi
--
import Prelude hiding ((<>))
import Numeric.LinearAlgebra
--
-- ピボット選択
--
findMaxRow :: [Vector R] -> Int -> Int
findMaxRow [] _ = 0
findMaxRow (z:zs) n =
iter zs 1 (abs (z ! n)) 0
where
iter [] _ _ b = b
iter (z:zs) i a b =
let c = abs (z ! n)
in if a < c
then iter zs (i + 1) c i
else iter zs (i + 1) a b
swapRow :: [a] -> Int -> Int -> [a]
swapRow zs m n =
iter zs 0
where
x = zs !! m
y = zs !! n
iter [] _ = []
iter (z:zs) i
| i == m = y : iter zs (i + 1)
| i == n = x : iter zs (i + 1)
| otherwise = z : iter zs (i + 1)
selectPivot' :: [Vector R] -> Int -> [Int] -> ([Vector R], [Int])
selectPivot' zs n idx =
let i = findMaxRow zs n
in if i == 0
then (zs, idx)
else (swapRow zs i 0, swapRow idx (n + i) n)
transfer :: [a] -> [Int] -> [a]
transfer vs xs = map (\k -> vs !! k) xs
-- 上三角行列 (ピボット選択)
makeUp' :: [Vector R] -> Int -> [Int] -> ([Vector R], [Int])
makeUp' [v] _ idx = ([v], idx)
makeUp' (v:vs) n idx =
let
a = map (\w -> let temp = (w ! n) / (v ! n)
in w - (v * scalar temp))
vs
(b, d) = selectPivot' a (n + 1) idx
(xs, ys) = makeUp' b (n + 1) d
in
(v : xs, ys)
-- 下三角行列
makeLow :: [Vector R] -> Int -> [[R]]
makeLow [_] _ = [[1]]
makeLow (v:vs) n =
let
a = map (\w -> let temp = (w ! n) / (v ! n)
in (w - (v * scalar temp), temp))
vs
b = makeLow (map fst a) (n + 1)
in
(1 : map snd a) : map (0:) b
-- LU 分解 (ピボット選択)
makeLU' :: Matrix R -> (Matrix R, Matrix R, [Int])
makeLU' xs =
let
a = toRows xs
(b, d) = selectPivot' a 0 [0..((rows xs) - 1)]
(ys, d1) = makeUp' b 0 d
zs = makeLow (transfer a d1) 0
in
(fromColumns $ map fromList zs, fromRows ys, d1)
-- 前進代入
lu1 :: [Vector R] -> Int -> [R] -> [R]
lu1 [] _ zs = zs
lu1 (v:vs) 0 [] =
lu1 vs 1 [(v ! (size v - 1)) / (v ! 0)]
lu1 (v:vs) n xs =
let
x = ((v ! (size v - 1)) - (subVector 0 n v <.> vector xs)) / (v ! n)
in
lu1 vs (n + 1) (xs ++ [x])
-- 後退代入
lu2 :: [Vector R] -> Int -> [R]
lu2 [v] n = [(v ! (n + 1)) / (v ! n)]
lu2 (v:vs) n =
let m = n + 1
x = lu2 vs m
y = subVector m (size v - 1 - m) v <.> vector x
z = ((v ! (size v - 1)) - y) / (v ! n)
in z : x
-- 連立方程式の解法 (ピボット選択)
luSolver' :: (Matrix R, Matrix R, [Int]) -> Vector R -> [R]
luSolver' (ls, us, idx) ys =
let
a = ls ||| (col $ transfer (toList ys) idx)
b = lu1 (toRows a) 0 []
c = us ||| col b
in
lu2 (toRows c) 0
--
-- 累乗法
--
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
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)
--
-- 逆累乗法
--
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
-- 逆反復法
powerInvShift :: Matrix R -> R -> Maybe (R, Matrix R, Int)
powerInvShift xs d =
powerInv(xs - (scalar d) * ident (rows xs)) >>= \(k, ys, i) -> return (k + d, ys, i)
-- 原点移動法
powerShift :: Matrix R -> R -> Maybe (R, Matrix R, Int)
powerShift xs d =
power(xs - (scalar d) * ident (rows xs)) >>= \(k, ys, i) -> return (k + d, ys, i)
-- ギブンス回転
givens :: Matrix R -> Int -> Int -> Matrix R
givens xs x y =
accum (ident n) const [((x,x), cos r), ((y,y), cos r), ((x,y), sin r), ((y,x), - sin r)]
where
n = rows xs
d = (xs ! y ! y) - (xs ! x ! x)
r = if d == 0 then pi / 4 else 0.5 * atan(2.0 * (xs ! x ! y) / d)
-- 非対角線の成分で絶対値最大の位置を求める
searchMaxIdx :: Matrix R -> (Int, Int)
searchMaxIdx xs = maxIndex $ cmap abs (xs - diag (takeDiag xs))
--
-- ヤコビ法
--
jacobi :: Matrix R -> Maybe (Vector R, Matrix R, Int)
jacobi xs =
iter 0 xs (ident (rows xs))
where
iter i xs ks
| i >= 512 = Nothing
| otherwise = let
(x, y) = searchMaxIdx xs
zs = givens xs x y
in
if abs (xs ! x ! y) < 1e-14
then
Just (takeDiag xs, ks, i + 1)
else
iter (i + 1) ((tr zs) <> xs <> zs) (ks <> zs)