Derivation of SGPR equations

James Hensman, March 2016. Corrections by Alex Matthews, December 2016

This notebook contains a derivation of the form of the equations for the marginal likelihood bound and predictions for the sparse Gaussian process regression model in GPflow, gpflow.models.SGPR.

The primary reference for this work is Titsias 2009 [1], though other works (Hensman et al. 2013 [2], Matthews et al. 2016 [3]) are useful for clarifying the prediction density.

Marginal likelihood bound

The bound on the marginal likelihood (Titsias 2009) is:

$$\log p(\mathbf{y})\ge \log \mathcal{N}(\mathbf{y}|\mathbf{0},{\mathbf{Q}}_{ff}+{\sigma }^2\mathbf{I})-\frac{1}{2}{\sigma }^{-2}\text{tr}({\mathbf{K}}_{ff}-{\mathbf{Q}}_{ff})\triangleq \mathcal{L}$$ where: $${\mathbf{Q}}_{ff}={\mathbf{K}}_{fu}{\mathbf{K}}_{uu}^{-1}{\mathbf{K}}_{uf}$$

The kernel matrices ${\mathbf{K}}_{ff}$, ${\mathbf{K}}_{uu}$, ${\mathbf{K}}_{fu}$ represent the kernel evaluated at the data points $\mathbf{X}$, the inducing input points $\mathbf{Z}$, and between the data and inducing points respectively. We refer to the value of the GP at the data points $\mathbf{X}$ as $\mathbf{f}$, at the inducing points $\mathbf{Z}$ as $\mathbf{u}$, and at the remainder of the function as $f^{\star }$.

To obtain an efficient and stable evaluation on the bound $\mathcal{L}$, we first apply the Woodbury identity to the effective covariance matrix:

$$[{\mathbf{Q}}_{ff}+{\sigma }^2\mathbf{I}{]}^{-1}={\sigma }^{-2}\mathbf{I}-{\sigma }^{-4}{\mathbf{K}}_{fu}[{\mathbf{K}}_{uu}+{\mathbf{K}}_{uf}{\mathbf{K}}_{fu}{\sigma }^{-2}{]}^{-1}{\mathbf{K}}_{uf}$$

Now, to obtain a better conditioned matrix for inversion, we rotate by $\mathbf{L}$, where $\mathbf{L}{\mathbf{L}}^{\mathrm{\top }}={\mathbf{K}}_{uu}$:

$$[{\mathbf{Q}}_{ff}+{\sigma }^2\mathbf{I}{]}^{-1}={\sigma }^{-2}\mathbf{I}-{\sigma }^{-4}{\mathbf{K}}_{fu}{\mathbf{L}}^{-\mathrm{\top }}{\mathbf{L}}^{\mathrm{\top }}[{\mathbf{K}}_{uu}+{\mathbf{K}}_{uf}{\mathbf{K}}_{fu}{\sigma }^{-2}{]}^{-1}\mathbf{L}{\mathbf{L}}^{-1}{\mathbf{K}}_{uf}$$

This matrix is better conditioned because, for many kernels, it has eigenvalues bounded above and below. For more details, see section 3.4.3 of Gaussian Processes for Machine Learning.

$$\phantom{[{\mathbf{Q}}_{ff}+{\sigma }^2\mathbf{I}{]}^{-1}}={\sigma }^{-2}\mathbf{I}-{\sigma }^{-4}{\mathbf{K}}_{fu}{\mathbf{L}}^{-\mathrm{\top }}[{\mathbf{L}}^{-1}({\mathbf{K}}_{uu}+{\mathbf{K}}_{uf}{\mathbf{K}}_{fu}{\sigma }^{-2}){\mathbf{L}}^{-\mathrm{\top }}{]}^{-1}{\mathbf{L}}^{-1}{\mathbf{K}}_{uf}$$

$$\phantom{[{\mathbf{Q}}_{ff}+{\sigma }^2\mathbf{I}{]}^{-1}}={\sigma }^{-2}\mathbf{I}-{\sigma }^{-4}{\mathbf{K}}_{fu}{\mathbf{L}}^{-\mathrm{\top }}[\mathbf{I}+{\mathbf{L}}^{-1}\mathbf{(}{\mathbf{K}}_{uf}{\mathbf{K}}_{fu}){\mathbf{L}}^{-\mathrm{\top }}{\sigma }^{-2}{]}^{-1}{\mathbf{L}}^{-1}{\mathbf{K}}_{uf}$$

For notational convenience, we’ll define ${\mathbf{L}}^{-1}{\mathbf{K}}_{uf}{\sigma }^{-1}\triangleq \mathbf{A}$, and $[\mathbf{I}+\mathbf{A}{\mathbf{A}}^{\mathrm{\top }}]\triangleq \mathbf{B}$:

