SVGP-Poor performance on testing data

[mbelalsh]

Hello there. I am following an SVGP example from GPyTorch docs: SVGP

This code works fine on elevators dataset, but it is not working on my problem as it performs very poorly on testing data. I am applying SVGP on an inverted pendulum problem wherein the input features and output are given below

• Input features: angle, angular velocity, velocity, force. X_train is (8000, 4)
• Output: angular acceleration. Y_train_theta is (8000, )

I am using a single-output SVGP and following the exact approach you described in the docs. I have also normalized the dataset.
Right now, I am playing with batch_size and num_inducing_points to get an improved performance on testing dataset. Also, in SVGP example, you are passing the first 500 features (X_train[:500, :]) as inducing points to the model, whereas in some other examples, you are randomly selecting inducing points from the training dataset. Any difference?

Please suggest me some tips to get good performance on testing data. I have tried a lot. FYI, simple GP works fine on this example.

I have some other related questions as well. It seems like the smaller number of inducing points (like 500) gives better results.

• How many inducing points to use? Thinking of 500 as a benchmark, when I increase the num_inducting_points, I don't necessarily get better results. If I reduce the num_inducting_points, it throws an error like matrix is not positive definite after adding jitter.
• How long do we need to train the model? I trained the model for 27k epochs with a batch size of 1024 (which takes a long time), but the loss value saturates after 3k epochs and doesn't improve. The lowest loss value I could get with SVGP was 0.675. This seems unusual to me. I tuned batch_size and num_inducing_points a lot many times.
• What is the expected lowest loss value for SVGP? I know that a simple GP can give us as low as -1.956 with pretty good results.
  The code is given below. This is a complete working example:

import tqdm
import torch
import gpytorch
import numpy as np
from matplotlib import pyplot as plt
import time

def cartPoleDyn(theta, theta_d, x, x_d, m1, m2, length, k_inertia, mu_c, mu_p, gravity, force, dt):
    theta_dd = ((m1+m2)*gravity*np.sin(theta)-np.cos(theta)*(force+m2*length*(theta_d**2)*np.sin(theta)-mu_c*x_d)-((m1+m2)*mu_p*theta_d)/(m2*length)) / ((m1+m2)*(k_inertia+1)*length - m2*length*(np.cos(theta))**2)
    x_dd = (force-m2*length*(theta_dd*np.cos(theta)-(theta_d**2)*np.sin(theta))-mu_c*x_d) / (m1+m2)
    # Euler Integrator 
    x_d = x_d + x_dd*dt
    x = x + x_d*dt
    theta_d = theta_d + theta_dd*dt
    theta = theta + theta_d*dt
    return theta_dd, theta_d, theta, x_dd, x_d,  x

m1 = 1 
m2 = 0.1
length = 0.5
k_inertia = 1/3
mu_c = 0.0005
mu_p = 0.000002
gravity = 9.8 
dt = 0.01

# Initialization 
theta = 0.3 # np.random.uniform(low=-0.05, high=0.05)
theta_d = 1 # np.random.uniform(low=-0.05, high=0.05)
x = 0 # np.random.uniform(low=-0.05, high=0.05)
x_d = 0 # np.random.uniform(low=-0.05, high=0.05)

# lists for storing results 
theta_dd_sim = []           # theta double derivative 
theta_d_sim = [theta_d]     # theta derivative 
theta_sim = [theta]         # theta  
x_dd_sim = []               # x double derivative 
x_d_sim = [x_d]             # x derivative 
x_sim = [x]                 # x

noOfSteps = 12000
ts = np.linspace(0, noOfSteps/(1/dt), num=noOfSteps, endpoint=False)

