Unified Analysis of Periodization-Based Sampling Methods for Mat\'ern Covariances

TL;DR

This paper compares periodization methods for Matérn covariances, revealing that smooth truncation offers significant advantages over nonsmooth approaches in terms of eigenvalue decay and required domain size.

Contribution

It provides a complete analysis of eigenvalue decay and domain size requirements for both nonsmooth and smooth periodization methods for Matérn covariances.

Findings

Smooth truncation outperforms nonsmooth periodization in eigenvalue decay.

The required torus size depends on the smoothness index and correlation length.

Numerical results confirm theoretical advantages of smooth truncation.

Abstract

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\'ern covariances with smoothness index $\nu$ and correlation length $\lambda$, 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 $\nu$-dependent asymptotic decay of eigenvalues of the resulting circulant in the nonsmooth case, the second concerning the required size in terms of $\nu$, $\lambda$ 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.

Figures

Equations

$$\mathbb{E}[(Z(x,\cdot )-\overline{Z}(x))(Z(x^{\prime },\cdot )-\overline{Z}(x^{\prime }))]=:\rho (x-x^{\prime }),x,x^{\prime }\in D,$$

$${\Sigma }_{i,j}=\rho (x_i-x_j),i,j=1,\dots ,n.$$

$$\Sigma =F{F}^{\top },$$

$${\Sigma }^{\mathrm{ext}}=Q^{\mathrm{ext}}{\Lambda }^{\mathrm{ext}}(Q^{\mathrm{ext}}{)}^{\top },$$

$$V^{\text{ext}}=\sum _{j=1}^m{y}_j\sqrt{{\Lambda }_j^{\mathrm{ext}}}{q}_j$$

$$Z=\sum _{j=1}^{\infty }{y}_j\sqrt{{\lambda }_j}{\phi }_j,y_j\sim \mathcal{N}(0,1)\text{ i.i.d.,}$$

$$\phi =1\text{on}[-1,1],\text{and}\phi =0\text{on}\mathbb{R}\backslash [-\kappa ,\kappa ],\text{where}\kappa :=2\gamma -1\ge \gamma .$$

$$\displaystyle\rho^{\mathrm{ext}}(x)=\sum_{n\in\mathbb{Z}}\bigl{(}\rho\varphi\bigr{)}(x+2\gamma n),\quad x\in\mathbb{R}\,.$$

$$\phi =\text{raise}1.29167pt{\chi }_{[-\gamma ,\gamma )},\text{the characteristic function of [−γ,γ);}$$

$$\gamma >1\text{and}\phi \in C_0^{\infty }(\mathbb{R}),\text{with }\mathrm{supp}(\phi )=[-\kappa ,\kappa ].$$

$$\gamma\geq C\bigl{(}1+\lambda\max\{\sqrt{\nu}(1+\lvert\log\nu\rvert),1/\sqrt{\nu}\}\bigr{)},$$

$$\gamma\geq C\bigl{(}1+\lambda\max\{\sqrt{\nu}(1+\lvert\log\nu\rvert),\sqrt{\nu}\,|\!\log(h/\lambda)|\}\bigr{)}$$

$${\lambda }_j\le C{j}^{-(1+2\nu /d)}.$$

$$\hat{v}(\omega )={\int }_{{\mathbb{R}}^d}\exp (-\mathrm{i}\omega \cdot x)v(x)\mathrm{d}x,\text{and}v(x)=\frac{1}{(2\pi {)}^d}{\int }_{{\mathbb{R}}^d}\exp (\mathrm{i}\omega \cdot x)\hat{v}(\omega )\mathrm{d}\omega .$$

$$f(x)=(2\gamma {)}^{-d}\sum _{k\in {\mathbb{Z}}^d}{\hat{f}}_k\exp (\mathrm{i}{\omega }_k\cdot x),\text{where}{\hat{f}}_k={\int }_{[-\gamma ,\gamma {]}^d}f(x)\exp \left(-\mathrm{i}{\omega }_k\cdot x\right)\mathrm{d}x,$$

$$x_n:=nh\text{for}n\in {\mathbb{Z}}^d.$$

$$\displaystyle{\overline{\mathbb{Z}}_{N}^{d}:=\big{\{}-N/2,\ldots,N/2-1\big{\}}^{d}},$$

$$(S_{N}f{)}_k:=h^d\sum _{n\in {\overline{\mathbb{Z}}}_N^d}f(x_n)\exp (-\mathrm{i}{\omega }_k\cdot x_n),k\in {\overline{\mathbb{Z}}}_N^d,$$

$$h^d\sum _{n\in {\overline{\mathbb{Z}}}_N^d}f(x_n)-{\int }_{[-\gamma ,\gamma {]}^d}f(x)\mathrm{d}x=\sum _{0\ne m\in {\mathbb{Z}}^d}{\hat{f}}_{mN},$$

$$(S_{N}f{)}_k-{\hat{f}}_k=\sum _{0\ne m\in {\mathbb{Z}}^d}{\hat{f}}_{k+mN}.$$

$$Q(f):=h^{d}\sum_{n\in\overline{\mathbb{Z}}_{N}^{d}}f(x_{n})=\,(2\gamma)^{-d}\sum_{m\in\mathbb{Z}^{d}}\widehat{f}_{m}\Bigg{(}h^{d}\sum_{n\in\overline{\mathbb{Z}}_{N}^{d}}\exp(\mathrm{i}\omega_{m}\cdot x_{n})\Bigg{)}.$$

$$\sum _{n\in {\overline{\mathbb{Z}}}_N^d}\exp (\mathrm{i}{\omega }_m\cdot x_n)$$

$$=\sum _{n\in {\overline{\mathbb{Z}}}_N^d}\prod _{j=1}^d\exp \left(\mathrm{i}\frac{2\pi }{N}{m}_j{n}_j\right)=\prod _{j=1}^d\sum _{n\in {\overline{\mathbb{Z}}}_N}\exp \left(\mathrm{i}\frac{2\pi }{N}{m}_{j}n\right).$$

$$(S_{N}f{)}_k-{\hat{f}}_k=\sum _{0\ne m\in {\mathbb{Z}}^d}{\hat{g}}_{mN}.$$

$$\rho(x):=\rho_{\lambda,\nu}(x):=\frac{2^{1-\nu}}{\Gamma(\nu)}\bigg{(}\frac{\sqrt{2\nu}|x|}{\lambda}\bigg{)}^{\nu}K_{\nu}\bigg{(}\frac{\sqrt{2\nu}|x|}{\lambda}\bigg{)}\,.$$

$$\displaystyle\widehat{\rho}_{\lambda,\nu}(\omega):=\int_{\mathbb{R}^{d}}\rho_{\lambda,\nu}(x)\,\exp(-\mathrm{i}\omega\cdot x)\,\mathrm{d}x=C_{\lambda,\nu}\bigg{(}\frac{2\nu}{\lambda^{2}}+|\omega|^{2}\bigg{)}^{-(\nu+d/2)},$$

$$C_{\lambda ,\nu }:=(2\sqrt{\pi }{)}^d\frac{\Gamma (\nu +d/2)(2\nu {)}^{\nu }}{\Gamma (\nu ){\lambda }^{2\nu }}.$$

$$K_{\nu }(t)={\int }_0^{\infty }{e}^{-t\cosh (s)}\cosh (\nu s)\mathrm{d}s,$$

$$K_{\nu }(t)\le e\frac{2^{2\nu }\Gamma (\nu )}{2\sqrt{2t}}{e}^{-t}.$$

