動作環境
本記事の記載内容は以下の環境で検証しています。
- macOS Sierra
- CPU: Intel Core i7 3615QM (2.3GHz, 物理4コア)
- RAM: DDR3 16GB
- Python 3.6.6
- NumPy 1.14.2
multivariate_normal関数が妙に遅いし不安定
多変量正規分布に従う乱数を使ったシミュレーションを行なうことが多い筆者。最近はPythonを使用しているため、ここのところnumpyのmultivariate_normalに頻繁にお世話になっていました。分散共分散関数を引数とするだけで勝手に乱数を生成してくれるスグレモノ。
参考: numpy.random.multivariate_normal
...が、これが妙に遅い。特に分散共分散行列の規模が大きくなった際の挙動が酷くて、数千×数千の行列を引数とすると1回の乱数生成に数十秒かかることもザラにありました。
しかも、時々以下のようなエラーが表示されプログラムが停止することも。
LinAlgError: SVD did not converge SVD
※数千回に1度程度。モンテカルロシミュレーションで繰り返し試行することが多いので...
np.linalg.choleskyを使って自力で書いた方が早そう
どうもこの関数、上記のエラーからもわかる通り内部で特異値分解(SVD)を行なっているそうで... ここで重くなっているようです。SVDによって、関数内で分散共分散行列が半正定値かどうかの判定も行なっているようですね。しかも高速化のため繰り返し処理による近似演算を行なっているらしく、計算結果が収束しない場合にSVD失敗と判定してしまっている模様。(この辺は割愛)
ただ、入力した行列の性質がわかっている身にとってこの処理は無駄of無駄。最も素朴な作法に従い、コレスキー分解に基づく方法を使って自力で書いてみました。
参考: 統計数理研究所
幸い、NumPyにはコレスキー分解を行なうための関数が用意されているためこれを使えばラクに書き直せそうです。
比較用の2つの関数
- 普通にmultivariate_normalを使った場合:
def func1(cov):
return np.random.multivariate_normal(np.zeros(len(cov)), cov)
※covはN*Nの分散共分散行列, np.zeros(len(cov))は期待値が0であることを明示するためのベクトル
※len(cov)はNとイコール
- コレスキー分解を使った場合:
def func2(cov):
L = np.linalg.cholesky(cov) #コレスキー分解
z = np.random.standard_normal(len(cov)) #標準正規乱数ベクトル
return np.dot(L, z) #行列Lとベクトルzの積
速度の比較
適当に生成したN×Nの分散共分散行列に従い、10回繰り返し乱数生成を行なった場合の計算時間を表示しました。行列はシミュレーション開始時に1つだけ生成し、以降はこれを使い回すことで、func1, func2の比較を行なっています。コードの全体は文末に記します。
| N | func1 | func2 |
|---|---|---|
| 64 | 0.04807114601135254s | 0.001219034194946289s |
| 128 | 0.07236790657043457s | 0.004015207290649414s |
| 1024 | 5.460850238800049s | 0.31621503829956055s |
いや、軽く1ケタ時間は高速化されてんぞ...
シンプルイズザベスト。
※計算時間は分散共分散行列の値に依存します。同じプログラムでも試行ごとに上記計算時間の関係とは異なる点にご注意ください。
使用したコード
以下に今回使用したコードの、上記2つの関数以外の箇所を示します。
分散共分散関数の生成時には、予め２次元座標上からランダムにN地点を決めました。地点間の距離に従って相関係数が指数的に減衰していき、かつ分散は場所に依存せず一定という仮定のもとで行列を生成しています。筆者は空間統計勢なので、よくこのような流れを踏んだシミュレーションを組んでいます。
(もっと具体的に言うと、電波伝搬におけるシャドウイングの空間相関に関するモデルです...が、どマニアックな話になるのでこれも割愛。とりあえず、N×Nの分散共分散行列が生成されていると考えていただければOKです)
import numpy as np
import time
#分散共分散関数の生成
def genCovMat(x, y, dcor, sigma):
dmat = np.zeros([len(x), len(y)]) #サンプル間の距離に関する行列
for i in range(len(x)):
dmat[i] = distance(x[i], y[i], x[:], y[:])
h = np.exp(-dmat * np.log(2.0) / dcor) #相関行列
return sigma * sigma * h #分散共分散行列
#2地点間のユークリッド距離の計算
def distance(x1, y1, x2, y2):
return np.sqrt((x1-x2)**2 + (y1-y2)**2)
#観測位置の取得(x, yとも一様分布)
def genMeasurementLocation(size, len_area):
x = np.random.uniform(0.0, len_area, size)
y = np.random.uniform(0.0, len_area, size)
return x, y
if __name__ == '__main__':
LEN_AREA = 5000.0 #２次元平面１辺の長さ[m]
DCOR = 100.0 #空間相関の距離[m]
N = 64 #観測地点の数(=生成するサンプル数)
STDEV = 8.0 #標準偏差
LOOP = 10 #繰り返し回数
x, y = genMeasurementLocation(N, LEN_AREA) #観測地点の決定
cov = genCovMat(x, y, DCOR, STDEV) #分散共分散関数の生成
print("== np.random.multivariate_normal ==")
start = time.time()
for i in range(LOOP):
func1(cov)
print(time.time() - start)
print("")
print("== np.linalg.cholesky ==")
start = time.time()
for i in range(LOOP):
func2(cov)
print(time.time() - start)
まとめ
引数とする分散共分散行列の性質がきちんとわかっているのであれば、コレスキー分解から自力で書いた方がよさそうです。(というか不明確な場合ってどれだけあるんだ...)
みんなでビッグデータしようぜ
お礼
本記事の内容の検証には、@kojiis さんにご協力いただきました。ありがとうございます。
参考
StackOverflowで同様の内容に関する議論が行われていました。
https://stackoverflow.com/questions/40100340/numpy-multivariate-normal-bug-when-dimension-too-high