●QR 分解と QR 法 (2)
前回は QR 分解と QR 法の基本を説明し、ギブンス回転を使った QR 分解と QR 法のプログラムを作りました。今回はギブンス回転よりも効率が良いといわれている「ハウスホルダー変換 (Householder transformation)」を使って、QR 分解と QR 法を実装してみましょう。
●ハウスホルダー変換の基本
ハウスホルダー変換は簡単に言うと「鏡映変換」のことになります。たとえば、ベクトル [x, y] を x 軸で鏡映変換すると [x, -y] になります。また、y 軸で鏡映変換すると [-x, y] になります。このとき、ベクトルの長さは変わらないことに注意してください。
一般に、ベクトル X = [x1, x2, ..., xn] をベクトル Y = [y1, y2, ..., yn] に鏡映変換するとき、ベクトル 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 は次式のように表すことができます。
Y = X - V = X - 2(V, X)V / |V|2
ここで、ベクトル X, Y, V を n 行 1 列の行列と考えると、 (V, X)V = V(V, X) = V(VTX) = (VVT)X になるので、Y は次のようになります。
Y = HX, H = I - 2VVT / |V|2
この行列 H を「ハウスホルダー行列」といい、H で表される線形変換を「ハウスホルダー変換」といいます。H は対称行列でかつ直交行列になります。
ハウスホルター変換で特に重要なのが X = [x1, x2, ..., xn] を Y = [y1, 0, ...], (y1 は大きさ |X| で x1 の符号を反転したもの) に移す変換です。このとき、ベクトル V は次のようになります。
V = c * [x1±|X|, x2, ..., xn], c != 0, ±|X| は x1 と同じ符号
V は単位ベクトルにしたほうが都合が良いので、c を次のように定義します。
|V|2 = c2 * ((x1±|X|)2 + x22 + ... + xn2) = 2c2|X|(|X|±x1) c = 1 / √(2|X|(|X|±x1)) H = I - 2VVT
ハウスホルダー行列 H を作る場合、c = 1 / √(|X|(|X|±x1)) とすれば、H = I - VVT になります。
●ハウスホルダー変換のプログラム
それではプログラムを作りましょう。ハウスホルダー行列を生成する関数 householer() は次のようになります。
リスト : ハウスホルダー変換
def householder(xs):
n = len(xs)
x = np.linalg.norm(xs)
if xs[0] < 0: x = -x
vs = xs.astype(np.float_)
vs[0] += x
vs /= np.sqrt(vs[0] * x)
return np.eye(n) - vs.reshape((n, 1)) @ vs.reshape((1, n))
ベクトル xs の大きさを linalg.norm() で求めて変数 x にセットします。xs[0] が負ならば x も負にします。あとはベクトル vs を作って、行列 I - vsvsT を作成して返すだけです。
それでは実行してみましょう。
>>> h = householder(np.array([0, 0, 1]))
>>> h
array([[ 0., 0., -1.],
[ 0., 1., 0.],
[-1., 0., 0.]])
>>> h.T
array([[ 0., 0., -1.],
[ 0., 1., 0.],
[-1., 0., 0.]])
>>> h @ h
array([[ 1., 0., 0.],
[ 0., 1., 0.],
[ 0., 0., 1.]])
householder() にベクトル [0, 0, 1] を渡すと、それを [-1, 0, 0] に変換するハウスホルダー行列を返します。ハウスホルダー行列は対称行列かつ直交行列なので、h.T と h は等しくなり、h.T と h-1 は等しくなるので、h @ h は単位行列になります。
ところで、行 (または列) を交換する置換行列 R は、行列 A の左側から掛け算 (RA) すると A の行を交換します。右側から掛け算 (AR) すると列を交換します。簡単な例を示しましょう。
>>> r = np.array([[0, 0, 1], [0, 1, 0], [1, 0, 0]])
>>> x = np.arange(9).reshape((3,3))
>>> x
array([[0, 1, 2],
[3, 4, 5],
[6, 7, 8]])
>>> r @ x
array([[6, 7, 8],
[3, 4, 5],
[0, 1, 2]])
>>> x @ r
array([[2, 1, 0],
[5, 4, 3],
[8, 7, 6]])
r は 1 行 (列) と 3 行 (列) を交換する置換行列です。r @ x は 1 行と 3 行を交換し、x @ r は 1 列と 3 列を交換しています。同様に、ハウスホルダー行列 H は行列 A の左側から掛け算 (HA) すると行に対してハウスホルダー変換を行い、右側から掛け算 (AH) すると列に対してハウスホルダー変換を行います。
簡単な例を示しましょう。
>>> np.array([[0, 0, 1]]) @ h
array([[-1., 0., 0.]])
>>> a = np.array([[0, 0, 1], [0, 0, 1], [0, 0, 1]])
>>> a
array([[0, 0, 1],
[0, 0, 1],
[0, 0, 1]])
>>> a @ h
array([[-1., 0., 0.],
[-1., 0., 0.],
[-1., 0., 0.]])
>>> b = a.T
>>> b
array([[0, 0, 0],
[0, 0, 0],
[1, 1, 1]])
>>> h @ b
array([[-1., -1., -1.],
[ 0., 0., 0.],
[ 0., 0., 0.]])
ハウスホルダー変換は鏡映なので、[-1, 0, 0] に H を適用すれば [0, 0, 1] に変換されます。つまり、H を二回適用すれば元に戻ります。H @ H = I になるので、これは当然ですね。
>>> a
array([[0, 0, 1],
[0, 0, 1],
[0, 0, 1]])
>>> a @ h
array([[-1., 0., 0.],
[-1., 0., 0.],
[-1., 0., 0.]])
>>> (a @ h) @ h
array([[ 0., 0., 1.],
[ 0., 0., 1.],
[ 0., 0., 1.]])
>>> b
array([[0, 0, 0],
[0, 0, 0],
[1, 1, 1]])
>>> h @ b
array([[-1., -1., -1.],
[ 0., 0., 0.],
[ 0., 0., 0.]])
>>> h @ (h @ b)
array([[ 0., 0., 0.],
[ 0., 0., 0.],
[ 1., 1., 1.]])
●ハウスホルダー変換による QR 分解
次はハウスホルダー変換を使って行列を QR 分解するプログラムを作りましょう。基本的な考え方は簡単です。n * n の行列 A を QR 分解する場合、左端の列が [x, 0, ..., 0] になるようにハウスホルダー行列 H を生成して、それを A に適用します。
次は、行列の次数を減じて n-1 * n-1 の行列 A' = A[1:, 1:] を考えて、A' の左端の列が [x', 0, ..., 0] になるようにハウスホルダー行列 H' を生成します。このとき、単位行列 I の I[1:, 1:] の部分を H' に置き換えたものが実際の変換行列 H になります。
A @ H の結果は、1 行目は A と同じになり、1 列目は A の 1 列 2 行目以降が 0 なので、A @ H の 1 列 2 行目以降も 0 になります。この行列も対称行列でかつ直交行列になることに注意してください。あとはこれを繰り返していけば QR 分解を行うことができます。
簡単な例を示しましょう。
>>> a = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 10]])
>>> a
array([[ 1, 2, 3],
[ 4, 5, 6],
[ 7, 8, 10]])
>>> h = householder(a[:,0])
>>> h
array([[-0.12309149, -0.49236596, -0.86164044],
[-0.49236596, 0.7841456 , -0.3777452 ],
[-0.86164044, -0.3777452 , 0.33894589]])
>>> a1 = h @ a
>>> a1
array([[ -8.12403840e+00, -9.60113630e+00, -1.19398746e+01],
[ 4.44089210e-16, -8.59655700e-02, -5.49676344e-01],
[ -4.44089210e-16, -9.00439748e-01, -1.46193360e+00]])
>>> h1 = np.eye(3)
>>> h1[1:,1:] = householder(a1[1:, 1])
>>> h1
array([[ 1. , 0. , 0. ],
[ 0. , -0.09503851, -0.9954736 ],
[ 0. , -0.9954736 , 0.09503851]])
>>> h1.T
array([[ 1. , 0. , 0. ],
[ 0. , -0.09503851, -0.9954736 ],
[ 0. , -0.9954736 , 0.09503851]])
>>> h1 @ h1
array([[ 1.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, 1.00000000e+00, 1.11022302e-16],
[ 0.00000000e+00, 1.11022302e-16, 1.00000000e+00]])
>>> a2 = h1 @ a1
>>> a2
array([[ -8.12403840e+00, -9.60113630e+00, -1.19398746e+01],
[ 3.99873504e-16, 9.04534034e-01, 1.50755672e+00],
[ -4.84284661e-16, -1.24900090e-16, 4.08248290e-01]])
>>> h @ h1
array([[-0.12309149, 0.90453403, 0.40824829],
[-0.49236596, 0.30151134, -0.81649658],
[-0.86164044, -0.30151134, 0.40824829]])
>>> h @ h1 @ a2
array([[ 1., 2., 3.],
[ 4., 5., 6.],
[ 7., 8., 10.]])
これをそのままプログラムすると、次のようになります。
リスト : QR 分解 (ハウスホルダー変換)
def qr_h(a):
n = len(a)
xs = np.eye(n)
for i in range(n - 1):
h = np.eye(n)
h[i:, i:] = householder(a[i:, i])
a = h @ a
xs = xs @ h
return xs, a
for ループで行列 a の次数を一つずつ減らしながら変換行列 h を生成して、それを a の左側から掛け算していくだけです。h は対称行列でかつ直交行列なので、直交行列 xs を更新するときは h を転置する必要はありません。
簡単な実行例を示します。
>>> a
array([[2, 1],
[1, 3]])
>>> q, r = qr_h(a)
>>> q
array([[-0.89442719, -0.4472136 ],
[-0.4472136 , 0.89442719]])
>>> r
array([[ -2.23606798e+00, -2.23606798e+00],
[ 2.22044605e-16, 2.23606798e+00]])
>>> q @ r
array([[ 2., 1.],
[ 1., 3.]])
>>> b
array([[2, 1],
[1, 2]])
>>> q, r = qr_h(b)
>>> q
array([[-0.89442719, -0.4472136 ],
[-0.4472136 , 0.89442719]])
>>> r
array([[ -2.23606798e+00, -1.78885438e+00],
[ 2.22044605e-16, 1.34164079e+00]])
>>> q @ r
array([[ 2., 1.],
[ 1., 2.]])
>>> c
array([[1, 4, 5],
[4, 2, 6],
[5, 6, 3]])
>>> q, r = qr_h(c)
>>> q
array([[-0.15430335, 0.80178373, -0.57735027],
[-0.6172134 , -0.53452248, -0.57735027],
[-0.77151675, 0.26726124, 0.57735027]])
>>> r
array([[ -6.48074070e+00, -6.48074070e+00, -6.78934740e+00],
[ 5.47064065e-16, 3.74165739e+00, 1.60356745e+00],
[ -3.08466142e-16, 0.00000000e+00, -4.61880215e+00]])
>>> 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 = qr_h(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.24499800e+00, -2.40192231e+00, -2.56205046e+00, -2.72217861e+00],
[ -1.04432106e-16, -6.79932123e+00, -1.59518185e+00, -1.68568863e+00],
[ -9.46397882e-17, -2.77402609e-16, -7.60863275e+00, -1.22711418e+00],
[ 9.39872435e-17, 2.01276466e-16, 0.00000000e+00, 8.49955224e+00]])
>>> 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 = qr_h(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.28010989e+00, -2.47249015e+00, -2.60985072e+00, -2.74721128e+00, -2.88457184e+00],
[ -1.96834211e-16, -7.86681590e+00, -1.72206519e+00, -1.80600987e+00, -1.88995455e+00],
[ -1.80645774e-16, 7.01728601e-17, -8.67312924e+00, -1.35131411e+00, -1.40861191e+00],
[ -1.78196962e-16, 8.50288251e-17, -2.40796181e-16, -9.55850975e+00, -1.12556437e+00],
[ 1.83667216e-16, -1.01765000e-16, 2.08700803e-16, 1.11022302e-16, 1.04812325e+01]])
>>> 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.]])
正常に動作していますね。ところで、前回作成した固有値を求める関数 qr_eig_shiftd() は QR 分解にギブンス回転 qr_g() を使っていますが、これを qr_h() に変えると実行速度がぐっと速くなります。
ギブンス回転 (qr_g() を使用)
>>> s = time.time(); qr_eig_shiftd(np.ones((50,50)) + np.diag(np.arange(51, 101))); print(time.time() - s)
array([ 99.74906781, 98.69950593, 97.66614768, 96.63999547,
95.61808475, 94.59902629, 93.58204187, 92.56664498,
・・・省略・・・
55.24442862, 54.23218455, 53.21782509, 52.19967428,
51.17236607, 129.59687693])
1.1123793125152588
ハウスホルダー変換 (qr_h() を使用)
>>> s = time.time(); qr_eig_shiftd(np.ones((50,50)) + np.diag(np.arange(51, 101))); print(time.time() - s)
array([ 99.74906781, 98.69950593, 97.66614768, 96.63999547,
95.61808475, 94.59902629, 93.58204187, 92.56664498,
・・・省略・・・
55.24442862, 54.23218455, 53.21782509, 52.19967428,
51.17236607, 129.59687693])
0.1171114444732666
- Ubunts 22.04 LTS (WSL2 + WSLg, Windows 10), Intel Core i5-6200U 2.30GHz
gr_h() のほうが 9.5 倍も速くなりました。もっとも、これは qr_g() の実装がナイーブだからで、工夫次第では qr_g() でも速くすることが可能です。これは「三重対角行列の 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 法を適用することが行われます。一般的には、このほうが速くなるといわれています。実対称行列の場合、三重対角化はハウスホルダー変換を使うと簡単です。プログラムは次のようになります。
リスト : 対称行列の三重対角化
def tridiag(a):
n = len(a)
for i in range(n - 2):
h = np.eye(n)
h[i+1:, i+1:] = householder(a[i+1:, i])
a = h @ a @ h
return a
QR 分解と同じように、行列 a の次数を一つずつ減らしながらハウスホルダー変換を行います。このとき、左上隅の成分はハウスホルダー変換を適用しないことに注意してください。左端の列であれば 2 行目以降の成分に、一番上の行であれば 2 列目以降の成分にハウスホルダー変換を適用します。たとえば、行 (または列) の成分が [x1, x2, ..., xn] とすると、[x2, ..., xn] を [x2', 0, ..., 0] に変換するわけです。
対称行列の性質から、1 行目と 1 列目の成分は同じになるので、ハウスホルダー行列 h は同じものになります。つまり、a の左右に h を掛け算 (h @ a @ h) すれば、行と列を変換することができます。また、h は対称行列でかつ直交行列なので、h @ a @ h は相似変換になります。したがって、三重対角化した行列の固有値は元の行列の固有値と一致します。
簡単な実行例を示します。
>>> c
array([[1, 4, 5],
[4, 2, 6],
[5, 6, 3]])
>>> tridiag(c)
array([[ 1.00000000e+00, -6.40312424e+00, 8.88178420e-16],
[ -6.40312424e+00, 8.46341463e+00, 8.29268293e-01],
[ 8.88178420e-16, 8.29268293e-01, -3.46341463e+00]])
>>> d
array([[ 6., 1., 1., 1.],
[ 1., 7., 1., 1.],
[ 1., 1., 8., 1.],
[ 1., 1., 1., 9.]])
>>> tridiag(d)
array([[ 6.00000000e+00, -1.73205081e+00, 1.43157667e-16, 1.05314450e-16],
[ -1.73205081e+00, 1.00000000e+01, 8.16496581e-01, 5.82867088e-16],
[ 1.43157667e-16, 8.16496581e-01, 7.00000000e+00, -5.77350269e-01],
[ 1.05314450e-16, 2.49800181e-16, -5.77350269e-01, 7.00000000e+00]])
>>> 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.]])
>>> tridiag(e)
array([[ 7.00000000e+00, -2.00000000e+00, -4.34443477e-16, 2.77555756e-16, -1.24126708e-16],
[ -2.00000000e+00, 1.25000000e+01, -1.11803399e+00, 6.01479669e-16, 1.17080690e-15],
[ -4.34443477e-16, -1.11803399e+00, 8.50000000e+00, 8.94427191e-01, 5.27355937e-16],
[ 2.77555756e-16, -1.94695753e-16, 8.94427191e-01, 8.50000000e+00, 6.70820393e-01],
[ -1.24126708e-16, 2.92152494e-16, 8.32667268e-17, 6.70820393e-01, 8.50000000e+00]])
●三重対角化の高速化
さきほどのプログラムは、変換行列 h の大きさを行列 a と同じ大きさのままで計算しました。変換行列 h にも減次 (デフレーション) を適用すると、処理を高速化することができます。
行列 a にハウスホルダー変換による相似変換を適用すると、a0,0, a0,1, a1,0 以外の成分 ai,0, a0,i は 0 になります。次に相似変換を適用すると a0,0, a0,1, a1,0 の値は変化せず、それ以外の成分 ai,0 と a0,i は 0 のままです。つまり、1 行目と 1 列目の値は変化しないので、次数を一つ減らした行列 a' = a[1:, 1:] に相似変換を適用すればいいことになります。これをプログラムすると次のようになります。
リスト : 三重対角化の高速化
def tridiag1(a):
n = len(a)
a = a.astype(np.float_)
for i in range(n - 2):
h = np.eye(n - i)
h[1:, 1:] = householder(a[i+1:, i])
a[i:, i:] = h @ a[i:, i:] @ h
return a
最初に行列 a の要素の型を float_ に変換します。要素が型が整数のままだと、a に値を代入するとき、整数に型変換されてしまいます。あとは a[i:, i:] の相似変換の結果を a[i:, i:] に代入していくだけです。行列 h の大きさが一つずつ減っていくので、大きさが n のまま演算するよりも効率的です。
それでは実際に試してみましょう。
>>> s = time.time(); tridiag(np.ones((100,100)) + np.diag(np.arange(101, 201))); print(time.time() - s)
array([[ 1.02000000e+02, -9.94987437e+00, 1.52418769e-15, ...,
1.03216015e-16, 7.87676636e-17, -3.72077168e-17],
[ -9.94987437e+00, 2.50000000e+02, -2.85773803e+01, ...,
9.24669584e-16, -3.06356128e-15, 3.92078764e-16],
[ 1.52418769e-15, -2.85773803e+01, 1.51000000e+02, ...,
6.81594520e-17, 8.49211710e-16, -8.33618831e-16],
...,
[ 1.03216015e-16, -1.39857673e-16, -1.56427235e-16, ...,
1.51000000e+02, -4.94981323e+00, 1.50990331e-14],
[ 7.87676636e-17, 1.80999785e-16, 1.69716002e-16, ...,
-4.94981323e+00, 1.51000000e+02, 3.50896288e+00],
[ -3.72077168e-17, 1.08778989e-16, 4.47843865e-16, ...,
-2.22044605e-16, 3.50896288e+00, 1.51000000e+02]])
0.034636497497558594
>>> s = time.time(); tridiag1(np.ones((100,100)) + np.diag(np.arange(101, 201))); print(time.time() - s)
array([[ 1.02000000e+02, -9.94987437e+00, 8.56953397e-16, ...,
9.67975700e-16, 9.83588211e-16, 9.99200722e-16],
[ -9.94987437e+00, 2.50000000e+02, -2.85773803e+01, ...,
1.08524301e-14, 1.06858966e-14, 1.24344979e-14],
[ 8.56953397e-16, -2.85773803e+01, 1.51000000e+02, ...,
-3.02813330e-14, -5.99520433e-15, 1.77635684e-15],
...,
[ 9.67975700e-16, 2.77555756e-16, 2.44249065e-15, ...,
1.51000000e+02, -4.94981323e+00, 1.50990331e-14],
[ 9.83588211e-16, 1.49880108e-15, 2.41473508e-15, ...,
-4.94981323e+00, 1.51000000e+02, 3.50896288e+00],
[ 9.99200722e-16, -8.88178420e-16, 0.00000000e+00, ...,
-2.22044605e-16, 3.50896288e+00, 1.51000000e+02]])
0.019008874893188477
100 * 100 の行列で約 1.8 倍ちょっと速くなりました。なお、参考 URL 6 によると、もっと効率の良い方法もあるようです。興味のある方はいろいろ試してみてください。
●三重対角行列の QR 分解
三重対角行列に変換すると、QR 分解も簡単になります。対角成分 ai,i の 1 行下の成分 ai+1,i を 0 にするだけです。ギブンス回転を使うと、プログラムは次のようになります。
リスト : 三重対角行列の QR 分解
def qr_tri(a):
n = len(a)
qs = np.eye(n) # 直交行列
for y in range(n - 1):
x = y + 1
d = np.sqrt(a[x, y] ** 2 + a[y, y] ** 2)
c = a[y, y] / d
s = a[x, y] / d
a1 = a[y].copy()
a2 = a[x]
a[y] = a1 * c + a2 * s
a[x] = a1 * (-s) + a2 * c
q1 = qs[:, y].copy()
q2 = qs[:, x]
qs[:, y] = q1 * c + q2 * s
qs[:, x] = q1 * (-s) + q2 * c
return qs, a
ギブンス回転で行列を QR 分解する関数 qr_g() は二重の for ループになりますが、三重対角行列は一重の for ループで済みます。それから、a[y + 1, y] を 0 にする場合、ギブンス回転で値が変わるのは y, y + 1 行だけ、直交行列 qs では y, y + 1 列だけです。ここだけ計算すれば、処理をさらに高速化することができます。たとえば、4 行 4 列の行列 a で a[0, 1] を 0 にする場合は下図のようになります。
s = sin(r), c = cos(r) とする [[ c, s, 0, 0], [[a11, a12, a13, a14], [[ ca11 + sa21, ca12 + sa22, ca13 + sa23, ca14 + sa24], [-s, c, 0, 0], @ [a21, a22, a23, a24], = [-sa11 + ca21, -sa12 + ca22, -sa13 + ca23, -sa14 + ca24], [ 0, 0, 1, 0], [a31, a32, a33, a34], [ ... ], [ 0, 0, 0, 1]] [a41, a42, a43, a44]] [ ... ]] [[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 の次元数)
関数 qr_tri() は、これをそのままプログラムしています。ただし、a[y], q[:, y] を書き換えると、a[x], q[: x] (x = y + 1) を計算できなくなるので、copy() で a[y], q[:, y] をコピーしていることに注意してください。なお、三重対角行列の特徴を使えば、さらに計算量を減らすことができると思います。興味のある方はプログラムを改良してみてください。
簡単な実行例を示します。
>>> x = np.ones((5, 5))
>>> x1 = x - np.triu(x, 2) - np.tril(x, -2)
>>> x1
array([[ 1., 1., 0., 0., 0.],
[ 1., 1., 1., 0., 0.],
[ 0., 1., 1., 1., 0.],
[ 0., 0., 1., 1., 1.],
[ 0., 0., 0., 1., 1.]])
>>> q, r = qr_tri(x1)
>>> q
array([[ 0.70710678, 0. , -0.40824829, 0.28867513, 0.5 ],
[ 0.70710678, 0. , 0.40824829, -0.28867513, -0.5 ],
[ 0. , 1. , 0. , 0. , 0. ],
[ 0. , 0. , 0.81649658, 0.28867513, 0.5 ],
[ 0. , 0. , 0. , 0.8660254 , -0.5 ]])
>>> r
array([[ 1.41421356e+00, 1.41421356e+00, 7.07106781e-01, 0.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 1.22474487e+00, 8.16496581e-01, 8.16496581e-01],
[ 0.00000000e+00, 0.00000000e+00, -5.55111512e-17, 1.15470054e+00, 1.15470054e+00],
[ 0.00000000e+00, 0.00000000e+00, -9.61481343e-17, 0.00000000e+00, 0.00000000e+00]])
>>> q @ r
array([[ 1.00000000e+00, 1.00000000e+00, -1.75121059e-16, 5.55111512e-17, 5.55111512e-17],
[ 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, -5.55111512e-17, -5.55111512e-17],
[ 0.00000000e+00, 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 1.00000000e+00, 1.00000000e+00, 1.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00, 1.00000000e+00]])
●三重対角行列による QR 法
最後に、三重対角化による QR 法のプログラムを作ります。まず前提として、三重対角行列を QR 分解して相似変換 (R @ Q) する場合、三重対角行列の性質は残ることに注意してください。つまり、主対角線と副対角線以外は 0 のままで、相似変換を繰り返すと副対角線が 0 に近づいていくことになります。したがって、QR 分解は qr_tri() で大丈夫です。収束判定も対角成分 ai,i の左隣 ai.i-1 が 0 になったかチェックすれば OK です。
プログラムは次のようになります。
リスト : 三重対角化による QR 法
def qr_eig_tri(a, max_iter = 1024, e = 1e-14):
a = tridiag1(a.astype(np.float_))
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_tri(a - dk)
a = r @ q + dk
#print(i, size, k)
#print(a)
i += 1
if np.fabs(a[-1, -2]) < e:
ks.append(a[-1, -1])
break
else:
raise Exception("Repeated Over!!", i)
ks.append(a[0, 0])
return np.array(ks)
前回作成した関数 qr_eig_shiftd() とほぼ同じですが、最初に引数の行列 a を tridiag1() で三重対角化するところと、QR 分解に関数 qr_tri() を使うところが異なります。QR 分解は while ループの中で行っているので、三重対角行列専用の関数 qr_tri() を使うことで高速化が期待できます。これはあとで試してみましょう。
簡単な実行例を示します。ご参考までに、収束の様子を表示しています。
>>> c
array([[1, 4, 5],
[4, 2, 6],
[5, 6, 3]])
>>> qr_eig_tri(c)
0 3 -3.5207972894
[[ 1.20145137e+01 -1.53143502e+00 3.33066907e-16]
[ -1.53143502e+00 -2.35012171e+00 -7.08263423e-02]
[ 4.59423293e-16 -7.08263423e-02 -3.66439201e+00]]
1 3 -3.66819783368
[[ 1.21750977e+01 -1.13239864e-01 4.33680869e-16]
[ -1.13239864e-01 -2.50641459e+00 2.95524544e-05]
[ 4.55054030e-16 2.95524544e-05 -3.66868310e+00]]
2 3 -3.66868309795
[[ 1.21759664e+01 -8.30084433e-03 4.35743712e-16]
[ -8.30084433e-03 -2.50728327e+00 -3.04511733e-16]
[ 4.55042407e-16 -1.00508522e-17 -3.66868310e+00]]
3 2 -2.50728796709
[[ 1.21759711e+01 1.15532584e-15]
[ 8.40335200e-20 -2.50728797e+00]]
array([ -3.6686831 , -2.50728797, 12.17597107])
>>> d
array([[ 6., 1., 1., 1.],
[ 1., 7., 1., 1.],
[ 1., 1., 8., 1.],
[ 1., 1., 1., 9.]])
>>> qr_eig_tri(d)
0 4 6.42264973081
[[ 8.97741269e+00 -2.58338262e+00 1.00586083e-15 -1.72535908e-16]
[ -2.58338262e+00 7.13358012e+00 2.32518880e-01 4.02455846e-16]
[ 1.13360881e-16 2.32518880e-01 7.49567838e+00 -3.45458067e-02]
[ 2.80048133e-16 1.03479346e-16 -3.45458067e-02 6.39332881e+00]]
1 4 6.39224726776
[[ 1.06395144e+01 -9.34697934e-01 9.28521421e-16 3.79441330e-16]
[ -9.34697934e-01 5.47712638e+00 1.96185628e-01 3.81639165e-17]
[ 1.22244881e-16 1.96185628e-01 7.49108393e+00 8.68091796e-07]
[ -1.28793409e-16 2.69446074e-16 8.68091796e-07 6.39227529e+00]]
2 4 6.39227529027
[[ 1.07934382e+01 -2.39021506e-01 1.02695658e-15 -3.62370381e-16]
[ -2.39021506e-01 5.32412798e+00 1.96764180e-01 2.55713397e-16]
[ 1.17515441e-16 1.96764180e-01 7.49015852e+00 7.53213355e-16]
[ 1.83695271e-16 2.31774684e-16 -4.15551556e-17 6.39227529e+00]]
3 3 7.50788764581
[[ 1.07992387e+01 -1.59927486e-01 -1.00613962e-15]
[ -1.59927486e-01 5.30073731e+00 1.24529391e-05]
[ -1.14266539e-16 1.24529391e-05 7.50774871e+00]]
4 3 7.50774870536
[[ 1.08017929e+01 -1.07358757e-01 1.03190117e-15]
[ -1.07358757e-01 5.29818309e+00 5.07133872e-16]
[ 1.14131741e-16 -4.15956942e-18 7.50774871e+00]]
5 2 5.29608964531
[[ 1.08038864e+01 1.94289029e-16]
[ -4.65869187e-17 5.29608965e+00]]
array([ 6.39227529, 7.50774871, 5.29608965, 10.80388636])
>>> 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_tri(e)
0 5 9.17082039325
[[ 7.53181722e+00 -2.68627173e+00 -1.32530587e-15 -5.81904209e-16 -6.17194761e-16]
[ -2.68627173e+00 1.21373060e+01 -3.52959217e-01 6.71816773e-16 4.64431551e-16]
[ -2.33813510e-16 -3.52959217e-01 7.55932190e+00 4.90314959e-01 2.77555756e-16]
[ -7.58193089e-17 2.38790622e-16 4.90314959e-01 8.72658473e+00 -5.86229152e-01]
[ -1.55014516e-16 1.42362413e-16 3.71282797e-17 -5.86229152e-01 9.04497016e+00]]
1 5 9.49323685745
[[ 7.97707959e+00 -3.02484815e+00 -1.59772657e-16 1.41624091e-16 -1.17349348e-16]
[ -3.02484815e+00 1.16943961e+01 -1.89021997e-01 -1.39670781e-15 2.91314074e-16]
[ -1.25842996e-16 -1.89021997e-01 7.39545218e+00 2.09306199e-01 -5.55111512e-17]
[ 1.04363197e-16 -3.17661052e-16 2.09306199e-01 8.39430281e+00 4.27716821e-02]
[ -6.57059826e-17 1.75300016e-17 -1.22697385e-17 4.27716821e-02 9.53876930e+00]]
2 5 9.54036556488
[[ 8.44329835e+00 -3.27011113e+00 -1.05162728e-15 -1.29168070e-15 1.55788979e-16]
[ -3.27011113e+00 1.12269720e+01 -1.10748713e-01 8.17557264e-16 -2.09460822e-16]
[ -3.46918859e-17 -1.10748713e-01 7.36567317e+00 1.08808987e-01 -1.42247325e-16]
[ -2.38048992e-16 2.30112873e-16 1.08808987e-01 8.42366202e+00 -1.10180441e-06]
[ -2.36925232e-17 -5.82261204e-17 -8.85746838e-18 -1.10180441e-06 9.54039443e+00]]
3 5 9.54039442569
[[ 8.97314326e+00 -3.44946821e+00 9.28206042e-17 3.69872171e-16 -1.49031720e-16]
[ -3.44946821e+00 1.06961840e+01 -6.62804112e-02 -1.53653111e-15 2.09893879e-16]
[ -1.81335123e-17 -6.62804112e-02 7.35840009e+00 5.54613329e-02 -2.00646389e-16]
[ 1.41717382e-16 -3.00431031e-16 5.54613329e-02 8.43187827e+00 1.03443242e-15]
[ -6.27382944e-17 -4.20966945e-18 -8.72245433e-18 -4.43516041e-19 9.54039443e+00]]
4 4 8.43473607998
[[ 1.17064098e+01 -3.02355802e+00 -9.73617563e-16 1.62861246e-15]
[ -3.02355802e+00 7.96219116e+00 -2.33331306e-02 -1.28803218e-16]
[ 1.91991166e-17 -2.33331306e-02 7.35626796e+00 3.02353106e-08]
[ 3.18125643e-16 -9.27726297e-17 3.02353108e-08 8.43473667e+00]]
5 4 8.4347366665
[[ 1.29958222e+01 -1.62841264e+00 6.62664404e-16 1.28348271e-15]
[ -1.62841264e+00 6.67251013e+00 -1.04860762e-02 -1.00874597e-15]
[ -1.09275157e-17 -1.04860762e-02 7.35653660e+00 2.12260792e-16]
[ 2.96599114e-16 -1.47775458e-16 1.43716986e-18 8.43473667e+00]]
6 3 7.35669731302
[[ 1.33772036e+01 -3.07718479e-01 -8.18355800e-16]
[ -3.07718479e-01 6.29103346e+00 -6.18875671e-07]
[ 1.09793540e-17 -6.18875671e-07 7.35663185e+00]]
7 3 7.35663185484
[[ 1.33901140e+01 -5.51238493e-02 7.91901214e-16]
[ -5.51238493e-02 6.27812305e+00 -5.38469432e-16]
[ -1.09650410e-17 5.59313840e-19 7.35663185e+00]]
8 2 6.27769581992
[[ 1.33905412e+01 1.79023463e-15]
[ -1.10361221e-17 6.27769582e+00]]
array([ 9.54039443, 8.43473667, 7.35663185, 6.27769582, 13.39054123])
正常に動作していますね。最後に、qr_eig_shiftd() (ハウスホルダー変換版) と qr_eig_tri() の実行時間を比較してみましょう。100 * 100 の行列の固有値を求めます。結果は次のようになりました。
>>> s = time.time(); qr_eig_shiftd(np.ones((100,100)) + np.diag(np.arange(101, 201))); print(time.time() - s)
array([ 199.78276406, 198.74325265, 197.71691823, 196.6963277 ,
・・・省略・・・
103.19130324, 102.17693717, 101.15470811, 258.69669139])
1.3814737796783447
>>> s = time.time(); qr_eig_tri(np.ones((100,100)) + np.diag(np.arange(101, 201))); print(time.time() - s)
array([ 154.41654998, 153.41274961, 155.42038716, 152.40898353,
・・・省略・・・
199.78276406, 102.17693717, 101.15470811, 258.69669139])
0.565026044845581
qr_eig_tri() のほうが約 2.5 倍高速になりました。行列が大きくなると、その差はもっと大きくなるでしょう。興味のある方はいろいろ試してみてください。
●参考文献・URL
- 奥村晴彦,『Ｃ言語による最新アルゴリズム事典』, 技術評論社, 1991
- 行列の固有値問題 (PDF), (桂田祐史さん)
- 固有値問題ノートの補足, (桂田祐史さん)
- 固有値問題 (1) 対称行列の固有値, (fussy さん)
- HouesHolder変換 - [物理のかぎしっぽ], (東條遼平さん)
- 三重対角化 - [物理のかぎしっぽ], (東條遼平さん)
- QR法 - [物理のかぎしっぽ], (東條遼平さん)
- ハウスホルダー変換 - Wikipedia
●プログラムリスト
#
# qr1.py : ハウスホルダー変換による QR 分解と QR 法
#
# Copyright (C) 2018-2023 Makoto Hiroi
#
import numpy as np
#
# ハウスホルダー変換
#
def householder(xs):
n = len(xs)
x = np.linalg.norm(xs)
if xs[0] < 0: x = -x
vs = xs.astype(np.float_)
vs[0] += x
vs /= np.sqrt(vs[0] * x)
return np.eye(n) - vs.reshape((n, 1)) @ vs.reshape((1, n))
#
# QR 分解
#
def qr_h(a):
n = len(a)
xs = np.eye(n)
for i in range(n - 1):
h = np.eye(n)
h[i:, i:] = householder(a[i:, i])
a = h @ a
xs = xs @ h
return xs, a
#
# 三重対角化
#
def tridiag(a):
n = len(a)
for i in range(n - 2):
h = np.eye(n)
h[i+1:, i+1:] = householder(a[i+1:, i])
a = h @ a @ h
return a
def tridiag1(a):
n = len(a)
a = a.astype(np.float_)
for i in range(n - 2):
h = np.eye(n - i)
h[1:, 1:] = householder(a[i+1:, i])
a[i:, i:] = h @ a[i:, i:] @ h
return a
#
# 三重対角行列の QR 分解
#
def qr_tri(a):
n = len(a)
qs = np.eye(n) # 直交行列
for y in range(n - 1):
x = y + 1
d = np.sqrt(a[x, y] ** 2 + a[y, y] ** 2)
c = a[y, y] / d
s = a[x, y] / d
a1 = a[y].copy()
a2 = a[x]
a[y] = a1 * c + a2 * s
a[x] = a1 * (-s) + a2 * c
q1 = qs[:, y].copy()
q2 = qs[:, x]
qs[:, y] = q1 * c + q2 * s
qs[:, x] = q1 * (-s) + q2 * c
return qs, a
#
# 行列 [[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
#
# 三重対角化による QR 法
#
def qr_eig_tri(a, max_iter = 1024, e = 1e-14):
a = tridiag1(a)
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_tri(a - dk)
a = r @ q + dk
# print(i, size, k)
# print(a)
i += 1
if np.fabs(a[-1, -2]) < e:
ks.append(a[-1, -1])
break
else:
raise Exception("Repeated Over!!", i)
ks.append(a[0, 0])
return np.array(ks)