行列の固有値の計算2020年3月6日 | |
はじめにPython で行列の固有値を計算する。 環境
行列の固有値の計算Python で行列の固有値の計算をするための関数として、以下のようなものがある。
たとえば、以下のようにする。 import numpy as np from scipy.linalg import eig from scipy.sparse.linalg import eigs 行列を用意。
# 正定値行列
A = np.array([
[1., 0., 0.,],
[0., 2., 1.,],
[0., 1., 1.],
])
固有値の計算。 value, vector = np.linalg.eig(A) value, vector 結果
(array([0.38196601, 2.61803399, 1. ]),
array([[ 0. , 0. , 1. ],
[ 0.52573111, -0.85065081, 0. ],
[-0.85065081, -0.52573111, 0. ]]))
あるいは value, vector = eig(A) value, vector 結果
(array([0.38196601+0.j, 2.61803399+0.j, 1. +0.j]),
array([[ 0. , 0. , 1. ],
[ 0.52573111, -0.85065081, 0. ],
[-0.85065081, -0.52573111, 0. ]]))
最大のものが 1 つだけ欲しい場合は、eigs を使う。 value, vector = eigs(A, 1) value, vector 結果
(array([2.61803399+0.j]), array([[-7.72116492e-17+0.j],
[ 8.50650808e-01+0.j],
[ 5.25731112e-01+0.j]]))
最小の場合。 value, vector = eigs(A, 1, which="SM") value, vector 結果
(array([0.38196601+0.j]), array([[-4.26671901e-16+0.j],
[-5.25731112e-01+0.j],
[ 8.50650808e-01+0.j]]))
一般化固有値問題 (行列が 2 つある) を解きたい場合は、次のようにすればよい。 value, vector = eig(A, B) value, vector = eigs(A, 1, B) 固有値を並び替えるたとえば、正の実数の固有値にしか興味がなく、それを小さい順に並び替えたい場合は、次のような関数を作ればよい。
def get_sorted_positive_value(value, vector):
index = np.where((value > 0.) & (np.imag(abs(value)) < 1e-5))[0]
value = np.real(value[index])
vector = np.real(vector[:, index])
index = np.argsort(value)
value = value[index]
vector = vector[:, index]
return value, vector
次のように使う。 value, vector = eig(A) value, vector = get_sorted_positive_value(value, vector) value, vector 結果
(array([0.38196601, 1. , 2.61803399]),
array([[ 0. , 1. , 0. ],
[ 0.52573111, 0. , -0.85065081],
[-0.85065081, 0. , -0.52573111]]))
固有ベクトルの正規化固有ベクトルをノルムで正規化するには次のような関数を作ればよい。
def normalize_eigenvector(vector):
vector = vector/np.linalg.norm(vector)
return vector
def normalize_eigenvectors(vector):
n = vector.shape[1]
for j in range(n):
vector[:, j] = normalize_eigenvector(vector[:, j])
return vector
実行 normalize_eigenvectors(vector) 結果
array([[ 0. , 1. , 0. ],
[ 0.52573111, 0. , -0.85065081],
[-0.85065081, 0. , -0.52573111]])
最大値で正規化したい場合は、次のようにする。
def max_normalize_eigenvector(vector):
max_index = np.where(abs(vector) == max(abs(vector)))[0][0]
vector = vector/vector[max_index]
return vector
def max_normalize_eigenvectors(vector):
n = vector.shape[1]
for j in range(n):
vector[:, j] = max_normalize_eigenvector(vector[:, j])
return vector
実行 max_normalize_eigenvectors(vector) 結果
array([[-0. , 1. , -0. ],
[-0.61803399, 0. , 1. ],
[ 1. , 0. , 0.61803399]])
固有値の検証固有値の検証は、次のような関数でできる。
def eig_check(A, value, vector):
return np.linalg.norm(A.dot(vector) - value*vector)
def eig_checks(A, value, vector):
r = np.zeros(value.shape[0])
i = 0
for l, v in zip(value, vector.T):
r[i] = eig_check(A, l, v)
i += 1
return r
次のように使う。 value, vector = eig(A) value, vector = get_sorted_positive_value(value, vector) eig_checks(A, value, vector) 結果 array([1.57009246e-16, 0.00000000e+00, 2.22044605e-16]) 本来は 0 にならないといけないので、0 に近ければよい。 減次行列から 1 つの固有値の成分を差し引くことができる。 # 減次 l = np.real(value[0]) v = np.array([vector[:, 0]]).T A2 = A - l*v.dot(v.T) value2, vector2 = eig(A2) value2, vector2 = get_sorted_positive_value(value2, vector2) value2, vector2 結果
(array([1. , 2.61803399]),
array([[ 1. , 0. ],
[ 0. , -0.85065081],
[ 0. , -0.52573111]]))
べき乗法あえて自分で固有値を計算してみるなら、べき乗法というものが使える。
def eig_power(A, eps=1e-5, maxiter=100):
x = np.ones((A.shape[0], 1))
x = x/np.linalg.norm(x)
v0 = np.zeros((A.shape[0], 1))
v = x
i = 0
while np.linalg.norm(v - v0) > eps:
v0 = v
v = A.dot(x)
l = v.T.dot(x)[0, 0]
x = v/np.linalg.norm(v)
i += 1
if i == maxiter:
break
x[np.abs(x) < eps] = 0.
return l, x
実行 value, vector = eig_power(A) value, vector 結果
(2.618033988738303,
array([[0. ],
[0.85065081],
[0.52573111]]))
最大の固有値が得られる。 複数の固有値が欲しい場合は、いい加減な方法としては、単純な減次を使う方法がある。
def eig_power_n(A, n=0, eps=1e-15, maxiter=100):
if n == 0 or n > A.shape[0]:
n = A.shape[0]
value = np.zeros(n)
vector = np.zeros((A.shape[0], n))
A2 = A
for i in range(n):
l, v = eig_power(A2, eps, maxiter)
value[i] = l
vector[:, i] = v[:, 0]
A2 = A2 - l*v.dot(v.T)
return value, vector
実行 value, vector = eig_power_n(A, 3) value, vector 結果
(array([2.61803399, 1. , 0.38196601]),
array([[ 0. , 1. , 0. ],
[ 0.85065081, 0. , -0.52573111],
[ 0.52573111, 0. , 0.85065081]]))
行列によってはデタラメな結果を出すときがあるので注意。結果をチェックすること。 逆べき乗法べき乗法は固有値の最大値を計算するが、最小値が欲しい場合もある。その場合は逆べき乗法を使う。
def eig_inverse_power(A, eps=1e-15, maxiter=100):
x = np.ones((A.shape[0], 1))
x = x/np.linalg.norm(x)
v0 = np.zeros((A.shape[0], 1))
v = x
i = 0
while np.linalg.norm(v - v0) > eps:
v0 = v
v = np.linalg.solve(A, x)
l = v.T.dot(x)[0, 0]
x = v/np.linalg.norm(v)
i += 1
if i == maxiter:
break
x[np.abs(x) < eps] = 0.
return 1./l, x
実行 value, vector = eig_inverse_power(A) value, vector 結果
(0.38196601125010515,
array([[ 0. ],
[-0.52573111],
[ 0.85065081]]))
最小の固有値が得られる。 複数の固有値が欲しい場合は、次のようにする。
def eig_inverse_power_n(A, n=0, eps=1e-15, maxiter=100):
if n == 0 or n > A.shape[0]:
n = A.shape[0]
value = np.zeros(n)
vector = np.zeros((A.shape[0], n))
A2 = A
for i in range(n-1):
l, v = eig_inverse_power(A2, eps, maxiter)
value[i] = l
vector[:, i] = v[:, 0]
A2 = A2 - (1./l)*v.dot(v.T)
l, v = eig_power(A2, eps, maxiter)
value[i+1] = l
vector[:, i+1] = v[:, 0]
return value, vector
計算 value, vector = eig_inverse_power_n(A, 3) value, vector 結果
(array([0.38196601, 1. , 2.61803399]),
array([[ 0. , 1. , 0. ],
[-0.52573111, 0. , 0.85065079],
[ 0.85065081, 0. , 0.52573114]]))
行列によってはデタラメな結果を出すときがあるので注意。結果をチェックすること。 | |
| PENGUINITIS | |