$$\phantom{[{\mathbf{Q}}_{ff}+{\sigma }^2\mathbf{I}{]}^{-1}}={\sigma }^{-2}\mathbf{I}-{\sigma }^{-2}{\mathbf{A}}^{\mathrm{\top }}[\mathbf{I}+\mathbf{A}{\mathbf{A}}^{\mathrm{\top }}{]}^{-1}\mathbf{A}$$

$$\phantom{[{\mathbf{Q}}_{ff}+{\sigma }^2\mathbf{I}{]}^{-1}}={\sigma }^{-2}\mathbf{I}-{\sigma }^{-2}{\mathbf{A}}^{\mathrm{\top }}{\mathbf{B}}^{-1}\mathbf{A}$$

We also apply the matrix determinant lemma to the same:

:nbsphinx-math:`begin{equation} |{\mathbf Q_{ff}} + \sigma^2 {\mathbf I}| = |{mathbf K_{uu}} +

mathbf K_{uf}mathbf K_{fu}sigma^{-2}| , |\mathbf K_{uu}^{-1}| , |\sigma^{2}\mathbf I|

end{equation}`

Substituting ${\mathbf{K}}_{uu}={\mathbf{L}\mathbf{L}}^{\mathrm{\top }}$: :nbsphinx-math:`begin{equation} |{\mathbf Q_{ff}} + \sigma^2 {\mathbf I}| = |{mathbf {L L}^top} +

mathbf K_{uf}mathbf K_{fu}sigma^{-2}| , |\mathbf L^{-\top}|,| mathbf L^{-1}| , |\sigma^{2}\mathbf I|

end{equation}`

:nbsphinx-math:`begin{equation} |{\mathbf Q_{ff}} + \sigma^2 {\mathbf I}| = |mathbf I +

mathbf L^{-1}mathbf K_{uf}mathbf K_{fu} mathbf L^{-top}sigma^{-2}| , |\sigma^{2}\mathbf I|

end{equation}`

$$|{\mathbf{Q}}_{ff}+{\sigma }^2\mathbf{I}|=|\mathbf{B}||{\sigma }^2\mathbf{I}|$$

With these two definitions, we’re ready to expand the bound:

$$\mathcal{L}=\log \mathcal{N}(\mathbf{y}|\mathbf{0},{\mathbf{Q}}_{ff}+{\sigma }^2\mathbf{I})-\frac{1}{2}{\sigma }^{-2}\text{tr}({\mathbf{K}}_{ff}-{\mathbf{Q}}_{ff})$$

$$=-\frac{N}{2}\log 2\pi -\frac{1}{2}\log |{\mathbf{Q}}_{ff}+{\sigma }^2\mathbf{I}|-\frac{1}{2}{\mathbf{y}}^{\mathrm{\top }}[{\mathbf{Q}}_{ff}+{\sigma }^2\mathbf{I}{]}^{-1}\mathbf{y}-\frac{1}{2}{\sigma }^{-2}\text{tr}({\mathbf{K}}_{ff}-{\mathbf{Q}}_{ff})$$

$$=-\frac{N}{2}\log 2\pi -\frac{1}{2}\log (|\mathbf{B}||{\sigma }^2\mathbf{I}|)-\frac{1}{2}{\sigma }^{-2}{\mathbf{y}}^{\mathrm{\top }}(\mathbf{I}-{\sigma }^{-2}{\mathbf{A}}^{\mathrm{\top }}{\mathbf{B}}^{-1}\mathbf{A})\mathbf{y}-\frac{1}{2}{\sigma }^{-2}\text{tr}({\mathbf{K}}_{ff}-{\mathbf{Q}}_{ff})$$

$$=-\frac{N}{2}\log 2\pi -\frac{1}{2}\log |\mathbf{B}|-\frac{N}{2}\log {\sigma }^2-\frac{1}{2}{\sigma }^{-2}{\mathbf{y}}^{\mathrm{\top }}\mathbf{y}+\frac{1}{2}{\sigma }^{-2}{\mathbf{y}}^{\mathrm{\top }}{\mathbf{A}}^{\mathrm{\top }}{\mathbf{B}}^{-1}\mathbf{A}\mathbf{y}-\frac{1}{2}{\sigma }^{-2}\text{tr}({\mathbf{K}}_{ff})+\frac{1}{2}\text{tr}({\mathbf{AA}}^{\mathrm{\top }})$$

where ${\sigma }^{-2}\text{tr}(\mathbf{Q})=\text{tr}({\mathbf{AA}}^{\mathrm{\top }})$.

