SVGP, SGPR - Inducing Points locations

[mxef]

Hello, I am trying to fit a SVGP model to time-series data but I get some inducing points outside the range of the training data on the x-axis. The same thing happens with the SGPR model. Is this expected behavior?

[gpleiss]

How did you initialize/optimize the inducing point locations? Did you follow the example in the SGPR/SVGP notebook?

The ELBO, which is used to optimize inducing point locations, is non-convex, and so inducing point locations are highly dependent on initialization. Try initializing the inducing points to be a random subset of data.

[mxef]

Ok, now I used four different initializations for the SGPR model (to optimize only the locations of the inducing points for more reliable results), which are:
inducing_points=torch.rand(inducing_p)
inducing_points=torch.FloatTensor(inducing_p).uniform_(torch.min(train_x), torch.max(train_x)
inducing_points=torch.linspace(torch.min(train_x), torch.max(train_x), inducing_p)
inducing_points = train_x[np.random.randint(0,train_x.shape, size=inducing_p)]
#inducing_p: number of inducing points
The last three of them behaved pretty much the same. Sometimes I get a few inducing points completely outside the training data, like the plot here (with the last initialization)

Is this only related to the non-convexity of the ELBO or is it also an implementation issue?

[gpleiss]

Can you post a code example?

[mxef]

The code below prints 5 plots of the UCR dataset 'Beef' with 5*log2(length+1) inducing points. I ran it a lot of times and I always get such irrelevant inducing points in 2 or 3 of the 5 plots each time.

import numpy as np
import matplotlib.pyplot as plt
from tslearn.datasets import UCR_UEA_datasets
import math
import torch
import gpytorch
from gpytorch.means import ConstantMean
from gpytorch.kernels import ScaleKernel, RBFKernel, MaternKernel, InducingPointKernel
from gpytorch.distributions import MultivariateNormal
from gpytorch.means import ConstantMean
from gpytorch.kernels import ScaleKernel, RBFKernel, InducingPointKernel
from gpytorch.distributions import MultivariateNormal
import random

inducing_p = 0

class GPRegressionModel(gpytorch.models.ExactGP, ):
    def __init__(self, train_x, train_y, likelihood):
        super(GPRegressionModel, self).__init__(train_x, train_y, likelihood)
        global inducing_p
        self.mean_module = ConstantMean()
        self.base_covar_module = MaternKernel(0.5)
        self.covar_module = InducingPointKernel(self.base_covar_module, 
                                                inducing_points = train_x[np.random.randint(0,train_x.shape, size=inducing_p)], 
                                                likelihood=likelihood)
    def forward(self, x):
        mean_x = self.mean_module(x)
        covar_x = self.covar_module(x)
        return MultivariateNormal(mean_x, covar_x)


def run_set(dataset, ind_points):
    global inducing_p
    X_train, y_train, X_test, y_test = UCR_UEA_datasets().load_dataset(dataset)

    train_y = X_train[0,:,0]
    for i in range(1,20):
        train_y = np.vstack([train_y, X_train[i,:,0]])

    train_x = torch.from_numpy(np.arange(0,train_y.shape[1])).float()
    train_y = torch.from_numpy(train_y).float()
    train_x = (train_x - train_x.mean(0)) / train_x.std(0)

    #Normalization [-1,1]
    for i in range(20):
        train_y[i] = train_y[i] - torch.min(train_y[i])
        train_y[i] = 2 * (train_y[i] / torch.max(train_y[i])) - 1 

    inducing_p = ind_points*(int(math.log2(train_y.shape[1]))+1)

    # initialize likelihood and model
    likelihood = []
    model = []
    for i in range(20):
        likelihood.append(gpytorch.likelihoods.GaussianLikelihood())
        model.append(GPRegressionModel(train_x, train_y[i], likelihood[i]))

    training_iter = 50
    for i in range(20):
        # Find optimal model hyperparameters
        model[i].train()
        likelihood[i].train()

        # Use the adam optimizer
        optimizer = torch.optim.Adam(model[i].parameters(), lr=0.1)  # Includes GaussianLikelihood parameters
        # "Loss" for GPs - the marginal log likelihood
        mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood[i], model[i])

        for it in range(training_iter):
            # Zero gradients from previous iteration
            optimizer.zero_grad()
            # Output from model
            output = model[i](train_x)
            # Calc loss and backprop gradients
            loss = -mll(output, train_y[i])
            loss.backward()        
            optimizer.step()
            

    # Get into evaluation (predictive posterior) mode
    for i in range(20):
        model[i].eval()
        likelihood[i].eval()
        
    preds = []
    with gpytorch.settings.max_preconditioner_size(10), torch.no_grad():
        for i in range(20):
            preds.append(model[i](train_x + 0.0001))

    # Test points are regularly spaced along [0,1]
    # Make predictions by feeding model through likelihood
    print("{:s}:\n".format(dataset))
    for i in range(5):
        with torch.no_grad():
            # Initialize plot
            f, ax = plt.subplots(1, 1, figsize=(12, 9))
            # Get upper and lower confidence bounds
            lower, upper = preds[i].confidence_region()
            # Plot training data as black stars
            ax.plot(train_x.numpy(), train_y[i].numpy(), 'k*')            
            ax.plot(model[i].covar_module.inducing_points.detach().numpy()[:,-1], 
                       model[i](model[i].covar_module.inducing_points).mean.detach().numpy(), 'ro')

            # Plot predictive means as blue line
            ax.plot(train_x.numpy(), preds[i].mean.numpy(), 'b')
            # Shade between the lower and upper confidence bounds
            ax.fill_between(train_x.numpy(), lower.numpy(), upper.numpy(), alpha=0.5)
            ax.set_ylim([-3, 3])
            ax.legend(['Observed Data', 'Inducing Points', 'Mean', 'Confidence'])            
            plt.show()


run_set('Beef', 5)