Spatial mapping with Gaussian processes and nonstationary Fourier features

The use of covariance kernels is ubiquitous in the field of spatial statistics. Kernels allow data to be mapped into high-dimensional feature spaces and can thus extend simple linear additive methods to nonlinear methods with higher order interactions. However, until recently, there has been a strong reliance on a limited class of stationary kernels such as the Matérn or squared exponential, limiting the expressiveness of these modelling approaches. Recent machine learning research has focused on spectral representations to model arbitrary stationary kernels and introduced more general representations that include classes of nonstationary kernels. In this paper, we exploit the connections between Fourier feature representations, Gaussian processes and neural networks to generalise previous approaches and develop a simple and efficient framework to learn arbitrarily complex nonstationary kernel functions directly from the data, while taking care to avoid overfitting using state-of-the-art methods from deep learning. We highlight the very broad array of kernel classes that could be created within this framework. We apply this to a time series dataset and a remote sensing problem involving land surface temperature in Eastern Africa. We show that without increasing the computational or storage complexity, nonstationary kernels can be used to improve generalisation performance and provide more interpretable results.


Introduction
The past decade has seen a tremendous and ubiquitous growth in both data collection and computational resources available for data analysis.In particular, spatiotemporal modelling has been brought to the forefront with applications in finance [1], weather forecasting [2], remote sensing [3][4][5] and demographic/disease mapping [6][7][8].The methodological workhorse of mapping efforts has been Gaussian process (GP) regression [9,10].The reliance on GPs stems from their convenient mathematical framework which allows the modelling of distributions over nonlinear functions.GPs offer robustness to overfitting, a principled way to integrate over hyperparameters, and provide uncertainty intervals.In a GP, every point in some continuous input space is associated with a normally distributed random variable, such that every finite collection of those random variables has a multivariate normal distribution -entirely defined through a mean function µ(•) and a covariance kernel function k(•, •).In many settings, µ(•) = 0 and modelling proceeds through selecting the appropriate kernel function which entirely determines the properties of the GP, and can have a significant influence on both the predictive performance and on the model interpretability [11,12].However, in practice, the kernel function is often (somewhat arbitrarily) set a priori to the squared exponential or Matérn class of kernels [12], justifying this choice by the fact that they model a rich class of functions [13].
While offering an elegant mathematical framework, performing inference with GP models is computationally demanding.Namely, evaluating the GP posterior involves a matrix inversion, which for n observations, requires O(n 3 ) time and O(n 2 ) storage complexity.This makes fitting full GP models prohibitive for any dataset that exceeds a few thousand observations (thereby limiting their use exactly in the regimes where a flexible nonlinear model is of interest).
In response to these limitations, many scalable approximations to GP modelling have been proposed.Examples include inducing points representations [14], the Nyström approximation [9], the Fully Independent Training Conditional (FITC) model [15], stochastic partial differential equation representations [16], and representations in the Fourier domain [17][18][19].Most of these approaches reduce the computational complexity to O(nm 2 ) time and O(nm) storage for m n playing the role of the number of dimensions of the functional space under consideration (e.g. the number of inducing points or the number of frequencies in Fourier representations).
In this contribution, we will focus on large-scale Fourier representations of GPs.These methods traditionally rely on the strong assumption of the stationarity (or shift-invariance) of kernel functions, which is made in the vast majority of applications (and is indeed satisfied by the most often used squared exponential and Matérn kernels).Stationarity in the spatio-temporal data means that the similarity between two responses in space and time does not depend on the location and time itself, but only on the difference (or lag) between them, i.e. kernel function can be written as k(x 1 , x 2 ) = κ(x 1 − x 2 ) for some function κ.Several recent works [18,20,21] consider flexible families of kernels based on Fourier representations, thus avoiding the need to choose a specific kernel a priori and allowing the kernel to be learned from the data, but these approaches are restricted to the stationary case.In many applications, particularly when data is rich, relaxing the assumption of stationarity can greatly improve generalisation performance [11].
To address this, recent work in [12,19] note that a more general spectral characterisation exists that includes nonstationary kernels [12,22] and uses it to construct nonstationary kernel families.In this paper, we build on the work of [18,19,23] and develop a simple and practicable framework for learning spatiotemporal nonstationary kernel functions directly from the data by exploiting the connections between Fourier feature representations, Gaussian processes and neural networks [9].Specifically, we directly learn frequencies in nonstationary spectral kernel representations using an appropriate neural network architecture, and adopt techniques used for deep learning regularisation [24] to prevent overfitting.We demonstrate the utility of the proposed method for learning nonstationary kernel functions in a time series example and in spatial mapping of land surface temperature in East Africa.
2 Methods and Theory

