# numerically stable sparse gasussian 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) $$$$\mathbf{y}^T \Gamma^{-1} \mathbf{y} + \log |\Gamma|$$$$ 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$$)?