Time series reconstruction with ModelList (Multi-Output) GP Regression

[dhenriques62]

I am trying to reconstruct a time series (OMP) from another as a proxy (ERA5) and adapting this example. Supposedly, the model should fill in the missing values at the instants where there are values in the proxy series, but that's not what I get:

Do I have to make any changes in the example code ?

[ekinugurel]

Do not know much about your domain specifics but you usually need to represent time in a variety of ways in order to encode enough structure. Otherwise, you'd just be using GPs to "smooth" the data. If you're representing time as a monotonically-increasing datetime object (i.e., Unix time), the GP will just revert to its mean function in the absence of any nearby data. I recommend one-hot encoding various temporal dimensions like the day-of-week (0-6), hour-of-day (0-23), and whatever else you think would capture variability for your domain. Then, you can use an ARD kernel to fit each input dimension w.r.t. to your dependent variable.

Alternatively, you can try messing around with the PeriodicKernel to see if you can capture any periodic trends in your data, but I recommend doing this only with a monotonically-increasing variable like Unix time. I also don't see a clear periodic trend in your data, so this may not be the best approach. The kernel cookbook (and David Duvenaud's dissertation) is an easy read if you want to better understand how to use and combine different kernels to achieve your goals in time-series analysis.

[dhenriques62]

Thanks @ekinugurel.
I am exploring the potential of GP in climate series. In fact the data are encoded with an increasing time ordinate. These series have basically a daily and an annual variability. There may be others frequencies in between, but of lesser magnitude. I had not thought to order the data with respect to another time variable, but your suggestion makes sense. I will give it a try.

[dhenriques62]

Hi @ekinugurel ,

Unfortunately it didn't work out. I encoded the times in hours (0 to 23) and got this:

In fact, if the inputs vary only between 0 and 23 they should repeat the same outputs.
Is my reasoning correct ?

[gpleiss]

As an alternative, consider also using the spectral mixture kernel (without one-hot encodings). The spectral mixture kernel is good for time series extrapolations (assuming your data are stationary, which may not be the case based on the images you have shared), although you will have to properly initialize the hyperparameters. See https://docs.gpytorch.ai/en/stable/examples/01_Exact_GPs/Spectral_Mixture_GP_Regression.html

If you are not doing a one-hot encoding of your time data, it is helpful to normalize your data (e.g. so that the time variable is between 0 and 1).

[dhenriques62]

Thank you @gpleiss.

I do not intend to extrapolate but rather to interpolate long periods. I did what you suggested and got this:

The data was normalised, and the model setup used was this:

class SpectralMixtureGPModel(gpytorch.models.ExactGP):

def __init__(self, train_x, train_y, likelihood):
    super(SpectralMixtureGPModel, self).__init__(train_x, train_y, likelihood)
    self.mean_module = gpytorch.means.ConstantMean()
    self.covar_module = gpytorch.kernels.SpectralMixtureKernel(num_mixtures=4)
    self.covar_module.initialize_from_data(train_x, train_y)

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

likelihood1 = gpytorch.likelihoods.GaussianLikelihood()
model1 = SpectralMixtureGPModel(train_x1, train_y1, likelihood1)

likelihood2 = gpytorch.likelihoods.GaussianLikelihood()
model2 = SpectralMixtureGPModel(train_x2, train_y2, likelihood2)

However, I had tried this other setup:

class ExactGPModel(gpytorch.models.ExactGP):

def __init__(self, train_x, train_y, likelihood):
  super().__init__(train_x, train_y, likelihood)
  self.mean_module = gpytorch.means.ConstantMean()
  self.mean_module = gpytorch.means.ConstantMean()
  constraint1 = gpytorch.constraints.Interval(1e-1, 1)
  constraint2 = gpytorch.constraints.Interval(1, 1e2)
  self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel(lengthscale_constraint = constraint1) +
                                                  gpytorch.kernels.RBFKernel(lengthscale_constraint = constraint2))

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

likelihood1 = gpytorch.likelihoods.GaussianLikelihood()
model1 = ExactGPModel(train_x1, train_y1, likelihood1)

likelihood2 = gpytorch.likelihoods.GaussianLikelihood()
model2 = ExactGPModel(train_x2, train_y2, likelihood2)

and got this:

However, in either case I cannot get a realistic recovery of periods with missing data.

[gpleiss]

It looks like you may have a problem that is challenging for optimization. You could try manually adjusting hyperparameters of the kernels. Also consider the Spectral Delta Kernel (https://docs.gpytorch.ai/en/stable/examples/01_Exact_GPs/Spectral_Delta_GP_Regression.html).

[dhenriques62]

Hi @gpleiss

Yes ! I'm just trying the Spectra Delta Kernel,:

class SpectralDeltaGP(gpytorch.models.ExactGP):

def __init__(self, train_x, train_y, num_deltas, noise_init=None):
    likelihood = gpytorch.likelihoods.GaussianLikelihood(noise_constraint=gpytorch.constraints.GreaterThan(1e-11))
    likelihood.register_prior("noise_prior", gpytorch.priors.HorseshoePrior(0.1), "noise")
    likelihood.noise = 1e-2

    super(SpectralDeltaGP, self).__init__(train_x, train_y, likelihood)
    self.mean_module = gpytorch.means.ConstantMean()
    base_covar_module = gpytorch.kernels.SpectralDeltaKernel(
        num_dims=train_x.size(-1),
        num_deltas=num_deltas,
    )
    base_covar_module.initialize_from_data(train_x[0], train_y[0])
    self.covar_module = gpytorch.kernels.ScaleKernel(base_covar_module)

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

model1 = SpectralDeltaGP(train_x1, train_y1, num_deltas=1500)

if torch.cuda.is_available():
    model1 = model1.cuda()

model2 = SpectralDeltaGP(train_x2, train_y2, num_deltas=1500)

if torch.cuda.is_available():
    model2 = model2.cuda()`

but got the following error:

I don't understand why train_x needs to have at least two dimensions.
It only happens with this kernel setup. Could it be a bug?

[gpleiss]

Just shape your data to be N x 1, and that will fix the problem.

[dhenriques62]

Thank you @gpleiss.

It works, but is still too softened. Needs further adjustment: