実対称行列の固有値・固有ベクトルをNumpyで計算した話

つまりなにしたの？

実対称行列の固有値・固有ベクトルをNumpyで計算した。
1つは numpy.linalg.eigによる方法
2つは Jacobi法を自前で実装する方法(Numpyのアシスト付き）
1を信用して、2があってるか確認した。
f:id:ensekitt:20180820012846j:plain
*1

固有値・固有ベクトルってなに？

正方行列Aに対して
$Ax = \lambda x (x \neq 0)$
が成り立つときのλが固有値でxが固有ベクトル。

行列Aに関して、ベクトルxに行列Aをかけてもベクトルxのλ倍のベクトルになるってことなので、
ベクトルの向きは変わらず大きさだけ変わる。
逆に行列Aを何も考えてないベクトル（Aの固有ベクトルではない大きさを持つベクトル）に掛けると普通は向きも変わる。

そんな面白い性質を持つxとλを求めてみようと思う。

どうやって求めるのか？

QR法とかべき乗法とかJacobi法とかが用いられるみたい。
ヤコビ法はAが実対称行列でそんなに大きくないときに使うこともあるくらいのやつで、実装が簡易なので今日はそれを試してみる。
numpy.linalg.eigは様々な行列に対応できるように、LAPACKのxGEEV関数が使われるらしい。
xGEEV関数の中身はよくわかっていないけど、よく最適化された固有値求める関数。

実行環境

>>> 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でやってみる

import numpy as np
# 例にした行列
val = np.array([
    [1,2,3],
    [2,4,5],
    [3,5,1]])
print("val\n", val)
eig_val, eig_vec = np.linalg.eig(val)
    
print("eigen_values", eig_val)
print("eigen_vectors\n", eig_vec.T)

val
 [[1 2 3]
 [2 4 5]
 [3 5 1]]
eigen_values [ 9.15554844  0.03424156 -3.18979   ]
eigen_vectors
 [[-0.38907681 -0.71466558 -0.58126788]
 [-0.84846243  0.52376469 -0.07603987]
 [-0.35879067 -0.4635986   0.81015159]]

Jacobi法でやってみる

Jacobi方の大まかな手順

初期値Rを設定
非対角成分の最大値を求める
収束判定する
最大値の行と列を取り出す
回転角度を計算する
行列を更新する
ギブンス回転行列を作って直交行列Rを更新する

コード

import numpy as np

val = np.array([
    [1,2,3],
    [2,4,5],
    [3,5,1]], dtype=np.float)

# 固定値
inveye = -(np.eye(val.shape[0]) - 1) # 非対角成分作るための行列
MAX = 100 # 最大繰り返し回数
TOL = 1e-10 # 収束条件

# 初期値設定
R = np.eye(val.shape[0])

for i in range(MAX):
    # 非対角成分の最大値を求める
    absval = np.abs(val * inveye).reshape(val.shape[0]*val.shape[1])
    amax = np.max(absval)
    aindex = np.argmax(absval)
    
    # 収束判定
    if amax < TOL:
        break

    # 最大値の行と列
    r = int(aindex / val.shape[0])
    c = int(aindex % val.shape[0])
    
    # 回転角度計算
    if val[r,r]-val[c,c] == 0:
        theta = 0.25*np.pi
    else:
        theta = 0.5*np.arctan(-2.0*val[r,c]/(val[r,r]-val[c,c]))
    
    # 行列を更新する
    new_val = val.copy()
    for k in range(val.shape[0]):
        new_val[r,k] = val[r,k]*np.cos(theta) - val[c,k]*np.sin(theta)
        new_val[k,r] = new_val[r,k]
        
        new_val[c,k] = val[r,k]*np.sin(theta) + val[c,k]*np.cos(theta)
        new_val[k,c] = new_val[c,k]

    new_val[r,r] = (val[r,r] + val[c,c] )/2 + ((val[r,r]-val[c,c])/2)*np.cos(2*theta)-val[r,c]*np.sin(2*theta)
    new_val[c,c] = (val[r,r] + val[c,c] )/2 - ((val[r,r]-val[c,c])/2)*np.cos(2*theta)+val[r,c]*np.sin(2*theta)
    new_val[r,c] = 0.0
    new_val[c,r] = 0.0
    val = new_val.copy()

    # ギブンス回転行列を作って直交行列Rを更新する
    G = np.identity(val.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)
    R = np.dot(R,G)

eig_val = np.diag(val)
eig_vec = R.copy()

print("eigen_values", eig_val)
print("eigen_vectors\n", eig_vec.T)

実行結果

eigen_values [ 0.03424156  9.15554844 -3.18979   ]
eigen_vectors
 [[ 0.84846243 -0.52376469  0.07603987]
 [ 0.38907681  0.71466558  0.58126788]
 [-0.35879067 -0.4635986   0.81015159]]