はじめに
多くの理・工系分野で線形代数の固有値問題に直面することは非常に多い。
ここではPythonのnumpyおよびscipyを用いて(一般化)固有値問題を解く方法をまとめる。
メソッドは補遺にまとめた。
一般化固有値問題とは，$A$と$B$を行列として，以下の方程式を満たすベクトル$\mathbf{x}$と(複素)数値$\lambda$を求める問いである。
$ A \ \mathbf{x} = \lambda \ B \ \mathbf{x} \ {\tag 1}$
$B$が単位行列$E$の場合，問題は
$ A \ \mathbf{x} = \lambda \ \mathbf{x} \ {\tag 2}$
となる。この問題を標準固有値問題という。
＃なお，固有値問題の数値計算手法としては，べき乗法やQR法などがあるが[1]，本稿ではそれらには触れずライブラリの基本的な利用法[2]だけ述べる。
追記: 2017 9/12
● べき乗法による数値解法に関する記事を投稿しました。
[Pythonによる科学・技術計算]べき乗法による行列の固有値問題の数値解法，線形代数
内容
まず，標準固有値問題である上記の式(2)を考える。
numpy.linalg.eigメソッドを用いる。
(1) 実行列の場合
(2) 複素行列の場合
次に，上記の式(1)型を考える。scipy.linalg.eigを用いる。
(3) 一般固有値問題の解法
なお，$B$を単位行列にすれば，(1)と(2)の問題も解ける。
☆ 実行列の計算コード(1)および(2)で固有ベクトルを正規化するため得られた固有ベクトルをそのノルムで割った。
つまり
for i in range(len(eig_vec)):
eig_vec [i] = eig_vec[i]/np.linalg.norm(eig_vec[i])
とした。
複素行列の場合(コード(2))はもとから正規化された解が出力される。
コード(1): 実行列の固有値・固有ベクトル
import numpy as np
"""
実行列の固有値・固有ベクトルの計算
"""
#A_matrix=np.matrix([[2,4],[0,4]]) # 行列Aを生成
A_matrix=np.array([[2,4],[0,4]]) # 行列Aを生成
eig_val, eig_vec =np.linalg.eig(A_matrix) # 固有値・固有ベクトルをそれぞれeig_val, eig_vecに格納する
for i in range(len(eig_vec)):
eig_vec [i] = eig_vec[i]/np.linalg.norm(eig_vec[i])
#計算結果の出力
print (eig_val)
print(eig_vec)
結果(1)
[ 2. 4.] #2つの固有値 (2と4)
[[ 0.74535599 0.66666667]
[ 0. 1. ]]
対応する固有ベクトル。1に規格化されている。
コード(2): 複素行列の固有値・固有ベクトル
実行列と全く同じように計算できる。複素数 $a+bi$ を，pythonでは iをjに変えてa+bjと表すことに注意する。
import numpy as np
"""
複素行列の固有値・固有ベクトルの計算
"""
#A_matrix=np.matrix([[1+4j,2j],[2,4j]]) # 複素数が行列要素にある
A_matrix=np.array([[1+4j,2j],[2,4j]]) # 複素数が行列要素にある
eig_val, eig_vec =np.linalg.eig(A_matrix) # 固有値・固有ベクトルの計算
print (eig_val)
print(eig_vec)
結果(2)
固有値
[ 1.95907589+5.37073062j -0.95907589+2.62926938j]
固有ベクトル (ノルムが1に規格化された複素ベクトル)
[[ 0.76703668+0.j -0.36782319-0.52570033j], [ 0.52570033-0.36782319j 0.76703668+0.j ]]
コード(3): 一般化固有値問題
import numpy as np
import scipy.linalg
A = np.array([[-32,16,0],[16,-32,16],[0,16,-32]])
B = np.array([[-2, 16, 13], [16, -5, 16], [0, 1, -32]]) # 行列B
#B=np.identity(3)
eig_val,eig_vec = scipy.linalg.eig(A,B) # 一般固化有値問題を解く
for i in range(len(eig_vec)): #正規化
eig_vec [i] = eig_vec[i]/np.linalg.norm(eig_vec[i])
print (eig_val)
print(eig_vec)
*結果(3):
固有値: 固有値は一般に複素数なので，純虚数としてjを用いて複素数表示されている
[-0.81267990+0.j 2.48526372+0.j 1.00000000+0.j]
array([-0.81267990+0.j, 2.48526372+0.j, 1.00000000+0.j])
固有ベクトル:
[[ -7.43016156e-01 -5.06563678e-01 4.37401683e-01]
[ -6.38461105e-01 7.69654089e-01 2.24401124e-16]
[ -2.11406258e-01 -2.50015624e-01 -9.44880724e-01]]
補遺: numpyおよびscipyのlinalgのメソッド
1. eig：一般行列の固有値・固有ベクトルの計算する
2. eigh：エルミート・実対称行列の固有値・固有ベクトルの計算
3. eigvals：一般行列の固有値の計算 (固有ベクトルは計算しない)
4. eigvalsh：エルミート・ 実対称行列の固有値の計算 (固有ベクトルは計算しない)
特に，対称行列に対するeighの計算は，一般行列に対して用いるeigよりも
非常に高速である。
参考文献
数値線形代数についてはいろいろなテキストが出ているようです。。
[1] F. シャトラン， 『行列の固有値』, シュプリンガーフェアラーク東京, 2003.,
森正武ほか， 『岩波講座 応用数学〈12 〔方法2〕 線形計算』, 岩波書店, 1997.
[2] scipyのwebサイト: https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.eig.html
https://docs.scipy.org/doc/scipy-0.18.1/reference/tutorial/linalg.html
なお，scipy.linalgは常にBLAS/LAPACKとリンクされているため(numpyはオプション)，numpy.linalgよりも基本的に高速である旨がかかれている。
以下，その部分を転載する。
"scipy.linalg vs numpy.linal
scipy.linalg contains all the functions in numpy.linalg. plus some other more advanced ones not contained in numpy.linalg
Another advantage of using scipy.linalg over numpy.linalg is that it is always compiled with BLAS/LAPACK support, while for numpy this is optional. Therefore, the scipy version might be faster depending on how numpy was installed.
Therefore, unless you don’t want to add scipy as a dependency to your numpy program, use scipy.linalg instead of numpy.linalg"