$$\frac{\mathrm{d}^{n}}{\mathrm{d}t^{n}}\big{(}t^{\nu}K_{\nu}(t)\big{)}=\sum_{j=0}^{\lfloor{n}/{2}\rfloor}a_{n,j}t^{\nu-j}K_{\nu-n+j}(t)\,,$$

$$\sum _{j=0}^{\lfloor n/2\rfloor }|a_{n,j}|\le n!.$$

Full text

Keywords. stationary Gaussian random fields, circulant embedding, periodization, Matérn covariances

Mathematics Subject Classification. 60G15, 60G60, 42B05, 65T40

M.B. and V.K.N. acknowledge support by the Hausdorff Center of Mathematics, University of Bonn. V.K.N. would like to thank Vietnam Institute for Advanced Study in Mathematics for hospitality during his visit in summer 2018 when part of the work on this paper was undertaken. All authors thank the Johann Radon Institute for Computational and Applied Mathematics, Linz, for support to attend the Workshop on “Frontier Technologies for High-Dimensional Problems and Uncertainty Quantification”, December 2018, during which part of this work was carried out. I.G.G. and R.S. thank Dirk Nuyens (Leuven) for useful discussions

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 spatial points sampled. While these algorithms are well-known 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 family of Matérn covariances. In particular we present new results related to their efficiency and to the preservation of spectral properties in the periodization.

To give more context, we shall be concerned with Gaussian random fields $Z=Z(x,\omega)$ which are assumed to be stationary, that is,

[TABLE]

where $\overline{Z}(x)=\mathbb{E}[Z(x,\cdot)]$ and where $D$ is a bounded spatial domain in $d$ dimensions. Given $n$ sampling points $\{x_{j}\}_{j=1}^{n}$, our task is to sample the multivariate Gaussian $V:=\{Z(x_{i})\}_{i=1}^{n}\in\mathcal{N}(0,\Sigma)$, where

[TABLE]

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

[TABLE]

with $F=\Sigma^{1/2}$ being one possible choice, from which the desired samples are provided by the product $FY$ where $Y\in\mathcal{N}(0,I)$. When $n$ 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 $\Sigma$. The methods considered in this work follow the latter approach by embedding into a problem with periodic structure.

The methods we consider are based on the observation that in many applications the sampling points can be placed uniformly on a rectangular grid, as justified theoretically, for instance, in [10]. Then the $n\times n$ symmetric positive semidefinite matrix $\Sigma$ is block Toeplitz with Toeplitz blocks due to (1.2). Via periodization, $\Sigma$ is embedded into an $m\times m$ block circulant matrix $\Sigma^{\mathrm{ext}}$ with circulant blocks (with $m>n$, but $m$ typically proportional to $n$), which is diagonalized using FFT (with log-linear complexity) to provide the spectral decomposition

[TABLE]

with $\Lambda^{\mathrm{ext}}$ diagonal and containing the eigenvalues of $\Sigma^{\mathrm{ext}}$ and $Q^{\mathrm{ext}}$ being a Fourier matrix. Provided these eigenvalues are non-negative, the required $F$ in (1.3) can be computed by taking $n$ appropriate rows of the square root of $\Sigma^{\mathrm{ext}}$. Then samples of the grid values of $Z$ can be drawn by first drawing a random vector $(y_{j})_{j=1,\ldots,m}$ with $y_{j}\sim{\mathcal{N}}(0,1)$ i.i.d., then computing

[TABLE]

using the FFT, with $q_{j}$ the columns of $Q^{\mathrm{ext}}$ and $\Lambda^{\mathrm{ext}}_{j}$ the corresponding eigenvalues of $\Sigma^{\mathrm{ext}}$. Note that (1.5) is the Karhunen–Loève (KL) expansion of the random vector $V^{\text{ext}}$.

Finally, a sample of $V$ is obtained by extracting from $V^{\text{ext}}$ the entries corresponding to the original grid points. Proceeding in this manner yields exact samples of $Z$ at the grid points, provided that $\Sigma^{\mathrm{ext}}$ is positive semidefinite. Positive definiteness can be verified a posteriori by checking non-negativity of the computed entries of $\Lambda^{\mathrm{ext}}$ and is guaranteed to be satisfied for sufficiently large $m$ under fairly general conditions. For instance, it was shown in [9, Thm. 2.3] that this holds true under weak integrability and regularity assumptions on $\rho$.

It thus suffices to iteratively increase $m$ until $\Sigma^{\mathrm{ext}}$ is positive semidefinite to obtain a reliable sampling method. This paper is concerned with the following two key questions concerning this procedure.

• I

How large does the extended size $m$ need to be, compared to the cardinality $n$ of the original grid? This completely determines the efficiency of the sampling scheme. The required $m$ depends both on the covariance function $\rho$ and on the type of periodization.

• II

Do the eigenvalues $\{\Lambda^{\mathrm{ext}}_{j}:j=1,\ldots,m\}$ maintain a consistent rate of decay as $n$ (and hence $m$) increases? This determines the efficiency of numerical methods that build on the decomposition (1.5), since faster decay of the eigenvalues reduces the number of independent variables which are effectively needed to describe $Z$.

The answer to Question II 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.5) plays a key rôle in the convergence theory of QMC for high-dimensional problems [10]. The rate of decay of the $\Lambda^{\mathrm{ext}}_{j}$ (as $m$ increases) is intricately linked to the rate of decay of the exact KL eigenvalues $\lambda_{j}$ of the covariance operator with kernel $\rho$ defined on $L^{2}(D)$, which appear in the KL expansion of the continuous field $Z$,

[TABLE]

where $\varphi_{j}$ are the $L^{2}(D)$-orthonormal eigenfunctions of the covariance operator. Thus, this is a question of how well the properties of the spectrum of the original covariance operator are preserved by the periodization.

While partial answers to Questions I and II have appeared elsewhere [3, 9], this paper gives a full answer to both questions 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. Our results provide a complete quantitative characterization in terms of the parameters of the Matérn covariances. Although we focus on this class of covariances to avoid further complication of the already substantial technical difficulties, the techniques based on suitable cutoff functions that we use in our main results may be more generally applicable. In our concrete arguments, the combination of the exponential decay of $\rho$ and the algebraic decay of its Fourier transform $\hat{\rho}$ towards infinity, which holds in particular for Matérn covariances, plays an important rôle.

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

[TABLE]

Then we define $\rho^{\mathrm{ext}}$ on $\mathbb{R}$ as the infinite sum of shifts of $\rho\varphi$:

[TABLE]

Clearly, $\rho^{\mathrm{ext}}$ is $2\gamma$-periodic. Moreover, when $x\in[-1,1]$ and $0\not=n\in\mathbb{Z}$, we have $x+2n\gamma\in\mathbb{R}\backslash(-\kappa,\kappa)$, and (1.8) collapses to a single term, yielding $\rho^{\mathrm{ext}}=\rho$ on $[-1,1]$. We shall discuss in detail two examples corresponding to different choices of $\varphi$:

[TABLE]

In (1.9), the periodization is obtained by simply repeating the function $\rho|_{[-\gamma,\gamma]}$ periodically. It is easy to implement but has the disadvantage that artificial non-smoothness is introduced at the points $\{2n\gamma:0\not=n\in\mathbb{Z}\}$. By contrast, in the smooth periodization (1.10), the function $\rho\varphi$ has the same smoothness properties on $\mathbb{R}$ as $\rho$ but is supported on $[-\kappa,\kappa]$. An illustration of (1.8) for the choices (1.9) and (1.10) is given in Figure 1.

