Asymptotic Behavior of the Likelihood Function of Covariance Matrices of Spatial Gaussian Processes

Abstract

The covariance structure of spatial Gaussian predictors (aka Kriging predictors) is generally modeled by parameterized covariance functions; the associated hyperparameters in turn are estimated via the method of maximum likelihood. In this work, the asymptotic behavior of the maximum likelihood of spatial Gaussian predictor models as a function of its hyperparameters is investigated theoretically. Asymptotic sandwich bounds for the maximum likelihood function in terms of the condition number of the associated covariance matrix are established. As a consequence, the main result is obtained: optimally trained nondegenerate spatial Gaussian processes cannot feature arbitrary ill-conditioned correlation matrices. The implication of this theorem on Kriging hyperparameter optimization is exposed. A nonartificial example is presented, where maximum likelihood-based Kriging model training is necessarily bound to fail.

1. Introduction

Spatial Gaussian processing, also known as best linear unbiased prediction, refers to a statistical data interpolation method, which is nowadays applied in a wide range of scientific fields, including computer experiments in modern engineering context; see, for example, [1–5]. As a powerful tool for geostatistics, it has been pioneered by Krige in 1951 [6], and to pay tribute to his achievements, the method is also termed Kriging; see [7, 8] for geostatistical background.

In practical applications, the data's covariance structure is modeled through covariance functions depending on the so-called hyperparameters. These, in turn, are estimated by optimizing the corresponding maximum likelihood function. It has been demonstrated by many authors that the accuracy of Kriging predictors relies both heavily on hyperparameter-based model training and, from the numerical point of view, on the condition number of the associated Kriging correlation matrix. In this regard, we relate to the following, nonexhaustive selection of papers: Warnes and Ripley [9] and Mardia and Watkins [10] present numerical examples of difficult-to-optimize covariance model functions. Ababou et al. [11] show that likelihood-optimized hyperparameters may correspond to ill-conditioned correlation matrices. Diamond and Armstrong [12] prove error estimates under perturbation of covariance models, demonstrating a strong dependence on the correlation matrix' condition number. In the same setting, Posa [13] investigates numerically the behavior of this precise condition number for different covariance models and varying hyperparameters. An extensive experimental study of the condition number as a function of all parameters in the Kriging exercise is provided by Davis and Morris [14]. Schöttle and Werner [15] propose Kriging model training under suitable conditioning constraints. Related is the work of Ying [16] and Zhang and Zimmerman [17], who prove asymptotic results on limiting distributions of maximum likelihood estimators when the number of sample points approaches infinity. Radial basis function interpolant limits are investigated in [18]. Modern textbooks covering recent results are [19, 20].

In this paper, the connection between hyper-parameter optimization and the condition number of the correlation matrix is investigated from a theoretical point of view. The setting is as follows. All sample data is considered as fixed. An arbitrary feasible covariance model function is chosen for good, so that only the covariance models' hyperparameters are allowed to vary in order to adjust the model likelihood. This is exactly the situation as it occurs in the context of computer experiments, where, based on a fixed set of sample data, predictor models have to be trained numerically. We prove that, under weak conditions, the limit values of the quantities in the model training exercise exist. Subsequently, by establishing asymptotic sandwich bounds on the model likelihood based on the condition number of the associated correlation matrix, it is shown that ill-conditioning eventually also decreases the model likelihood. This result implies a strategy for choosing good starting solutions for hyperparameter-based model training. We emphasize that all covariance models applied in the papers briefly reviewed above subordinate to the theoretical setting of this work.

The paper is organized as follows. In the next section, a short review of the basic theory behind Kriging is given. The main theorem is stated and proved in Section 3. In Section 4, an example of a Kriging data set is presented, which illustrates the limitations of classical model training.

2. Kriging in a Nutshell

Kriging is a statistical approach for estimating an unknown (scalar) function𝑦∶ℝ^𝑑⊇𝑈⟶ℝ,𝑥⟼𝑦(𝑥), (2.1) based on a finite data set of sample locations 𝑥¹,…,𝑥^𝑛∈𝑈⊂ℝ^𝑑 with corresponding responses 𝑦₁∶=𝑦(𝑥¹),…,𝑦_𝑛∶=𝑦(𝑥^𝑛)∈ℝ obtained from measurements or numerical computations. The collection of responses is denoted by𝑌^𝑇=𝑦₁,…,𝑦_𝑛∈ℝ^𝑛. (2.2) The function 𝑦∶𝑈→ℝ to be estimated is assumed to be the realization of an underlying random process given by a regression model and a random error function 𝜖(𝑥) with zero mean. More precisely𝑦(𝑥)=𝑓(𝑥)𝛽+𝜖(𝑥)=𝑓₀(𝑥),…,𝑓_𝑝(𝑥)𝛽₀⋮𝛽_𝑝+𝜖(𝑥), (2.3) where the components of the row vector function 𝑓∶ℝ^𝑑→ℝ^𝑝+1 are the basis functions of the regression model and 𝛽=(𝛽₀,…,𝛽_𝑝) is the corresponding vector of regression coefficients. By assumption,𝐸[𝜖(𝑥)]=0∀𝑥. (2.4) The component functions of 𝑓 can be chosen arbitrarily, yet they should form a function basis suitable to the specific application. The most common choices for practical applications are (1)constant regression (ordinary Kriging): 𝑝=0, 𝑓∶ℝ^𝑑→ℝ,𝑥↦1, 𝑓(𝑥)𝛽=𝛽∈ℝ,(2)linear regression (universal Kriging): 𝑝=𝑑, 𝑓∶ℝ^𝑑→ℝ^𝑑+1,𝑥↦(1,𝑥₁,…,𝑥_𝑑), 𝑓(𝑥)𝛽=𝛽₀+𝛽₁𝑥₁+⋯+𝛽_𝑑𝑥_𝑑∈ℝ,

