おっぱいそん！

pythonを数値計算に使うときのテクニックとかをまとめていきたい。

固有値問題

numpy scipy

Pythonで固有値問題を解く方法についてメモしておく。
メジャーな方法として、以下の3つがある

  1. numpy.linalgの関数を使う。
  2. scipy.linalgの関数を使う。
  3. scipy.sparse.linalgの関数を使う。


numpy.linalgとscipy.linalgには以下の4つの関数がある。

関数名のhはHermitianの略。Scipyだと一般化固有値問題*1もオプションで出来る。
ちなみに、scipy.linalgはnumpy.linalgに含まれる関数はすべて含んでいて、さらに追加で他の関数も含んでいる。また、scipy.linalgが常にBLAS/LAPACKでコンパイルされるのに対し、numpy.linalgはそうとは限らないので、scipyの方が早い場合もある。なので、特別な理由（コードが既にnumpy.linalgで書かれていて、書きなおすのがめんどいとか）がない限り、scipy.linalgを使うほうがよい。たぶん。
http://docs.scipy.org/doc/scipy-0.14.0/reference/tutorial/linalg.html


scipy.sparse.linalgには以下の2つの関数がある。

関数名のsはsparseの略。一般化固有値問題も出来る。

詳細

以下で、簡単に使い方を解説する。

w, v = numpy.linalg.eigh(A, UPLO='L')

  • エルミート（or 実対称）行列Aの固有値・固有ベクトルを求める。
  • UPLOは行列Aの下三角の部分を使うか(‘L’)、上三角の部分を使うか(‘U’)を指定する（デフォルトでは‘L’）。Aはエルミートなので、どちらを使っても同じ結果が得られるはず。
  • 返り値は固有値w(エルミートより常に実)と固有ベクトルからなる行列v(ユニタリ−行列)。ただし、固有値wの順番は特に決まっていない*2。vの列v[:, i]が固有値w[i]に対応する規格化された固有ベクトル。
  • 固有値だけでいいなら、 numpy.linalg.eigvalsh(A, UPLO=‘L’)を使えばいい。
  • ちなみに、中身はLAPACKの_ssyevd, _heevd*3が使われているらしい。

行列がエルミートでない一般の行列の時は

numpy.linalg.eig(A) or numpy.linalg.eigvals(A)

を使う。

http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eigh.html
http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eig.html


scipy.linalg.eig(a, left=False, right=True, overwrite_a=False, check_finite=True)

http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.eig.html

scipy.linalg.eigh(a, lower=True, eigvals_only=False, overwrite_a=False, eigvals=None, check_finite=True)

  • 固有値wはnumpy.linalg.eighと違って、昇順に並んでいる。
  • lowerはaの下三角からとるか、上三角から取るか。
  • eigvals_onlyは固有ベクトルを計算するかどうか。
  • eigvalsは求める固有値・固有ベクトルの範囲（指定がなければ全て求める）。

http://docs.scipy.org/doc/scipy-dev/reference/generated/scipy.linalg.eigh.html

注意

import scipy as spでインポートして、sp.linalg.eighで使おうとすると、

Traceback (most recent call last):
File "eigen.py", line 16, in
eval, evec = sp.linalg.eigh(A)
AttributeError: 'module' object has no attribute 'linalg'

のようなエラーがでる。
これは単にScipyをインポートしただけでは、Numpyと違ってサブパッケージまではインポートしないため（Scipyはかなり多くのサブパッケージを含んでいるので、このような仕様になっているみたい）。
たとえば、from scipy import linalgでインポートして、scipy.linalg.eighで使うようにすればいい。
http://stackoverflow.com/questions/9819733/scipy-special-import-issue

注意

左固有ベクトルの定義がややこしい（普通と違う）ので注意。
http://stackoverflow.com/questions/15560905/is-scipy-linalg-eig-giving-the-correct-left-eigenvectors


scipy.sparse.linalg.eigs(A, k=6, sigma=None, which='LM', v0=None, ncv=None, maxiter=None, tol=0, return_eigenvectors=True)

  • scipy.sparseの場合はAにndarrayだけでなく、疎行列やLinearOperatorも使える。
  • k: 求める固有値・固有ベクトルの数。Aの次元Nより小さくないとだめ（すべての固有値は求められない）。
  • sigmaはshift-invert mode。この値のあたりの固有値を求める時に、指定すると有効。
  • v0はiterationのInitial vectorの指定。
  • ncv:Lanczosベクトルの数。kより大きくないとダメ。
  • which:どの固有値を求めるか。（最大固有値を求める時'LM'はなくてもいいが、最小固有値'SM'を求める時はsigmaを指定したほうがよい。）
  • maxiter: Arnoldiのiterationの回数の最大値。
  • tol: 固有値の相対精度（iterationをストップする値）。
  • return_eigenvectors: 固有値を求めるかどうか。
  • 返り値はk個の固有値wとk個の固有ベクトルv[:,i]
  • ARPACKのSNEUPD, DNEUPD, CNEUPD, ZNEUPDという関数のラッパー。中身はthe Implicitly Restarted Arnoldi Methodを使っている。