for i in range(noOfSteps):
    force = -0.2+0.5*np.sin(6*ts[i])
    #force = -1+np.sin(ts[i])+3*np.sin(2*ts[i])
    #force = -0.5+0.2*np.sin(ts[i])+0.3*np.sin(2*ts[i])
    theta_dd, theta_d, theta, x_dd, x_d, x = cartPoleDyn(theta, theta_d, x, x_d, m1, m2, length, k_inertia, mu_c, mu_p, gravity, force, dt)
    theta_dd_sim.append(theta_dd)
    theta_d_sim.append(theta_d)
    theta_sim.append(theta)
    x_dd_sim.append(x_dd)
    x_d_sim.append(x_d)
    x_sim.append(x)

# Converting to Numpy Arrays for Plotting
theta_dd_sim = np.array(theta_dd_sim)
theta_d_sim = np.array(theta_d_sim)
theta_sim = np.array(theta_sim)
x_dd_sim = np.array(x_dd_sim)
x_d_sim = np.array(x_d_sim)
x_sim = np.array(x_sim)

# Normalizing Training and Testing data
st_Traj = 0
en_Traj = 8000 # 2307 # 9405

timeSteps = np.linspace(0, noOfSteps/(1/dt), num=noOfSteps, endpoint=False)
forceIn = -0.2+0.5*np.sin(6*timeSteps) # control input function 

# Training Input (Standardized)
theta_sim_mean = theta_sim[st_Traj:en_Traj].mean()
theta_sim_stdev = theta_sim[st_Traj:en_Traj].std()

theta_d_sim_mean = theta_d_sim[st_Traj:en_Traj].mean()
theta_d_sim_stdev = theta_d_sim[st_Traj:en_Traj].std()

x_d_sim_mean = x_d_sim[st_Traj:en_Traj].mean()
x_d_sim_stdev = x_d_sim[st_Traj:en_Traj].std()

forceIn_mean = forceIn[st_Traj:en_Traj].mean()
forceIn_stdev = forceIn[st_Traj:en_Traj].std()  

# Training Output (Standardized)
theta_dd_sim_mean = theta_dd_sim[st_Traj:en_Traj].mean()
theta_dd_sim_stdev = theta_dd_sim[st_Traj:en_Traj].std()  

x_dd_sim_mean = x_dd_sim[st_Traj:en_Traj].mean()
x_dd_sim_stdev = x_dd_sim[st_Traj:en_Traj].std() 

X_train = np.stack(((theta_sim[st_Traj:en_Traj]-theta_sim_mean)/theta_sim_stdev,(theta_d_sim[st_Traj:en_Traj]-theta_d_sim_mean)/theta_d_sim_stdev,(x_d_sim[st_Traj:en_Traj]-x_d_sim_mean)/x_d_sim_stdev,(forceIn[st_Traj:en_Traj]-forceIn_mean)/forceIn_stdev), axis=1)

Y_train_theta = ((theta_dd_sim[st_Traj:en_Traj]-theta_dd_sim_mean)/theta_dd_sim_stdev)  # RMSE(A)
Y_train_x = ((x_dd_sim[st_Traj:en_Traj]-x_dd_sim_mean)/x_dd_sim_stdev)  # RMSE(V)

# Testing Output (Standardized)

Y_test_theta = ((theta_dd_sim[en_Traj:]-theta_dd_sim_mean)/theta_dd_sim_stdev)  # RMSE(A)
Y_test_x = ((x_dd_sim[en_Traj:]-x_dd_sim_mean)/x_dd_sim_stdev)  # RMSE(V)

Y_theta = ((theta_dd_sim-theta_dd_sim_mean)/theta_dd_sim_stdev)
Y_x = ((x_dd_sim-x_dd_sim_mean)/x_dd_sim_stdev)
X_prior = np.stack(((theta_sim[:-1]-theta_sim_mean)/theta_sim_stdev,(theta_d_sim[:-1]-theta_d_sim_mean)/theta_d_sim_stdev,(x_d_sim[:-1]-x_d_sim_mean)/x_d_sim_stdev,(forceIn-forceIn_mean)/forceIn_stdev), axis=1) # Prior