and higher-order polynomials.

Introducing the regression design matrix𝑓𝑥𝐹∶=⋮¹𝑓𝑥^𝑛=𝑓₀𝑥¹𝑓₁𝑥¹⋯𝑓_𝑝𝑥¹𝑓⋮⋮⋯⋮₀𝑥^𝑛𝑓₁𝑥^𝑛⋯𝑓_𝑝𝑥^𝑛∈ℝ^{𝑛×(𝑝+1)}, (2.5) the vector of errors at the sampled sites can be written as𝜖Σ∶=⋮¹𝜖_𝑛=𝜖𝑥¹⋮𝜖𝑥^𝑛=𝑦₁−𝑓𝑥¹𝛽⋮𝑦_𝑛−𝑓𝑥^𝑛𝛽=𝑌−𝐹𝛽. (2.6) Note that the first column of 𝐹 equals 𝟏∈ℝ^𝑛 for all polynomial regression models.

The Kriging predictor 𝑦 estimates 𝑦 at an untried site 𝑥 as a linear combination of the sampled data𝑦(𝑥)=⟨𝜔(𝑥),𝑌⟩=^𝑛_𝑖=1𝜔_𝑖(𝑥)𝑦_𝑖. (2.7) For each 𝑥∈ℝ^𝑑, the unique vector of weights 𝜔(𝑥)=(𝜔₁(𝑥),…,𝜔_𝑛(𝑥)) that leads to an unbiased prediction minimizing the mean squared error is given by the solution of the Kriging equation system𝐹𝐶𝐹_𝑇0𝜔(𝑥)12𝜇(𝑥)=𝑐(𝑥)𝑓(𝑥)^𝑇∈ℝ^{𝑛+(𝑝+1)}. (2.8) Here,𝐶∶=Cov𝜖𝑥^𝑖,𝜖𝑥^𝑗_{𝑖,𝑗≤𝑛}∈ℝ^𝑛×𝑛,𝑐(𝑥)∶=Cov(𝜖(𝑥^𝑖),𝜖(𝑥))_𝑖≤𝑛∈ℝ^𝑛 (2.9) denote the covariance matrix and the covariance vector, respectively, and the entries of the vector 𝜇=(𝜇₀,…,𝜇_𝑝) are Lagrange multipliers. Solving (2.8) by Schur matrix complement inversion and substituting in (2.7) leads to the Kriging predictor formula𝑦(𝑥)=𝑓(𝑥)𝛽+𝑐^𝑇(𝑥)𝐶⁻¹𝑌−𝐹𝛽, (2.10) where 𝛽=(𝐹^𝑇𝑅⁻¹𝐹)⁻¹𝐹^𝑇𝑅⁻¹𝑌 is the generalized least squares solution to the regression problem 𝐹𝛽≃𝑌. For details, see, for example, [20, 21].

For setting up Kriging predictors, it is therefore mandatory to estimate the covariances based on the sampled data set. The two most popular approaches to tackle this problem are variogram fitting (the geostatistical literature, see [7]) and application of spatial correlation functions (computer experiments, see [4, 20]). The latter ones are usually of the formCov𝜖𝑥^𝑖,𝜖𝑥^𝑗=𝜎²𝑟𝜃,𝑥^𝑖,𝑥^𝑗=𝜎^2𝑑_𝑘=1scf_𝑘𝜃,𝑥^𝑖,𝑥^𝑗. (2.11) Here 𝜃=(𝜃₁,…,𝜃_𝑑)∈ℝ^𝑑 is a vector of distance weights, which models the influence of the coordinate-wise spatial correlation on the prediction. The correlation matrix is defined by𝑅(𝜃)=𝑟𝜃,𝑥^𝑖,𝑥^𝑗_𝑖,𝑗∈ℝ^𝑛×𝑛. (2.12) In order to avoid ambiguity due to different parameterizations of the correlation models, we fix for the rest of the paper the following.

Convention 1
Large distance weight values correspond to weak spatial correlation, and small distance weight values correspond to strong spatial correlation. More precisely, we assume that feasible spatial correlation functions are always parameterized such that𝑟𝜃,𝑝,𝑞⟶1,for‖𝜃‖⟶0,0,for‖𝜃‖⟶∞, (2.13) at distinct locations 𝑝≠𝑞.

