numpy的矩陣掛羊頭賣狗肉 - 逆時代矩陣的總和

Question

我試圖做到以下幾點，並重復直至收斂：

其中X _我是n x p，並有他們的r在r x n x p數組中被稱爲samples。 U是n x n,V是p x p。 （我得到了matrix normal distribution的MLE。） 這些尺寸都可能很大，我期待的事情至少在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的Cholesky因子，但是我不知道如何去做，因爲總和。

我可以做這樣的事情

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

（或利用的PSD岬類似的東西）避免逆，但後來有一個Python的循環，這使得該numpy的仙女哭泣。

我也想象不出這樣一種方式，我能得到使用solve爲每x的A^-1 x數組，而不必做一個Python循環重塑samples，但能造成很大的輔助陣列是浪費內存。

是否存在一些線性代數或無節制的技巧，我可以做到最好的三種：沒有明確的逆，沒有Python循環，沒有大的輔助數組？或者，我最好的辦法是用更快的語言來實現Python循環並向其發出呼籲？ （直接將它直接移植到Cython可能會有所幫助，但仍然會涉及到很多Python方法調用;但是，直接製作相關的blas/lapack例程並不會太麻煩。）

（事實證明，實際上我並不需要矩陣U和V - 只是它們的決定因素，或者實際上只是Kronecker產品的決定因素。所以如果有人對如何減少工作有一個聰明的想法，仍然得到了決定了，那將是非常讚賞。）

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多倍多390倍。如果循環和求和是免費的，你會得到少於5％的改進。可能會有一些調用python函數的開銷，但對於實際解決方程的時間來說，它可能仍然可以忽略不計，無論您使用哪種語言編寫。