http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigs.html

scipy.sparse.linalg.eigsh(A, k=6, M=None, sigma=None, which='LM', v0=None, ncv=None, maxiter=None, tol=0, return_eigenvectors=True)

  • ARPACKのSSEUPDとDSEUPDという関数のラッパー。中身はRestarted Lanczos Methodが動いている。
注意：最小固有値を求める場合

ARPACKは最大固有値を求める時に最適化されているようなので、最小固有値の場合は工夫した方がいい（そのままだと収束に失敗したりする）。

  • tolerance (tol) を増やして、はやく収束させる。→ 収束には成功するが、精度が悪くなる。
  • maximum number of iterations (maxiter)を増やす。→ 精度はいいが、時間がかかる。
  • shift-invert modeを(sigma)使う。→ 短い時間で精度も良い。

らしいので、sigmaを指定するのがおすすめ。
http://students.mimuw.edu.pl/~pbechler/scipy_doc/tutorial/arpack.html


使い方

エルミート行列の固有値・右固有ベクトルをscipy.linalgとscipy.sparse.linalgで。

import numpy as np
import scipy.linalg
import scipy.sparse.linalg

N = 3; np.random.seed(5)
A = np.random.rand(N,N); A = np.dot(A, A.conj().T) #make Hermitian matrix

w1, v1 = scipy.linalg.eigh(A)
print 'eigenvalues =', w1
print 'A v =', np.dot(A, v1); print 'v w =', np.dot(v1, np.diag(w1))
w2, v2 = scipy.sparse.linalg.eigsh(A, k=2)
print 'eigenvalues =', w2
i = 0
print 'A v[:,i] =', np.dot(A, v2[:, i]); print 'w[i] v[:,i] =', w2[i] * v2[:, i]

結果は

eigenvalues = [ 0.02340572 0.34973568 2.87713176]
A v = [[-0.00366442 0.30775003 1.29047002]
[-0.0131621 -0.15834243 1.99081348]
[ 0.01900418 -0.05032536 1.62764737]]
v w = [[-0.00366442 0.30775003 1.29047002]
[-0.0131621 -0.15834243 1.99081348]
[ 0.01900418 -0.05032536 1.62764737]]
eigenvalues = [ 0.34973568 2.87713176]
A v[:,i] = [-0.30775003 0.15834243 0.05032536]
w[i] v[:,i] = [-0.30775003 0.15834243 0.05032536]

scipy.linalgでは、固有値が昇順に並んでいる。
scipy.sparse.linalgでは、全ての固有値を求めることは出来ず、求めたい固有値の数kを設定する（kは行列Aの次元より小さくないといけない）。
scipy.sparse.linalgも固有値は昇順に並ぶ。
行列Aの固有値wと右固有ベクトルからなる行列V（固有値w[i]に対応する右固有ベクトルがv[:,i]）は
A V = V W ⇔ A v[:,i] = w[i] v[:,i]
を満たす。

左固有ベクトル

エルミートの時（eig(s)h）は、右固有ベクトルのエルミート共役をとれば、左固有ベクトルが作れるので、自明。
そこで、一般の行列の右固有ベクトルを考える。
scipy.linalg.eigは右固有ベクトルを求めるかどうかのオプション（left=True）があるので、それを使えばいい（ただし、右固有ベクトルの定義が普通と違うので、注意する）。
scipy.sparse.linalg.eigsにはそのようなオプションはないので、エルミート共役を取った行列に対して、右固有ベクトルを求めることで、元の行列の左固有ベクトルを求める。:

import numpy as np
import scipy.linalg
import scipy.sparse.linalg

N = 3; np.random.seed(5)
B = np.random.rand(N,N)

w1, v1_left = scipy.linalg.eig(B, left=True, right=False)
print np.dot(B.conj().T, v1_left) - np.dot(v1_left, np.diag(w1.conj()))
print np.dot(v1_left.conj().T, B) - np.dot(np.diag(w1), v1_left.conj().T)

w2, v2_left = scipy.sparse.linalg.eigs(B.conj().T, k=1)
w2= w2.conj(); v2_left = v2_left.conj().T
i = 0
print np.dot(v2_left[i,:], B) - w2[i] * v2_left[i,:]

結果

