Unified Analysis of Periodization-Based Sampling Methods for Matérn Covariances

05/31/2019 ∙ by Markus Bachmayr, et al. ∙ University of Bath University of Heidelberg University of Bonn University of Mainz 0

The periodization of a stationary Gaussian random field on a sufficiently large torus comprising the spatial domain of interest is the basis of various efficient computational methods, such as the classical circulant embedding technique using the fast Fourier transform for generating samples on uniform grids. For the family of Matérn covariances with smoothness index ν and correlation length λ, we analyse the nonsmooth periodization (corresponding to classical circulant embedding) and an alternative procedure using a smooth truncation of the covariance function. We solve two open problems: the first concerning the ν-dependent asymptotic decay of eigenvalues of the resulting circulant in the nonsmooth case, the second concerning the required size in terms of ν, λ of the torus when using a smooth periodization. In doing this we arrive at a complete characterisation of the performance of these two approaches. Both our theoretical estimates and the numerical tests provided here show substantial advantages of smooth truncation.



There are no comments yet.


page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1. Introduction

The simulation of Gaussian random fields with specified covariance is a fundamental task in computational statistics with a wide range of applications – see, for example, [4, 19, 6]. Since these fields provide natural tools for spatial modeling under uncertainty, they have played a fundamental rôle in the modern field of uncertainty quantification (UQ). Thus the construction and analysis of efficient methods for sampling such fields is of widespread interest.

In this paper we consider a class of algorithms for sampling stationary fields based on (artificial) periodization of the field, and fast Fourier transform. These algorithms enjoy log-linear complexity in the number of samples computed. While they are classical in computational statistics e.g., [4, 19] and are widely applied in UQ [8, 15, 9, 3, 2], so far they have been subjected to relatively little rigorous numerical analysis. This paper provides a detailed novel analysis of these algorithms in the case of the important Matérn covariance, providing new results related to their efficiency and their use in obtaining separable expansions of the random fields.

We consider Gaussian random fields which are assumed stationary, i.e.,


where and where is a bounded spatial domain in dimensions. Given sampling points , our task is to sample the multivariate Gaussian , where


Since is symmetric positive semidefinite, this can in principle be done by performing a factorisation


with being one possible choice, from which the desired samples are provided by the product where . When is large a naive direct factorisation is prohibitively expensive. A variety of methods for drawing samples have been developed that either perform approximate sampling, for instance by using a less costly inexact factorisation (e.g., [11, 5]), or maintain exactness of the distribution by using an exact factorisation that exploits additional structure in . The methods considered in this work follow the latter approach by embedding into a problem with periodic structure.

Fast methods for sampling are based on the observation that in many applications the sampling points can be placed uniformly on a rectangular grid, as justified theoretically, e.g., in [10]. Then the symmetric positive semidefinite matrix is block Toeplitz with Toeplitz blocks. Via periodization, can be embedded into an block circulant matrix with circulant blocks (and ). This can be diagonalised using FFT (with log-linear complexity) leading to the spectral decomposition , with diagonal and containing the eigenvalues of . Provided these eigenvalues are non-negative, the required in (1.3) can be computed by taking appropriate rows of the square root of . Then samples of the grid values of can then be drawn by computing


with the columns of (suitably truncated) and the corresponding eigenvalues of . Note that (1.4

) is the Karhunen–Loève (KL) expansion of the random vector


Two key questions of numerical analysis arise concerning this procedure. The first is to describe conditions which guarantee the positive semidefiniteness of the extension , since without this property the formula (1.4) will fail to produce samples of . The second is to determine the rate of decay of the eigenvalues as (and hence ) increases, This rate of decay is crucial in understanding the asymptotic properties of the KL expansion (1.4) and in particular in determining the convergence properties of many methods for uncertainty quantification (see, e.g., [8]). In this paper these questions are answered in detail in the case where the covariance is of Matérn type and in the context of two different methods of periodization, both in use in practice, partial answers having previously been obtained in the recent papers [3, 9].

Before describing our main results, we briefly describe the periodization procedure in the case of physical dimension ; higher dimensions are entirely analogous. Since the sampling domain is bounded, without loss of generality we assume that it is contained in , and so the covariance is only evaluated on the domain . For any we construct a -periodic extension of as follows: First choose a cut-off function with the property that


