EnsekiTT Blog

EnsekiTTが書くブログです。

QR法をNumpyで実装して固有値を計算してみた話

つまりなにしたの?

固有値にはQR法を使うのが良いらしいのでNumpyで計算した。

  1. numpy.linalg.eigを使う方法(固有値の確認)
  2. numpy.linalg.qrを使う方法(QR分解はNumpyに頼る)
  3. QR分解を自前で実装する方法(Numpyのアシスト付き)

1と2を信用して、3の妥当性を確認した。
f:id:ensekitt:20180822032058j:plain
*1

前回に引き続いて

ensekitt.hatenablog.com
前回はJacobi法で実対称行列に対して固有値を求めた。

今回はQR法で実正方行列の固有値が求められるので、それをやってみる。
固有ベクトルは必要があれば固有値から求められるので、今回は対象としない。

QR法とは

QR法はQR分解を用いた固有値などの計算方法
QR分解はある行列(A)を直交行列Qと上三角行列Rの積に分解すること。
QR分解にも色々な方法があるらしいが、今回はギブンス回転を用いた方法を試すことにした。
ギブンス回転は直交行列で、主に相似変換によって行列に0の要素を増やすために使われる。
これをつかって行列を上三角行列にしているみたい。

Qを求めることができたら、Aを更新していくことで固有値を得ることができる。
 A_k = Q_k R
 A_{k+1} = R_k Q_k = Q_{k-1}^{-1}A_k Q_k = Q_k^T A_k Q_k
 (k = 1,2,...,n)

実行環境

>>> 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]

ちゃんと一致を確認できた。

クリエイティブ・コモンズ・ライセンス
この 作品 は クリエイティブ・コモンズ 表示 4.0 国際 ライセンスの下に提供されています。