# numerically stable sparse Gaussian process regression (matrix inversion)

In sparse approximations of GP for large data set $$(X,\mathbf{y})$$ with $$n$$ samples, usually $$m$$ inducing points are chosen such that the true covariance matrix is approximated by $$K_{nn}\to K_{nm}K_{mm}^{-1}K_{mn}$$, where $$K_{mn}=K_{nm}^T$$ is the cross-covariance between inducing point and data points. This results in a lower bound for the marginal likelihood (negative loss function) of (see Titsias 2009) $$\begin{gather} F_V = \log\mathcal{N}(0,\sigma^2 I_{nn}+K_{nm}K_{mm}^{-1}K_{mn}) - \frac{1}{2\sigma^2} Tr({K_{nn}-K_{nm}K_{mm}^{-1}K_{mn}}) \end{gather}$$ I have already coded this and the code is fast but it sometimes (50% actually!) fails. Let $$\Gamma=\sigma^2 I_{nn}+K_{nm}K_{mm}^{-1}K_{mn}$$. For evaluating this, we need to calculate (among other things) $$\begin{equation} \mathbf{y}^T \Gamma^{-1} \mathbf{y} + \log |\Gamma| \end{equation}$$ I used Woodbury equations to convert $$\Gamma^{-1} = \sigma^{-2}I_{nn} - \sigma^{-4} K_{nm} \Sigma_{mm}^{-1} K_{mn}$$ where $$\Sigma_{mm} = K_{mm} + \sigma^{-2} K_{mn}K_{nm}$$. But apparently the Cholesky factorization (or inversion) of $$\Sigma$$ is not well conditioned. Do you know any numerically stable and efficient algorithm to do this (invert $$\Gamma$$)?

It is some time ago that this question was asked, but I now also came across this problem and first implemented the inverse like you did. However, this is apparently numerical extremely unstable.

Following Gaussian Processes for Machine Learning I did a more stable implementation, where I first decompose the approximate co-variance matrix in a symmetric way as

$$K_{nm}K_{mm}^{-1}K_{mn}=K_{nm}(L_{mm}L_{mm}^T)^{-1}K_{mn}=K_{nm}(L_{mm}^T)^{-1} L_{mm}^{-1}K_{mn}=\bar{K}_{nm}\bar{K}_{mn}$$

where I used the Cholesky decomposition $$K_{mm}=L_{mm}L_{mm}^T$$ and inverted it to form the intermediate $$\bar{K}_{nm}=L_{mm}^{-1}K_{mn}$$, but you can also use an eigen decomposition to construct a similar intermediate.

Afterwards I could use the matrix inversion lemma to get

$$(\hat{K}_{nm}\hat{K}_{mn}+\sigma^2I)^{-1}=\sigma^{-2}I-\sigma^{-2}\hat{K}_{nm}(\sigma^2I+\hat{K}_{mn}\hat{K}_{nm})^{-1}\hat{K}_{mn}$$

which is at least for my cases numerical more stable. Still, for co-variance matrices with large condition number this still might be not so robust.