Periodic Gaussian Process - covariance matrix not positive semidefinite


#1

Hi all,
I am trying to model a periodic function as a Gaussian process. The modelled function has not other constraints than the need to be periodic and smooth. My initial idea was to use squared exponential kernel, and only redefine the distance metric, so that it treats the time variable as a circle. However, for certain combinations of parameters this fails to yield a positive semi-definite matrix. Any idea what might be wrong? Or a pointer to some literature on GPs using different distance metrics than Euclidean distance?

Thanks a lot!

Here is R code that shows the problem (I encountered this why preparing simulated data for a Stan model in R)

generate_periodic_gp <- function(time, scale, length) {
  # Construct the squared exponential covariance matrix
  cov_matrix = array(0,c(length(time), length(time)));
  maxTime = max(time)
  minTime = min(time)
  for(i in 1:length(time))  {
    cov_matrix[i,i] = scale^2 + 0.00000001; #Adding jitter to ensure positive semi-definiteness
    if (i < length(time)) {
      for(j in (i+1):length(time)) {
        #Take either the normal distance or the distance wrapped around the edge, whichever is shorter
        distance = min(abs(time[i] - time[j]), min(time[i], time[j]) - minTime + maxTime - max(time[i], time[j]))

        covariance = (scale^2) * exp( (-0.5 / (length ^ 2)) * (distance ^ 2) );
        cov_matrix[i,j] = covariance
        cov_matrix[j,i] = covariance
      }
    }
  }
  # Draw from the multinormal distribution using the cholesky decomposition of the cov. matrix
  chol_cov = t(chol(cov_matrix));
  values = chol_cov %*% rnorm(length(time));
  return(values)
}

times = seq(1,10, by =0.1)
values = generate_periodic_gp(times,1,0.7)
plot(times, values, type ="l") #Works as expected, gives a smooth curve and values[1] == values[length(values)]
values = generate_periodic_gp(times,1,0.8) #Fails, not positive definite

#2

@flaxter has committed the crime of talking about periodic kernels in the last few days, so maybe he can give better comments here, but I’ll try to just provide pointers.

Check out: Approximate GPs with Spectral Stuff

The Stan model at the top has Fourier in the name. Not sure what it is, but anything Fourier seems like it’d be good for periodic stuff. The arxiv paper on the Poisson intensity estimation he links there has another kernel defined on a periodic domain.

I just Googled periodic RBF and this page came up: http://www.cs.toronto.edu/~duvenaud/cookbook/index.html . It has some sorta strange periodic kernel there that might work!


#3

Thanks for the help - periodic kernel from http://www.cs.toronto.edu/~duvenaud/cookbook/index.html does the trick. Actually, I played before with the periodic kernel from Wikipedia (very slightly different formulation) and it gave weird results, but it must have been an implementation bug on my side, because now it works :-)

For future travellers, the formula for the periodic kernel is:

covariance[i,j] = (scale^2) * exp(-2*(sin(PI * abs(x[i] - x[j]) / period)^2) / (length ^ 2))

Thanks again!


#4

Sorry for the late reply–you could also use the spectral mixture kernel (https://arxiv.org/pdf/1302.4245.pdf). Here is Stan code for a two component mixture (i.e. an additive kernel with two different periodic components) in case it’s helpful to someone.

  matrix[n1, n1] Sigma1;

  for (i in 1:n1) {
    Sigma1[i, i] <- var1 + var2;
    for (j in (i+1):n1) {
      Sigma1[i, j] <- var1 * exp(-(x1[i]-x1[j])^2*bw1) * cos(twopi * fabs(x1[i] - x1[j]) * period1) +
                      var2 * exp(-(x1[i]-x1[j])^2*bw2) * cos(twopi * fabs(x1[i] - x1[j]) * period2);
      Sigma1[j, i] <- Sigma1[i, j];
    }
  }

Adding Gaussian Process Covariance Functions