Returning to domains of general dimension $d$, we first highlight our results concerning Question I. In the case of smooth periodization, it was previously shown in [3, Thm. 2.3] for Matérn and various related (e.g., anisotropic Matérn) covariances that by taking $\gamma$ 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 $\gamma$ was not quantified in [3]. In our first main result – Theorem 10 below – we show that for Matérn covariances, for any Matérn smoothness parameter $\nu>0$ and correlation length $\lambda>0$, it suffices to take

[TABLE]

with a constant $C>0$ that is independent of $\nu$ and $\lambda$. Due to the existence of the periodic random field on the continuous level, for sampling on a discrete grid the result applies to any grid size $h$, defined as the (uniform) spacing between the sampling points $x_{i}$.

This should be compared with the corresponding result for classical periodization [9]. There, assuming in addition $\nu\geq\frac{1}{2}$, a sufficient condition for positive definiteness of the form

[TABLE]

was obtained (again with $C$ independent of $\nu$ and $\lambda$, as well as $h$). This bound was seen to be essentially sharp in numerical experiments in [9]; in particular, the required $\gamma$ indeed diverges as $h\to 0$, and one thus obtains, for any fixed $h>0$, a positive definite covariance matrix on a finite grid, but no underlying periodic random field.

For smooth periodization, we thus obtain in (1.11) essentially the same qualitative behaviour with respect to $\lambda$, $\nu$ as in (1.12), but without the dependence on the grid size $h$, and including the regime $\nu\in(0,1/2)$. The absence in (1.11) of the logarithmic term in $h$ that was present in (1.12) leads to a substantial reduction of the computational complexity of the resulting sampling scheme. Moreover, the periodic extension of the original random field $Z$ on the continuous level can be used to obtain computationally attractive series expansions of $Z$ that provide alternatives to the classical KL expansion. These further conclusions are explained in detail in Section 6.

The condition (1.11) for all $\nu>0$ can also be regarded as a generalisation of the previous works [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 $\nu<2$. Methods using smooth cutoff functions similar to the ones analysed here have also been tested computationally in [13, 12].

Our second main result concerns the answer to Question II in the case of non-smooth periodization. In the Matérn case it is well-known (see, e.g. [7, Corollary 5], [3, eq.(64)]) that the exact KL eigenvalues $\lambda_{j}$ of $Z$ in $L^{2}(D)$ decay with the rate

[TABLE]

Due to the preservation of covariance regularity enjoyed by the smooth periodization, the ordered eigenvalues of $\Sigma^{\mathrm{ext}}$ in the smooth periodization case decay at the same optimal rate as in (1.13), 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.12), the eigenvalues of $\Sigma^{\mathrm{ext}}$ also exhibit, uniformly in $h$, the same asymptotic decay rate (1.13) in the classical case. In Theorem 6, we prove this conjecture, up to a minor modification by a multiplicative factor of order $\mathcal{O}(|\!\log h|^{\nu})$.

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 $\gamma$ of the periodization cell is substantially more favorable in the smooth periodization case.

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 $\Sigma^{\mathrm{ext}}$ in the classical case; in Section 4, we establish the quantitative condition (1.11) for a smooth truncation; and in Section 5, we illustrate our findings by numerical tests. In Section 6, we conclude with a discussion of the computational implications of our findings and of further applications.

We use the following notational conventions: $\lvert x\rvert$ is the Euclidean norm of $x\in\mathbb{R}^{d}$; $B_{r}(x)$ is the Euclidean ball of radius $r>0$ with center $x$. We use $C>0$ 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 $v:\mathbb{R}^{d}\rightarrow\mathbb{C}$, the Fourier transform on $\mathbb{R}^{d}$ and its inverse are defined for $\omega,x\in\mathbb{R}^{d}$ by

[TABLE]

When $f:\mathbb{R}^{d}\rightarrow\mathbb{C}$ is $2\gamma-$periodic in each coordinate direction and f\in L_{2}\big{(}[-\gamma,\gamma]^{d}\big{)} then $f$ can be represented as its Fourier series:

[TABLE]

for all $k\in\mathbb{Z}^{d}$, with $\omega_{k}:=\pi k/\gamma$. Moreover, if $f$ belongs to a Hölder space C^{0,\alpha}\big{(}[-\gamma,\gamma]^{d}\big{)} for some $\alpha>0$ then the sum in (2.2) converges uniformly.

Let $N\geq 2$ be an even integer, set $h=2\gamma/N$ and introduce the infinite uniform grid of points on $\mathbb{R}^{d}$:

[TABLE]

By restricting $n$ to lie in

[TABLE]

we obtain a uniform grid on $[-\gamma,\gamma]^{d}$. The (appropriately scaled) discrete Fourier transform of the corresponding grid values of $f$, given by

[TABLE]

yields an approximation of the Fourier coefficients $\widehat{f}_{k}$ for the same values of $k$. The approximation error can be quantified uniformly in $k\in\overline{\mathbb{Z}}_{N}^{d}$ 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 $f$ be $2\gamma-$periodic in each coordinate direction and f\in C^{0,\alpha}\big{(}[-\gamma,\gamma]^{d}\big{)} for some $\alpha>0$. Then

[TABLE]

and in particular, for any $k\in{\overline{\mathbb{Z}}}_{N}^{d}$,

[TABLE]

Proof.

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

[TABLE]

Moreover,

[TABLE]

The last term vanishes unless $m_{j}=0(\mathrm{mod}\ N)$ for each $j=1,\ldots,d$, in which case it takes the value $N^{d}$. Since $h^{d}N^{d}=(2\gamma)^{d}$, we have $Q(f)\ =\ \sum_{m\in N\mathbb{Z}^{d}}\widehat{f}_{m}$. Now (2.5) is obtained by noting that $\int_{[-\gamma,\gamma]^{d}}f(x)\,\mathrm{d}x=\widehat{f}_{{0}}$.

For fixed $k\in{\overline{\mathbb{Z}}}_{N}^{d}$, we introduce the function $g(x)=f(x)\exp(-\mathrm{i}\omega_{k}\cdot x)$, for $x\in[-\gamma,\gamma]^{d}$. Then $\widehat{f}_{k}=\int_{[-\gamma,\gamma]^{d}}g(x)\,\mathrm{d}x$ and $(S_{N}f)_{k}=Q(g)$. From (2.5) we conclude that

[TABLE]

Now (2.6) follows since $\widehat{g}_{m}=\widehat{f}_{k+m}$. ∎

2.2. Covariance functions

On a computational domain $D\subset\mathbb{R}^{d}$, we consider the fast evaluation of a Gaussian random field $Z(x,\omega)$ with covariance function given by (1.1). In this paper we consider the important case of Matérn covariance kernels with correlation length $\lambda>0$ and smoothness exponent $\nu>0$, given by

[TABLE]

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

[TABLE]

where

[TABLE]

The modified Bessel functions of the second kind $K_{\nu}$ have (see [1, 9.6.24]) the integral representations

[TABLE]

which shows in particular that their values for fixed $t$ are monotonically increasing with respect to $\nu$ for $\nu\geq 0$, and also that $K_{-\nu}=K_{\nu}$. We will also use the following results that directly imply exponential decay of $\rho_{\lambda,\nu}(x)$ as $\lvert x\rvert\to\infty$. Their proofs are given in Appendix A.

Lemma 2**.**

Let $\nu\geq 0$ and $t\geq 1/2$. Then we have

[TABLE]

Lemma 3**.**

Let $n\in\mathbb{N}_{0}$. Then \frac{\mathrm{d}^{n}}{\mathrm{d}t^{n}}\big{(}t^{\nu}K_{\nu}(t)\big{)} is of the form

[TABLE]

with coefficients $a_{n,j}$ satisfying

[TABLE]