Then we define on as the infinite sum of shifts of :


Clearly is -periodic. Moreover, when and , we have , and (1.6) collapses to a single term, yielding on . We shall discuss in detail two examples corresponding to different choices of :

(1.7) classical periodization:
(1.8) smooth periodization:

In (1.7), the periodization is obtained by simply repeating the function periodically. It is simple to implement but has the disadvantage that articial non-smoothness is introduced at the points . By contrast, in the smooth periodization (1.8), the function , has the same smoothness properties on as but is supported on . Thus the summands in (1.6) do not have disjoint supports. For example, the supports of and overlap on ).

Returning to domains of general dimension , to guarantee positive definiteness of the circulant , needs to be chosen sufficiently large. For the Matérn family of covariances defined in (2.8), for classical periodization, [9] gave a sufficient condition on of the form


when , where is the Matérn smoothness parameter, is the length scale and is the sampling rate, i.e., the spacing between the sampling points . This bound was seen to be essentially sharp in numerical experiments in [9]; in particular, the required indeed diverges as , and one thus obtains a positive definite covariance matrix on each finite grid, but no underlying periodic random field.

By contrast, for smooth periodization, it was shown in [3] for Matérn and various related covariances that by taking sufficiently large, one can always obtain a positive definite periodized covariance function, and hence a grid-independent periodic random field. However, the required size of was not quantified. In one of our main results (proved in Theorem 10), we show that for Matérn covariances, there exists such that for any , it suffices to take


We thus obtain essentially the same qualitative behaviour with respect to , as in (1.9), but without the dependence on the grid size , and including the regime . This means that in terms of , whereas in the classical case, one sample on costs , in the smooth case the cost is only (with the remaining logarithmic factor stemming from the costs of the FFT). Our numerical tests below demonstrate that, especially for , this indeed translates to substantial gains in quantitative efficiency.

The condition (1.10) for all can also be regarded as a generalisation of the previous contributions [17, 6, 16] where periodization was considered with smooth truncations specifically designed for particular types of covariances, including the Matérn case, but with the corresponding analysis limited to certain ranges of . Methods using smooth cutoff functions similar to the ones analysed here have also been tested computationally in [13, 12].

One further useful consequence of the new result (1.10) covering the full range of , is that smooth periodization also provides an attractive way of sampling with hyperpriors on these two parameters: in this case, one needs to first randomly select and

according to some probability distribution, and then draw a sample of the Gaussian random field with the corresponding Matérn covariance. This problem has been addressed, for instance, in

[14] by approximate sampling using a reduced basis. In this context, it is relevant that for both types of periodization, setting up the factorisation for a new covariance and drawing a sample come at the same cost of one FFT each: thus by (1.10), provided that the distributions of , have compact supports inside , smooth periodization is guaranteed to provide each sample – without further approximation – at near-optimal cost .

Our second main result concerns the decay of the ordered eigenvalues of the circulant matrix . This question is particularly relevant in several areas of UQ. For example, in the analysis of Quasi-Monte Carlo integration methods, the rate of decay of the terms in the discrete KL expansion (1.4) plays a key rôle in the convergence theory of QMC for high-dimensional problems. This rate needs to be compared to the KL expansion of on ,

where are the

-orthonormal eigenfunctions of the covariance operator. In the Matérn case it is well-known (see, e.g.

[7, Corollary 5], [3, eq.(64)]) that the exact KL eigenvalues of on decay with the rate


Due to the preservation of covariance regularity enjoyed by the smooth periodization, the ordered eigenvalues of in the smooth periodization case decay at the same optimal rate as in (1.11), see [3, eq. (64)]. However, in the case of classical periodization, the associated loss of regularity means that the decay rate is not obvious. Supported by numerical evidence, it was conjectured in [9] that under the condition (1.9), the eigenvalues of in the classical case also exhibit, uniformly in , the same asymptotic decay rate (1.11). In Theorem 6 we prove this conjecture, up to a minor modification by a multiplicative factor of order .

In summary, this paper provides a complete characterisation of the performance of the two types of periodization in the case of Matérn covariances. Both lead to optimal decay of covariance matrix eigenvalues, whereas the required size parameter of the periodization cell is substantially more favorable in the smooth periodization case.