All correlation models applied in all the papers briefly reviewed in the introduction can be parameterized accordingly. A collection of spatial correlation functions is given in several publications on Cokriging/Kriging, including [21, Table 2.1]. For example, the Gaussian correlation function parameterized with respect to the convention above is given byscf_𝑘𝜃,𝑥^𝑖,𝑥^𝑗=exp−𝜃_𝑘𝑥^𝑖_𝑘−𝑥^𝑗_𝑘²,𝑟𝜃,𝑥^𝑖,𝑥^𝑘=exp−_𝑘𝜃_𝑘𝑥^𝑖_𝑘−𝑥^𝑗_𝑘². (2.14) The results and proofs presented below hold true without change for every admissible spatial correlation model, assuming that the sample errors are normally distributed and that the process variance is stationary; that is, 𝜎 is independent of the locations 𝑥^𝑖, 𝑥^𝑗. In this setting, hyper-parameter training for Kriging models consists of optimizing the corresponding maximum likelihood function =1MaxL𝜎,𝛽,𝜃2𝜋det𝜎²𝑅(𝜃)−1exp2𝜎²𝑌−𝐹𝛽.^𝑇𝑅(𝜃)⁻¹𝑌−𝐹𝛽.. (2.15) For 𝜃 fixed, optima for 𝜎=𝜎(𝜃) and 𝛽=𝛽(𝜃) can be derived analytically, see [20], so that hyper-parameter training for Kriging models is reduced to the following optimization problem:min_{{𝜃∈ℝ𝑑,𝜃𝑗>0}}ML(𝜃)∶=𝜎²(𝜃)^𝑛det(𝑅(𝜃)), (2.16) where the dependency on 𝜃 is as follows:𝜎²(𝜃)=1𝑛Σ^𝑇(𝜃)𝑅(𝜃)⁻¹Σ(𝜃),Σ(𝜃)[𝑌](𝜃)𝛽(𝜃)=𝐹=𝑌−𝐸=𝑌−𝐹𝛽,_𝑇𝑅⁻¹(𝜃)𝐹⁻¹𝐹^𝑇𝑅⁻¹(𝜃)𝑌. (2.17)(2.18)(2.19) Because the logarithm is monotonic, this is equivalent to minimizing the so-called condensed log-likelihood function min_{𝜃∈ℝ𝑑,𝜃𝑗>0}LogML(𝜃)∶=𝑛ln𝜎²(𝜃)+ln(det(𝑅(𝜃))), (2.20) often encountered in the literature.

3. Asymptotic Behavior of the Maximum Likelihood Function—Why Kriging Model Training Is Tricky: Part I

The condition number of a regular matrix 𝑅∈ℝ^𝑛×𝑛 with respect to a given matrix norm ‖⋅‖ is defined as cond(𝑅)∶=‖𝑅‖⋅‖𝑅⁻¹‖. For the matrix norm induced by the euclidean vector norm, one can show that(𝑅)=𝜆cond^max𝜆_min, (3.1) where 𝜆_max, 𝜆_min are the largest, the smallest eigenvalue of 𝑅, respectivley. In order to prevent the solution of the Kriging equation system (2.8) from being spoiled by numerical errors, it is important to prevent the covariance matrix from being severely ill-conditioned. However, the next theorem shows that eventually, when the condition number approaches infinity, so does the associated likelihood function. (Keep in mind that we have formulated likelihood estimation as a minimization problem; see (2.16).) Throughout this section, we will assume that the regression design matrix 𝐹 from (2.5) features the vector 𝟏∈ℝ^𝑛 as first column. This is the case of the highest practical relevance and, in fact, is of particular difficulty, since in this case the first column of 𝐹 coincides with a limit eigenvector of the correlation matrix as will be seen in the following.

Theorem 3.1. Let 𝑥¹,…,𝑥^𝑛∈ℝ^𝑑, 𝑌∶=(𝑦(𝑥¹),…,𝑦(𝑥^𝑛))∈ℝ^𝑛 be a data set of sampled sites and responses. Let 𝑅(𝜃) be the associated spatial correlation matrix, and let Σ(𝜃) be the vector of errors with respect to the chosen regression model. Furthermore, let cond(𝑅(𝜃)) be the condition number of 𝑅(𝜃). Suppose that the following conditions hold true: (1)the eigenvalues 𝜆_𝑖(𝜃) of 𝑅(𝜃) are mutually distinct for small 0<‖𝜃‖≤𝜀, 𝜃∈ℝ^𝑑,(2)the derivatives of the eigenvalues do not vanish in the limit: (𝑑/𝑑𝜏)|_𝜏=0𝜆_𝑗(𝜏𝟏)=𝜆^′_𝑗(0)≠0 for all 𝑗=2,…,𝑛,(3)Σ(0)∉span{𝟏}, Σ(0)∉𝟏^⟂. Then, cond(𝑅(𝜃))⟶∞for‖𝜃‖⟶0, (3.2) and there exist constants 𝑐₁,𝑐₂∈ℝ such that 𝑐₁cond(𝑅(𝜃))≤ML(𝜃)≤𝑐₂cond(𝑅(𝜃)), (3.3) for 𝜃∈ℝ^𝑑,0<‖𝜃‖≤𝜀. The constants 𝑐₁,𝑐₂ are independent of 𝜃.

