numpy oszustwo matrix - suma macierzy odwrotnych razy

Question

Próbuję wykonać następujące czynności i powtarzać aż do zbieżności:

gdzie każdy X _i jest n x p i istnieje r z nich w macierzy r x n x p o nazwie samples. U jest n x n, jest p x p. (Dostaję MLE od matrix normal distribution.) Rozmiary są potencjalnie duże-ish; Oczekuję rzeczy przynajmniej na zamówienie: r = 200, n = 1000, p = 1000.

Mój obecny kod robi

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) 

Działa to dobrze, ale oczywiście nigdy nie powinniśmy się rzeczywiście znaleźć odwrotność i pomnożyć przez niego rzeczy. Byłoby dobrze, gdybym mógł w jakiś sposób wykorzystać fakt, że U i V są symetryczne i pozytywne. Chciałbym móc po prostu obliczyć współczynnik Cholesky U i V w iteracji, ale nie wiem jak to zrobić z powodu sumy.

mogłem uniknąć odwrotność wykonując coś jak

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

(lub coś podobnego, że wykorzystał PSD-ności), ale to nie jest to pętla Python, a to sprawia, że ​​NumPy wróżki płakać.

mogłem również wyobrazić przekształcanie samples w taki sposób, że mogę dostać tablicę A^-1 x korzystając solve dla każdego x bez konieczności wykonywania pętli Pythona, ale to robi dużą tablicę pomocniczą, że to strata pamięci.

Czy istnieje algebra liniowa lub sztuczka numpy, którą mogę zrobić, aby uzyskać najlepsze z wszystkich trzech: brak wyraźnych odwrotności, brak pętli w Pythonie i brak dużych tablic Aux? A może mój najlepszy zakład implementacji tego z pętlą Python w szybszym języku i wywoływanie go? (Po prostu przeniesienie go bezpośrednio do Cythona mogłoby pomóc, ale nadal wymagałoby wielu wywołań w języku Python, ale może nie byłoby zbyt wielkim problemem, aby odpowiednie procedury blas/lapack były wykonywane bezproblemowo.)

(Jak się okazuje, naprawdę nie potrzebuję w końcu matryc U i V - tylko ich wyznaczniki, a właściwie tylko wyznacznik ich produktu Kroneckera. Więc jeśli ktoś ma sprytny pomysł na to, jak wykonać mniej pracy i wciąż dostać się uwarunkowania, które byłyby mile widziane.)

Answer 1

Dopóki ktoś nadjeżdża z bardziej natchnionym odpowiedź, gdybym był tobą, bym pozwolił wróżki płakać ...

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 

A więc posiadanie pętli w języku Pythona i sumowanie wszystkich wyników zajmuje 390 ms więcej niż 200 razy czasu potrzebnego na rozwiązanie każdego z 200 systemów, które należy rozwiązać. Dostaniesz mniej niż 5% poprawy, jeśli pętla i suma będą darmowe. Możliwe, że niektóre z nich wywołują również funkcję Pythona, ale prawdopodobnie nadal będzie ona pomijalna w stosunku do rzeczywistego czasu rozwiązywania równań, bez względu na język, w którym się je koduje.