# Numpy to tensors 
X_train = torch.from_numpy(X_train).float()
Y_train_theta = torch.from_numpy(Y_train_theta).float()
Y_train_x = torch.from_numpy(Y_train_x).float()
X_prior = torch.from_numpy(X_prior).float()
Y_theta = torch.from_numpy(Y_theta).float()
Y_x = torch.from_numpy(Y_x).float()

# Putting tensors on GPU
if torch.cuda.is_available():
    X_train, Y_train_theta, Y_train_x, X_prior, Y_theta, Y_x = X_train.cuda(), Y_train_theta.cuda(), Y_train_x.cuda(), X_prior.cuda(), Y_theta.cuda(), Y_x.cuda()

train_dataset = TensorDataset(X_train, Y_train_theta)
train_loader = DataLoader(train_dataset, batch_size=1024, shuffle=True)

test_dataset = TensorDataset(X_prior, Y_theta)
test_loader = DataLoader(test_dataset, batch_size=1024, shuffle=False)

class StandardApproximateGP(gpytorch.models.ApproximateGP):
    def __init__(self, inducing_points):
        variational_distribution = gpytorch.variational.CholeskyVariationalDistribution(inducing_points.size(0))
        variational_strategy = gpytorch.variational.VariationalStrategy(self, inducing_points, variational_distribution, learn_inducing_locations=True)
        super().__init__(variational_strategy)
        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)

indPts = np.sort(np.random.randint(low=0, high=8000, size=(500,)))
inducing_points = X_train[indPts,:]
#inducing_points = X_train[:500, :]
#inducing_points = torch.randn(4000, X_train.size(-1), dtype=X_train.dtype, device=X_train.device)
model = StandardApproximateGP(inducing_points=inducing_points)
likelihood = gpytorch.likelihoods.GaussianLikelihood()

if torch.cuda.is_available():
    model = model.cuda()
    likelihood = likelihood.cuda()

model.train()
likelihood.train()

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

# Our loss object. We're using the VariationalELBO
mll = gpytorch.mlls.VariationalELBO(likelihood, model, num_data=Y_train_theta.size(0))

start_time = time.time()
losses = []
num_epochs = 18000
epochs_iter = tqdm.notebook.tqdm(range(num_epochs), desc="Epoch")
for i in epochs_iter:
    # Within each iteration, we will go over each minibatch of data
    minibatch_iter = tqdm.notebook.tqdm(train_loader, desc="Minibatch", leave=False, disable=True)
    for x_batch, y_batch in minibatch_iter:
        optimizer.zero_grad()
        output = model(x_batch)
        loss = -mll(output, y_batch)
        minibatch_iter.set_postfix(loss=loss.item())
        losses.append(loss.detach().cpu().numpy())
        loss.backward()
        optimizer.step()
    if i%1000==0:
        print('Loss: %.3f' % (loss.item()))    
 
print('Total time is %s:' %(time.time() - start_time)) 

model.eval()
likelihood.eval() 

means = np.zeros(X_prior.shape[0])
lowers = np.zeros(X_prior.shape[0])
uppers = np.zeros(X_prior.shape[0])

with torch.no_grad():
    i = 0
    #for x_batch, y_batch in test_loader:
    preds = model(X_prior)
    lower, upper = preds.confidence_region()
    means = preds.mean.cpu().numpy()
    lowers = lower.cpu().numpy()
    uppers = upper.cpu().numpy()
    i += 1  

start_traj = 0
end_traj = 12000

fig, ax = plt.subplots(1)

fig.set_figheight(5)
fig.set_figwidth(16)

ax.plot(timeSteps[start_traj:end_traj], Y_theta[start_traj:end_traj].cpu().numpy(), 'k')
ax.plot(timeSteps[start_traj:end_traj], means[start_traj:end_traj], 'b')
ax.fill_between(timeSteps[start_traj:end_traj], lowers[start_traj:end_traj], uppers[start_traj:end_traj], alpha=0.3)
ax.legend(['Observed Data', 'Mean', 'Confidence'])
#plt.savefig('withInducPts(500).png')

