On the consistency of inversion-free parameter
estimation for Gaussian random fields ✩
Hossein Keshavarza , Clayton Scottb,a , XuanLong Nguyena,b
a

arXiv:1601.03822v2 [math.ST] 21 Jun 2016

b

Department of Statistics, University of Michigan
Department of Electrical Engineering and Computer Science, University of Michigan

Abstract
Gaussian random fields are a powerful tool for modeling environmental processes. For high dimensional samples, classical approaches for estimating the covariance parameters require highly
challenging and massive computations, such as the evaluation of the Cholesky factorization or solving linear systems. Recently, Anitescu, Chen and Stein [2] proposed a fast and scalable algorithm
which does not need such burdensome computations. The main focus of this article is to study the
asymptotic behavior of the algorithm of Anitescu et al. (ACS) for regular and irregular grids in
the increasing domain setting. Consistency, minimax optimality and asymptotic normality of this
algorithm are proved under mild differentiability conditions on the covariance function. Despite
the fact that ACS’s method entails a non-concave maximization, our results hold for any stationary
point of the objective function. A numerical study is presented to evaluate the efficiency of this
algorithm for large data sets.
Keywords: Inversion-free estimation, Covariance function, Stationary Gaussian process,
Asymptotic analysis

1. Introduction
Gaussian processes have plethora of applications, ranging from the modeling of environmental
processes in geostatistics (e.g., [9, 11]) to supervised regression and classification in machine learning
[5, 8, 16], and the simulation of complex computer models [18]. The versatility of the correlation
structure of Gaussian processes provides a tractable and powerful tool for the modeling of large and
highly dependent environmental variables. As a common approach in the field of spatial statistics,
the covariance functions of Gaussian processes are assumed to belong to a parametric family. High
precision estimates of the covariance parameters are pivotal for interpolating Gaussian processes
which is the ultimate goal in many geostatistical problems [9, 19].
In the last two decades, there has been extensive research regarding the statistical and computational facets of the estimation of the Gaussian processes’ covariance parameters. Maximum
likelihood estimation (MLE) was the earliest favored algorithm in the geostatistics community, e.g.,
✩
This research is partially supported by NSF grant ACI-1047871. Additionally, CS is partially supported by NSF
grants 1422157, 1217880, and 0953135, and LN by NSF CAREER award DMS-1351362, NSF CNS-1409303, and NSF
CCF-1115769.
Email: {hksh,clayscot,xuanlong}@umich.edu

Preprint submitted to Journal of Multivariate Analysis

June 22, 2016

Mardia et al. [13] and Ying [25]. However, solving systems of linear equations is inevitable to evaluate the Gaussian likelihood. Notwithstanding the recent advances toward scalable and efficient
solution of the system of linear equation ( e.g., iterative Krylov subspace method or block preconditioned conjugate gradient algorithm [15]) which moderately reduces the computational and memory
costs of the direct evaluation of the precision matrix, obtaining the MLE of unknown covariance
parameters using such linear systems solvers is still a strenuous task, especially for a generic Gaussian spatial process observed at numerous and possibly irregularly spaced locations. Approximating
the likelihood function by tapering the covariance matrix is another class of algorithms aiming to
reduce the numerical burden of MLE (see Kaufman et al. [12]). Due to the sparsity of the tapered
covariance matrix, its inverse can be computed in a faster and more stable way. Recent studies
[10, 12, 22] demonstrate the consistency and asymptotic normality of this algorithm under some
mild conditions on the taper function.
Because of the obstacles of solving system of linear equations for massive data, which is necessary
for tapered and exact MLE, it is of great interest to develop estimation techniques without requiring
such extensive computations. Such class of algorithms, which will be referred to as inversion-free,
are based upon optimizing a loss function whose form (and its derivatives) is independent of the
precision matrix of data. The first attempt toward such a goal has been done by Anitescu, Chen
and Stein [2] in 2014 (referred here to as ACS). Their proposed procedure is faster and more stable
than likelihood based algorithms. In [2], the covariance parameters are estimated by computing
the global maximizer of a non-concave program. Simulation studies verify the efficiency of ACS’s
approach in the case that the covariance matrix has a bounded condition number. The main
purpose of this paper is to appraise the asymptotic properties of the ACS’s algorithm such as
consistency, minimax optimality and asymptotic normality. The developed theory in this paper
shows that ACS’s algorithm has the same asymptotic rate of convergence as the MLE. In practice,
the solution of ACS’s optimization problem may also serve as the starting point (initial guess) of a
likelihood maximization procedure.
In geostatistics, there are two common asymptotic regimes: increasing-domain and fixeddomain, the latter sometimes referred to as infill asymptotics (see [19], Section 3.3 or [9], p. 480).
In the former setting, the minimum distance among the sampling points is bounded away from
zero and more samples are collected by increasing the diameter of the spatial domain. In the latter
regime, the data are sampled in a fixed and bounded domain, and the observations get denser as
the sample size n increases.
Zhang [26] showed that not all the covariance parameters are consistently estimable in the fixed
domain regime. Strictly speaking, there is no asymptotically consistent algorithm for estimating
the non-micro ergodic covariance parameters, which do not asymptotically affect the interpolation
mean square error (see [19] for a precise definition). On the other hand, it is known in the literature
that subject to some mild regularity conditions, maximizing the likelihood provides a strongly
consistent and asymptotically normal estimate for all the covariance parameters in the increasing
domain setting [4, 13].
Increasing domain asymptotic analysis of covariance estimation has two significant benefits.
First, unlike the infill asymptotic setting, the geometry of the spatial sampling has a crucial impact
on the asymptotic distribution of the parameter estimate. Thus, this regime is a natural asymptotic
framework for assessing the role of irregularity of spatial sampling on the covariance parameter
estimation [4]. This claim can be verified by a deeper look at the asymptotic distribution of
the microergodic parameter estimates in the fixed domain (see e.g., [19, 25] for MLE and [10,

2

12, 22] for tapered MLE). Another significant characteristic of increasing domain regime is that
the covariance matrix has a universally bounded condition number as n grows under some mild
regularity conditions. This feature of the covariance matrix plays a major role in our asymptotic
analysis. Although in many geostatistical applications in a fixed bounded domain the condition
number of the covariance matrix increases at least linearly with respect to n, preconditioning filters
is commonly used to uniformly control the condition number independent of n [7, 19]. Therefore,
we believe that our developed increasing domain asymptotics can be useful for the fixed domain
analysis of preconditioned inversion-free algorithms.
Outline of main results.. This paper studies the increasing domain asymptotic behaviour of ACS’s
estimation algorithm introduced in [2]. Specifically, suppose that G is a zero mean stationary
Gaussian process in Rd with covariance function cov (G (s) , G (s′ )) = R (s − s′ , η) in which η ∈ Ω
denotes the vector of unknown covariance parameters. One realization of G has been observed on a
d−dimensional perturbed regular lattice of n = N d points, which will be formally defined in Section
2. The specific contributions of this work are given as follow:
(a) Assuming the polynomial decay of R (s, η) and its gradient (with respect to η) in terms of
the Euclidean norm of s, and under some mild identifiability condition on R, we prove that
the global maximizer
√ of ACS’s method consistently estimates η. Furthermore, the
√ estimation
error is of order n−1 ln n which is shown to be minimax optimal up to some ln n term.
(b) As the proposed loss function in Anitescu et al. [2] is not jointly concave in η, finding its
global maximizer is challenging. For a large enough sample size and under an additional
condition regarding the polynomial decay of the second derivative of R (s, η) with respect to
η, we show that any stationary
√ point of this non-concave program is concentrated around the
true η with radius of order n−1 ln n.
(c) The asymptotic normality of the stationary points of the aforementioned algorithm will be
substantiated under some mild restriction on the third derivative of R with respect to η.
Furthermore, the Appendix contains several easy-to-reference nonasymptotic results on the
global and local behaviour of the quadratic forms of Gaussian processes which may come in handy
for the analysis of certain problems in statistics and machine learning.
Plan of the paper.. In Section 2, we formulate ACS’s inversion-free estimation method and precisely
introduce the geometry of the sampling points. Section 3 expresses the necessary assumptions and
studies the asymptotic properties of the estimation algorithm. Section 3.1 presents the convergence rate of the global and local maximizers of the optimization problem introduced in Section
2. We investigate the minimax optimality and the asymptotic normality of the local maximizers
in Section 3.2. The objective of Section 4 is to assess the performance of ACS’s algorithm and
verify the developed theory using simulation studies on synthetic data. Section 5 serves as the
conclusion and discusses the future directions. Section 6 presents the proof of the main results.
Finally, the Appendix contains some auxiliary technicalities on the nonasymptotic behaviour of the
quadratic forms of Gaussian processes and of large covariance matrices with polynomially decaying
off-diagonal entries, which are essential in Section 6.

3

Notation.. For any m ∈ N, Im and 0m respectively denote the m by m identity matrix and all
zeros column vector of length m. Moreover, ∧ and ∨ stand for theP
minimum and maximum
operators. For two matrices of the same size M and M ′ , hM, M ′ i := i,j Mij Mij′ denotes their
usual inner product. We use the following matrix norms on M ∈ Rm×n . For any 1 ≤ p < ∞,
P
1/p
p
kM kℓp :=
|M
|
stands for the element-wise p−norm of M . The usual operator norm
ij
i,j

d
(largest singular value of M ) is represented by kM k2→2 . Associated to any
 ′ finite set ′ D ⊂ R and
′
′
s ∈ D, we define D−s := {s − s : s ∈ D}. Moreover, we write D (s, r) := s ∈ D : ks − skℓ2 ≤ r
and D c (s, r) = D \ D (s, r), for any non-negative r. S m stands for the m−dimensional unit sphere
with respect to the Euclidean norm, i.e., S m := v ∈ Rm+1 : kvkℓ2 = 1 . For a random sequence
xn and a deterministic positive sequence an , we write xn = OP (an ) when xn is bounded below
by an asymptotically, i.e., lim Pr (|xn | ≥ Can ) = 0 for some C > 0. For two sets Ω1 , Ω2 ⊂
n→∞

Rm , dist (Ω1 , Ω2 ) := inf ωi ∈Ωi , i=1,2 kω1 − ω2 kℓ2 represents their mutual distance with respect to the
Euclidean norm. Moreover, for A ⊂ Rm and r > 0, Nr (A) denotes a subset of A (of minimal size)
such that for any a ∈ A, dist ({a} , Nr (A)) ≤ r. The cardinality of such set will be referred as the
covering number of A. Given spatial points {s1 , . . . , sn } ∈ Rd and the covariance function R (·, η)
parametrized by η = (η1 , . . . , ηm ), the associated covariance matrix and its derivatives are defined
as

n
∂
∂
n
, ∀ r = 1, . . . , m.
Rn (η) =
R (si − sj , η)
Rn (η) = [R (si − sj , η)]i,j=1 ,
∂ηr
∂ηr
i,j=1
The higher order derivatives can be defined in an analogous way. For two random vectors v1 and v2 ,
d
the expression v1 = v2 means that they have the same distribution. Lastly, D (P1 k P2 ) indicates
the Kullback-Leibler divergence of two distributions Pi , i = 1, 2.
2. Problem set up and ACS’s estimation algorithm
Consider a mean zero and stationary (real valued) Gaussian process G : Rd 7→ R whose covariance function belongs to a parametric family CR,Ω := {R (·, η) : η ⊂ Ω}. In other words, there
exists η0 ∈ Ω for which


(2.1)
EG (s) G s′ = R s − s′ , η0 , ∀ s, s′ ∈ Rd .

Moreover, there is m ∈ N such that Ω is a compact (m + 1) dimensional subset of Rm+1 with
respect to the Euclidean topology. Thus, CR,Ω is assumed to be a finite dimensional class. For
analytical convenience, we consider an alternative formulation for the unknown parameters of the
covariance function given as
η0 = (φ0 , θ0 ) , φ0 ∈ I, θ0 ∈ Θ.
In this new representation, φ0 is a strictly positive scalar denoting the variance of G and the
m−dimensional vector θ0 stands for the other parameters of R. Moreover, I ⊂ (0, ∞) is a bounded
interval and Θ ⊂ Rm is compact. For instance in isotropic Matern or powered exponential classes,
θ0 is a positive vector representing the range parameter and fractal index. Finally, (2.1) can be
rewritten as R (s − s′ , η0 ) = φ0 K (s − s′ , θ0 ), in which K (·, θ0 ) indicates the correlation function
parametrized by θ0 .
The objective is to estimate η0 observing one realization of G at a deterministic set of spatial
locations Dn = {s1 , . . . , sn } ⊂ Rd . It is beneficial to emphasize that our asymptotic analysis lies
4

in the increasing domain regime in which the diameter of Dn tends to infinity as n → ∞. The
collected samples form a zero mean Gaussian column vector Y = [G (s1 ) , . . . , G (sn )]⊤ of length n.
Before proceeding further, let us precisely introduce the geometric structure of Dn .
Assumption 2.1. Suppose that there is N ∈ N such that n = N d . There exists δ ∈ [0, 1/2) for
which Dn is a d−dimensional δ−perturbed regular lattice (with unit grid size). Namely,
on
n
,
Dn = vi + δpi : vi ∈ VN,d , pi ∈ [−1, 1]d
i=1

in which VN,d := {v1 , . . . , vn } = {1, . . . , N }d denotes the d−dimensional regular lattice.

The condition δ ∈ [0, 1/2) guarantees the existence of a strictly positive minimum distance
(1 − 2δ) between the distinct points in Dn . The scalar δ quantifies the amount of irregularity in
Dn . In the case of δ = 0, Dn forms a regular lattice and the irregularity can be more apparent as δ
increases. Although the absence of randomness in Assumption 2.1 may appear problematic at first
sight, our theoretical contributions are not restricted to any further set of strong conditions on pi ’s.
For instance, the presented results in the next section hold almost surely if pi ’s are independent (or
even dependent) draws of a distribution supported on [−1, 1]d which is absolute continuous with
respect to the Lebesgue measure..
Now we present ACS’s estimation algorithm introduced in [2]. Define,


1
1
2
⊤
(2.2)
Y Rn (η) Y − kRn (η)kℓ2 .
η̂n = argmax Fn (Y, η) , where Fn (Y, η) :=
n
2
η∈Ω
Note that Fn (Y, η) does not depend on the Cholesky factorization
of Rn (η) and regardless of the

choice of covariance function, it can be evaluated in O n2 operations, even for irregularly spaced
samples which is an improvement over the conventional likelihood function. The optimization
algorithm in (2.2) can be reformulated as




1
φ2
2
⊤
kKn (θ)kℓ2 .
φY Kn (θ) Y −
φ̂n , θ̂n = argmax Fn (Y, φ, θ) , where Fn (Y, φ, θ) :=
n
2
(φ,θ)∈I×Θ
(2.3)
Despite the fact that Fn (Y, φ, θ) has a simple quadratic (concave) form of φ, its dependence to
θ is fairly complicated . For instance Fn is not a concave function of θ even for the classic case
of isotropic exponential covariance. So, accurate approximation of its global maximizer can be
computationally expensive.
Remark 2.1. We conclude this section mentioning two characteristics of Fn (Y, φ, θ) that can
provide a theoretical clue for generalizing ACS’s loss function to a broader class of inversion-free
losses. The first property is also critical for the theoretical analysis in the next section.
1. As stated in [2], the true parameter η0 is a stationary point of the expected value of Fn (Y, η).
That is,


∂
Fn (Y, η) |η=η0 = 0, ∀ j = 1, . . . , (m + 1) .
E
∂ηj
Roughly speaking, η0 is located in a small neighborhood of a stationary point of (2.2) if the
gradient of Fn (Y, η) is smooth enough and concentrated around its expected value. The next
fact, which has not been stated in [2], reveals a profound connection of (2.2) to MLE.
5

1/2

1/2

2. For brevity, define Hn (Y, η) := Fn (Y, η)−Fn (Y, η0 ) and Ln (η, η0 ) := Rn (η0 ) Rn−1 (η) Rn (η0 ).
Also, let η̃n denotes the MLE of η. Obvious calculations lead to
η̂n = argmax Hn (Y, η) ,
η∈Ω

η̃n = argmax Hn′ (Y, η) ,
η∈Ω

where Hn′ (Y, η) :=

o
1n
− log det Ln (η, η0 ) − n + Y ⊤ Rn−1 (η) Y .
n