Inspection of (2.8) shows that changing $\lambda$ amounts to a rescaling of the computational domain $D$. Hence without loss of generality, we may assume our computational domain $D$ to be contained in the box $[-\frac{1}{2},\frac{1}{2}]^{d}$, so that any difference of two points in $D$ is contained in $[-1,1]^{d}$. In what follows, we subsequently embed this box into a torus $[-\gamma,\gamma]^{d}$ with $\gamma\geq 1$, on which we define a $2\gamma$-periodic covariance $\rho^{\mathrm{ext}}$ such that $\rho^{\mathrm{ext}}(x)=\rho(x)$ for $x\in[-1,1]^{d}$, which means that the covariance between any pair of points in $D$ is preserved in replacing $\rho$ by $\rho^{\mathrm{ext}}$.

3. Classical Periodization

In this section, we treat the classical periodization (1.9). We prove in Theorem 6 that for Matérn covariances, the asymptotic decay of the eigenvalues of the extended matrix $\Sigma^{\mathrm{ext}}$ is the same as that of the underlying KL eigenvalues (1.13), up to a multiplicative factor which grows logarithmically in $\lvert\log(h)\rvert^{\nu}$. This confirms a recent conjecture [9, eq. (3.9)].

In this case $\Sigma^{\mathrm{ext}}$ is given by

[TABLE]

with $\rho^{\mathrm{ext}}$ defined in (1.8), (1.9). $\Sigma^{\mathrm{ext}}$ is a circulant extension of the covariance matrix $\Sigma$ of the form (1.2), obtained when sampling $Z$ at those points $\{x_{n}\}$ which lie in $[-1/2,1/2]^{d}$. Sampling on more general $d$-dimensional rectangles can be treated in the same manner with the obvious modifications. If the index set $n$ is given lexicographical ordering, then $\Sigma$ is a nested block Toeplitz matrix where the number of nested levels is the physical dimension $d$ and $\Sigma^{\mathrm{ext}}$ is a nested block circulant extension of it.

To analyse the eigenvalues of $\Sigma^{\mathrm{ext}}$ it is useful to also consider the continuous periodic covariance integral operator

[TABLE]

Then the scaled circulant matrix $h^{d}\Sigma^{\mathrm{ext}}$ can be identified as a Nyström approximation of ${\mathcal{R}}^{\rm ext}$, using the composite trapezoidal rule with respect to the uniform grid on $[-\gamma,\gamma]^{d}$ given by the points (2.3) with $n\in\overline{\mathbb{Z}}_{N}^{d}$. The operator ${\mathcal{R}}^{\rm ext}$ is a compact operator on the space of $2{\gamma}$-periodic continuous functions on $\mathbb{R}^{d}$, 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**.**

• (i)

The eigenvalues of ${\mathcal{R}}^{\rm ext}$ are $\widehat{\rho}_{k},\,k\in\mathbb{Z}^{d}$, as defined in (2.2), with corresponding eigenfunctions normalized in $L^{2}([-\gamma,\gamma]^{d})$:

[TABLE]

• (ii)

The eigenvalues of $h^{d}\Sigma^{\mathrm{ext}}$ are $(S_{N}\rho)_{k}$, $k\in{\overline{\mathbb{Z}}}_{N}^{d}$, as defined in (2.4), with corresponding eigenvectors normalized with respect to the Euclidean norm:

[TABLE]

Note that $S_{N}\rho^{\mathrm{ext}}=S_{N}\rho$ since $\rho$ coincides with $\rho^{\mathrm{ext}}$ on $[-\gamma,\gamma]^{d}$. It was convenient to introduce the scaling factor $h$ in (ii) above, since we can then identify (2.4) as the trapezoidal rule approximation of the ($2\gamma$-periodic) Fourier transform defining $\widehat{\rho}_{k}$.

Our analysis for estimating the decay of $(S_{N}\rho)_{k}$ as $|k|\rightarrow\infty$ will be based on the formula (2.6) and estimating the rate of decay of $\widehat{\rho}_{k}$. To this end, we now restrict our consideration to the Matérn covariances, that is, for the remainder of this section we assume

[TABLE]

with some $\lambda,\nu>0$, where $\rho_{\lambda,\nu}$ is defined in (2.8). In the recent paper [9] the following theorem was proven for this case.

Theorem 5**.**

Let $1/2\leq\nu<\infty$, $\lambda\leq 1$, and $h/\lambda\leq e^{-1}$. Then there exist $C_{1},C_{2}>0$ which may depend on $d$ but are independent of $\gamma,h,\lambda,\nu$, such that $\Sigma^{\mathrm{ext}}$ is positive definite if

[TABLE]

We are now in a position to formulate our result on the rate of decay of eigenvalues of $\Sigma^{\mathrm{ext}}$ for Matérn covariances, which proves the conjecture made in [9, eq. (3.9)] up to an additional factor of order $|\!\log h|^{\nu}$ in the constant.

Theorem 6**.**

Let $\nu$, $\lambda$, and $h$ be as in Theorem 5. Let $\gamma^{*}=\gamma^{*}(\lambda,\nu,h)$ denote the smallest value of $\gamma\geq 1$ such that condition (3.2) holds true and, adjusting $C_{2}$ if necessary, assume that $C_{2}>2\sqrt{2}(2d-2)$. Suppose $\gamma$ is chosen in the range $\gamma^{*}\leq\gamma\leq a\gamma^{*}$ for some $a\geq 1$ independent of $h$ and let $\Lambda^{\mathrm{ext}}_{j}$ denote eigenvalues of $\Sigma^{\mathrm{ext}}$ in non-increasing order. Then there exists $C>0$ such that for all $N$,

[TABLE]

The analysis which follows will be explicit in $j$, $h$, and $\gamma$ but not in the Matérn parameters $\lambda,\nu$ or in the dimension $d$. Note that in applications, one is typically interested in the $h\ll\lambda\leq\operatorname{diam}(D)$, and thus our assumptions on $\lambda$ and $h$ only exclude practically irrelevant cases. Recall that $C$ denotes a generic constant which may change from line to line and may depend on $\lambda,\nu$ or $d$.

The proof of Theorem 6 uses Proposition 4(ii), which tells us that $h^{d}\Lambda_{j}^{\mathrm{ext}}=(S_{N}\rho)_{k(j)}$ for some $k(j)\in\overline{\mathbb{Z}}_{N}^{d}$. Since Theorem 5 and Lemma 1 give us

[TABLE]

the proof proceeds by obtaining suitable estimates for the Fourier coefficients $\widehat{\rho}_{k}$ of the periodized covariance, as defined in (2.2). To this end, we use a cut-off function to isolate the artificial nonsmooth part of $\rho^{\mathrm{ext}}$ created by the classical periodization. We thus define an even smooth univariate cut-off function $\phi:\mathbb{R}\rightarrow\mathbb{R}$ by requiring that $\phi$ is supported on $[-3/4,3/4]$ and

[TABLE]

For any $\gamma>0$, we can scale this to a cut-off function supported on $[-3\gamma/4,3\gamma/4]$  by defining

[TABLE]

Using this function, we now write

[TABLE]

Thus $\beta$ coincides with $\rho$ in a neighbourhood of the origin and vanishes in a neighbourhood of the interface where $\rho^{\mathrm{ext}}$ has undergone its (nonsmooth) extension, while the support of $\sigma$ covers exactly this interface. In the following two lemmas, we separately estimate $\widehat{\beta}_{k}$ and $\widehat{\sigma}_{k}$, defined as in the right-hand side of (2.2).

Lemma 7**.**

For $r\in\mathbb{N}$, there exists $C>0$ independent of $\gamma\geq 1$ such that

[TABLE]

where $\widehat{\rho}$ is given by (2.9).

Proof.