Remark 3.2. (1) The conditions given in the above theorem cannot be proven to hold true in general, since they depend on the data set in question. However, they hold true for nondegenerate data set. In Appendix A, a relationship between condition 2 and the regularity of 𝑅′(0) is established, giving strong support that condition 2 is generally valid. Concerning the third condition, it will be shown in Lemma 3.4, that the limit Σ(0) exists, given conditions 1 and 2. Note that the set span{𝟏}∪𝟏^⟂ is of Lebesgue measure zero in ℝ^𝑛. In all practical applications known to the author, these conditions were fulfilled.
(2) It holds that lim_{‖𝜃‖→∞}(𝑅(𝜃))_𝑖,𝑗=𝐼∈ℝ^𝑛×𝑛. Hence, the likelihood function approaches a constant limit for ‖𝜃‖→∞ and lim_{‖𝜃‖→∞}cond(𝑅(𝜃))=1. The corresponding predictor behavior is investigated in Section 4.
(3) Even though Theorem 3.1 shows that the model likelihood becomes arbitrarily bad for hyperparameters ‖𝜃‖→0, the optimum might lie very close to the blowup region of the condition number, leading to still quite ill-conditioned covariance matrices [11]. This fact as well as the general behavior of the likelihood function as predicted by Theorem 3.1 is illustrated in Figure 1.
(4) Figure 2(b), provides an additional illustration of Theorem 3.1.
(5) Theorem 3.1 offers a strategy for choosing starting solutions for the optimization problem (2.16): take each 𝜃_𝑘, 𝑘∈1,…,𝑑 as small as possible such that the corresponding correlation matrix is still (numerically) positive definite.
(6) A related investigation of interpolant limits has been performed in [18] but for standard radial basis functions.

(a)

(b)

(c)

(d)

Figure 1 

Ordinary Kriging estimation (a) of a two-dimensional analytical test function (b) based on 15 samples points. ((c) and (d)) Two views of the associated LogML, scaled by a factor of 1/100. This example shows that hyperparameter optima might lie very close to the blowup of the LogML due to ill-conditioning, that is proved to occur by Theorem 3.1. Model function and sample locations are listed in Appendix B.

(a)

(b)

Figure 2 

Kriging of data stemming from reduced-order modeling of solutions to the Navier-Stokes equations via POD coefficient interpolation; see, for example, [23]. (a) Kriging predictors for different choices of the distance weight 𝜃. (b) Corresponding condensed log-likelihood function, showing the limit behavior as predicted by Theorem 3.1. Numerically, the function features no local minimum. Thus, it is impossible to say which one of the estimator functions in the left picture is the most likely.

In order to support readability, we divide the proof of Theorem 3.1 into smaller units, organized as follows. As a starting point, we establish two auxiliary lemmata on the existence of limits of eigenvalue quotients and of errors vectors. Subsequently, the proof of the main theorem is conducted relying on the lemmata.

Lemma 3.3. In the setting of Theorem 3.1, let 𝜆_𝑖(𝜃), 𝑖=1,…,𝑛 be the eigenvalues of 𝑅(𝜃), ordered by size. Then, lim_{‖𝜃‖→0}𝜆_𝑖(𝜃)𝜆_𝑗(𝜃)=const.>0∀𝑖,𝑗=2,…,𝑛. (3.4)

Proof. Because of (2.13), it holds that (𝑅(0))_𝑖,𝑗=𝟏𝟏^𝑇∈ℝ^𝑛×𝑛 for every admissible spatial correlation function. Since (𝟏𝟏^𝑇)𝟏=𝑛𝟏∈ℝ^𝑛 and (𝟏𝟏^𝑇)𝑊=0 for all 𝑊∈𝟏^⟂⊂ℝ^𝑛, the limit eigenvalues of the correlation matrix ordered by size are given by 𝜆₁(0)=𝑛>0=𝜆₂(0)=⋯=𝜆_𝑛(0).
Under the present conditions, the eigenvalues 𝜆_𝑖 are differentiable with respect to 𝜃. Hence, it is sufficient to proof the lemma for ℝ∋𝜏↦𝜃(𝜏)=𝜏𝟏∈ℝ^𝑑 and 𝜏→0. Now, condition 2 and L'Hospital's rule imply the result.

Lemma 3.4. In the setting of Theorem 3.1, let Σ(𝜃) be defined by (2.18) and (2.19). Then, lim_{‖𝜃‖→0}Σ(𝜃)=Σ(0) (3.5) exists.

