NumPy matrice inganno - somma dei tempi inversi matrici

Question

sto cercando di effettuare le seguenti operazioni, e ripetere fino convergenza:

dove ciascun X _i è n x p, e ci sono r di loro in un array r x n x p chiamato samples. U è n x n, V è p x p. (Sto ricevendo il MLE di un . Le dimensioni sono tutte potenzialmente di grandi dimensioni; Mi aspetto cose almeno nell'ordine di r = 200, n = 1000, p = 1000.

mio codice corrente fa

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) 

Questo funziona bene, ma naturalmente non si è mai dovuto trovare in realtà l'inverso e moltiplicare roba da essa. Sarebbe anche bello se potessi in qualche modo sfruttare il fatto che U e V sono simmetrici e positivi-definiti. Mi piacerebbe essere in grado di calcolare il fattore Cholesky di U e V nell'iterazione, ma non so come farlo a causa della somma.

ho potuto evitare l'inverso facendo qualcosa di simile

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

(o qualcosa di simile che sfruttava il PSD-ness), ma poi c'è un ciclo di Python, e che rende le fate numpy piangono.

I può anche ipotizzare rimodellare samples in modo tale da poter ottenere un array di A^-1 x usando solve indipendentemente dai x senza dover fare un loop Python, ma che rende un array ausiliario grande che una perdita di memoria.

C'è qualche trucco lineare algebrico o numpico che posso fare per ottenere il meglio da tutti e tre: nessuna inversione esplicita, nessun loop Python e nessun grande array aux? Oppure la mia scommessa migliore è implementare quella con un loop Python in un linguaggio più veloce e chiamarla? (Portarlo direttamente su Cython potrebbe essere d'aiuto, ma comporterebbe comunque molte chiamate al metodo Python, ma forse non sarebbe troppo difficile rendere le routine blas/lapack direttamente senza troppi problemi.)

(Come risulta, in realtà non ho bisogno delle matrici U e V alla fine - solo i loro determinanti, o in realtà solo il determinante del loro prodotto Kronecker. Quindi, se qualcuno ha un'idea intelligente di come fare meno lavoro e ancora ottenere i determinanti fuori, che sarebbe molto apprezzato.)

Answer 1

fino a quando qualcuno esce con una risposta più ispirato, se fossi in te, mi piacerebbe lasciare le fate piangere ...

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 

Quindi avere un loop python e dover sommare tutti i risultati insieme, sta prendendo 390 ms più di 200 volte quello che serve per risolvere ciascuno dei 200 sistemi che devono essere risolti. Otterresti meno di un miglioramento del 5% se il ciclo e la somma fossero gratuiti. Potrebbero esserci alcuni richiami di overhead della funzione python, ma probabilmente sarà ancora trascurabile rispetto al tempo reale di risoluzione delle equazioni, indipendentemente dalla lingua in cui viene codificato.