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 loglinear 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.,
(1.1) 
where and where is a bounded spatial domain in dimensions. Given sampling points , our task is to sample the multivariate Gaussian , where
(1.2) 
Since is symmetric positive semidefinite, this can in principle be done by performing a factorisation
(1.3) 
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 loglinear complexity) leading to the spectral decomposition , with diagonal and containing the eigenvalues of . Provided these eigenvalues are nonnegative, 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
(1.4) 
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 cutoff function with the property that
(1.5) 
Then we define on as the infinite sum of shifts of :
(1.6) 
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 nonsmoothness 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
(1.9) 
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 gridindependent 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
(1.10) 
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 nearoptimal 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 QuasiMonte 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 highdimensional 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 wellknown (see, e.g.
[7, Corollary 5], [3, eq.(64)]) that the exact KL eigenvalues of on decay with the rate(1.11) 
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
(1.12) 
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, wavelettype representations
are constructed, with summation over and , and where
have the same multileveltype 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
[2].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
(2.1) 
When is periodic in each coordinate direction and then can be represented as its Fourier series:
(2.2) 
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 :
(2.3) 
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
(2.4) 
yields an approximation of the Fourier coefficients for the same values of . The approximation error can be quantified uniformly in by the following wellknown 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
(2.5) 
and in particular, for any ,
(2.6) 
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
(2.8) 
For its Fourier transform, we have (see, e.g., [9, eq. (2.22)]
where
(2.9) 
The modified Bessel functions of the second kind have (see [1, 9.6.24]) the integral representations
(2.10) 
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
(2.11) 
Lemma 3.
Let . Then is of the form
(2.12) 
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 nonsmooth 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:
(3.1) 
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) above^{1}^{1}1Note 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
(3.2) 
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.
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
(3.4) 
It therefore remains to obtain suitable estimates for the Fourier coefficients of the periodised covariance. To this end, we use a cutoff function to separate the artificial nonsmooth part created by the periodization. We thus define an even smooth univariate cutoff function by requiring that is supported on and
For any , we can scale this to a cutoff function supported on by defining
Using this function, we now write
(3.5) 
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
Proof.
We have , since vanishes outside . By the convolution theorem,
(3.6) 
The second term on the right can be estimated as
with constants independent of , where we have used
(3.7) 
and the decay properties of for Matérn covariances.
Lemma 8.
Let be a Matérn covariance function with . Then there exists independent of such that
Proof.
When we use the estimate
(3.8) 
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
(3.9)  
Denoting
we get . If , integrating by parts with respect to then gives
(3.10)  
We use the following auxiliary statement, which is proved in the appendix.
Lemma 9.
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.
Comments
There are no comments yet.