Proof. We prove Lemma 3.4 by showing that lim_{‖𝜃‖→0}𝛽(𝜃) exists.
Remember that 𝛽=(𝛽₀,…,𝛽_𝑝)^𝑇∈ℝ^𝑝+1, with 𝑝∈ℕ₀ depending on the chosen regression model. As in the above lemma, we can restrict the considerations to the direction 𝜏↦𝜃(𝜏)=𝜏𝟏.
Let 𝜆_𝑖(𝜃),𝑖=1,…,𝑛 be the eigenvalues of 𝑅(𝜃) ordered by size with corresponding eigenvector matrix 𝑄(𝜃)=(𝑋₁,…,𝑋_𝑛)(𝜃) such that 𝑄(𝜃)𝑅(𝜃)𝑄^𝑇(𝜃)=Λ(𝜃)=diag(𝜆₁,…,𝜆_𝑛). For brevity, define 𝑋_𝑖(𝜏)∶=𝑋_𝑖(𝜃(𝜏))=𝑋_𝑖(𝜏𝟏) and so forth.
It holds that 𝑋₁(0)=(1/√𝑛)𝟏; see Lemma 3.3. Hence, ⟨𝟏,𝑋_𝑗(𝜏)⟩→0 for 𝜏→0 and 𝑗=2,…,𝑛. In the present setting, the derivatives of eigenvalues and (normalized) eigenvectors exist and can be extended to 0; see, for example, [22]. By another application of L'Hospital's rule, lim_𝜏→0⟨𝟏,𝑋_𝑗(𝜏)⟩𝜆_𝑗(𝜏)existsfor𝑗=2,…,𝑛. (3.6) Introducing 𝐿(𝜏)=𝜆₁0⋯00𝜆₂⋮⋮⋱00⋯0𝜆₂(𝜏)∈ℝ^{(𝑝+1)×(𝑝+1)}, (3.7) we can restate (2.19) as 𝛽=𝐹^𝑇𝑄Λ⁻¹𝑄^𝑇𝐹⁻¹𝐿⁻¹𝐿𝐹^𝑇𝑄Λ⁻¹𝑄^𝑇𝑌=𝐿𝐹^𝑇𝑄Λ⁻¹𝑄^𝑇𝐹⁻¹𝐿𝐹^𝑇𝑄Λ⁻¹𝑄^𝑇𝑌. (3.8) It is sufficient to show that lim_𝜏→0(𝐿𝐹^𝑇𝑄Λ⁻¹)(𝜏) exists.
Writing columnwise (𝐹₀,𝐹₁,…,𝐹_𝑝)∶=𝐹, a direct computation shows 𝐿𝐹^𝑇𝑄Λ⁻¹(𝜏)=𝜆₁𝜆₁⟨𝐹₀,𝑋₁⟩𝜆₁𝜆₂⟨𝐹₀,𝑋₂𝜆⟩⋯¹𝜆_𝑛⟨𝐹₀,𝑋_𝑛⟩𝜆₂𝜆₁⟨𝐹₁,𝑋₁⟩𝜆₂𝜆₂⟨𝐹₁,𝑋₂𝜆⟩⋯²𝜆_𝑛⟨𝐹₁,𝑋_𝑛⟩𝜆⋮⋮⋮₂𝜆₁⟨𝐹_𝑝,𝑋₁⟩𝜆₂𝜆₂⟨𝐹_𝑝,𝑋₂𝜆⟩⋯²𝜆_𝑛⟨𝐹_𝑝,𝑋_𝑛⟩(𝜏)∈ℝ^{(𝑝+1)×𝑛}. (3.9) Note that 𝐹₀=𝟏=√𝑛𝑋₁(0) for the default choices of regression basis functions, such that ⟨𝐹₀,𝑋₁(0)⟩≠0 and ⟨𝐹₀,𝑋_𝑗(0)⟩=0 for 𝑗=2,…,𝑛. The desired result follows from (3.6) and Lemma 3.3.
Remark 3.5. Actually, one cannot prove for (𝐿𝐹^𝑇𝑄Λ⁻¹) to be regular in general, since this matrix depends on the chosen sample locations. It might be possible to artificially choose samples such that, for example, 𝐹 has not full rank. Yet if so, the whole Kriging exercise cannot be performed, since (2.19) is not well defined in this case. For constant regression, that is, 𝐹=𝟏, this is impossible. Note that 𝐹 is independent of 𝜃.

Now, let us prove Theorem 3.1 using notation as introduced above.

