Edit Question - KnowledgeHub

Title Body

I'm trying to do the following, and repeat until convergence: where each X i is n x p , and there are r of them in an r x n x p array called samples . U is n x n , V is p x p . (I'm getting the MLE of a matrix normal distribution .) The sizes are all potentially large-ish; I'm expecting things at least on the order of r = 200 , n = 1000 , p = 1000 . My current code does 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) This works okay, but of course you're never supposed to actually find the inverse and multiply stuff by it. It'd also be good if I could somehow exploit the fact that U and V are symmetric and positive-definite. I'd love to be able to just calculate the Cholesky factor of U and V in the iteration, but I don't know how to do that because of the sum. I could avoid the inverse by doing something like V = sum(np.dot(x.T, scipy.linalg.solve(A, x)) for x in samples) (or something similar that exploited the psd-ness), but then there's a Python loop, and that makes the numpy fairies cry. I could also imagine reshaping samples in such a way that I could get an array of A^-1 x using solve for every x without having to do a Python loop, but that makes a big auxiliary array that's a waste of memory. Is there some linear algebra or numpy trick I can do to get the best of all three: no explicit inverses, no Python looping, and no big aux arrays? Or is my best bet implementing the one with a Python loop in a faster lan

Tags (comma-separated) Cancel