Gaussian Process Regression
Gaussian process regression (GPR) takes a training dataset D = {(x i , y i )} n i=1 where y i ∈ R is real-valued response/output and x i ∈ R D is a D-dimensional input vector.The response y i and the input x i are connected via the observation model GPR is a Bayesian non-parametric approach that imposes a prior distribution on functions f , namely a GP prior, such that any vector f = [f (x 1 ), .., f (x m )] of a finite number of evaluations of f follows a multivariate normal distribution f ∼ N (0, K xx ), where the covariance matrix Throughout this paper we will assume that the mean function of the GP prior is µ = 0, however, all the approaches in this paper can be easily extended to include a mean function [25].
In stationary settings, k(x i , x j ) = κ(x i − x j ) for some function κ(δ).A popular choice is the automatic relevance determination (ARD) kernel [9], given by κ(δ) = τ 2 exp(−δ Λδ) where τ 2 > 0 and Λ = diag(λ 1 , . . ., λ D ).Kernel k will typically have hyperparameters θ (e.g.θ = [τ, λ 1 , . . ., λ D ] for the ARD kernel) and one can thus consider a Bayesian hierarchical model: The posterior predictive distribution is straightforward to obtain from the conditioning properties of multivariate normal distributions.For a new input x * , we can find the posterior predictive distribution of the associated response where and it is understood that the dependence on θ is through the kernel k = k θ .The computational complexity in prediction stems from the matrix inversion ( The marginal likelihood (also called model evidence) of the vector of outputs y = [y 1 , . . ., y n ], is given by p(y|θ) = p(y|f , θ)p(f |θ)df , is obtained by integrating out the GP evaluations f from the likelihood of the observation model.Maximising the marginal likelihood over hyperparameters allows for automatic regularisation and hence for selecting an appropriate model complexity.For a normal observation model in (1), the log marginal likelihood is available in closed form Computing the inverse and determinant in (6) are computationally demanding -moreover, they need to be computed for every hyperparameter value θ under consideration.To allow for computational tractability, we will use an approximation of K xx based on Fourier features (see section 2.3 and 2.4).
Alternative representations can easily be used such as the primal/dual representations in a closely related frequentist method, kernel ridge regression (KRR) [26].In contrast to KRR, optimising the marginal likelihood as above retains the same computational complexity while providing uncertainty bounds and automatic regularisation without having to tune a regularisation hyperparameter.However, the maximisation problem of ( 6) is non-convex thereby limiting the chance of finding a global optimum, but instead relying on reasonable local optima [9].

Random Fourier Feature mappings
The Wiener-Khintchine theorem states that the power spectrum and the autocorrelation function of a random process constitute a Fourier pair.Given this, random Fourier feature mappings and similar methodologies [17-19, 21, 23] appeal to Bochner's theorem to reformulate the kernel function in terms of its spectral density.
Theorem 1. (Bochner's Theorem) A stationary continuous kernel k(x i , x j ) = κ(x i − x j ) on R d is positive definite if and only if κ(δ) is the Fourier transform of a non-negative measure.
Hence, for an appropriately scaled shift invariant complex kernel κ(δ), i.e. for κ(0) = 1, Bochner's Theorem ensures that its inverse Fourier Transform is a probability measure: Thus, Bochner's Theorem introduces the duality between stationary kernels and the spectral measures P(dω).Note that the scale parameter of the kernel, i.e. σ 2 f = κ(0) can be trivially added back into the kernel construction by rescaling.Table 1 shows some popular kernel functions and their respective spectral densities.

Summary table of kernels and their spectral densities
By taking the real part of equation 7 (since we are commonly interested only in real-valued kernels in the context of GP modelling) and performing standard Monte Carlo integration, we can derive to a finite-dimensional, reduced rank approximation of the kernel function where ∼ P and we denoted For a covariate design matrix X ∈ R n×D (with rows corresponding to data vectors x 1 , . . ., x n ), and frequency matrix Ω ∈ R m×D (with rows coresponding to frequencies ω 1 , . . ., ω m ), we let Φ x = cos(XΩ ) sin(XΩ ) be a n × 2m matrix referred to as the feature map of the dataset.
The estimated covariance matrix can be computed as T which has rank at most 2m.Substituting K xx into (6) now allows rewriting the determinant and the inverse in terms of the 2m × 2m matrix Φ x T Φ x , thereby reducing the computational cost of inference from O(n 3 ) to O(nm 2 ), where m is the number of Monte Carlo samples/frequencies.Typically, m n.
In particular, by defining where σ 2 n is the observation noise variance and σ 2 f = κ(0) is the kernel scale parameter, and taking R = chol(A) to be the Cholesky factor of A, we first calculate vectors α 1 , α 2 solving the linear systems of equations Rα 1 = Φ x T y and The log marginal likelihood can then be computed efficiently as: Additionally, the posterior predictive mean and variance can be estimated as There are two important disadvantages of standard random Fourier features as proposed by [17]: firstly, only stationary (shift invariant) kernels can be approximated, and secondly we have to select a priori a specific class of kernels and their corresponding spectral distributions (e.g. Table 1).In this paper, we address both of these limitations, with a goal to construct methods to learn a nonstationary kernel from the data, while preserving the computational efficiency of random Fourier features.
While we can think about the quantities in (16) as giving approximations to the full GP inference with a given kernel k, they are in fact performing exact GP calculations for another kernel k defined using the explicit feature map Φ x defined through frequencies sampled from the spectral measure of k.We can thus think about these feature maps as parametrizing a family of kernels in their own right and treat frequencies ω 1 , . . ., ω m as kernel parameters to be optimized, i.e. learned from the data by maximizing the log marginal likelihood.It should be noted that dropping the imaginary part of our kernel symmetrizes the spectral measure allowing us to use any P(dω) -regardless of its symmetry properties, we will still have a real-valued kernel.In particular, one can use an empirical spectral measure defined by any finite set of frequencies.

Nonstationary random Fourier features
Contrary to stationary kernels, which only depend on the lag vector i.e. δ = x i −x j , nonstationary kernels depend on the inputs themselves.A simple example of a nonstationary kernel would be the polynomial kernel defined as: To extend the stationary random feature mapping to nonstationary kernels, following [12,19,22], we will need to use a more general spectral characterisation of positive definite functions which encompasses stationary and nonstationary kernels.
From the above theorem, a nonstationary kernel can be characterized by a spectral measure µ(dω 1 , dω 2 ) on the product space R D × R D .Again, without loss of generality we can assume that µ is a probability measure.If µ is concentrated along the diagonal, ω 1 = ω 2 , we recover the spectral representation of stationary kernels in the previous section.However, exploiting this more general characterisation, we can construct feature mappings for nonstationary kernels.
Just like in the stationary case, we can approximate (19) using Monte Carlo integration.In order to ensure a valid positive semi-definite spectral density we first have to symmetrize f (ω 1 , ω 2 ) by ensuring f (ω 1 , ω 2 ) = f (ω 2 , ω 1 ) and including the diagonal components f (ω 1 , ω 1 ) and f (ω 2 , ω 2 ) [23].We can take a general form of density g on the product space and "symmetrize": Once again using Monte Carlo integration we get: where Hence, by denoting Ω l ∈ R m×D (with rows coresponding to frequencies ω l 1 , . . ., ω l m ) for l = 1, 2 as before, we obtain the corresponding feature map for the approximated kernel as an n × 2m matrix and can be condensed to an identical form as in the stationary case.
The non-stationarity in equation ( 21) arises from the product of differing locations x 1 and x 2 by different frequencies ω 1 k , ω 2 k , hence making the kernel dependent on the values of x 1 and x 2 and not only the lag vector.If the frequencies were exactly the same we just revert back to the stationary case.The complete construction of random Fourier feature approximation is summarized in the algorithm below.

Algorithm 1 Random Fourier features for nonstationary kernels
Input: spectral measure µ, dataset X, number of frequencies m However, just like in the stationary case, we can think about nonstationary Fourier feature maps as parametrizing a family of kernels and treat frequencies {(ω 1 k , ω 2 k )} m k=1 as kernel parameters to be learned by maximizing the log marginal likelihood, which is an approach we pursue in this work.Again, symmetrization due to dropping imaginary parts implies that any empirical spectral measure is valid (there are no constraints on the frequencies).

On the choice of spectral measure in non-stationary case
Using the characterisation in equation ( 21) one only requires the specification of the (Lebesgue-Stieltjes measurable) distribution f (ω 1 , ω 2 ) in order to construct a nonstationary kernel.This very general formulation allows us to create the full spectrum encompassing both simple and highly complex kernels.
In the simplest case, f (ω 1 , ω 2 ) = f 1 (ω 1 )f 2 (ω 2 ), i.e. it can be a product of popular spectral densities listed in Table 1.Furthermore, one could consider cases where these individual spectral densities factorize further across dimensions.This corresponds to a notion of separability.
In spatio-temporal data, separability can be very useful as it enables interpretation of the relationship between the covariates as well as computationally efficient estimations and inferences [27].Practical implementation is straightforward; consider the classic spatio-temporal setting with 3 covariates -longitude, latitude and time.When drawing random samples of ω l = (ω 1 l , ω 2 l , ω 3 l ) where l ∈ {1, 2}, we could define the ω i l to come from different distributions, allowing us to individually model each input dimension.If the distribution on frequencies are independent across dimensions then we see that if A practical example for spatio-temporal modelling of a nonstationary separable kernel would be generating a four dimensional (ω 2 ) , sample from independent Gaussian distributions (whose spectral density corresponds to a squared exponential kernel) representing nonstationary spatial coordinates, and a two dimensional (ω 3  1 , ω 3 2 ) from a Student-t distribution with 0.5 degrees of freedom (whose spectral density corresponds to a Matérn 1/2 kernel or exponential kernel) representing temporal coordinates.
To move from separable to non-separable, nonstationary kernels one only needs to introduce some dependence structure within ω 1 or ω 2 i.e. across feature dimensions, such as for example using the multivariate normal distribution in R D , in order to prevent the factorization in equation (22).The correlation structure in these multivariate distributions are what creates the non-separability.
To create non-separable kernels with different spectral densities along each feature dimension copulas can be used.An example in a spatial (latitude, longitude feature dimensions) setting using the Gaussian copula, would involve generating samples for ω 1 or ω 2 ∈ R 2 (or both) from a multivariate normal distribution {ω 1 k } m k=1 i.i.d.
∼ N (0, Σ), pass these through the Gaussian cumulative distribution function, and then used in the quantile function of another distribution )).This transformation can also be done using different Λs along different feature dimensions.Alternative copulas can be readily used, including the popular Archimedian Copulas: Clayton, Frank and Gumbel [28].Additionally, mixtures of multivariate normals can be used [21,23] to create arbitrarily complex non-separable and nonstationary kernels.Given sufficient components any probability density function can be approximated to the desired accuracy.
In this paper, we focus on the most general case where the frequencies {(ω 1 k , ω 2 k )} m k=1 are treated as kernel parameters and are learnt directly from the data by optimising the marginal likelihood, i.e. they are not associated to any specific family of joint distributions.This approach allows us to directly learn nonstationary kernels of arbitrary complexity as m increases.However, a major problem with such a heavily overparametrized kernel is the possibility of overfitting.
Stationary examples of learning frequencies directly from the data [18,29,30] have been known to overfit despite the regularisation due to working with marginal likelihood.This problem is further exacerbated in high-dimensional settings, such as those in spatio-temporal mapping with covariates.In this paper, we include an additional regularisation inspired by dropout [24] which prevents the co-adaptation of the learnt frequencies ω 1 , ω 2 .