Notice that EHn (Y, η0 ) = EHn′ (Y, η0 ) = 0. Under Assumption 2.1 and using similar techniques
as Appendix A, one can guarantee the existence of a universal scalar C ∈ (0, ∞) such that
EHn (Y, η) ≤ C EHn′ (Y, η) ,

∀η ∈ Ω.

(2.4)

Namely, in the increasing domain regime, the objective function proposed in [2] can be viewed
as an approximate minorizing surrogate of the likelihood function in the expected value sense
(it forms a perfect minorizer whenever C = 1 in (2.4)).
3. Main results
We establish the asymptotic characteristics of the estimation algorithm in (2.2). Section 3.1
examines the consistency of the global maximizer and the stationary points of (2.2) under some
sufficient conditions on Ω and the correlation function K (·, θ). The near minimax optimality and
the asymptotic normality of the stationary points will be covered in Section 3.2
3.1. Consistency and the convergence rate
The following assumptions are assumed on the parameter space Ω = I × Θ and the correlation
function K (·, θ) for studying the asymptotic behaviour of the global maximizer of (2.2). Similar
but slightly stronger conditions have been used in [4] for the increasing domain asymptotic analysis
of MLE.
Assumption 3.1. The following conditions are satisfied by Ω and K.
(A1) Θ and I are compact connected subsets of Rm and (0, ∞), respectively.
(A2) There are bounded scalars M > 0 and r1 > 1 such that for any s ∈ Dn ,


max
K s′ − s, θ2 − K s′ − s, θ1 ≥ M kθ2 − θ1 kℓ2 , ∀ θ1 , θ2 ∈ Θ.
s′ ∈Dn (s,r1 )

(A3) For some q ∈ {1, 2, 3}, there exists a positive scalar CK,Θ such that


CK,Θ
∂
∂
K (s, θ) ≤
...
,
max |K (s, θ)| ∨
θ∈Θ
∂θj1
∂θjq
1 + kskd+1
ℓ
2

(3.1)

∀ s ∈ Rm ,

for any j1 , . . . , jq ∈ {1, . . . , m}.
Condition (A2), assuring the identifiability of θ from the K, holds for the standard class of
correlation functions such as Matern, powered exponential and rational quadratic. A detailed look
at (A2) is postponed to the end of this section. Before commenting on (A3), let us define the family
of geometric anisotropic covariance functions.
6

Definition 3.1. Let G : Rd 7→ R be a zero mean stationary Gaussian process in Rd . Then G is
called geometric anisotropic if

q


(s − s′ )⊤ A0 (s − s′ ) , ∀ s, s′ ∈ Rd ,
(3.2)
R s − s′ , η0 := EG (s) G s′ = φ0 K

for φ0 > 0, symmetric positive definite matrix A0 ∈ Rd×d , η0 = (φ0 , A0 ) and a correlation function
K. Specifically if A0 = θ0−1 Id for some strictly positive θ0 , then G is said to be an isotropic Gaussian
process.

For geometric anisotropic processes, K is either assumed to be a fully known function (in this
case η0 = (φ0 , A0 ) in which φ0 ∈ I and A0 ∈ Θ, denotes the unknown parameters of covariance
function) or known up to some strictly positive scalar ν0 , usually refers to as the fractal index. In
the latter case, η0 = {φ0 , θ0 = (A0 , ν0 )}. Now, we mention some commonly used class of covariance
functions, with unknown fractal index, satisfying (A3) with q = 1 (appearing in the statement of
the first main result in this section). It is supposed in the following Remark that
Λmin,Θ ≤ min

A0 ∈Θ

1
A−1
0
2→2

≤ max kA0 k2→2 ≤ Λmax,Θ ,

(3.3)

A0 ∈Θ

for strictly positive and bounded scalars Λmin,Θ and Λmax,Θ . Namely, all eigenvalues of A0 are
universally bounded away from zero and infinity.
Remark 3.1. Any compactly supported correlation function, such as spherical or Wendland family
[23] on Rd trivially admits (A3). Assumption (A3) with q = 1 also holds for some classical families
of geometric anisotropic covariances such as:
(a) Matern: The Gaussian process G has Matern covariance function if it fulfills (3.2) with
K (r) =

21−ν0 ν0
r Kν0 (r) ,
Γ (ν0 )

(3.4)

in which ν0 is an unknown, strictly positive scalar lies in a compact space. Moreover, Γ (·) and
Kν0 (·) represent the Gamma function and the modified Bessel function of the second kind,
respectively. The parametric Matern family satisfies (A3) provided condition (3.3).
ν

(b) Powered exponential: A covariance function in this class satisfies (3.2) with K (r) = e−r 0 and
ν0 ∈ (0, 2). Like Matern class, assuming (3.3), any member of powered exponential family
fulfills (A3) with q = 1.
− d +ν
(c) Rational quadratic: The elements of this class are of the form (3.2) with K (r) = 1 + r 2 ( 2 0 )
and ν0 > 0. For the case of known fractal index, (A3) with q = 1 is valid, if (3.3) holds. Note
that for unknown ν0 the exact same statement is satisfied under a slightly stronger condition
of ν0 > 1/2
Parts (b) and (c) of Remark 3.1 are verifiable by straightforward algebra and differentiation
rules. In order to demonstrate part (a), see [1] for the derivative properties of the Bessel function
of the second kind (with respect to the entries of A0 ) and see Lemma A.9 for the asymptotic
behaviour of the partial derivatives of the Matern covariance with respect to ν0 . Now, we state
the first significant result of this section regarding the consistency of the global maximizer of (2.3)
under Assumption 3.1 and perturbed regular lattice sampling.
7

Theorem 3.1. Suppose that Assumptions 2.1 and 3.1 with q = 1 hold for Dn , Ω and K. Then
the maximizer of (2.3) satisfies
!
r
ln n
φ̂n
Pr
→ 0, as n → ∞,
(3.5)
∨
θ̂n − θ0
−1 ≥C
φ0
n
ℓ2
for some constant C (which depends on Dn , Ω and K).
Remark 3.2. Let φmin and φmax denote the smallest and largest element in I. Obviously, φmin
and φmax are well defined and finite due to (A1). Moreover,
!
!
p
φ̂
φ̂n
n
∨
θ̂n − θ0
(1 ∧ φmin )
∨
θ̂n − θ0
− 1 ≤ kη̂n − η0 kℓ2 ≤ 1 + φ2max
−1 .
φ0
φ0
ℓ2
ℓ2
Thus (3.5) is a stronger statement than kη̂n − η0 kℓ2 = OP
when φmin > 0 (which is true under A1).

√


n−1 ln n and they are equivalent

An analogous consistency result has been proved recently by Bachoc [4] for the MLE and cross
validation estimator. Based upon Theorem 3.1, the asymptotic rate of ACS’s algorithm has not
been sacrificed for increasing the speed and memory efficiency comparing to the MLE.
Finally, we concisely address the role of the identifiability condition (A2) in Theorem 3.1.
Actually, (A2) plays a decisive role in translating consistent estimation of the correlation matrix
(in the relative sense) to η0 . Strictly speaking, (A2) is required to deduce (3.5) from the probabilistic
statement

√
1
√ Kn (θ̂n ) − Kn (θ0 )
= OP
n−1 ln n .
n
ℓ2

The rest of this section is devoted to the analysis of the stationary points of (2.3). Solving the
unique root of the derivative of Fn (Y, φ, θ) with respect to φ, yields a closed form formula for φ̂n
in terms of θ̂n , Y and the correlation function, namely
φ̂n =

Y ⊤ Kn (θ̂n )Y
Kn (θ̂n )

2

.

(3.6)

ℓ2

Moreover, θ̂n can be obtained using
θ̂n = argmax Gn (Y, θ) ,
θ∈Θ

where Gn (Y, θ) =

Y ⊤ Kn (θ) Y
.
kKn (θ)kℓ2

(3.7)

Note that for large n, computing the global maximizer of (3.7) can be less intensive than (2.2)
due to searching over a smaller space Θ. We first visually assess the key properties of Fn in
some simple scenarios. In Figure 1, Gn (Y, θ) (which is a univariate function of scalar θ) has been
plotted versus θ for the two dimensional (d = 2) isotropic Matern covariance function in two
different scenarios. In the left panel the isotropic Gaussian process G has been generated with the
parameters (φ0 , θ0 , ν0 ) = (1, 4, 0.5) and has been sampled in a randomly perturbed regular lattice
with δ = 0.2 and of size N = 100. In the right panel, the covariance parameters are given by
(φ0 , θ0 , ν0 ) = (1, 6, 1.5) and the GP is sampled at a randomly generated perturbed regular lattice
8

θ0 = 4, ν0 = 0.5, δ = 0.3

θ0 = 6, ν0 = 1.5, δ = 0.2

3.5

12.0
10.0

2.5

Gn√
(Y,θ)
n

Gn√
(Y,θ)
n

3.0

2.0
1.5

8.0
6.0
4.0

1.0
0

2

4
θ

6

8

2.0

0

2

4

6
θ

8

10

12

Figure 1: The above figures exhibit n−1/2 Gn (Y, θ) for the isotropic Matern covariance function
(with known ν0 ). In the left panel, θ0 = 4, ν0 = 0.5 and the spatial samples form a two dimensional
randomly perturbed regular lattice of size N = 100 with δ = 0.3. In the right panel, θ0 = 6,
ν0 = 1.5 and Dn is a randomly chosen two dimensional perturbed regular lattice with N = 100 and
δ = 0.3.
with N = 100 and δ = 0.2. As is apparent from Figure 1, for these two parsimonious scenarios
Gn (Y, θ) is not a concave function of θ and has a single inflection point. However, Gn (Y, θ) has
only one stationary point which coincides with its global maximizer. In the following, we rigorously
study the large sample behaviour of the stationary points of Gn (Y, θ) (as well as Fn ) in a generic
case. We initiate our analysis by stating the sufficient conditions on K and Ω.
Assumption 3.2. (A1) holds for Ω, and K fulfills (A2) and (A3) with q = 2 in Assumption 3.1.
Remark 3.3. The analysis of the stationary points of (2.3) requires a slightly stronger conditions
than that of the global maximizer in Assumption 3.1. The main distinction is the polynomial decay
of the second order derivative of K with respect to θ. Note that the new condition on the second
derivative of K is not too restrictive. For instance the same analysis as Remark 3.1 validates this
condition for all covariance families introduced in Remark 3.1 (with a larger constant CK,Θ ).
Theorem 3.2. Suppose that Dn admits Assumption 2.1 and Assumption 3.2 holds for Ω and K.
Then any stationary point of the optimization problem (2.3) satisfies
!
r
ln n
φ̂n
= 0,
(3.8)
−1 ≥C
∨
θ̂n − θ0
lim Pr
n→∞
φ0
n
ℓ2
for an appropriately chosen constant C > 0 depending on Dn , Ω and K.
Theorem 3.2 shows that any stationary point of Fn is concentrated in a small neighborhood of
(φ0 , θ0 ), with high probability. In other words, Fn (Y, φ, θ) shows a similar behaviour as Figure 1
in the general case. In addition, the comparison between (3.5) and (3.8) reveals that stationary
points converge to (φ0 , θ0 ) with the same rate as the global maximizer.
9

We conclude this section by illustrating how restrictive the identifiability assumption (A2) can
be for the frequently used classes of the covariance functions. We first introduce a slightly stronger
identifiability condition than (A2), which will be referred to as (A4). Note that a slightly modified
version of (A4) has been first introduced in [4] for studying the increasing domain asymptotics of the
maximum likelihood and cross validation algorithms. That is to say, these identifiability conditions
are not exclusive to ACS method and pop up in the asymptotic analysis of other algorithms.
Proposition 3.1. (A2) is satisfied, whenever
(A4) (a) There are positive scalars r2 > 1 and M2 such that for any η ∈ Ω and λ ∈ S m ,
min

max

s∈Dn s′ ∈Dn (s,r2 )

m+1
X

λj

j=1


∂
R s − s′ , η ≥ M2 .
∂ηj

(b) The following inequality holds for any distinct pair of points η1 , η2 ∈ Ω


min max
R s − s′ , η2 − R s − s′ , η1 > 0.
s∈Dn s′ ∈Dn (s,r2 )

Clearly, (A4.b) is necessary for any algorithm consistently estimating η and it can be verified
for all typical classes of geometrical anisotropic covariance function. However, understanding the
role and restrictiveness of (A4.a) is more subtle than that of (A4.b). Note that unlike (A3), all the
introduced identifiability conditions not only depend on the choice of the covariance function but
also to the observed locations Dn . It may be excessive to seek the class of covariances satisfying
(A4.a) for any perturbed lattice Dn . So, a more pertinent question is: which class of covariance
functions do almost surely satisfy (A4.a) for a randomly generated perturbed lattice? The following result responds to question by rigorously characterizing a broad subclass of the geometrically
anisotropic covariances (as defined in Definition 3.1) fulfilling (A4.a).
Proposition 3.2. Let P be a distribution in [−1, 1]d which is absolutely continuous with respect to
the Lebesgue measure. Suppose that R (·, η) : Rd 7→ [0, ∞) is a geometrically anisotropic covariance
function with a known ν0 (if exists). Then, (A4) almost surely holds if
(a) K : [0, ∞) 7→ [0, ∞) is a nonzero, differentiable and strictly decreasing function (K may only
have right derivative at zero).
(b) Dn is a randomly generated δ−perturbed lattice associated to P. That is, pi are independent
draws of P in Assumption 2.1.
Corollary 3.1. Let Dn be d−dimensional regular lattice (associated to δ = 0) and assume that
R (·, η) : Rd 7→ [0, ∞) is a geometrically anisotropic covariance function with known ν0 . Then, (A4)
holds if K : [0, ∞) 7→ [0, ∞) admits the condition (a) in Proposition 3.2.
Although the conditions of Proposition 3.2 trivially hold for the non-compactly supported covariance function introduced in Remark 3.1, deploying analogous proof techniques can lead to a
similar result for compactly supported covariance function.
Proposition 3.3. Suppose that P, R (·, η) and Dn satisfy the same conditions as Proposition 3.2.
Then, (A4) almost surely holds if there exists a large enough positive scalar r0 for which
10

• K : [0, ∞) 7→ [0, ∞) is a nonzero, differentiable and strictly decreasing function in the interval
[0, r0 ] and K (r) = 0 for any r > r0 .
Although the required conditions on the covariance function’s formulation, in Propositions 3.2
and 3.3, are very minimal, we assume that the fractal index ν0 (if exists) is fully known. However
in the following result, ν0 is one of the unknown parameters to be estimated. Here the central
emphasis is on the powered exponential and rational quadratic classes, as their partial derivative
with respect to the fractal index have a somewhat simple closed form that can be handled without
great difficulty.
Proposition 3.4. Let P be a distribution in [−1, 1]d which is absolutely continuous with respect to
the Lebesgue measure. Suppose that R (·, η) : Rd 7→ [0, ∞) is a geometrically anisotropic covariance
function. Then, (A4) almost surely holds if
(a) K : [0, ∞) 7→ [0, ∞) is either a powered exponential or rational quadratic covariance functions
in Remark 3.1 with unknown ν0 > 0.
(b) Dn is a randomly generated δ−perturbed lattice associated to P. That is, pi are independent
draws of P in Assumption 2.1.
Remark 3.4. A prudent look at the proof of Proposition 3.4 reveals that the following property
(which is satisfied by the powered exponential and rational quadratic families) has the crucial role.
∂K
∂ν

{r=0d , θ}

= 0,

∀ θ ∈ Θ.

(3.9)

For the Matern class, not only ∂K
∂ν not satisfy (3.9), it does not have a tractable algebraic form.
We believe that (A4) holds true for the geometric anisotropic Matern family with unknown ν, even
though it is beyond the reach of our current proof technique.
3.2. Minimax optimality and Asymptotic normality
Now, we further investigate the asymptotic statistical properties of ACS’s algorithm. Near
minimax optimality and asymptotic normality are respectively presented in Theorems 3.3 and 3.4.
Theorem 3.3. Suppose that Assumptions 2.1 and 3.1 hold for Dn , Ω and K. Then there exist
n0 ∈ N and a bounded scalar C > 0 such that


