pythonで特異値分解 - おっぱいそん！

pythonで特異値分解(singular value decomposition,SVD)をする時のメモ。

一般の密行列のSVD

あまり選択肢がないみたい。とりあえず、Numpy or ScipyのSVDを使っとけば間違いなさそう。

numpy.linalg.svd(a, full_matrices=1, compute_uv=1)
scipy.linalg.svd

http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.svd.html
http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.svd.html
行列AをSVDしたければ、

import numpy as np
X, Y, Z = np.linalg.svd(A)

import scipy.linalg
X, Y, Z = scipy.linalg.svd(A)

のように使います。
full_matrices=0を指定すれば、上の行列X,Zの次元が小さくなります。
compute_uv=1を指定すれば、X,Yを求めなくなります。
どちらも、LAPACKを使っているみたいです。
この2つに違いがあるのかは知りませんが、僕がテストしてみた範囲では速さや結果は変わらないみたいです。

ちなみに、疎行列用の関数scipy.sparse.linalg.svdsなどでも、密行列のSVDを出来るみたいです。

scipy.sparse.linalg.svds(A, k=Aの次元)

もちろん、疎行列のSVDをするなら、scipy.linalg.svdでなく、scipy.sparse.linalg.svdsや下の記事の関数を使う方が早いです。
一方、密行列のSVDはscipy.sparse.linalg.svdsより、scipy.linalg.svdの方が速いようなのですが、メモリはscipy.sparse.linalg.svdsの方が節約出来るみたいです(下のリンクを参照)。
http://fa.bianp.net/blog/2012/singular-value-decomposition-in-scipy/

疎行列のSVD

いろいろあるみたいですが、僕は使ってないので、下の記事を参考に。
幾つかの関数の速さを比較しています。
http://jakevdp.github.io/blog/2012/12/19/sparse-svds-in-python/
http://fa.bianp.net/blog/2012/singular-value-decomposition-in-scipy/

収束しないエラー

numpy.linalg.svdやscipy.linalg.svdを使っていると、

raise LinAlgError, ’SVD did not converge’
numpy.linalg.linalg.LinAlgError: SVD did not converge SVD

とエラーが出る時があります。
SVDは原理的には任意の行列に対して(エルミート行列の対角化を使って)実行できますが、numpy.linalg.svdでは高速化のために、"ある操作"を繰り返してSVDして求められる行列(に近い行列?)に収束させているようです。
その際、収束に必要な"ある操作"の繰り返しのステップ数が一定数n以上になると、諦めて上のようなエラーを吐くようです。
対策としては

numpy.linalg.svdモジュールの中の"ある操作"の繰り返しの上限nを書き換えて、収束が遅くてもエラーが出にくくする。
SVDする行列を2乗してからSVDする。
対称（orエルミート）行列の対角化を使う。

2つ目の方法のように、2乗してからSVDすると、なぜエラーがでなくなるかは知りません。
（そもそも、モジュールの中でどういうアルゴリズムでSVDしているのか知らないので…とりあえず、使う分にはエラーさえ出なければ困らない）
1つ目の方法は、モジュールの中身を書き換えられない環境の時には使えない。
2つ目の方法は以下のようにする。

import numpy as np
try:
    X, Y, Z = np.linalg.svd(A)
except:
    A2 = np.dot(A.T, A)
    X2, Y2, Z = np.linalg.svd(A2)
    Y = np.sqrt(Y2)
    X = np.dot(A, Z.T); X = np.dot(X, np.linalg.inv(np.diag(Y)))

ここで、何回もSVDする場合には、毎回2乗していると遅くなるので(エラーはたまにしかでないと想定している)、普通にSVDを試みてから、エラーが出たら2乗してSVDするようにしている。
1つ目の方法と2つ目の方法のどちらが高速化は分からない。

3つ目の方法はエルミート行列np.dot(A.T ,A)の対角化に落とす方法です。
（任意の行列Aに対して、特異値分解ができることの証明としてよく使われる方法をそのままコードにした）

try:
	X, Y, Z = np.linalg.svd(A)
except:
	w,v = np.linalg.eigh(np.dot(A.T ,A))
	w = w[::-1]; v = v[:,::-1]
	Y = np.sqrt(w)
	X = np.dot(A,v); X = np.dot(X,np.diag(Y**(-1))); Z = v.T

ただし、Aの形が正方行列出ない時には、形に応じて、上のnp.dot(A.T ,A) or np.dot(A ,A.T)を使わなくてはいけません。適当に書き換えてください。

2つ目の方法より、3つ目の方法の方が安定するようです。