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$)?
 A: 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.