You can see the performance of the model below. It performs well on training data (i.e., first 8k points), then it is messed up after that which is testing data.

[gpleiss]

The function that you are modeling has low intrinsic lengthscale (it changes rapidly). This is known to be a setting that is notoriously difficult for SVGP - you would need nearly N inducing points (where N is the amount of training data) to accurately model the function everywhere. See for example this paper: https://arxiv.org/pdf/2202.01694.pdf

I am playing with batch_size and num_inducing_points to get an improved performance on testing dataset.

Increasing num_inducing_points should improve performance.

Also, in SVGP example, you are passing the first 500 features (X_train[:500, :]) as inducing points to the model, whereas in some other examples, you are randomly selecting inducing points from the training dataset. Any difference?

This affects optimization. Using k_means to initialize inducing points tends to work the best from my experience, but initializing to the first n training data points also works well.

What is the expected lowest loss value for SVGP?

It is larger than the loss for the exact GP, since the loss for SVGP is an ELBO.

[mbelalsh]

Geoff. Thanks for the insightful answers.

• I admit that my control input function is time-varying, as you already noticed, and it changes continuously after every 0.01s. That being said, I would still like to see the good performance of SVGP on this problem because Deep GP and DSPP are using inducing points-based implementation. So, I want to see SVGP working on this simple toy problem. Also, I want to try Multi-Output SVGP later on and compare its performance with Multi-Output GP in terms of speed and accuracy.
• I just want to say that I also tried with num_inducing_points = N earlier, but it did not work. Also, it tends to be super slow in terms of training compared to simple GP. This defeats the purpose of SVGP because I saw some of your answers wherein you said that SVGP tends to be faster. I do admit that the inference time on testing data is super fast with SVGP compared to simple GP, which I appreciate.
• I am again gonna try SVGP with num_inducing_points = N, and then I will let you know.
• If you say so, I can give Variational Nearest Neighbour GP a try with num_inducing_points = N. I will let you know how it goes.

[mbelalsh]

@gpleiss I was able to get better performance with SVGP. Previously, I did not pass ard_num_dims=4 to gpytorch.kernels.RBFKernel and it turned out that this has caused a major performance improvement. Without passing this as an argument, I was not getting good performance despite having nearly N inducing points. Now, you can see the performance in the attached figure. You can see that the prediction are good on some testing data points and then all of a sudden goes to zero. Any thoughts as to why the predictions would suddenly go to zero?

Predictions are black, and training/testing data are green/red.
Please click on the figure to open it in a new tab for better visualization

[gpleiss]

This is the expected behavior of an RBF GP. The posterior GP reverts to the mean when it is not near observed data.
If you want better extrapolation, I suggest looking at the spectral mixture kernel: https://arxiv.org/pdf/1302.4245.pdf

[mbelalsh]

Thanks, @gpleiss. An amazing answer like always.

I have been looking into kernels recently as to how they affect predictions on a given dataset. They seem to be important as per my experience, but I don't see discussions on selecting kernels in GP-related papers published recently.

• Do you know of any systematic way to select kernels for a given dataset other than heuristic approaches?
• For a given dataset, how do we select a kernel that would give the best performance? Or Do we just follow the trial-and-error approach?

Any suggestions on related papers would mean a lot to me.

Thanks :)

[gpleiss]

Do you know of any systematic way to select kernels for a given dataset other than heuristic approaches?

This question is the fundamental question of Bayesian model selection. Chapter 5 of Rasmussen and Williams provides some ideas for kernel selection. The standard approach is comparing the log marginal likelihoods of GP models with different kernels. However - this is not always the best approach (see https://arxiv.org/abs/2106.06510 and https://arxiv.org/abs/2202.11678 for recent discussions).