# numerical linear algebra tricks for repeated sums and inversions with symmetric positive-definite matrices

I'm doing the following procedure to get the max-likelihood estimate of a matrix-variate normal distribution from $r$ samples of matrices in $\mathbb{R}^{n \times p}$ (algorithm from Dutilleul (1999)):

\begin{align} V &= \frac{1}{r n} \sum_{i=1}^r X_i^T U^{-1} X_i \\ U &= \frac{1}{r p} \sum_{i=1}^r X_i V^{-1} X_i^T \end{align}

and repeating until convergence. ($U \in \mathbb{R}^{p \times p}$, $V \in \mathbb{R}^{n \times n}$). I need to do this for moderately large $n$ and $p$ (on the order of 1,000 each), with somewhat smaller $r$ (order of 200), and I need to do it over and over and over again, so any algorithmic speed increases would be very handy.

That paper didn't really talk about the practical implementation of the algorithm, though. Since $U$ and $V$ are symmetric positive-definite, it seems that there should be some trick smarter than just doing the sums as written (using a solve call rather than an inverse and multiply for the right half, of course).

Currently, I'm not exploiting the SPD-ness at all. I can Cholesky-factorize $U$ and $V$ before doing each sum and then using a solver that exploits that, which will presumably save some work in the sum over $i$.

1. Is it possible to just directly get a Cholesky-style factor of $U$ and $V$ in the sum, rather than the full matrix? The outer-product-with-an-spd-matrix-in-the-middle form is super tempting for that, but of course the sum means you can't just add together the Cholesky factors of each summand.

2. Can I exploit some structure here, more than just doing $X_i$ times the result of a LU solver for $U^{-1} X_i$? Would it maybe be better to do $L_U^{-1} X_i$ and outer product that with itself, for example, or something smarter?

3. As it turns out, I don't actually need $U$ and $V$, just $\det(V \otimes U) = \det(U)^n \det(V)^p$. Is there some shortcut I can take that'll get the determinants out without having to do as much work? Seems unlikely, because of the iterative nature, but people out there are smarter than me. :)

(Also posted with some additional programming-related components on StackOverflow.)

-
Did you try computing $U^{-1}X$ by using an iterative method such as Gauss-Seidel? Depending on the accuracy that you need, it may be faster than solving directly. Gauss-Seidel is guaranteed to converge for SPD matrices. – Victor May Jan 20 '14 at 14:03
@VictorMay Thanks for the suggestion. It ended up that just caching the Cholesky factors made it fast enough to no longer be the bottleneck in my algorithm, so I left it at that; it may well be that an approximate iterative solution would speed things up more than that, though. – Dougal Jan 20 '14 at 18:20