We have $\beta=\rho\phi_{\gamma}$, since $\phi_{\gamma}$ vanishes outside $[-\gamma,\gamma]^{d}$. By the convolution theorem,

[TABLE]

The second term on the right can be estimated as

[TABLE]

with generic constants $C>0$ independent of $\gamma$, where we have used

[TABLE]

and the decay properties of $\widehat{\rho}$ given in (2.9).

For the first term on the right of (3.6), we use two separate bounds, suitable for for small and large $\lvert\omega_{k}\rvert$, to obtain

[TABLE]

Now note that $\lVert\widehat{\rho}\rVert_{L^{\infty}(\mathbb{R}^{d})}\lVert\widehat{\phi}_{\gamma}\rVert_{L^{1}(\mathbb{R}^{d})}$ is uniformly bounded with respect to $\gamma$ by (3.8) and (2.9). Since $\phi_{\gamma}$ is smooth, integration by parts yields

[TABLE]

Thus, since $\widehat{\rho}\in L^{1}(\mathbb{R}^{d})$, we have, for any $r\in\mathbb{N}$,

[TABLE]

Combining (3.7) and (3.9) with (3.6) completes the proof. ∎

For estimating $\lvert\widehat{\sigma}_{k}\rvert$, we use the following auxiliary result, which is proved in the appendix.

Lemma 8**.**

Let $\alpha\in\mathbb{Z}^{d}$ with $|\alpha|_{\infty}\leq 2$ and $|x|\geq\gamma/2$ with $\gamma\geq 1$. Then, with $\sigma$ as given in (3.5), there exists $C$ independent of $\gamma$ such that

[TABLE]

Lemma 9**.**

Let $1/2\leq\nu<\infty$. Then there exists $C>0$ independent of $\gamma\geq 1$ such that

[TABLE]

Proof.

First we bound $\sigma_{k}$ in terms of $\sigma$ and its derivatives. When $\max_{i}\lvert\omega_{k,i}\rvert\leq 1$ we use the estimate

[TABLE]

For the case where $\lvert\omega_{k,i}\rvert>1$ for at least one value of $i$, we integrate by parts dimensionwise to best exploit the limited smoothness of $\sigma$ across the boundary of $[-\gamma,\gamma]^{d}$. We assume without loss of generality that $\lvert\omega_{k,1}\rvert>1$ and we just give the proof for $d=2$ to simplify the exposition; higher dimensions are analogous. Integration by parts in (3.11) twice with respect to $x_{1}$ gives

[TABLE]

Now denoting

[TABLE]

we get $\tilde{\sigma}_{1}(\gamma)=\tilde{\sigma}_{1}(-\gamma)$ by the periodicity of $\sigma(x)$ as a function of $x_{2}$. If, in addition, $\lvert\omega_{k,2}\rvert\geq 1$, then integrating by parts with respect to $x_{2}$ gives

[TABLE]

To estimate the right-hand sides of (3.11), (3.12) and (3.13), we use Lemma 8, which directly gives the desired bound for the last term on the right-hand side of (3.13). Moreover,

[TABLE]

which is the required bound for (3.11) and the first terms on the right-hand sides of (3.12) and (3.13). Finally, for the second term on the right hand side of (3.13) we have with $\alpha=(1,2)$

[TABLE]

Carrying this out in the same way for further other terms we obtain the desired result for $d=2$. The proof for $d>2$ can be done analogously, giving rise to $2^{d}$ terms in (3.13). ∎

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). For the first term, combining Lemmas 7 and 9 yields

[TABLE]

with $C>0$ independent of $\gamma$ and $r\in\mathbb{N}$ arbitrary. We choose fixed $r\geq\nu+d/2$. For the first term on the right of (3.14), using (2.9) with $\omega_{k}=\pi k/\gamma$ gives

[TABLE]

Moreover, the second term of (3.14) can be estimated by

[TABLE]

To see the first inequality in (3.15), consider $k={0}$ and $k\not={0}$ separately. In the latter case the inequality follows from the elementary estimate $\pi|k|\geq(1+|k|$). Estimating the remaining term of (3.14) in a similar way, we obtain

[TABLE]

The second term on the right-hand side of (3.16) will turn out to be dominated by the first. But to finish the argument we also have to estimate the second term on the right-hand side of (3.4). To do this note that since $\max_{i=1,\ldots,d}\lvert k_{i}\rvert\leq N/2$ for $k\in{\overline{\mathbb{Z}}}_{N}^{d}$, and $m\in\mathbb{Z}^{d}$, we have

[TABLE]

for $m\in\mathbb{Z}^{d}$ with $C$ independent of $k$ and $N$. Thus, we get from (3.16)

[TABLE]

Now, by elementary arguments,

[TABLE]

Inserting this into the second term of (3.17), and estimating the first term similarly, we obtain

[TABLE]

We now show that the first term in (3.16) is dominant in both estimates (3.16), (3.18). First note that, by elementary arguments, $\prod_{i=1}^{d}(1+\lvert k_{i}\rvert)^{2}\ \geq\ 1+|k|^{2}\ \geq\ \frac{1}{2}(1+|k|)^{2}$. Then, for $k\in\overline{\mathbb{Z}}_{N}^{d}$, we have

[TABLE]

Now by choice of $\gamma$ and (3.2), we have $\gamma\,\geq\,\gamma^{*}\,>\,C_{2}\lambda\sqrt{\nu}\log(\lambda/h)$, and using the definition of $L$ in (3.10), we obtain

[TABLE]

Then, combining (3.19) and (3.20) and recalling that $h=2\gamma/N$, we obtain

[TABLE]

By choice of $C_{2}$, and using $\nu\geq 1/2$, we have $C_{2}>2\sqrt{2}(2+(d-2)/\nu)$ and so the exponent of $h$ in (3.21) is positive. Also since $\gamma\leq a\gamma^{*}$, $\gamma$ grows at most logarithmically in $h$ with a multiplicative constant which grows at most linearly in $a$. This yields a bound on the second term in (3.16) and thus

[TABLE]

Turning to the second term on the right-hand side of (3.18) we obtain, similarly,

[TABLE]

when $k\in\overline{\mathbb{Z}}_{N}^{d}$. Inserting (3.21) and (3.22) into (3.16) and (3.18), we obtain

[TABLE]

Now, to finish the proof, let $\{\lvert k^{*}_{j}\rvert:\ j=1,\ldots,N^{d}\}$, be a non-decreasing ordering of the numbers $\{\lvert k\rvert:\,k\in\overline{\mathbb{Z}}_{N}^{d}\}$. As shown in [9, Theorem 3.4], $\lvert k^{*}_{j}\rvert$ is then uniformly proportional to $j^{1/d}$, with constants that depend only on $d$. Thus by (3.23),

[TABLE]

Now by the hypothesis of the theorem, the numbers $\lambda^{*}_{j}:=N^{-d}\Lambda^{\mathrm{ext}}_{j}$, $j=1,\ldots,N^{d}$, are non-increasing and, by Proposition 4(ii), provide a non-increasing ordering of the values $(2\gamma)^{-d}(S_{N}\rho)_{k}$, $k\in{\overline{\mathbb{Z}}}_{N}^{d}$. Then we claim that, for any integer $0<J<N^{d}$,

[TABLE]

with $C$ as in (3.24). If this were not true then, by (3.24), and for some $J$,

[TABLE]

Since the terms in the right-hand sum also provide an ordering for the eigenvalues $\lambda_{j}^{*}$ this contradicts the assumed non-increasing property of $\lambda^{*}_{j}$. As a consequence of (3.25), we then have

[TABLE]

With the condition (3.2) on $\gamma$, implying that $\gamma\leq Ca\log(\lambda/h)$, this completes the proof. ∎