Gaussian dropout regularisation
Dropout [24] is a regularisation technique introduced to mitigate overfitting in deep neural networks.In its simplest form, dropout involves setting features/matrix entries to zero with probability q = 1 − p, i.e. according to a Bernoulli(p) for each feature.The main motivation behind the algorithm is to prevent co-adaptation by forcing features to be robust and rely on population behaviour.This prevents individual features from overfitting to idiosyncrasies of the data.
Using standard dropout, where zeros are introduced into the frequencies {(ω 1 k , ω 2 k )} m k=1 can be problematic due to the trigonometric transformations in the projected features.An alternative to dropout that has been shown to be just as effective if not better is Gaussian dropout [24,31].Regularisation via Gaussian dropout involves augmenting our sample distribution as ] (see figure 1).As with dropout, this approach prevented our population Monte Carlo sample from coadapting, and ensured that the learnt frequencies are robust and not overfitting noise in the data.An additional benefit of this procedure over improving generalisation error and preventing overfitting was to speed up the convergence of gradient descent optimisers through escaping saddle points more effectively [32].The noise parameter σ p defines the degree of regularisation and is a hyperparameter that needs to be tuned.However, we found in practice when coupled with an early stopping procedure, learning the frequencies is robust to sensible choices of σ p .

Google daily high stock price
To demonstrate the use of the developed method and the utility of nonstationary modelling, we consider time series data of the daily high stock price of Google spanning 3295 days from 19th August 2004 to 20th September 2017.We set x ∈ {1, . . ., 3295} and y = log(Stock high ).For the stationary case we use vanilla random Fourier features [17,18] with the squared exponential kernel (Gaussian spectral density) and m = 600 fixed frequencies.For the nonstationary case we use m = 300 frequencies for each ω 1 and for ω 2 .We performed a sensitivity analysis to check that no improvements in either the log marginal likelihood or testing error resulted from using more features.Cross-validation was performed using a 70 − 30 training-testing split averaged over 20 independent runs and testing performance evaluated via the mean squared error and correlation.
Optimisation of the log marginal likelihood was performed using ADAM [33] gradient descent in TensorFlow [34] with early stopping and patience [35].
Figure 2 (top left) shows the comparison in the optimisation paths of the negative log marginal likelihood between the two methods.It is clear that the nonstationary approach reaches a lower minima than the vanilla random Fourier features approach.This is also mirrored in the testing accuracy over the 20 independent runs where our approach achieves a mean squared error and correlation of 3.29 × 10 −5 and 0.999, while the vanilla Fourier features approach achieves a mean squared error and correlation of 5.69 × 10 −5 and 0.987.Of note is the impact of our Gaussian dropout regularisation, which, through the injection of noise, appears to converge faster and avoid plateaus.This is entirely in keeping with previous experiences using dropout variants [31] and highlights an added benefit over using only ridge (weight decay) regularisation.Figure 3 (top) shows the overall fits compared to the raw data.Both methods appear to fit the data very well, as reflected in the testing statistics, but when examining a zoomed-in transect it is clear that the learnt nonstationary features fit the data better than the vanilla random features by allowing variable degree of smoothness in time.The combination of nonstationarity and kernel flexibility allowed us to learn a much better characterisation of the data patterns without overfitting.The covariance matrix comparisons (Figure 3 bottom) further highlight this point where the learnt nonstationary covariance matrix shares some similarities with the vanilla random features covariance matrix, such as the concentration on the diagonal, but exhibits a much greater degree of texture.The histograms in Figure 3 provide another perspective on the covariance structure, where the vanilla features are by design Gaussian distributed, but learnt nonstationary frequencies are far from Gaussian (Kolmogorov-Smirnov test p-value< 10 −16 ) exhibiting skewness and heavy tails.Additionally, the differences between the learnt frequencies ω 1 and ω 2 show that, not only is the learnt kernel far from Gaussian, but that it is indeed also nonstationary.This simple example also highlights the potential defficiencies of choosing kernels/frequencies a priori.

