Computing posterior variance at noisy training samples using Gaussian Process regression

Question

Sorry for a possible naive question but this has been unclear to me for awhile...

Consider the data model $y = f(x) + \epsilon, \;\; \epsilon \sim \mathcal{N}(0, \sigma_v^2)$. They show in Rasmussen and Williams (Eq. 2.26) for noisy observational data, the posterior variance of a Gaussian Process at a new sample $x_*$ is equal to,

$\sigma^2(x_*) = K_{ss} - K_s^T (K + \sigma_v^2I)^{-1} K_{s} \qquad \qquad (1)$

where $K$ is the covariance matrix at the training samples, $K_s$ is the covariance matrix of new test samples with the training samples, $K_{ss}$ is the covariance matrix of the new test samples, and $\sigma_v^2$ is the observation noise.

Several lecture slides use a different equation where an additional $\sigma^2_v$ is added,

$\sigma^2(x_*) = K_{ss} + \sigma_v^2 I - K_{s}^T (K + \sigma_v^2I)^{-1} K_{s}\qquad \qquad (2)$

If I wanted to compute the variance at the training locations and test locations, which equation would give the correct variance?

Is Eq. 1 just the variance associated with the GP model, and neglects the variance from the inferred noise model $\sigma^2_v$? Then, when making a prediction I would add the variances, e.g., $\sigma^2(x) = \sigma^2_{gp}(x) + \sigma^2_v(x)$ where $\sigma^2_{gp}$ is equal to equation 1?

Answer 1

$\newcommand{\f}{\mathbf f}\newcommand{\e}{\varepsilon}\newcommand{\y}{\mathbf y}\newcommand{\0}{\mathbf 0}$If $f \sim \mathcal {GP}(0, k)$ and we have sample points $x_1,\dots,x_n$ then letting $\f = (f(x_1), \dots, f(x_n))^T$ we have $$ \f \sim \mathcal N(\0, K) $$ with $K_{ij} = k(x_i,x_j)$ as the kernel matrix. If we observe $$ y_i = f(x_i) + \e_i $$ so $$ \y = \f + \e $$ where $\e \sim \mathcal N(\0, \sigma^2_\e I) \perp f$ then $$ \y \sim \mathcal N(\0, K + \sigma^2_\e I). $$ Now for a new point $x_s$ we have two options. If we are just interested in the posterior distribution of a point of the Gaussian process itself, i.e. $f_s := f(x_s)$, we'll have $$ {\f \choose f_s} \sim \mathcal N\left({\0 \choose 0}, \begin{bmatrix}K & k_s^T \\ k_s & k_{ss}\end{bmatrix}\right) $$ so $$ {\y \choose f_s} = {\f \choose f_s} + {\e \choose 0} \sim \mathcal N\left({\0 \choose 0}, \begin{bmatrix}K + \sigma^2_\e I & k_s^T \\ k_s & k_{ss}\end{bmatrix} \right) $$ which means the posterior variance of $f_s\mid \y$ is $$ \text{Var}(f_s\mid \y) = k_{ss} - k_s^T(K + \sigma^2_\e I)^{-1}k_s. $$

If however we are interested in the posterior distribution of a new noisy data point $y_s = f_s + \e_s$ then we'll have $$ {\y \choose y_s} = {\f \choose f_s} + {\e \choose \e_s} \sim \mathcal N\left({\0 \choose 0}, \begin{bmatrix}K & k_s^T \\ k_s & k_{ss}\end{bmatrix} + \sigma^2_\e I_{n+1}\right). $$ so posterior variance of $y_s\mid \y$ is $$ \text{Var}(y_s\mid \y) = (k_{ss} + \sigma^2_s) - k_s^T(K + \sigma^2_\e I)^{-1}k_s. $$

Thus it depends on whether or not we're talking about the posterior variance of the GP itself versus new data points coming from the noisy process. If you want to get a sense of the variation of the GP $f$ then you can look at $f_s \mid \y$ but if you're interested in what noisy data will look like coming from it then $y_s\mid \y$ will be more interesting.

In linear regression this is analogous to looking at the estimated mean for a new point $\text{Var}(x_s^T\hat\beta)$ versus looking at the variance of a new data point noise and all $\text{Var}(x_s^T\hat\beta + \e_s)$. Both can be reasonable but it depends what we're trying to do. Often though in my experience we care more about the mean, i.e. $f$ or $\hat\beta$, so usually that's what I go with, rather than the predictive distribution for a new data point.

Here's an example showing the posterior mean with 1.96 times the standard deviation of the variance for just the GP itself vs for new data. The underlying GP is well captured by the red line while the data is well described by the blue line, which just has the extra $\sigma^2_\e$ term.