[[ 4.44089210e-16+0.j 2.77555756e-17+0.j -6.59194921e-17+0.j]
[ 8.88178420e-16+0.j 1.11022302e-16+0.j -1.94289029e-16+0.j]
[ 3.33066907e-16+0.j -5.55111512e-17+0.j -2.77555756e-17+0.j]]
[[ 4.44089210e-16+0.j 8.88178420e-16+0.j 2.22044605e-16+0.j]
[ 2.77555756e-17+0.j 1.11022302e-16+0.j -8.32667268e-17+0.j]
[ -6.59194921e-17+0.j -1.94289029e-16+0.j -2.77555756e-17+0.j]]
[ 2.22044605e-16+0.j 2.22044605e-16+0.j 0.00000000e+00+0.j]

マニュアルにあるように、scipy.linalg.eig(B, left=True)で求められる左固有ベクトルv_leftは
B.H v_left[:,i] = w[i].conj() v_left[:,i] ⇔ B.H V_left = V_left W
⇔ V_left.H B = W V_left.H
を満たすもの(.Hはエルミート共役の意味)。
普通、右固有ベクトルは
V_left B = W V_left ⇔ v_left[i,:] B = w[i] v_left[i,:]
⇔ B.H V_left.H = V_left.H W.conj()
で定義されるので、scipy.linalg.eigで求められる左固有ベクトルは普通の定義のエルミート共役になっている。

scipy.sparse.linalg.eigsでは、行列のエルミート共役を取ってものの固有値・右固有ベクトルを求め、それぞれ複素共役・エルミート共役をとれば、求めたい固有値・左固有ベクトルが得られる。


比較

3つの方法の違いをまとめておく。

  1. numpy.linalg
  2. scipy.linalg
  3. scipy.sparse.linalg
  • 1は一般化固有値問題は出来ない。
  • 1は右固有ベクトルは求められない。
  • 1,2はLAPACK(Linear Algebra PACKage)、3はARPACK(ARnoldi PACKage)が中で動いている。
  • 3はnumpy.ndarrayだけでなく、scipy.sparse.csr_matrixやscipy.sparse.linalg.LinearOperatorにも使える。

linalg.eig(s)では

  • 1,2は常に全ての固有値を求めるのに対し、3は求める固有値の数とどういう固有値のみを求めるか（whichを使う）を指定できる（ただし、全ての固有値は求められない）。また、1,2では固有値の並びは決まっていないのに対し、3では固有値の並びは昇順。

[2015.7.1 訂正] scipy.sparseでもeigsの固有値の並びは決まっていませんでした。実際、適当にランダムな行列を作って計算を回してみると、昇順になる時や、降順になる時や、バラバラになる時があるようなので、固有値の順を指定したいときは、ソートしないといけない。
エルミート行列のlinalg.eig(s)hでは

  • 1は常に全ての固有値を求めるのに対し、2,3は求める固有値の数を指定できる。さらに、3はどういう固有値のみを求めるか（whichを使う）を指定できる（ただし、3は全ての固有値は求められない）。固有値の並びは、1は決まっておらず、2,3は昇順。

※ 3.spaeseが全ての固有値を求められないのはLanczos or Arnoldiに基づいているので。

以上より、2.scipy.linalgは1.numpy.linalgの上位互換なので、1を使う理由はない(早さは1と2はほとんど同じ)。
2.scipy.linalgと3.scipy.sparse.linalgのどちらを選ぶべきかは、固有値を求めたい行列の性質や何を求めたいかに依って変わる（早さは行列の性質に依る）。
例えば、密な行列の全ての固有値を求めたい時は2を、疎な行列の最大固有値のみを求めたい時は3を選ぶ（ただし、最大固有値がほぼ縮退している時は、3はiterative methodsなので遅い）。
密行列でデモしてみたところ、最大固有値のみを求める場合、縮退がある時は2の方が2倍くらい速く、縮退がない時は3の方が100倍くらい早かった。
http://stackoverflow.com/questions/11083660/python-eigenvectors-differences-among-numpy-linalg-scipy-linalg-and-scipy-spar
http://stackoverflow.com/questions/24418534/least-eigenvalue-in-python
また、上の3つ以外にもMPIを使ったSLEPcとかもあるらしいので、並列化する人はありかも。
http://stackoverflow.com/questions/6684238/whats-the-fastest-way-to-find-eigenvalues-vectors-in-python

*1:普通の固有値問題がAx=axに対し、一般化固有値問題はAx=aMxのような固有値問題のこと。

*2:numpy.linalg.eighは固有値の順を揃えないので、numpy.linalg.eighでソートしたいなら、工夫しないとだめ。例えば、以下をリンクを参考。http://stackoverflow.com/questions/7839413/python-numpy-compute-first-eigenvalue-and-eigenvector、http://d.hatena.ne.jp/nohzen/20140407/1396852091 ただし、scipyを使えば、固有値の順を最初から揃えてくれるので、Scipyを使った方が楽です。

*3:それぞれ、実対称行列、エルミート行列のすべての固有値を求め、固有ベクトルを分割統治法を使って求める関数 http://d.hatena.ne.jp/blono/20080925/1222296801