Gpytorch: Different lengthscales in different domains

[ankushaggarwal]

I am new to GPytorch and am trying to regress GP to some data which has fast varying and slow varying parts. As a result a single lengthscale cannot capture both. I was wondering if it is possible to have a kernel with two different lengthscales depending.

[dhruvbalwada]

I am working on a similar problem. Using the option ard_num_dims and has_lengthscale seems to work, at least for the simple case.

Here is some example code:

nsamples = 100
x= np.random.rand(nsamples, 2)
y = np.sin(x[:,0] + 2*x[:,1]) + np.random.randn(nsamples) * np.sqrt(0.04)

train_x = torch.from_numpy(x).float()
train_y = torch.from_numpy(y).float()

class ExactGPModel(gpytorch.models.ExactGP):

    def __init__(self, train_x, train_y, likelihood):
        super(ExactGPModel, self).__init__(train_x, train_y, likelihood)  
        self.mean_module = gpytorch.means.ConstantMean()
        self.covar_module = gpytorch.kernels.ScaleKernel(
                            gpytorch.kernels.RBFKernel(ard_num_dims=2,
                                                      has_lengthscale=True))

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

likelihood = gpytorch.likelihoods.GaussianLikelihood()
model = ExactGPModel(train_x, train_y, likelihood)

# just showing how initial values to hyper parameters can be given
hypers = {
    'likelihood.noise_covar.noise': torch.tensor(.01),
    'covar_module.base_kernel.lengthscale': torch.tensor([1.,1.]),
    'covar_module.outputscale': torch.tensor(1.),
     }

model.initialize(**hypers);  

print(model.covar_module.base_kernel.lengthscale[0])

model.train() 
likelihood.train();

optimizer = torch.optim.Adam([
                {'params': model.parameters()},  
                              ], lr=0.1)

mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)


training_iter = 500 # set as many stels as we want to take 

for i in range(training_iter):

        # Zero gradients from previous iteration
    optimizer.zero_grad()

    # Output from model
    output = model(train_x)

    # Calc loss and backprop gradients
    # loss is defined is -log(p(Y(x)))

    # the way the PyTorch optimizer works is that it looks for a variable called loss,
    # which also has a gradient (calc by calling backward), and tries to reduce it.
    loss = -mll(output, train_y)

    loss.backward()

    if np.mod(i,10)==0:
        print('Iter %d/%d - Loss: %.3f   lengthscale1: %.3f lengthscale2: %.3f  noise: %.3f' % (
            i + 1, training_iter, loss.item(),
            model.covar_module.base_kernel.lengthscale[0][0].item(),
            model.covar_module.base_kernel.lengthscale[0][1].item(),
            model.likelihood.noise.item()
              ))

    optimizer.step()

[dhruvbalwada]

^Oops found a bug in older code, edited above code with working version now.

However, if you are thinking of a multiscale problem - slow trend + fast variation, you could refer to the example that is discussed here: https://george.readthedocs.io/en/latest/tutorials/hyper/ , which includes a long term trend + fast seasonal variation, and also refer to chapter 5 in R&W (http://www.gaussianprocess.org/gpml/chapters/RW5.pdf).

[jacobrgardner]

@dhruvbalwada's suggestion to use ARD is on the money assuming you're just talking about different lengthscales in different dimensions. One very small point: has_lengthscale is true by default for RBFKernel, so specifying it should be unnecessary.

If your goal is to achieve different lengthscales in different regions of the space, this effect can often best be achieved by using a small neural net feature extractor before your GP.

If you have a mixture of seasonal trends, the references suggested are good. You could also try spectral kernel learning, which we have a few options for, which can help automate that process (although the initialization can be a bit tricky).

[ankushaggarwal]

@jacobrgardner You are right. I am looking for different lengthscales in different regions of the space. Is there an example on how to use neural net feature extractor before GP?

[dhruvbalwada]

I am curious to know about the neural network feature extractor too.
One approach I am aware of is to use a moving window and assuming local stationarity (used here - https://royalsocietypublishing.org/doi/10.1098/rspa.2018.0400 )

[jacobrgardner]

If you want to try it out, just have your forward method pass data through a NN first:

def forward(self, x):
    x = self.feature_extractor(x)  # any PyTorch NN
    covar_x = self.covar_module(x)
    ...

The NN parameters will automatically get gradient updates during training, since it's all just PyTorch code.