Finally, we define $\mathbf{c}\triangleq {\mathbf{L}}_{\mathbf{B}}^{-1}\mathbf{A}\mathbf{y}{\sigma }^{-1}$, with $\text{Double subscripts: use braces to clarify}$, so that:

$${\sigma }^{-2}{\mathbf{y}}^{\mathrm{\top }}{\mathbf{A}}^{\mathrm{\top }}{\mathbf{B}}^{-1}\mathbf{A}\mathbf{y}={\mathbf{c}}^{\mathrm{\top }}\mathbf{c}$$

The SGPR code implements this equation with small changes for multiple concurrent outputs (columns of the data matrix Y), and also a prior mean function.

Prediction

At prediction time, we need to compute the mean and variance of the variational approximation at some new points ${\mathbf{X}}^{\star }$.

Following Hensman et al. (2013), we know that all the information in the posterior approximation is contained in the Gaussian distribution $q(\mathbf{u})$, which represents the distribution of function values at the inducing points $\mathbf{Z}$. Remember that:

$$q(\mathbf{u})=\mathcal{N}(\mathbf{u}|\mathbf{m},{\mathbf{\Lambda }}^{-1})$$

with:

$$\mathbf{\Lambda }={\mathbf{K}}_{uu}^{-1}+{\mathbf{K}}_{uu}^{-1}{\mathbf{K}}_{uf}{\mathbf{K}}_{fu}{\mathbf{K}}_{uu}^{-1}{\sigma }^{-2}$$

$$\mathbf{m}={\mathbf{\Lambda }}^{-1}{\mathbf{K}}_{uu}^{-1}{\mathbf{K}}_{uf}\mathbf{y}{\sigma }^{-2}$$

To make a prediction, we need to integrate:

$$p({\mathbf{f}}^{\star })=\int p({\mathbf{f}}^{\star }|\mathbf{u})q(\mathbf{u})\text{d}\mathbf{u}$$

with:

$$p({\mathbf{f}}^{\star }|\mathbf{u})=\mathcal{N}({\mathbf{f}}^{\star }|{\mathbf{K}}_{\star u}{\mathbf{K}}_{uu}^{-1}\mathbf{u},{\mathbf{K}}_{\star \star }-{\mathbf{K}}_{\star u}{\mathbf{K}}_{uu}^{-1}{\mathbf{K}}_{u\star })$$

The integral results in:

$$p({\mathbf{f}}^{\star })=\mathcal{N}({\mathbf{f}}^{\star }|{\mathbf{K}}_{\star u}{\mathbf{K}}_{uu}^{-1}\mathbf{m},{\mathbf{K}}_{\star \star }-{\mathbf{K}}_{\star u}{\mathbf{K}}_{uu}^{-1}{\mathbf{K}}_{u\star }+{\mathbf{K}}_{\star u}{\mathbf{K}}_{uu}^{-1}{\mathbf{\Lambda }}^{-1}{\mathbf{K}}_{uu}^{-1}{\mathbf{K}}_{u\star })$$

Note from our above definitions we have:

$${\mathbf{K}}_{uu}^{-1}{\mathbf{\Lambda }}^{-1}{\mathbf{K}}_{uu}^{-1}={\mathbf{L}}^{-\mathrm{\top }}{\mathbf{B}}^{-1}{\mathbf{L}}^{-1}$$

and further:

$${\mathbf{K}}_{uu}^{-1}\mathbf{m}={\mathbf{L}}^{-\mathrm{\top }}{\mathbf{L}}_{\mathbf{B}}^{-\mathrm{\top }}\mathbf{c}$$

substituting:

$$p({\mathbf{f}}^{\star })=\mathcal{N}({\mathbf{f}}^{\star }|{\mathbf{K}}_{\star u}{\mathbf{L}}^{-\mathrm{\top }}{\mathbf{L}}_{\mathbf{B}}^{-\mathrm{\top }}\mathbf{c},{\mathbf{K}}_{\star \star }-{\mathbf{K}}_{\star u}{\mathbf{L}}^{-\mathrm{\top }}(\mathbf{I}-{\mathbf{B}}^{-1}){\mathbf{L}}^{-1}{\mathbf{K}}_{u\star })$$

The code in SGPR implements this equation, with an additional switch depending on whether the full covariance matrix is required.

References

[1] Titsias, M: Variational Learning of Inducing Variables in Sparse Gaussian Processes, PMLR 5:567-574, 2009

[2] Hensman et al: Gaussian Processes for Big Data, UAI, 2013

[3] Matthews et al: On Sparse Variational Methods and the Kullback-Leibler Divergence between Stochastic Processes, AISTATS, 2016