# full rank update to cholesky decomposition for multivariate normal distribution

This question is a specialization of full rank update to cholesky decomposition, to which I hope to get a more positive answer.

When calculating the minus log of the multivariate normal distribution, the most difficult part is evaluating $$f(\epsilon,\Sigma) \equiv \epsilon^\top \Sigma^{-1} \epsilon + \log \det \Sigma$$

where $\epsilon$ is a vector, and $\Sigma$ is a symmetric positive definite matrix. When calculating its gradient, we are also interested in $$\frac{\partial f}{\partial \epsilon}(\epsilon,\Sigma) \equiv 2\Sigma^{-1}\epsilon \qquad \frac{\partial f}{\partial \Sigma}(\epsilon,\Sigma) \equiv -(\Sigma^{-1}\epsilon)(\Sigma^{-1}\epsilon)^\top + \Sigma^{-1}$$

In this question, $\epsilon$ is constant, and $\Sigma$ is such that $$\exists (x_1,\ldots,x_N) \quad \Sigma_{ij} = \sigma^2(x_i,x_j; \theta)$$

$\sigma^2$ is twice differentiable along $\theta$, which is a parameter, or a parameter vector, depending on the application. Therefore, an infinitesimal change in $\theta$ results in an infinitesimal change in $\Sigma$.

In trying to optimize the value of $\theta$, I perform minute changes in its value, possibly following the gradient of $f$, and recompute $f$ at each step to see if the value of $\theta$ got better. I would like to make this step faster by not Cholesky-decomposing $\Sigma(\theta+d\theta)$ from scratch, but instead use the decomposition of $\Sigma(\theta)$ as a starting point. If it makes things easier, I would simply like to update the following quantities

• $\epsilon\Sigma^{-1}\epsilon$
• $\log \det \Sigma$
• $\Sigma^{-1}\epsilon$
• $\Sigma^{-1}$

Here's my personal intuition as to why there should be a solution. We focus on the simplest case, updating $g(\theta) \equiv \epsilon^\top\Sigma^{-1}\epsilon$ and with $\theta$ a scalar parameter.

The question is, therefore, "what is $g(\theta+\delta\theta)$ knowing $g(\theta)$ ?" We know that $g(\theta+d\theta) = g(\theta) + g'(\theta)d\theta$, and in this particular case, we have $$g'(\theta) \equiv -(\Sigma^{-1}\epsilon)^\top \frac{\partial\Sigma}{\partial \theta}(\Sigma^{-1}\epsilon) \qquad \text{with} \quad \Sigma \equiv \Sigma(\theta)$$

Once the Cholesky decomposition of $\Sigma(\theta)$ is known, the calculation of $g'(\theta)$ is $O(N^2)$ and therefore by simple iteration of this gradient descent, we should be able to reach $\Sigma(\theta+\delta\theta)$ in a small number of iterations, making it $O(kN^2)$ with $k\ll N \equiv \dim \Sigma$. So, what should it be, and how to approach this? Steepest descent ? Conjugate Gradients ? And how do I set $k$?

• To clarify, you have two questions: (1) Can you make the computation of $f$ cheaper by not having to perform a Cholesky decomposition of $\Sigma$ for every different parameter value $\theta$ (2) Absent any issues of how fast you can compute $f$, what is the best way to go about finding $\max_\theta f$. Is that about right? – Daniel Shapero Jan 30 '14 at 18:04
• I agree with (1) but (2) is incorrect. It is indeed what I'm trying to do, but in this question I am only asking about (1). Conjugate gradients is O(kN^2) if I'm correct, so that could be a way to find a solution to (1). – yannick Jan 31 '14 at 10:11