cornellius-gp / gpytorch

A highly efficient implementation of Gaussian Processes in PyTorch
MIT License
3.59k stars 562 forks source link

[Question] Constraints for inducing points in SGPR #1675

Closed DanielH1994 closed 3 years ago

DanielH1994 commented 3 years ago

Hi,

I am using a simple example to implement the SGPR model. I find that the selected inducing points located sometimes outside the input domain. It is unacceptable in my case. In fact, in my application, the point has physical meanings and can not be negative.

I find that, by registering the raw parameters, those raw parameters raw_inducing_points will have the same value as the original parameter inducing_point but some of them will still locate outside the constraint interval. Could you please help me for how to set the constraint for the inducing_points? Here are the result figure I get and the code snippet.

SGPt

import math, torch, gpytorch, time, os
from matplotlib import pyplot as plt

def testfunction(x):
    return torch.sin((x+1) * x * math.pi) / (x + 0.5)

train_x = torch.linspace(-0.1, 1.1, 121)
train_y = testfunction(train_x) 
inducingnum = 20
learningrate = 0.05
training_iterations = 150

ind_ini = train_x[torch.randint(low=11, high=train_x.size()[0]-11, size=(inducingnum,))]

class  SGPRModel(gpytorch.models.ExactGP):
    def __init__(self, train_x, train_y, likelihood) -> None:
        super(SGPRModel, self).__init__(train_x, train_y, likelihood)
        self.mean_module = gpytorch.means.ConstantMean()
        self.base_covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.MaternKernel(nu=1.5))
        self.covar_module = gpytorch.kernels.InducingPointKernel(
            base_kernel=self.base_covar_module, likelihood=likelihood,
            inducing_points=ind_ini
        )
    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(noise_constraint=gpytorch.constraints.Positive())
model = SGPRModel(train_x, train_y, likelihood)

lower_bound, upper_bound = torch.zeros_like(ind_ini).reshape(-1,1), torch.ones_like(ind_ini).reshape(-1,1)
constrained_tensor = gpytorch.constraints.Interval(lower_bound=lower_bound, upper_bound=upper_bound)
model.covar_module.register_parameter(name="raw_inducing_points", parameter=torch.nn.Parameter(ind_ini.reshape(-1,1)))
model.covar_module.register_constraint("raw_inducing_points", constrained_tensor)

for param_name, param in model.named_parameters():
    print(f'Parameter name: {param_name:42} value = {param}')
for constraint_name, constraint in model.named_constraints():
    print(f'Constraint name: {constraint_name:55} constraint = {constraint}')

model.train()
likelihood.train()

optimizer = torch.optim.Adam(model.parameters(), lr=learningrate)

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

def train():
    loss_pre = math.inf
    for i in range(training_iterations):
        optimizer.zero_grad()
        output = model(train_x)
        loss = -mll(output, train_y)
        loss.backward()
        if i%15 == 14:
            print('Iter %d/%d - Loss: %.5f lengthscale: %.5f  noise: %.5f' % (
                i+1, training_iterations, loss.item(), 
                model.covar_module.base_kernel.base_kernel.lengthscale.item(),
                model.likelihood.noise.item()
                ))
        loss_pre = loss.item()
        optimizer.step()

start = time.time()
train()
end = time.time() - start
print('The training time: %.5f s' % end)

model.initialize(**{'likelihood.noise_covar.noise': torch.tensor(1.e-8),})

print(model.covar_module.raw_inducing_points)
print(model.covar_module.inducing_points)
print(constrained_tensor.transform(model.covar_module.raw_inducing_points))

model.eval()
likelihood.eval()

with torch.no_grad(), gpytorch.settings.fast_pred_var():
    test_x = torch.linspace(0, 1, 21)
    true_y = testfunction(test_x)
    observed_pred = likelihood(model(test_x)) 
    inducingpoints = model.covar_module.inducing_points
    inducinglabels = model(inducingpoints).mean
with torch.no_grad():
    f, ax = plt.subplots(1, 1, figsize=(8, 4))
    lower, upper = observed_pred.confidence_region()
    ax.plot(inducingpoints.detach().numpy(), inducinglabels.detach().numpy(), 'ro')
    ax.plot(test_x.numpy(), true_y.numpy(), 'r')
    ax.plot(test_x.numpy(), observed_pred.mean.numpy(), 'b')
    ax.fill_between(test_x.numpy(), lower.numpy(), upper.numpy(), alpha=0.5)
    ax.set_ylim([-1.0, 1.5])
    ax.legend(['Inducing Points', 'True', 'Mean', 'Confidence'])
    f.savefig('SGPt.png') 
wjmaddox commented 3 years ago

Hi, assigning the constraint over the inducing points on the module level is the correct thing to do here. Here's a working example:

from gpytorch.constraints import Interval
from gpytorch.kernels import InducingPointKernel

class ConstrainedInducingPointKernel(InducingPointKernel):
    def __init__(self, base_kernel, inducing_points, likelihood, constraint, active_dims=None):
        super().__init__(
            base_kernel=base_kernel, 
            inducing_points=inducing_points, 
            likelihood=likelihood,
            active_dims=None,
        )
        # define the raw inducing points in the kernel initialization method
        del self.inducing_points # not sure if this is strictly necessary but probably good for book-keeping
        self.register_parameter(name="raw_inducing_points", parameter=torch.nn.Parameter(inducing_points))
        self.register_constraint("raw_inducing_points", constraint)

    @property
    def inducing_points(self):
        return self.raw_inducing_points_constraint.transform(
            self.raw_inducing_points
        )

# now the model looks something like this
class  SGPRModel(gpytorch.models.ExactGP):
    def __init__(self, train_x, train_y, likelihood) -> None:
        super(SGPRModel, self).__init__(train_x, train_y, likelihood)
        self.mean_module = gpytorch.means.ConstantMean()
        self.base_covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.MaternKernel(nu=1.5))

        # only change is we swapped out the kernel class and put in a constraint over the inducing pts
        self.covar_module = ConstrainedInducingPointKernel(
            base_kernel=self.base_covar_module, likelihood=likelihood,
            inducing_points=ind_ini, constraint=Interval(0.0, 2.0),
        )

    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(noise_constraint=gpytorch.constraints.Positive())
model = SGPRModel(train_x, train_y, likelihood)

With this model and the rest of your code, I'm able to get something like this: image

DanielH1994 commented 3 years ago

It works very well! Thank you very much!