numpyの行列の策略 - 逆回行列の合計

Question

私は次の操作を行い、収束するまで繰り返すようにしようとしている：

各X _は私 n x pあり、そしてそれらの rがある

samplesと呼ばれるr x n x pの配列に格納されます。 Uはn x n、Vはp x pである。 （私はmatrix normal distributionのMLEを取得しています） サイズはすべて潜在的に大きいishです。私は少なくともr = 200、n = 1000、p = 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-ネスを活用似たような）のような何かをすることによって、逆を避けることができますが、そこPythonのループがだ、それがnumpyの妖精たちが泣きます。

また、私はPythonのループを行うことなく、すべてのxためsolveを使用してA^-1 xの配列を得ることができるようにsamplesを整形想像できるが、それはメモリの無駄だ大きな補助配列になります。

明示的な逆も、Pythonのループも、大きなaux配列もありません。この3つの中で最良のものを得るためにできる線形代数またはnumpyのトリックはありますか？それとも、Pythonループを高速言語で実装して呼び出すのが最善でしょうか？ （これを直接Cythonに移植するだけで、多くのPythonメソッド呼び出しが必要になりますが、あまり問題なく、関連するblas/lapackルーチンを直接作成するのはあまり問題にならないでしょう）。

（結局のところ、私は実際に行列UとVを最終的には必要としません - その決定要因、または実際にはクロネッカー製品の決定要因になります。決定を出す、それははるかに高く評価されるだろう。）誰かがより多くのインスピレーションの答えを思いつくまでは私があなただったら

Answer 1

が、私は妖精が泣いてみましょうと思います...

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ミリ秒以上かかることがあります。ルーピングと加算が自由であれば、5％未満の改善が得られます。 Python関数のオーバーヘッドもあるかもしれませんが、あなたがコードしている言語にかかわらず、方程式を解く実際の時間に対して無視できるほどです。