# Choosing kernel for Gaussian process regression (kriging) with strong directional correlation

My problem is very simple.

I intend to use kriging (for various reasons) to perform regression between the trajectories marked with red dots on the below plot, where a simple linear-regression surface has been superimposed as well.

As you can see there is a strong correlation in a particular direction of the four 'curves'.

How might one build in such strong prior-knowledge in the choice of kernel? Or indeed, what might an appropriate choice of kernel be in and off itself?

Thanks!

• can you explain why are you ruling out multivariate splines? Your case seems a perfect fit for the method. If you're choosing GPR because it comes with handy confidence/prediction intervals, you can compute them with splines too, of course. Also: this is regression, not interpolation, right? Your red dots are observations with error. I'm asking just to be sure, because they look extremely smooth, so I was wondering if they were results from a computer simulation (deterministic) or experimental observation (with random error). – DeltaIV Dec 13 '16 at 12:47
• @DeltaIV thank you for you incisive comment. I am using both. I am comparing them, and my reasoning for using GPR is my response variable values that I am trying to model are highly correlated and the nature of the correlation is quite complex (or so I believe) - so other techniques for handling this (like ARMA), may not be a good fit. – Astrid Dec 13 '16 at 14:43
• @DeltaIV in response to your questions: a) I picked GPR because of the confidence intervals that come with it, I was not aware that you could get these with splines too. Will investigate! b) That was something I had not thought of. So the red dots are computer simulations (deterministic) which are computed from experimental observation (which naturally have error). I am really interested in the region between the red-dot curves. – Astrid Dec 13 '16 at 14:47
• for point a) yes, and it's very easy! Spline regression is a linear model, meaning that it's linear in the model parameters (if you're using fixed knots). Thus it's actually very easy to compute confidence/prediction intervals. The same is true for Fourier regression, which may be more indicated in your case (periodic latent function). b) I'm not sure I understand your setting...if you have experimental errors on the inputs, but no error on the outputs, this is quite different from the usual GPR setting. Suppose you could acquire another data set at the same input points [1/2] – DeltaIV Dec 13 '16 at 16:45
• [2/2], i.e., if the $(x_1,x_2)$ coordinates of the red dots are the same. In practice this is not possible, but suppose for a moment that you could. Then, if the $(x_1,x_2)$ coordinates of each point in the training set are the same, would the response values (the $y$ coordinates of the red dots) be the same or not? – DeltaIV Dec 13 '16 at 16:46

You have very smooth trajectories, with different correlation lengths between trajectories and along a trajectory. Simple Squared Exponential Kernel with a different correlation length for each coordinate $x_1,x_2$ might suffice:

$$K(\mathbf{x},\mathbf{x}')=\sigma \text{exp}\left(-\left(\frac{(x_1-x_1')^2}{2l_1^2}+\frac{(x_2-x_2')^2}{2l_2^2}\right)\right)=\sigma \text{exp}\left(-\frac{(x_1-x_1')^2}{2l_1^2}\right)\text{exp}\left(-\frac{(x_2-x_2')^2}{2l_2^2}\right)$$

where $\sigma,l_1,l_2$ and the error standard deviation $\sigma_n$ would be estimated from data, either with MLE or by computing their posterior distribution in a Bayesian setting.

I said that this might work, because it depends also on the use you want to make of this Gaussian Process. Given a GP prior with zero mean and this kernel, then, outside from the region where you have data (for example, in your picture the interval $I=[0,2750]\text{x}[0.6,1.6]$ or something like that), the posterior conditional mean will invariably go to 0. This will happen more quickly in the direction with shorter correlation length, and over a longer distance in the direction with longer correlation length, but sooner or later it will happen nonetheless (more sooner than later with an SE kernel). Thus, if you want to use the GP for interpolation inside $I$ or not too much outside, you'll be ok, but for extrapolation well outside $I$ you'll get in trouble.

Fixing this issue requires prior knowledge about the behavior of your ideal function $f(\mathbf{x})$ outside $I$. For example, if you expect $f$ to asymptotically become linear, then use a linear mean function in your GP prior, whose parameters you will estimate from data. Another possibility: by looking at your data, you might think that $f$ becomes asymptotically periodic. Based on the red dots, it might seem a plausible hypothesis, but you cannot say for sure: you must have some knowledge of the physical process which generated your data. In this case, you may either resort to a periodic kernel, or again use a non-zero mean function, in this case a truncated Fourier series whose coefficients you would, as always, estimate from data.

An important practical note: $x_1$ and $x_2$ have very different magnitudes. To improve the conditioning of your covariance matrix, it would be good to rescale $x_1$ and $x_2$ so that they have similar orders of magnitude, and since you're there already, you could also center them. Of course, if you do this, then when you use the GP for prediction at a new point $\mathbf{x}^*=(x_1^*,x_2^*)$, you need to rescale its coordinates before evaluating the posterior mean and variance of the GP. This rescaling is not strictly necessary, because judging by your plot you seem to have an "easy" problem for GPR (relatively few data points, and no pairs of data points which are much closer than the average distance between data points). Anyway, the Cholesky decomposition can always be a dodgy step in GPR, so I'd rescale if I were in you, expecially if the physics of the problem gives you reasonable measures of scale and location for $x_1$ and $x_2$ (if it doesn't, you could estimate them from your data set).

EDIT: from the comments it's clear that the latent function is periodic in $x_1$, and that $x_1$ is a time, thus I'm changing my notation: $w=f(t,x)$, where $f$ is the latent function, $t$ is time (it was $x_1$ in my old notation) and $x$ is space ($x_2$ in my old notation). In this case I definitely don't recommend the SE kernel. This periodic-SE kernel would probably be a better idea:

$$K((t,x),(t',x'))=\sigma \text{exp}\left(-2\frac{\text{sin}^2\left(\pi\frac{|t-t'|^2}{T}\right)}{l_t^2}\right)\text{exp}\left(-\frac{(x-x')^2}{2l_x^2}\right)$$

If you know already know $T$ (the period) from physics, you could plug it in in the above kernel, but it may be more accurate to estimate it from data. Of course, if you end up estimating a $T$ value which is consistently different from what physics tells you, there is a problem either with the model, the data or the estimation algorithm.

• The data is highly periodic! Shoot, I forgot to add that in the question; this is but one snapshot, of one complete cycle of the sinusoidal-looking latent function. Given that this is actually the case; a periodic kernel in the domain of each red-dot curve (of which there are four) ought then be the correct choice? – Astrid Dec 13 '16 at 14:50
• it would seem Duvenaud suggests using a locally periodic kernel, for functions that do not exactly repeat themselves. This would basically be my case. – Astrid Dec 13 '16 at 14:55
• Do you mean $x_2$ and not $y_2$ in the above answer? – Astrid Dec 13 '16 at 15:01
• @Astrid yes, of course, $x_2$, not $y_2$, it was a typing mistake. A very important question: let's indicate with $x_1$ the coordinate along the trajectory of each curve (thus, $x_1 \in [0,2750]$ in your plot, or something like that). Is your latent function $f(x_1,x_2)$ periodic only along $x_1$, or also along $x_2$, obviously with a different period? BTW please read the new version of my answer, I added an important advice. – DeltaIV Dec 13 '16 at 16:27
• My function is $only$ periodic along $x_1$ in your notation, and not, say $x_2 \in [0.5,1.8]$. I expect there to be little prior knowledge about the behaviour of $x_2$; if anything I expect the ridge in the middle (see plot) to become wider (in $x_1$) and deeper (in $x_3 \in [10,-25]$) as $x_2 \rightarrow 0$. – Astrid Dec 13 '16 at 22:20