Spatial temperature anomaly for East Africa in 2016
We next consider MOD11A2 Land Surface Temperature (LST) 8-day composite data [36], gap filled for cloud cover images [5] and averaged to a synoptic yearly mean for 2016.To replicate common situations for spatial mapping, such as interpolation from sparse remote sensed sites or cluster house hold survey locations [25] we randomly sample 4000 LST locations (only ∼ 4% of the total) from the East Africa region (see Figures 4 and 5).We set x ∈ R 2 = {Latitude, Longitude} i.e. using only the spatial coordinates as covariates, and use the LST temperature anomaly as the response.We apply our nonstationary approach, learning m = 600 frequencies for ω 1 and ω 2 each.Cross validation was evaluated over all pixels excluding the training set (∼ 83, 000) and averaged over 10 independent runs with testing performance evaluated via the mean squared error and correlation.Our final performance estimates were 4.23 and 0.91 for the mean squared error and correlation respectively.Figure 4 shows our predicted surface (left) compared to the actual data (right).Our model shows strong correspondence to the underlying data and highlights the suitability of using our approach in settings where no relevant covariates exist outside of the spatial coordinates.Figure 5 shows 3 randomly sampled points and the covariance patterns around those points.For comparison Figure 6 shows the equivalent plot when using a stationary squared exponential kernel.In stark contrast to the stationary covariance function, which has an identical structure for all three points, the nonstationary kernel shows considerable heterogeneity in both patterns and shapes.Interestingly the learnt lengthscale/bandwidth seems to be much smaller in the stationary case than the nonstationary case, we hypothesise that this is due to the inability of the stationary kernel to learn the rich covariance structure needed to accurately model temperature anomaly.Intuitively, nonstationarity allows locally dependent covariance structures which conform to the properties of a particular location and imply (on average) larger similarity of nearby outputs and better generalisation ability.In contrast, stationary kernels are trying to fit one covariance structure to all locations and as a result end up with a much shorter lengthscale as it needs to apply to all directions from all locations.Our results are in concordance with other studies showing that temperature anomaly data is nonstationary [19,23,37].This interpolation problem can readily be expanded into multiple dimensions including time and other covariates.

