{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SGPR notes\n",
    "--\n",
    "*James Hensman, March 2016*\n",
    "*Small corrections Alex Matthews, December 2016*\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The GPflow method SGPR performs Sparse Gaussian Process Regression. It is based mainly on Titsias (2009), though other works (Hensman et al. (2013), Matthews et. al. 2016) are useful for clarifying the prediction density.\n",
    "\n",
    "This notebook contains a derivation of the form of the equations used in GPflow for the marginal likelihood bound and for the prediction equations. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Marginal Likelihood bound\n",
    "The bound on the marginal likelihod is (Titsias 2009):\n",
    "$$ \\log p(\\mathbf y) \\geq \\log \\mathcal N(\\mathbf y\\,|\\,\\mathbf 0,\\, \\mathbf Q_{ff} + \\sigma^2 \\mathbf I) - \\tfrac{1}{2} \\sigma^{-2}\\textrm{tr}(\\mathbf K_{ff} - \\mathbf Q_{ff}) \\triangleq \\mathcal L$$\n",
    "with \n",
    "$$\\mathbf Q_{ff} = \\mathbf K_{fu}\\mathbf K_{uu}^{-1}\\mathbf K_{uf}$$.\n",
    "\n",
    "The kernel matrices $\\mathbf K_{ff}$, $\\mathbf K_{uu}$, $\\mathbf K_{fu}$ represent the kernel evaulated 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$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To obtain an efficient and stable evaluation os the bound $\\mathcal L$, we first apply the Woodbury identity to the effective covariance matrix:\n",
    "\n",
    "$$[\\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}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, to obtain a better conditioned matrix for inversion, we rotate by $\\mathbf L$, where $\\mathbf L\\mathbf L^\\top = \\mathbf K_{uu}$:\n",
    "$$[\\mathbf Q_{ff} + \\sigma^2 \\mathbf I ]^{-1} = \\sigma^{-2} \\mathbf I - \\sigma^{-4} \\mathbf K_{fu}\\color{red}{\\mathbf L^{-\\top} \\mathbf L^\\top}[\\mathbf K_{uu} + \\mathbf K_{uf}K_{fu}\\sigma^{-2}]^{-1}\\color{red}{\\mathbf L \\mathbf L^{-1}}\\mathbf K_{uf}$$\n",
    "\n",
    "This matrix is better conditioned because for many kernels it has eigenvalues which are bounded above and below. For more detail see GPML section 3.4.3 ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\phantom{[\\mathbf Q_{ff} + \\sigma^2 \\mathbf I ]^{-1}} = \\sigma^{-2} \\mathbf I - \\sigma^{-4} \\mathbf K_{fu}\\color{red}{\\mathbf L^{-\\top}} [\\color{red}{\\mathbf L^{-1}}\\mathbf (K_{uu} + \\mathbf K_{uf}K_{fu}\\sigma^{-2})\\color{red}{\\mathbf L^{-\\top}}]^{-1}\\color{red}{ \\mathbf L^{-1}}\\mathbf K_{uf}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\phantom{[\\mathbf Q_{ff} + \\sigma^2 \\mathbf I ]^{-1} }= \\sigma^{-2} \\mathbf I - \\sigma^{-4} \\mathbf K_{fu}\\color{red}{\\mathbf L^{-\\top}} [\\mathbf I + \\color{red}{\\mathbf L^{-1}}\\mathbf (\\mathbf K_{uf}K_{fu})\\color{red}{\\mathbf L^{-\\top}}\\sigma^{-2}]^{-1}\\color{red}{ \\mathbf L^{-1}}\\mathbf K_{uf}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For notational convenience, we'll define $\\mathbf L^{-1}\\mathbf K_{uf}\\sigma^{-1} \\triangleq \\mathbf A$, and  $[\\mathbf I + \\mathbf A\\mathbf A^\\top]\\triangleq \\mathbf B$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\phantom{[\\mathbf Q_{ff} + \\sigma^2 \\mathbf I ]^{-1} }= \\sigma^{-2} \\mathbf I - \\sigma^{-2} \\mathbf A^{\\top} [\\mathbf I + \\mathbf A\\mathbf A^\\top]^{-1}\\mathbf A\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\phantom{[\\mathbf Q_{ff} + \\sigma^2 \\mathbf I ]^{-1}}= \\sigma^{-2} \\mathbf I - \\sigma^{-2} \\mathbf A^{\\top} \\mathbf B^{-1}\\mathbf A\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also apply the [matrix determinant lemma](https://en.wikipedia.org/wiki/Matrix_determinant_lemma) to the same:\n",
    "$$\n",
    "|{\\mathbf Q_{ff}} + \\sigma^2 {\\mathbf I}| = |{\\mathbf K_{uu}} + \n",
    " \\mathbf K_{uf}\\mathbf K_{fu}\\sigma^{-2}| \\, |\\mathbf K_{uu}^{-1}| \\, |\\sigma^{2}\\mathbf I|\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Substituting $\\mathbf K_{uu} = {\\mathbf {L L}^\\top}$:\n",
    "$$\n",
    "|{\\mathbf Q_{ff}} + \\sigma^2 {\\mathbf I}| = |{\\mathbf {L L}^\\top} + \n",
    " \\mathbf K_{uf}\\mathbf K_{fu}\\sigma^{-2}| \\, |\\mathbf L^{-\\top}|\\,| \\mathbf L^{-1}| \\, |\\sigma^{2}\\mathbf I|\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "|{\\mathbf Q_{ff}} + \\sigma^2 {\\mathbf I}| = |\\mathbf I + \n",
    " \\mathbf L^{-1}\\mathbf K_{uf}\\mathbf K_{fu} \\mathbf L^{-\\top}\\sigma^{-2}| \\, |\\sigma^{2}\\mathbf I|\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "|{\\mathbf Q_{ff}} + \\sigma^2 {\\mathbf I}| = |\\mathbf B| \\, |\\sigma^{2}\\mathbf I|\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With these two definitions, we're ready to expand the bound:\n",
    "$$ \\mathcal L = \\log \\mathcal N(\\mathbf y\\,|\\,\\mathbf 0,\\, \\mathbf Q_{ff} + \\sigma^2 \\mathbf I) - \\tfrac{1}{2} \\sigma^{-2}\\textrm{tr}(\\mathbf K_{ff} - \\mathbf Q_{ff}) $$\n",
    "\n",
    "$$ = -\\tfrac{N}{2}\\log{2\\pi} -\\tfrac{1}{2}\\log|\\mathbf Q_{ff}+\\sigma^2\\mathbf I| - \\tfrac{1}{2}\\mathbf y^\\top [ \\mathbf Q_{ff} + \\sigma^2 \\mathbf I]^{-1}\\mathbf y - \\tfrac{1}{2} \\sigma^{-2}\\textrm{tr}(\\mathbf K_{ff} - \\mathbf Q_{ff}) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "$$ = -\\tfrac{N}{2}\\log{2\\pi} -\\tfrac{1}{2}\\log(|\\mathbf B||\\sigma^{2}\\mathbf I|)  - \\tfrac{1}{2}\\sigma^{-2}\\mathbf y^\\top (\\mathbf I - \\sigma^{-2} \\mathbf A^{\\top} \\mathbf B^{-1}\\mathbf A)\\mathbf y - \\tfrac{1}{2} \\sigma^{-2}\\textrm{tr}(\\mathbf K_{ff} - \\mathbf Q_{ff}) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ = -\\tfrac{N}{2}\\log{2\\pi} \n",
    "-\\tfrac{1}{2}\\log|\\mathbf B|\n",
    "-\\tfrac{N}{2}\\log\\sigma^{2}\n",
    "-\\tfrac{1}{2}\\sigma^{-2}\\mathbf y^\\top\\mathbf y\n",
    "+\\tfrac{1}{2}\\sigma^{-2}\\mathbf y^\\top\\mathbf A^{\\top} \\mathbf B^{-1}\\mathbf A\\mathbf y\n",
    "-\\tfrac{1}{2}\\sigma^{-2}\\textrm{tr}(\\mathbf K_{ff})\n",
    "+ \\tfrac{1}{2}\\textrm{tr}(\\mathbf {AA}^\\top) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where we have used $\\sigma^{-2}\\textrm{tr}(\\mathbf Q) = \\textrm{tr}(\\mathbf {AA}^\\top)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we define $\\mathbf c \\triangleq \\mathbf L_{\\mathbf B}^{-1}\\mathbf A\\mathbf y \\sigma^{-1}$, with $\\mathbf {L_BL_B}^\\top = \\mathbf B$, so that:\n",
    "$$\n",
    "\\sigma^{-2}\\mathbf y^\\top\\mathbf A^{\\top} \\mathbf B^{-1}\\mathbf A\\mathbf y = \n",
    "\\mathbf c^\\top \\mathbf c\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prediction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At precition 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 on function values at the inducing points $\\mathbf Z$. Recall:\n",
    "\n",
    "$$\n",
    "q(\\mathbf u) = \\mathcal N(\\mathbf u\\,|\\,  \\mathbf m, \\mathbf \\Lambda)\n",
    "$$\n",
    "\n",
    "with \n",
    "\n",
    "$$\\mathbf \\Lambda = \\mathbf K_{uu}^{-1} + \\mathbf K_{uu}^{-1}\\mathbf K_{uf}\\mathbf K_{fu}\\mathbf K_{uu}^{-1} \\sigma^{-2}\n",
    "$$\n",
    "$$\n",
    "\\mathbf m = \\mathbf \\Lambda^{-1} \\mathbf K_{uu}^{-1}\\mathbf K_{uf}\\mathbf y\\sigma^{-2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make a prediction, we need to integrate:\n",
    "    $$ p(\\mathbf f^\\star) = \\int p(\\mathbf f^\\star \\,|\\, \\mathbf u) q (\\mathbf u) \\textrm d \\mathbf u\n",
    "    $$\n",
    "    with \n",
    "    $$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})\n",
    "    $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The integral results in \n",
    "$$\n",
    "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})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note from our above definitions that we have\n",
    "$$\\mathbf K_{uu}^{-1}\\mathbf \\Lambda^{-1}\\mathbf K_{uu}^{-1} = \n",
    "\\mathbf L^{-\\top}\\mathbf B^{-1}\\mathbf L^{-1}$$\n",
    "\n",
    "and further\n",
    "$$\\mathbf K_{uu}^{-1}\\mathbf m = \\mathbf L^{-\\top}\\mathbf L_{\\mathbf B}^{-\\top}\\mathbf c$$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "substituting:\n",
    "$$\n",
    "p(\\mathbf f^\\star) = \\mathcal N(\\mathbf f^\\star\\,|\\, \\mathbf K_{\\star u}\\mathbf L^{-\\top}\\mathbf L_{\\mathbf B}^{-\\top}\\mathbf c, \\,\\mathbf K_{\\star\\star} - \\mathbf K_{\\star u}\\mathbf L^{-1}(\\mathbf I - \\mathbf B^{-1})\\mathbf L^{-1}\\mathbf K_{u\\star})\n",
    "$$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code in `sgpr.py` implements this equation, with an additional switch depending on whether the full covariance matrix is required."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}