Non-Stationary Spectral Kernels

Introduction

Gaussian processes are a flexible method for non-linear regression . They define a distribution over functions, and their performance depends heavily on the covariance function that constrains the function values. Gaussian processes interpolate function values by considering the value of functions at other similar points, as defined by the kernel function. Standard kernels, such as the Gaussian kernel, lead to smooth neighborhood-dominated interpolation that is oblivious of any periodic or long-range connections within the input space, and can not adapt the similarity metric to different parts of the input space.

Two key properties of covariance functions are stationarity and monotony. A stationary kernel $K(x,x^{\prime})=K(x+a,x^{\prime}+a)$ is a function only of the distance $x-x^{\prime}$ and not directly the value of $x$. Hence it encodes an identical similarity notion across the input space, while a monotonic kernel decreases over distance. Kernels that are both stationary and monotonic, such as the Gaussian and Matérn kernels, can encode neither input-dependent function dynamics nor long-range correlations within the input space. Non-monotonic and non-stationary functions are commonly encountered in realistic signal processing , time series analysis , bioinformatics , and in geostatistics applications .

Recently, several authors have explored kernels that are either non-monotonic or non-stationary. A non-monotonic kernel can reveal informative manifolds over the input space by coupling distant points due to periodic or other effects. Non-monotonic kernels have been derived from the Fourier decomposition of kernels , which renders them inherently stationary. Non-stationary kernels, on the other hand, are based on generalising monotonic base kernels, such as the Matérn family of kernels , by partitioning the input space , or by input transformations .

We propose an expressive and efficient kernel family that is – in contrast to earlier methods – both non-stationary and non-monotonic, and hence can infer long-range or periodic relations in an input-dependent manner. We derive the kernel from first principles by solving the more expressive generalised Fourier decomposition of non-stationary functions, than the more limited standard Fourier decomposition exploited by earlier works. We propose and solve the generalised spectral density as a mixture of Gaussian process density surfaces that model flexible input-dependent frequency patterns. The kernel reduces to a stationary kernel with appropriate parameterisation. We show the expressivity of the kernel with experiments on time series data, image-based pattern recognition and extrapolation, and on climate data modelling.

Non-stationary spectral kernels

This section introduces the main contributions. We employ the generalised spectral decomposition of non-stationary functions and derive a practical and efficient family of kernels based on non-stationary spectral components. Our approach relies on associating input-dependent frequencies for data inputs, and solving a kernel through the generalised spectral transform.

where $\mu_{S}$ is a Lebesgue-Stieltjes measure associated to some positive semi-definite (PSD) spectral density function $S(s,s^{\prime})$ with bounded variations , which we denote as the spectral surface since it considers the amplitude of frequency pairs (See Figure 1a).

The generalised Fourier transform (27) specifies that a spectral surface $S(s,s^{\prime})$ generates a PSD kernel $K(x,x^{\prime})$ that is non-stationary unless the spectral measure mass is concentrated only on the diagonal $s=s^{\prime}$. We design a practical, efficient and flexible parameterisation of spectral surfaces that, in turn, specifies novel non-stationary kernels with input-dependent characteristics and potentially long-range non-monotonic correlation structures.

Next, we introduce spectral kernels that remove the restriction of stationarity of earlier works. We start by modeling the spectral density as a mixture of $Q$ bivariate Gaussian components

with parameterization using the correlation $\rho_{i}$, means $\mu_{i},\mu_{i}^{\prime}$ and variances $\sigma_{i}^{2},{\sigma_{i}^{\prime}}^{2}$. To produce a PSD spectral density $S_{i}$ as required by equation (27) we need to include symmetries $S_{i}(s,s^{\prime})=S_{i}(s^{\prime},s)$ and sufficient diagonal components $S_{i}(s,s)$, $S_{i}(s^{\prime},s^{\prime})$. To additionally result in a real-valued kernel, symmetry is required with respect to the negative frequencies as well, i.e., $S_{i}(s,s^{\prime})=S_{i}(-s,-s^{\prime})$. The sum $\sum_{\boldsymbol{\mu}_{i}\in\pm\{\mu_{i},\mu_{i}^{\prime}\}^{2}}$ satisfies all three requirements by iterating over the four permutations of $\{\mu_{i},\mu_{i}^{\prime}\}^{2}$ and the opposite signs $(-\mu_{i},-\mu_{i}^{\prime})$, resulting in eight components (see Figure 1a).