Proof. As shown in the proof of Lemma 3.3(𝑅(𝜃=𝜆cond))^max(𝜃)𝜆_min(𝜃)⟶∞for‖𝜃‖⟶0. (3.10) Because the correlation matrix is symmetric and positive definite, a decomposition 𝑅(𝜃)⁻¹=𝑄(𝜃)Λ(𝜃)⁻¹𝑄^𝑇(𝜃), (3.11) with 𝑄 orthogonal and Λ⁻¹=Diag(1/𝜆_𝑖)_{𝑖={1,…,𝑛}}, exists.
If necessary, renumber such that 𝜆_max=𝜆₁≥⋯≥𝜆_𝑛=𝜆_min. Let 𝑊(𝜃)=(𝑊₁(𝜃),…,𝑊_𝑛(𝜃))∶=𝑄(𝜃)Σ(𝜃). By Lemma 3.4, 𝑊(0)∶=𝑄(0)^𝑇Σ(0) exists. Condition 3 insures that 𝑊₁(0)≠0.
Case 1. Suppose that 𝑊_𝑖(0)≠0 for all 𝑖=1,…,𝑛.
By continuity, 𝑊_𝑖(𝜃)≠0 for 0≤‖𝜃‖≤𝜀 and 𝜀>0 small enough. Since {𝜃∈ℝ^𝑑,0≤‖𝜃‖≤𝜀} is a compact set, 𝑊²_𝑚∶=min_{0≤‖𝜃‖≤𝜀}𝑊²_𝑖(𝜃),𝑖=1,…,𝑛 (3.12) exists. Then, for ‖𝜃‖∈[0,𝜖], (𝜃)=1ML𝑛^𝑛Σ^𝑇(𝜃)𝑅(𝜃)⁻¹Σ(𝜃)^𝑛(𝑅(𝜃=1det))𝑛^𝑛_𝑖𝑊²_𝑖(𝜃)𝜆_𝑖(𝜃)^𝑛_𝑖𝜆_𝑖(𝜃)≥𝑊²_𝑚𝑛^𝑛𝑛−1𝜆_max+1(𝜃)𝜆_min(𝜃)^𝑛𝜆^𝑛−1_min(𝜃)𝜆_max(𝜃)≥𝑊²_𝑚𝑛^𝑛𝑛−1𝜆_max(𝜃)^𝑛+1𝜆_min(𝜃)^𝑛𝜆^𝑛−1_min(𝜃)𝜆_max(𝜃)≥𝑊²_𝑚𝑛^𝑛(𝑛−1)^𝑛cond(𝑅(𝜃))^𝑛−1(𝑅(𝜃≥𝑊+cond))_2𝑚𝑛^𝑛cond(𝑅(𝜃)). (3.13)Case 2. Suppose that Case 1 does not hold true.
From Σ(0)∉span{𝟏}, it follows that 𝑊(𝜃)=𝑄(𝜃)^𝑇Σ(𝜃)∉span(1,0,…,0)^𝑇, (3.14) for ‖𝜃‖ sufficiently small. Let 𝐽∶={𝑖∈{1,…,𝑛}∣𝑊_𝑖(0)≠0}. Then, 𝑛_𝐽∶=|𝐽|≥2.
Define 𝑊²_𝑚∶=min_{0≤‖𝜃‖≤𝜀,𝑗∈𝐽}𝑊²_𝑗(𝜃). (3.15) For the index 𝑚 defined by 𝜆_𝑚∶=min_𝑗∈𝐽{𝜆_𝑗}, it holds that 𝑚∈{2,…,𝑛}.
By Lemma 3.3, 𝐿∶=min_{0≤‖𝜃‖≤𝜀}𝜆_min(𝜃)𝜆_𝑚(𝜃)>0 (3.16) exists. Using ^𝑛_𝑖=1𝑊²_𝑖(𝜃)𝜆_𝑖(𝜃)≥_𝑗∈𝐽𝑊²_𝑗(𝜃)𝜆_𝑗(𝜃)≥𝑊²_𝑚𝑛_𝐽−1𝜆_max(𝜃)+1𝜆_𝑚(𝜃), (3.17) together with 𝜆_max𝜆_𝑚(𝜃𝜆)^min𝜆_𝑚^𝑛−1(𝜃)=𝜆_max𝜆_min(𝜃𝜆)^min𝜆_𝑚^𝑛(𝜃)≥𝐿^𝑛cond(𝑅(𝜃)), (3.18) the result can be established as in Case 1.The estimate of the upper bound of ML is obtained in an analogous manner. Let 𝑊²_𝑀∶=max_{0≤‖𝜃‖≤𝜀}𝑊²_𝑖(𝜃)∣𝑊_𝑖(𝜃)𝐿≠0,𝑖=1,…,𝑛>0,_𝑀∶=max_{0≤‖𝜃‖≤𝜀}𝜆₂(𝜃)𝜆_min(𝜃)>0. (3.19) Then, (𝜃)=1ML𝑛^𝑛Σ^𝑇(𝜃)𝑅(𝜃)⁻¹Σ(𝜃)^𝑛(𝑅(𝜃=1det))𝑛^𝑛_𝑖𝑊²_𝑖(𝜃)𝜆_𝑖(𝜃)^𝑛_𝑖𝜆_𝑖(𝜃)≤𝑊²_𝑀𝑛^𝑛𝑛−1𝜆_min+1(𝜃)𝜆_max(𝜃)^𝑛𝜆_max(𝜃)𝜆₂^𝑛−1(𝜃)=𝑊²_𝑀𝑛^𝑛𝑛−1𝜆_min(𝜃)^𝑛11+𝜆𝑛−1^min(𝜃)𝜆_max(𝜃)^𝑛𝜆_max(𝜃)𝜆₂^𝑛−1(𝜃)_(∗)≤𝑊²_𝑀^𝑛11−𝑛^𝑛(𝑅(𝜃𝜆cond))²(𝜃)𝜆_min(𝜃)^𝑛−1𝑛1+()(𝑅(𝜃≤𝑊𝑛−1cond))_𝑀2𝑛𝑒𝐿_𝑀^𝑛−1((𝑅(𝜃)𝑊cond))+2≤3^𝑀2𝑛𝑒𝐿_𝑀^𝑛−1cond(𝑅(𝜃)), (3.20) where we used Bernoulli's inequality at (∗), Lemma 3.3, and the fact that 𝑛/(𝑛−1)≤2.

