As noted here
Why is the mean function in Gaussian Process uninteresting?
the mean function is usually not the main focus of the modeling effort, for Gaussian Processes. However, there are cases, such as extrapolation, where we need to use something better than a constant mean function, because otherwise the response of a Gaussian Process with a constant mean function $C$ will revert to just $C+\bar{y}$ "sufficiently far away" from the training data. And "sufficiently far away" can be "very close", if we use a Squared Exponential covariance function, and/or the length-scales which best fit the training data are very small with respect to the "diameter" of the training set.
Excluding the trivial case where the mean function is deterministic (i.e., it's a known function of the inputs but it doesn't depend on the training data, such as for example $\mathbf{c}^T\cdot\mathbf{x}$, with $\mathbf{c}$ a predetermined vector), we have basically two cases:
The mean function is a linear model
This means that the mean function is
$$g(\mathbf{x}\vert\boldsymbol{\beta})=\boldsymbol{\beta}^T\cdot\mathbf{b}(\mathbf{x})$$
where $\boldsymbol{\beta}$ is a vector of unknown parameters, and $\mathbf{b}(\mathbf{x})$ is a fixed set of basis functions, such as for example:
- the monomials of maximum degree $p$, i.e., $\{x_1^{\alpha_1}\dots x_d^{\alpha_d}\vert\sum_{i=1}^d\alpha_i\le p\}$
- the Fourier monomials in $\mathbb{R}^d$, i.e., $\{\exp(i\mathbf{m}\cdot \mathbf{x})\vert \ \Vert\mathbf{m}\Vert_1\le M\}$
- splines (I leave you the pleasure of writing out the multivariate expression)
- etc.
In this case, if we choose a Gaussian prior $\boldsymbol{\beta}\sim\mathcal{N}(\mathbf{b},\boldsymbol{\Sigma})$, then the predictive mean vector and covariance matrix still have an analytical expression, just like in the case of the constant mean function. The expression is a bit cumbersome: you can find it as equations (2.41) of C. E. Rasmussen & C. K. I. Williams, Gaussian Processes for Machine Learning. Note: as in the constant mean function case, these analytical expressions are only exact if the covariance function (kernel) is pre-determined (except at most for a nugget term to accommodate noisy observations). If (as it's nearly always the case) the covariance function contains free hyperparameters whose posterior distribution you need to estimate, based on the training data, then you need to use simulation (MCMC) if you want to perform "exact" Bayesian inference.
The mean function is a nonlinear model
This is the case when, for example:
- the basis functions are themselves functions of the training data
- the number of basis functions depends on the training data
- the mean function is not a linear combination of basis functions (e.g., rational functions)
- etc.
In this case, the only way to compute the predictive mean and covariance is through simulation, even if the covariance function is prespecified. However, I've never seen a practical application of nonlinear mean functions. I guess that, when the data generating process is so complicated that a linear model for the mean function is inappropriate, you either focus on improving/complicating the covariance function, or you use another statistical model instead of a Gaussian Process (for example, a Bayesian Neural Network).
Criteria for selection
Now that you have the expression, you can either perform the selection based on purely heuristic criteria (e.g., WAIC or cross-validation), or based on prior knowledge. For example, if you know from Physics that for $\Vert \mathbf{x} \Vert_2\to\infty$, your response should be a linear function of the inputs, you will select a mean function which is a linear polynomial, if you know that it must become periodic, you will choose a Fourier basis, etc.
Another possible criterion is interpretability: for obvious reasons, a GP is not the most immediately interpretable model, but if you use a linear mean function, then at least asymptotically, when the effects of the kernel have "died out", you can interpret the coefficients of the linear model as a sort of effect size.
Finally, nonconstant mean functions can be used to show the strict relationship between spline models, Generalized Additive Models (GAMs) and Gaussian Processes.
