title: numpy matrix trickery - sum of inverse times matrices
body: <p>I'm trying to do the following, and repeat until convergence:</p>
  
  <p><img src="https://i.imgur.com/JbF9ech.png" /></p>
  
  <p>where each X<sub>i</sub> is <code>n x p</code>, and there are <code>r</code> of them in an <code>r x n x p</code> array called <code>samples</code>. <code>U</code> is <code>n x n</code>, <code>V</code> is <code>p x p</code>. (I'm getting the MLE of a <a href="http://en.wikipedia.org/wiki/Matrix_normal_distribution" rel="noreferrer">matrix normal distribution</a>.)
  The sizes are all potentially large-ish; I'm expecting things at least on the order of  <code>r = 200</code>, <code>n = 1000</code>, <code>p = 1000</code>.</p>
  
  <p>My current code does</p>
  
  <pre><code>V = np.einsum('aji,jk,akl-&gt;il', samples, np.linalg.inv(U) / (r*n), samples)
  U = np.einsum('aij,jk,alk-&gt;il', samples, np.linalg.inv(V) / (r*p), samples)
  </code></pre>
  
  <p>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.</p>
  
  <p>I could avoid the inverse by doing something like</p>
  
  <pre><code>V = sum(np.dot(x.T, scipy.linalg.solve(A, x)) for x in samples)
  </code></pre>
  
  <p>(or something similar that exploited the psd-ness), but then there's a Python loop, and that makes the numpy fairies cry.</p>
  
  <p>I could also imagine reshaping <code>samples</code> in such a way that I could get an array of <code>A^-1 x</code> using <code>solve</code> for every <code>x</code> without having to do a Python loop, but that makes a big auxiliary array that's a waste of memory.</p>
  
  <p>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 language and calling out to it? (Just porting it directly to Cython might help, but would still involve a lot of Python method calls; but maybe it wouldn't be too much trouble to make the relevant blas/lapack routines directly without too much trouble.)</p>
  
  <p>(As it turns out, I don't actually need the matrices <code>U</code> and <code>V</code> in the end - just their determinants, or actually just the determinant of their Kronecker product. So if anyone has a clever idea for how to do less work and still get the determinants out, that would be much appreciated.)</p>