Student T Likelihood for outlier detection

[jduerholt]

Hi,

first of all, sorry for creating so many discussions and issues in the moment. We are currently evaluating which steps to do next in a our BO endeavour ;) I promise, that I will also implement some of the proposed solutions and you will see PRs.

I came across this paper today http://proceedings.mlr.press/v84/martinez-cantin18a/martinez-cantin18a.pdf, which I found very interesting. They are using a GP with an approximate Student T likelihood to detect outliers in the BO loop. The predictive distribution of the approximate student T likelihood which allows for exact inference is shown in equation 8.

So far, I never touched likelihoods in Gpytorch, but from what I saw this would amount to implementing a new likelihood or? Could you give me some hints how to best implement the proposed approach, or do you know if it is already available?

Best,

Johannes

[Balandat]

Thanks for the suggestion, support for outlier-robust inference would be quite useful in many places.

So far, I never touched likelihoods in Gpytorch, but from what I saw this would amount to implementing a new likelihood or?

Indeed. Looking at the Laplace approximation derived in Sec 3.1 of https://papers.nips.cc/paper/2009/file/13fe9d84310e77f13a6d184dbf1232f3-Paper.pdf that was used in the paper, this should be pretty straightforward to implement. I haven't looked into this in detail, but one should be able to implement this pretty straightforwardly by just following the Gaussian Likelihood implementation and modifying the marginal code by replacing the noise_covar term with the W_ii from eq (4) in https://papers.nips.cc/paper/2009/file/13fe9d84310e77f13a6d184dbf1232f3-Paper.pdf

I'd be happy to help review a gpytorch PR.

There is a related gpytorch issue about Student-t processes (which are somewhat different than using a Student-t likelihood): cornellius-gp/gpytorch#1858. cc @qingfeng10 had thought about looking into this in the past as well.

[jduerholt]

Thanks for your suggestions. I had a look at the marginal code and for me it is not yet clear how to calculate w_ii inside. To calculate w_ii, I have to calculate $r_i$ which is defined as $y_i - f_i$. But I do not see how I can access $y_i$ and $f_i$ within the marginal code. Can you give me a hint on this?

[Balandat]

Hmm so I took another look at the paper, and I guess I should have spent more time on it the first time. Things aren't as straightforward as I suggested above.

To compute the MLL for a fixed set of hyperparameters you will need to run the EM algorithm (section 4.1) to find the posterior mode \hat{f} that you can then use in the MLL computation of eq (6) [or (15) for the approximate one]. Then fitting the model (i.e. its hyperparameters) then amounts to optimizing that MLL.

The y_i are the observation data. The f_i IIUC are the posterior mean of the latent variable from eq (7) under the current mode estimate \hat{f}.

Implementing this in gpytorch is possible, I'm sure, but I'm not sure how much work it is; we'd have to take a much more detailed look at this. Sorry about the confusion here :|

[jduerholt]

I did a bit more literature research on this topic, I just note it down here:

• Concerning, the paper from above and the Laplace approximation to the Student-t process, I found the implementation from the orginal author Jarna Vanhatalo. Code: https://github.com/gpstuff-dev/gpstuff (unfortunately matlab) Paper:https://arxiv.org/pdf/1206.5754.pdf
• This paper (https://arxiv.org/abs/2011.11057) heavily cites the student-t models from Vanhatalo but claims to have found a better method. I have to read it in detail ... Code can be found here: https://github.com/syrte/robustgp/tree/master

[jduerholt]

@Balandat you even implemented the Wishart priors from @andrewgordonwilson's paper on student t processes into gpytorch five years ago: https://github.com/cornellius-gp/gpytorch/blob/master/gpytorch/priors/wishart_prior.py 👍

[jduerholt]

But it could also be that I am mixing up the student T process from this paper with t-likelihood papers from Vanhatalo ;)

[gpleiss]

It sounds like you're asking for a student-T observation model (on top of a latent Gaussian process), rather than a student T process model (what is proposed in cornellius-gp/gpytorch#1858). GPyTorch has a Student-T likelihood; see the documentation here: https://docs.gpytorch.ai/en/stable/likelihoods.html#studenttlikelihood

Note that the student T likelihood has to be used in conjunction with a variational GP model since exact inference is intractable.

[jduerholt]

Yep, in the original paper from this thread, they are using a student-t likelihood on top of a latent GP. But they are using there some laplace transformation to to make it analytically tractable. This one is also implemented in Gpy (as linked above). The question would be how to implment this in the best way into Gpytorch.

[jduerholt]

After @gpleiss hint I played a bit around with variational GPs with lik-t likelihood and produced a minimal working example:

# Imports
import gpytorch
import torch
import numpy as np
import matplotlib.pyplot as plt
from botorch.models import SingleTaskGP
from botorch.fit import fit_gpytorch_mll
from botorch.utils import standardize
from gpytorch.mlls import ExactMarginalLogLikelihood

# Create the training data
real_std = 0.1
X = np.linspace(0.0, np.pi * 2, 100)[:, None]
Y = np.sin(X) + np.random.randn(*X.shape) * real_std
Y = Y / Y.max()
Yc = Y.copy()
X = X/X.max()
# Slightly noisy data
Yc[75:80] += 1



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

training_iterations = 100

train_x = torch.from_numpy(X.ravel()).to(dtype=torch.float64)
train_y = torch.from_numpy(Yc.ravel()).to(dtype=torch.float64)



# Our testing script takes in a GPyTorch MLL (objective function) class
# and then trains/tests an approximate GP with it on the supplied dataset

#def train_and_test_approximate_gp(objective_function_cls):
model = ApproximateGPModel(train_x).to(dtype=torch.float64)
#likelihood = gpytorch.likelihoods.GaussianLikelihood()
likelihood = gpytorch.likelihoods.StudentTLikelihood()
objective_function = gpytorch.mlls.VariationalELBO(likelihood, model, num_data=train_y.numel())
optimizer = torch.optim.Adam(list(model.parameters()) + list(likelihood.parameters()), lr=0.1)

# Train
model.train()
likelihood.train()
for _ in range(training_iterations):
    output = model(train_x)
    loss = -objective_function(output, train_y)
    loss.backward()
    optimizer.step()
    optimizer.zero_grad()

# Test
model.eval()
likelihood.eval()

with gpytorch.settings.num_likelihood_samples(128):
    observed_pred = likelihood(model(train_x))
samples = observed_pred.sample()


# fit GP
gp = SingleTaskGP(train_x.reshape(100,1).to(dtype=torch.float64), train_y.reshape(100,1).to(dtype=torch.float64))
mll = ExactMarginalLogLikelihood(gp.likelihood, gp)
fit_gpytorch_mll(mll)

mean = gp.posterior(train_x.reshape(100,1).to(dtype=torch.float64)).mean.reshape(100).detach()


fig, ax = plt.subplots(1, 1, figsize=(5, 3))
ax.plot(train_x, samples.mean(dim=0), label="lik-t GP")
ax.plot(train_x, mean, label="GP")

#ax.plot(train_x, samples)

ax.scatter(train_x, train_y, c='k', marker='.', label="Data")
ax.legend(loc="best")
ax.set(xlabel="x", ylabel="y")

plt.show()

This results in a behavior like this:

I will now test it also in a BO setup as in the original paper and maybe also try to implement the Laplace approximation. I keep you posted ;)