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!  
 A: 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.