C
1
sup Pr kη̂n − η0 kℓ2 ≥ √
≥ ,
n
4
η0 ∈Ω
for any estimator η̂n and any n ≥ n0 .
√ Theorem 3.3 reveals that the established bounds in Theorems 3.1 and 3.2 are sharp up to order
ln n. This means that for the perturbed regular lattice sampling scheme, no algorithm can achieve
a significantly better rate than the one considered in this paper.
Theorem 3.4. Suppose that Dn is a perturbed lattice introduced in Assumption 2.1. Furthermore,
(A1), (A3) with q = 3 and (A4) are fulfilled by Ω and R. There is a positive definite matrix
Σ ∈ R(m+1)×(m+1) with bounded operator norm such that
√
d
n (η̂n − η0 ) → N (0m+1 , Σ) .
11

(3.10)

The exact formulation of Σ has been omitted in this section due to its complicated algebraic
form. We refer the reader to the proof of Theorem 3.4 in Section 6 for further details. It is
worthwhile to mention that the entries of Σ heavily depend on the configuration of points in Dn ,
which is a major disparity between fixed and increasing domain asymptotics. Comparing Theorems
3.2 and 3.4, here we impose a slightly stronger differentiability condition (polynomially decaying
of the third derivative) for establishing asymptotic normality. This condition has been formerly
introduced in [4] and holds for the geometrically anisotropic covariances in Remark 3.1.
4. Simulation studies
The relatively large-scaled numerical studies in this section give a fairly comprehensive appraisal
of the statistical and computational performance of the optimization problem (2.3). Despite the
popularity of R language among the statisticians, running the iterative programs such as loops
in R is much slower than that of C++ (around 250 times slower according to some studies [3]).
Taking advantage of the Rcpp package and hybrid programming techniques in R can considerably
expedite the execution time (up to 50 times in our simulation studies). In order to get the maximum
speed, the open MP application programming interface has been used to exploit the multi-threaded
programing technology. All the numerical experiments in this section have been executed on 12
processors, except for the second simulation study (n = 106 ) which has been implemented on 60
cores.
Generating high dimensional samples from a Gaussian process on an irregularly spaced grid
is the foremost challenge that we confronted in our synthetic data simulations. Applying the
traditional method based upon the Cholesky decomposition of the covariance matrix is almost
infeasible in the case that n ≈ 105 or larger. Hence, we use the considerably faster spectral
method (pp. 203 − 205, [9]) for generating stationary Gaussian processes. For completeness, this
algorithm will be concisely presented here. Strictly speaking, the objective is to simulate a real
valued zero mean stationary Gaussian process G in Rd with the covariance function φ0 K (·, θ0 ) over
a δ−perturbed lattice Dn = {s1 , . . . , sn }. For the purpose of generating a realization of G on a
perturbed grid, without loss of generality we can assume that the samples are all of unit variance,
i.e., φ0 = 1. We also assume that G is geometric anisotropic. Recalling from Definition 3.2, there
d×d which represents the symmetric square root of
is a symmetric positive definite matrix
 B0 ∈ R
A0 , such that K (r, θ0 ) = K kB0 rkℓ2 . Throughout this section d = 2 and K is either the Matern
or rational quadratic covariance function which have been previously introduced in Remark 3.1.
Let p ∈ N be a large enough number and {ξk }pk=1 be i.i.d. uniform random variables on [−π, π].
Let K̂ : Rd 7→ R denotes the spectral density of G defined by
Z
Z

−d
−d
:=
K kB0 rkℓ2 cos (hω, ri) dr.
K (r, θ0 ) cos (hω, ri) dr = (2π)
K̂ (ω)
(2π)
Rd

Rd

The non-negative mapping K̂ (·) is a density function in Rd (since it integrates to K (0, θ0 ) = 1).
Furthermore, let {ωk }N
k=1 be independent draws from the density K̂ (·). Now, define
r p
2X
cos (hωk , si + ξk ) , ∀ s ∈ Rd .
(4.1)
G (s) =
p
k=1


It is known that G is a geometric anisotropic process with cov (G (s) , G (s′ )) = K kB0 (s − s′ )kℓ2
for any pair s, s′ ∈ Rd (p. 204, [9]), converging in distribution to a Gaussian process as p tends to
12

infinity. Next, we explain how to generate the random variables {ωk }pk=1 . The following fact which
can be proved using the integration by substitution plays a principal role in our algorithm.
Remark 4.1. Let ω ′ ∈ Rd be a draw from the following density function
Z

K krkℓ2 cos (hu, ri) dr.
K̂I (u) = (2π)−d
Rd

d

Then ω and B0 ω ′ have the same distribution, i.e., ω = B0 ω ′ . Note that K̂I is an isotropic function.

d
Namely, there is a function Φ : R 7→ [0, ∞) for which K̂I (u) = Φ kukℓ2 . Moreover ω ′ = r kψψddk in
ℓ2

which ψd is a standard d− dimensional Gaussian vector and r is a non-negative random variable
with the density function
fr (r) =

2π d/2 d−1
dπ d/2
r d−1 Φ (r) =
r Φ (r) .
Γ (d/2 + 1)
Γ (d/2)
d

For instance in the case of d = 2, we have fr (r) = 2πrΦ (r). Hence, ω =

r
kψd kℓ

B0 ψd .
2

For Matern covariance function in two dimensional plane (d = 2), generating independent
samples of the random variable r is a straightforward task. In this case
fr (r) =

−(ν+d/2)
2π d/2 d−1 π −d/2 Γ (d/2 + ν)
2rν
2π d/2 d−1
.
r Φ (r) =
r
1 + r2
=
Γ (d/2)
Γ (d/2)
Γ (ν)
(1 + r 2 )1+ν

Thus the cumulative distribution is of the form Pr (r ≤ r) := Fr (r) = 1 − 1 + r 2
d

r=

Fr−1 (u)

=

q

1 − (1 − u)−1/ν ,

−ν

. So

in which u is a uniform random variable in [0, 1]. One can find a closed from expression
for r in

terms of u for the rational quadratic covariance function, in the case that τ + 12 ∈ N (Recall τ
from Remark 3.1). In this case, Φ (·) has a form of the Matern covariance function (3.4) (with
different constants) due to the duality principle of the Fourier transform.
Throughout this section, G is assumed to be a zero mean Gaussian process in R2 , whose
covariance function is a member of either Matern or rational quadratic families with a known
fractal index. In the first experiment, G is an isotropic spatial process. In other words, we set
A0 = θ0−2 I2 in the Definition 3.2. θ0 is a strictly positive scalar known as the range parameter.
Furthermore, Dn is a randomly generated δ−perturbed lattice of size 3202 , i.e., n = 102400 ≈ 105 .
The approximated realizations of G are generated using (4.1) with p = 1.5 × 105 . To investigate the
role of spatial irregularity in the computational and statistical performance of ACS’s algorithm, we
vary δ in √
the set {0.1, 0.3}. The range parameter and the standard deviation, which is represented
by σ0 = φ0 , are respectively estimated solving the optimization problem (3.7) and closed form
formula (3.6). The range parameter space is chosen as Θ = [0.1, 15]. The single variable constrained
optimization problem (3.7) is solved using the optimize function in R ,which exploits a combination
of golden section search and successive parabolic interpolation. We stop the iteration of the solver
when the relative change in the objective is below 10−3 . Table 1 displays the summary of our first
simulation study.
13

ν0 = 0.5
Matern

ν0 = 1.5
ν0 = 2.5
ν0 = 0.5

Rational quadratic
ν0 = 1.5

η̂
η̂
η̂
η̂
η̂

δ = 0.1
η0 = (1, 4)
η0 = (1, 7)
= (0.993, 4.420) η̂ = (0.978, 7.565)
η0 = (1, 4)
η0 = (1, 7)
= (0.973, 3.618) η̂ = (1.023, 6.568)
η0 = (1, 4)
η0 = (1, 7)
= (1.014, 4.186) η̂ = (1.010, 7.428)
η0 = (1, 4)
η0 = (1, 7)
= (1.017, 4.100) η̂ = (0.976, 7.415)
η0 = (1, 4)
η0 = (1, 7)
= (1.000, 4.001) η̂ = (1.014, 7.210)

η̂
η̂
η̂
η̂
η̂

δ = 0.3
η0 = (1, 4)
η0 = (1, 7)
= (1.005, 3.941) η̂ = (1.028, 6.553)
η0 = (1, 4)
η0 = (1, 7)
= (1.026, 4.343) η̂ = (1.005, 7.102)
η0 = (1, 4)
η0 = (1, 7)
= (1.044, 4.037) η̂ = (0.993, 6.505)
η0 = (1, 4)
η0 = (1, 7)
= (1.013, 3.916) η̂ = (1.001, 6.639)
η0 = (1, 4)
η0 = (1, 7)
= (1.017, 3.920) η̂ = (0.994, 7.063)

Table 1: Estimation of η0 = (σ0 , θ0 ) for the isotropic Matern and rational quadratic covariance
functions, where Dn is a perturbed lattice of size 3202 with associated δ ∈ {0.1, 0.3}.
The required CPU times for the numerical experiments in Table 1 are approximately 30 and 60
minutes for the rational quadratic and Matern kernels, respectively. However evaluating the full
MLE for a such large sample size is intractable. As is apparent from Table 1, the estimated
parameters, η̂, are in a close neighborhood of η0 . Moreover, estimating σ0 has a significantly higher
precision than that of the range parameters, since as the distant samples in Dn carry negligible
information about θ0 . Lastly, the condition number of the covariance matrix increases with the
value of range parameter, leading to a higher estimation error kη0 − η̂kℓ2 for larger θ0 . In the second
simulation study which has the same set up as the first experiment, Dn is a irregular grid of size
10002 = 106 . We also set p = 5 × 105 in (4.1). Table 2 encapsulates the results of this experiment.
The evaluation of η̂ for this very high dimensional numerical study takes 8 hours on 60 cores with
4GB RAM.

Rational quadratic

δ = 0.1, ν0 = 0.5
η0 = (1, 4)
η̂ = (1.004, 4.053)

δ = 0.3, ν0 = 1.5
η0 = (1, 7)
η̂ = (0.999, 6.948)

Table 2: Estimation of η0 = (σ0 , θ0 ) for the isotropic rational quadratic covariance functions, where
Dn is a perturbed lattice of size 10002 with associated δ ∈ {0.1, 0.3}.
In the next set of experiments featuring Gaussian processes with isotropic covariance functions, we
set Dn to be a two dimensional perturbed lattice of size 1002 . For such a scenario, η̂ can be estimated
in a few minutes. Thus, we simulated T = 100 independent realizations and η0 is estimated using
the same procedure as previous studies, for each realization. The mean and root mean squared
error (RMSE) have been computed across T experiments. Table 3 displays the average and RMSE
for the standard deviation and range parameters for different values of η, δ and covariance kernels.
The instances for which η̂ hits the boundary points of Θ have been excluded in the procedure of
calculating the mean and RMSE of the estimates.
Looking at the left and right panels of the Table 3 reveals that the RMSE of σ̂ and θ̂ are slightly
larger for the higher values of δ. Moreover, as we have discussed before, the RMSE and the average
norm of (η̂ − η0 ) directly depends on the range parameter. It is immediately clear that there is a
considerable reduction in RMSE for rational quadratic kernel comparing to the Matern class. This
observation may look surprising for the reader, as the condition number of the covariance matrix
associated to Matern kernel is smaller than that of rational quadratic due to its faster decay.
Thus at the first glance, it may not corroborate our developed theory regarding the consistency
14

θ0 = 4, ρ0 = 6, ν0 = 0.5, δ = 0.1

θ0 = 3, ρ0 = 7, ν0 = 0.5, δ = 0.3

4
4

3
Gn (Y,θ)
√
n

Gn√
(Y,θ)
n

5

2
5
θ

6
2

0

10

5
ρ

0

4
10

5

2
00

θ

ρ

Figure 2: The above figures exhibit n−1/2 Gn (Y, θ) for geometric anisotropic Matern covariance
function with ν0 = 0.5. The spatial samples form a two dimensional randomly δ−perturbed regular
lattice of size N = 100. In the left panel, (θ0 , ρ0 ) = (4, 6) and δ = 0.1. In the right panel,
(θ0 , ρ0 ) = (3, 7) and δ = 0.3.
of ACS’s estimation algorithm for covariance matrices with bounded condition number. However,
obtaining a highly accurate estimate of the dependence parameters are more difficult for a rapidly
decaying covariance function, as more samples are almost independent. In the extreme case θ0 is
unidentifiable if K (·, θ0 ) is a compactly supported covariance function whose support size is strictly
less than (1 − 2δ) (In this case all the samples are independent).
Matern covariance
(ν0 = 0.5)
Matern covariance
(ν0 = 1.5)
Rational quadratic
covariance (ν0 = 0.5)
Rational quadratic
covariance (ν0 = 1.5)

δ = 0.1
(σ0 , θ0 ) = (1, 4)
(σ0 , θ0 ) = (1, 7)
θ̂ ± RSME = 4.107 ± 1.224 θ̂ ± RSME = 7.259 ± 2.462
σ̂ ± RSME = 0.999 ± 0.067 σ̂ ± RSME = 0.991 ± 0.089
(σ0 , θ0 ) = (1, 4)
(σ0 , θ0 ) = (1, 7)
θ̂ ± RSME = 3.936 ± 1.127 θ̂ ± RSME = 6.588 ± 2.060
σ̂ ± RSME = 1.002 ± 0.070 σ̂ ± RSME = 0.992 ± 0.096
(σ0 , θ0 ) = (1, 4)
(σ0 , θ0 ) = (1, 7)
θ̂ ± RSME = 3.889 ± 0.599 θ̂ ± RSME = 6.855 ± 1.507
σ̂ ± RSME = 1.002 ± 0.062 σ̂ ± RSME = 0.986 ± 0.082
(σ0 , θ0 ) = (1, 4)
(σ0 , θ0 ) = (1, 7)
θ̂ ± RSME = 3.984 ± 0.342 θ̂ ± RSME = 7.160 ± 1.010
σ̂ ± RSME = 0.999 ± 0.028 σ̂ ± RSME = 0.994 ± 0.074

δ = 0.3
(σ0 , θ0 ) = (1, 4)
(σ0 , θ0 ) = (1, 7)
θ̂ ± RSME = 3.982 ± 0.980 θ̂ ± RSME = 6.814 ± 2.233
σ̂ ± RSME = 0.995 ± 0.062 σ̂ ± RSME = 1.003 ± 0.096
(σ0 , θ0 ) = (1, 4)
(σ0 , θ0 ) = (1, 7)
θ̂ ± RSME = 4.180 ± 1.181 θ̂ ± RSME = 6.519 ± 2.127
σ̂ ± RSME = 0.995 ± 0.072 σ̂ ± RSME = 1.018 ± 0.107
(σ0 , θ0 ) = (1, 4)
(σ0 , θ0 ) = (1, 7)
θ̂ ± RSME = 4.032 ± 0.647 θ̂ ± RSME = 6.793 ± 1.373
σ̂ ± RSME = 0.992 ± 0.046 σ̂ ± RSME = 0.990 ± 0.069
(σ0 , θ0 ) = (1, 4)
(σ0 , θ0 ) = (1, 7)
θ̂ ± RSME = 4.016 ± 0.348 θ̂ ± RSME = 7.127 ± 1.116
σ̂ ± RSME = 0.994 ± 0.026 σ̂ ± RSME = 0.995 ± 0.049

Table 3: Mean and RMSE of η̂ over 100 independent experiments for the isotropic Matern and
rational quadratic covariance functions, where Dn is a perturbed lattice of size 1002 with associated
δ ∈ {0.1, 0.3}.
Now we turn to investigate the precision and RMSE of estimation algorithm (2.3) for the geometric
anisotropic covariance structure. Same as before, G is a zero mean stationary Gaussian process in
R2 observed on a perturbed lattice of size 1002 . G has a geometric anisotropic covariance kernel

15

(Matern or rational quadratic) with

 −1
θ0
0
,
B0 =
0 ρ−1
0

θ0 = 4, and ρ0 = 6.