The generalised Fourier transform (27) can be solved in closed form for a weighted spectral surface mixture $S(s,s^{\prime})=\sum_{i=1}^{Q}w_{i}^{2}S_{i}(s,s^{\prime})$ using Gaussian integral identities (see the appendix):

We immediately notice that the BSM kernel vanishes rapidly outside the origin $(x,x^{\prime})=(0,0)$. We would require a huge number of components centered at different points $x_{i}$ to cover a reasonably-sized input space.

2 Generalised Spectral Mixture (GSM) kernel

To overcome the deficiencies of the kernel derived in Section 2.1, we extend it further by parameterizing the frequencies, length-scales and mixture weights as a Gaussian processesSee the appendix for a tutorial on Gaussian processes., that form a smooth spectrogram (See Figure 2l):

We accommodate the input-dependent lengthscale by replacing the exponential part of (3) by the Gibbs kernel

which is a non-stationary generalisation of the Gaussian kernel . We propose a non-stationary generalised spectral mixture (GSM) kernel with a simple closed form (see the appendix):

We have presented the proposed kernel (7) for univariate inputs for simplicity. The kernel can be extended to multivariate inputs in a straightforward manner using the generalised Fourier transform with vector-valued inputs . However, since in many applications multivariate inputs have a grid-like structure, for instance in geostatistics, image analysis and temporal models. We exploit this assumption and propose a multivariate extension that assumes the inputs to decompose across input dimensions :

Inference

with a standard predictive GP posterior $f(\mathbf{x}_{\star}|\mathbf{y})$ for a new input point $\mathbf{x}_{\star}$ . The posterior can be efficiently computed using Kronecker identities (see the appendix).

Related Work

Bochner’s theorem for stationary signals, whose covariance can be written as $k(\tau)=k(x-x^{\prime})=k(x,x^{\prime})$, implies a Fourier dual

The dual is a special case of the more general Fourier transform (27), and has been exploited to design rich, yet stationary kernel representations and used for large-scale inference . Lazaro-Gredilla et al. proposed to directly learn the spectral density as a mixture of Dirac delta functions leading to a sparse spectrum (SS) kernel $k_{\text{SS}}(\tau)=\frac{1}{Q}\sum_{i=1}^{Q}\cos(2\pi s_{i}^{T}\tau)$ . Wilson et al. derived a stationary spectral mixture (SM) kernel by modelling the univariate spectral density using a mixture of normals $S_{\text{SM}}(s)=\sum_{i}w_{i}[\mathcal{N}(s|\mu_{i},\sigma_{i}^{2})+\mathcal{N}(s|-\mu_{i},\sigma_{i}^{2})]/2$ , corresponding to the kernel function $k_{\text{SM}}(\tau)=\sum_{i}w_{i}\exp(-2\pi^{2}\sigma_{i}^{2}\tau)\cos(2\pi\mu_{i}\tau)$, which we generalized to the non-stationary case. Kernels derived from the spectral representation are particularly well suited to encoding long-range, non-monotonic or periodic kernels; however, they have so far been unable to handle non-stationarity.

Non-stationary kernels, on the other hand, have been constructed by non-stationary extensions of Matérn and Gaussian kernels with input-dependent lengthscales , input space warpings , and with local stationarity with products of stationary and non-stationary kernels . The simplest non-stationary kernel is arguably the dot product kernel , which has been used as a way to assign input-dependent signal variances . Non-stationary kernels are a good match for functions with transitions in their dynamics, yet are unsuitable for modelling non-monotonic properties.

Our work can also be seen as a generalisation of wavelets, or time-dependent frequency components, into general and smooth input-dependent components. In signal processing, Hilbert-Huang transforms and Hilbert spectral analysis explore input-dependent frequencies, but with deterministic transform functions on the inputs .