Remark 3.6. The extension of the main theorem to Cokriging prediction [7, 20] is a straight-forward exercise, since the limit eigenvectors of the Cokriging correlation matrix corresponding to nonzero eigenvalues can also be determined explicitly.

4. Why Kriging Model Training Is Tricky: Part II

The following simple observation illustrates Kriging predictor behavior for large-distance weights 𝜃. Notation is to be understood as introduced in Section 3.

Observation 1. Suppose that sample locations {𝑥¹,…,𝑥^𝑛}⊂ℝ^𝑑 and responses 𝑦_𝑖=𝑦(𝑥^𝑖)∈ℝ, 𝑖=1,…,𝑛 are given. Let 𝑦 be the corresponding Kriging predictor according to (2.10). Then, for ℝ^𝑑∋𝑥∉{𝑥¹,…,𝑥^𝑛} and distance weights ‖𝜃‖→∞, it holds that 𝑦(𝑥)⟶𝑓(𝑥)𝛽. (4.1)

Put in simple words: if too large distance weights are chosen, then the resulting predictor function has the shape of the regression model, 𝑥↦𝑓(𝑥)𝛽, with peaks at the sample sites, compare Figure 2, dashed curve.

Proof. According to (2.10) it holds that 𝑦(𝑥)=𝑓(𝑥)𝛽+𝑐^𝑇_𝜃(𝑥)𝐶_𝜃⁻¹𝑌−𝐹𝛽, (4.2) where 𝐶=𝜎²𝑅. By (2.13), it holds that 𝑅(𝜃)→𝐼, for ‖𝜃‖→∞ for every admissible spatial correlation model of the form of (2.11). By Cauchy-Schwartz' inequality, 𝑐^𝑇_𝜃(𝑥)𝐶_𝜃⁻¹𝑌−𝐹𝛽=𝑐_𝜃(𝑥),𝐶_𝜃⁻¹𝑌−𝐹𝛽≤‖𝑐_𝜃(𝑥)‖‖𝐶_𝜃⁻¹𝑌−𝐹𝛽‖, (4.3) where ‖𝑐_𝜃(𝑥)‖→0 and ‖𝐶_𝜃⁻¹(𝑌−𝐹𝛽)‖→const. for ‖𝜃‖→∞.

Remark 4.1. The same predictor behavior arises at locations far away from the sampled sites, that is, for dist(𝑥,{𝑥¹,…,𝑥^𝑛})→∞. This has to be considered, when extrapolating beyond the sample data set.

Figure 2 shows an example data set for which the Kriging maximum likelihood function is constant over a large range of 𝜃 values. This example was not constructed artificially but occured in the author's daily work of computing approximate fluid flow solutions based on proper orthogonal decomposition (POD) followed by coefficient interpolation as described in [23, 24].

The sample data set is given in Table 1. The Kriging estimator given by the dashed line shows a behavior as predicted by Observation 1. Note that from the model training point of view, all distance weights 𝜃>1 are equally likely, yet lead to quite different predictor functions. Since the ML features no local minimum, classical hyperparameter estimation is impossible.

Table 1 

Sampled sites corresponding to the example displayed in Figure 2.

Appendices

A. On the Validity of Condition 2 in Theorem 3.1

The next lemma strongly indicates that the second condition in the main Theorem 3.1 is given in nondegenerate cases.

Lemma A.1. Let ℝ^𝑑∋𝜃↦𝑅(𝜃)∈ℝ^𝑛×𝑛 be the correlation matrix function corresponding to a given set of Kriging data and a fixed spatial correlation model.
Let 𝜆_𝑖(𝜃), 𝑖=1,…,𝑛 be the eigenvalues of 𝑅 ordered by size with corresponding eigenvector matrix 𝑄=(𝑋₁,…,𝑋_𝑛), and define 𝜃∶ℝ→ℝ^𝑑, 𝜏↦𝜃(𝜏)∶=𝜏𝟏. Suppose that the eigenvalues are mutually distinct for 𝜏>0 close to zero.
Denote the directional derivative in the direction 𝟏 with respect to 𝜏 by a prime ′, that is, (𝑑/𝑑𝜏)𝑅(𝜏𝟏)=𝑅′(𝜏) and so forth. Then, it holds that 𝜆^′_𝑖(0)=𝑋^𝑇_𝑖(0)𝑅′(0)𝑋_𝑖. (A.1) If 𝑅′(0) is regular, then 𝜆^′_𝑖(0)≠0, (A.2) for all 𝑖=1,…,𝑛, with at most one possible exception.

