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

Authors: Ralf Zimmermann
Journal: Journal of Applied Mathematics
 (2010)

---

## 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.


---

## Body

## 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-calledhyperparameters. 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 Section3. 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)y:Rd⊇U⟶R,x⟼y(x),
based on a finite data set of sample locations x1,…,xn∈U⊂ℝd with corresponding responses y1:=y(x1),…,yn:=y(xn)∈ℝ obtained from measurements or numerical computations. The collection of responses is denoted by(2.2)YT=(y1,…,yn)∈Rn.
The function y:U→ℝ to be estimated is assumed to be the realization of an underlying random process given by a regression model and a random error function ϵ(x) with zero mean. More precisely(2.3)y(x)=f(x)β+ϵ(x)=(f0(x),…,fp(x))(β0⋮βp)+ϵ(x),
where the components of the row vector function f:ℝd→ℝp+1 are the basis functions of the regression model and β=(β0,…,βp) is the corresponding vector of regression coefficients. By assumption,(2.4)E[ϵ(x)]=0∀x.
The component functions of f 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): p=0, f:ℝd→ℝ,x↦1, f(x)β=β∈ℝ,(2)
linear regression (universal Kriging): p=d, f:ℝd→ℝd+1,x↦(1,x1,…,xd), f(x)β=β0+β1x1+⋯+βdxd∈ℝ,and higher-order polynomials.Introducing the regression design matrix(2.5)F:=(f(x1)⋮f(xn))=(f0(x1)f1(x1)⋯fp(x1)⋮⋮⋯⋮f0(xn)f1(xn)⋯fp(xn))∈Rn×(p+1),
the vector of errors at the sampled sites can be written as(2.6)Σ:=(ϵ1⋮ϵn)=(ϵ(x1)⋮ϵ(xn))=(y1-f(x1)β⋮yn-f(xn)β)=Y-Fβ.
Note that the first column of F equals 1∈ℝn for all polynomial regression models.The Kriging predictorŷ estimates y at an untried site x as a linear combination of the sampled data(2.7)ŷ(x)=〈ω(x),Y〉=∑i=1nωi(x)yi.
For each x∈ℝd, the unique vector of weights ω(x)=(ω1(x),…,ωn(x)) that leads to an unbiased prediction minimizing the mean squared error is given by the solution of the Kriging equation system(2.8)(CFFT0)(ω(x)12μ(x))=(c(x)f(x)T)∈Rn+(p+1).
Here,(2.9)C:=(Cov(ϵ(xi),ϵ(xj)))i,j≤n∈Rn×n,c(x):=(Cov(ϵ(xi),ϵ(x)))i≤n∈Rn
denote the covariance matrix and the covariance vector, respectively, and the entries of the vector μ=(μ0,…,μp) are Lagrange multipliers. Solving (2.8) by Schur matrix complement inversion and substituting in (2.7) leads to the Kriging predictor formula(2.10)ŷ(x)=f(x)β+cT(x)C-1(Y-Fβ),
where β=(FTR-1F)-1FTR-1Y is the generalized least squares solution to the regression problem Fβ≃Y. 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 form(2.11)Cov(ϵ(xi),ϵ(xj))=σ2r(θ,xi,xj)=σ2∏k=1dscfk(θ,xi,xj).
Here θ=(θ1,…,θd)∈ℝd 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)R(θ)=(r(θ,xi,xj))i,j∈Rn×n.
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(2.13)r(θ,p,q)⟶{1,for∥θ∥⟶0,0,for∥θ∥⟶∞,at distinct locationsp≠q.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 by(2.14)scfk(θ,xi,xj)=exp⁡(-θk|xki-xkj|2),r(θ,xi,xk)=exp⁡(-∑kθk|xki-xkj|2).
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 xi, xj. In this setting, hyper-parameter training for Kriging models consists of optimizing the corresponding maximum likelihood function(2.15)MaxL(σ,β,θ)=12πdet⁡(σ2R(θ))exp⁡(-12σ2(Y-Fβ.)TR(θ)-1(Y-Fβ.)).
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:(2.16)min⁡{θ∈Rd,θj>0}ML(θ):=(σ2(θ))ndet⁡(R(θ)),
where the dependency on θ is as follows:(2.17)σ2(θ)=1nΣT(θ)R(θ)-1Σ(θ),(2.18)Σ(θ)=Y-E[Y]=Y-Fβ(θ),(2.19)β(θ)=(FTR-1(θ)F)-1FTR-1(θ)Y.
Because the logarithm is monotonic, this is equivalent to minimizing the so-called condensed log-likelihood function(2.20)min⁡{θ∈Rd,θj>0}LogML(θ):=nln⁡(σ2(θ))+ln⁡(det⁡(R(θ))),
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 matrixR∈ℝn×n with respect to a given matrix norm ∥·∥ is defined as cond(R):=∥R∥·∥R-1∥. For the matrix norm induced by the euclidean vector norm, one can show that(3.1)cond(R)=λmax⁡λmin⁡,
where λmax, λmin are the largest, the smallest eigenvalue of R, 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 F from (2.5) features the vector 1∈ℝn 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 F coincides with a limit eigenvector of the correlation matrix as will be seen in the following.Theorem 3.1.
Letx1,…,xn∈ℝd, Y:=(y(x1),…,y(xn))∈ℝn be a data set of sampled sites and responses. Let R(θ) be the associated spatial correlation matrix, and let Σ(θ) be the vector of errors with respect to the chosen regression model. Furthermore, let cond(R(θ)) be the condition number of R(θ). Suppose that the following conditions hold true: (1)
the eigenvaluesλi(θ) of R(θ) are mutually distinct for small 0<∥θ∥≤ɛ, θ∈ℝd,(2)
the derivatives of the eigenvalues do not vanish in the limit:(d/dτ)|τ=0λj(τ1)=λj′(0)≠0 for all j=2,…,n,(3)
Σ(0)∉span{1}, Σ(0)∉1⊥.Then,
(3.2)cond(R(θ))⟶∞for∥θ∥⟶0,
and there exist constants c1,c2∈ℝ such that
(3.3)c1cond(R(θ))≤ML(θ)≤c2cond(R(θ)),
for θ∈ℝd,0<∥θ∥≤ɛ. The constants c1,c2 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 AppendixA, a relationship between condition 2 and the regularity of R′(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{1}∪1⊥ is of Lebesgue measure zero in ℝn. In all practical applications known to the author, these conditions were fulfilled.
(2) It holds thatlim∥θ∥→∞(R(θ))i,j=I∈ℝn×n. Hence, the likelihood function approaches a constant limit for ∥θ∥→∞ and lim∥θ∥→∞cond(R(θ))=1. The corresponding predictor behavior is investigated in Section 4.
(3) Even though Theorem3.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) Figure2(b), provides an additional illustration of Theorem 3.1.
(5) Theorem3.1 offers a strategy for choosing starting solutions for the optimization problem (2.16): take each θk, k∈1,…,d 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.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 of1/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)(c)(d)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.
(a)(b)In order to support readability, we divide the proof of Theorem3.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 Theorem3.1, let λi(θ), i=1,…,n be the eigenvalues of R(θ), ordered by size. Then,
(3.4)lim⁡∥θ∥→0λi(θ)λj(θ)=const.>0∀i,j=2,…,n.Proof.
Because of (2.13), it holds that (R(0))i,j=11T∈ℝn×n for every admissible spatial correlation function. Since (11T)1=n1∈ℝn and (11T)W=0 for all W∈1⊥⊂ℝn, the limit eigenvalues of the correlation matrix ordered by size are given by λ1(0)=n>0=λ2(0)=⋯=λn(0).
Under the present conditions, the eigenvaluesλi are differentiable with respect to θ. Hence, it is sufficient to proof the lemma for ℝ∋τ↦θ(τ)=τ1∈ℝd and τ→0. Now, condition 2 and L'Hospital's rule imply the result.Lemma 3.4.
In the setting of Theorem3.1, let Σ(θ) be defined by (2.18) and (2.19). Then,
(3.5)lim⁡∥θ∥→0Σ(θ)=Σ(0)
exists.Proof.
We prove Lemma3.4 by showing that lim∥θ∥→0β(θ) exists.
Remember thatβ=(β0,…,βp)T∈ℝp+1, with p∈ℕ0 depending on the chosen regression model. As in the above lemma, we can restrict the considerations to the direction τ↦θ(τ)=τ1.
Letλi(θ),i=1,…,n be the eigenvalues of R(θ) ordered by size with corresponding eigenvector matrix Q(θ)=(X1,…,Xn)(θ) such that Q(θ)R(θ)QT(θ)=Λ(θ)=diag(λ1,…,λn). For brevity, define Xi(τ):=Xi(θ(τ))=Xi(τ1) and so forth.
It holds thatX1(0)=(1/n)1; see Lemma 3.3. Hence, 〈1,Xj(τ)〉→0 for τ→0 and j=2,…,n. 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,
(3.6)lim⁡τ→0〈1,Xj(τ)〉λj(τ)existsforj=2,…,n.
Introducing
(3.7)L(τ)=(λ10⋯00λ2⋮⋮⋱00⋯0λ2)(τ)∈R(p+1)×(p+1),
we can restate (2.19) as
(3.8)β=(FTQΛ-1QTF)-1(L-1L)FTQΛ-1QTY=(LFTQΛ-1QTF)-1(LFTQΛ-1QT)Y.
It is sufficient to show that limτ→0(LFTQΛ-1)(τ) exists.
Writing columnwise(F0,F1,…,Fp):=F, a direct computation shows
(3.9)(LFTQΛ-1)(τ)=(λ1λ1〈F0,X1〉λ1λ2〈F0,X2〉⋯λ1λn〈F0,Xn〉λ2λ1〈F1,X1〉λ2λ2〈F1,X2〉⋯λ2λn〈F1,Xn〉⋮⋮⋮λ2λ1〈Fp,X1〉λ2λ2〈Fp,X2〉⋯λ2λn〈Fp,Xn〉)(τ)∈R(p+1)×n.
Note that F0=1=nX1(0) for the default choices of regression basis functions, such that 〈F0,X1(0)〉≠0 and 〈F0,Xj(0)〉=0 for j=2,…,n. The desired result follows from (3.6) and Lemma 3.3.
Remark 3.5.
Actually, one cannot prove for(LFTQΛ-1) 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, F 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, F=1, this is impossible. Note that F is independent of θ.Now, let us prove Theorem3.1 using notation as introduced above.Proof.
As shown in the proof of Lemma3.3(3.10)cond(R(θ))=λmax⁡(θ)λmin⁡(θ)⟶∞for∥θ∥⟶0.
Because the correlation matrix is symmetric and positive definite, a decomposition
(3.11)R(θ)-1=Q(θ)Λ(θ)-1QT(θ),
with Q orthogonal and Λ-1=Diag(1/λi)i={1,…,n}, exists.
If necessary, renumber such thatλmax=λ1≥⋯≥λn=λmin. Let W(θ)=(W1(θ),…,Wn(θ)):=Q(θ)Σ(θ). By Lemma 3.4, W(0):=Q(0)TΣ(0) exists. Condition 3 insures that W1(0)≠0.
Case 1.
Suppose thatWi(0)≠0 for all i=1,…,n.
By continuity,Wi(θ)≠0 for 0≤∥θ∥≤ɛ and ɛ>0 small enough. Since {θ∈ℝd,0≤∥θ∥≤ɛ} is a compact set,
(3.12)Wm2:=min⁡0≤∥θ∥≤ɛ{Wi2(θ),i=1,…,n}
exists. Then, for ∥θ∥∈[0,ϵ],
(3.13)ML(θ)=1nn(ΣT(θ)R(θ)-1Σ(θ))ndet⁡(R(θ))=1nn(∑iWi2(θ)λi(θ))n∏iλi(θ)≥(Wm2n)n(n-1λmax⁡(θ)+1λmin⁡(θ))nλmin⁡n-1(θ)λmax⁡(θ)≥(Wm2n)n((n-1λmax⁡(θ))n+(1λmin⁡(θ))n)λmin⁡n-1(θ)λmax⁡(θ)≥(Wm2n)n((n-1)ncond(R(θ))n-1+cond(R(θ)))≥(Wm2n)ncond(R(θ)).Case 2.
Suppose that Case1 does not hold true.
FromΣ(0)∉span{1}, it follows that
(3.14)W(θ)=Q(θ)TΣ(θ)∉span⁡{(1,0,…,0)T},
for ∥θ∥ sufficiently small. Let J:={i∈{1,…,n}∣Wi(0)≠0}. Then, nJ:=|J|≥2.
Define(3.15)Wm2:=min⁡0≤∥θ∥≤ɛ,j∈J{Wj2(θ)}.
For the index m̃ defined by λm̃:=minj∈J{λj}, it holds that m̃∈{2,…,n}.
By Lemma3.3,
(3.16)L:=min⁡0≤∥θ∥≤ɛ{λmin⁡(θ)λm̃(θ)}>0
exists. Using
(3.17)∑i=1nWi2(θ)λi(θ)≥∑j∈JWj2(θ)λj(θ)≥Wm2(nJ-1λmax⁡(θ)+1λm̃(θ)),
together with
(3.18)λmax⁡λm̃(θ)(λmin⁡λm̃)n-1(θ)=λmax⁡λmin⁡(θ)(λmin⁡λm̃)n(θ)≥Lncond(R(θ)),
the result can be established as in Case 1.The estimate of the upper bound of ML is obtained in an analogous manner. Let
(3.19)WM2:=max⁡0≤∥θ∥≤ɛ{Wi2(θ)∣Wi(θ)≠0,i=1,…,n}>0,LM:=max⁡0≤∥θ∥≤ɛ{λ2(θ)λmin⁡(θ)}>0.
Then,
(3.20)ML(θ)=1nn(ΣT(θ)R(θ)-1Σ(θ))ndet⁡(R(θ))=1nn(∑iWi2(θ)λi(θ))n∏iλi(θ)≤(WM2n)n(n-1λmin⁡(θ)+1λmax⁡(θ))nλmax⁡(θ)λ2n-1(θ)=(WM2n)n(n-1λmin⁡(θ))n(1+1n-1λmin⁡(θ)λmax⁡(θ))nλmax⁡(θ)λ2n-1(θ)≤(*)(WM2)n(1-1n)ncond(R(θ))(λ2(θ)λmin⁡(θ))n-1(1+n(n-1)cond(R(θ)))≤WM2neLMn-1(cond(R(θ))+2)≤3WM2neLMn-1cond(R(θ)),
where we used Bernoulli's inequality at (*), Lemma 3.3, and the fact that n/(n-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{x1,…,xn}⊂ℝd and responses yi=y(xi)∈ℝ, i=1,…,n are given. Let ŷ be the corresponding Kriging predictor according to (2.10). Then, for ℝd∋x∉{x1,…,xn} and distance weights ∥θ∥→∞, it holds that
(4.1)ŷ(x)⟶f(x)β.Put in simple words: if too large distance weights are chosen, then the resulting predictor function has the shape of the regression model,x↦f(x)β, with peaks at the sample sites, compare Figure 2, dashed curve.Proof.
According to (2.10) it holds that
(4.2)ŷ(x)=f(x)β+cθT(x)Cθ-1(Y-Fβ),
where C=σ2R. By (2.13), it holds that R(θ)→I, for ∥θ∥→∞ for every admissible spatial correlation model of the form of (2.11). By Cauchy-Schwartz' inequality,
(4.3)|cθT(x)Cθ-1(Y-Fβ)|=|〈cθ(x),Cθ-1(Y-Fβ)〉|≤∥cθ(x)∥∥Cθ-1(Y-Fβ)∥,
where ∥cθ(x)∥→0 and ∥Cθ-1(Y-Fβ)∥→const. for ∥θ∥→∞.Remark 4.1.
The same predictor behavior arises at locations far away from the sampled sites, that is, fordist(x,{x1,…,xn})→∞. This has to be considered, when extrapolating beyond the sample data set.Figure2 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 Table1. 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 Figure2.
x=α:0.02.04.06.08.0y(x):−0.2293340.2770180.455534−0.7695580.26634



---

*Source: 494070-2011-01-13.xml*