Experiments

We apply our proposed kernel first on simple simulated time series, then on texture images and lastly on a land surface temperature dataset. With the image data, we compare our method to two stationary mixture kernels, specifically the spectral mixture (SM) and sparse spectrum (SS) kernels , and the standard squared exponential (SE) kernel. We employ the GPML Matlab toolbox, which directly implements the SM and SE kernels, and the SS kernel as a meta kernel combining simple cosine kernels. The GPML toolbox also implements Kronecker inference automatically for these kernels. We implemented the proposed GSM kernel and inference in Matlab.

For optimizing the log posterior (60) we employ the L-BFGS algorithm. For both our method and the comparisons, we restart the optimization from 10 different initialisations, each of which is chosen as the best among 100 randomly sampled hyperparameter values as evaluating the log posterior is cheap compared to evaluating gradients or running the full optimisation.

2 Image data

We applied our kernel to two texture images. The first image of a sheet of metal represents a mostly stationary periodic pattern. The second, a wood texture, represents an example of a very non-stationary pattern, especially on the horizontal axis. We use majority of the image as training data (the non-masked regions of Figure 3a and 3f) , and use the compared kernels to predict a missing cross-section in the middle, and also to extrapolate outside the borders of the original image.

Figure 4 shows the two texture images, and extrapolation predictions given by the proposed GSM kernel, with a comparison to the spectral mixture (SM), sparse spectrum (SS) and standard squared exponential (SE) kernels. For GSM, SM and SS we used $Q=5$ mixture components for the metal texture, and $Q=10$ components for the more complex wood texture.

The GSM kernel gives the most pleasing result visually, and fills in both patterns well with consistent external extrapolation as well. The stationary SM kernel does capture the cross-section, but has trouble extrapolation outside the borders. The SS kernel fails to represent even the training data, it lacks any smoothness in the frequency space. The gaussian kernel extrapolates poorly.

3 Spatio-Temporal Analysis of Land Surface Temperatures

NASAhttps://neo.sci.gsfc.nasa.gov/view.php?datasetId=MOD11C1_M_LSTDA provides a land surface temperature dataset that we used to demonstrate our kernel in analysis of spatio-temporal data. Our primary objective is to demonstrate the capability of the kernel in inferring long-range, non-stationary spatial and temporal covariances.

We took a subset of four years (February 2000 to February 2004) of North American land temperatures for training data. In total we get 407,232 data points, constituting 48 monthly temperature measurements on a $84\times 101$ map grid. The grid also contains water regions, which we imputed with the mean temperature of each month. We experimented with the data by learning a generalized spectral mixture kernel using $Q=5$ components.

Figure 5 presents our results. Figure 5b highlights the training data and model fits for a winter and summer month, respectively. Figure 5a shows the non-stationary kernel slices at two locations across both latitude and longitude, as well as indicating that the spatial covariances are remarkably non-symmetric. Figure 5c indicates five months of successive training data followed by three months of test data predictions.

Discussion

In this paper we have introduced non-stationary spectral mixture kernels, with treatment based on the generalised Fourier transform of non-stationary functions. We first derived the bivariate spectral mixture (BSM) kernel as a mixture of non-stationary spectral components. However, we argue it has only limited practical use due to requiring an impractical amount of components to cover any sufficiently sized input space. The main contribution of the paper is the generalised spectral mixture (GSM) kernel with input-dependent Gaussian process frequency surfaces. The Gaussian process components can cover non-trivial input spaces with just a few interpretable components. The GSM kernel is a flexible, practical and efficient kernel that can learn both local and global correlations across the input domains in an input-dependent manner. We highlighted the capability of the kernel to find interesting patterns in the data by applying it on climate data where it is highly unrealistic to assume the same (stationary) covariance pattern for every spatial location irrespective of spatial structures.

Even though the proposed kernel is motivated by the generalised Fourier transform, the solution to its spectral surface