Proof. For every admissible spatial spatial correlation function 𝑟(𝜃,⋅,⋅) of the form (2.11) and 𝑥≠𝑧∈ℝ^𝑑, it holds that 𝑟(𝜃(𝜏),𝑥,𝑧)⟶1,for𝜏⟶0. (A.3) Thus, (𝑅(0))_𝑖,𝑗=𝟏𝟏^𝑇∈ℝ^𝑛×𝑛.
It holds that (𝟏𝟏^𝑇)𝟏=𝑛𝟏∈ℝ^𝑛 and (𝟏𝟏^𝑇)𝑊=0 for all 𝑊∈𝟏^⟂⊂ℝ^𝑛; therefore, the limits of the eigenvalues of the correlation matrix ordered by size are given by 𝜆₁(0)=𝑛>𝜆₂(0)=⋯=𝜆_𝑛(0)=0. The assumption, that no multiple eigenvalues occur, ensures that the eigenvalues 𝜆_𝑖 and corresponding (normalized, oriented) eigenvectors 𝑋_𝑖 are differentiable with respect to 𝜏. Let 𝑄(𝜏)=(𝑋₁(𝜏),…,𝑋_𝑛(𝜏))∈ℝ^𝑛×𝑛 be the (orthogonal) matrix of eigenvectors, such that Λ(𝜏)∶=diag𝜆_𝑖(𝜏)^𝑛₁=𝑄^𝑇(𝜏)𝑅(𝜏)𝑄(𝜏). (A.4) Then, Λ^′(𝜏)=𝑄(𝜏)^𝑇𝑅^′(𝜏)𝑄(𝜏)+𝑄^′(𝜏)^𝑇𝑅^′(𝜏)𝑄(𝜏)+𝑄(𝜏)^𝑇𝑅(𝜏)𝑄^′(𝜏),Λ′(0)=𝑄(0)^𝑇𝑅^′(0)𝑄(0)+0𝜆₁(0)⟨𝑋^′₂(0),𝑋₁(0)⟩⋯𝜆₁(0)⟨𝑋^′_𝑛(0),𝑋₁(0)⟩𝜆₁(0)⟨𝑋^′₂(0),𝑋₁(0)𝜆⟩0⋯0⋮⋮⋮₁(0)⟨𝑋^′_𝑛(0),𝑋₁(0),⟩0⋯0 (A.5) where 𝑄(0) is the continuous extension of 𝑄(𝜏); see, for example, [22].
Hence, 𝜆^′_𝑖(0)=𝑒^𝑇_𝑖Λ^′(0)𝑒_𝑖=𝑋^𝑇_𝑖(0)𝑅^′(0)𝑋_𝑖(0). (A.6)
Let us assume, that there exist two indices 𝑗₀,𝑘₀,𝑗₀≠𝑘₀ such that 𝜆^′_𝑘0(0)=0=𝜆^′_𝑗0(0).
Let 𝑊∶=(0,…,−𝜆₁⟨𝑋^′_𝑗0(0),𝑋₁(0)⟩,…,−𝜆₁⟨𝑋_𝑘0′(0),𝑋₁(0),…,0)^𝑇∈ℝ^𝑛. Then, 0⋮0=Λ^′(0)𝑊=𝑄(0)^𝑇𝑅^′(0)𝑄(0)𝑊, (A.7) contradicting the regularity of 𝑅′(0), if 𝑊≠0.
If 𝑊=0, replace 𝑊 by 𝑊∶=(0,…,1^𝑗0,…,1^𝑘0,…,0) and repeat the above argument.

For most correlation models, the derivative 𝑅′(0) can be computed explicitly.

B. Test Setting Corresponding to Figure 1

In order to produce Figure 1, the following test function has been applied:𝑦∶[0,100]×[0,100]⟶ℝ,𝑥₁,𝑥₂𝑥⟼5+²¹⋅𝑥₂𝑥10,000sin².10 (B.1)

The Kriging predictor function displayed in this figure has been constructed based on the fifteen (randomly chosen) sample points shown in Table 2.

Table 2 

Nomenclature

𝑑∈ℕ:	Dimension of parameter space
𝑛∈ℕ:	(Fixed) number of sample points
𝑥𝑖∈ℝ𝑑:	𝑖th sample location
𝑦𝑖∈ℝ:	Sample value at sample location 𝑥𝑖
𝐼∈ℝ𝑛×𝑛:	Unit matrix
𝟏∶=(1,…,1)𝑇∈ℝ𝑛:	Vector with all entries equal to 1
𝑒𝑖=(0,…,0,1𝑖,0,…,0)𝑇:	𝑖th standard basis vector
𝑉⟂⊂ℝ𝑛:	Subspace of all vectors orthogonal to 𝑉∈ℝ𝑛
𝑅∈ℝ𝑛×𝑛:	Correlation matrix
𝐶∈ℝ𝑛×𝑛:	Covariance matrix
cond(𝑅):	Condition number of 𝑅∈ℝ𝑛×𝑛
⟨⋅,⋅⟩:	Euclidean scalar product
𝑒=exp(1):	Euler’s number
𝑝+1∈ℕ:	Dimension of regression model
𝛽∈ℝ𝑝+1:	Vector of regression coefficients
𝑓∶ℝ𝑑→ℝ𝑝+1:	Regression model
𝜖∶ℝ𝑑→ℝ:	Random error function
𝐸[⋅]:	Expectation value
𝜎=√Var:	Standard deviation
𝜃∈ℝ𝑑:	Distance weights vector, model hyperparameters.

Acknowledgments

This research was partly sponsored by the European Regional Development Fund and Economic Development Fund of the Federal German State of Lower Saxony Contract/Grant no. W3-80026826.