4. Smooth Periodization

We now establish a sufficient criterion on the periodization cell size $\gamma$ ensuring a positive definite periodic covariance function in the case of periodization with a smooth cutoff function as in (1.10). This amounts to proving a quantitative version of [3, Theorem 2.3] for the case of Gaussian random fields with Matérn covariance as in (2.8).

We explicitly construct a suitable even cutoff function $\varphi_{\kappa}$ which vanishes outside $B_{\kappa}({0})$ such that $\varphi_{\kappa}=1$ on $[-1,1]^{d}$ and $\rho_{\lambda,\nu}\varphi_{\kappa}$ is a positive definite function. In this case, for $\gamma$ sufficiently large, there exists a periodic Gaussian random field on the torus $[-\gamma,\gamma]^{d}$ with the periodized covariance kernel

[TABLE]

such that $\rho^{\mathrm{ext}}=\rho_{\lambda,\nu}$ on $[-1,1]^{d}$, so that the corresponding random fields have the same law on the domain of interest contained in $[-\frac{1}{2},\frac{1}{2}]^{d}$. As shown in [3, §3], if $\varphi_{\kappa}$ is sufficiently smooth, the eigenvalues of the covariance operator of the periodized random field then have the same asymptotic decay as those of the corresponding Matérn covariance operator.

The cutoff function defined here is different from $\phi$ used in Section 3, which served only as a tool in the proof of Theorem 6. By contrast, the cutoff function $\varphi_{\kappa}$ on $\mathbb{R}^{d}$ derived from a univariate cutoff function $\varphi$ constructed here is used numerically in the computation of the random field. In order to cover the full range of Matérn smoothness parameters $\nu>0$ in our analytical results, we need precise control of high-order derivatives of $\varphi$. Specifically, for $p:=\lceil\nu+\frac{d}{2}\rceil$ we require bounds of the form

[TABLE]

with some $c_{1},c_{2}>0$. The use of such $\varphi$ mainly allows us to circumvent some further major technicalities in our proofs, and as the numerical tests in Section 5 show, one still observes similar results for cutoff functions for which no bound of the form (4.2) is available.

Our concrete choice of $\varphi$ is as follows: let $N_{P}$ be the B-spline function with nodes $\{-P,\ldots,-1,0\}$, where $P:=2p+1$. For $\kappa>0$ we define the even function $\varphi\in C^{2p}(\mathbb{R})$ by

[TABLE]

It is easy to see that $\varphi(t)=0$ if $|t|\geq\kappa$. This choice of $\varphi$ provides us with explicit bounds of the form (4.2) on all required derivatives. From $N^{\prime}_{r+1}(t)=N_{r}(t)-N_{r}(t-1)$ we infer that for $0\leq\alpha\leq 2p$,

[TABLE]

We now define

[TABLE]

With this choice of $\varphi_{\kappa}$ in (4.1), we have $\rho^{\mathrm{ext}}=\rho_{\lambda,\nu}$ on $[-1,1]^{d}$ provided that

[TABLE]

which reduces to the condition in (1.7) when $d=1$.

In our following main result, we pursue the basic strategy of [3, Theorem 2.3] to establish sufficient conditions in terms of $\nu,\lambda$ on the required value of $\kappa>0$ such that

[TABLE]

In [3], for a more general class of covariance functions than considered here, only the existence of such $\kappa$ is established without further information on its size. The proof uses the integrability of derivatives of $\rho_{\lambda,\nu}\theta_{\kappa}$ to show that $\widehat{\rho_{\lambda,\nu}\theta_{\kappa}}$ decays at least as fast as $\widehat{\rho_{\lambda,\nu}}$, and then uses the exponential spatial decay of these derivatives to show that $\widehat{\rho_{\lambda,\nu}\theta_{\kappa}}$ can indeed be bounded by $\widehat{\rho_{\lambda,\nu}}$ if $\kappa$ is chosen sufficiently large. Here, we follow the same basic strategy, but extract information on the required size of $\kappa$. This needs detailed information on higher-order derivatives of $\rho_{\lambda,\nu}$ and $\theta_{\kappa}$, where the order increases with the value of $\nu$. The proof of Theorem 5 for the classical periodization relies in a similar manner on using spatial decay of $\rho_{\lambda,\nu}$ to control a perturbation term, but does not require derivative information.

Theorem 10**.**

For $d\in\{1,2,3\}$ and $\varphi_{\kappa}$ as defined above, there exist constants $C_{1},C_{2}$ such that for any $0<\lambda,\nu<\infty$, we have $\widehat{\rho_{\lambda,\nu}\varphi_{\kappa}}>0$ provided that $\kappa>1$ and

[TABLE]

Note that the condition (4.8), together with (4.6), is similar to the one in Theorem 5. There are two key differences: the restriction $\nu\geq\frac{1}{2}$ does not appear, and since there is no discretization involved in Theorem 10, there is no dependence on a grid size $h$ as in Theorem 5. The restriction to the dimensionalities $d\in\{1,2,3\}$ that are relevant in applications is not essential, but allows us to avoid some further technicalities in the proof.

Proof.

Note that since $\widehat{\rho_{\lambda,\nu}\varphi_{\kappa}}(\omega)=\lambda^{d}\widehat{\rho_{1,\nu}\varphi_{\frac{\kappa}{\lambda}}}(\lambda\omega)$ for any $\lambda>0$, if (4.7) holds with some $\kappa_{1}$ for $\rho_{1,\nu}$, then it also holds with $\kappa:=\lambda\kappa_{1}$ for $\rho_{\lambda,\nu}$. Consequently, it suffices to consider the case $\lambda=1$ in what follows, and we write $\rho=\rho_{1,\nu}$ and $\theta=\theta_{\kappa}$. We consider the cases $\frac{1}{2}\leq\nu<\infty$ and $0<\nu<\frac{1}{2}$ separately.

Step 1. The case $\nu\geq\frac{1}{2}$. With $r=|\omega|$, for the Fourier transform of the radial function $\rho\theta$ we have

[TABLE]

where $J_{\alpha}$ is the classical Bessel function of order $\alpha$. Now condition (4.7) with $\lambda=1$ is equivalent to

[TABLE]

where $C_{1,\nu}$ is given in (2.10). Since $\nu\geq\frac{1}{2}$ and $p\geq\nu+\frac{d}{2}$, it thus suffices to choose $\kappa$ such that

[TABLE]

In what follows, as a consequence of (4.8) we can assume without loss of generality that

[TABLE]

with $d_{1}=1+2(d-1)$. We now proceed to estimate

[TABLE]

with

[TABLE]

We consider the case $\ell\geq 1$. Using \frac{\mathrm{d}}{\mathrm{d}z}\big{[}z^{\alpha}J_{\alpha}(z)\big{]}=z^{\alpha}J_{\alpha-1}(z), see [18, page 45], we can write

[TABLE]

Integrating by parts, we get from $[\rho\theta](\frac{\kappa}{2})=0$ and the exponential decay of the Matérn covariance function

[TABLE]

where in the last step we use \frac{\mathrm{d}}{\mathrm{d}z}\big{[}z^{-\alpha}J_{\alpha}(z)\big{]}=-z^{-\alpha}J_{\alpha+1}(z) (see [18, page 45]). Integrating by parts again we arrive at

[TABLE]

Repeating this argument we conclude that

[TABLE]

where derivatives are taken $2\ell$ times. Employing Lommel’s expression of $J_{\alpha}$, see [18, page 47],

[TABLE]

and $J_{-1/2}(z)=\sqrt{\frac{2}{\pi}}\frac{\cos z}{\sqrt{z}}$, we can bound

