pythonによる固有値固有ベクトル計算
pythonで固有値固有ベクトルを計算しようと思うと numpy.linalg (線形代数の便利関数のモジュール)以下にある numpy.linalg.eigs を使うことになります。
numpy.linalg.eig — NumPy v1.11 Manual
この挙動が個人的な感覚とちょっとずれていて、そこでかなり消耗してしまったのでメモとして残しておきます。
nump.linalg.eigs でハマる
numpy にはnumpy.linalg.eigs
という関数が用意されています。
返り値として第一変数に固有値(eigenvalues)、第二変数に正規化された(すなわち \( \|u_i\|=1 \) )固有ベクトル(eigenvector)を返します。
例
X = # 適当な(n,n)のnp.array
eig_vals,eig_vec = np.linalg.eigs(X)
ここで1点注意点があります。
http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eig.html
Returns:
w : (..., M) array The eigenvalues, each repeated according to its multiplicity. (以下略
v : (..., M, M) array The normalized (unit “length”) eigenvectors, such that the column v[:,i] is the eigenvector corresponding to the eigenvalue w[i].
重要なのは、index を i とした時に固有値 eig_vals[i]
に対応する固有ベクトルは eig_vec[i]
ではないということ!!
正しくはeig_vec[:,i]
が対応します。
この問題のやっかいな点がデバックすることが困難である点です。
なぜなら、インデックスを入れる次元を間違えて eig_val[i] としても、その配列の長さは正しい場合と同じ M になり、ベクトルの長さでは判定できません。また「固有ベクトル同士は直行する」という条件を使ってデバッグをしようとしても、間違えたベクトルに対しても成り立ってしまいます。すなわちi=jの時だけ
eig_vec[i].dot(eif_vec[j])=1
が成り立ちます。ですから間違えたindexを与えた際にも固有ベクトルの特徴である「お互いが直交する」という性質も満たしてしまっています。
したがって、計算した結果がなんだかおかしい、ということ以外でのデバッグが不可能になります。僕はこれで3時間ぐらい持って行かれました。おとなしく最初にドキュメントを読めばいいのですが。ちょっと慣れてくると適当にやってもそれっぽく出来てしまうのはダメだなと改めて感じました。ドキュメント大事。
までも普通に考えて、i番目の固有ベクトルだからeig_vec[i]
でいいと思うしソッチのほうが自然な定義な気がしてならないんですが、なぜそうなっていないのか……