ハウスホルダー変換は簡単に言うと「鏡映変換」のことになります。たとえば、ベクトル [x, y] を x 軸で鏡映変換すると [x, -y] になります。また、y 軸で鏡映変換すると [-x, y] になります。このとき、ベクトルの大きさは変わらないことに注意してください。
一般に、ベクトル \(X = [x_1, x_2, \ldots, x_n]\) をベクトル \(Y = [y_1, y_2, \ldots, y_n]\) に鏡映変換するとき、ベクトル \(V \ (= X - Y)\) に垂直で、\(X\) と \(Y\) の角度を二等分する超平面が鏡になります。
Ｘ
／│
／ │
／ │
Ｏ───Ａ── 鏡映面
＼ │
＼ │
＼│
Ｙ
\(V\) の中点を A, 原点を O とすると、O-X-A は直角三角形になるので、X から A までの長さは \(|X|\cos r\), (r は角 O-X-A) になります。内積 \((V, X)\) の値は \(|X||V| \cos r\) になるので、 ベクトル \(Y\) は次式のように表すことができます。
ここで、ベクトル \(X, Y, V\) を n 行 1 列の行列と考えると、 \((V, X)V = V(V, X) = V(V^{\mathrm{T}}X) = (VV^{\mathrm{T}})X\) になるので、\(Y\) は次のようになります。
\(H\) を「ハウスホルダー行列」といい、\(H\) で表される線形変換を「ハウスホルダー変換」といいます。\(H\) は対称行列でかつ直交行列になります。
ハウスホルター変換で特に重要なのが \(X = [x_1, x_2, \ldots, x_n]\) を \(Y = [y_1, 0, \ldots]\), (\(y_1\) は大きさ \(|X|\) で \(x_1\) の符号を反転したもの) に移す変換です。このとき、ベクトル \(V\) は次のようになります。
V は単位ベクトルにしたほうが都合が良いので、c を次のように定義します。
ハウスホルダー行列 \(H\) を作る場合、\(c = \dfrac{1}{\sqrt{|X|(|X| \pm x_1)}}\) とすれば、\(H = I - VV^{\mathrm{T}}\) になります。
それでは実際にプログラムを作ってみましょう。ハウスホルダー行列を生成する関数 houseHolder は次のようになります。
リスト : ハウスホルダー変換
houseHolder :: Vector R -> Matrix R
houseHolder xs =
let
x = norm_2 xs
(y:ys) = toList xs
z = if y < 0 then (y - x) else (y + x)
k = sqrt(x * abs z)
v = z/k : map (/k) ys
in
ident (size xs) - (col v <> row v)
ベクトル xs の大きさを norm_2 で求めて変数 x にセットし、xs をリストに変換して y : ys と照合します。y が負ならば先頭要素 z を y - x に、そうでなければ y + x にします。あとはベクトル V のかわりにリスト v を作り、行列 I - (col v <> row v) を計算して返すだけです。
それでは実行してみましょう。
ghci> h = houseHolder $ vector [0,0,1] ghci> h (3><3) [ 0.0, 0.0, -1.0 , 0.0, 1.0, 0.0 , -1.0, 0.0, 0.0 ] ghci> tr h (3><3) [ 0.0, 0.0, -1.0 , 0.0, 1.0, 0.0 , -1.0, 0.0, 0.0 ] ghci> h <> h (3><3) [ 1.0, 0.0, 0.0 , 0.0, 1.0, 0.0 , 0.0, 0.0, 1.0 ]
houseHolder にベクトル [0, 0, 1] を渡すと、それを [-1, 0, 0] に変換するハウスホルダー行列を返します。ハウスホルダー行列は対称行列かつ直交行列なので、h.T と h は等しくなり、h.T と h-1 は等しくなるので、h * h は単位行列になります。
ところで、行 (または列) を交換する置換行列 R は、行列 A の左側から掛け算 (RA) すると A の行を交換します。右側から掛け算 (AR) すると列を交換します。簡単な例を示しましょう。
ghci> r = matrix 3 [0,0,1,0,1,0,1,0,0] ghci> r (3><3) [ 0.0, 0.0, 1.0 , 0.0, 1.0, 0.0 , 1.0, 0.0, 0.0 ] ghci> x = matrix 3 [0..8] ghci> x (3><3) [ 0.0, 1.0, 2.0 , 3.0, 4.0, 5.0 , 6.0, 7.0, 8.0 ] ghci> r <> x (3><3) [ 6.0, 7.0, 8.0 , 3.0, 4.0, 5.0 , 0.0, 1.0, 2.0 ] ghci> x <> r (3><3) [ 2.0, 1.0, 0.0 , 5.0, 4.0, 3.0 , 8.0, 7.0, 6.0 ]
r は 1 行 (列) と 3 行 (列) を交換する置換行列です。r <> x は 1 行と 3 行を交換し、x <> r は 1 列と 3 列を交換しています。同様に、ハウスホルダー行列 H は行列 A の左側から掛け算 (HA) すると行に対してハウスホルダー変換を行い、右側から掛け算 (AH) すると列に対してハウスホルダー変換を行います。
簡単な例を示しましょう。
ghci> row [0,0,1] <> h (1><3) [ -1.0, 0.0, 0.0 ] ghci> row [-1,0,0] <> h (1><3) [ 0.0, 0.0, 1.0 ] ghci> a = matrix 3 [0,0,1,0,0,1,0,0,1] ghci> a (3><3) [ 0.0, 0.0, 1.0 , 0.0, 0.0, 1.0 , 0.0, 0.0, 1.0 ] ghci> a <> h (3><3) [ -1.0, 0.0, 0.0 , -1.0, 0.0, 0.0 , -1.0, 0.0, 0.0 ] ghci> a <> h <> h (3><3) [ 0.0, 0.0, 1.0 , 0.0, 0.0, 1.0 , 0.0, 0.0, 1.0 ] ghci> b = tr a ghci> b (3><3) [ 0.0, 0.0, 0.0 , 0.0, 0.0, 0.0 , 1.0, 1.0, 1.0 ] ghci> h <> b (3><3) [ -1.0, -1.0, -1.0 , 0.0, 0.0, 0.0 , 0.0, 0.0, 0.0 ] ghci> h <> h <> b (3><3) [ 0.0, 0.0, 0.0 , 0.0, 0.0, 0.0 , 1.0, 1.0, 1.0 ]
ハウスホルダー変換は鏡映なので、[-1, 0, 0] に H を適用すれば [0, 0, 1] に変換されます。つまり、H を二回適用すれば元に戻ります。H * H = I になるので、これは当然ですね。
次はハウスホルダー変換を使って行列を QR 分解するプログラムを作りましょう。基本的な考え方は簡単です。n * n の行列 \(A\) を QR 分解する場合、左端の列が \([x, 0, \ldots, 0]\) になるようにハウスホルダー行列 \(H\) を生成して、それを \(A\) に適用します。
次は、行列 \(A\) の次数を減じて n-1 * n-1 の行列 \(A'\) = \(A\) ?? (Drop 0, Drop 0) を考えて、\(A'\) の左端の列が \([x', 0, \ldots, 0]\) になるようにハウスホルダー行列 \(H'\) を生成します。このとき、単位行列 \(I\) の右下の部分を \(H'\) に置き換えたものが実際のハウスホルダー行列 \(H\) になります。
\(A \times H\) の結果は、1 行目は \(A\) と同じになり、1 列目は \(A\) の 1 列 2 行目以降が 0 なので、\(A \times H\) の 1 列 2 行目以降も 0 になります。この行列も対称行列でかつ直交行列になることに注意してください。あとはこれを繰り返していけば QR 分解を行うことができます。
ハウスホルダー行列の作成には関数 diagBlock を使うと便利です。
diagBlock :: (Element t, Num t) => [Matrix t] -> Matrix t
diagBlock はリストに格納された行列を対角線上に配置した行列を返します。
ghci> diagBlock [ident 2, matrix 2 [1..4], matrix 3 [1..9]] (7><7) [ 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 , 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 , 0.0, 0.0, 1.0, 2.0, 0.0, 0.0, 0.0 , 0.0, 0.0, 3.0, 4.0, 0.0, 0.0, 0.0 , 0.0, 0.0, 0.0, 0.0, 1.0, 2.0, 3.0 , 0.0, 0.0, 0.0, 0.0, 4.0, 5.0, 6.0 , 0.0, 0.0, 0.0, 0.0, 7.0, 8.0, 9.0 ] ghci> diagBlock [ident 2, matrix 2 [1..6], matrix 3 [1..9]] (8><7) [ 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 , 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 , 0.0, 0.0, 1.0, 2.0, 0.0, 0.0, 0.0 , 0.0, 0.0, 3.0, 4.0, 0.0, 0.0, 0.0 , 0.0, 0.0, 5.0, 6.0, 0.0, 0.0, 0.0 , 0.0, 0.0, 0.0, 0.0, 1.0, 2.0, 3.0 , 0.0, 0.0, 0.0, 0.0, 4.0, 5.0, 6.0 , 0.0, 0.0, 0.0, 0.0, 7.0, 8.0, 9.0 ]
簡単な例を示しましょう。
ghci> a = matrix 3 [1,2,3,4,5,6,7,8,10] ghci> a (3><3) [ 1.0, 2.0, 3.0 , 4.0, 5.0, 6.0 , 7.0, 8.0, 10.0 ] ghci> flatten $ takeColumns 1 a [1.0,4.0,7.0] ghci> h = houseHolder $ flatten $ takeColumns 1 a ghci> disp 3 h 3x3 -0.123 -0.492 -0.862 -0.492 0.784 -0.378 -0.862 -0.378 0.339 ghci> a1 = h <> a ghci> disp 3 a1 3x3 -8.124 -9.601 -11.940 0.000 -0.086 -0.550 -0.000 -0.900 -1.462 ghci> b = a1 ?? (Drop 1, Drop 1) ghci> disp 3 b 2x2 -0.086 -0.550 -0.900 -1.462 ghci> h1 = diagBlock[ident 1, houseHolder $ flatten $ takeColumns 1 b] ghci> disp 3 h1 3x3 1.000 0.000 0.000 0.000 -0.095 -0.995 0.000 -0.995 0.095 ghci> disp 3 $ tr h1 3x3 1.000 0.000 0.000 0.000 -0.095 -0.995 0.000 -0.995 0.095 ghci> disp 3 $ h1 <> h1 3x3 1.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 1.000 ghci> a2 = h1 <> a1 ghci> disp 3 $ a2 3x3 -8.124 -9.601 -11.940 0.000 0.905 1.508 -0.000 0.000 0.408 ghci> disp 3 $ h <> h1 3x3 -0.123 0.905 0.408 -0.492 0.302 -0.816 -0.862 -0.302 0.408 ghci> disp 3 $ h <> h1 <> a2 3x3 1.000 2.000 3.000 4.000 5.000 6.000 7.000 8.000 10.000
これをそのままプログラムすると、次のようになります。
リスト : QR 分解 (ハウスホルダー変換)
qrHouse :: Matrix R -> (Matrix R, Matrix R)
qrHouse xs =
foldl (\(ys, xs) i ->
let
h0 = houseHolder $ flatten $ takeColumns 1 (xs ?? (Drop i, Drop i))
h1 = diagBlock [ident i, h0]
in
(ys <> h1, h1 <> xs))
(ident n, xs)
[0 .. n-1]
where
n = rows xs
畳み込みを行う関数 foldl で、行列 xs の次数を一つずつ減らしながらハウスホルダー行列 h1 を生成して、それを ys の右側から掛け算し、xs の左側から掛け算していくだけです。h1 は対称行列でかつ直交行列なので、直交行列 ys を更新するときは h1 を転置する必要はありません。
簡単な実行例を示します。
ghci> a = matrix 2 [2,1,1,3] ghci> a (2><2) [ 2.0, 1.0 , 1.0, 3.0 ] ghci> (q, r) = qrHouse a ghci> disp 3 q 2x2 -0.894 0.447 -0.447 -0.894 ghci> disp 3 r 2x2 -2.236 -2.236 -0.000 -2.236 ghci> disp 3 $ q <> r 2x2 2.000 1.000 1.000 3.000 ghci> b = matrix 2 [2,1,1,2] ghci> b (2><2) [ 2.0, 1.0 , 1.0, 2.0 ] ghci> (q, r) = qrHouse b ghci> disp 3 q 2x2 -0.894 0.447 -0.447 -0.894 ghci> disp 3 r 2x2 -2.236 -1.789 -0.000 -1.342 ghci> disp 3 $ q <> r 2x2 2.000 1.000 1.000 2.000 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) = qrHouse c ghci> disp 3 q 3x3 -0.154 0.802 0.577 -0.617 -0.535 0.577 -0.772 0.267 -0.577 ghci> disp 3 r 3x3 -6.481 -6.481 -6.789 0.000 3.742 1.604 0.000 0.000 4.619 ghci> disp 3 $ q <> r 3x3 1.000 4.000 5.000 4.000 2.000 6.000 5.000 6.000 3.000
正常に動作していますね。ところで、前回作成した固有値を求める関数 qrEigShiftd は QR 分解にギブンス回転 qrGivens を使っていますが、これを qrHouse に変えると実行速度がぐっと速くなります。
ghci> a = (scalar 1) + (diag $ vector [51..100]) ghci> :set +s ghci> qrEigShiftd a Just ... 省略 ... (10.75 secs, 8,663,207,152 bytes) ghci> qrHouseEigShiftd a Just ... 省略 ... (0.51 secs, 516,018,696 bytes) 実行環境 : Ubuntu 22.04 (WSL2), Intel Core i5-6200U 2.30GHz
行列 a のようなテストデータを生成する場合、関数 diagRect を使うと便利です。
diagRect :: Foreign.Storable.Storable t => t -> Vector t -> Int -> Int -> Matrix t
diagRect は対角成分以外の成分を第 1 引数で、対角成分を第 2 引数のベクトルで指定します。そして、行列の大きさを第 3, 4 引数で指定します。簡単な使用例を示します。
ghci> diagRect 1 (vector [4,5,6]) 3 3 (3><3) [ 4.0, 1.0, 1.0 , 1.0, 5.0, 1.0 , 1.0, 1.0, 6.0 ] ghci> diagRect 1 (vector [4,5,6]) 3 4 (3><4) [ 4.0, 1.0, 1.0, 1.0 , 1.0, 5.0, 1.0, 1.0 , 1.0, 1.0, 6.0, 1.0 ] ghci> diagRect 1 (vector [4,5,6]) 4 3 (4><3) [ 4.0, 1.0, 1.0 , 1.0, 5.0, 1.0 , 1.0, 1.0, 6.0 , 1.0, 1.0, 1.0 ] ghci> diagRect 1 (vector [4,5,6]) 4 4 (4><4) [ 4.0, 1.0, 1.0, 1.0 , 1.0, 5.0, 1.0, 1.0 , 1.0, 1.0, 6.0, 1.0 , 1.0, 1.0, 1.0, 1.0 ]
qrHouseEigShiftd は qrGivens を qrHouse に変更しただけの関数です。qrHouse のほうが約 15 倍ちょっと速くなりました。もっとも、これは qrGivens の実装がナイーブだからで、工夫次第では qrGivens でも速くすることが可能です。これは「三重対角行列の QR 分解」で取り上げます。
「三重対角行列 (tridiagonal matrix)」は、対角成分とそれに隣接する成分以外が 0 の行列のことです。簡単な例を下図に示します。
[ a, b, 0, 0, 0 , c, d, e, 0, 0 , 0, f, g, h, 0 , 0, 0, i, j, k , 0, 0, 0, l, m ] 図 : 三重対角行列
実対称行列の固有値を求める場合、三重対角行列に相似変換してから QR 法を適用することが行われます。一般的には、このほうが速くなるといわれています。実対称行列の場合、三重対角化はハウスホルダー変換を使うと簡単です。プログラムは次のようになります。
リスト : 対称行列の三重対角化
tridiag :: Matrix R -> Matrix R
tridiag xs = foldl (\a i ->
let
h0 = houseHolder $ flatten $ takeColumns 1 (a ?? (Drop (i + 1), Drop i))
h1 = diagBlock [ident (i + 1), h0]
in
h1 <> a <> h1)
xs
[0 .. ((rows xs) - 3)]
QR 分解と同じように、行列 a の次数を一つずつ減らしながらハウスホルダー変換を行います。このとき、左上隅の成分はハウスホルダー変換を適用しないことに注意してください。左端の列であれば 2 行目以降の成分に、一番上の行であれば 2 列目以降の成分にハウスホルダー変換を適用します。たとえば、行 (または列) の成分が [x1, x2, ..., xn] とすると、[x2, ..., xn] を [x2', 0, ..., 0] に変換するわけです。
対称行列の性質から、1 行目と 1 列目の成分は同じになるので、ハウスホルダー行列 h は同じものになります。つまり、a の左右に h を掛け算 (h * a * h) すれば、行と列を変換することができます。また、h は対称行列でかつ直交行列なので、h * a * h は相似変換になります。したがって、三重対角化した行列の固有値は元の行列の固有値と一致します。
簡単な実行例を示します。
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> disp 3 $ tridiag c 3x3 1.000 -6.403 0.000 -6.403 8.463 0.829 0.000 0.829 -3.463 ghci> d = diagRect 1 (vector [6..9]) 4 4 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> disp 3 $ tridiag d 4x4 6.000 -1.732 0.000 0.000 -1.732 10.000 0.816 -0.000 0.000 0.816 7.000 -0.577 0.000 -0.000 -0.577 7.000 ghci> e = diagRect 1 (vector [7..11]) 5 5 ghci> e (5><5) [ 7.0, 1.0, 1.0, 1.0, 1.0 , 1.0, 8.0, 1.0, 1.0, 1.0 , 1.0, 1.0, 9.0, 1.0, 1.0 , 1.0, 1.0, 1.0, 10.0, 1.0 , 1.0, 1.0, 1.0, 1.0, 11.0 ] ghci> disp 3 $ tridiag e 5x5 7.000 -2.000 -0.000 0.000 -0.000 -2.000 12.500 -1.118 -0.000 -0.000 -0.000 -1.118 8.500 0.894 -0.000 0.000 -0.000 0.894 8.500 0.671 -0.000 0.000 0.000 0.671 8.500
さきほどのプログラムは、変換行列 h の大きさを行列 a と同じ大きさのままで計算しました。変換行列 h にも減次 (デフレーション) を適用すると、処理を高速化することができます。
行列 a にハウスホルダー変換による相似変換を適用すると、\(a_{0,0}, a_{0,1}, a_{1,0}\) 以外の成分 \(a_{i,0}, a_{0,i}\) は 0 になります。次に相似変換を適用すると \(a_{0,0}, a_{0,1}, a_{1,0}\) の値は変化せず、それ以外の成分 \(a_{i,0}\) と \(a_{0,i}\) は 0 のままです。つまり、1 行目と 1 列目の値は変化しないので、次数を一つ減らした行列 a1 = a ?? (Drop i, Drop i) に相似変換を適用すればいいことになります。これをプログラムすると次のようになります。
リスト : 三重対角化の高速化
tridiag1 :: Matrix R -> Matrix R
tridiag1 xs = foldl (\a i ->
let
h0 = houseHolder $ flatten $ takeColumns 1 (a ?? (Drop (i + 1), Drop i))
h1 = diagBlock [ident 1, h0]
a1 = a ?? (Drop i, Drop i)
in
a ?? (Take i, All) === (a ?? (Drop i, Take i) ||| h1 <> a1 <> h1))
xs
[0 .. ((rows xs) - 3)]
a ?? (Drop i, Drop i) の成分を相似変換の結果に置き換えていくだけです。ハウスホルダー行列の大きさが一つずつ減っていくので、大きさが n のまま演算するよりも効率的です。
それでは実際に試してみましょう。
ghci> :set +s ghci> a = diagRect 1 (vector [101..300]) 200 200 (0.00 secs, 30,952 bytes) ghci> (tridiag a ! 199 ! 199) 199.99999999999594 (4.05 secs, 909,299,288 bytes) ghci> (tridiag1 a ! 199 ! 199) 199.99999999999594 (1.52 secs, 773,984,896 bytes) 実行環境 : Ubuntu 22.04 (WSL2), Intel Core i5-6200U 2.30GHz
200 * 200 の行列で約 2.5 倍ちょっと速くなりました。なお、参考 URL 6『三重対角化 - [物理のかぎしっぽ]』によると、もっと効率の良い方法もあるようです。興味のある方はいろいろ試してみてください。
三重対角行列に変換すると、QR 分解も簡単になります。対角成分 \(a_{i,i}\) の 1 行下の成分 \(a_{i+1,i}\) を 0 にするだけです。ギブンス回転を使うと、プログラムは次のようになります。
リスト : 三重対角行列の QR 分解
qrTri :: Matrix R -> (Matrix R, Matrix R)
qrTri xs =
foldl (\(qs, a) y ->
let
x = y + 1
d = sqrt((a ! x ! y) ** 2 + (a ! y ! y) ** 2)
c = (a ! y ! y) / d
s = (a ! x ! y) / d
a1 = a ?? (Pos (idxs [y]), All)
a2 = a ?? (Pos (idxs [x]), All)
a3 = (scalar c) * a1 + (scalar s) * a2
a4 = (scalar (-s)) * a1 + (scalar c) * a2
q1 = qs ?? (All, Pos (idxs [y]))
q2 = qs ?? (All, Pos (idxs [x]))
q3 = (scalar c) * q1 + (scalar s) * q2
q4 = (scalar (-s)) * q1 + (scalar c) * q2
in
(qs ?? (All, Take y) ||| q3 ||| q4 ||| qs ?? (All, TakeLast (n - 1 - x)),
a ?? (Take y, All) === a3 === a4 === a ?? (TakeLast (n - 1 - x), All)))
(ident n, xs)
[0 .. n-2]
where
n = rows xs
行列 a の (y + 1, y) を 0 にする場合、ギブンス回転で値が変わるのは y, y + 1 行だけ、直交行列 qs では y, y + 1 列だけです。ここだけ計算すれば、処理を高速化することができます。たとえば、4 行 4 列の行列 a で (0, 1) を 0 にする場合は下図のようになります。
s = sin(r), c = cos(r) とする
[[ c, s, 0, 0], [[a11, a12, a13, a14],
[-s, c, 0, 0], * [a21, a22, a23, a24],
[ 0, 0, 1, 0], [a31, a32, a33, a34],
[ 0, 0, 0, 1]] [a41, a42, a43, a44]]
[[ ca11 + sa21, ca12 + sa22, ca13 + sa23, ca14 + sa24],
= [-sa11 + ca21, -sa12 + ca22, -sa13 + ca23, -sa14 + ca24],
[ ... ],
[ ... ]]
[[q11, q12, q13, q14], [[c, -s, 0, 0], [[q11c+q12s, -q11s+q12c, ... ],
[q21, q22, q23, q24], * [s, c, 0, 0], = [q21c+q22s. -q21s+q22c, ... ],
[q31, q32, q33, q34], [0, 0, 1, 0], [q31c+q32s. -q31s+q32c, ... ],
[q41, q42, q43, q44]] [0, 0, 0, 1]] [q41c+q42s. -q41s+q42c, ... ]],
一般に、行列 a の (i+1, i) を 0 にする場合、a の要素と直交行列 q の要素は次のようになります。
ai,j = cos(r) * ai,j + sin(r) * ai+1,j ai+1,j = -sin(r) * ai,j + cos(r) * ai+1,j qj,i = cos(r) * qj,i + sin(r) * qj,i+1 qj,i+1 = -sin(r) * qj,i + cos(r) * qj,i+1 j = 0, 1, ..., n-1 (n は行列 a の次元数)
関数 qrTri は、これをそのままプログラムしています。ただし、Matrix は immutable なデータ型なので、配列の要素を書き換えることはできません。このため、プログラムはちょっと複雑になってしまいました。mutable な配列を使ったほうが効率よくプログラムできるかもしれません。実際に NUmPy や Julia で作成したプログラムでは効果がありました。興味のある方はプログラムを改良してみてください。
簡単な実行例を示します。
ghci> a = diagRect 1 (vector [6,7,8,9,10]) 5 5 ghci> a (5><5) [ 6.0, 1.0, 1.0, 1.0, 1.0 , 1.0, 7.0, 1.0, 1.0, 1.0 , 1.0, 1.0, 8.0, 1.0, 1.0 , 1.0, 1.0, 1.0, 9.0, 1.0 , 1.0, 1.0, 1.0, 1.0, 10.0 ] ghci> b = tridiag1 a ghci> disp 3 b 5x5 6.000 -2.000 -0.000 -0.000 -0.000 -2.000 11.500 -1.118 0.000 0.000 -0.000 -1.118 7.500 0.894 0.000 -0.000 0.000 0.894 7.500 0.671 -0.000 0.000 0.000 0.671 7.500 ghci> (q, r) = qrTri b ghci> disp 3 q 5x5 0.949 0.314 0.034 -0.004 0.000 -0.316 0.943 0.102 -0.012 0.001 0.000 -0.108 0.987 -0.120 0.011 0.000 0.000 0.121 0.989 -0.090 0.000 0.000 0.000 0.091 0.996 ghci> disp 3 r 5x5 6.325 -5.534 0.354 -0.000 -0.000 0.000 10.338 -1.866 -0.097 0.000 -0.000 0.000 7.396 1.790 0.081 -0.000 0.000 0.000 7.368 1.346 -0.000 0.000 0.000 0.000 7.408 ghci> disp 3 $ q <> r 5x5 6.000 -2.000 -0.000 -0.000 -0.000 -2.000 11.500 -1.118 0.000 0.000 -0.000 -1.118 7.500 0.894 0.000 -0.000 0.000 0.894 7.500 0.671 -0.000 0.000 0.000 0.671 7.500
最後に、三重対角化による QR 法のプログラムを作ります。まず前提として、三重対角行列を QR 分解して相似変換 (R * Q) する場合、三重対角行列の性質は残ることに注意してください。つまり、主対角線と副対角線以外は 0 のままで、相似変換を繰り返すと副対角線が 0 に近づいていくことになります。したがって、QR 分解は qrTri で大丈夫です。収束判定も対角成分 \(a_{i,i}\) の左隣 \(a_{i.i-1}\) が 0 になったかチェックすれば OK です。
プログラムは次のようになります。
リスト : 三重対角化による QR 法
qrEigTri :: Matrix R -> Maybe ([R], Int)
qrEigTri xs =
iter2 0 $ tridiag1 xs
where
iter1 i xs
| i >= 1024 = Nothing
| otherwise =
let
n = rows xs
k = eig22 xs
dk = (scalar k) * ident n
(q, r) = qrTri (xs - dk)
xs1 = r <> q + dk
in
if abs (xs1 ! (n-1) ! (n-2)) < 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)
前回作成した関数 qrEigShiftd とほぼ同じですが、最初に引数の行列 xs を tridiag1 で三重対角化するところと、QR 分解に関数 qrTri を使うところが異なります。QR 分解は繰り返し (再帰呼び出し) の中で行っているので、三重対角行列専用の関数 qrTri を使うことで高速化が期待できます。これはあとで試してみましょう。
簡単な実行例を示します。
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> qrEigTri c Just ([-3.668683097953265,-2.5072879670936414,12.175971065046898],4) ghci> d = diagRect 1 (vector [6,7,8,9]) 4 4 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> qrEigTri d Just ([6.39227529027298,7.507748705363651,5.29608964531212,10.803886359051257],6) ghci> e = diagRect 1 (vector [7 .. 11]) 5 5 ghci> e (5><5) [ 7.0, 1.0, 1.0, 1.0, 1.0 , 1.0, 8.0, 1.0, 1.0, 1.0 , 1.0, 1.0, 9.0, 1.0, 1.0 , 1.0, 1.0, 1.0, 10.0, 1.0 , 1.0, 1.0, 1.0, 1.0, 11.0 ] ghci> qrEigTri e Just ([8.434736666495787,7.3566318548442196,9.540394425688136,6.277695819922921, 13.39054123304896],10)
正常に動作していますね。最後に、qrHouseEigShiftd と qrEigTri の実行時間を比較してみましょう。100 * 100 の行列の固有値を求めます。結果は次のようになりました。
ghci> a = (scalar 1) + diag (vector [101 .. 200]) ghci> :set +s ghci> qrHouseEigShiftd a Just ... 省略 ... (11.42 secs, 5,977,575,896 bytes) ghci> qrEigTri a Just ... 省略 ... (6.14 secs, 7,757,437,960 bytes) 実行環境 : Ubuntu 22.04 (WSL2), Intel Core i5-6200U 2.30GHz
qrEigTri の方が速くなりました。NumPy や Julia で作成したプログラムでは 150 * 150 の行列で約 20 倍ほど高速化できたのですが、そこまでの効果はありませんでした。プログラムをコンパイルして、もう少し大きな行列で試してみるとよいかもしれません。ただし、メモリをかなり消費しているので、immutable な Matrix ではこれ以上の高速化は難しいのではないか、と思っています。興味のある方はいろいろ試してみてください。
--
-- qr.hs : QR 分解と QR 法
--
-- Copyright (C) 2021 Makoto Hiroi
--
import Prelude hiding ((<>))
import Numeric.LinearAlgebra
-- ギブンス回転
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
-- 後退代入
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
--
-- 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
-- eig22 xs => 右下の 2 * 2 行列 [[a, b], [c, d]]
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 :: 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)
-- ハウスホルダー変換
houseHolder :: Vector R -> Matrix R
houseHolder xs =
let
x = norm_2 xs
(y:ys) = toList xs
z = if y < 0 then (y - x) else (y + x)
k = sqrt(x * abs z)
v = z/k : map (/k) ys
in
ident (size xs) - (col v <> row v)
qrHouse :: Matrix R -> (Matrix R, Matrix R)
qrHouse xs =
foldl (\(ys, xs) i ->
let
h0 = houseHolder $ flatten $ takeColumns 1 (xs ?? (Drop i, Drop i))
h1 = diagBlock [ident i, h0]
in
(ys <> h1, h1 <> xs))
(ident n, xs)
[0 .. n-1]
where
n = rows xs
qrHouseEigShiftd :: Matrix R -> Maybe ([R], Int)
qrHouseEigShiftd 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) = qrHouse (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)
-- 対称行列の三重化
tridiag :: Matrix R -> Matrix R
tridiag xs = foldl (\a i ->
let
h0 = houseHolder $ flatten $ takeColumns 1 (a ?? (Drop (i + 1), Drop i))
h1 = diagBlock [ident (i + 1), h0]
in
h1 <> a <> h1)
xs
[0 .. ((rows xs) - 3)]
-- 高速化
tridiag1 :: Matrix R -> Matrix R
tridiag1 xs = foldl (\a i ->
let
h0 = houseHolder $ flatten $ takeColumns 1 (a ?? (Drop (i + 1), Drop i))
h1 = diagBlock [ident 1, h0]
a1 = a ?? (Drop i, Drop i)
in
a ?? (Take i, All) === (a ?? (Drop i, Take i) ||| h1 <> a1 <> h1))
xs
[0 .. ((rows xs) - 3)]
-- 三重対角行列の QR 分解
qrTri :: Matrix R -> (Matrix R, Matrix R)
qrTri xs =
foldl (\(qs, a) y ->
let
x = y + 1
d = sqrt((a ! x ! y) ** 2 + (a ! y ! y) ** 2)
c = (a ! y ! y) / d
s = (a ! x ! y) / d
a1 = a ?? (Pos (idxs [y]), All)
a2 = a ?? (Pos (idxs [x]), All)
a3 = (scalar c) * a1 + (scalar s) * a2
a4 = (scalar (-s)) * a1 + (scalar c) * a2
q1 = qs ?? (All, Pos (idxs [y]))
q2 = qs ?? (All, Pos (idxs [x]))
q3 = (scalar c) * q1 + (scalar s) * q2
q4 = (scalar (-s)) * q1 + (scalar c) * q2
in
(qs ?? (All, Take y) ||| q3 ||| q4 ||| qs ?? (All, TakeLast (n - 1 - x)),
a ?? (Take y, All) === a3 === a4 === a ?? (TakeLast (n - 1 - x), All)))
(ident n, xs)
[0 .. n-2]
where
n = rows xs
-- 三重対角化による QR 法
qrEigTri :: Matrix R -> Maybe ([R], Int)
qrEigTri xs =
iter2 0 $ tridiag1 xs
where
iter1 i xs
| i >= 1024 = Nothing
| otherwise =
let
n = rows xs
k = eig22 xs
dk = (scalar k) * ident n
(q, r) = qrTri (xs - dk)
xs1 = r <> q + dk
in
if abs (xs1 ! (n-1) ! (n-2)) < 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)