remains unknown due to having multiple GP functions inside the integral. Figure 2h highlights a numerical integration of the surface equation (12) on an example GP frequency surface. Furthermore, the theoretical work of Kom Samo and Roberts on generalised spectral transforms suggests that the GSM kernel may also be dense in the family of non-stationary kernels, that is, to reproduce arbitrary non-stationary kernels.

References

Appendix A A tutorial on Gaussian processes

We summarise here Gaussian process regression for completeness. For an interested reader, we refer to the excellent and comprehensive book by Rasmussen and Williams .

Gaussian processes (GP) are a Bayesian nonparameteric machine learning framework for regression, classification and unsupervised learning . A Gaussian process is a collection of random variables, any finite combination of which has a Multivariate normal distribution. A GP prior defines a distribution over functions, denoted as

Furthermore, the GP prior determines that for any finite collection of input points $x_{1},\ldots,x_{N}$, the corresponding function values follow a Multivariate normal distribution

Assume a dataset $\mathcal{D}=(x_{i},y_{i})_{i=1}^{N}$ and an additive Gaussian likelihood

where $K(x_{*},X)=K(X,x_{*})^{T}$ is a row kernel.

has a closed form as well. The marginal log likelihood is related to the amount of functions compatible with the prior and matching the data. Hence, the marginal log likelihood automatically promotes priors that induce functions matching the data while penalising model complexity. The marginal log likelihood can be directly maximised using standard gradient ascent techniques to infer optimal hyperparameters $\theta$.

Appendix B Deriving the bivariate spectral mixture kernel

where $\mu_{S}$ is a Lebesgue-Stieltjes measure associated to some positive semi-definite (PSD) spectral density function $S(s,s^{\prime})$ with bounded variations, which we denote as the spectral surface since it considers the amplitude of frequency pairs.

We define a spectral density $S(s,s^{\prime})$ as a mixture of $Q$ bivariate Gaussian components

with parameterization using the correlation $\rho_{i}$, means $\mu_{i},\mu_{i}^{\prime}$ and variances $\sigma_{i}^{2},{\sigma_{i}^{\prime}}^{2}$. To ensure the PSD property of spectral density $S_{i}(s,s^{\prime})$ it must hold that $S_{i}(s,s^{\prime})=S_{i}(s^{\prime},s)$ and sufficient diagonal components $S_{i}(s,s)$, $S_{i}(s^{\prime},s^{\prime})$ exist. In addition to retrieve a real-valued kernel we require symmetry with respect to the negative frequencies as well, i.e. $S_{i}(s,s^{\prime})=S_{i}(-s,-s^{\prime})$. The sum $\sum_{\boldsymbol{\mu}_{i}\in\pm\{\mu_{i},\mu_{i}^{\prime}\}^{2}}$ satisfies all three requirements by iterating over four permutations of $\{\mu_{i},\mu_{i}^{\prime}\}^{2}$ and the opposite signs $(-\mu_{i},-\mu_{i}^{\prime})$, resulting in eight components

The full $Q$-component spectral density is

Next, we compute the generalised Fourier transform in closed form by exploiting Gaussian integral identities

The $i$’th component of the kernel mixture is then

where the complex part cancels out. Now by defining a function

we can express the sum of the 8 exponentials in (37) as $\Psi_{\mu,\mu^{\prime}}(x)^{T}\Psi_{\mu,\mu^{\prime}}(x^{\prime})$. The final kernel thus takes the form

where we introduced mixture weights $w_{i}$ for each component.

Now, we immediately notice that the kernel vanishes rapidly outside the origin $(x,x^{\prime})=(0,0)$; we would require a huge number of components centered at different points $x_{i}$ to cover a reasonably-sized input space. One simple fix would be to change the exponential part to e.g. a Gaussian kernel $\exp(-\frac{1}{2}\sigma^{2}||x-x^{\prime}||^{2})$ to prevent the component from vanishing but this still would not allow us to account for non-stationary frequencies, which is what we address next.

Appendix C Deriving the generalised spectral mixture (GSM) kernel

The generalised spectral mixture kernel defines Gaussian process frequencies, lengthscales and mixture weights:

