2
$\begingroup$

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

$\endgroup$
1
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.