Prescribing asymptotic shape in Gaussian process regression

Question

A feature of the Bayesian approach to Gaussian processes is the fact that asymptotic behaviour far away from data is controlled by the prior mean function. In most cases covered in the literature the prior mean is either a function independent of the data (such as a constant, often zero) or a function in a low dimensional function space, such as a polynomial, whose posterior parameters are specified by data.

Both of these approaches (constant or few parameters) do not work in my application. Prior knowledge tells me the mean function is asymptotically linear, i.e. its second derivative vanishes, but the gradient will depend on the data. Hence constant priors will not work: The gradient is not necessarily zero. Parametrising the mean is not feasible, since my problem is high dimensional (multivariate time series data, dimension on the order of hundreds) and a reasonable parametrisation would require hundreds of parameters, leading to overfitting.

Being twice differentiable almost everywhere is not a problem in my application, hence it would be great if, instead of specifying mean zero or $f=0$, I could specify $f''=0$ as a prior mean. Theoretically this should work since $f''$ is also a Gaussian process and I could condition $f''$ on observations of $f$. Practically. I have a problem, since the correlations of $f$ and $f''$ involve high dimensional double integrals of the kernel, which seem infeasible to calculate.

So my question is: What is the best way to approach Gaussian regression, in case your prior belief about the process is that it is asymptotically linear or $f''=0$?

Answer 1

I'm not sure why overfitting would be an issue here. Suppose you have $d=100$ where $d$ is the number of the input variables. How many training samples do you have? If you were using a Squared Exponential kernel, then unless you're using an isotropic one (in $d=100$ dimensions?) you would already have (at least) $d+2$ hyperparameters, which add way more flexibility to your model than the $d+1$ parameters of the linear mean function. Thus I would expect that you have already enough training data to be able to fit (or to compute a posterior distribution for) $d+2$ hyperparameters. It doesn't seem evident to me that the limited extra flexibility added by a linear trend would make your nonparametric model much more flexible.

Of course I could be wrong: maybe you already tried fitting a linear mean function, and you found that it led to a sharp increase in the roll forward cross-validation. Or maybe you are using a kernel with way less than $d$ parameters, because in your specific application that makes sense (data could be living in a low-dimensional manifold of the high-dimensional input space, and you have some prior knowledge about the manifold). In that case, what about using regularization? In regression, the usual way to add a linear mean function to a GP is to put a Gaussian prior $\boldsymbol{\beta}\sim N(\mathbf{b},B)$ on the $d+1$ coefficients of the mean function. Then, people usually let $B^{-1}\to O$, where $O$ is the matrix of zeros, i.e., the prior becomes vague, in order to obtain prediction equations which don't depend on $\mathbf{b}$ and $B$. However, you're not forced to use a vague prior: you could just use the prior

$$\boldsymbol{\beta}\sim N(\mathbf{0},\lambda I)$$

where of course $\lambda>0$. In other words, you put independent Gaussian priors centered at 0 and with a common variance on the coefficients. This prior will shrink the coefficients towards 0 and prevent overfitting. As a matter of fact, if you were doing Bayesian Linear Regression, instead than Bayesian Gaussian Process Regression, this would correspond to what's called Bayesian Ridge Regression, which allows to fit models even in the $n<d$ case (i.e., when you have less training samples than parameters).

You could easily generalize this to the case in which the covariance is not diagonal, if you have a priori knowledge of the correlation among the coefficients, or even to the case where they are not normally distributed. For example, putting double-exponential, independent priors on each coefficient would again shrink coefficients to 0. Again, in Bayesian Linear Regression, instead than Bayesian Gaussian Process Regression, this would correspond to the Bayesian Lasso. The drawback here is that, once $\boldsymbol{\beta}$ is not normally distributed, then you don't have anymore an analytical expression for the posterior mean of the coefficients: you need an MCMC sampler which works well in high dimensions (for example something like Hamiltonian Monte Carlo).

Answer 2

Caveat: This may be a wacky idea

My suggestion is that you place a GP on $f''$ with zero mean and given variance. Then the observations can be integral observations (simply integrate the kernel as necessary) and so too are the predicted points,

$$ \begin{bmatrix} k(a,b) & \int \int k(a,b) db \\ \int \int k(a,b) da & \int \int \int \int k(a,b)da db \\ \end{bmatrix} $$

where aa and b are replaced with the corresponding x values.