前回は固有値と固有ベクトルを求める簡単な方法として「累乗法」と「ヤコビ法」を紹介しました。今回は行列の固有値を求める方法の一つである「QR 法 (QR algorithm)」と、その基礎になる「QR 分解 (QR decomposition)」について取り上げます。
QR 法の基本的な考え方はシンプルですが、これをそのままプログラムするのは効率的ではありません。このため、対称行列であれば「三重対角行列」に、それ以外の行列は「ヘッセンベルグ行列」に相似変換してから QR 法を適用するのが主流のようです。本稿では「実対称行列」の固有値を求める QR 法の基本について簡単に説明します。
まずは最初に QR 分解を説明しましょう。QR 分解は行列 \(A\) を直交行列 \(Q\) と上三角行列 \(R\) に分解することです。
基本的な考え方は簡単です。\(A\) の左側から適当な直交行列 \(Q_i\) をいくつか掛け算したら上三角行列 \(R\) になったとしましょう。これを式で示すと次のようになります。
直交行列の逆行列は転置行列なので、上の式は次のように変換できます。
直交行列の掛け算は直交行列になるので、行列 \(A\) は直交行列 \(Q\) と上三角行列 \(R\) に分解することができました。
QR 分解は「ギブンス回転」、「ハウスホルダー変換 (Householder transformation)」、「グラム・シュミット分解」などの手法を使って行うことができます。ここでは実装が簡単なギブンス回転を使って説明します。ギブンス回転の場合、直交行列 \(Q_i\) を掛け算するとき、対角成分よりも下の成分の一つ \(a_{ij} \ (i \lt 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 = atan(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 になります。
それでは実際に試してみましょう。
ghci> import Prelude hiding ((<>)) ghci> :m + Numeric.LinearAlgebra ghci> givens a c = let d = sqrt(a * a + c * c) in (2><2) [a/d, c/d, -c/d, a/d] ghci> a = matrix 2 [2,1,1,3] ghci> a (2><2) [ 2.0, 1.0 , 1.0, 3.0 ] ghci> q = givens 2 1 ghci> q (2><2) [ 0.8944271909999159, 0.4472135954999579 , -0.4472135954999579, 0.8944271909999159 ] ghci> q <> a (2><2) [ 2.23606797749979, 2.23606797749979 , 0.0, 2.23606797749979 ]
確かにギブンス回転を使って正方行列を上三角行列に変換することができました。
それではギブンス回転を使って正方行列を QR 分解するプログラムを作りましょう。次のリストを見てください。
リスト : QR 分解 (ギブンス回転, 効率は悪い)
-- ギブンス回転
mkGivens :: Matrix R -> Int -> Int -> Matrix R
mkGivens xs x y =
-- x 列, y 行
accum (ident n) const [((x,x), a/d), ((x,y), c/d), ((y,x), -c/d), ((y,y), a/d)]
where
n = rows xs
a = xs ! x ! x
c = xs ! y ! x
d = sqrt(a * a + c * c)
mkIndex :: Int -> [(Int, Int)]
mkIndex n = [(i, j) | i <- [0 .. n-1], j <- [i+1 .. n-1]]
qrGivens :: Matrix R -> (Matrix R, Matrix R)
qrGivens xs =
foldl (\(qs, ys) (i, j) ->
let zs = mkGivens ys i j in (qs <> tr zs, zs <> ys))
(ident n, xs)
(mkIndex n)
where
n = rows xs
関数 mkGivens は x 列 y 行目の成分を 0 にするギブンス回転行列を生成します。基本的な考え方は前回作成した関数 givens と同じです。関数 mkIndex は 0 にする成分の位置 (添字) を生成します。これはリスト内包表記を使えば簡単です。関数 qrGivens はギブンス回転行列を使って行列 xs を QR 分解します。mkIndex で添字を生成し、畳み込みを行う関数 foldl に渡して、位置 (i, j) の成分をギブンス回転で 0 にしていきます。
それでは簡単な実行例を示します。
ghci> (x, y) = qrGivens $ matrix 2 [2,1,1,3] ghci> disp 3 x 2x2 0.894 -0.447 0.447 0.894 ghci> disp 3 y 2x2 2.236 2.236 0.000 2.236 ghci> disp 3 $ x <> y 2x2 2 1 1 3 ghci> (x, y) = qrGivens $ matrix 2 [2,1,1,2] ghci> disp 3 x 2x2 0.894 -0.447 0.447 0.894 ghci> disp 3 y 2x2 2.236 1.789 0.000 1.342 ghci> disp 3 $ x <> y 2x2 2.000 1.000 1.000 2.000 ghci> (x, y) = qrGivens $ matrix 3 [1,4,5,4,2,6,5,6,3] ghci> disp 3 x 3x3 0.154 0.802 0.577 0.617 -0.535 0.577 0.772 0.267 -0.577 ghci> disp 3 y 3x3 6.481 6.481 6.789 0.000 3.742 1.604 -0.000 0.000 4.619 ghci> disp 3 $ x <> y 3x3 1.000 4.000 5.000 4.000 2.000 6.000 5.000 6.000 3.000 ghci> (x, y) = qrGivens $ matrix 4 [6,1,1,1,1,7,1,1,1,1,8,1,1,1,1,9] ghci> disp 3 x 4x4 0.961 -0.192 -0.152 -0.130 0.160 0.973 -0.126 -0.108 0.160 0.091 0.979 -0.093 0.160 0.091 0.059 0.981 ghci> disp 3 y 4x4 6.245 2.402 2.562 2.722 0.000 6.799 1.595 1.686 0.000 -0.000 7.609 1.227 -0.000 0.000 0.000 8.500 ghci> disp 3 $ x <> y 4x4 6.000 1.000 1.000 1.000 1.000 7.000 1.000 1.000 1.000 1.000 8.000 1.000 1.000 1.000 1.000 9.000
正常に動作していますね。なお、0 の精度が足りない (もっと 0 に近づけたい) 場合は、上三角行列に qrGivens を繰り返し適用すれば OK です。
ghci> (x, y) = qrGivens $ matrix 4 [6,1,1,1,1,7,1,1,1,1,8,1,1,1,1,9] ghci> x (4><4) [ 0.9607689228305228, -0.1923268901207253, -0.15176735883886408, -0.13000043480943765 ,0.16012815380508713, 0.972947797081316, -0.1264727990323867, -0.10833369567453138 , 0.1601281538050871,9.050677182151778e-2, 0.9785422790840111,-9.285745343531263e-2 , 0.1601281538050871,9.050677182151777e-2,5.8534672981559975e-2, 0.9811937579664701] ghci> y (4><4) [ 6.2449979983984, 2.401922307076307, 2.562050460881394, 2.722178614686481 ,1.0451130034545904e-17, 6.799321233091524,1.5951818533542508,1.6856886251757683 ,1.1027165024386951e-16,-1.1082526910889261e-16, 7.608632747782399,1.2271141780468002 ,-7.542458324249159e-18, 6.611457556382143e-18, 0.0, 8.499552237778948] ghci> x <> y (4><4) [ 6.000000000000002, 1.0,1.0000000000000002, 1.0 ,1.0000000000000002, 7.0,0.9999999999999998,0.9999999999999997 ,1.0000000000000002,0.9999999999999999, 8.0, 1.0 ,1.0000000000000002, 1.0, 1.0, 9.0] ghci> (x1, y1) = qrGivens y ghci> x1 (4><4) [ 1.0, -1.67352015760e-18, -1.765759577059e-17, 1.207759926613e-18 , 1.673520157608e-18, 1.0, 2.253716761223e-17, -1.399022452331e-18 , 1.765759577059e-17, -2.253716761223e-17, 1.0, -1.133773539631e-19 , -1.207759926613e-18, 1.399022452331e-18, 1.13377353963e-19, 1.0 ] ghci> y1 (4><4) [ 6.2449979983984, 2.401922307076307, 2.562050460881394, 2.722178614686481 , 2.777920381828957e-49, 6.799321233091524, 1.5951818533542506, 1.6856886251757683 , -1.232595164407831e-32, 0.0, 7.608632747782399, 1.2271141780468002 , 1.397483782483125e-51, 0.0, 0.0, 8.499552237778948 ] ghci> x <> x1 <> y1 (4><4) [ 6.000000000000002, 1.0, 1.0000000000000002, 1.0 , 1.0000000000000002, 7.0, 0.9999999999999998, 0.9999999999999997 , 1.0000000000000002, 0.9999999999999999, 8.0, 1.0 , 1.0000000000000002, 1.0, 1.0, 9.0 ]
実用的には、ギブンス回転よりも効率が良い「ハウスホルダー変換」を用いて QR 分解を行うのが主流のようです。ハウスホルダー変換は次回以降に取り上げたいと思います。なお、hmatrix には行列を QR 分解する関数 qr が用意されているので、私たちが作る必要はありません。
qr :: Field t => Matrix t -> (Matrix t, Matrix t)
簡単な使用例を示します。
ghci> c = matrix 3 [1,4,5,4,2,6,5,6,3] ghci> c (3><3) [ 1.0, 4.0, 5.0 , 4.0, 2.0, 6.0 , 5.0, 6.0, 3.0 ] ghci> (q, r) = qr c ghci> q (3><3) [ -0.15430334996209183, 0.8017837257372727, -0.577350269189626 , -0.6172133998483675, -0.5345224838248492, -0.5773502691896253 , -0.7715167498104593, 0.2672612419124245, 0.5773502691896258 ] ghci> r (3><3) [ -6.480740698407861, -6.480740698407859, -6.789347398332042 , 0.0, 3.7416573867739387, 1.603567451474541 , 0.0, 0.0, -4.618802153517005 ] ghci> disp 3 $ q <> r 3x3 1.000 4.000 5.000 4.000 2.000 6.000 5.000 6.000 3.000 ghci> d = matrix 4 [6,1,1,1,1,7,1,1,1,1,8,1,1,1,1,9] ghci> d (4><4) [ 6.0, 1.0, 1.0, 1.0 , 1.0, 7.0, 1.0, 1.0 , 1.0, 1.0, 8.0, 1.0 , 1.0, 1.0, 1.0, 9.0 ] ghci> (q, r) = qr d ghci> q (4><4) [ -0.960768922830, 0.192326890120, 0.1517673588388, -0.130000434809 , -0.1601281538050, -0.972947797081, 0.1264727990323, -0.108333695674 , -0.1601281538050, -9.05067718215e-2, -0.978542279084, -9.28574534353e-2 , -0.1601281538050, -9.05067718215e-2, -5.85346729815e-2, 0.981193757966 ] ghci> r (4><4) [ -6.244997998398398, -2.4019223070763065, -2.5620504608813937, -2.7221786146864804 , 0.0, -6.799321233091523, -1.5951818533542512, -1.685688625175769 , 0.0, 0.0, -7.608632747782399, -1.2271141780468002 , 0.0, 0.0, 0.0, 8.49955223777895 ] ghci> disp 3 $ q <> r 4x4 6.000 1.000 1.000 1.000 1.000 7.000 1.000 1.000 1.000 1.000 8.000 1.000 1.000 1.000 1.000 9.000
ところで、行列 \(A\) が正則行列であれば、QR 分解を使って連立一次方程式を解くことができます。次の式を見てください。
\(Q\) は直交行列なので、y の値は \(Q^{\mathrm{T}} b\) で求めることができます。\(R\) は上三角行列なので、ガウスの消去法と同様に後退代入で x を求めることができます。ガウスの消去法はピボット選択が必要になる場合がありますが、QR 分解だとその必要がないところが利点です。
プログラムと実行例を示します。
リスト : 連立一次方程式の解法
-- 後退代入 (LU 分解のプログラムと同じ)
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
-- 連立一次方程式の解法
qrSolver :: Matrix R -> Matrix R -> [R]
qrSolver xs ys =
let
(q, r) = qrGivens xs
zs = r ||| ((tr q) <> ys)
in
lu2 (toRows zs) 0
ghci> qrSolver (matrix 2 [1,1,2,4]) (col [100, 272]) [64.00000000000001,36.0] ghci> qrSolver (matrix 3 [1,1,1,2,4,6,2,0,4]) (col [10,38,14]) [3.0,5.0,2.0] ghci> qrSolver (matrix 3 [0,2,4,1,1,1,4,2,6]) (col [14,10,38]) [5.0,3.0,2.0] ghci> qrSolver (matrix 4 [1,1,1,1,-1,1,-1,1,8,4,2,1,-8,4,-2,1]) (col [-5,-7,-31,-35]) [-3.8949200278671966e-16,-9.0,1.0000000000000009,3.0] ghci> qrSolver (matrix 5 [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]) (col [1,0,8,0,1]) [0.3125000000000003,-1.4665169976505503e-16,-1.8750000000000013,3.5000000000000013, 6.062500000000001]
次は QR 法の基本を簡単に説明します。上 (下) 三角行列の固有値は、対角行列と同様に対角成分と一致します。QR 分解は相似変換ではないので、行列 A を QR 分解したとしても、R の対角成分は A の固有値とは一致しません。そこで、次に示すような相似変換を考えます。
ここで \(A_i\) を QR 分解して、\(A_i = Q_i R_i\) とおけば、上式は次のようになります。
つまり、\(A_i\) を QR 分解して \(Q_i\) と \(R_i\) を求め、それを逆に掛け算したものが \(A_{i+1}\) になります。これを繰り返すと \(A_{i+1}\) の対角成分より下側の成分は 0 に収束することが知られています。つまり、\(A_{i+1}\) は上三角行列になるので、その対角成分が固有値になります。これは相似変換なので、\(A_{i+1}\) の固有値は \(A\) の固有値と一致します。
ただし、ヤコビ法と違って直交行列が固有ベクトルと一致するとは限りません。固有ベクトルを求めるときは他の方法 (たとえば、逆反復法など) を使用することになります。
これをそのままプログラムすると次のようになります。
リスト : QR 法 (ナイーブ版)
qrEig :: Matrix R -> Maybe (Vector R, Int)
qrEig xs =
iter 0 xs (takeDiag xs)
where
iter i a d
| i >= 512 = Nothing
| otherwise =
let
(q, r) = qrGivens a
a1 = r <> q
d1 = takeDiag a1
in
if sumElements (cmap abs ((d1 - d) / d1)) < 1e-8
then Just (d1, i + 1)
else iter (i + 1) a1 d1
対角成分 d1 と d の差分が許容誤差 1e-8 内に入ったら、収束したとして対角成分 (固有値) d1 と繰り返し回数を返します。あとは、QR 法の原理をそのままプログラムしただけなので、特に難しいところはないと思います。
それでは実行してみましょう。
ghci> qrEig $ matrix 2 [2,1,1,3] Just ([3.618033988205326,1.3819660117946744],12) ghci> qrEig $ matrix 2 [2,1,1,2] Just ([2.999999999426406,1.0000000005735947],10) ghci> qrEig $ matrix 3 [1,4,5,4,2,6,5,6,3] Just ([12.175971065046914,-3.6686830908604544,-2.5072879741864527],30) ghci> qrEig $ matrix 4 [6,1,1,1,1,7,1,1,1,1,8,1,1,1,1,9] Just ([10.803886359051218,7.507748637722436,6.392275357812125,5.296089645414204],65)
ナイーブな QR 法の場合、固有値は対角線上に絶対値の大きい順に並びます。結果を見ればお分かりのように、ナイーブな QR 法は収束回数が多くなるという欠点があります。これを改良する方法に「原点移動法」と「減次 (デフレーション)」があります。
QR 法の原点移動法は次式のように行います。
\(A_{i}\) の対角成分から移動量 \(\mu_{i}\) を引き算して QR 分解を行い、\(A_{i+1}\) を求めるとき \(R_{i} @ Q_{i}\) の対角成分に \(\mu_{i}\) を足し算します。移動量 \(\mu\) は、右下隅 2 行 2 列の行列の固有値で、右下の要素に近い値を選びます。これを「ウィルキンソンの移動法」といいます。
プログラムは簡単に修正できるので、さっそく試してみたのですが、ウィルキンソンの移動法だけでは QR 法の収束回数を大幅に減らすことはできませんでした。右下隅の対角成分の収束は速くなるのですが、それ以外の対角成分の収束は速くなりません。
そこで「減次 (デフレーション)」という方法と組み合わせることにします。デフレーションは n * n の行列 \(A\) の固有値 \(\lambda_{n}\) をひとつ求めたら、行列の次数を一つ減らして n-1 * n-1 の行列 \(A'\) で固有値 \(\lambda_{n-1}\) を求めるという方法です。
相似変換で n * n の行列 \(A\) が上三角行列に近づいていくとき、対角成分は固有値に近づいていきます。一番下の行の成分 \(a_{n,0}, \ldots, a_{n,n-1}\) が 0 に収束すると、対角成分 \(a_{n,n}\) は固有値 \(\lambda_{n}\) に収束します。他の非対角成分が固有値 \(\lambda_{n}\) に影響を与えることはありません。残りの固有値 \(\lambda_{1}, \ldots, \lambda_{n-1}\) は n-1 * n-1 の行列 \(A'\) を上三角行列に相似変換すれば求めることができるはずです。
そしてここがポイントですが、最初にウィルキンソンの移動法を適用し、そのあとデフレーションしたときも、その行列での右下隅 2 * 2 行列の固有値を求めてウィルキンソンの移動法を適用します。これで右下隅の対角成分 \(a_{n-1,n-1}\) の収束が速くなります。一番下の行の成分 \(a_{n-1,0}, \ldots, a_{n-1,n-2}\) が 0 に収束したら、対角成分 \(a_{n-1,n-1}\) の値が固有値 \(\lambda_{n-1}\) になります。あとはこれを繰り返すことで、すべての固有値を高速に求めることができます。
それではプログラムを作りましょう。最初に右下隅 (2 行 2 列) の 行列 [[a, b], [c, d]] の固有値を求め、d に近い値を返す関数 eig22 を作ります。これは二次方程式を解くだけなので簡単です。
リスト : シフト値を求める
eig22 :: Matrix R -> R
eig22 xs =
if abs(d - k1) < abs(d - k2) then k1 else k2
where
n = rows xs
a = xs ! (n-2) ! (n-2)
b = xs ! (n-2) ! (n-1)
c = xs ! (n-1) ! (n-2)
d = xs ! (n-1) ! (n-1)
e = a + d
f = sqrt(e * e - 4 * (a * d - b * c))
k1 = (e + f) / 2
k2 = (e - f) / 2
次はウィルキンソンの移動法とデフレーションを行う関数 qrEigShiftd を作ります。
リスト : ウィルキンソンの移動法とデフレーション
qrEigShiftd :: Matrix R -> Maybe ([R], Int)
qrEigShiftd xs =
iter2 0 xs
where
iter1 i xs
| i >= 1024 = Nothing
| otherwise =
let
n = rows xs
k = eig22 xs
dk = (scalar k) * ident n
(q, r) = qrGivens (xs - dk)
xs1 = r <> q + dk
in
if maxElement (cmap abs (subVector 0 (n-1) (xs1 ! (n-1)))) < 1e-14
then Just (xs1 ! (n-1) ! (n-1), xs1, i + 1)
else iter1 (i + 1) xs1
iter2 i xs =
iter1 i xs >>=
\(k, xs1, j) ->
if rows xs1 == 2
then return ([k, xs1 ! 0 ! 0], j)
else iter2 j (xs1 ?? (DropLast 1, DropLast 1)) >>= \(ks, j1) -> return (k:ks, j1)
引数 xs が固有値を求める行列です。実際の処理は局所関数 iter1, iter2 で行います。iter2 の再帰呼び出し (繰り返し) で、行列 xs の次元数を一つずつ減らしていきます。iter1 は再帰呼び出し (繰り返し) で相似変換を行い、求めた固有値、相似変換した行列、繰り返した回数をを返します。
iter2 は行列の次元数が 2 であれば、iter1 で求めた固有値 k ともう一つの固有値 xs1 ! 0 ! 0 をリストに格納して返します。そうでなければ、xs1 の最後の行と列を削除して iter2 を再帰呼び出します。hmatrix の場合、部分行列は演算子 ?? で簡単に求めることができます。
iter1 は最初に繰り返しの回数 (引数 i) をチェックします。1024 以上になったら収束しなかったと判定して Nothing を返します。そうでなければ、シフト値を関数 eig22 で求めて変数 k にセットし、変数 dk に対角成分が k の対角行列をセットします。それから、関数 grGivens で行列 xs - dk を QR 分解して、r <> q + dk で相似変換します。
最後の if 文で、最終行の対角成分以外の絶対値が 1e-14 未満になったかチェックします。ここでは maxElement で最大値を求め、それが 1e-14 未満であれば収束したと判定しています。そうであれば、xs1 の右下隅の成分と xs1 と繰り返し回数 i + 1 を返します。なお、対角成分の左隣の成分 (xs1 ! (n-1) ! (n-2)) だけで収束を判定する方法もあるようです。興味のある方は試してみてください。
プログラムはこれで完成です。簡単な実行例を示します。
ghci> qrEigShiftd $ matrix 2 [2,1,1,3] Just ([3.618033988749895,1.381966011250105],1) ghci> qrEigShiftd $ matrix 2 [2,1,1,2] Just ([1.0,2.9999999999999996],1) ghci> qrEigShiftd $ matrix 3 [1,4,5,4,2,6,5,6,3] Just ([12.175971065046904,-3.6686830979532563,-2.5072879670936308],5) ghci> qrEigShiftd $ matrix 4 [6,1,1,1,1,7,1,1,1,1,8,1,1,1,1,9] Just ([10.803886359051248,7.5077487053636505,6.392275290272985,5.2960896453121205],7)
結果を見ればお分かりのように、ウィルキンソンの移動法とデフレーションを適用することで、QR 法の収束回数を大幅に減少させることができました。ここまで効果が高いとは M.Hiroi も大変驚きました。なお、実用的には行列 A をそのまま QR 分解するのではなく、「三重対角行列」に相似変換してから QR 法を適用したほうが良いようです。これは次回以降に取り上げたいと思います。