12

私は次のことを行おうとしています、そして収束するまで繰り返します：

ここで、各X iはn x pであり、 。と呼ばれる配列rにそれらがあります。です、です。（行列正規分布のMLEを取得しています。）サイズはすべて潜在的に大きいです。私は少なくとも、、、のオーダーで物事を期待して います。r x n x psamplesUn x nVp x pr = 200n = 1000p = 1000

私の現在のコードは

V = np.einsum('aji,jk,akl->il', samples, np.linalg.inv(U) / (r*n), samples)
U = np.einsum('aij,jk,alk->il', samples, np.linalg.inv(V) / (r*p), samples)

これは問題なく機能しますが、もちろん、実際に逆を見つけてそれを乗算することは想定されていません。UとVが対称で正定値であるという事実をどうにかして利用できれば、それも良いことです。反復でUとVのコレスキー係数を計算できるようにしたいと思いますが、合計のため、その方法がわかりません。

私は次のようなことをすることで逆を避けることができます

V = sum(np.dot(x.T, scipy.linalg.solve(A, x)) for x in samples)

（またはpsd-nessを利用した類似のもの）が、Pythonループがあり、それは厄介な妖精を泣かせます。

また、Pythonループを実行せずに、すべてに使用samplesする配列を取得できるような方法で再形成することも想像できますが、それはメモリの浪費である大きな補助配列になります。A^-1 xsolvex

3つすべてを最大限に活用するために実行できる線形代数または厄介なトリックはありますか？明示的な逆関数、Pythonループ、および大きなaux配列はありませんか？それとも、Pythonループを備えたものをより高速な言語で実装し、それを呼び出すのが最善の策ですか？（Cythonに直接移植するだけでも役立つ場合がありますが、それでも多くのPythonメソッド呼び出しが必要になります。ただし、関連するblas / lapackルーチンをそれほど問題なく直接作成するのはそれほど問題ではないでしょう。）

（結局のところ、私は実際には行列式を必要とせずU、V最終的には行列式だけ、または実際にはクロネッカー積の行列式だけを必要とします。したがって、誰かがより少ない作業を行い、それでも決定要因が出て、それは大いにありがたいです。）

4

1 に答える 1

8

誰かがもっと刺激的な答えを思い付くまで、もし私があなたなら、妖精たちを泣かせたでしょう...

r, n, p = 200, 400, 400

X = np.random.rand(r, n, p)
U = np.random.rand(n, n)

In [2]: %timeit np.sum(np.dot(x.T, np.linalg.solve(U, x)) for x in X)
1 loops, best of 3: 9.43 s per loop

In [3]: %timeit np.dot(X[0].T, np.linalg.solve(U, X[0]))
10 loops, best of 3: 45.2 ms per loop

したがって、Pythonループがあり、すべての結果を合計する必要がある場合、解決する必要のある200のシステムのそれぞれを解決するのにかかる時間の200倍以上の390ミリ秒かかります。ループと合計が無料の場合、5％未満の改善しか得られません。Python関数のオーバーヘッドを呼び出すこともあるかもしれませんが、どの言語でコーディングしても、方程式を解く実際の時間に対してはおそらく無視できるでしょう。

于 2013-02-09T01:38:32.033 に答える