前回に引き続いて
ensekitt.hatenablog.com
前回はJacobi法で実対称行列に対して固有値を求めた。
今回はQR法で実正方行列の固有値が求められるので、それをやってみる。
固有ベクトルは必要があれば固有値から求められるので、今回は対象としない。
QR法とは
QR法はQR分解を用いた固有値などの計算方法
QR分解はある行列(A)を直交行列Qと上三角行列Rの積に分解すること。
QR分解にも色々な方法があるらしいが、今回はギブンス回転を用いた方法を試すことにした。
ギブンス回転は直交行列で、主に相似変換によって行列に0の要素を増やすために使われる。
これをつかって行列を上三角行列にしているみたい。
Qを求めることができたら、Aを更新していくことで固有値を得ることができる。
実行環境
>>> import sys >>> sys.version 3.6.4 (default, Mar 2 2018, 00:45:54) [GCC 4.2.1 Compatible Apple LLVM 9.0.0 (clang-900.0.39.2)] >>> import numpy >>> numpy.__version__ 1.14.1
numpy.linalg.eigを使う方法
import numpy as np A = np.array([[6, 4, 1], [1, 8, -2], [3, 2, 0]]) val, vec = np.linalg.eig(A) print(val) # [-0.42968559 6.29086844 8.13881715]
numpy.linalg.qrを使う方法
import numpy as np A = np.array([[6, 4, 1], [1, 8, -2], [3, 2, 0]]) MAX = 100 # 最大繰り返し回数 TOL = 1e-10 # 収束条件 for i in range(MAX): absval = np.abs(A*np.tri(A.shape[0], k=-1)).reshape(A.shape[0]*A.shape[1]) amax = np.max(absval) if amax < TOL: break Q, R = np.linalg.qr(A) A = np.dot(np.dot(Q.T, A), Q) print(np.diag(A)) # [ 8.13881715 6.29086844 -0.42968559]
QR分解から自前で実装する方法
import numpy as np A = np.array([[6, 4, 1], [1, 8, -2], [3, 2, 0]], dtype=np.float) MAX = 100 # 最大繰り返し回数 TOL = 1e-10 # 収束条件 for i in range(MAX): # 下三角要素の絶対値のうち最大値とその行、列を求める absval = np.abs(A*np.tri(A.shape[0], k=-1)).reshape(A.shape[0]*A.shape[1]) amax = np.max(absval) # 最大値が収束条件未満なら終了 if amax < TOL: break # Rを順次求めていく Q = np.identity(A.shape[0]) R = A.copy() for i in range(R.shape[0]*R.shape[1]): # Rの下三角要素の絶対値のうち最大値とその行、列を求める absval = np.abs(R*np.tri(R.shape[0], k=-1)).reshape(R.shape[0]*R.shape[1]) amax = np.max(absval) # 最大値が収束条件未満なら終了 if amax < TOL: break aindex = np.argmax(absval) # 最大値の行と列 r = int(aindex / R.shape[0]) c = int(aindex % R.shape[0]) # ギブンス回転行列の角度を計算 theta = np.arctan2(R[r,c],R[c,c]) # ギブンス回転行列を作って直交行列Pを更新する G = np.identity(R.shape[0]) G[r,r] = np.cos(theta) G[r,c] = -np.sin(theta) G[c,r] = np.sin(theta) G[c,c] = np.cos(theta) # QとRを更新する R = np.dot(G,R) Q = np.dot(Q,G.T) # Aを更新する A = np.dot(np.dot(Q.T, A), Q) val = np.diag(A) print(val) # [ 8.13881715 6.29086844 -0.42968559]
ちゃんと一致を確認できた。