Variance seems to be off by 1

[ToennisStef]

I wanted to validate the calculations of posterior varaince of a GP conditioned on one datapoint and noticed while validating that the calculated variance of GP with gpytorch is off by 1.
To better explain my problem:
I created a basic squared-exponential-kernel(SE-Kernel) class.
with the squared exponential kernel function beeing:

$k(x,x') = \sigma_{f} exp(-\frac{\parallel x-x'\parallel ^2}{2l^2})$

I set the signal variance $\sigma_f$ , and the lenthscale $l$ to 1

Introducing the notation of the data
$y=f(x) + \epsilon, \quad \epsilon \sim N(0,\sigma_n^2) $

the code for the (SE-Kernel) is very basic:

import numpy as np
class RBFKernel():
    signal_variance = 1
    lengthscale = 1

    def CalcCovar(self, x1, x2):
        return (self.signal_variance**2)*np.exp(-np.linalg.norm((x1-x2))**2/(2*self.lengthscale**2))
 myKernel = RBFKernel()

I instantiated a random datapoint at a random location and set an evaluation range:

#datapoint
x_hat,y_hat = np.random.normal(0,1,2)
#evaluation range
res = 10000
x = np.linspace(-5,5,res)

Then i computed the posterior mean and posterior covariance of the GP following this equation on the evaluation range:
$m_{post}(x) = m(x) + K(x,X)(K(X,X)+ \sigma_{n}^{2}I)^{-1}(y-m(X))$
$k_post(x,x') = k(x,x') - K(x,X)(K(X,X)+\sigma_{n}^{2})^{-1}K(X,x')$

The code for this calculation is:

data_noise = 1

m_post = np.empty(res)
covar_post = np.empty(res)
for i in range(0,x.size):
    m_post[i] =  0 + myKernel.CalcCovar(x[i],x_hat)*((myKernel.CalcCovar(x_hat,x_hat)+data_noise**2)**(-1))*(y_hat-0)
    covar_post[i] = myKernel.CalcCovar(x[i],x[i])-myKernel.CalcCovar(x_hat, x[i])*((myKernel.CalcCovar(x_hat,x_hat)+data_noise**2)**(-1))*myKernel.CalcCovar(x[i], x_hat)

As you can see i set the data noise $\sigma_n$ to 1 as well.

next i wanted to validate my calculation of covar_post with the same gaussian process calculated with gpytorch
As I am also using BoTorch for y Thesis i used a custom GP for use in BoTorch, but the implementation is also super basic so here is the code i used for generating the posterior variance:

import gpytorch
import torch
import botorch

class ExactGP(gpytorch.models.ExactGP, botorch.models.gpytorch.GPyTorchModel):

    _num_outputs = 1  # to inform GPyTorchModel API

    def __init__(self, train_X, train_Y, Likelihood): 
        super().__init__(train_X, train_Y, Likelihood)
        self.mean_module = gpytorch.means.ConstantMean()
        self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())

    def forward(self, x):
        mean_x = self.mean_module(x)
        covar_x = self.covar_module(x)
        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)

X_train = torch.tensor([x_hat])
Y_train = torch.tensor([y_hat])

#Model Initialization
likelihood = gpytorch.likelihoods.GaussianLikelihood()
model = ExactGP(X_train, Y_train, likelihood)

model.likelihood.noise_covar.noise = 1
model.covar_module.outputscale = 1
model.covar_module.base_kernel.lengthscale = 1
model.mean_module.constant = 0

test_x = torch.linspace(-5,5,res)
observed_pred = model.likelihood(model(test_x))

The posterior variance can be extracted with observed_pred.variance

But i noticed that the posterior variance calculated with the gpytorch is exactly the varance that i calculated but offset by 1.
When plotting the variance on the evaluation range i get the follwing results:

The code i used for plotting is:

from matplotlib import pyplot as plt

with torch.no_grad():
    f, ax = plt.subplots(1, 1, figsize=(10, 10))
    ax.plot(x,covar_post)
    ax.plot(test_x.numpy(),observed_pred.variance)

Can somebody make sense of this? i can't find a mistake in my implementation of the calculation of the posterior variance. Maybe there is sense to this offset. Maybe i am just to stupid to understand the problem.
I am really confused by this offset and i am struggleing to undertand why this offset exists.
I hope somebody can help me with this.

[Balandat]

The issue is the very last line here:

observed_pred = model.likelihood(model(test_x))

What this does is add the observation noise variance to the posterior. So you're computing the posterior predictive over y, not the posterior over f. model(test_x) will return the posterior mvn that you're looking for. Calling likelihood() on that will add the noise variance to the diagonal, which is why you're getting the offset. If you want the posterior over f just use model(test_x).

[ToennisStef]

Thank you very much! I appologize for my seemingly stupid questions. I still have to learn the basic usage of GPyTorch and BoTorch but i am sometimes struggeling to find the correct reference in the API Reference.

I still don't quite get the interpretation of the posterior predictive over y.
What is the use of:
likelihood(model(test_x)) # Returns the (approximate) predictive posterior distribution at test_x
then?

Is it correct that if i now wanted to sample another datapoint from my exact function f at new_x, the posterior predictive distribution that i got from model(new_x) represents the value range where i expect the new sample value?
If that is correct what again is the use of likelihood(model(new_x))?

[Balandat]

model(new_x) returns the distribution over f, and likelihood(model(new_x)) returns the distribution over y. The latter just adds the observation noise sigma_n^2. So if you want to reason about the next realized function value is you'd look at model(new_x), if you want to reason about what observations y you can expect you'd look at likelihood(model(new_x)).

[ToennisStef]

Ok thank very much this really helps my understanding!