full rank update to cholesky decomposition for multivariate normal distribution

Question

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$?

Answer 1

in general there is no formula for updating the cholesky factorization of the covariance matrix when the hyperparameter is changed. There are other techniques for speeding up the computations though. For instance the method SKI (http://proceedings.mlr.press/v37/wilson15.pdf) can be used to compute matrix-vector products with the covariance matrix. Lanczos can also be used to compute fast estimates of the posterior covariance matrix (https://arxiv.org/pdf/1803.06058.pdf).

The Gpytorch website has a lot of references. https://docs.gpytorch.ai/en/v1.2.1/examples/02_Scalable_Exact_GPs/index.html