Before proceeding we make a few more remarks about the other advantages enjoyed by the smooth periodization method. The periodic random field on obtained using smooth truncation also provides a tool for deriving series expansions of the original random field on . As mentioned above the KL eigenvalues of the periodic field in the smooth periodization case have the decay rate (1.11). Moreover, in contrast to the KL eigenfunctions on , which are typically not explicitly known, the corresponding eigenfunctions of the periodic covariance are explicitly known trigonometric functions, and, restricting this back to , one obtains an expansion


of the original random field on . Thus no approximate computation of eigenfunctions is required, in contrast to the the standard KL expansion in , and the eigenvalues can again be approximated efficiently by FFT as described in [3, Sec. 5.1]. The asymptotic decay of is the same as that of . What is more, the scaled -norms also decay at the same rate. This decay of -norms is important for applications to random PDEs. It can actually be faster than for the standard KL expansion, where is in general not uniformly bounded, see [3, Sec. 3].

The KL expansion of also enables the construction of alternative expansions of of the basic form (1.12), but with the spatial functions having additional properties. In [3], applying the square root of the covariance operator in the corresponding factorisation to periodic Meyer wavelets, wavelet-type representations

are constructed, with summation over and , and where

have the same multilevel-type localisation as the Meyer wavelets. This feature yields improved convergence estimates for tensor Hermite polynomial approximations of solutions of random diffusion equations with lognormal coefficients


The outline of the paper is as follows: in Section 2, we introduce some notions and basic results that are relevant to both the classical and smooth periodizations; in Section 3, we prove the conjecture from [9] (slightly modified) on the rate of eigenvalue decay for in the classical case; in Section 4, we establish the quantitative condition (1.10) for a smooth truncation; and in Section 5, we illustrate our findings by numerical tests.

We use the following notational conventions: vectorial quantities are denoted by boldface letters; is the Euclidean norm of ; is the Euclidean ball of radius with center . We use as a generic constant whose definition can change whenever it is used in a new inequality.

2. Preliminaries

2.1. Fourier transforms

For a suitably regular function , the Fourier transform on and its inverse are defined for by


When is periodic in each coordinate direction and then can be represented as its Fourier series:


for all , with . Moreover, if belongs to a Hölder space for some then the sum in (2.2) converges uniformly.

Let be an even integer, set and introduce the infinite uniform grid of points on :


By restricting to lie in

we obtain a uniform grid on . The (appropriately scaled) discrete Fourier transform of the corresponding grid values of , given by


yields an approximation of the Fourier coefficients for the same values of . The approximation error can be quantified uniformly in by the following well-known result on the error of the trapezoidal rule applied to periodic functions, whose proof we include for convenience of the reader.

Lemma 1.

Let be periodic in each coordinate direction and for some . Then


and in particular, for any ,


Using (2.2) and uniform convergence of the Fourier series by Hölder continuity of , we have



The last term vanishes unless for each , in which case it takes the value . Since , we have

Now (2.5) is obtained by noting that .

For fixed , we introduce the function , for . Then and . From (2.5) we conclude that

Now (2.6) follows since . ∎

2.2. Covariance functions

On a computational domain , we consider the fast evaluation of a Gaussian random field with covariance function given by (1.1). We focus on the Matérn covariance kernel with correlation length and smoothness exponent , which is given by