Discussion
We have shown that nonstationary kernels of arbitrary complexity are as easy to model as stationary ones, and can be learnt with sufficient efficiency to be applicable to datasets of all sizes.The qualitative superiority of predictions when using nonstationary kernels has previously been noted [11].In many applications, such as in epidemiology where data can be noisy, generalisation accuracy is not the only measure of model performance, and there is a need for models that conform to known biological constraints and external field data.The flexibility of nonstationary kernels allows for more plausible realities to be modelled without the assumption of stationarity limiting the expressiveness of the predictions.It has also been noted that while nonstationary GPs give more sensible results than stationary GPs, they often show little generalisation improvement [11].For the examples in this work we show clear improvements in generalisation accuracy when using nonstationary kernels.We conjecture that the differences in generalisation performance are likely due to the same reasons limiting neural network performance a decade ago [38] -namely, a combination of small, poor quality data and a lack of generality in the underlying specification.Given more generalised specifications, such as those introduced in this paper, coupled with the current trend of increasing quantities of high quality data [39] we believe nonstationary approaches will be more and more relevant in spatio-temporal modelling.
There has long been codes of practice on which kernel to use on which spatial dataset [10] based on a priori assumptions about the roughness of the underlying process.Using the approach introduced in this paper, ad hoc choices of kernel and decisions on stationarity versus nonstationarity may no longer be needed as it may be possible to learn the kernel automatically from the data.For example, our approach can be easily modified to vary the degree of nonstationarity according to patterns found in the data.
In this work we have focused on optimising the marginal likelihood in Gaussian Process regression and added extra regularisation via Gaussian dropout.However, for non-Gaussian observation models the marginal likelihood cannot be obtained in a closed form.In these settings, one may resort to frequentist methods instead and resort to variational [30], approximate [40,41] or suitable MCMC [42] approaches in order to provide uncertainty measures.For very large models with non-Gaussian observation models, stochastic gradient descent in mini-batches [43] or stochastic gradient Bayesian methods [44] can be used.
The matrix Φ that results from the Fourier feature expansion can be thought of as a hidden layer in a single layer neural network.This formulation explicitly connects the random Fourier feature approach and single layer neural networks using learnable (or random) nodes and specific trigonometric activation functions.Our approach can therefore be used in deep learning settings while retaining an explicit link to kernel methods and their interpretability [45].

Figure 1 .
Figure 1.Histogram of the Euclidean norm of a covariance matrix ΦΦ T with Gaussian dropout of σ p = 0.05.The black line is the norm of ΦΦ T without noise

Figure 5 .Figure 6 .
Figure 5. Covariance matrix images for 3 random points showing different covariance structures due to nonstationarity