Frequency parameter $\operatorname{logit}\mu(x)$ to limit the learned frequencies between zero and the Nyquist frequency $F_{N}$, which can be defined as half of the sampling rate of the signal (or for non-equispaced signals as the inverse of the smallest time interval between the samples).

To accommodate lengthscale functions we replace the exponential part of the BSM kernel by the Gibbs kernel

The cosine part (38) is replaced by a function

The non-stationary generalised spectral mixture (GSM) kernel has a closed form

due to identity $\cos\alpha\cos\beta+\sin\alpha\sin\beta=\cos(\alpha-\beta)$. The kernel is a product of three kernels, namely a linear kernel, a Gibbs kernel and a novel cosine kernel with a feature mapping $\Psi_{i}(x)$. The full kernel is PSD due to all of its product kernels being PSD. The cosine kernel is PSD due to a dot product.

We show that the proposed non-stationary GSM kernel reduces to the stationary SM kernel with appropriate parameterisation. We show this identity for univariate inputs for simplicity, with the same result being straightforward to derive for multivariate kernel variants as well.

The proposed generalised spectral mixture (GSM) kernel for univariate inputs is

where the parameters are the weights $w_{i}$, mean frequencies $\mu_{i}$ and variances $\sigma_{i}^{2}$. Now if we assign the following constant functions for the GSM kernel to match the parameters of the SM kernel on the right-hand side,

Appendix D Inference

In many applications multivariate inputs have a grid-like structure, for instance in geostatistics, image analysis and temporal models. We exploit this assumption and propose a multivariate extension that assumes the inputs to decompose across input dimensions :

with a standard predictive GP posterior $f(\mathbf{x}_{\star}|\mathbf{y})$ for a new input point $\mathbf{x}_{\star}$ . The posterior can be efficiently computed using Kronecker identities .

The marginal likelihood (60) can be evaluated using the eigen decomposition $\mathbf{K}=\boldsymbol{Q}\boldsymbol{V}\boldsymbol{Q}^{T}$. Using known results for Kronecker products we can compute the eigen decomposition as $\boldsymbol{Q}=\bigotimes_{p}{Q}_{p}$, $\boldsymbol{V}=\bigotimes_{p}{V}_{p}$ and $\boldsymbol{Q}^{T}=\bigotimes_{p}{Q}_{p}^{T}$ using the decompositions of the smaller kernels $K_{p}={Q}_{p}{V}_{p}{Q}_{p}^{T}$. Thus we can decompose the computation of the first term in (60) as

where the inversion is taken only of the diagonal matrix of eigenvalues and matrix-vector products with a Kronecker matrix can be computed efficiently. The second term of (60) can be computed using the eigenvalues $\boldsymbol{\lambda}=\operatorname{diag}(\boldsymbol{V})=\bigotimes_{p}\operatorname{diag}({V}_{p})$ as $\log|\mathbf{K}+\sigma_{n}^{2}\mathbf{I}|=\sum_{i}\log(\lambda_{i}+\sigma_{n}^{2})$.

The gradient of the marginal likelihood is given by

where $\boldsymbol{\alpha}=(\mathbf{K}+\sigma_{n}^{2}\mathbf{I})^{-1}\mathbf{y}$ is computed as in (62). The gradient of the Kronecker product kernel can be computed as

assuming that $\frac{\partial\mathbf{K}_{p}}{\partial\theta_{i}}=\boldsymbol{0}$ for $i\neq p$. As this is a Kronecker product, the first term in (63) can be computed efficiently. The trace term in (63) can be computed by exploiting the cyclic property and the eigen decomposition as

where the latter term can be computed efficiently as

and its diagonal as a Kronecker product of the diagonals of each factor in the product. For the noise parameter $\sigma_{n}$ we get $\frac{\partial(\mathbf{K}+\sigma_{n}^{2}\mathbf{I})}{\partial\log\sigma_{n}}=2\sigma_{n}^{2}\mathbf{I}$ which makes both terms in (63) easy to compute.

Kronecker methods are also easily extensible for non-complete grids and non-Gaussian likelihoods .