For its Fourier transform, we have (see, e.g., [9, eq. (2.22)]



The modified Bessel functions of the second kind have (see [1, 9.6.24]) the integral representations


which shows in particular that their values for fixed are monotonically increasing with respect to for , and also that . We will also use the following results that directly imply exponential decay of as . Their proofs are given in Appendix A.

Lemma 2.

Let and . Then we have

Lemma 3.

Let . Then is of the form


with coefficients satisfying

Inspection of (2.8) shows that changing amounts to a rescaling of the computational domain . Hence without loss of generality, we may assume our computational domain to be contained in the box , so that any difference of two points in is contained in . In what follows, we subsequently embed this box into a torus with , on which we define a -periodic covariance such that for , which means that the covariance between any pair of points in is preserved in replacing by .

3. Classical Periodization

In this section, we prove that for Matérn covariances, the asymptotic decay of eigenvalues is preserved in the covariance matrix obtained by classical circulant embedding with non-smooth periodization as in (1.7). This confirms a conjecture made recently in [9].

We consider the extension defined as in (1.7). Then the matrix is defined by:


For any even positive integer , is a circulant extension of the covariance matrix given in (1.2), obtained when sampling on a uniform grid on any cube . More general -dimensional rectangles can be treated in the same manner with the obvious modifications. If the index set is given lexicographical ordering, then is a nested block Toeplitz matrix where the number of nested levels is the physical dimension and is a nested block circulant extension of it.

To analyse the eigenvalues of it is useful to also consider the continuous periodic covariance integral operator:

Then the scaled circulant matrix can be identified as a Nyström approximation of using the composite trapezoidal rule with respect to the uniform grid on given by the points (2.3) with . The operator is a compact operator on the space of -periodic continuous functions on , and so it has a discrete spectrum with the only accumulation point at the origin.

The following result is standard (see for example [9]).

Proposition 4.
  • The eigenvalues of are , as defined in (2.2), with corresponding eigenfunctions normalised in :

  • The eigenvalues of are , , as defined in (2.4

    ), with corresponding eigenvectors normalized with respect to the Euclidean norm:

It was convenient to introduce the scaling factor in (ii) above111Note that in the notation of [9], is and is ., since we can then identify (2.4) as the trapezoidal rule approximation of the (-periodic) Fourier transform defining .

Our analysis for estimating the decay of as will be based on the formula (2.6) and estimating the rate of decay of . To prepare for the main result we now restrict our consideration to the Matérn covariances, that is, for the remainder of this section we assume

with some , where is defined in (2.8). In the recent paper [9] the following theorem was proven for this case.

Theorem 5.

Let , , and . Then there exist which may depend on but are independent of , such that is positive definite if


We are now in a position to formulate our result on the rate of decay of eigenvalues of for Matérn covariances, which proves the conjecture made in [9, eq. (3.9)] up to an additional factor of order in the constant. In what follows, let us denote by the smallest value of such that the condition (3.2) holds true.

Theorem 6.

Let , , and be as in Theorem 5, and let be the decreasingly ordered eigenvalues of . Let be chosen such that with an independent of , where in (3.2) is sufficiently large with respect to . Then there exists such that for all ,


where with depending only on , , and .

In the proof of this result, we rely on the following observation. Noting that , for each , by Proposition 4 we have for some . Using Lemma 1, we thus obtain


It therefore remains to obtain suitable estimates for the Fourier coefficients of the periodised covariance. To this end, we use a cut-off function to separate the artificial nonsmooth part created by the periodization. We thus define an even smooth univariate cut-off function by requiring that is supported on and

For any , we can scale this to a cut-off function supported on  by defining

Using this function, we now write


In the following two lemmas, we separately estimate and , defined as in (2.2).

Lemma 7.

Let be a Matérn covariance function. Then for , there exists independent of such that


We have , since vanishes outside . By the convolution theorem,


The second term on the right can be estimated as

with constants independent of , where we have used


and the decay properties of for Matérn covariances.

For the first term on the right of (3.6), we use two separate bounds for small and large , that is,

where is uniformly bounded with respect to by (3.7). Moreover, since and since is smooth, for any ,

Combining these bounds completes the proof. ∎

Lemma 8.

Let be a Matérn covariance function with . Then there exists independent of such that


When we use the estimate


For the case where for at least one value of , we integrate by parts dimensionwise to best exploit the limited smoothness of across the boundary of . We assume without loss of generality that and first give the proof for to simplify notation. Integration by parts of (3.8) with respect to gives



we get . If , integrating by parts with respect to then gives


We use the following auxiliary statement, which is proved in the appendix.

Lemma 9.

Let with and with . Then, with as given in (3.5),

where is independent of  .

This gives the desired bound for the last term in (3.10). Moreover,

which is the bound for (3.8) and the first terms on right hand side of(3.9) and (3.10). Finally, for the second term on the right hand side of (3.10) we have with

Carrying this out in the same way for further other terms we obtain the desired result for . The proof for can be done analogously, giving rise to terms in (3.10). ∎

Based on these preparations, we now complete the proof of Theorem 6.

Proof of Theorem 6..

We estimate the two terms on the right hand side of (3.4). Combining Lemmas 7, with , and 8 yields