[TABLE]

where in the last equality we have used the relation between Gamma and Beta functions, see [1, Section 6.2]. Consequently

[TABLE]

with $C_{0}=\frac{(1/2)^{d/2-1}}{\Gamma(d/2)}$ for $d\in\{1,2,3\}$.

To finish the proof we need a technical lemma. For $\ell\in\mathbb{N}$ and $f$ having derivatives of sufficiently high order, we denote

[TABLE]

Lemma 11**.**

For $\ell\geq 1$ and $f$ having $2\ell$-th derivative, the term $B_{2\ell,d}(f,t)$ has the form

[TABLE]

with $B_{2\ell,1}(f,t)=f^{(2\ell)}(t)$, $B_{2\ell,3}(f,t)=f^{(2\ell)}(t)t^{2}+2\ell f^{(2\ell-1)}(t)t$, and when $d=2$ we have

[TABLE]

The proof of this lemma is given in Appendix A. We continue the proof of Theorem 10 by using the above lemma to obtain the estimate

[TABLE]

Thus, from (4.4) we get

[TABLE]

which with (4.14) leads to

[TABLE]

From Lemma 3, we obtain

[TABLE]

with $a_{n,j}$ as in (2.13), and consequently

[TABLE]

Since $\max\{2\nu-n+j,n-j\}\leq 2p$ we get $|\nu-n+j|\leq 2p-\nu$ which implies $K_{\nu-n+j}(t)\leq K_{2p-\nu}(t)$ by the representation (2.11). Again with the assumption $\tau:=\frac{\kappa\sqrt{2\nu}}{2}=\frac{\kappa\sqrt{2\nu}}{2\lambda}\geq 2P$ and $j+2\ell-\alpha\geq 0$ we can estimate

[TABLE]

where in the second step we have used Lemma 3. Inserting this into (4.16) we obtain

[TABLE]

These estimates hold form all $\ell\geq 1$. Moreover, using (4.12) and $|\theta(t)|\leq 1$ for all $t\geq 0$ we have from (4.11)

[TABLE]

Consequently

[TABLE]

Since $\nu+d-1\leq\nu(1+2(d-1))$, as a consequence of (2.12) with $d_{1}=1+2(d-1)$,

[TABLE]

Using the assumption $\tau=\frac{\kappa\sqrt{2\nu}}{2}\geq 2\nu d_{1}$ as well as

[TABLE]

we obtain

[TABLE]

Combining (4.19) and (4.20), we arrive at

[TABLE]

Now the required bound (4.10) follows from

[TABLE]

Since $(2\nu)^{p-\nu-\frac{d}{2}}\leq 2\nu$, a sufficient condition for (4.21) is

[TABLE]

with $C>0$ independent of $\kappa$ and $\nu$.

Taking logarithms and using the Stirling bounds

[TABLE]

as well as $p\leq\nu+\frac{d}{2}+1$, for general $\lambda>0$ shows that the condition

[TABLE]

where $C_{1},C_{2}$ depend only on $d$ (or more precisely, $C_{1},C_{2}=\mathcal{O}(d\log d)$) is sufficient to ensure (4.22)

Step 2. The case $\nu<\frac{1}{2}$. Since we restrict ourselves to $d\leq 3$, we have $p=\lceil\nu+\frac{d}{2}\rceil\leq 2$. Repeatedly using the identities $K_{\nu}^{\prime}(t)=-K_{\nu-1}(t)-\frac{\nu}{t}K_{\nu}(t)$ and $K_{\nu}^{\prime}=\frac{\nu}{t}K_{\nu}-K_{\nu+1}$, we obtain

[TABLE]

Next, we note that $K_{\nu-1}=K_{1-\nu}$, and by [1, 10.2.17]

[TABLE]

where the monotonicity in $\nu$ can be seen from the explicit representation (2.11). Continuing from (4.16), we now estimate for $0\leq\ell\leq p$

[TABLE]

By (4.8) we may use the assumption $\kappa\sqrt{2\nu}/2\geq 1$, and thus

[TABLE]

if $t\geq\kappa/2$. As a consequence,

[TABLE]

and

[TABLE]

For $0<\nu<\frac{1}{2}$, the condition for (4.9) to hold becomes

[TABLE]

or equivalently $e^{\frac{\kappa\sqrt{2\nu}}{4}}\geq C\left(1+\lvert\omega\rvert^{2}/(2\nu)\right)^{\nu+\frac{d}{2}-p}$ with $C>0$ independent of $\kappa$, $\nu$, which is implied by $\kappa\geq(2\sqrt{2}\ln C)\,\nu^{-1/2}$. This completes the proof. ∎

5. Numerical Experiments

The eigenvalue decay established in Theorem 6 has already been studied numerically in [9]. Note that the results given there are also consistent with the presence of the extra logarithmic factor in (3.3). Here, we thus focus on a numerical study of the required extension size $\gamma$. In order to assess the sharpness of the necessary conditions (3.2) and (4.8), we use a simple bisection scheme to find the minimum value of $\gamma$ that is actually required in each case to ensure that the obtained covariance matrix is positive definite. In all tests, we assume the box $[-1,1]^{d}$ as the computational domain, and we show results only for $\lambda=\frac{1}{2}$ since the resulting values of $\gamma$ exhibit an approximately linear scaling with respect to $\lambda$.

In the case of smooth periodization, in addition to the cutoff function using integrated B-splines defined in (4.3), (4.5) we also test a standard infinitely differentiable cutoff function as used in [3], which is simpler to implement in practice: let

[TABLE]

One can then replace the definition of $\varphi$ in (4.3) by

[TABLE]

and again define $\varphi_{\kappa}(x):=\varphi(\lvert x\rvert)$.

First, we compare the extension sizes $\gamma$ in terms of the grid size $h$ that are needed for classical circulant embedding and for smooth periodization; Figure 2 shows the resulting ratio of the number of grid points in the extension torus $\mathbb{T}=[-\gamma,\gamma]^{d}$ to the number of sampling grid points in the original domain. Whereas, as expected in view of Theorem 10, the minimum required values of $\gamma$ are indeed independent of $h$ in the case of the smooth truncation, in the case of the classical circulant embedding the corresponding values of $\gamma$ indeed exhibit a dependence of order $\lvert\log h\rvert^{d}$ on $h$. In this sense, we observe the result of Theorem 5 to be sharp. Especially for $d=3$ and smaller values of $h$, the smooth periodization leads to substantially more favorable extension sizes. The results shown are for the $C^{\infty}$-cutoff function (5.1), and one obtains very similar results with (4.3).

In addition, in Figure 3 we consider the dependence of the minimum required $\gamma$ on $\nu$ for the periodization with smooth truncation and compare to the asymptotics in the condition (4.8) of Theorem 10. We show results for both cutoff function constructions (4.3) and (5.1). We observe that in the particular case $d=1$, $\gamma$ remains bounded as $\nu\to 0$ (indeed, we observe $\gamma\to 1$ in this limit), whereas for $d>1$ we find an increase in $\gamma$ both as $\nu\to 0$ and $\nu\to\infty$. The actual required increase of $\gamma$ as $\nu\to 0$ appears to be slightly slower than the order $\nu^{-1/2}$ in (4.8). The observed behaviour of $\gamma$ for larger $\nu$ is consistent with the sufficient condition of order $\nu^{1/2}\log\nu$ in (4.8). Note that the B-spline cutoff of limited smoothness leads to a slower increase of $\gamma$ as $\nu\to\infty$, whereas $\gamma$ in this case increases slightly faster as $\nu\to 0$ for $d>1$.

6. Conclusions

