行列の固有値の計算

2020年3月6日

はじめに

Python で行列の固有値を計算する。

環境

  • Anaconda3 (Python 3.7)

行列の固有値の計算

Python で行列の固有値の計算をするための関数として、以下のようなものがある。

  • numpy.linalg.eig
  • scipy.linalg.eig
  • scipy.sparse.linalg.eigs

たとえば、以下のようにする。

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

行列によってはデタラメな結果を出すときがあるので注意。結果をチェックすること。