I'm trying to do the following, and repeat until convergence:
where each Xi 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