We have seen that both classical and smooth periodization preserve the asymptotic decay rate of covariance eigenvalues. Concerning the factor $\gamma\geq 1$ by which the sampling grid needs to be extended to ensure positive definiteness, there is a difference: the smooth periodization requires an extension factor $\gamma$ that is independent of the sample grid size $h$, whereas classical periodization requires $\gamma\sim\lvert\log h\rvert$ according to Theorem 5. For the extended grid size $N^{d}=(2\gamma h^{-1})^{d}$ we thus have $N^{d}\sim h^{-d}$ in the case of smooth periodization, and $N^{d}\sim h^{-d}\lvert\log h\rvert^{d}$ in the case of classical periodization. This is exactly what we observe for the numerically determined minimum extension sizes in Figure 3. This directly translates to efficiency advantages of the smooth periodization that are more pronounced for larger $d$. The total cost is dominated by the FFT on the extended grid, which requires $\mathcal{O}(N^{d}\log N)$ operations. Whereas in the case of the classical circulant embedding, one sample on $D$ thus costs $\mathcal{O}(h^{-d}\lvert\log h\rvert^{d+1})$ operations, with the smooth periodization the cost is only $\mathcal{O}(h^{-d}\lvert\log h\rvert)$ in terms of $h$.

One further useful consequence of the new result (1.11) covering the full range of $\lambda$, $\nu$ 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 $\lambda$ and $\nu$ 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 factorization for a new covariance and drawing a sample come at the same cost of one FFT each: thus by (1.11), provided that the distributions of $\lambda$, $\nu$ have compact supports inside $(0,\infty)$, smooth periodization is guaranteed to provide each sample – without any further approximation – at near-optimal cost $\mathcal{O}(h^{-d}\lvert\log h\rvert)$.

The periodic random field $Z^{\mathrm{p}}$ on $\mathbb{T}$, obtained using smooth truncation, also provides a tool for deriving series expansions of the original random field. As mentioned above, in the smooth periodization case the KL eigenvalues of the periodic field have the decay rate (1.13). Moreover, in contrast to the KL eigenfunctions on $D$, which are typically not explicitly known, the corresponding eigenfunctions $\varphi^{\mathrm{p}}_{j}$ of the periodic covariance are explicitly known trigonometric functions and one has the following KL expansion for the periodized random field:

[TABLE]

with $\lambda^{\mathrm{p}}_{j}$ denoting the eigenvalues of the periodized covariance and the $\varphi^{\mathrm{p}}_{j}$ are normalized in $L^{2}(\mathbb{T})$. Restricting this expansion back to $D$, one obtains an exact expansion of the original random field on $D$ in terms of independent scalar random variables

[TABLE]

This provides an alternative to the standard KL expansion (1.6) of $Z$ in terms of eigenvalues $\lambda_{j}$ and eigenfunctions $\varphi_{j}$ normalized in $L^{2}(D)$. The main difference is that the functions \varphi^{\mathrm{p}}_{j}\big{|}_{D} in (6.1) are not $L^{2}(D)$-orthogonal. However, these functions are given explicitly, and thus no approximate computation of eigenfunctions is required.

The eigenvalues $\lambda^{\mathrm{p}}_{j}$ can be approximated efficiently by FFT as described in [3, Sec. 5.1] and the asymptotic decay of $\lambda^{\mathrm{p}}_{j}$ is the same as that of $\lambda_{j}$ in (1.6). Moreover, since the eigenfunctions of the periodized covariance are explicitly known trigonometric functions, the $L^{\infty}$-norms of the scaled eigenfunctions in the expansion (6.1) also decay at the same rate as $\sqrt{\lambda_{j}}$, that is,

[TABLE]

This decay of the $L^{\infty}$-norms is important for applications, e.g., to random PDEs. Since $\lVert\varphi_{j}\rVert_{L^{\infty}(D)}$ is in general not uniformly bounded, see [3, Sec. 3], this constitutes a marked advantage of the expansion (6.1) over the standard KL expansion (1.6). In addition, on geometrically complicated $D$, (6.1) is significantly easier to handle numerically.

The KL expansion of $Z^{p}$ also enables the construction of alternative expansions of $Z$ of the basic form (6.1), but with the spatial functions having additional properties. In [3], wavelet-type representations

[TABLE]

are constructed by applying the square root of the covariance operator in the corresponding factorisation to periodic Meyer wavelets, with the summation running over $\ell\geq 0$ and $k\in\{0,\ldots,2^{\ell}-1\}^{d}$. The functions $\psi_{\ell,k}$ 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 [2].

Appendix A Proofs of Auxiliary Results

Proof of Lemma 2..

If $0\leq\nu\leq\frac{1}{2}$ we use the estimate $K_{\nu}(t)\leq K_{1/2}(t)=\sqrt{\frac{\pi}{2t}}e^{-t}$. When $\nu>1/2$ we have (see [18, Page 206])

[TABLE]

Changing variable we obtain

[TABLE]

Now using the duplication formula $\Gamma(\nu)\Gamma(\nu+1/2)=2^{1-2\nu}\sqrt{\pi}\Gamma(2\nu)$ we get the desired result. ∎

Proof of Lemma 3..

We proceed by induction. The statement in Lemma 3 holds for $n=0$ and $n=1$ since \frac{\mathrm{d}}{\mathrm{d}t}\big{(}t^{\nu}K_{\nu}(t)\big{)}=-t^{\nu}K_{\nu-1}(t) which follows from $tK_{\nu}^{\prime}(t)+\nu K_{\nu}(t)=-tK_{\nu-1}(t)$, see [18, Section 3.71]. Assume that the statement holds for some $n\geq 1$. We write

[TABLE]

Taking derivatives on both sides and using \frac{\mathrm{d}}{\mathrm{d}t}\big{(}t^{\tau}K_{\tau}(t)\big{)}=-t^{\tau}K_{\tau-1}(t) we get

[TABLE]

This is of the form given in (2.13), moreover

[TABLE]

where we have used the induction hypothesis. The proof is completed. ∎

Proof of Lemma 8.

For $|\alpha|_{\infty}\leq 2$ and $|x|\geq\gamma/2$, we have

[TABLE]

which implies, since $\gamma\geq 1$,

[TABLE]

for some positive constant $C$ depending on $\alpha$. Since $r=|x|\geq\gamma/2$ implies in particular $r\geq 1/2$, we have $|\partial^{\alpha}r|\leq C$ (which in general depends only on $d$) and therefore obtain from Lemma 3

[TABLE]

Note that, using $K_{\nu}=K_{-\nu}$, by the monotonicity of $K_{\nu}$ in $\nu$ for $\nu\geq 0$ and the series representation of $K_{\nu}$ [18, Page 80], we can estimate, for all $t\geq 0$,

[TABLE]

Now, if $|x|\geq\gamma/2$, then $t:=\sqrt{2\nu}|x|/\lambda\geq C$, (A.1) implies that $K_{\nu}(t)\leq Ct^{-1/2}e^{-t}$ and so

[TABLE]

Proof of Lemma 11.

The assertion is obvious for $d=1$. With $\ell=1$ we have $B_{2,d}=f^{\prime\prime}(t)t^{d-1}+(d-1)f^{\prime}(t)t^{d-2}\,.$ Assuming that (4.13) is true for $\ell$, we prove the this statement for $\ell+1$. We have

[TABLE]

and

[TABLE]

When $d=2$ we obtain $\sum_{\alpha=1}^{2}\lvert a_{2,\alpha}\rvert=d=2$ and

[TABLE]

This proves (4.13). When $d=3$, by the induction assumption we have

[TABLE]

Using (A.2) again we obtain

[TABLE]

This completes the proof of the lemma. ∎