The parameter space Θ = {(θ0 , ρ0 ) ∈ Θ} is a two dimensional box chosen as [0.1, 15]2 . The range
parameters θ0 and ρ0 are estimated by solving the optimization problem (3.7). The Figure 2 exhibits
the objective function in (3.7) and its contours for a Gaussian process with geometric anisotropic
Matern covariance with ν0 = 0.5 which has been sampled on a δ−perturbed regular grid. In the left
and right panels of Figure 2, the other parameters are respectively given by (θ0 , ρ0 , δ) = (4, 6, 0.1)
and (θ0 , ρ0 , δ) = (3, 7, 0.3). It can be seen from Figure 2 that Gn is a unimodal function with
only one stationary point. We perform the maximization using the optim function in R and with
L-BFGS-B algorithm [6] (box constrained BFGS). The maximum iteration and the initial guess of
the L-BFGS-B method are respectively 100 and (2, 2). The components of the gradient function are
computed using the finite difference approximation with the step size of 10−3 . We cease the iteration
when the relative change in the objective function is below 10−5 . The computation procedure of
the average and RMSE of η̂ = (θ̂, ρ̂, σ̂) is exactly the same as the former simulation study. Table
4 presents a summary of the final results of this simulation study. It is clear from Table 4 that
the RMSE for Matern covariance is significantly larger in comparison to rational quadratic class.
Furthermore, increasing ν0 for each covariance kernel leads to a slightly larger RMSE.

Matern covariance (ν0 = 0.5)

Matern covariance (ν0 = 1.5)

Rational quadratic covariance (ν0 = 0.5)

Rational quadratic covariance (ν0 = 1.5)

δ = 0.1
(σ0 , ρ0 , θ0 ) = (1, 6, 4)
σ̂ ± RSME = 0.988 ± 0.096
ρ̂ ± RSME = 6.042 ± 1.885
θ̂ ± RSME = 4.091 ± 1.110
(σ0 , ρ0 , θ0 ) = (1, 6, 4)
σ̂ ± RSME = 0.993 ± 0.108
ρ̂ ± RSME = 05.965 ± 1.981
θ̂ ± RSME = 3.740 ± 1.146
(σ0 , ρ0 , θ0 ) = (1, 6, 4)
σ̂ ± RSME = 0.992 ± 0.071
ρ̂ ± RSME = 5.978 ± 1.241
θ̂ ± RSME = 4.092 ± 0.843
(σ0 , ρ0 , θ0 ) = (1, 6, 4)
σ̂ ± RSME = 0.996 ± 0.036
ρ̂ ± RSME = 6.116 ± 0.821
θ̂ ± RSME = 4.045 ± 0.543

δ = 0.3
(σ0 , ρ0 , θ0 ) = (1, 6, 4)
σ̂ ± RSME = 0.993 ± 0.097
ρ̂ ± RSME = 6.478 ± 1.908
θ̂ ± RSME = 4.038 ± 1.272
(σ0 , ρ0 , θ0 ) = (1, 6, 4)
σ̂ ± RSME = 0.984 ± 0.104
ρ̂ ± RSME = 6.160 ± 1.890
θ̂ ± RSME = 3.970 ± 1.243
(σ0 , ρ0 , θ0 ) = (1, 6′ 4)
σ̂ ± RSME = 0.989 ± 0.076
ρ̂ ± RSME = 5.921 ± 1.208
θ̂ ± RSME = 4.037 ± 1.064
(σ0 , ρ0 , θ0 ) = (1, 6, 4)
σ̂ ± RSME = 0.998 ± 0.036
ρ̂ ± RSME = 6.158 ± 0.766
θ̂ ± RSME = 4.150 ± 0.524

Table 4: Mean and RMSE of η̂ over 100 independent experiments for the geometric anisotropic
Matern and rational quadratic covariance functions, where Dn is a perturbed lattice of size 1002
with associated δ ∈ {0.1, 0.3}.
5. Discussion
Investigation of the asymptotic properties of the non-likelihood based optimization algorithms for
estimating covariance parameters has remained relatively intact. To our knowledge, this paper is the
first asymptotic analysis of ACS’s algorithm in the increasing domain regime. Notwithstanding the
thorough study of the consistency, minimax optimality and asymptotic normality of the stationary
16

points of ACS’s loss function, there is much future work to be done to determine the computational
and statistical strengths and weaknesses of this algorithm in either of the two frequently used
asymptotic regimes. Here we mention a few among the many future directions which were beyond
the scope of this paper.
(a) As indicated in Remark 2.1, the inversion-free loss function can be viewed as an approximate
minorizer for the likelihood loss (in the expected value sense) in the increasing domain setting.
However, more work needs to be done to know how to precisely characterize a rich class of
minorizers for the likelihood loss. We believe that responding to this question will provide a
flexible class of fast and consistent estimators of covariance parameters.
(b) Spatial statisticians usually cast doubt upon the benefits of increasing domain asymptotics as
spatial processes are unlikely to be stationary over a large domain. However the developed
theory in this manuscript has persuaded us that inversion free algorithm can consistently
estimate microergodic covariance parameters, when applied on the preconditioned samples of
a Gaussian random field in a fixed domain. In this regard, a modified version of the estimator
studied in this paper can provide a powerful tool for interpolation of Gaussian processes.
6. Proofs of the main results
We first introduce a few notation to simplify the algebra in the forthcoming sections. For any
strictly positive scalar r and any θ0 ∈ Θ, the ball of radius r (with respect to the Euclidean norm)
centered at θ0 and its complement are defined by

Θθ0 (r) := θ ∈ Θ : kθ − θ0 kℓ2 ≤ r , Θcθ0 (r) := Θ \ Θθ0 (r) .
Furthermore, for any θ1 , θ2 ∈ Θ define
1/2

Mθ1 ,θ2 :=

1/2

Kn (θ1 ) Kn (θ2 ) Kn (θ1 )
,
kKn (θ2 )k2ℓ2

Hθ1 ,θ2 := kKn (θ2 )kℓ2 Mθ1 ,θ2 .

(6.1)

Proof of Theorem 3.1. Our proof has two major parts. In the first part, the consistency of θ̂n
(correlation function’s parameters) will be substantiated. In the second part, we establish the
consistency of φ̂n , which has a closed form solution in terms of Y , correlation function and θ̂n ,
by conditioning on the consistency of θ̂n . To this end, various types of concentration inequalities
regarding the quadratic forms (and their supremum over a bounded space) of Gaussian processes
are of the indispensable importance. Such results will be presented in the Appendix A.
√
1/2
Let Z be a standard Gaussian vector in Rn . As Y and φ0 Kn (θ0 )Z have the same distribution,
(2.3) can be equivalently written by




φ2
2
⊤ 1/2
1/2
φ̂n , θ̂n = argmax φφ0 Z Kn (θ0 )Kn (θ) Kn (θ0 )Z −
kKn (θ)kℓ2 .
(6.2)
2
(φ,θ)∈I×Θ
The objective function in (6.2) is quadratic in terms of φ and its maximizer φ̂n has a simple closed
form. Replacing φ̂n to (6.2) gives a surrogate form for θ̂n . Omitting the cumbersome algebra, the
final results are given by
φ̂n
= Z ⊤ Mθ0 ,θ̂n Z,
φ0

θ̂n = argmax Z ⊤ Hθ0 ,θ Z.
θ∈Θ

17

(6.3)

We first show (as Claim 1) the consistency of θ̂n , which is the supremum of a generalized chi-square
random variable. The purpose of Claim 2 is to find the estimation rate of φ0 .
q
Claim 1. Choose ξ > (m − 1) and let rn := Cmin lnnn for some bounded positive scalar Cmin (See
Lemma A.5 for its exact form). Then,
n
o
Pr θ̂n ∈ Θcθ0 (rn ) → 0, as n → ∞.

√
Proof of Claim 1. Consider the sequence rn′ = n−1 ln n, ∀n. The boundedness of Θ guarantees
the existence of some R0 > 0 such that a ball of radius R0 contains Θ. So, the classical volume
argument implies that
 m

R0
c
= o (nm ) .
(6.4)
Nrn′ Θθ0 (rn ) ≤ Nrn′ (Θ) .
rn′

So, there is n1 ∈ N such that Nrn′ Θcθ0 (rn ) ≤ nm for any n ≥ n1 . It follows from (6.3) that


n
o
Pr θ̂n ∈ Θcθ0 (rn ) ≤ RHS := Pr (An ) := Pr Z ⊤ Hθ0 ,θ0 Z ≤ sup Z ⊤ Hθ0 ,θ Z  .
θ∈Θcθ (rn )
0

Thus, it suffices to control RHS from above. For a properly chosen positive scalar C2 , define the
event πn by


r



ln n 
⊤
⊤
πn :=
Z Hθ0 ,θ Z + C2
.
sup
sup Z Hθ0 ,θ Z ≤



n 

θ∈Θcθ0 (rn )
θ∈Nr ′ Θc (rn )
θ0

n

√

√
Notice that rn n = O
n−1 ln n . According to Lemma A.2, there is a bounded C2 > 0 for
which τn := Pr (πnc ) → 0. We refer the reader to Lemma A.2 for the closed form expression of C2 .
An upper bound on RHS is obtained by conditioning An on πn .
RHS

=
(A)

≤

(B)

≤

Pr (An ∩ πn ) + τn Pr (An | πnc ) ≤ τn + Pr (An ∩ πn )



τn + Pr Z ⊤ Hθ0 ,θ0 Z ≤
τn + n m

sup




θ∈Nr ′ Θcθ (rn )
n
0

θ∈Nr ′

n

Z ⊤ Hθ0 ,θ Z + C2

sup


Θcθ (rn )
0



r



ln n 

n

Pr Z ⊤ Hθ0 ,θ0 Z ≤ Z ⊤ Hθ0 ,θ Z + C2

r

ln n
n

!

(6.5)

The way that πn and An have been defined trivially justifies inequality (A). Furthermore (B) is
inferred from the combination of (6.4) and
 the union bound. Applying Lemma A.5 guarantees the
following result for any θ ∈ Nrn′ Θcθ0 (rn ) .
⊤

⊤

Pr Z Hθ0 ,θ0 Z ≤ Z Hθ0 ,θ Z + C2
18

r

ln n
n

!

≤ n−(1+ξ) .

(6.6)

Finally, substituting (6.6) into (6.5) yields that


RHS ≤ τn + nm−(1+ξ) → 0,

as n → ∞.

Claim 2. There exists a bounded scalar C > 0, depending on Dn , K and Ω, such that
!
r
ln
n
φ̂
n
πn′ := Pr
→ 0, as n → ∞.
− 1 ≥ rn′′ := C
φ0
n
√
(Here C = 2 (Dmax Cmin + C ′ Λmax m) in which C ′ is a large enough positive universal constant
and Cmin has been defined in Claim 1. Λmax and Dmax are given in the Proposition A.1.)
Proof of Claim 2. Recall rn and rn′ form Claim 1. Obviously,
!
)
(
r

n
o
ln n \ 
φ̂n
′
c
θ̂n ∈ Θθ0 (rn )
.
−1 ≥C
πn ≤ T1 + T2 := Pr θ̂n ∈ Θθ0 (rn ) + Pr
φ0
n
Since T1 tends to zero (via Claim 1), it suffices to show that T2 is a diminishing sequence as
n → ∞. Let βθ̂n be the closest point in Nrn′ (Θθ0 (rn )) (which is a deterministic set) to θ̂n . Based
upon identity (6.3), we have
φ̂n
(6.7)
− 1 = Z ⊤ Mθ0 ,θ̂n Z − 1
φ0
Given that θ̂n belongs to Θθ0 (rn ), applying the triangle inequality on the right hand side of the
identity (6.7) yields that, almost surely
φ̂n
−1
φ0

≤
(D)





+ tr Mθ0 ,βθ̂ − 1
Z ⊤ Mθ0 ,θ̂n Z − Z ⊤ Mθ0 ,βθ̂ Z + Z ⊤ Mθ0 ,βθ̂ Z − tr Mθ0 ,βθ̂

≤

θ∈Θθ0 (rn )

:=

sup

sup

θ∈Θθ0 (rn )

n

n

n

n

n

|tr (Mθ0 ,βθ − 1)| + Z ⊤ Mθ0 ,βθ Z − tr (Mθ0 ,βθ ) + Z ⊤ Mθ0 ,θ Z − Z ⊤ Mθ0 ,βθ Z

n
o
T21 (θ) + T22 (θ) + T23 (θ) .

(6.8)

Replacing the random quantities θ̂n and βθ̂n with the nonrandom parameters θ and βθ is the key
advantage of (D). Now we control the terms T21 , T22 and T23 from above, uniformly over Θθ0 (rn ).
Lemma A.3 guarantees the existence of a scalar C0 , for which
!
lim Pr

n→∞

sup

θ∈Θθ0 (rn )

T23 (θ) ≤ C0 rn′

→ 1.

(6.9)

C0 depends on Λmax and Dmax . See Lemma A.3 for its exact formulation. For large enough n, we
have
r

 1
√
ln n
′
−1
C0 rn = O n
.
(6.10)
ln n < Dmax Cmin
2
n
19

o

Now we control T21 (θ) uniformly from above using Proposition A.1. The goal is to show that
r
ln n
3
sup T21 (θ) < Dmax Cmin
.
(6.11)
2
n
θ∈Θθ (rn )
0

Applying the Cauchy-Schwartz inequality shows that (recalling Mθ1 ,θ2 from (6.1))
sup

T21 (θ) =

θ∈Θθ0 (rn )

sup
θ∈Θθ0 (rn )

≤

hKn (βθ ) , Kn (βθ ) − Kn (θ0 )i
hKn (βθ ) , Kn (θ0 )i
− 1 = sup
2
kKn (βθ )kℓ2
kKn (βθ )k2ℓ2
θ∈Θθ0 (rn )

kKn (βθ ) − Kn (θ0 )kℓ2
kKn (βθ ) − Kn (θ0 )kℓ2
√
≤ sup
.
kKn (βθ )kℓ2
n
θ∈Θθ (rn )
(rn )

sup
θ∈Θθ0

(6.12)

0

Furthermore, using the part (b) of the Proposition A.1, we get
kKn (βθ ) − Kn (θ0 )kℓ2
√
n
(rn )

sup
θ∈Θθ0

≤ Dmax
≤ Dmax
≤

sup
θ∈Θθ0 (rn )

kθ0 − βθ kℓ2

sup
θ∈Θθ0 (rn )

Dmax rn + rn′



kθ − βθ kℓ2 + kθ0 − θkℓ2
3
< Dmax Cmin
2

r

ln n
.
n


(6.13)

Note that the last inequality holds for large enough n. So (6.10) follows from replacing (6.13)
into (6.12).
In the sequel we achieve a uniform upper bound on T22 . For brevity define un :=
√
Λmax mn−1 ln n and select a large enough universal constant C ′ . Recall that βθ , by its definition,
is an element of the finite set Nrn′ (Θθ0 (rn )). Thus,
!

′
sup Pr T22 (θ) ≥ C ′ un .
Pr
sup T22 (θ) ≥ C un ≤ Nrn′ (Θθ0 (rn ))
θ∈Θθ0 (rn )

θ∈Θθ0 (rn )

The same trick as (6.4) leads to Nrn′ (Θθ0 (rn )) = o (nm ). So, it is adequate to show that

Pr T22 (θ) ≥ C ′ un ≤ n−m , ∀ θ ∈ Θθ0 (rn ) .

(6.14)

We employ Hanson-Wright inequality (Theorem 1.1, [17]) for obtaining a probabilistic upper bound
on T22 (θ) (for a fixed θ).

o
n
√
√
≤ n−m , ∀ θ ∈ Θθ0 (rn ) .
Pr T22 (θ) ≥ C ′ m ln n kMθ0 ,βθ kℓ2 ∨ kMθ0 ,βθ k2→2 m ln n
For simplifying the upper bound on T22 (θ), we control the operator and Frobenius norms of Mθ0 ,βθ
from above. The following inequalities can be easily justified by Proposition A.1.
kMθ0 ,βθ kℓ2 ≤ Λmax n−1/2 ,

kMθ0 ,βθ k2→2 ≤ Λ2max n−1 .

Replacing (6.15) into Hanson-Wright inequality justifies (6.14). Hence
!
Pr

sup

θ∈Θθ0 (rn )

T22 (θ) ≥ C ′ un

≤ Nrn′ (Θθ0 (rn )) n−m → 0,

as n → ∞.

(6.15)

(6.16)

Substituting inequalities (6.9), (6.10), (6.11) and (6.16) into (6.8) concludes the proof by confirming
that T2 goes to zero as n → ∞.
20

Combining Claims 1 and 2 ends the proof.
√
Proof of Theorem 3.2. Let rn = C n−1 ln n for some strictly positive C whose exact form will be
given shortly. Let (φ̂n , θ̂n ) be any stationary point of the optimization problem (2.3). Due to the
space constraint, we just show that θ̂n ∈ Θθ0 (rn ). The same technique as Claim 2 in the proof
of Theorem 3.1 attains the convergence rate of φ̂n . Notice that θ̂n is a stationary point of the
optimization problem (6.3). We show that with a high probability there is no θ ∈ Θcθ0 (rn ) for
which the gradient of the objective function in (6.3) be exactly zero. In order to substantiate our
claim, we prove that the absolute value of the inner product of the gradient and a fixed non-zero
vector is uniformly greater than zero on Θcθ0 (rn ).
Let us first give the closed form of the gradient function in (6.3), which will be denoted by [Gl (θ)]m
l=1 .
A few lines of direct algebra leads to
∂ ⊤
Z Hθ0 ,θ Z
∂θl
(
)
∂
K
(θ)
∂
n
Z ⊤ Kn1/2 (θ0 )
Kn (θ) − h
Kn (θ) , Kn (θ)i
Kn1/2 (θ0 ) Z,
∂θl
∂θl
kKn (θ)k2ℓ

Gl (θ) := Z ⊤ Pθl0 ,θ Z :=
=

2

where Z ∈ Rn is a standard Gaussian vector and θl , l = 1, . . . , m is the lth component of θ.
Trivially, Gl (θ̂n ) = 0 for any l = 1, . . . , m. Choose an arbitrary unit norm vector λ ∈ Rm and define
1/2
Y := Kn (θ0 ) Z. Observe that


m
m
m
X
X
X
∂
Kn (θ) 
∂
λj
Kn (θ) − h
Kn (θ) , Kn (θ)i
λj Gj (θ) = Y ⊤
λj
W (θ) :=
2  Y.

∂θj
∂θj
kK
(θ)k
n
ℓ2
j=1
j=1
j=1
(6.17)

We can conclude that θ̂n ∈ Θθ0 (rn ) in probability, if we can prove that
!
Pr

inf

θ∈Θcθ (rn )
0

|W (θ)| > 0

→ 1,

as n → ∞.

(6.18)

(6.18) is immediate after we show the two following claims.
Claim 1. There exists a positive finite constant C0 (depending on K, Θ and Dn ) such that


m


X
√
lim Pr sup W (θ) −
(6.19)
λj tr Pθj0 ,θ ≥ C0 n ln n = 0.
n→∞

θ∈Θ

j=1

Claim 2. The succeeding inequality holds for large enough n (C0 is from the previous claim).
inf
c

θ∈Θθ (rn )
0

m
X
j=1



√
λj tr Pθj0 ,θ > C0 n ln n

The Claim 1 provides a uniform concentration inequality regarding the random function W (θ). In
the Claim 2, we obtain a uniform lower bound on the expected value of W (θ) over Θcθ0 (rn ).
21

Proof of Claim 1. For simplicity, define


m
m
X
X
∂
Kn (θ) 
∂
.
λj
Kn (θ) − h
Kn (θ) , Kn (θ)i
λj
Tn (θ) :=

∂θj
∂θj
kKn (θ)k2ℓ 
j=1

j=1

(6.20)

2

The matrix Tn (θ) has appeared in the defining formula for W (θ) in (6.17). The basic facts of the
directional derivative yield
(
)
K
(θ)
n
Tn (θ) = lim γ −1 Kn (θ + γλ) − Kn (θ) − hKn (θ + γλ) − Kn (θ) , Kn (θ)i
γ→0
kKn (θ)k2ℓ2
(
)
Kn (θ)
−1
Kn (θ + γλ) − hKn (θ + γλ) , Kn (θ)i
= lim γ
.
(6.21)
γ→0
kKn (θ)k2ℓ2
Proposition A.3 is the essential tool in our proof. As inequality (6.19) is the same as (A.26),
it suffices to justify the three conditions of Proposition A.3. We verify conditions (a) and (b)
simultaneously. Choose θ, θ ′ ∈ Θ in an arbitrary way. The following inequalities are obvious
implications of Proposition A.1 (Λmax , Λmin and Dmin are introduced there).
m
X

λj Pθj0 ,θ

j=1

m
X
j=1

m
X

≤ Λmax kTn (θ)k2→2 ,
2→2

λj Pθj0 ,θ′ −

m
X
j=1

λj Pθj0 ,θ

j=1

λj Pθj0 ,θ

2→2

m
X


≤ Λmax Tn θ ′ − Tn (θ)

2→2
−1

λj Pθj0 ,θ

j=1

≤
ℓ2

2→2

,

Λmax kTn (θ)k2→2
.
Λmin kTn (θ)kℓ2

Thus, it is sufficient to justify the conditions of Proposition A.3 for {Tn (θ)}θ∈Θ . According to
Lemma A.6, Tn (θ) admits conditions (a) and (b). As the operator√norm of Tn (θ) is uniformly
bounded over Θ, condition (c) holds if kTn (θ)kℓ2 grows faster than ln n. We aim to show that
√
kTn (θ)kℓ2 = O ( n).


2 1/2

hK
(θ)
,
K
(θ
+
γλ)i
(A)
n
n
kTn (θ)kℓ2 = lim γ −1 kKn (θ + γλ)k2ℓ2 −


γց0
kKn (θ)kℓ2
=

≥
≥
(B)

≥


2 1/2




kKn (θ + γλ)kℓ2
Kn (θ + γλ)
Kn (θ)
√
,
i
2−2 h

γց0
kKn (θ)kℓ2 kKn (θ + γλ)kℓ2 
γ 2
)1/2
(
kKn (θ + γλ)kℓ2
Kn (θ)
Kn (θ + γλ)
√
2−2 h
,
i
lim
γց0
kKn (θ)kℓ2 kKn (θ + γλ)kℓ2
γ 2
 √ −1
kKn (θ + γλ)kℓ2
lim γ 2
Kn (θ)
Kn (θ + γλ) −
γց0
kKn (θ)kℓ2
ℓ2
r
(C)
kKn (θ + γλ) − Kn (θ)kℓ2
Dmin n
√
lim
.
≥
γց0
Λmax 2
Λmax γ 2
lim

22

(6.22)

Trivial calculations can prove identity (A). Moreover, (B) and (C) are immediate consequences of
Lemma A.4 and Proposition A.1, respectively.
Proof of Claim 2. For notational convenience, define vn (θ) := kKn (θ)k−1
ℓ2 Kn (θ) for any θ ∈ Θ.
√
−1
Also, recall that rn = C n ln n for some C to be chosen. We first derive a lower bound on
|E {W (θ)}| in terms of vn vectors. Regrouping the terms in (6.21) leads to
|E {W (θ)}| = |hKn (θ0 ) , Tn (θ)i|
kKn (θ + γλ)kℓ2
= lim
hvn (θ0 ), vn (θ + γλ)i − hvn (θ + γλ) , vn (θ)ihvn (θ) , vn (θ0 )i
γ→0 γ kKn (θ0 )k−1
ℓ2
≥ n lim

γ→0

hvn (θ0 ), vn (θ + γλ)i − hvn (θ + γλ) , vn (θ)ihvn (θ) , vn (θ0 )i
.
γ

Hence, it suffices to show that
A :=

inf
c

lim

θ∈Θθ (rn ) γ→0
0

>

rn′ := C0

r

hvn (θ0 ), vn (θ + γλ)i − hvn (θ + γλ) , vn (θ)ihvn (θ) , vn (θ0 )i
γ

ln n
.
n

(6.23)

Define two unit norm vectors u1 , u2 ∈ Rm and positive scalars τ1 , τ2 by
vn (θ0 ) = hvn (θ0 ), vn (θ)ivn (θ) + τ1 u1 ,

vn (θ + γλ) = hvn (θ + γλ), vn (θ)ivn (θ) + τ2 u2 .
Using the fact that both vectors u1 and u2 are perpendicular to vn (θ), we get
hvn (θ0 ), vn (θ + γλ)i − hvn (θ + γλ) , vn (θ)ihvn (θ) , vn (θ0 )i
τ1 τ2
=
hu1 , u2 i.
γ
γ

(6.24)

Achieving tight lower bounds on τ1 and τ2 is a crucial step of the proof. Without loss of generality
we can assume that hvn (θ0 ), vn (θ)i is non-negative. Thus, for any θ ∈ Θcθ0 (rn ).
τ1

1 p
1 − hvn (θ0 ), vn (θ)i2 ≥ √
2 − 2hvn (θ0 ), vn (θ)i
2
kKn (θ)kℓ2
1
Kn (θ0 )
=√
Kn (θ) −
kKn (θ0 )kℓ2
2 kKn (θ)kℓ

=

vn (θ0 ) − hvn (θ0 ), vn (θ)ivn (θ) =

=

1
√ kvn (θ0 ) − vn (θ)kℓ2
2

(D)

≥

p

2

kKn (θ) − Kn (θ0 )kℓ2 (E) Dmin
Dmin
√
kθ − θ0 kℓ2 ≥ √
rn .
≥ √
2
2 kKn (θ)kℓ2 Λmax
2Λmax
2Λ2max

ℓ2

(6.25)

In above inequality, (D) and (E) has been concluded from Lemma A.4 and Proposition A.1.
Applying similar tricks as (6.25) leads to
Dmin
Dmin
τ2 ≥ √
kθ + γλ − θkℓ2 = γ √
.
2
2Λmax
2Λ2max

23

(6.26)

Replacing the two inequalities (6.25) and (6.26) into (6.24) shows that
hvn (θ0 ), vn (θ + γλ)i − hvn (θ + γλ) , vn (θ)ihvn (θ) , vn (θ0 )i
≥C
γ



D
√ min
2Λ2max

2

|hu1 , u2 i|

r

ln n
.
n

Using Assumption 3.2, one can show that hu1 , u2 i is nonzero, i.e., |hu1 , u2 i| ≥ α0 for some α0 > 0,
√ −1 2 2
2Dmin Λmax implies that
universally over θ ∈ Θ and θ0 ∈ Θcθ (rn ). Finally, choosing C > C0 /α0
r

2 r
Dmin
ln n
ln n
′
A ≥ Cα0 √
> rn = C 0
,
2
n
n
2Λmax
which ends the proof.

Proof of Theorem 3.3. We follow the standard techniques presented in the Chapter 2 of [21] for
bounding the minimax risk from below. For any θ ∈ Θ, Pθ stands for the associated distribution
to a zero mean Gaussian vector with the covariance function Kn (θ). Finding far enough (with
respect to the Euclidean distance) pair of the correlation parameters, θi ∈ Θ, i = 1, 2, for which
D (Pθ1 k Pθ2 ) < α = 1/2 lies at the heart of our proof. The two bounded positive scalars Dmax and
min√
(Choose
Λmin appearing here are defined in Proposition A.1. To ease notation let rn := 8DΛmax
n
n large enough so that 4rn ≤ diam (Θ)). Choose θ1 , θ2 ∈ Θ with 2rn ≤ kθ2 − θ1 kℓ2 ≤ 4rn . The
connectedness of Θ guarantees the existence of such pair of points. We first use the Proposition
A.4 to show that D (Pθ1 k Pθ2 ) ≤ α.
D (Pθ1 k Pθ2 ) ≤ 2n



Dmax
kθ2 − θ1 kℓ2
Λmin

2

≤ 32n



Dmax
rn
Λmin

2

1
=α= .
2

As α ≥ D (Pθ1 k Pθ2 ), Theorem 2.2 of [21] yields
inf sup Pr
θ̂n θ0 ∈Θ



θ̂n − θ0

ℓ2

≥ rn



≥




1 −α
e
∨
4

p !
1 − α2
1
= .
2
4

The desired statement follows from the fact that kη̂n − η0 kℓ2 ≥ θ̂n − θ0

ℓ2

(6.27)

.

Proof of Proposition 3.1. Recall our notation from Assumptions 3.1 and 3.2. The main objective
is to verify the inequality (3.1) (Assumption (A2)). Let φ be any point in I and let φmax :=
sup {φ′ : φ′ ∈ I}. We also choose r1 = r2 = r. Select small enough δ > 0, which its value will be
chosen later and an arbitrary s ∈ Dn . First, choose two distinct θ1 , θ2 ∈ Θ whose mutual Euclidean
distance is at least δ and define ηi := (θi , φ) , i = 1, 2. Clearly,




= φ K s′ − s, θ1 − K s′ − s, θ2
R s′ − s, η1 − R s′ − s, η2


≤ φmax K s′ − s, θ1 − K s′ − s, θ2 ,
for any s′ ∈ Dn (s, r). Thus, according (A4.b),
max

h∈Dn (s,r)−s

|K (h, θ1 ) − K (h, θ2 )| ≥ CΩ,K := φ−1
max
24

max

h∈Dn (s,r)−s

|R (h, η1 ) − R (h, η2 )| > 0.

(6.28)

In other words, (A2) holds by choosing M = δ−1 CΩ,K . So, we need to verify condition (3.1) when
the distance θ1 and θ2 is smaller than δ. In this case, there is small τ < δ and λ0 ∈ S m−1 for which
θ2 − θ1 = τ λ0 . Moreover, define the unit norm vector λ ∈ Rm+1 by λ = (λ0 , 0). Because of the
continuity of K, we can choose δ small enough for which there exists M ′ ∈ (0, 1) such that
|K (h, θ1 ) − K (h, θ2 )|
kθ2 − θ1 kℓ2
h∈Dn (s,r)−s
max

≥ M′

max

lim

h∈Dn (s,r)−s r→0

= ϑ :=

M′
φmax

K (h, θ1 + rλ) − K (h, θ1 )
r

max

h∈Dn (s,r)−s

m+1
X
j=1

λj

∂
R (h, η)
∂ηj

η=η1

.

′
:=
Condition (A4.a) assures the existence of a positive universal constant M2 , for which ϑ ≥ CΩ,K
′
′
−1
M2 M /φmax > 0. Hence, (A2) holds by taking M = CΩ,K ∧ δ CΩ,K .
d×d and h ∈ Rd , the norm khk repreProof of Proposition
A
√ 3.2. For a positive definite matrix A ∈ R

sents the quantity h⊤ Ah. Furthermore, define Br (Dn ) := h = s′ − s : s, s′ ∈ Dn , khkℓ2 ≤ r , for
any r > 0. Any geometrical anisotropic covariance function associated to the correlation function
K are of the form



R s′ − s, η = φK s′ − s A , ∀ s, s′ ∈ Rd .

Notice that η = (φ, A) is a vector of size m + 1 = d2 + 1. Consider any λ ∈ S m and η ∈ Ω. The
chain rule for derivative leads to the following identity for any η ∈ Ω and h ∈ Rd
J (h) :=

m+1
X
j=1

d
X
K ′ (khkA )
∂
λij hi hj
,
R (h, η) = λ1 K (khkA ) + φ
λj
∂ηj
2 khkA
i,j=1

h 6= 0d .

(6.29)

Slightly abusing the notation, λij denotes the entry of λ corresponding to Aij . Furthermore, let
J˜ (h) denotes the second term in the right hand side of (6.29) and define Λ := [λij ]di,j=1 . One can
easily show that lim J˜ (h) = 0. Thus, identity (6.29) can be rewritten as
h→0d

J (0d ) = λ1 K (0) = λ1 ,
J (h) = λ1 K (khkA ) + φK ′ (khkA )

h⊤ Λh
,
2 khkA

h 6= 0d .

(6.30)

It suffices to prove that J is not a zero function of h over Br (Dn ) for some appropriately chosen r.
Assume toward contradiction that J is a zero function over Br (Dn ). As 0d ∈ Br (Dn ), λ1 equals
zero from the first identity in (6.30). Therefore, J˜ (h) = 0 for any h ∈ Br (Dn ). Since K is a
differentiable and strictly decreasing function in (0, ∞), we have K ′ (u) 6= 0 for any strictly positive
u. Thus,
J˜ (h) = 0, ∀ h ∈ Br (Dn ) ⇒ hΛ, hh⊤ i = 0, h ∈ Br (Dn ) .
(6.31)
Now choose r in such a way that |Br (Dn )| > (d + 1). Due to the non-atomicity of P, there are
almost surely d non-zero and linearly independent vectors, {h1 , . . . , hd }, in Br (Dn ). According to
(6.31), hΛ, hi h⊤
i i = 0 for i = 1, . . . , d. Hence Λ is a zero matrix, which leads to a clear contradiction
to the fact that λ ∈ S m (recall that kλk2ℓ2 = λ21 + kΛk2ℓ2 ).

25

Proof of Proposition 3.4. As the proof is akin to that of Proposition 3.2, the details are dropped for
avoiding redundancy. We also keep the notation used there, for simplicity. Notice that η = (φ, θ)
and here θ = (ν, A) is of size m = d2 + 1. Choose any unit norm λ ∈ S m . We first consider the
rational quadratic case. It is easy to verify the following identity for any h ∈ Br (Dn ).


d
m+1




X
X
d/2 + ν
λ1
∂
λ
h
h
.
R (h, η) = φKν (khkA )
− λ2 ln 1 + khk2A −
λj
J (h) :=
ij
i
j

φ
∂ηj
1 + khk2A
i,j=1

j=1

It suffices to show that J is a nonzero function on Br (Dn ) unless all λ coefficients are zero. Assume
that the premise is not true and J vanishes in Br (Dn ). Thus, J (0d ) = λ1 = 0 and so J˜ defined by
!


d/2
+
ν
2
J˜ (h) := λ2 ln 1 + khkA +
h⊤ Λh,
2
1 + khkA
should be a zero function in Br (Dn ). If λ2 = 0, then the rest of the proof is similar to the argument
of Proposition 3.2. So assume that λ2 is nonzero. Hence for C = φ/λ2 (d/2 + ν) we have

 

1 + khk2A ln 1 + khk2A = C khk2Λ , ∀ h ∈ Br (Dn ) ,

which is almost surely absurd for large enough r and so terminates the proof. Apart from the
algebraic form of J, our proof for the powered exponential family goes through without any change.
In the sequel, we just express J and skip the further details for page constraints.
λ1
ν
J (h)
−(2−ν)
khk2Λ ,
=
− λ2 khkνA ln (khkA ) − khkA
φKν (khkA )
φ
2

J (0d ) = λ1 .

Proof of Theorem 3.4. Let g : Ω 7→ Rm+1 represents the gradient of the objective function in (2.2)
with respect to η. Here gj , j = 11, . . . , (m + 1) stands for the j th entry of g. Analyzing the exact
√
second order Taylor expansion of ng (η) around η0 at η = η̂n is the integral part of the proof. We
argue that the second order term of the expansion, which involves the third order derivatives of the
covariance function, converges to zero in probability as n grows to infinity. We also show that the
first term (zeroth order term) in the expansion converges weakly to a Gaussian random variable.
These two ingredients lead to the desirable result by showing the asymptotic normality of the first
√
order term in the expansion, which directly depends on n (η̂n − η0 ).
∂Rn (η)
for any q ∈ {1, 2} and J ∈ {1, . . . , m + 1}q . In fact
For simplicity, define RnJ (η) = ∂η
j ∂ηj
q

1

ngj (η) = Y ⊤ Rnj (η) Y − hRnj (η) , Rn (η)i.

(6.32)

Let η̂n be an arbitrary stationary point of optimization problem (2.2). Clearly, g (η̂n ) = 0m+1 . The
second order approximation of gj around η̂n yields
√
√
√
√
i + n∆j (η0 , η̂n ) ,
ngj (η̂n ) = ngj (η0 ) + h n (η̂n − η0 ) , ∇η gj (η)
η=η0

for some residual function ∆j (·, ·). Note that ∆j (η0 , η̂n ) is given by

m+1
∂gj (η)
⊤
∆j (η0 , η̂n ) = (η̂n − η0 )
(η̂n − η0 )
∂ηl1 ∂ηl2 η=zj l1 ,l2 =1
26

in which zj lies on the line segment between η0 and η̂n . Proposition D.10 of [4] guarantees the
−1
statement (3.10) for Σ = Σ−1
2 Σ1 Σ2 and hence concludes the proof, if
(a) The matrix Σ2 defined as the following, is well defined and positive definite.

m+1
∂
Pr
n
→ Σ2 , as n → ∞.
gk (η)
V := −
∂ηl
η=η0 l,k=1
√

d

ng (η0 ) → N (0m+1 , Σ1 ) for some positive semidefinite matrix Σ1 ∈ R(m+1)×(m+1) .


√
(c) Pr lim n∆j (η0 , η̂n ) = 0 = 1, for any j ∈ {1, . . . , m + 1}.

(b)

n→∞

The reminder of the proof hinges on the following technicalities which verify conditions (a) − (c).
Validating condition (a). The entries of V n has the following explicit form.
o
1 n ⊤ lk
1
Y Rn (η0 ) Y − hRnlk (η0 ) , Rn (η0 )i .
Vlkn = hRnl (η0 ) , Rnk (η0 )i +
n
n
Now define

m+1
hRnl (η0 ) , Rnk (η0 )i
Σ2 := lim
.
n→∞
n
l,k=1
Notice that the entries of Σ2 are well defined and bounded due to part (a) of Proposition A.2. The
proof will be presented
steps: First we show that Σ2 is a positive definite matrix. Second,
 in⊤ two
lk
n
we prove that Φlk := Y Rn (η0 ) Y − hRnlk (η0 ) , Rn (η0 )i /n converges to zero in probability. To
substantiate the first claim, consider an arbitrary λ ∈ S m . It is required to show that λ⊤ Σ2 λ > c,
for some constant c > 0. The condition (A4.a) guarantees the existence of positive scalars M2 , r2
such that for any s ∈ Dn ,
max

s′ ∈Dn (s,r2 )

m+1
X
l=1

Thus,
⊤

λ Σ2 λ =

lim

n→∞

m+1
X

l,k=1

λl


∂
R s − s′ , η
∂ηl

η=η0

≥ M2 .

(6.33)

m+1
1 X
hRnl (η0 ) , Rnk (η0 )i
= lim
λl Rnl (η0 )
λl λk
n→∞ n
n

"m+1
X

∂
1
λl
R s′ − s, η
= lim
n→∞ n
∂ηl
l=1

η=η0

l=1
2

#

2

ℓ2

(A)

s,s′ ∈Dn ℓ2

≥ M22 .

(6.34)

Here, inequality (A) is an easy consequence of (6.33). The rest of the proof is devoted to prove
the second claim. Choose an arbitrary strictly positive ǫ. As Φnlk is a zero mean random variable,
using Chebyshev’s inequality we get
Pr (|Φnlk |

≥ ǫ)

≤
(C)

=

var (Φnlk )
ǫ2
O n


−1

1/2

1/2

2 Rn (η0 ) Rnlk (η0 ) Rn (η0 )
=

n 2 ǫ2

→ 0,
27

2
ℓ2

(B)

.



φ0
Rnlk (η0 )
nǫ

ℓ2

2

in which (B) and (C) are implied by Propositions A.1 and A.2, respectively (See Appendix A
for more details about the constants).
1/2
1/2
Validating condition (b). Define Qjn := n−1/2 Rn (η0 ) Rnj (η0 ) Rn (η0 ) for 1 ≤ j ≤ m + 1, and
for any λ = (λ1 , . . . , λm+1 ) ∈ S m . Rewriting (6.32) yields
write Ψn,λ := λ1 Q1n + . . . + λm+1 Qm+1
n

√
d
ngj (η0 ) = Z ⊤ Qjn Z − tr Qjn .

√
The asymptotic normality of ng (η0 ) is justified if there is a positive semi-definite Σ1 such that

√
d
hλ, ng (η0 )i → N 0, λ⊤ Σ1 λ for any λ ∈ S m . This statement trivially holds for zero Ψn,λ . So,
without loss of generality assume that Ψn,λ is non-zero. Observe that
 ⊤
√
Z Ψn,λ Z − tr (Ψn,λ )
hλ, ng (η0 )i =
kΨn,λ kℓ2 .
kΨn,λ kℓ2
We claim that lim 2 kΨn,λ k2ℓ2 = λ⊤ Σ1 λ for a covariance matrix Σ1 . The construction of Ψn,λ yields
n→∞

h
im+1
λ.
2 kΨn,λ k2ℓ2 = 2λ⊤ hQkn , Qln i

(6.35)

l,k=1

m+1

Thus, it is enough to show that the matrix defined by Σ1 := lim 2 hQkn , Qln i l,k=1 is well defined
n→∞

(with bounded entries) and positive semi-definite. Well definiteness of Σ1 can be proved using
the same techniques as the proof of Claim 1 and by employing Propositions A.1 and A.2. The
positive semi-definite property of Σ1 is an immediate consequence of (6.35). We conclude the proof
by showing that
 ⊤
Z Ψn,λ Z − tr (Ψn,λ ) d
√
→ N (0, 1) .
2 kΨn,λ kℓ2
According to Lemma A.8, this statement is valid if the following claim holds.
√
Claim 1. kΨn,λ k−1
ℓ2 kΨn,λ k2→2 ≤ C/ n for some positive scalar C.
Proof of Claim 1. We show that C depends on m, Λmax , Λmin and Λ′max (Except m, all the constant
are introduced in Propositions A.1 and A.2). Obviously Ψn,λ can be rewritten as,


m+1


X
1
Ψn,λ = √ Rn1/2 (η0 )
λj Rnj (η0 ) Rn1/2 (η0 )


n
j=1

Applying Propositions A.1, we get

Pm+1
j
m+1
kRn (η0 )k2→2
j=1 λj Rn (η0 )
kΨn,λ k2→2
Λmax X
2→2
≤
λj Rnj (η0 )
≤ P
j
m+1
kΨn,λ kℓ2
Λ
min
λmin {Rn (η0 )}
j=1
j=1 λj Rn (η0 )
ℓ2

2→2

m+1
X

−1

λj Rnj (η0 )

j=1

Furthermore, using Proposition A.2 leads to
m+1
X
j=1

λj Rnj (η0 )

≤
2→2

m+1
X
j=1

|λj | Rnj (η0 )
28

2→2

√
≤ Λ′max kλkℓ1 ≤ Λ′max m + 1.

ℓ2

(6.36)

From (6.34) we know that there is a scalar C0 ∈ (0, ∞) for which
Replacing the last two inequalities into (6.36) ends the proof.

Pm+1
j=1

λj Rnj (η0 )

ℓ2

√
≥ C0 n.

In conclusion we state that the condition (c) can be proved using akin techniques as the proof of
Proposition 3.2 of [4]. We omit the technical details due to the space constraints.

A

Appendix: Auxiliary results

We first state an elementary lemma regarding the basic properties of the operator norm. The proof
is skipped since it is straightforward.
Lemma A.1. For any matrix A ∈ Rn×n , define its absolute value by |A| := [|Aij |]ni,j=1 . Then,
kAk2→2 ≤ k|A|k2→2 . Moreover, the largest eigenvalue of |A| can be written as
k|A|k2→2 = max v ⊤ |A| v,
n−1
v∈S+

n−1
denotes the collection of unit norm vectors with non-negative entries in Rn .
in which S+

The next result examines the perturbation of some norms of Kn (θ) with respect to θ. It appears
to be of great importance for proving Theorems 3.1–3.4 in the Section 6.
Proposition A.1. Suppose that Dn admits Assumption 2.1. Moreover, Assumption 3.1 holds for
Θ and K. Construct n × n correlation matrix Kn (θ) := [K (s − s′ , θ)]s,s′ ∈Dn for any θ ∈ Θ.
(a) There are bounded positive scalars Λmin and Λmax (depending on K, Θ, d and δ) such that
Λmin ≤ min min
n∈N θ∈Θ

1

,
Kn−1 (θ) 2→2

max max kKn (θ)k2→2 ≤ Λmax .
n∈N θ∈Θ

(b) There exist scalars Dmin , Dmax ∈ (0, ∞) (depending on K, Θ, d and δ) such that

and

kKn (θ2 ) − Kn (θ1 )k2→2 ≤ Dmax kθ2 − θ1 kℓ2 ,

(A.1)

1
√ kKn (θ2 ) − Kn (θ1 )kℓ2 ≤ Dmax kθ2 − θ1 kℓ2 ,
n

(A.2)

1
√ kKn (θ2 ) − Kn (θ1 )kℓ2 ≥ Dmin kθ2 − θ1 kℓ2 ,
n

(A.3)

for any θ1 , θ2 ∈ Θ.
Proof. Part (a) is similar to Proposition D.4 and Lemma D.5 of [4] and can be substantiated in an
analogous way. So we skip its proof due to the space limit. Furthermore, Eq. (A.2) is an immediate
consequence of Eq. (A.1). Thus we only need to prove inequalities (A.1) and (A.3). For brevity,
define two n × n matrices A and B by
Aij = |K (si − sj , θ2 ) − K (si − sj , θ1 )| ,
29


−1
Bij = 1 + ksi − sj kd+1
.
ℓ2

For distinct θ1 , θ2 ∈ Θ, we have
kKn (θ2 ) − Kn (θ1 )k2→2
kθ2 − θ1 kℓ2

(A)

≤

(C)

=

≤
(D)

≤

(B)

−1
kθ2 − θ1 k−1
ℓ2 kAk2→2 = kθ2 − θ1 kℓ2 max

n
X

n−1
v∈S+
i,j=1

max

n−1
v∈S+

max

X

vi vj

X

vi vj

Z1

i,j=1

n−1
v∈S+
i,j=1

0

Z1
0

vi Aij vj

h∇θ K (si − sj , θ1 + t (θ2 − θ1 )) ,

θ2 − θ1
idt
kθ2 − θ1 kℓ2

k∇θ K (si − sj , θ1 + t (θ2 − θ1 ))kℓ2 dt

CK,Θ kBk2→2 .

(A.4)

Notice that (A) and (B) follow from Lemma A.1. Identity (C) is implied by the fundamental
theorem of calculus for line integrals and inequality (D) is a consequence of (A3) in the Assumption
3.1. In the sequel we introduce an upper bound on kBk2→2 . Using Gersgorin’s circle theorem
(Theorem 6.1.1, [? ]), we get
kBk2→2 ≤ 1 + max

i=1,...,n

X
j6=i

(E)

1
1 + ksi −

sj kd+1
ℓ2

≤

Dmax
.
CK,Θ

(A.5)

Because of Assumption 2.1, we can apply Lemma D.1 of [4] to guarantee the existence of a bounded
positive scalar Dmax for which inequality (E) holds. Combining the inequalities (A.4) and (A.5)
concludes the proof of this part.
Now we turn to prove Eq. (A.3). Recall M and r from (A2) in Assumption 3.1 and set Dmin = M .
Applying inequality (3.1) terminates the proof by
X

2
kKn (θ2 ) − Kn (θ1 )k2ℓ2 =
K s − s′ , θ2 − K s − s′ , θ1
s,s′ ∈Dn

≥

=

X

max

h∈Dn (s,r)

s∈Dn
D2min n kθ2

|K (h, θ2 ) − K (h, θ1 )|2 ≥ M 2 n kθ2 − θ1 k2ℓ2

− θ1 k2ℓ2 .

For ease of reference, we present the following result as a standalone Proposition. Note that its
proof is akin to that of Proposition A.1 and so it will be skipped to avoid repetition.
Proposition A.2. Suppose that Dn admits Assumption 2.1. Moreover, Θ and K satisfy Assumption 3.2. Construct the matrix ∂Kn (θ) /∂θj := [∂K (s − s′ , θ) /∂θj ]s,s′ ∈Dn ,for θ ∈ Θ and
j = 1, . . . , m.
(a) There is a bounded strictly positive scalar Λ′max (depending on K, Θ, d and δ) such that
max max
n∈N θ∈Θ

∂
Kn (θ)
∂θj
30

2→2

≤ Λ′max ,

(b) There is D′max > 0 such that for any θ1 , θ2 ∈ Θ
∂
∂
Kn (θ2 ) −
Kn (θ1 )
∂θj
∂θj

2→2

≤ D′max kθ2 − θ1 kℓ2 .

(A.6)

The bounded positive scalars Dmax , Dmin and Λmax , which have been introduced in the Proposition
A.1, become frequently apparent in the subsequent results in this section. It is also proper to
remind the reader that Nǫ (A) stands for the ǫ−net of A with respect to the Euclidean distance.
Furthermore, the matrices Hθ1 ,θ2 and Mθ1 ,θ2 have been formerly defined in (6.1) for any pair of the
correlation function parameters θ1 , θ2 . The succeeding two Lemmas ( A.2 and A.3), which come
in handy in the proof of Theorem 3.1, establish a probabilistic upper bound on the maximum of a
quadratic Gaussian expression over a uncountable set Θ in terms of its largest value over one of its
finite subset. To avoid repetition, we omit the proof of the latter Lemma. Since it is akin to that
of the former one (with a slightly different algebra).
Lemma A.2. Let Z ∈ Rn be a standard Gaussian vector and suppose that Dn satisfies Assumption
2.1. Furthermore, assume that Θ and K admit Assumption 3.1. For any vanishing positive sequence
{rn }n∈N , any non-empty Θ̄ ⊂ Θ and each θ0 ∈ Θ,




√ 
(A.7)
LHS := lim Pr sup Z ⊤ Hθ0 ,θ Z − sup Z ⊤ Hθ0 ,θ Z  ≥ Crn n = 0,
n→∞

 θ∈Θ̄
θ∈Nr (Θ̄)
n

where C = 2Λmax (1 + Dmax ).

1/2

Proof. For simplicity, define Y = Kn (θ0 ) Z. Choose θ ∈ Θ̄ arbitrarily. Trivially, there is βθ ∈
Nrn Θ̄ such that kθ − βθ kℓ2 ≤ rn . Thus, with probability one
(
! )


Kn (θ)
Kn (βθ )
⊤
⊤
−
Y
LHS ≤ sup Z Hθ0 ,θ Z − Z Hθ0 ,βθ Z = sup Y
kKn (θ)kℓ2
kKn (βθ )kℓ2
θ∈Θ̄
θ∈Θ̄
≤ kY k2ℓ2 sup
θ∈Θ̄

Kn (βθ )
Kn (θ)
−
kKn (θ)kℓ2
kKn (βθ )kℓ2

.

(A.8)

2→2

Next we obtain an upper bound on the operator norm term in the inequality (A.8). Observe that
)
(
(A)
kKn (θ) − Kn (βθ )k2→2 kKn (βθ )k2→2 kKn (θ) − Kn (βθ )kℓ2
Kn (βθ )
Kn (θ)
−
+
≤
kKn (θ)kℓ2
kKn (βθ )kℓ2
kKn (θ)kℓ2
kKn (θ)kℓ2 kKn (βθ )kℓ2
2→2


(B)
kKn (θ) − Kn (βθ )k2→2 kKn (βθ )k2→2 kKn (θ) − Kn (βθ )kℓ2
√
+
≤
n
n
√ !
(C)
Λmax Dmax kθ − βθ kℓ2 n
Dmax kθ − βθ kℓ2
√
+
≤
n
n
≤

Dmax (1 + Λmax )
√
rn .
n

Here (A) follows from the triangle inequality alongside some straightforward algebra. Also (B) can
be deduced from the fact that Kn (·) is correlation matrix and (C) is implied by Proposition A.1.
31

(A.9)

Substituting Eq. (A.9) into Eq. (A.8) leads to
LHS ≤

Dmax (1 + Λmax )
√
rn kY k2ℓ2 .
n

(A.10)

We finish the proof by finding a tight probabilistic upper bound on kY k2ℓ2 . Obvious calculations
show that E kY k2ℓ2 = n. Hanson-Wright inequality (Theorem 1.1, [17]) ensures the existence of a
bounded universal constant C > 0 such that

 


√
√
kKn (θ0 )kℓ2
2
Pr kY kℓ2 − n ≥ C kKn (θ0 )kℓ2 ln n
≤ exp − ln n ∨ ln n
kKn (θ0 )k2→2
(
!)
√
(D)
1
n ln n
≤ exp − ln n ∨
= .
Λmax
n
Notice that (D) is implied by the part (a) of Proposition A.1 and the last inequality holds whenever
n ≥ n0 for some appropriately chosen n0 . Thus, with probability at least 1 − n−(1+ξ) , we have
p
p
kY k2ℓ2 ≤ n + C kKn (θ0 )kℓ2 (1 + ξ) ln n ≤ n + CΛmax (1 + ξ) n ln n ≤ 2n.
(A.11)

The last inequality in (A.11) holds if n ≥ n′0 for some n′0 . Assuming that n is at least n0 ∨ n′0 ,
combining the inequalities (A.10) and (A.11) concludes the proof.

Lemma A.3. Let Z ∈ Rn be a standard Gaussian vector. Suppose that Assumption 2.1 and
Assumption 3.1 hold for Dn , Θ and K. For any strictly positive vanishing sequence {rn }∞
n=1 , any
non-empty Θ̄ ⊂ Θ and arbitrary θ0 ∈ Θ,
(
)
Pr sup Z ⊤ (Mθ0 ,θ − Mθ0 ,βθ ) Z ≥ Crn
θ∈Θ̄

→ 0,

as n → ∞.



Here βθ represents the nearest element of Nrn Θ̄ to θ and C = 2Dmax 1 + 2Λ2max .

The next result not only appears in the proof of Theorem 3.2, but is also required for substantiating
some of the subsequent results in this section. Although it can be easily verified by Lemma A.1 of
[4], we provide its proof for the sake of completeness.
Lemma A.4. Under the same notation and conditions as Proposition A.1,
Kn (θ2 ) −

kKn (θ2 )kℓ2
Kn (θ1 )
kKn (θ1 )kℓ2

ℓ2

≥

kKn (θ2 ) − Kn (θ1 )kℓ2
.
Λmax

(A.12)

Proof. Let Ξ denote the left hand side of (A.12) and define λ := kKn (θ1 )k−1
ℓ2 kKn (θ2 )kℓ2 . Applying
Lemma A.1. of [4] along obvious calculations imply that


n 

2 
X
X
(λ − 1)2 +
K (si − sj , θ2 ) − λK (si − sj , θ1 )
Ξ2 =


i=1

≥

j6=i

n
P

n
X
j=1
i=1


2
K (si − sj , θ2 ) − K (si − sj , θ1 )
n
P

j=1

K 2 (si − sj , θ1 )

kKn (θ2 ) − Kn (θ1 )k2ℓ2
≥
.
n
P
maxi=1,...,n
K 2 (si − sj , θ1 )
j=1

32

Notice that the operator norm of any matrix is no smaller than the largest ℓ2 norm of its columns.
That is,
n
X
K 2 (si − sj , θ1 ) ≤ kKn (θ1 )k22→2 ≤ Λ2max .
max
i=1,...,n

j=1

The combination of these two inequalities ends the proof.
Now we state a lemma which plays a crucial role in the proof of Theorem 3.1.
Lemma A.5. Let Z ∈ Rn be a standard Gaussian vector and let C, ξ > 0. Suppose that Θ and
K satisfy Assumption 3.1. Select θ1 , θ2 ∈ Θ such that
r
ln n
kθ2 − θ1 kℓ2 ≥ Cmin
,
(A.13)
n
p
2
′
in which Cmin := 4D−1
min Λmax C (1 + ξ) (Recall Dmin and Λmax , from the Proposition A.1), for
some appropriately chosen universal constant C ′ > 0. There exists n0 = O (1) (depending on C, ξ,
K, Dn and Θ) such that for any n ≥ n0
)
(
r
ln n
⊤
≤ n−(1+ξ) .
(A.14)
p := Pr Z (Hθ2 ,θ2 − Hθ2 ,θ1 ) Z ≤ C
n
Refer to the identity (6.1) for the definition of Hθ2 ,θ1 .
Proof. We first obtain a lower bound on kKn (θ2 ) − Kn (θ1 )kℓ2 . Observe that,
(A)
(B)
p
√
kKn (θ2 ) − Kn (θ1 )kℓ2 ≥ Dmin n kθ2 − θ1 kℓ2 ≥ 4Λ2max C ′ (1 + ξ) ln n,

(A.15)

where (A) follows from inequality (A.3) in Proposition A.1 and (B) is implied by condition (A.13).
For brevity, let A := Hθ2 ,θ2 − Hθ2 ,θ1 . It is trivial from the identity (6.1) that
!
K
(θ
)
K
(θ
)
n
2
n
1
A = Kn1/2 (θ2 )
−
Kn1/2 (θ2 ) .
(A.16)
kKn (θ2 )kℓ2
kKn (θ1 )kℓ2
Obtaining a closed form for the trace and Frobenius norm of A, which express the expected value
and variance of Z ⊤ AZ, plays a significant role in the proof. Straightforward algebra leads to
kKn (θ2 )kℓ2
hKn (θ1 ) , Kn (θ2 )i
=
tr (A) = kKn (θ2 )kℓ2 −
kKn (θ1 )kℓ2
2

Kn (θ1 )
Kn (θ2 )
−
kKn (θ2 )kℓ2
kKn (θ1 )kℓ2

2

. (A.17)
ℓ2

Applying Lemma A.4 to the identity (A.17) shows that
tr (A)

≥
(C)

≥

1
2 kKn (θ2 )kℓ2
8C ′ Λmax



kKn (θ2 ) − Kn (θ1 )kℓ2
Λmax

(1 + ξ) ln n
√
.
n

2

≥

1

√

2Λmax n



kKn (θ2 ) − Kn (θ1 )kℓ2
Λmax

2

(A.18)

33

Note that (C) follows from
√ the condition (A.15), and (A.18) assures the existence
 ⊤ of some n2 ∈ N
−1
for which tr (A) ≥ 2C n ln n, whenever n ≥ n2 . Hence, p ≤ q := Pr Z AZ ≤ 1/2 tr (A) .
We aim to show that q is smaller than n−(1+ξ) . Based upon Hanson-Wright inequality there is a
bounded universal constant C ′ > 0 for which
n

o
p
p
Pr Z ⊤ AZ ≤ tr (A) − C ′ (1 + ξ) ln n kAkℓ2 ∨ kAk2→2 (1 + ξ) ln n
≤ n−(1+ξ) .
(A.19)

Thus it suffices to show that for some n3 ∈ N


p
p
tr (A) ≥ 2 C ′ (1 + ξ) ln n kAkℓ2 ∨ kAk2→2 (1 + ξ) ln n ,

Claim 1. Define a0 := Λmin Dmin Λ3max Dmax (1 + Λmax )

−1

∀ n ≥ n3 .

(A.20)

√
. Then, kAkℓ2 ≥ a0 n kAk2→2 .

Proof of Claim 1. In order to substantiate the claim, we require a lower bound on kAkℓ2 and an
upper bound on kAk2→2 . Applying the part (a) of Proposition A.1 (alongside the Cauchy-Schwartz
inequality) on the identity (A.16) shows that
Λ−1
max kAk2→2 ≤ Lup :=

Kn (θ1 )
Kn (θ2 )
−
kKn (θ2 )kℓ2
kKn (θ1 )kℓ2

Λ−1
min kAkℓ2 ≥ Llw :=

Kn (θ2 )
Kn (θ1 )
−
kKn (θ2 )kℓ2
kKn (θ1 )kℓ2

,
2→2

.

(A.21)

ℓ2

We achieve an upper bound on Lup by taking advantage of the same technique as (A.9). Furthermore, a tight lower bound on Llw can be obtained from Lemma A.4. The algebraic details are
skipped to avoid redundancy.
Lup ≤ Dmax (1 + Λmax )

kθ2 − θ1 kℓ2
√
,
n

Llw ≥

Dmin kθ2 − θ1 kℓ2
.
Λ2max

(A.22)

The combination of (A.21) and (A.22) finishes the proof.
With Claim 1 in hand, proving (A.20) is equivalent to showing that
p
tr (A) ≥ 2C ′ (1 + ξ) ln n kAkℓ2 ,

(A.23)

for large enough n. Obviously, kAkℓ2 ≤ Λmax Llw . Thus, identity (A.17) implies that
tr (A)
kAkℓ2

≥

kKn (θ2 )kℓ2
2Λmax Llw

= b0 :=

Kn (θ1 )
Kn (θ2 )
−
kKn (θ2 )kℓ2
kKn (θ1 )kℓ2

2

ℓ2

kKn (θ2 )kℓ2
1
Kn (θ1 )
Kn (θ2 ) −
2Λmax
kKn (θ1 )kℓ2

(A.24)
ℓ2

In the sequel, notice that
(D)

b0 ≥

kKn (θ2 ) − Kn (θ1 )kℓ2 (E) p ′
≥ 2 C (1 + ξ) ln n,
2Λ2max

in which (D) and (E) respectively follow from Lemma A.4 and condition (A.15). Substituting the
last inequality into (A.24) verifies (A.23) and concludes the proof.
34

The next proposition rigorously expresses the uniform concentration of the Euclidean squared norm
of Gaussian vectors with the covariance matrix Kn (θ) , θ ∈ Θ around their mean. It is worthwhile
to mention that such inequality is crucial for proving the Theorem 3.2.
Proposition A.3. Let Θ ⊂ Rm be a bounded set. Consider the class of n by n matrices
{Πn (θ)}θ∈Θ parametrized by θ ∈ Θ. Suppose that the following conditions hold
(a) The operator norm of Πn (θ) is uniformly bounded in Θ. Namely,
M := sup sup kΠn (θ)k2→2 < ∞.
n θ∈Θ


(b) The mapping θ, k·kℓ2 →
7 (Πn (θ) , k·k2→2 ) is Lipschitz. Namely, there is C > 0 for which
kΠn (θ2 ) − Πn (θ1 )k2→2 ≤ C kθ2 − θ1 kℓ2 ,

∀ θ1 , θ2 ∈ Θ.

(A.25)

(c)
kΠn (θ)k2→2
=o
kΠn (θ)kℓ2



1
√
ln n



,

∀ θ ∈ Θ.

Then, there is a constant C ′ > 0 such that


√
⊤
′
lim Pr sup Z Πn (θ) Z − tr {Πn (θ)} ≥ C n ln n = 0.
n→∞

(A.26)

θ∈Θ

p
Proof. Let rn = C −1 ln n/n in which C has been defined in (A.25) and let Nrn (Θ) denote the
rn -covering set of Θ. As before, for any θ let βθ represents the closest element of Nrn (Θ) to θ.
Observe that,
RHS := Z ⊤ Πn (θ) Z − tr {Πn (θ)} − Z ⊤ Πn (βθ ) Z + tr {Πn (βθ )}
= hΠn (θ) − Πn (βθ ) , ZZ ⊤ + In i ≤ kΠn (θ) − Πn (βθ )k2→2 ZZ ⊤ + In
S1
r
(A)

ln n
=
n + kZkℓ2 ,
≤ Crn ZZ ⊤ + In
≤ C kθ − βθ kℓ2 ZZ ⊤ + In
n
S1
S1

in which k·kS1 stands for the nuclear norm (absolute sum of eigenvalues). Note that the obtained
upper bound does not depend on θ (uniform upper
bound). Moreover, based upon Hanson-Wright

inequality there is c > 0 for which n + kZkℓ2 ≤ 3n with probability at least 1 − exp (−cn). Thus,
√
RHS ≥ 3 n ln n with probability at most exp (−cn). Hence, as n → ∞ we get
!
√
⊤
⊤
Pr sup Z Πn (θ) Z − tr {Πn (θ)} ≥ sup
Z Πn (θ) Z − tr {Πn (θ)} + 3 n ln n → 0.
θ∈Θ

θ∈Nrn (Θ)

(A.27)
− tr {Πn (θ)} . Applying
In the sequel we find a tight upper bound on
condition (c) on Hanson-wright inequality and using union bound leads to
!
√
Z ⊤ Πn (θ) Z − tr {Πn (θ)} ≥ C0 sup kΠn (θ)kℓ2 ln n ≤ |Nrn (Θ)| n−m ,
Pr
sup
supθ∈Nrn (Θ) Z ⊤ Πn (θ) Z

θ∈Θ

θ∈Nrn (Θ)

35

√
for some constant C0 > 0 depending on m. Notice that supθ∈Θ kΠn (θ)kℓ2 ≤ M n according to
condition (a). Moreover, as we argued in (6.14), |Nrn (Θ)| = o (nm ). Thus,
!
√
lim Pr
sup
Z ⊤ Πn (θ) Z − tr {Πn (θ)} ≥ C0 M n ln n = 0.
n→∞

θ∈Nrn (Θ)

Replacing the last inequality into (A.27) concludes the proof.
Lemma A.6. Under the same notation and assumptions as Propositions A.1 and A.2 and for
any unit norm vector λ, the class of matrices defined by
)
(
K
(θ)
n
Tn (θ) := lim γ −1 Kn (θ + γλ) − hKn (θ + γλ) , Kn (θ)i
γց0
kKn (θ)k2ℓ2
satisfy the conditions (a) and (b) in Proposition A.3.
Proof. The following inequality uniformly holds for any θ and n.
(

kTn (θ)k2→2

(A)

≤

(B)

≤

lim γ

−1

lim γ

−1

lim γ

−1

γց0

γց0

(C)

≤

γց0

kKn (θ + γλ) − Kn (θ)k2→2 + kKn (θ)k2→2

(

kKn (θ + γλ) − Kn (θ)k2→2 +



Λmax √
nDmax γ
Dmax γ + √
n



hKn (θ + γλ) , Kn (θ)i
1−
kKn (θ)k2ℓ2
)

kKn (θ)k2→2
kKn (θ + γλ) − Kn (θ)kℓ2
kKn (θ)kℓ2

= Dmax (1 + Λmax ) < ∞.

In above, (A) and (B) are simple applications of triangle and Cauchy-Schwartz inequalities, respectively. Furthermore, (C) follows from Propositions A.1. Namely, supn∈N supθ∈Θ kTn (θ)k2→2 < ∞
which is same as condition (a).
The next phase of the proof is devoted to justify condition (b). We have introduced an equivalent
representation for Tn (θ) in (6.20). Define
m
X
∂
λj
Kn (θ) , Kn (θ)i
an (θ) := kKn (θ)k−2
h
ℓ2
∂θj
j=1

It follows from Propositions A.1 and A.2 that both matrices Kn (θ) and

m
P

λj ∂Kn (θ) /∂θj admit

j=1

conditions (a) and (b) in Proposition A.3. Thus, according to Lemma A.7, it is sufficient to prove
that an (θ) is a uniformly bounded and Lipschitz function in Θ. Notice that,

sup |an (θ)| ≤ sup
θ∈Θ

j=1

λj ∂θ∂ j Kn (θ)
ℓ2

≤ sup

r

m

∂
1X
|λj |
Kn (θ)
n
∂θj

kKn (θ)kℓ2
θ∈Θ
j=1
√
′
Λmax n kλkℓ1
√
√
≤ Λ′max m, ∀ n ∈ N.
n

θ∈Θ

≤

m
P

36

ℓ2

(A.28)

)

In the sequel, we justify the Lipschitz property. Taking advantage of Proposition A.2 and inequality
(A.28), one show that for any θ ∈ Θ,
m
X

λj

∂
Kn (θ)
∂θj

λj



j=1

m
X

√
≤ Λ′max mn,
ℓ2

j=1

∂
∂
Kn (θ2 ) −
Kn (θ1 )
∂θj
∂θj



√
≤ D′max mn kθ2 − θ1 kℓ2 .

(A.29)

ℓ2

Thus
|an (θ2 ) − an (θ1 )|

≤



m
X
∂
∂
Kn (θ1 )
i
λj
h
Kn (θ2 ) −
Kn (θ1 ) ,
∂θj
∂θj
kKn (θ1 )k2ℓ

+

m
X
Kn (θ1 )
∂
Kn (θ2 )
−
i
λj
h
Kn (θ2 ) ,
2
∂θj
kKn (θ2 )kℓ2
kKn (θ1 )k2ℓ2
j=1

(D)

≤
+

j=1

m
P

≤

(F )

≤

λj

j=1

m
X
j=1

(E)

2

√

λj



∂
∂
∂θj Kn (θ2 ) − ∂θj Kn (θ1 )

kKn (θ1 )kℓ2

∂
Kn (θ2 )
∂θj
ℓ2



ℓ2

Kn (θ1 )
Kn (θ2 )
2 −
kKn (θ2 )kℓ2
kKn (θ1 )k2ℓ2

ℓ2

mnD′max kθ2 − θ1 kℓ2 √
√
Kn (θ1 )
Kn (θ2 )
√
+ mnΛ′max n
−
2
n
kKn (θ2 )kℓ2
kKn (θ1 )k2ℓ2

2→2

√
√
D′max m kθ2 − θ1 kℓ2 + 3Λ′max mΛmax Dmax kθ2 − θ1 kℓ2 = C0 kθ2 − θ1 kℓ2 ,

where (D) is implied by triangle inequality, (E) follows from (A.29). Lastly, (F ) can be proved
using akin tricks as (A.9). Thus, an (·) is a Lipschitz function with some constant C0 .
The following technical result, which comes in handy in the proof of Lemma A.6, is an easy
consequence of triangle inequality. The proof is skipped for brevity.
Lemma A.7. Suppose that a1 , a2 : Θ →
7 R are two uniformly bounded and Lipschitz functions
with constant L1a and L2a , respectively. Consider matrices A1 , A2 : Θ 7→ Rn×n parametrized by
θ ∈ Θ such that
kAi (θ2 ) − Ai (θ1 )k2→2 ≤ LiA kθ2 − θ1 kℓ2 ,

∀ θ1 , θ2 ∈ Θ, i = 1, 2

Let A (θ) := a1 (θ) A1 (θ) + a2 (θ) A2 (θ) for any θ. There is a scalar L, which depends on La , Lb ,
LA and LB , for which the following inequality holds.
kA (θ2 ) − A (θ1 )k2→2 ≤ L kθ2 − θ1 kℓ2 ,
37

∀ θ1 , θ2 ∈ Θ.

The following result gives an upper bound on the Kullback-Leibler divergence of two zero mean
multivariate Gaussian distributions respectively associated with the two covariance matrices
Kn (θi ) , i = 1, 2. Such upper bound is extremely useful for establishing Theorem 3.3.
Proposition A.4. Choose θ1 , θ2 ∈ Θ in such a way that kθ2 − θ1 kℓ2 ≤ Λmin / (2Dmax ). Let
Pi , i = 1, 2, denotes the associated probability distribution to a zero mean Gaussian vector with
the covariance matrix Kn (θi ) ∈ Rn×n , i = 1, 2. Then,
2

Dmax
kθ2 − θ1 kℓ2 .
D (P1 k P2 ) ≤ 2n
Λmin
Proof. For any symmetric matrix A ∈ Rn×n , let λi (A) , i = 1, . . . , n, denotes its ith eigenvalue in
decreasing order. The von Neumann’s trace inequality [14] yields


det Kn (θ2 )
−1
D (P1 k P2 ) = hKn (θ2 ) , Kn (θ1 )i − n + ln
det Kn (θ1 )


n
X λj (Kn (θ1 ))
λj (Kn (θ1 ))
≤ Q :=
− 1 − ln
.
λj (Kn (θ2 ))
λj (Kn (θ1 ))
j=1

We finish the proof by acquiring a proper upper bound on Q. Define f : (0, ∞) 7→ R by f (x) =
|x − 1 − ln x|. Applying the second order Taylor’s expansion around x = 1 shows that f (x) ≤
2 (x − 1)2 for |x − 1| ≤ 1/2.
Claim 1. The succeeding inequality holds for any j = 1, . . . , n.
λj (Kn (θ1 ))
1
Dmax
−1 ≤
kθ2 − θ1 kℓ2 ≤ .
λj (Kn (θ2 ))
Λmin
2
Claim 1 provides the key tool to control Q from above.
2
2


n 
n 
n
X
X
X
λj {Kn (θ1 )}
λj {Kn (θ1 )}
Dmax
−1 ≤2
kθ2 − θ1 kℓ2
≤2
f
Q =
λj {Kn (θ2 )}
λj {Kn (θ2 )}
Λmin
j=1
j=1
j=1

2
Dmax
= 2n
kθ2 − θ1 kℓ2 .
Λmin
In conclusion, we substantiate Claim 1. The first inequality can be established using Proposition
A.1 and the second one is obvious.
Dmax kθ2 − θ1 kℓ2
kKn (θ2 ) − Kn (θ1 )k2→2
λj {Kn (θ1 )} − λj {Kn (θ2 )}
λj {Kn (θ1 )}
≤
−1 =
≤
.
λj {Kn (θ2 )}
λj {Kn (θ2 )}
λn {Kn (θ2 )}
Λmin
Now we demonstrate the asymptotic normality of the normalized quadratic Gaussian forms. We
exploit this fact in the proof of Theorem 3.4.
Lemma A.8. For n ∈ N, let Zn ∈ Rn be a standard Gaussian vector and let An ∈ Rn×n . Then,
(
)
Zn⊤ An Zn − tr (An ) d
Ψn :=
→ N (0, 2) ,
kAn kℓ2
provided that lim kAn k−1
ℓ2 kAn k2→2 = 0.
n→∞

38

Proof. Let Ψ∞ be a zero mean Gaussian random variable with variance 2. So, ln E exp (tΨ∞ ) = t2
for any t ∈ R. The basic properties of the quadratic forms of Gaussian vectors yields
!
t tr (An )
An
1
−
ln E exp (tΨn ) = − ln det In − 2t
2
kAn kℓ2
kAn kℓ2
!
)
(
n
2tλj (An )
1X
2tλj (An )
+
= −
ln 1 −
2
kAn kℓ2
kAn kℓ2
j=1


!2
!2  
n

X
tλj (An )
tλj (An ) 
(A)

=
+o
→ t2 , as n → ∞.

 kAn kℓ2
kAn kℓ2
j=1

Here (A) follows from expanding ln (1 − x) around 1 for infinitesimal x (since λj (An ) / kAn kℓ2
vanishes as n → ∞). Consequently, Ψn converges in distribution to Ψ∞ by the continuity theorem
of moment generating functions.

The last result of this section studies the shrinkage behaviour of the partial derivatives of Matern
covariance function with respect to its fractal index. In turns out to be useful for corroborating
the part (a) of Remark 3.1.
Lemma A.9. Let Kν : Rd 7→ R be a geometric anisotropic (Recall from Definition 3.1) Matern
correlation function given by
√

21−ν
r ⊤ Ar ,
Kν
Kν (r) =
Γ (ν)

in which Kν stands for the modified Bessel function of the second kind and A satisfies the condition
3.3. Then for any β ∈ N and m ∈ N, there is a bounded constant Cβ,A such that
Cβ,A
∂m
Kν (r) ≤
,
∂ν m
1 + krk2β
ℓ2

∀ r = (r1 , . . . , rd ) ∈ Rd .

(A.30)

Proof. Choose β, m ∈ N arbitrarily. The inequality (A.30) obviously holds if we can show that
!
d
X
∂m
2β
ri
sup 1 +
Kν (r) < ∞.
(A.31)
∂ν m
r∈Rd
i=1

Let F stands for the Fourier operator. For proving (A.31) it suffices to show that
)
!
(
d
m
X
∂
Kν (r) ∈ L1 (Rd ),
ri2β
F
1+
∂ν m

(A.32)

i=1

where L1 (Rd ) stands for the class of absolutely integrable d−variate functions. For doing so we
need a closed form expression for the Matern spectral density K̂ν , which can be done using the
integration by substitution technique.
K̂ν (ω) =

−(ν+d/2)
π −d/2 
1 + ω ⊤ A−1 ω
|det (A)|
39

The condition 3.3 says that K̂ν decays polynomially in terms of kωkℓ2 . We now evaluate the Fourier
m
ν
transform of ∂∂νK
m . Observe that
Z
Z


∂ m Kν
∂m
m
⊤ −1
−jhr,ωi
m
−jhr,ωi
1
+
ω
A
ω
dω.
e
K̂
(ω)
ln
e
K̂
(ω)
dω
=
(−1)
=
ν
ν
∂ν m
∂ν m Rd
Rd


m
ν
= Ĝν (ω) := (−1)m K̂ν (ω) lnm 1 + ω ⊤ A−1 ω . Thus,
In other words, F ∂∂νK
m
F

(

1+

d
X
i=1

ri2β

!

)
d
X
∂m
∂ 2β Ĝν
β
.
K
(r)
=
Ĝ
+
(−1)
ν
ν
∂ν m
∂ωi2β
i=1

The absolutely integrability of Ĝν is clear. The reader can verify the same property for the functions
∂ 2β Ĝν
2β , i = 1, . . . , d using simple differentiation. So (A.32) holds by the triangle inequality.
∂ωi

References
References
[1] M. Abramowitz, and I.A. Stegun. “Handbook of mathematical functions”, Applied Mathematics Series 55 (1966): 62.
[2] M. Anitescu, J. Chen, and M.L. Stein, “An inversion-free estimating equation approach for
Gaussian process models”, submitted for publication (2014).
[3] S.B. Aruoba, and J. Fernndez-Villaverde, “A comparison of programming languages in economics”, No. w20263. National Bureau of Economic Research, 2014.
[4] F. Bachoc, “Asymptotic analysis of the role of spatial sampling for covariance parameter
estimation of Gaussian processes”, Journal of Multivariate Analysis 125 (2014): 1-35.
[5] J. Bernardo, J. Berger, A. Dawid, and A. Smith. “Regression and classification using Gaussian
process priors”, Bayesian statistics 6 (1998): 475.
[6] R.H. Byrd, P. Lu, J. Nocedal, and C. Zhu, “A limited memory algorithm for bound constrained
optimization”, SIAM Journal on Scientific Computing 16, no. 5 (1995): 1190-1208.
[7] J. Chen, “On the use of discrete Laplace operator for preconditioning kernel matrices”, SIAM
Journal on Scientific Computing 35, no. 2 (2013): A577-A602.
[8] W. Chu, and Z. Ghahramani, “Gaussian processes for ordinal regression”, Journal of Machine
Learning Research, pp. 1019-1041. 2005.
[9] N. Cressie, “Statistics for spatial data”, John Wiley & Sons, 2015.
[10] J. Du, H. Zhang and V.S. Mandrekar, “Fixed-domain asymptotic properties of tapered maximum likelihood estimators”, The Annals of Statistics 37, no. 6A (2009): 3330-3361.
[11] A.E. Gelfand, P. Diggle, P. Guttorp, and M. Fuentes, “Handbook of spatial statistics”, CRC
press, 2010.
40

[12] C.G. Kaufman, M.J. Schervish, and D.W. Nychka, “Covariance tapering for likelihood-based
estimation in large spatial data sets”, Journal of the American Statistical Association 103, no.
484 (2008): 1545-1555.
[13] K.V. Mardia, and R.J. Marshall, “Maximum likelihood estimation of models for residual covariance in spatial regression”, Biometrika 71, no. 1 (1984): 135-146.
[14] L. Mirsky, “A trace inequality of John von Neumann”, Monatshefte fr Mathematik 79, no. 4
(1975): 303-306.
[15] D.P. O’Leary, “The block conjugate gradient algorithm and related methods”, Linear algebra
and its applications 29 (1980): 293-322.
[16] C.K. Williams, and C.E. Rasmussen, “Gaussian Processes for Machine Learning”, The MIT
Press, 2(3), p.4.
[17] M. Rudelson, and R. Vershynin, “Hanson-Wright inequality and sub-Gaussian concentration”,
arXiv preprint arXiv:1306.2872 (2013).
[18] J. Sacks, W.J. Welch, T.J. Mitchell, and H.P. Wynn, “Design and analysis of computer experiments”, Statistical science (1989): 409-423.
[19] M.L. Stein, “Interpolation of spatial data: some theory for kriging”, Springer Science & Business Media, 2012.
[20] M.L. Stein, J. Chen, and M. Anitescu, “Difference filter preconditioning for large covariance
matrices”, SIAM Journal on Matrix Analysis and Applications 33, no. 1 (2012): 52-72.
[21] A.B. Tsybakov, “Introduction to Nonparametric Estimation”, Springer Series in Statistics,
2009.
[22] D. Wang, and W.L. Loh, “On fixed-domain asymptotics and covariance tapering in Gaussian
random field models”, Electronic Journal of Statistics 5 (2011): 238-269.
[23] H. Wendland, “Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree”, Advances in computational Mathematics 4, no. 1 (1995): 389-396.
[24] W.Y. Wu, C.Y. Lim, and Y. Xiao, “Tail estimation of the spectral density for a stationary
Gaussian random field”, Journal of Multivariate Analysis 116 (2013): 74-91.
[25] Z. Ying, “Asymptotic properties of a maximum likelihood estimator with data from a Gaussian
process”, Journal of Multivariate Analysis 36, no. 2 (1991): 280-296.
[26] H. Zhang, “Inconsistent estimation and asymptotically equal interpolations in model-based
geostatistics”, Journal of the American Statistical Association 99, no. 465 (2004): 250-261.

41

