Fixed-domain Posterior Contraction Rates for Spatial Gaussian Process Model with Nugget

1]Cheng Li<EMAIL_ADDRESS>1]Saifei Sun<EMAIL_ADDRESS>[1]Department of Statistics and Data Science, National University of Singapore
2]Yichen Zhu<EMAIL_ADDRESS>[2]Department of Statistical Science, Duke University

Spatial Gaussian process regression models typically contain finite dimensional covariance parameters that need to be estimated from the data. We study the Bayesian estimation of covariance parameters including the nugget parameter in a general class of stationary covariance functions under fixed-domain asymptotics, which is theoretically challenging due to the increasingly strong dependence among spatial observations. We propose a novel adaptation of the Schwartz's consistency theorem for showing posterior contraction rates of the covariance parameters including the nugget. We derive a new polynomial evidence lower bound, and propose consistent higher-order quadratic variation estimators that satisfy concentration inequalities with exponentially small tails. Our Bayesian fixed-domain asymptotics theory leads to explicit posterior contraction rates for the microergodic and nugget parameters in the isotropic Matérn covariance function under a general stratified sampling design. We verify our theory and the Bayesian predictive performance in simulation studies and an application to sea surface temperature data.

Keywords: Bayesian inference, Matérn covariance function, Evidence lower bound, Higher-order quadratic variation.

§ INTRODUCTION

Gaussian process has been widely used in spatial statistics for modeling spatial and spatiotemporal correlations. This paper studies the parameter estimation in the following spatial Gaussian process regression model:

\begin{align} \label{eq:obs.model}
& Y(\bfs) = \ff(\bfs)^\T  \beta + X(\bfs) + \varepsilon(\bfs), \text{ for } \bfs\in [0,1]^d,
\end{align}

where $d$ is the dimension of domain, $\ff(\cdot)=(\ff_1(\cdot),\ldots,\ff_p(\cdot))^\T  \in \RR^p$ is a $p$-dimensional vector of deterministic known functions on $[0,1]^d$, and $\beta\in \RR^p$ is the regression coefficient vector. In applications, $\ff_1,\ldots,\ff_p$ can include the constant function $1$, and hence $\beta$ can include an intercept term. The term $X(\cdot)$ is a Gaussian process $X=\left\{X(\bfs): \bfs\in [0,1]^d \right\}$, and $\{\varepsilon(\bfs): \bfs\in [0,1]^d\}$ is a Gaussian white noise process independent of $X(\cdot)$ that satisfies $\varepsilon(\bfs)\sim \Ncal(0,\tau)$ for all $\bfs\in [0,1]^d$. In practice, we observe the process $Y(\cdot)$ and the functions $\ff(\cdot)$ on a set of distinct sampling points $S_n=\{\bfs_1,\ldots,\bfs_n\}\subseteq [0,1]^d$. The observed data from the model (<ref>) are $Y_n=(Y(\bfs_1),\ldots,Y(\bfs_n))^\T $ and $ F_n  = (\ff(\bfs_1)^\T ,\ldots,\ff(\bfs_n)^\T ) \in \RR^{n\times p}$. All statistical inference is based on the data $(Y_n, F_n)$.

Model (<ref>) is one of the most important models in spatial statistics and can be used as building blocks for other more sophisticated models. Besides the regression term $\ff(\bfs)^\T  \beta$, the Gaussian process term $X(\cdot)$ captures the spatial association and explains the spatial random effects from unmeasured or unobserved covariates with spatial pattern (6). The noise $\varepsilon(\cdot)$ captures the measurement error, whose variance parameter $\tau$ is also known as the nugget parameter.

It remains to specify the Gaussian process for $X(\cdot)$. We assume that the mean function of $X(\cdot)$ is zero, and the covariance function of $X(\cdot)$ can be written as $\Cov(X(\bfs),X(\bft)) = \theta K_{\alpha,\nu}(\bfs-\bft)$ for $\bfs,\bft\in [0,1]^d$, where $\theta>0$ is the microergodic parameter that will be defined later, $K_{\alpha,\nu}(\cdot):\RR^d\to\RR$ is a positive definite function, $\alpha$ is the spatial range parameter, and $\nu$ is the smoothness parameter. We will give several examples of covariance functions with this form in Section <ref>. For example, one popular choice in spatial statistics is the Matérn covariance function (42)

\begin{align} \label{eq:MaternCov}
& \Cov(X(\bfs),X(\bft)) = \sigma^2\frac{2^{1-\nu}}{\Gamma(\nu)} \left(\alpha\|\bfs-\bft\|\right)^{\nu} \Kcal_{\nu}\left(\alpha\|\bfs-\bft\|\right),
\end{align}

where $\nu>0$ is the smoothness parameter, $\sigma^2>0$ is the variance (partial sill) parameter, $\alpha>0$ is the inverse range (or length-scale) parameter, $\Gamma(\cdot)$ is the gamma function, $\Kcal_{\nu}(\cdot)$ is the modified Bessel function of the second kind, and $\|\cdot\|$ is the Euclidean norm. If we let $\theta=\sigma^2\alpha^{2\nu}$  and $K(x)=\{2^{1-\nu}/\Gamma(\nu)\} \left(\|x\|/\alpha\right)^{\nu} \Kcal_{\nu}\left(\alpha\|x\|\right)$ for $x\in \RR^d$, then (<ref>) can be written as $\theta K_{\alpha,\nu}(\bfs-\bft)$ for $\bfs,\bft\in [0,1]^d$. For the Gaussian processes $X$ and $Y$ in Model (<ref>), we write $X\sim \gp(0,\theta K_{\alpha,\nu})$ and $Y\sim \gp( \beta^\T  \ff, \theta K_{\alpha,\nu} + \tau \delta_0)$, respectively, where $\delta_0$ denotes the Dirac delta function at zero.

We study the Bayesian inference of Model (<ref>) (15, 11), where the common practice is to assign prior distributions on the regression coefficients $\beta$, the nugget parameter $\tau$, and the parameters in the covariance function of $X(\cdot)$ such as $\theta$ and $\alpha$. Bayesian posterior inference, including the prediction of $Y(s^*)$ at a new location $s^*\in [0,1]^d$, is fully based on the posterior distribution of these model parameters. In application, these parameters are randomly drawn from their posterior using sampling algorithms such as Markov chain Monte Carlo. Therefore, the behavior of their posterior distributions will heavily affect the Bayesian posterior predictive performance. However, despite the routine practice of using them in spatial statistics, there is a severe lack of theoretical understanding of Bayesian parameter estimation in the general spatial Gaussian process regression model with nugget (<ref>), especially on the asymptotic behavior of the posterior distribution of model parameters when the number of observations becomes large. Such theory becomes increasingly important given the advance in high-resolution remote sensing technology, and the research on Gaussian process with massive spatial datasets is prevalent in Bayesian spatial statistics, such as 6, 36, 10, 17, 21, 34, 14, etc.

We consider the Bayesian fixed-domain asymptotics (or infill asymptotics) framework (42, 51). For the first time in the literature, we derive the fixed-domain Bayesian posterior contraction rates of model parameters in the general spatial Gaussian process regression model (<ref>) with both the regression term $\ff(\bfs)^\T \beta$ and the measurement error (with the nugget $\tau$). In the fixed-domain asymptotics regime, the spatial domain $[0,1]^d$, as we have assumed, always remains fixed and bounded. This implies that as $n$ increases, the sampling locations in $S_n$ become increasingly dense in the domain $[0,1]^d$, leading to increasingly stronger dependence between adjacent observations in the data $Y_n$. The fixed-domain asymptotics has several advantages over other regimes such as increasing-domain asymptotics (32). First, a fixed domain matches up with the reality in many spatial applications. For instance, the advances in remote sensing technology enable routine collection of spatial data in larger volume and higher resolution from a given region (44). Second, the stationarity assumption on the Gaussian process $X(\cdot)$ is more likely to hold on a fixed domain rather than an expanding domain. Therefore, the fixed-domain asymptotics regime is more suitable for interpolation of spatial processes (Section 3.3 of 42). Third, fixed-domain asymptotics has better parameter estimation performance than the increasing-domain asymptotics (52).

On the other hand, fixed-domain asymptotics also poses several significant theoretical challenges. We take the isotropic Matérn covariance function for example. The first challenge comes from the lack of identification for the covariance parameters $(\sigma^2,\alpha)$ in (<ref>) when the domain dimension $d=1,2,3$, due to the increasingly stronger dependence in the data $Y_n$. We briefly review some basic results in equivalence of Gaussian processes. If two Gaussian processes are defined on a fixed domain, then by the result in 20 and Chapter 4 of 42, they must be either equivalent or orthogonal to each other, where the equivalence means that their Gaussian measures are mutually absolutely continuous. In particular, 51 has shown the well-known result that when the spatial domain has the dimension $d=1,2,3$, for two isotropic Matérn with parameters $(\sigma^2_1,\alpha_1)$ and $(\sigma^2_2,\alpha_2)$ and the same smoothness $\nu$, their Gaussian measures are equivalent if $\sigma_1^2\alpha_1^{2\nu}=\sigma_2^2\alpha_2^{2\nu}$, and they are orthogonal otherwise. This equivalence result for $d=1,2,3$ is a consequence of the integral test for the spectral densities of two covariance functions; see for example, Theorem A.1 of 43. On the other hand, when $d\geq 5$, 1 has proposed consistent moment estimators for both $\sigma^2$ and $\alpha$ with a given $\nu$ under fixed-domain asymptotics. Recently 8 have shown that for any $d\geq 4$, two Gaussian measures are equivalent if and only if they have the same parameters $(\sigma^2,\alpha,\nu)$, though finding the consistent estimators of $(\sigma^2,\alpha)$ for $d=4$ under fixed-domain asymptotics requires further study.

A direct consequence of the theory on equivalence of Gaussian processes is that for the domain dimension $d=1,2,3$, which is of primary interest to spatial statistics, only the microergodic parameter $\theta=\sigma^2\alpha^{2\nu}$ and continuous functions of $\theta$ can be consistently estimated from the data under fixed-domain asymptotics. However, the individual variance parameter $\sigma^2$ and range parameter $\alpha$ have no consistent estimators when $d=1,2,3$, and therefore they cannot have Bayesian posterior consistency. In general, the microergodic parameter is defined to be the parameter that uniquely determines the equivalence class of Gaussian processes; see Section 6.2 of 42 for a detailed explanation on the microergodic parameter. Besides Matérn, the exact form of microergodic parameter $\theta$ for $d=1,2,3$ has been identified in several other families, such as generalized Wendland (7) and confluent hypergeometric covariance functions (31); see Examples <ref> and <ref> in Section <ref>. The microergodic parameter has been further studied for covariance functions on non-Euclidean domains, such as spheres (2) and Riemannian manifolds (26), where 26 have discovered a similar phenomenon for the Matérn covariogram when $d=1,2,3$. Consequently, statistical inference for fixed-domain asymptotics when $d=1,2,3$ differs completely from that of regular parametric models with independent or weakly dependent data. For example, instead of increasing with the sample size, the Fisher information for the non-microergodic parameters $(\sigma^2,\alpha)$ has a finite limit in the case $d=1,2,3$ under fixed-domain asymptotics (52, 27).

The second main challenge comes from the impact of the regression terms $\ff(\bfs)^\T \beta$ and the nugget parameter $\tau$ on the asymptotics of other parameters, in particular the microergodic parameter $\theta$ when $d=1,2,3$. In real spatial applications, typically one assumes that $Y_n$ from Model (<ref>) contains some regression terms and measurement error, though most existing frequentist fixed-domain asymptotics works have only considered estimation in a Gaussian process without regression terms and measurement error (50, 51, 12, 47, 23, 7, 31, etc.) In these works, the maximum likelihood estimator (MLE) of the microergodic parameter $\theta$ is asymptotically normal with a parametric $n^{-1/2}$ convergence rate. The cross-validation estimator (5) and the composite likelihood estimator (3) have the same parametric rate for Matérn with $\nu=1/2$ and $d=1$. This is no longer true for the microergodic parameter $\theta$ in Model (<ref>). Although the nugget $\tau$ is also identifiable in Model (<ref>) (42, 46), it will dramatically change the convergence rate for the estimators of $\theta$. For example, 9 have studied a special case of Model (<ref>) without regression terms, with $d=1$, $\nu=1/2$, and $S_n$ being the equispaced grid in $[0,1]$. The convergence rate for the MLE of $\theta$ slows down from $n^{-1/2}$ in the model without $\tau$ to $n^{-1/4}$ in the model with $\tau$. Under the assumptions of equispaced grid and exact decaying rate of eigenvalues of the Matérn covariance matrix, 46 have shown that the asymptotic normality for the MLE of $\theta$, whose convergence rate is $n^{-1/(2+4\nu/d)}$. This deterioration is mainly caused by the convolution of $X(\cdot)$ with the Gaussian noise $\varepsilon(\cdot)$ in Model (<ref>). As a consequence, the Bayesian fixed-domain asymptotic theory for the model with nugget (<ref>) becomes much more challenging than the Bayesian theory for the model without nugget (24).

The main contribution of this work is to devise a general theoretical framework for studying the Bayesian fixed-domain posterior contraction rates of the microergodic parameter $\theta$ and the nugget parameter $\tau$ in Model (<ref>) where $X(\cdot)$ has the covariance function $\theta K_{\alpha,\nu}(\bfs-\bft)$, which covers a wide range of stationary covariance functions in spatial statistics; see our Examples <ref>-<ref> in Section <ref>. In Bayesian asymptotic theory, the posterior contraction rate describes how fast the posterior distribution of certain parameters contracts to their true values. General Bayesian theory of posterior consistency and contraction rates has been extensively studied; see Chapters 6–9 in 13 for a thorough treatment. The main idea in establishing such Bayesian asymptotic results is to adapt the posterior consistency theorem from 37. That is, for models with independent data, one only needs to verify two conditions: (a) the prior distribution contains a Kullback-Leibler support of the target parameter or nonparametric function, and (b) the existence of an exponentially consistent test. Our theoretical framework resembles the Schwartz's consistency framework, but with two key adaptations to the spatial Gaussian process regression model (<ref>) under fixed-domain asymptotics.

First, we establish a new polynomial lower bound for the denominator of the posterior density given the strongly dependent data $Y_n$ under fixed-domain asymptotics, based on the spectral analysis of covariance function. This is vastly different from the lower bound of denominator for independent data where the lower bound is established by the law of large numbers in the original Schwartz's theorem. Based on this new evidence lower bound, we establish a general posterior contraction rate theorem for the microergodic parameter and the nugget parameter in Model (<ref>) under fixed-domain asymptotics.

Second, we propose new exponentially consistent tests for the microergodic parameter $\theta$ and the nugget parameter $\tau$ using higher-order quadratic variation estimators. We illustrate this using the isotropic Matérn covariance function in Example <ref>. Higher-order quadratic variation has been successfully applied to parameter estimation in frequentist fixed-domain asymptotics for Gaussian process models without nugget (28, 30). We further generalize this powerful tool to Model (<ref>) with nugget, similar to the recent work 29.

Our new techniques allow us to bypass the difficult problem of expanding the likelihood function for $(\theta,\tau)$ in the presence of nonidentifiable parameters such as the range parameter $\alpha$. For a fixed value of range parameter $\alpha$, the MLE of $(\theta,\tau)$ in Model (<ref>) without the regression terms $\ff(\bfs)^\T \beta$ satisfies the frequentist asymptotic normality (Theorem 5 of 46). However, these MLEs do not have any closed-form expression or even first-order approximation except for a few special cases (such as $\nu=1/2$, $d=1$, $S_n$ being an equispaced grid in 9), which makes it technically highly challenging to construct exponentially consistent tests based on these MLEs. Our higher-order quadratic variation estimators are constructive, explicit, and completely circumvent the difficulty caused by the complicated likelihood function. They can be directly applied to Model (<ref>) with the regression terms $\ff(\bfs)^\T \beta$. More importantly, they satisfy concentration inequalities with exponentially small tails which work uniformly well for the range parameter $\alpha$ in a wide expanding interval and for a large class of stratified sampling designs without requiring the sampling points $S_n$ to be on an equispaced grid. Details will be discussed in Section <ref>.

The rest of the paper is organized as follows. Section <ref> presents a new general framework for deriving posterior contraction rates for model parameters under fixed-domain asymptotics. Section <ref> presents the details of our high-order quadratic variation estimators and the explicit posterior contraction rates for the isotropic Matérn covariance function. Section <ref> contains simulation studies and a real data example that verify our theory. Section <ref> includes some discussion. The technical proofs of all theorems and propositions as well as additional simulation results can be found in the Supplementary Material.

§ GENERAL THEORY FOR FIXED-DOMAIN POSTERIOR
CONTRACTION RATES

We first provide some examples of the assumed covariance function $\theta K_{\alpha,\nu}(\bfs-\bft)$ for $\bfs,\bft\in [0,1]^d$. This format of covariance function is very general and includes several classes of covariance functions used in spatial statistics.

The isotropic Matérn covariance function defined in (<ref>) takes the form $\theta K_{\alpha,\nu}$ with $\theta=\sigma^2\alpha^{2\nu}$ and $K_{\alpha,\nu}(x)=\{2^{1-\nu}/\Gamma(\nu)\} \left(\|x\|/\alpha\right)^{\nu} \Kcal_{\nu}\left(\alpha\|x\|\right)$ for $x\in \RR^d$.

The tapered isotropic Matérn covariance function (22) is $\theta \tilde K_{\alpha,\nu}(\bfs-\bft)= \theta K_{\alpha,\nu}(\bfs-\bft) \times K_{\tap}(\bfs-\bft)$ for any $\bfs,\bft\in [0,1]^d$, where $\theta K_{\alpha,\nu}(\bfs-\bft)$ is the isotropic Matérn covariance function in Example <ref>, and $K_{\tap}: \RR^d \to \RR$ is a positive definite function, whose spectral density $f_{\tap}(w) = (2\pi)^{-d}\int_{\RR^d} \exp\left(-\imath w^\T  x \right) K_{\tap}(x) \ud x$ satisfies $0< f_{\tap}(w)\leq C_{\tap} (1+\|w\|^2)^{-(\nu+d/2+\eta_{\tap})}$ for some constants $C_{\tap}>0,\eta_{\tap}>0$ and all $w\in \RR^d$, where $\imath^2=-1$.

The isotropic generalized Wendland covariance function (7) is $\theta K_{\alpha,\nu}(\bfs-\bft)$ for any $\bfs,\bft\in [0,1]^d$, where
$\theta=\sigma^2 \alpha^{2\nu}$,

\begin{align} \label{eq:GWCov}
& K_{\alpha,\nu}(x) = \frac{\Ical(0\leq \|x\|<1/\alpha)}{\alpha^{2\nu} B(2\nu,\mu)}\int_{\alpha\|x\|}^1 \big(t^2-\alpha^2\|x\|^2 \big)^{\nu-1/2} (1-t)^{\mu -1} \ud t, \quad \text{ for } x\in \RR^d,
\end{align}

and $\sigma^2>0$, $\alpha>0$, $\nu \geq 1/2$, $\mu > \nu + d$, $B(\cdot,\cdot)$ is the Beta function, $\Ical(\cdot)$ is the indicator function. We assume that $\mu$ is known and suppress the dependence of $K_{\alpha,\nu}(\cdot)$ on $\mu$.

The isotropic confluent hypergeometric covariance function (31) is $\theta K_{\alpha,\nu}(\bfs-\bft)$ for any $\bfs,\bft\in [0,1]^d$, where $\theta=\sigma^2 \alpha^{2\nu} \Gamma(\nu+\mu)/\Gamma(\mu)$,

\begin{align} \label{eq:CHCov}
& K_{\alpha,\nu}(x) = \frac{1}{\Gamma(\nu)} \int_0^{\infty} t^{\nu-1} \big(\alpha^2t+1\big)^{-(\nu+\mu)} \exp\left(-\nu \|x\|^2/t \right) \ud t , \quad \text{ for } x\in \RR^d,
\end{align}

where $\sigma^2>0$, $\alpha>0$, $\nu>0$ and $\mu>0$. We assume that $\mu$ is known and suppress the dependence of $K_{\alpha,\nu}(\cdot)$ on $\mu$, though $\mu$ can also be estimated as in 31.

Among these examples, the Matérn covariance function is of central importance since the other three can be viewed as its generalizations. The tapered Matérn and the generalized Wendland covariance functions have compact supports and hence computationally more efficient, while the confluent hypergeometric is suitable for modeling polynomially decaying spatial dependence.

Throughout the paper, we assume that the smoothness parameter $\nu>0$ is fixed and known. Estimation of $\nu$ is an important and technically challenging problem for spatial statistics, with some recent advances in the frequentist fixed-domain asymptotics (28, 30, 29); see our detailed discussion in Section <ref>. Our main goal is to perform Bayesian inference on the parameters $(\theta,\alpha,\tau,\beta)$ in Model (<ref>) based on the observed data $(Y_n, F_n)$. We assume that the true parameters are $(\theta_0,\alpha_0,\tau_0,\beta_0)$. We use $\PP_{(\theta,\alpha,\tau,\beta)}$ to denote the probability measure of $\gp(\beta^\T  \ff,\sigma_0^2K_{\alpha_0,\nu}+\tau \delta_{0})$, and hence $\PP_{(\theta_0,\alpha_0,\tau_0,\beta_0)}$ is the true probability measure of the observed process $Y(\cdot)$.

The log-likelihood function based on the data $(Y_n, F_n )$ is

\begin{align} \label{eq:loglik}
\Lcal_n(\theta,\alpha,\tau,\beta) &= - \frac{1}{2} (Y_n-  F_n  \beta)^\T  \left\{\theta K_{\alpha,\nu}(S_n) + \tau I_n \right\}^{-1} (Y_n-  F_n  \beta)  \nonumber \\
&\quad  - \frac{1}{2}\log
\dett\left\{\theta K_{\alpha,\nu}(S_n) + \tau I_n \right\},
\end{align}

where $K_{\alpha,\nu}(S_n)$ is the $n\times n$ covariance matrix whose $(i,j)$-entry is $K_{\alpha,\nu}(\bfs_i-\bfs_j)$, $I_n$ is the $n\times n$ identity matrix, $\dett(A)$ is the determinant of a matrix $A$, and $cA$ denotes the matrix of $A$ with all entries multiplied by the number $c$.

For Bayesian inference, we impose a prior distribution with the density $\pi(\theta,\alpha,\tau,\beta)$ on the parameters. Let $\RR_+=(0,+\infty)$ and $\ZZ_+$ be the set of all positive integers. Then the posterior density of $(\theta,\alpha,\tau,\beta)$ is

\begin{align}\label{jointpost1}
\pi(\theta,\alpha,\tau,\beta~|~ Y_n, F_n ) &= \frac{\exp\left\{ \Lcal_n(\theta,\alpha,\tau,\beta) \right\} \pi(\theta,\alpha,\tau,\beta)}{\int_{\RR_+^{3}\times \RR^p} \exp\left\{\Lcal_n(\theta',\alpha',\tau',\beta') \right\} \pi(\theta',\alpha',\tau',\beta')\ud \theta' \ud \tau' \ud \alpha' \ud \beta' } .
\end{align}

We use $\Pi(\cdot)$ and $\Pi(\cdot \mid Y_n, F_n )$ to denote the prior and posterior measures, respectively.

Our general framework requires a few assumptions on the model setup.

The data $(Y_n, F_n)$ are observed on the distinct sampling points $S_n=\{\bfs_i:i=1,\ldots,n\}$ in a fixed domain $[0,1]^d$ following the model (<ref>). The sequence of $S_n$ is getting dense in $[0,1]^d$ as $n\to\infty$, in the sense that $\sup_{\bfs^*\in [0,1]^d} \min_{1\leq i\leq n} \|\bfs^*-\bfs_i\| \to 0$ as $n\to\infty$. There exists a constant $C_{\ff}>0$ such that $|\ff_l(\bfs)|\leq C_{\ff}$ for all $l=1,\ldots,p$ and all $\bfs\in [0,1]^d$.

Let $f_{\theta,\alpha,\nu}(w)= (2\pi)^{-d}\int_{\RR^d} \exp(-\imath w^\T  x) \theta K_{\alpha,\nu}(x) \ud x$ be the spectral density of the covariance function $\theta K_{\alpha,\nu}$, where $\imath^2 = -1$ and $\alpha \in \RR_+ $. Then

(i) There exist constants $L>0,r_0\in (0,1/2),\kappa>0$ that may depend on $\theta_0,\alpha_0,\nu,d$, such that for all $\alpha$ that satisfies $|\alpha/\alpha_0 -1|\leq r_0$,

\begin{align*}
\sup_{w\in \RR^d} \left|f_{\theta_0,\alpha,\nu}(w)/f_{\theta_0,\alpha_0,\nu}(w)-1\right| \leq L |\alpha/\alpha_0-1|^{\kappa}.
\end{align*}

(ii) For any given $\theta>0$ and $w\in \RR^d$, $f_{\theta,\alpha,\nu}(w)$ is a non-increasing function in $\alpha$.

All the four covariance functions in Examples <ref>–<ref> satisfy Assumption <ref>.

Assumption <ref> imposes mild and general conditions on the spectral density $f_{\theta,\alpha,\nu}$. For example, the spectral density of the isotropic Matérn in Example <ref> is $f_{\theta,\nu}(w)=\Gamma(\nu+d/2)/\Gamma(\nu) \times \pi^{-d/2}\theta (\alpha^2+\|w\|^2)^{\nu+d/2}$, which is both Lipschitz continuous ($\kappa=1$) and decreasing in $\alpha$. The same can be proved for Examples <ref>–<ref>. Assumption <ref> makes our Bayesian fixed-domain asymptotic theory much more general and works far beyond the popular Matérn covariance function. The proof of Proposition <ref> is in the Supplementary Material.

The prior probability measure $\Pi(\cdot)$ does not depend on $n$ and has a proper density $\pi(\theta,\alpha,\tau,\beta)$ that is continuous in a neighborhood of $(\theta_0,\alpha_0,\tau_0,\beta_0)$, and satisfies $\pi(\theta_0,\alpha_0,\tau_0,\beta_0)>0$. The marginal prior of $\tau$ satisfies $\int_0^{\infty} \tau^{-n/2}\pi(\tau) \ud \tau <\infty$ for any $n\in \ZZ_+$.

Assumption <ref> requires that the prior distribution has some positive probability mass around the true parameter $(\theta_0,\alpha_0,\tau_0,\beta_0)$. This is a mild condition and necessary for showing posterior consistency. The finite integral condition on $\pi(\tau)$ is to ensure the propriety of the posterior density defined in (<ref>), and is satisfied by the commonly used inverse gamma prior on $\tau$; see Proposition <ref> in Section <ref> below.

To formulate the posterior contraction, for any $\epsilon_1>0,\epsilon_2>0$, we define the neighborhood

$$\Bcal_0(\epsilon_1,\epsilon_2)=\{(\theta,\alpha,\tau,\beta): |\theta/\theta_0-1| < \epsilon_1 , |\tau/\tau_0-1|< \epsilon_2, \alpha \in \RR_+, \beta \in \RR^p \} .$$

This is an open neighborhood for the microergodic parameter $\theta$ and the nugget parameter $\tau$, but with the range parameter $\alpha$ and the regression coefficient $\beta$ unconstrained. We will focus on showing that the posterior probability of this neighborhood converges to one and characterize the dependence of $\epsilon_1,\epsilon_2$ on the sample size $n$ under the fixed-domain asymptotics regime.

Our posterior contraction results will exclude the range parameter $\alpha$ and the regression coefficient $\beta$. This is because they are in general inconsistent under fixed-domain asymptotics following the general frequentist theory of 51 when the domain dimension $d=1,2,3$. This will be proved rigorously in Section S3 of the Supplementary Material. In general, their inconsistency is not likely to affect the asymptotic prediction performance of Model (<ref>) as long as $(\theta,\tau)$ are consistent (40).

In the original Schwartz's consistency theorem, one important step is to show that the denominator of the posterior density has a lower bound. In the regular parametric models with independent observations, this is achieved by applying concentration inequalities to the log-likelihood ratio together with Jensen's inequality and Fubini's theorem; see for example, the proof of Theorem 6.17 and Lemma 8.10 in 13. This lower bound of denominator is also known as the evidence lower bound, which plays a key role in determining the posterior contraction rates.

Our first result is the following theorem that characterizes the evidence lower bound based on the strongly dependent data $Y_n$ observed from Model (<ref>) under fixed-domain asymptotics.

Suppose that Assumptions <ref>, <ref> and <ref> hold. There exists a constant $D>0$ and a large integer $N_1$ that only depend on $\nu,d,p$, $\theta_0,\tau_0,\alpha_0$, $C_{\ff},L,r_0,\kappa$ and $\pi(\theta_0,\tau_0,\alpha_0,\beta_0)$, such that for all $n>N_1$, with probability at least $1-3\exp(-\log^2 n)$,

\begin{align}
& \int \exp\left\{\Lcal_n(\theta,\alpha,\tau,\beta) - \Lcal_n(\theta_0,\alpha_0,\tau_0,\beta_0) \right\} \cdot \pi(\theta,\alpha,\tau,\beta)~ \ud \theta \ud \tau \ud \alpha \ud \beta \geq  Dn^{-(3p+2+2/\kappa)}.  \nonumber
\end{align}

Theorem <ref> shows a polynomially decaying evidence lower bound for the posterior distribution, which can be used for showing posterior contraction together with some exponentially consistent tests. The major difference between the proof of Theorem <ref> and that of regular parametric models with independent data is that our data $Y_n$ are strongly dependent and our likelihood ratio depends heavily on the covariance matrix $\theta K_{\alpha,\nu}(S_n)+\tau I_n$, where $\theta K_{\alpha,\nu}(S_n)$ becomes increasingly singular as $S_n$ becomes denser in the domain $[0,1]^d$. Thus the standard techniques for showing evidence lower bound with independent data such as Lemma 8.10 of 13 cannot be directly applied here. Fortunately, the structure of this covariance matrix can be related to the spectral density $f_{\theta,\alpha,\nu}$, such that under Assumption <ref>, the likelihood ratio can be lower bounded when $(\theta,\alpha,\tau,\beta)$ is in a shrinking neighborhood of $(\theta_0,\alpha_0,\tau_0,\beta_0)$.

The next two assumptions are on the existence of consistent estimators for $(\theta,\tau)$ and the prior distribution. We define the function $\varphi(x) = \min(x^2,x)$ for $x>0$, and $\Ecal^c$ to be the complement of a generic set $\Ecal$.

There exist estimators $\widehat\theta_n$ and $\widehat\tau_n$, positive constants $b_1,b_2,c_1,c_2$, and a parameter set $\Ecal_n\subseteq \RR_+^{3} \times \RR^p$, such that for any $\epsilon_1,\epsilon_2 \in (0,1/2)$, for all $n>n_0$ where $n_0$ depends on $\epsilon_1,\epsilon_2$,

\begin{align}
\PP_{(\theta_0,\alpha_0,\tau_0,\beta_0)} \left(|\widehat\theta_n/\theta_0-1|\geq \epsilon_1/2 \right) & \leq \exp\left\{- c_1 \varphi\left( n^{b_1} \epsilon_1 \right)\right\}, \label{eq:test.theta1} \\
\PP_{(\theta_0,\alpha_0,\tau_0,\beta_0)} \left( |\widehat\tau_n/\tau_0-1|\geq \epsilon_2/2 \right) & \leq \exp\left\{ - c_2 \varphi\left( n^{b_2} \epsilon_2 \right) \right\}, \label{eq:test.tau1} \\
\sup_{\Bcal_0(\epsilon_1,\epsilon_2)^c \cap \Ecal_n} \PP_{(\theta,\alpha,\tau,\beta)} \left( |\widehat\theta_n/\theta_0-1|\leq \epsilon_1/2 \right) & \leq \exp\left\{ - c_1 \varphi\left(n^{b_1} \epsilon_1\right)\right\}, \label{eq:test.theta2} \\
\sup_{\Bcal_0(\epsilon_1,\epsilon_2)^c \cap \Ecal_n} \PP_{(\theta,\alpha,\tau,\beta)} \left(  |\widehat\tau_n/\tau_0-1|\leq\epsilon_2/2 \right) & \leq \exp\left\{- c_2 \varphi\left( n^{b_2} \epsilon_2\right)\right\}, \label{eq:test.tau2}
\end{align}

where the supremum is taken over $(\theta,\alpha,\tau,\beta)$.

For the parameter set $\Ecal_n$ in Assumption <ref>, the prior satisfies $\Pi(\Ecal_n^c ) \leq n^{-(3p+4+2/\kappa)}$ for all sufficiently large $n$, where $\kappa$ is given in Assumption <ref>.

Assumption <ref> assumes the existence of consistent estimators $\widehat\theta_n$ and $\widehat\tau_n$ with exponential tail probabilities for large values of $\epsilon_1,\epsilon_2$ and sub-Gaussian tails for small values of $\epsilon_1,\epsilon_2$, which are reasonable given the normality assumption on both $X(\cdot)$ and $\varepsilon(\cdot)$ and can usually be derived from the Hanson-Wright inequality (35). The four inequalities are used for constructing uniformly consistent tests in the Schwartz's theorem; see for example, Proposition 6.22 in 13. In particular, the set $\Ecal_n$ is a sieve parameter space, which typically expands as $n$ increases and eventually covers the whole parameter space. Furthermore, we require in Assumption <ref> that the prior probability of the sieve space $\Ecal_n$ converges to one faster than the evidence lower bound in Theorem <ref>. This puts a mild restriction on the tail behavior of the prior distribution, which is generally weaker than assuming a prior tail probability exponentially small in $n$ as used in the Bayesian nonparametrics literature; see for example, Condition (iii) in Theorems 8.9 and 8.11 in 13. In Section <ref>, we will show that the higher-order quadratic variation estimators of $\theta$ and $\tau$ satisfy Assumption <ref>, and will discuss the concrete prior tail conditions that satisfy Assumption <ref>.

The following theorem provides the general posterior contraction rates for $\theta$ and $\tau$.

Suppose that Assumptions <ref>, <ref>, <ref>, <ref>, and <ref> hold. Then the posterior distribution $\Pi(\cdot \mid Y_n, F_n )$ with the density given in (<ref>) exists and contracts to the true parameters $(\theta_0,\tau_0)$ at the rate $(n^{-b_1} \log n, n^{-b_2}\log n)$ with $b_1$ and $b_2$ given in Assumption <ref>, in the sense that for any positive sequence $M_n\to \infty$ as $n\to\infty$,

\begin{align}
\Pi \left( |\theta/\theta_0 - 1| < M_n n^{-b_1} \log n \text{ and } |\tau/\tau_0 - 1|< M_n n^{-b_2}\log n \mid Y_n, F_n  \right) \to 1,  \nonumber
\end{align}

almost surely $\PP_{(\theta_0,\alpha_0,\tau_0,\beta_0)}$ as $n\to\infty$.

Theorem <ref> is an adapted version of the Schwartz's theorem and is established based on our new evidence lower bound in Theorem <ref>, as well as the uniformly consistent tests constructed from the estimators $\widehat\theta_n$ and $\widehat\tau_n$ in Assumption <ref>. Theorem <ref> essentially translates the complicated problem of finding Bayesian posterior contraction rates under fixed-domain asymptotics to the more direct problem of finding the frequentist consistent estimators $\widehat\theta_n$ and $\widehat\tau_n$ that satisfy Assumption <ref>. The provable rates $n^{-b_1} \log n$ for the microergodic parameter $\theta$ and $n^{-b_2} \log n$ for the nugget parameter $\tau$ are completely determined by how efficient the frequentist estimators of $\theta$ and $\tau$ are inside the four inequalities of Assumption <ref>. We find $b_1$ and $b_2$ explicitly in Section <ref>.

§ HIGHER-ORDER QUADRATIC VARIATION ESTIMATORS FOR
ISOTROPIC MATÉRN

For the isotropic Matérn covariance function in Example <ref>, we verify Assumptions <ref> and <ref> by proposing the higher-order quadratic variation estimators of the microergodic parameter $\theta$ and the nugget parameter $\tau$. We explicitly find the rates $n^{-b_1} \log n$ and $n^{-b_2} \log n$ in Theorem <ref>, and thus provide the posterior contraction rates for $\theta$ and $\tau$ in Model (<ref>) under fixed-domain asymptotics. To focus on the main idea and also due to the space limit, we will only show the rates for isotropic Matérn which can be compared directly with the existing frequentist fixed-domain asymptotic results. Our techniques can also be extended to the three covariance functions in Examples <ref>, <ref> and <ref> with additional technical adjustment.

Higher-order quadratic variation has been adopted in recent theoretical works of spatial Gaussian processes (28, 30, 29). The basic idea is to use a carefully designed series of constants ($c_{\bfi,d,\ell}^{(k_1,\ldots,k_d)}$ and $c_{d,\ell}^{(k_1,\ldots,k_d)}$ in Lemma <ref>), such that finite differencing the observations $Y_n$ weighted by these constants can approximately solve for both the microergodic and nugget parameters. This method has at least two major technical advantages. First, as a frequentist estimator, it is simple to implement in practice and does not involve any numerical optimization such as the maximum likelihood estimator. Second, it does not require the observations $Y_n$ to be made strictly on an equispaced grid and hence is applicable to a much wider range of sampling designs. The previous works 28 and 30 have considered such estimators for the model without measurement error and nugget parameter $\tau$, while our version of higher-order quadratic variation estimators are similar to those in 29 which work for the general spatial model (<ref>) with measurement error and the nugget $\tau$.

§.§ Construction of Estimators

We consider the stratified sampling design, in which the domain $[0,1]^d$ is divided into cells with the same size and observations are made inside each cell, but not necessarily on the grid points. Without loss of generality, we let $n=m^d$ for an integer $m$ since we only study the asymptotics in $n$. We rewrite the set $S_n=\{\bfs_1,\ldots,\bfs_n\}$ using the ordering in each coordinate:

\begin{align} \label{eq:Sn}
S_n=\left\{\bfs(\bfi)=(s_1(\bfi),\ldots,s_d(\bfi))^\T :~ \bfi=(i_1,\ldots,i_d), 1\leq i_1,\ldots,i_d\leq m\right\},
\end{align}

and for $k=1,\ldots,d$, $s_k(\bfi)$ satisfies

\begin{align} \label{eq:si}
s_k(\bfi) &= \frac{i_k-1}{m} + \frac{\delta_{\bfi;k}}{m} ,
\end{align}

where $0\leq \delta_{\bfi;k}<1$ is arbitrary for all $k=1,\ldots,d$. We make the following assumption on the stratified sampling design.

Suppose that the sampling points $S_n$ satisfies (<ref>) and (<ref>), where all $\delta_{\bfi;k}$'s can take arbitrary values in $[0,1)$. Then for any $1\leq i_1,\ldots,i_d\leq m$,

\begin{align*}
& \bfs(\bfi) \in \left[\frac{i_1-1}{m}, \frac{i_1}{m} \right) \times \ldots \times \left[\frac{i_d-1}{m}, \frac{i_d}{m} \right) .
\end{align*}

This stratified sampling design is much more relaxed than assuming sampling locations strictly on equispaced grid points. We allow the perturbations $\delta_{\bfi;k}$ to be arbitrary in $[0,1)$, so there is little restriction other than that the sampling locations need to be roughly evenly distributed across the domain.

For the regression functions $\ff_1,\ldots,\ff_p$, there exists a constant $C_{\ff}'>0$ such that their partial derivatives satisfy $\left|\sD^{\sj}~\ff_l(\bfs)\right| \leq C_{\ff}'$ for all $l=1,\ldots,p$, all $\bfs\in [0,1]^d$, and all index vector $\sj=(j_1,\ldots,j_d)\in \NN^d$ that satisfies $j_1+\ldots +j_d\leq \lceil \nu+d/2\rceil$, where $\lceil x\rceil$ is the smallest integer greater than or equal to $x$, and $\sD^{\sj}$ is the partial differentiation operator of order $\sj$.

Assumption <ref> requires the regression functions $\ff_1,\ldots,\ff_p$ to have bounded derivatives up to the order of $\lceil \nu+d/2\rceil$. Because the sample path $X(\cdot)$ of the Matérn covariance function (<ref>) are only mean square differentiable up to the order of $\nu$, we essentially assume that the functions $\ff_1,\ldots,\ff_p$ are smoother than the sample path $X(\cdot)$. Furthermore, the order $\lceil \nu+d/2\rceil$ implies that the functions $\ff_1,\ldots,\ff_p$ all lie in the reproducing kernel Hilbert space of $\theta K_{\alpha,\nu}$, which is known to be norm equivalent to the Sobolev space of order $\nu+d/2$ (Corollary 10.48 of 49), and therefore satisfy Assumption <ref>.

Let $\lfloor x\rfloor$ be the largest integer less than or equal to $x$. Define $\omega_m=\lfloor m^{\gamma} \rfloor$ for a positive constant $\max\big\{1-d/(4\nu),0\big\} < \gamma < 1$ and assume that $\omega_m$ is an even integer without loss of generality. The higher-order quadratic variation method relies on the following sequence of $m$-dependent (and so $n$-dependent) constants $c_{\bfi,d,\ell}^{(k_1,\ldots,k_d)}$ and $c_{d,\ell}^{(k_1,\ldots,k_d)}$ described in 30.

(Corollary 1 and Lemma 2 of 30) Let $d,\ell\in \ZZ_+$. Let $\bfi=(i_1,\ldots,i_d)^\T $ where $1\leq i_1,\ldots,i_d\leq m - \ell \omega_m$. Then there exists a sequence of constants

$$\left\{c_{\bfi, d, \ell}^{(k_1,\ldots,k_d)}: 0\leq k_1,\ldots,k_d\leq \ell, \text{and } \bfi= (i_1,\ldots,i_d), 1\leq i_1,\ldots,i_d\leq m \right\}, $$

such that for any $\bfs(\bfi) = (s_1(\bfi),\ldots,s_d(\bfi))^\T  \in [0,1]^d$, for all integers $l_1,\ldots,l_d$ satisfying $0\leq l_1,\ldots,l_{d-1}\leq \ell$, $0\leq l_d\leq \ell-1$, and $0\leq l_1+\ldots + l_d\leq \ell$,

\begin{align*}
& \textstyle \sum_{0\leq k_1,\ldots,k_d\leq \ell} c_{\bfi,d,\ell}^{(k_1,\ldots,k_d)} \cdot s_d(i_1+k_1\omega_m,\ldots,i_d+k_d\omega_m)^{\ell} = \ell! \left(\omega_m/m\right)^{\ell}, \\
\text{and } & \textstyle \sum_{0\leq k_1,\ldots,k_d\leq \ell} c_{\bfi,d,\ell}^{(k_1,\ldots,k_d)} \cdot \textstyle \prod_{j=1}^d s_j(i_1+k_1\omega_m,\ldots,i_d+k_d\omega_m)^{l_j} = 0 .
\end{align*}

Furthermore, there exists a sequence of $m$-independent constants $c_{d,\ell}^{(k_1,\ldots,k_d)}$, such that

\begin{align*}
& c_{\bfi,d,\ell}^{(k_1,\ldots,k_d)} = c_{d,\ell}^{(k_1,\ldots,k_d)} + O(\omega_m^{-1}), \quad \text{for all } 0\leq k_1,\ldots,k_d\leq \ell,
\end{align*}

as $m\to\infty$, where the term $O(\omega_m^{-1})$ is uniform over all $1\leq i_1,\ldots,i_d\leq m- \ell\omega_m$ and any $\delta_{\bfi;k}\in [0,1)$ for $k=1,\ldots,d$.

For $d, \ell \in \ZZ_+$ and $\bfi=(i_1,\ldots, i_d)^\T$ such that $1\leq i_1, \ldots, i_d\leq m-\ell \omega_m$, we
define the following differencing operator

\begin{align*}
& \nabla_{ d, \ell} Y \left(\bfs(\bfi) \right)
=  \textstyle \sum_{0\leq k_1, \ldots, k_d\leq \ell}  c_{\bfi, d, \ell}^{(k_1, \ldots, k_d)}  Y \left(\bfs (i_1+k_1 \omega_m , \ldots, i_d+ k_d \omega_m ) \right),
\end{align*}

and we define the $\ell$th-order quadratic variation $V_{u,d, \ell} = V_{u,d, \ell}(m)$ by

\begin{align} \label{eq:3.1}
& V_{u, d, \ell} = \textstyle \sum_{\bfi\in \Xi_{u,m}} \left\{ \nabla_{d, \ell} Y \left(\bfs(\bfi)\right) \right\}
\left \{ \nabla_{d, \ell} Y \left(\bfs(\bfi + u\bfe_1) \right) \right\},
\end{align}

where $\bfe_1 = (1,0,\ldots,0)^\T  \in \RR^d$, $u\in\{0,1\}$ and
$\Xi_{u,m}=\left\{\bfi:~ 1\leq i_1+ u,i_1, \ldots, i_d \leq m- 2 \ell \omega_m \right\}$. The cardinality of the set $\Xi_{u,m}$, denoted by $|\Xi_{u,m}|$, has the same order as $m^d$.

We set $\ell_{\star}=\lceil \nu+d/2 \rceil$. Our higher-order quadratic variation estimators for $\tau$ and $\theta$ are defined as follows:

\begin{align}
\widehat{\tau}_n &= V_{0,d,\ell_{\star}} / C_{V,0},  \qquad \widehat{\theta}_n = V_{1, d,\ell_{\star}} / g_{\ell_{\star},\nu} ,\label{eq:tau.theta.es} \\
\text{where }
C_{V,0} &= \textstyle\sum_{\bfi \in \Xi_{0,m}} \textstyle\sum_{0\leq k_1, \ldots, k_{d}\leq \ell_{\star}} \left\{c_{\bfi, d, \ell_{\star}}^{(k_1, \ldots, k_d)} \right\}^2, \quad g_{\ell_{\star},\nu} = \left(\omega_m/m\right)^{2 \nu} |\Xi_{1, m}| \xi^*_{\nu} H_{\ell_{\star},\nu},  \nonumber \\
H_{\ell_{\star},\nu} &= \textstyle\sum_{0\leq k_1,\ldots, k_{2d}\leq \ell_{\star}} c_{d,\ell_{\star}}^{(k_1,\ldots, k_d)} c_{d,\ell_{\star}}^{(k_{d+1},\ldots, k_{2d})}
G_\nu \left( \big\|  (k_1,\ldots, k_d)^\T  - (k_{d+1},\ldots, k_{2d} )^\T  \big\| \right), \nonumber \\
\text{and } G_\nu (t)& = t^{2 \nu}, ~~ \xi^*_{\nu} = -\pi/\{2^{2\nu}\Gamma(\nu+1)\Gamma(\nu)\sin(\nu\pi) \}  , \text{ if } \nu\notin \ZZ, \nonumber \\
G_\nu (t)& = t^{2 \nu} \log(t), ~~ \xi^*_{\nu} = (-1)^{\nu+1}/\{2^{2\nu-1}\nu!(\nu-1)!\} , \text{ if } \nu\in \ZZ_+ . \label{eq:CgH}
\end{align}

We also define the set $\Ecal_{n}$ in Assumption <ref> as

\begin{align} \label{eq:E1n}
\Ecal_{n}=\left\{ (\theta,\alpha,\tau,\beta)\in \RR_+^{3}\times \RR^p: \frac{\|\beta\|^2}{\theta}  \leq  n^{\rho_1}, \frac{\tau}{\theta} \in [n^{-\rho_{21}},n^{\rho_{22}}], \alpha\in [n^{-\rho_{31}},n^{\rho_{32}}] \right\} ,
\end{align}

where the constants $\rho_1,\rho_{21},\rho_{22},\rho_{31},\rho_{32}$ satisfy the following relations:

\begin{align} \label{eq:rho.relation1}
& 0< \rho_1 < 1-\gamma,  \quad
0<\rho_{21} < 2\nu(1-\gamma)/d, \nonumber \\
& 0<\rho_{22} < 1/2 - 2(1-\gamma)\nu/d, \quad \rho_{31}>0, \quad 0<\rho_{32} < (1-\gamma)/d.
\end{align}

In the next theorem, we show that the estimators of $\widehat{\tau}_n$ and $\widehat\theta_n$ in (<ref>) together with the set $\Ecal_{n}$ in (<ref>) satisfy Assumption <ref>.

Suppose that $X(\cdot)$ in Model (<ref>) has the isotropic Matérn covariance function in (<ref>). Suppose that Assumptions <ref>, <ref> and <ref> hold. Let $\omega_m=\lfloor m^{\gamma} \rfloor$ be an even number for a constant $\max\big\{1-d/(4\nu),0\big\} < \gamma < 1$. Let $\Ecal_n$ in Assumption <ref> be defined in (<ref>) with $\rho_1,\rho_{21},\rho_{22},\rho_{31},\rho_{32}$ satisfying (<ref>). Then uniform over all $\delta_{\bfi;k}\in [0,1)$ with $\bfi=(i_1,\ldots,i_d)^\T $ for all $1\leq i_1,\ldots,i_d\leq m$ and $1\leq k \leq d$, for any $d\in \ZZ_+$,

(i) $\widehat{\theta}_n$ in (<ref>) satisfies (<ref>) and (<ref>) in Assumption <ref>, where $b_1$ is

\begin{align} \label{eq:E1n.b1}
b_1& = \min\Bigg\{\frac{1}{2}  - \frac{2(1-\gamma)\nu}{d}- \rho_{22} , ~ \frac{1-\gamma}{2} - \varsigma , ~\frac{1}{4} + \frac{(1-\gamma)(\ell_{\star}-2\nu)}{d}- \frac{\rho_{1} +\rho_{22}}{2} , \nonumber \\
& \qquad \qquad \frac{1-\gamma}{4} + \frac{(1-\gamma)(\ell_{\star}-\nu)}{d} - \frac{\rho_{1}}{2} - \varsigma  \Bigg\}  ,
\end{align}

for an arbitrarily small constant $\varsigma>0$.

(ii) $\widehat{\tau}_n$ in (<ref>) satisfies (<ref>) and (<ref>) in Assumption <ref>, where $b_2$ is

\begin{align} \label{eq:E1n.b2}
b_2& =  \min\Bigg\{\frac{1}{2} , ~ \frac{(1-\gamma)(4\nu+d)}{2d} - \rho_{21} - \varsigma , ~ \frac{1}{4} +\frac{ (1-\gamma)\ell_{\star}}{d} - \frac{\rho_{1} + \rho_{21}}{2}, ~ \nonumber \\
& \qquad \qquad \frac{(1-\gamma)(4\nu+d +4\ell_{\star})}{4d} - \frac{\rho_{1} + 2\rho_{21} }{2}  - \varsigma  \Bigg\} ,
\end{align}

for an arbitrarily small constant $\varsigma>0$.

Theorem <ref> shows that the higher-order quadratic variation estimators $\widehat{\theta}_n$ and $\widehat{\tau}_n$ satisfy Assumption <ref>. We explicitly give the values of $b_1$ and $b_2$ in the four exponential tail inequalities in Assumption <ref>, which directly determine the posterior contraction rates for $\theta$ and $\tau$ in Theorem <ref>. Further simplified values of $b_1$ and $b_2$ will be given in Section <ref>.

Theorem <ref> demonstrates several advantages of the higher-order quadratic variation estimators $\widehat{\theta}_n$ and $\widehat{\tau}_n$ for constructing the exponentially consistent tests needed for our Bayesian theory. First, they satisfy Assumption <ref> with the sieve $\Ecal_n$ defined in (<ref>), where the range parameter $\alpha$ lies in an expanding interval $[n^{-\rho_{31}},n^{\rho_{32}}]$ which eventually covers the entire $\RR_+$ as $n\to\infty$. Though varying $\alpha$ is technically challenging in studying frequentist estimators (9 and 46), we do not need to fix $\alpha$ at a given value and we can assign a general prior on $\alpha$. Second, the inequalities in Assumption <ref> hold uniformly over all possible stratified sampling designs in Assumption <ref>, for arbitrary values of $\delta_{\bfi;k}\in [0,1)$. This generality significantly broadens the applicability of our results to real-world spatial data, as we do not need the sampling locations $S_n$ to be exactly on equispaced grids as in many frequentist fixed-domain asymptotics works. Third, Theorem <ref> holds for any domain dimension $d\in \ZZ_+$, which includes both the case of $d\in \{1,2,3\}$ when the range parameter $\alpha$ cannot be consistently estimated (51), and the case of $d\geq 5$ when $\alpha$ can be consistently estimated (1).

§.§ Explicit Posterior Contraction Rates for Isotropic Matérn

We recall that if Assumptions <ref>-<ref> hold, then Theorem <ref> shows that the posterior contraction rates for $\theta$ and $\tau$ are $n^{-b_1}\log n$ and $n^{-b_2}\log n$, respectively. Theorem <ref> provides the explicit values for $b_1$ and $b_2$ based on the higher-order quadratic variation estimators in (<ref>). In the following, we simplify their general expressions in (<ref>) and (<ref>). The positive constants $\rho_1,\rho_{21},\rho_{22},\rho_{31},\rho_{32}$ need to satisfy (<ref>), but can be all taken as arbitrarily small and close to zero. The condition on $\gamma$ in Theorem <ref> is equivalent to $1-\gamma\in \big(0,\min\{d/(4\nu),1\}\big)$, and $\gamma$ is related to the cell size $\omega_m=\lfloor m^{\gamma} \rfloor$ of the higher-order quadratic variation estimators. We also emphasize that the value of $\gamma$ can be chosen differently for the two estimators $\widehat\theta_n$ and $\widehat\tau_n$, since in Assumption <ref>, the inequalities (<ref>) and (<ref>) for $\widehat\theta_n$ and (<ref>) and (<ref>) for $\widehat\tau_n$ are fully separate and do not affect one another.

First we consider $b_1$ defined in (<ref>), which is related to the posterior contraction rate of the microergodic parameter $\theta$. When the constants $\rho_1,\rho_{21},\rho_{22},\rho_{31},\rho_{32},\varsigma$ are very close to zero, we can show that the first two terms in the minimum are smaller than the last two terms, i.e., $b_1\approx \min\big\{0.5 - 2(1-\gamma)\nu/d, (1-\gamma)/2 \big\}$. We can choose $\gamma$ to balance the two terms by setting $1-\gamma=1/(4\nu/d+1)$, such that $b_1\approx 1/\{2(4\nu/d+1)\}$. From our general Theorem <ref>, the posterior contraction rate of $\theta$ then becomes $n^{-b_1}\log n\approx n^{-1/\{2(4\nu/d+1)\}} \log n$. This rate seems to be slightly slower than the rate $n^{-1/\{2(2\nu/d+1)\}}$ for the MLE in Theorem 5 of 46, where the observations $Y_n$ are assumed to be made on an equispaced grid. Our stratified sampling design in Assumption <ref> does not require the equispaced grid. Meanwhile, it is completely unknown whether the faster rate $n^{-1/\{2(2\nu/d+1)\}}$ for $\theta$ can still hold under our stratified sampling design.  We emphasize that the implied posterior contraction rate for $\theta$ (Theorem <ref> below) is the first of its kind in the Bayesian literature, and the optimal convergence rate for such stratified sampling design in Assumption <ref> remains an open problem.

Next, we consider $b_2$ defined in (<ref>), which is related to the posterior contraction rate of the nugget parameter $\tau$. Again, we set the constants $\rho_1,\rho_{21},\rho_{22},\rho_{31},\rho_{32},\varsigma$ to be very close to zero. For $\gamma$, we choose it to be very close to its lower bound $\max\{0,1-d/(4\nu)\}$, such that the minimum in (<ref>) is equal to $1/2$. As a result, from our general Theorem <ref>, the largest possible posterior contraction rate for the nugget parameter $\tau$ implied by our Theorem <ref> is $n^{-1/2}\log n$. This is almost the same rate (up to the logarithm factor) as the parametric rate $n^{-1/2}$ for the MLE of $\tau$ as shown for the special case of $d=1,\nu=1/2$ in 9 and for the general cases in Theorem 5 of 46.

We summarize the analysis above for the posterior contraction rates in the following theorem, which works for the general stratified sampling design in Assumption <ref> and for any domain dimension $d\in \ZZ_+$.

Suppose that $X(\cdot)$ in Model (<ref>) has the isotropic Matérn covariance function in Example <ref>. Suppose that Assumptions <ref>, <ref>, <ref>, <ref> and <ref> hold with $\kappa=1$, $\Ecal_n$ defined in (<ref>) and $\rho_1,\rho_{21},\rho_{22},\rho_{31},\rho_{32}$ satisfying

\begin{align} \label{eq:rho.relation2}
& 0< \rho_1 < \frac{d}{4\nu+d},   ~~ 0< \rho_{21} <
\min\left(\frac{2\nu}{4\nu+d}, \frac{d}{8\nu} \right) ,
~~ 0<\rho_{22} < \frac{d}{2(4\nu+d)},   \nonumber \\
&  \rho_{31}>0, ~~ 0<\rho_{32} < \frac{1}{4\nu+d},~~ \rho_{1} + \rho_{21} < \min\left(\frac{1}{2}+\frac{2\nu}{d}, \frac{d}{4\nu}\right) .
\end{align}

Then uniform over all $\delta_{\bfi;k}\in [0,1)$ with $\bfi=(i_1,\ldots,i_d)^\T $ for all $1\leq i_1,\ldots,i_d\leq m$ and $1\leq k \leq d$, for any positive sequence $M_n\to \infty$ as $n\to\infty$,

\[\Pi \left(  |\theta/\theta_0 - 1| < M_n n^{-\frac{1}{2(4\nu/d+1)}+\varrho} \log n, \text{ and } |\tau/\tau_0 - 1|< M_n n^{-\frac{1}{2}}\log n  \mid  Y_n, F_n  \right) \to 1,  \]

almost surely $\PP_{(\theta_0,\alpha_0,\tau_0,\beta_0)}$ as $n\to\infty$, where
$\varrho  = \max(\rho_1/2, \rho_{22})$ and $\rho_1,\rho_{22}$ are as defined in (<ref>) and (<ref>).

We remark on the prior Assumptions <ref> and <ref>. They can be verified by many commonly used priors on $(\theta,\alpha,\tau,\beta)$, such as the priors in the following proposition. Let $\text{IG}(a,b)$ be the inverse gamma distribution with shape parameter $a>0$ and rate parameter $b>0$. Let $\text{IGauss}(\mu,\lambda)$ be the inverse Gaussian distribution with mean parameter $\mu>0$ and shape parameter $\lambda>0$. Then we have the following proposition.

Suppose that $X(\cdot)$ in Model (<ref>) has the isotropic Matérn covariance function in Example <ref>. Suppose that the independent priors are assigned on $\theta,\alpha,\tau,\beta$, where $\beta\sim \Ncal(0,a_0I_p)$, $\theta\sim \textup{IG}(a_1,b_1)$, $\tau\sim \textup{IG}(a_2,b_2)$, $\alpha\in \textup{IGauss}(\mu,\lambda)$, for some positive constant hyperparameters $a_0,a_1,b_1,a_2,b_2,\mu,\lambda$ that satisfy $a_1 > 2(3p+6)/\rho_{21}$ and $a_2 > 2(3p+6)/\rho_{22}$. Then Assumption <ref> is satisfied, and Assumption <ref> is satisfied by the set $\Ecal_{n}$ defined in (<ref>) with $\kappa=1$.

§ NUMERICAL EXPERIMENTS

We investigate the posterior contraction behavior for the posterior distribution of parameters in Model (<ref>). In particular, we focus on how fast the marginal posteriors of the microergodic parameter $\theta$ and the nugget parameter $\tau$ contracts towards their true values as we increase the sample size $n$. Simulation studies for the frequentist properties of the proposed higher-order quadratic variation estimators $\widehat\theta_n$ and $\widehat\tau_n$ in Section <ref> can be found in the PhD thesis 45. The simulations below will focus exclusively on the Bayesian posterior contraction for $\theta$ and $\tau$ as well as the Bayesian posterior predictive performance.

§.§ Simulations

We present the results for Model (<ref>) with domain dimension $d=2$ with the isotropic Matérn covariance function in (<ref>). We also have additional simulation results for $d=1$ in Section S8 of the Supplementary Material. We consider two values of the smoothness parameter $\nu=1/2$ and $\nu=1/4$, characterizing different smoothness of the Gaussian process sample paths. The true covariance parameters are set to be $\theta_0=5,\alpha_0=1,\tau_0=0.5$ for both $\nu=1/2$ and $\nu=1/4$. For $d=2$ and the domain $[0,1]^2$, we choose the sampling points $S_n$ to be the regular grid $\left((2i-1)/(2m),(2j-1)/(2m)\right)$ for $i,j=1,\ldots,m$, where we choose $m=20, 22, 25, 28, 31, 35, 39, 44, 49, 55$ such that the sample size roughly follows the geometric sequence $n= m^2\approx 400\times 1.25^{k-1}$ for $k=1,\ldots,10$. For the regression functions, we let $\ff(\bfs)=\left(1,s_1,s_2,s_1^2,s_1s_2,s_2^2\right)^\T$ for $\bfs=(s_1,s_2)\in [0,1]^2$ and $\beta_0=(1,-1.5,-1.5,2,1,2)^\T$.

For Bayesian inference, we assign the prior specified in Proposition <ref> on $(\theta,\alpha,\tau,\beta)$, with hyperparameters $a_0=10^6$, $a_1=b_1=a_2=b_2=0.1$, and $\mu=\lambda=1$. Such values of $a_1,a_2$ are small and do not satisfy the sufficient conditions in Proposition <ref>, but we show that this does not affect the convergence results. We integrate out $\beta$ given the conjugate normal prior and then use the random walk Metropolis algorithm to draw $2000$ samples of $(\theta,\alpha,\tau)$ after $1000$ burnins from the posterior density $\pi(\theta,\alpha,\tau|Y_n)$. We simulate Model (<ref>) for 40 independent copies of the dataset $(Y_n,F_n)$ and find their posterior distributions. The results are summarized in Figure <ref> for $\nu=1/2$ and Figure <ref> for $\nu=1/4$. The boxplots are the marginal posterior distributions of $\theta$ and $\tau$, obtained by averaging over the 40 macro replications of posterior distributions using the Wasserstein-2 barycenter (25). Clearly in both $\nu=1/2$ and $\nu=1/4$ cases, the posterior distribution contracts to the true parameter $\theta_0=5,\tau_0=0.5$ as $n$ increases. The right panels of Figures <ref> and Figure <ref> are the means of absolute differences from all posterior draws of $(\theta,\tau)$ to the true parameters $(\theta_0,\tau_0)$ versus the sample size on the logarithm scale. We can see that they approximately decrease in straight lines on the logarithm scale after $n$ becomes larger. These linear trends indicate that the posterior contractions for both $\theta$ and $\tau$ happen at polynomial rates and hence corroborate our theory.

Posterior contraction for $\nu=1/2$. Left column: Boxplots for the marginal posterior densities of $\theta$ and $\tau$ versus the increasing sample size $n$. The grey dashed lines are the true parameters $\theta_0=5$ and $\tau_0=0.5$. Right column: Posterior means of $|\theta-\theta_0|$ and $|\tau-\tau_0|$ versus the increasing sample size $n$, on the logarithm scale. The dashed lines are the linear regression fits, and the grey shaded areas are the 95% confidence bands. All posterior summaries are averaged over 40 macro Monte Carlo replications.

Posterior contraction for $\nu=1/4$. Left column: Boxplots for the marginal posterior densities of $\theta$ and $\tau$ versus the increasing sample size $n$. The grey dashed lines are the true parameters $\theta_0=5$ and $\tau_0=0.5$. Right column: Posterior means of $|\theta-\theta_0|$ and $|\tau-\tau_0|$ versus the increasing sample size $n$, on the logarithm scale. The dashed lines are the linear regression fits, and the grey shaded areas are the 95% confidence bands. All posterior summaries are averaged over 40 macro Monte Carlo replications.

We further investigate the Bayesian posterior prediction performance and compare with the best possible prediction. We draw another $N=2500$ points $\bfs_1^*,\ldots,\bfs_N^*$ uniformly from the domain $[0,1]^2$ as the testing locations, and compute the prediction mean squared error

$$\EE_{Y_n}  \Bigg\{N^{-1}\sum_{l=1}^N M_{\text{post}}(\bfs_l^*)\Bigg\} =\EE_{Y_n} \Bigg\{ N^{-1}\sum_{l=1}^N \EE_{\theta,\alpha,\tau,\beta\mid Y_n,F_n} \big\{(\tilde Y(\bfs_l^*)-Y_0(\bfs_l^*)\big\}^2\Bigg\}, $$

where $\tilde Y(\bfs_l^*)$ is a draw from the Bayesian predictive posterior distribution of $Y(\cdot)$ at $\bfs_l^*$ given a random draw $(\theta,\alpha,\tau,\beta)$ from the posterior distribution $\Pi(\cdot\mid Y_n,F_n)$, $Y_0(\bfs^*) = \beta_0^\T \ff(\bfs^*) + X(\bfs^*)$ is the true mean function at $\bfs^*$ without measurement error, and the two layers of expectations are with respect to both the posterior distribution and the distribution of $Y_n$ given the true parameters $(\theta_0,\alpha_0,\tau_0,\beta_0)$. The oracle prediction mean squared error for Model (<ref>) is calculated as $\EE_{Y_n}  \big\{N^{-1}\sum_{l=1}^N M_{0}(\bfs_l^*)\big\} = \EE_{Y_n} \big[N^{-1}\sum_{l=1}^N \big\{(\tilde Y_0(\bfs_l^*)-Y_0(\bfs_l^*)\big\}^2\big]$ where $\tilde Y_0(\cdot)$ is the best linear unbiased predictor of $Y(\cdot)$ given the true parameters $(\theta_0,\alpha_0,\tau_0,\beta_0)$ (42). The exact formulas to calculate $\tilde Y(\bfs_l^*),M_{\text{post}}(\bfs_l^*),\tilde Y_0(\bfs_l^*),M_{0}(\bfs_l^*)$ for Model (<ref>) can be found in Section S8 of the Supplementary Material.

In the left panel of Figure <ref>, we plot the prediction mean squared errors for both $\nu=1/2$ and $\nu=1/4$ together with those from $(\theta_0,\alpha_0,\tau_0,\beta_0)$ on the logarithm scale. Clearly the Bayesian predictive posterior has almost the same mean square errors as the oracle prediction. We further take the ratio of the two prediction mean squared errors and plot $\EE_{Y_n} \big\{N^{-1}\sum_{l=1}^N M_{\text{post}}(\bfs_l^*)/M_{0}(\bfs_l^*)\big\}$ in the right panel of Figure <ref>, which characterizes the relative efficiency of Bayesian prediction. 38, 39, 40, 41 have developed the frequentist theory of posterior asymptotic efficiency for Gaussian processes without regression terms and nugget. The decreasing trend of the ratio towards 1 in the right panel of Figure <ref> indicates a possibly similar phenomenon that the Bayesian posterior prediction is asymptotically efficient even for the more general spatial model (<ref>) under the fixed-domain asymptotics framework.

Prediction mean squared errors for $\nu=1/2$ and $\nu=1/4$. Left panel: The prediction mean squared errors under both the Bayesian posterior prediction and the oracle prediction based on the true parameters. Right panel: Ratios of the Bayesian prediction mean squared error and the oracle prediction mean squared error. All posterior summaries are averaged over 2500 testing locations in $[0,1]^2$ and 40 macro Monte Carlo replications.

§.§ Sea Surface Temperature Data

For a real data analysis, we fit Model (<ref>) and the isotropic Matérn covariance function (<ref>) with $\nu=1/2$ to the sea surface temperature data on the Pacific Ocean between $45^\circ$–$48^\circ$ north latitudes and $149^\circ$–$152^\circ$ west longitudes on August 16, 2016. The data were collected from remote sensing satellites with a high-resolution on a $0.025^\circ\times 0.025^\circ$ grid, and can be obtained from National Oceanographic Data Centres (NODC) World Ocean Database (<https://www.ncei.noaa.gov/products/world-ocean-database>). For Bayesian estimation, we use a total of 3600 observations on the $0.05^\circ\times 0.05^\circ$ grid, and set the sample size $n = 700, 1050, 1600, 2400, 3600$. For each $n$ smaller than 3600, we randomly choose 10 subsets of data $Y_n$. We set the regressors $\ff(\bfs)=(1,s_1,s_2)^\T$ for the latitude $s_1$ and the longitude $s_2$. We assign the same prior as in Proposition <ref> with hyperparameters $a_0=10^6$, $a_1=b_1=a_2=b_2=0.1$, and $\mu=\lambda=1$. We draw $4000$ posterior samples of $(\theta,\alpha,\tau)$ after $1000$ burnins for each of the 10 subsets, and then average the 10 marginal posterior distributions into one summary posterior distribution using the Wasserstein-2 barycenter (25). Figure <ref> presents the summary marginal posterior densities of $\theta$ and $\tau$ for $n\in \{700, 1050, 1600, 2400, 3600\}$. As $n$ increases, we can see the clear trend of posterior contraction for both $\theta$ and $\tau$, and the posterior of $\tau$ seems to contract faster than that of $\theta$.

Sea surface temperature data. Left panel: Temperature in the target region. Middle panel: Marginal posterior densities of $\theta$ with various sample size $n$. Right panel: Marginal posterior densities of $\tau$ with various sample size $n$.

§ DISCUSSION

The general fixed-domain posterior contraction theory developed in this paper can be potentially applied to many stationary covariance functions in spatial statistics. With the new evidence lower bound, we have effectively transformed the problem of finding posterior contraction rates to the problem of finding efficient frequentist estimators that satisfy the concentration inequalities with exponentially small tails as in Assumption <ref>.

There are many potential directions to extend the current work. First, for isotropic Matérn, our higher-order quadratic variation estimators deliver the explicit posterior contraction rates, though it remains unclear what is the optimal rate for the microergodic parameter $\theta$ under our flexible stratified sampling design. It would be of further interest to find the limiting Bayesian posterior distribution of $(\theta,\alpha,\tau,\beta)$ and establish the posterior asymptotic normality for both $\theta$ and the nugget $\tau$, similar to the Bayesian fixed-domain asymptotic theory for the model without nugget in 24.

Second, we have assumed that the smoothness parameter $\nu$ is fixed and known in our theory. Because the smoothness parameter $\nu$ determines the degree of mean square differentiability of the random field, estimation of $\nu$ has been an important and meanwhile challenging problem in the spatial literature. Recently, using the new higher-order quadratic variation method, 28, 30, and 29 have proposed consistent estimators of $\nu$ for irregularly spaced spatial data under fixed-domain asymptotics. In fact, our estimator of the microergodic parameter $\theta$ in Section <ref> is the same as the estimator in 30 and 29 when $\nu$ is assumed to be known. It requires further study whether we can include $\nu$ as part of the unknown parameters in the  Bayesian framework, and establish the posterior contraction for $\nu$ based on the new estimators.

Third, the isotropic covariance function considered in this paper leads to the simplest spatial Gaussian process model for real data. To extend our method to anisotropic covariance functions or even nonstationary spatial processes, we need to both establish the evidence lower bound similar to Theorem <ref> and find some consistent frequentist estimators for the microergodic parameters in these more general covariance functions under fixed-domain asymptotics. From the technical perspective, it seems that finding consistent estimators will be a more challenging task than deriving the evidence lower bound, though both of them may require a case-by-case analysis for different spatial covariance functions. Meanwhile, it is possible to adapt the current proof techniques to a general nonparametric class of covariance functions based the principal irregular terms in their Taylor expansions (see Section 2.7 of 42), in a similar spirit to the recent work 4.

Fourth, given the frequentist asymptotic efficiency in 39, 40 for Gaussian processes without nugget as well as our promising simulation results, it would be of theoretical interest to show that the asymptotic efficiency can still be preserved even in the Bayesian posterior prediction for the true mean function $Y_0(\bfs) = \beta_0^\T \ff(\bfs) + X(\bfs)$ in the more general model (<ref>) with both regression terms and nugget.

These potential developments will together provide strong theoretical justification for the existing Bayesian spatial inference based on the Gaussian process regression model, including parameter estimation, uncertainty quantification, and prediction. We leave these directions for future research.

§ ACKNOWLEDGEMENTS

The authors thank Professor Wei-Liem Loh for helpful discussion. This work was supported by Singapore Ministry of Education Academic Research Funds Tier 1 Grant A-0004822-00-00.

§ SUPPLEMENTARY MATERIAL

Supplementary Material for “Fixed-domain Posterior Contraction Rates for Spatial Gaussian Process Model with Nugget”: Technical proofs of all theorems, propositions, and additional simulation results.

Supplementary Material for “Fixed-domain Posterior Contraction Rates for Spatial Gaussian Process Model with Nugget"

This supplementary material contains the technical proofs of the theorems and propositions in the main paper as well as additional simulation results. Section <ref> contains the proof of Theorem <ref>. Section <ref> contains the proof of Theorem <ref> and auxiliary technical results on spectral analysis. Section <ref> proves a proposition on the posterior inconsistency of the range parameter $\alpha$ and the regression coefficient vector $\beta$. Sections <ref>, <ref>, and <ref> include the proofs of Proposition <ref>, Theorem <ref>, and Proposition <ref>, respectively. Section <ref> contains the lengthy proof of Theorem <ref>. Section <ref> includes the formulas for calculating the Bayesian and oracle prediction mean squared errors for Model (<ref>) in the main text and additional simulation results for the case of domain dimension $d=1$.

We define some universal notation that will be used throughout the proof. Let $\RR_+=(0,+\infty)$, $\ZZ_+$ be the set of all positive integers, and $\NN=\ZZ_+\cup \{0\}$. For any $x\in \RR$, $\lceil x\rceil$ and $\lfloor x\rfloor$ denote the smallest integer $\geq x$ and the largest integer $\leq x$. For any $x=(x_1,\ldots,x_d)^\T  \in \RR^d$, we let $\|x\|=\big(\sum_{i=1}^d x_i^2\big)^{1/2}$, $\|x\|_1=\sum_{i=1}^d |x_i|$, and $\|x\|_{\infty}=\max(x_1,\ldots,x_d)$. For a generic set $\Ecal$, we use $\Ecal^c$ to denote its complement. If $\Ecal$ contains finitely many elements, then $|\Ecal|$ denotes its cardinality.

For two positive sequences $a_n$ and $b_n$, we use $a_n\prec b_n$ and $b_n\succ a_n$ to denote the relation $\lim_{n\to\infty} a_n/b_n=0$. $a_n=O(b_n)$, $a_n\preceq b_n$ and $b_n\succeq a_n$ denote the relation $\limsup_{n\to\infty} a_n/b_n<+\infty$. $a_n\asymp b_n$ denotes the relation $a_n\preceq b_n$ and $a_n\succeq b_n$.

For any integers $k,m$, we let $I_k$ be the $k\times k$ identity matrix, $0_k$ be the $k$-dimensional column vectors of all zeros, $0_{k\times m}$ be the $k\times m$ zero matrix, and $\diag\{c_1,\ldots,c_k\}$ be the diagonal matrix with diagonal entries $c_1,\ldots,c_k\in \RR$. For any generic matrix $A$, $cA$ denotes the matrix of $A$ with all entries multiplied by the number $c$, and $\dett(A)$ denotes the determinant of $A$. For a square matrix $A$, $\tr(A)$ denotes the trace of $A$. If $A$ is symmetric positive semidefinite, then $\mathsf{s}_{\min}(A)$ and $\mathsf{s}_{\max}(A)$ denote the smallest and largest eigenvalues of $A$, and $A^{1/2}$ denotes a symmetric positive semidefinite square root of $A$. For two symmetric matrices $A$ and $B$, we use $A\preceq B$ and $B\succeq A$ to denote the relation that $B-A$ is symmetric positive semidefinite, and use $A\prec B$ and $B\succ A$ to denote the relation that $B-A$ is symmetric positive definite. For any matrix $A$, $\|A\|_{\op}=\big\{\mathsf{s}_{\max}(A^\T  A)\big\}^{1/2} = \big\{\mathsf{s}_{\max}(AA^\T )\big\}^{1/2}$ denotes the operator norm of $A$, and $\|A\|_F= \big\{\tr\left(A^\T  A\right)\big\}^{1/2}$ denotes the Frobenius norm of $A$. We use $\pr(\cdot)$ to denote the probability under the probability measure $\PP_{(\theta_0,\alpha_0,\tau_0,\beta_0)}$.

§ PROOF OF THEOREM <REF>

We first show the propriety of the posterior distribution defined in Equation (<ref>) of the main text, which is equivalent to showing that the integral in the denominator on the right-hand side of (<ref>) is finite. Since for two positive definite matrices $A$ and $B$, $\dett(A+B)\geq \dett(A)$, we have that for any given $n\in \ZZ_+$,

\begin{align*}
&\quad~ \int_{\RR_+^{3}\times \RR^p} \exp\left\{\Lcal_n(\theta,\alpha,\tau,\beta) \right\} \pi(\theta,\alpha,\tau,\beta) \ud \theta \ud \tau \ud \alpha \ud \beta  \\
&= \int_{\RR_+^{3}\times \RR^p} \frac{1}{\left[\dett\left\{\theta K_{\alpha,\nu}(S_n)+\tau I_n\right\}\right]^{1/2}} \exp\left\{ - \frac{1}{2} (Y_n-  F_n  \beta)^\T  \left\{\theta K_{\alpha,\nu}(S_n) + \tau I_n \right\}^{-1} (Y_n-  F_n  \beta) \right\} \\
&\qquad \times \pi(\theta,\alpha,\tau,\beta) \ud \theta \ud \tau \ud \alpha \ud \beta  \\
&\leq \int_{\RR_+^{3}\times \RR^p} \frac{1}{\left\{\dett(\tau I_n)\right\}^{1/2}} \cdot 1 \cdot \pi(\theta,\alpha,\tau,\beta) \ud \theta \ud \tau \ud \alpha \ud \beta  \\
&= \int_0^{\infty} \tau^{-n/2} \pi(\tau) \ud \tau < \infty,
\end{align*}

where the last integral is finite following Assumption <ref>. This proves that the joint posterior density $\pi(\theta,\alpha,\tau,\beta~|~ Y_n, F_n)$ in Equation (<ref>) is well defined.

Next, we prove the posterior contraction rates for $\theta$ and $\tau$ in Theorem <ref>. We follow the classical proof of the Schwartz's theorem for posterior consistency (37); see for example, Theorem 6.17 in 13. Let $\varepsilon_{1n}= M_n n^{-b_1} \log n$ and $\varepsilon_{2n}= M_n n^{-b_2} \log n$. Define the testing function

\begin{align} \label{eq:test}
\phi_n & = \Ical \left( |\widehat\theta_n/\theta_0-1|\geq \varepsilon_{1n}/2, \text{ or } |\widehat\tau_n/\tau_0-1|\geq \varepsilon_{2n}/2 \right),
\end{align}

where $\widehat\theta_n$ and $\widehat\tau_n$ are from Assumption <ref>.

We start with the decomposition

\begin{align}\label{eq:cons1}
& \Pi \left(\Bcal_0(\varepsilon_{1n},\varepsilon_{2n})^c ~|~ Y_n, F_n  \right)
= \frac{\int_{\Bcal_0(\varepsilon_{1n},\varepsilon_{2n})^c}\exp\left\{ \Lcal_n(\theta,\alpha,\tau,\beta) \right\} \pi(\theta,\alpha,\tau,\beta) \ud \theta \ud \tau \ud \alpha \ud \beta }
{\int_{\RR_+^{3} \times \RR^p} \exp\left\{\Lcal_n(\theta',\alpha',\tau',\beta') \right\} \pi(\theta',\alpha',\tau',\beta')\ud \theta' \ud \tau' \ud \alpha' \ud \beta' } \nonumber \\
&\leq \phi_n + \frac{(1-\phi_n)\int_{\Bcal_0(\varepsilon_{1n},\varepsilon_{2n})^c}\exp\left\{ \Lcal_n(\theta,\alpha,\tau,\beta)\right\} \pi(\theta,\alpha,\tau,\beta) \ud \theta \ud \tau \ud \alpha \ud \beta}
{\int_{\RR_+^{3} \times \RR^p} \exp\left\{\Lcal_n(\theta',\alpha',\tau',\beta') \right\} \pi(\theta',\alpha',\tau',\beta')\ud \theta' \ud \tau' \ud \alpha' \ud \beta' } .
\end{align}

By Assumption <ref>, we have that for some large positive constant $N_0$, for all $n>N_0$,

\begin{align*}
{\EE}_{(\theta_0,\alpha_0,\tau_0,\beta_0)} (\phi_n) & \leq \exp\left\{- c_1 \varphi\left(n^{b_1} \varepsilon_{1n} \right)\right\} + \exp\left\{- c_2 \varphi\left(n^{b_2} \varepsilon_{2n} \right)\right\} \\
& \leq \exp\left(- c_1 M_n \log n \right) + \exp\left(- c_2 M_n \log n  \right)  \to 0,
\end{align*}

where the last step follows from $M_n \to +\infty$ as $n\to\infty$. Notice that for any constant $c>0$ and an arbitrary $c'>1$, $\exp\left(-c M_n \log n\right)\leq \exp(-c' \log n) = n^{-c'}$ for all sufficiently large $n$ and hence the sequence $\exp\left(-c M_n \log n\right)$ is summable over $n$. Therefore, by the Markov's inequality and the Borel-Cantelli lemma, $\phi_n\to 0$ almost surely $\PP_{(\theta_0,\alpha_0,\tau_0)}$ as $n\to\infty$.

We now turn to the second term in (<ref>). It can be further decomposed into two terms:

\begin{align}\label{eq:cons2}
& \frac{(1-\phi_n)\int_{\Bcal_0(\varepsilon_{1n},\varepsilon_{2n})^c} \exp\left\{ \Lcal_n(\theta,\alpha,\tau,\beta)\right\} \pi(\theta,\alpha,\tau,\beta) \ud \theta \ud \tau \ud \alpha \ud \beta}
{\int_{\RR_+^{3} \times \RR^p} \exp\left\{\Lcal_n(\theta',\alpha',\tau',\beta') \right\} \pi(\theta',\alpha',\tau',\beta')\ud \theta' \ud \tau' \ud \alpha' \ud \beta'}  = T_1+T_2, \quad \text{ where }\nonumber \\
& T_1 = \frac{(1-\phi_n)\int_{\Bcal_0(\varepsilon_{1n},\varepsilon_{2n})^c \cap \Ecal_n} \exp\left\{ \Lcal_n(\theta,\alpha,\tau,\beta)-\Lcal_n(\theta_0,\alpha_0,\tau_0,\beta_0) \right\} \pi(\theta,\alpha,\tau,\beta) \ud \theta \ud \tau \ud \alpha \ud \beta }
{\int_{\RR_+^{3} \times \RR^p} \exp\left\{\Lcal_n(\theta',\alpha',\tau',\beta') -\Lcal_n(\theta_0,\alpha_0,\tau_0,\beta_0)\right\} \pi(\theta',\alpha',\tau',\beta')\ud \theta' \ud \tau' \ud \alpha' \ud \beta' } ,  \nonumber \\
&T_2 =  \frac{(1-\phi_n)\int_{\Bcal_0(\varepsilon_{1n},\varepsilon_{2n})^c \cap \Ecal_n^c} \exp\left\{ \Lcal_n(\theta,\alpha,\tau,\beta)-\Lcal_n(\theta_0,\alpha_0,\tau_0,\beta_0) \right\} \pi(\theta,\alpha,\tau,\beta) \ud \theta \ud \tau \ud \alpha \ud \beta  }
{\int_{\RR_+^{3}\times \RR^p} \exp\left\{\Lcal_n(\theta',\alpha',\tau',\beta')-\Lcal_n(\theta_0,\alpha_0,\tau_0,\beta_0) \right\} \pi(\theta',\alpha',\tau',\beta')\ud \theta' \ud \tau' \ud \alpha' \ud \beta'} .
\end{align}

For the first term $T_1$ in (<ref>), by Assumption <ref> and the Fubini's theorem, the numerator has expectation

\begin{align}\label{eq:cons2.1}
& {\EE}_{(\theta_0,\alpha_0,\tau_0,\beta_0)} (1-\phi_n)\int_{\Bcal_0(\varepsilon_{1n},\varepsilon_{2n})^c \cap \Ecal_n} \exp\Big\{ \Lcal_n(\theta,\alpha,\tau,\beta) \nonumber \\
&~~ -\Lcal_n(\theta_0,\alpha_0,\tau_0,\beta_0)  \Big\} \pi(\theta,\alpha,\tau,\beta) \ud \theta \ud \tau \ud \alpha \ud \beta \nonumber \\
={} & \int_{\Bcal_0(\varepsilon_{1n},\varepsilon_{2n})^c \cap \Ecal_n} {\EE}_{(\theta,\alpha,\tau,\beta)} (1-\phi_n) \pi(\theta,\alpha,\tau,\beta) \ud \theta \ud \tau \ud \alpha \ud \beta \nonumber \\
\leq{}& \sup_{\Bcal_0(\varepsilon_{1n},\varepsilon_{2n})^c \cap \Ecal_n} {\EE}_{(\theta,\alpha,\tau,\beta)} \left\{(1-\phi_n) \cdot \int_{\Bcal_0(\epsilon_1,\epsilon_2)^c \cap \Ecal_n} \pi(\theta,\alpha,\tau,\beta) \ud \theta \ud \tau \ud \alpha \ud \beta \right\} \nonumber \\
\leq{}& \sup_{\Bcal_0(\varepsilon_{1n},\varepsilon_{2n})^c \cap \Ecal_n} {\EE}_{(\theta,\alpha,\tau,\beta)} (1-\phi_n) \nonumber \\
\leq{}& \sup_{\Bcal_0(\varepsilon_{1n},\varepsilon_{2n})^c \cap \Ecal_n} \PP_{(\theta,\alpha,\tau,\beta)} \left( |\widehat\theta_n/\theta_0-1|\leq \varepsilon_{1n}/2 \right)
+ \sup_{\Bcal_0(\varepsilon_{1n},\varepsilon_{2n})^c \cap \Ecal_n} \PP_{(\theta,\alpha,\tau,\beta)} \left( |\widehat\tau_n/\tau_0-1|\leq \varepsilon_{2n}/2 \right) \nonumber \\
\leq{}& \exp\left\{- c_1 \varphi\left(n^{b_1}\varepsilon_{1n} \right)\right\} + \exp\left\{- c_2 \varphi\left(n^{b_2}\varepsilon_{2n} \right)\right\} \nonumber \\
={}& \exp\left(-c_1 M_n \log n\right) + \exp\left(-c_2 M_n \log n\right) .
\end{align}

Therefore, by applying the Markov's inequality and the Borel-Cantelli Lemma to (<ref>), the numerator of $T_1$ in (<ref>) is smaller than $\exp\left\{-\min(c_1,c_2) M_n \log n / 2\right\}$ almost surely $\PP_{(\theta_0,\alpha_0,\tau_0,\beta_0)}$ as $n\to\infty$. On the other hand, Theorem <ref> shows that almost surely $\PP_{(\theta_0,\alpha_0,\tau_0,\beta_0)}$ as $n\to\infty$, for all $d\in \ZZ_+$, the denominator in $T_1$ in (<ref>) is lower bounded by

\begin{align}\label{eq:cons2.2}
& \int_{\RR_+^{3}\times \RR^p} \exp\left\{\Lcal_n(\theta',\alpha',\tau',\beta')-\Lcal_n(\theta_0,\alpha_0,\tau_0,\beta_0)  \right\} \nonumber \\
& ~~ \times  \pi(\theta',\alpha',\tau',\beta')\ud \theta' \ud \tau' \ud \alpha' \ud \beta' \geq D n^{-(3p+2+2/\kappa)},
\end{align}

for some constant $D>0$. Hence for all sufficiently large $n$, the first term $T_1$ in (<ref>) is upper bounded by

\begin{align}\label{eq:cons2.3}
T_1 & \leq D^{-1} n^{(3p+2+2/\kappa)} \exp\left\{-\min(c_1,c_2) M_n \log n / 2\right\}\leq \exp(-c' \log n)  \to 0,
\end{align}

for an arbitrary $c'>1$, almost surely $\PP_{(\theta_0,\alpha_0,\tau_0,\beta_0)}$.

For the second term $T_2$ in (<ref>), by Assumption <ref> and the Fubini's theorem, the numerator has expectation

\begin{align}\label{eq:cons2.4}
& {\EE}_{(\theta_0,\alpha_0,\tau_0,\beta_0)} (1-\phi_n)\int_{\Bcal_0(\varepsilon_{1n},\varepsilon_{2n})^c \cap \Ecal_n^c} \exp\Big\{ \Lcal_n(\theta,\alpha,\tau,\beta) \nonumber \\
&\quad -\Lcal_n(\theta_0,\alpha_0,\tau_0,\beta_0) \Big\} \pi(\theta,\alpha,\tau,\beta) \ud \theta \ud \tau \ud \alpha \ud \beta \nonumber \\
={} & \int_{\Bcal_0(\varepsilon_{1n},\varepsilon_{2n})^c \cap \Ecal_n^c} {\EE}_{(\theta_0,\alpha_0,\tau_0,\beta_0)} (1-\phi_n)\exp\Big\{ \Lcal_n(\theta,\alpha,\tau,\beta) \nonumber \\
&\quad -\Lcal_n(\theta_0,\alpha_0,\tau_0,\beta_0) \Big\} \pi(\theta,\alpha,\tau,\beta) \ud \theta \ud \tau \ud \alpha \ud \beta  \nonumber \\
={} & \int_{\Bcal_0(\varepsilon_{1n},\varepsilon_{2n})^c \cap \Ecal_n^c} {\EE}_{(\theta,\alpha,\tau,\beta)} (1-\phi_n) \pi(\theta,\alpha,\tau,\beta) \ud \theta \ud \tau \ud \alpha \ud \beta \nonumber \\
\leq{}& \int_{\Ecal_n^c} 1\cdot \pi(\theta,\alpha,\tau,\beta) \ud \theta \ud \tau \ud \alpha \ud \beta
= \Pi(\Ecal_n^c) \leq n^{-(3p+4+2/\kappa)} .
\end{align}

Therefore, by applying the Markov's inequality to (<ref>), the numerator of $T_2$ in (<ref>) is smaller than $n^{-\left(3p+\frac{5}{2}+\frac{2}{\kappa}\right)}$ with probability at least $1-n^{-3/2}$. Since the sequence $\{n^{-3/2}:n\geq 1\}$ is summable over $n$, by the Borel-Cantelli Lemma, we have that the numerator of $T_2$ in (<ref>) is smaller than $n^{-\left(3p+\frac{5}{2}+\frac{2}{\kappa}\right)}$ almost surely $\PP_{(\theta_0,\alpha_0,\tau_0,\beta_0)}$ as $n\to\infty$.
We combine this with the lower bound in (<ref>) to conclude that almost surely $\PP_{(\theta_0,\alpha_0,\tau_0,\beta_0)}$ as $n\to\infty$, the second term $T_2$ in (<ref>) is upper bounded by

\begin{align}\label{eq:cons2.5}
T_2 &\leq  D^{-1} n^{3p+2+2/\kappa} \cdot n^{-\left(3p+\frac{5}{2}+\frac{2}{\kappa}\right)} = D^{-1} n^{-1/2} \to 0.
\end{align}

Finally, we combine (<ref>), (<ref>), (<ref>), and (<ref>) to conclude that the right-hand side of (<ref>) converges to zero almost surely $\PP_{(\theta_0,\alpha_0,\tau_0,\beta_0)}$ as $n\to\infty$, which completes the proof of Theorem <ref>.

§ PROOF OF THEOREMS <REF> AND AUXILIARY TECHNICAL LEMMAS

§.§ Proof of Theorem <ref>

Define the sets $\Acal_{n}$ and $\Dcal_{n}$ as

\begin{align}
& \Acal_{n} = \left\{(\theta,\alpha,\tau,\beta)\in \RR_+^{3} \times \RR^p:~ |\alpha/\alpha_0-1|\leq n^{-2/\kappa} \right\} , \nonumber \\
& \Dcal_{n} = \left\{(\theta,\alpha,\tau,\beta)\in \RR_+^{3} \times \RR^p:~ \|\beta-\beta_0\|\leq n^{-3} \right\}. \nonumber
\end{align}

By the continuity of the prior density in Assumption <ref>, we have that for all sufficiently large $n$, $\pi(\theta,\tau,\alpha,\beta)>\pi(\theta_0,\tau_0,\alpha_0,\beta_0)/2>0$ for all $(\theta,\tau,\alpha,\beta)\in \Bcal_0(n^{-1}, n^{-1}) \cap \Acal_{n} \cap \Dcal_{n}$. Therefore,

\begin{align} \label{eq:prior.lower1}
\Pi \left(\Bcal_0(n^{-1}, n^{-1}) \cap \Acal_{n} \cap \Dcal_{n}\right)
&\geq \frac{\pi(\theta_0,\tau_0,\alpha_0,\beta_0)}{2} \cdot (2\theta_0n^{-1}) \cdot (2\tau_0n^{-1}) \cdot (2\alpha_0 n^{-2/\kappa}) \cdot v_0 n^{-3p}  \nonumber \\
& = 4 v_0 \theta_0\tau_0\alpha_0 \pi(\theta_0,\tau_0,\alpha_0,\beta_0) n^{-(3p+2+2/\kappa)},
\end{align}

where $v_0 = \pi^{p/2}/\Gamma(p/2+1) $. Using the definition of log-likelihood in (<ref>), the exponent in Theorem <ref> can be written as

\begin{align}
& \Lcal_n(\theta,\alpha,\tau,\beta) - \Lcal_n(\theta_0,\alpha_0,\tau_0,\beta_0)  \nonumber \\
= {}& -\frac{1}{2} \left\{\tr\left[\left\{\theta K_{\alpha,\nu}(S_n) + \tau I_n\right\}^{-1} \widetilde Y_n \widetilde Y_n^\T  \right] - \tr\left[\left\{\theta_0 K_{\alpha_0,\nu}(S_n) + \tau_0 I_n\right\}^{-1} \widetilde Y_n \widetilde Y_n^\T  \right] \right\} \nonumber \\
& - \frac{1}{2} \log \frac{\dett\left\{\theta K_{\alpha,\nu}(S_n) +\tau I_n\right\}}{\dett\left\{\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right\}} \nonumber \\
& + \widetilde Y_n^\T  \left[\left\{\theta K_{\alpha,\nu}(S_n) +\tau I_n\right\}^{-1}-\left\{\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right\}^{-1}\right] F_n (\beta-\beta_0) \nonumber \\
& - \frac{1}{2} (\beta-\beta_0)^\T   F_n ^\T  \left[\left\{\theta K_{\alpha,\nu}(S_n) +\tau I_n\right\}^{-1}-\left\{\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right\}^{-1} \right] F_n (\beta-\beta_0). \nonumber
\end{align}

Using the lower bounds in Part (i) of Lemma <ref> and Lemma <ref> in Section <ref> below, we obtain that there exists a large integer $N_1$ that depends on $\nu,d,C_{\ff},\alpha_0,\theta_0$, $\tau_0,L,r_0,\kappa$, such that for any $\epsilon_1,\epsilon_2\in (0,1/2)$, for all $n>N_1$, with probability at least $1-3\exp(-\log^2 n)$,

\begin{align} \label{eq:denom.lower1}
& \int_{\Bcal_0(\epsilon_1,\epsilon_2)\cap \Acal_{n}\cap\Dcal_{n}} \exp\left\{\Lcal_n(\theta,\alpha,\tau,\beta) - \Lcal_n(\theta_0,\alpha_0,\tau_0,\beta_0) \right\} \cdot \pi(\theta,\alpha,\tau,\beta)~ \ud \theta \ud \tau \ud \alpha \ud \beta  \nonumber \\
& \geq  \Pi(\Bcal_0(\epsilon_1,\epsilon_2) \cap \Acal_{n} \cap \Dcal_{n}) \cdot \exp\Big\{-  \left(\frac{5}{2}n \epsilon_{1\vee 2} + 5Ln^{-1} \right)  \nonumber \\
& \quad - \left[8C_{\ff}p^{1/2} \left\{\theta_0 K_{\alpha_0,\nu}(0) + \tau_0 \right\}^{1/2} \tau_0^{-1} n^{-1} + 2C_{\ff}^2 p \tau_0^{-1} n^{-4} \right] \Big\} ,
\end{align}

where $\epsilon_{1\vee 2} = \max(\epsilon_1,\epsilon_2)$. We can take $\epsilon_1=\epsilon_2=n^{-1}$ in (<ref>), so $\epsilon_{1\vee 2}=n^{-1}$. We then combine (<ref>) and (<ref>) to obtain that for all sufficiently large $n>N_1$, with probability at least $1-3\exp(-\log^2 n)$,

\begin{align} \label{eq:post.lower1}
& \int \exp\left\{\Lcal_n(\theta,\alpha,\tau,\beta) - \Lcal_n(\theta_0,\alpha_0,\tau_0,\beta_0) \right\} \cdot \pi(\theta,\alpha,\tau,\beta)~ \ud \theta \ud \tau \ud \alpha \ud \beta  \nonumber \\
\geq{}& \int_{\Bcal_0(n^{-1},n^{-1})\cap \Acal_{n} \cap \Dcal_{n}} \exp\left\{\Lcal_n(\theta,\alpha,\tau,\beta) - \Lcal_n(\theta_0,\alpha_0,\tau_0,\beta_0) \right\} \cdot \pi(\theta,\alpha,\tau,\beta)~ \ud \theta \ud \tau \ud \alpha \ud \beta  \nonumber \\
\geq{}& \Pi\left(\Bcal_0\big(n^{-1},n^{-1}\big) \cap \Acal_{n}\cap \Dcal_{n}\right) \cdot \exp\Big\{- \Big(\frac{5}{2}n \cdot n^{-1} + 5Ln^{-1} \Big)  \nonumber \\
& \quad - \left[8C_{\ff}p^{1/2} \left\{\theta_0 K_{\alpha_0,\nu}(0) + \tau_0 \right\}^{1/2} \tau_0^{-1} n^{-1} + 2C_{\ff}^2 p \tau_0^{-1} n^{-4} \right] \Big\} \nonumber \\
\geq{}& \Pi\left(\Bcal_0\big(n^{-1},n^{-1}\big) \cap \Acal_{n}\cap \Dcal_{n}\right)   \exp\left[- 2.5 - 5L - 8C_{\ff}p^{1/2} \left\{\theta_0 K_{\alpha_0,\nu}(0) + \tau_0 \right\}^{1/2} \tau_0^{-1} - 2C_{\ff}^2 p \tau_0^{-1}  \right] \nonumber \\
\geq{}& Dn^{-(3p+2+2/\kappa)},
\end{align}

$$D=4v_0\theta_0\tau_0\alpha_0 \pi(\theta_0,\tau_0,\alpha_0,\beta_0)\cdot \exp\left[-2.5-5L -8C_{\ff}p^{1/2} \left\{\theta_0 K_{\alpha_0,\nu}(0) + \tau_0 \right\}^{1/2} \tau_0^{-1} -2C_{\ff}^2 p \tau_0^{-1}  \right]$$

is a positive constant. This completes the proof of Theorem <ref>.

§.§ Spectral Analysis of Covariance Functions

We present a series of results for the spectral analysis of the covariance function $\theta K_{\alpha,\nu}$ for the Gaussian process $X(\cdot)$, which is used in the proof of Theorem <ref>.

For $w \in \RR^d$, let

\begin{align}\label{f.specden}
f_{\theta,\alpha,\nu}(w) &= \frac{1}{(2\pi)^d}\int_{\RR^d} \exp\left(-\imath w^\T  x\right) \theta K_{\alpha,\nu}(x) \ud x
\end{align}

be the spectral density of the covariance function $\theta K_{\alpha,\nu}$.  For any given pair $(\theta,\alpha)$, let $\|h\|_{f_{\theta,\alpha,\nu}}^2 = \langle h,h \rangle_{f_{\theta,\alpha,\nu}} =\int_{\RR^d} |h(w)|^2 f_{\theta,\alpha,\nu}(w)\ud w$ be the norm of a generic function $h$ in the Hilbert space $L_2(f_{\theta,\alpha,\nu})$, with inner product $\langle h_1,h_2\rangle_{f_{\theta,\alpha,\nu}} =\int_{\RR^d} h_1(w) \overline{h_2(w)} f_{\theta,\alpha,\nu}(w)\ud w$ for any functions $h_1,h_2\in L_2(f_{\theta,\alpha,\nu})$.

The following lemma is important for our spectral analysis.

Let $K_{\alpha,\nu}(S_n)$ be the $n\times n$ covariance matrix whose $(i,j)$-entry is $K_{\alpha,\nu}(\bfs_i-\bfs_j)$. Then for any $\alpha \in \RR_+$, there exists an $n\times n$ invertible matrix $U_{\alpha}$ that depends on $\alpha,\alpha_0,\theta_0,\nu,S_n$, such that

\begin{align}\label{diagonalize}
& \theta_0  U_{\alpha}^\T  K_{\alpha_0,\nu}(S_n) U_{\alpha}= I_n, \qquad \theta_0 U_{\alpha}^\T  K_{\alpha,\nu}(S_n) U_{\alpha}= \Lambda_n(\alpha) = \diag\{\lambda_{k,n}(\alpha):k=1,\ldots,n\},
\end{align}

where $\{\lambda_{k,n}(\alpha),k=1,\ldots,n\}$ are the positive diagonal entries of the diagonal matrix $\Lambda_n(\alpha)$.

Furthermore, there exist orthonormal basis functions $\psi_1,\ldots,\psi_n \in L_2(f_{\theta_0,\alpha_0,\nu})$, such that for any $j,k\in \{1,\ldots,n\}$,

\begin{align}\label{wangloh.eq13}
& \langle \psi_j,\psi_k \rangle_{f_{\theta_0,\alpha_0,\nu}} =  \Ical(j=k), \qquad \langle \psi_j,\psi_k \rangle_{f_{\theta,\alpha,\nu}} =  \lambda_{j,n}(\alpha) \Ical(j=k),
\end{align}

where $\Ical(\cdot)$ is the indicator function.

The existence of such an invertible $U_{\alpha}$ is guaranteed by Theorem 7.6.4 and Corollary 7.6.5 on page 465–466 of 18. For completeness, we directly prove the existence of such an invertible matrix in the following general claim.

Claim: Suppose that $A$ and $B$ are two generic $n\times n$ symmetric positive definite matrices. Then there always exists an invertible matrix $U$, such that

\begin{align} \label{eq:UAU}
& U^\T  A U = I_n,\qquad U^\T  B U = \Lambda,
\end{align}

where $I_n$ is the $n\times n$ identity matrix and $\Lambda$ is an $n\times n$ diagonal matrix whose diagonal entries are all positive.

Proof of the Claim:
Since $B$ is symmetric positive definite, let $B=LL^\T $ be the Cholesky decomposition of $B$, where $L$ is an $n\times n$ lower triangular matrix with all positive diagonal entries and $L$ is invertible. Let $G=L^{-1} A L^{-\T}$. Then obviously $G$ is also a symmetric positive definite matrix with $G^\T =G$. Suppose that $G$ has the spectral decomposition $G=PDP^{-1}$ where $P$ is an $n\times n$ orthogonal matrix ($P^{-1}=P^\T $) and $D$ is a $n\times n$ diagonal matrix whose diagonal entries are all eigenvalues of $G$ and they are all positive. Then $P^\T  G P = D$. We let $U=L^{-\T} P D^{-1/2}$. It follows that

\begin{align*}
U^\T  A U &= D^{-1/2} P^\T  L^{-1} A L^{-\T} P D^{-1/2} \\
&= D^{-1/2} P^\T  G P D^{-1/2} = D^{-1/2} D D^{-1/2} = I_n,\\
U^\T  B U &= D^{-1/2} P^\T  L^{-1} B L^{-\T} P D^{-1/2} \\
&= D^{-1/2} P^\T  L^{-1} LL^\T  L^{-\T} P D^{-1/2} = D^{-1/2} P^\T  P  D^{-1/2} = D^{-1}.
\end{align*}

We set $\Lambda=D^{-1}$ which is an $n\times n$ diagonal matrix whose diagonal entries are all positive. This proves the claim.

Based on the claim, if we set $A=\theta_0 K_{\alpha_0,\nu}(S_n)$ and $B=\theta_0 K_{\alpha,\nu}(S_n)$, then we can find an invertible matrix $U$ such that (<ref>) holds. Because $\theta_0,\alpha_0,\nu$ are assumed to be fixed numbers, we can see that $U$ only changes with $\alpha$ and therefore we can write it as $U_{\alpha}$. Similarly, we can write $\Lambda_n(\alpha)$ to highlight its dependence on $\alpha$ and $n$. Correspondingly, we have $\theta_0  U_{\alpha}^\T  K_{\alpha_0,\nu}(S_n) U_{\alpha}= I_n$ and $\theta_0 U_{\alpha}^\T  K_{\alpha,\nu}(S_n) U_{\alpha}= \diag\{\lambda_{k,n}(\alpha):k=1,\ldots,n\} =\Lambda_n(\alpha)$.  This proves (<ref>). The existence of orthonormal basis functions $\psi_1,\ldots,\psi_n \in L_2(f_{\theta_0,\alpha_0,\nu})$ is proved in Section 4 of 47, which does not involve the specific forms of covariance functions. This completes the proof of Lemma <ref>.

Suppose that Assumption <ref> holds. Then for $\{\lambda_{k,n}(\alpha):k=1,\ldots,n\}$ defined in Lemma <ref> and any $\alpha\in \RR_+$,

\begin{align}
& \sup_{|\alpha/\alpha_0 -1|\leq r_0} \lambda_{k,n}(\alpha) \leq 1 + L|\alpha/\alpha_0-1|^{\kappa}, \label{lambda.upper1} \\
& \inf_{|\alpha/\alpha_0 -1|\leq r_0}\lambda_{k,n}(\alpha) \geq 1 - L|\alpha/\alpha_0-1|^{\kappa}, \label{lambda.lower1}
\end{align}

for all $k=1,\ldots,n$. Furthermore, if $\alpha\leq \alpha_0$, then $\lambda_{k,n}(\alpha)\geq 1$ for all $k=1,\ldots,n$; if  $\alpha\geq \alpha_0$, then $0\leq \lambda_{k,n}(\alpha)\leq 1$ for all $k=1,\ldots,n$.

For (<ref>) and (<ref>), we use the relation

$$\lambda_{k,n}(\alpha)= \int_{\RR^d} |\psi_k(w)|^2 f_{\theta_0,\alpha_0,\nu}(w) \cdot \frac{f_{\theta_0,\alpha,\nu}(w)}{f_{\theta_0,\alpha_0,\nu}(w)} \ud w$$

for $k=1,\ldots,n$ together with the bounds in Assumption <ref> (i) to obtain that for all $\alpha$ satisfying $|\alpha/\alpha_0 -1|\leq r_0$,

\begin{align}\label{lambda30}
\lambda_{k,n}(\alpha)&\leq \sup_{w \in \RR^d} \frac{f_{\theta_0,\alpha,\nu}(w)}{f_{\theta_0,\alpha_0,\nu}(w)}  \cdot \int_{\RR^d} |\psi_k(w)|^2 f_{\theta_0,\alpha_0,\nu}(w) \ud w \leq 1 + L|\alpha/\alpha_0-1|^{\kappa}, \nonumber \\
\lambda_{k,n}(\alpha)&\geq \inf_{w \in \RR^d} \frac{f_{\theta_0,\alpha,\nu}(w)}{f_{\theta_0,\alpha_0,\nu}(w)}  \cdot \int_{\RR^d} |\psi_k(w)|^2 f_{\theta_0,\alpha_0,\nu}(w) \ud w \geq 1 - L|\alpha/\alpha_0-1|^{\kappa}.
\end{align}

This proves (<ref>) and (<ref>).

By the non-increasing property of $f_{\theta,\alpha,\nu}(w)$ in $\alpha$ in Assumption <ref> (ii), we have that if $\alpha\leq \alpha_0$, $f_{\theta_0,\alpha,\nu}(w)\geq f_{\theta_0,\alpha_0,\nu}(w) $ for all $w\in \RR^d$, and for all $k=1,\ldots,n$,

\begin{align*}
\lambda_{k,n}(\alpha)& = \int_{\RR^d} |\psi_k(w)|^2 f_{\theta_0,\alpha,\nu}(w) \ud w \geq \int_{\RR^d} |\psi_k(w)|^2 f_{\theta_0,\alpha_0,\nu}(w) \ud w = 1.
\end{align*}

The conclusion for $\alpha\geq \alpha_0$ follows similarly.

Suppose that Assumption <ref> (ii) holds. For any $\alpha,\alpha'\in \RR_+$ such that $\alpha\leq \alpha'$, we have that $K_{\alpha,\nu}(S_n)\succeq K_{\alpha',\nu}(S_n)$, for any set $S_n=\{s_1,\ldots,s_n\}$ of $n$ distinct locations in $[0,1]^d$.

The proof of Lemma <ref> is motivated by the proof of Lemma 1 in 23. We define the matrix $\Omega^{\dagger} = K_{\alpha,\nu}(S_n) - K_{\alpha',\nu}(S_n)$. Then the entries of $\Omega^{\dagger}$ can be expressed in terms of a function $\widetilde K_{\Omega^{\dagger}}: \RR^d\to \RR$, with its $(i,j)$-entry

\begin{align*}
\Omega^{\dagger}_{ij} &= \widetilde K_{\Omega^{\dagger}}(\bfs_i-\bfs_j) = K_{\alpha,\nu}(\bfs_i-\bfs_j) - K_{\alpha',\nu}(\bfs_i-\bfs_j),
\end{align*}

for $i,j\in \{1,\ldots,n\}$. The matrix $\Omega^{\dagger}$ is positive definite if $\widetilde K_{\Omega^{\dagger}}$ is a positive semidefinite function.

By the non-increasing property of $f_{\theta,\alpha,\nu}(w)$ in $\alpha$ for any given $\theta>0$ and $w\in \RR^d$ as in Assumption <ref> (ii), we have that for $\alpha\leq \alpha'$,

\begin{align} \label{eq:f2.diff}
f_{\theta,\alpha,\nu}(w) & \geq f_{\theta,\alpha',\nu}(w).
\end{align}

Therefore, we can compute the spectral density of the function $\widetilde K_{\Omega^{\dagger}}$:

\begin{align} \label{eq:fOmega1}
f_{\Omega^{\dagger}}(w) &= \frac{1}{(2\pi)^d} \int_{\RR^d} \exp\left(-\imath w^\T x \right) \widetilde K_{\Omega^{\dagger}} (x) \ud x \nonumber \\
&= f_{1,\alpha,\nu}(w) - f_{1,\alpha',\nu}(w) \geq 0, \text{ for all } w \in \RR^d,
\end{align}

where the last step follows from (<ref>). This has shown that $\widetilde K_{\Omega^{\dagger}}$ is indeed a positive semidefinite function. Therefore, $\Omega^{\dagger}$ is a positive semidefinite matrix and the conclusion follows.

Suppose that Assumption <ref> holds. Let $\Acal_{n}^{\dagger}=\left\{\alpha\in \RR_+:|\alpha/\alpha_0-1|\leq n^{-2/\kappa} \right\}$ where $\kappa$ is as defined in Assumption <ref>. There exists a large integer $N_2$ that only depends on $\nu,d,\alpha_0,L,r_0,\kappa$, such that for all $n>N_2$,

\begin{align}\label{eq:lambda.3bound}
& \sup_{\alpha \in \Acal_{n}^{\dagger}}  \left\{ \sum_{k=1}^n \left| \lambda_{k,n}(\alpha) - 1 \right| \right\} \leq L n^{-1}, \qquad \sup_{\alpha \in \Acal_{n}^{\dagger}}  \left\{ \sum_{k=1}^n \left| \lambda_{k,n}(\alpha) - 1 \right|^2 \right\} \leq L^2 n^{-3}, \nonumber \\
& \sup_{\alpha \in \Acal_{n}^{\dagger}}  \left\{ \sum_{k=1}^n \left| \lambda_{k,n}(\alpha)^{-1} - 1 \right| \right\} \leq 2L n^{-1}, \qquad \sup_{\alpha \in \Acal_{n}^{\dagger}}  \left\{ \sum_{k=1}^n \left| \lambda_{k,n}(\alpha)^{-1} - 1 \right|^2 \right\} \leq 4L^2 n^{-3}.
\end{align}

For all sufficiently large $n$, $n^{-2/\kappa} < r_0$ for $r_0$ in Assumption <ref>. By Lemma <ref>, we have that for any $k=1,\ldots,n$, for all $\alpha\in \Acal_{n}^{\dagger}$,

\begin{align}
& \quad~ \left|\lambda_{k,n}(\alpha) - 1 \right| \leq L|\alpha/\alpha_0-1|^{\kappa} \nonumber.
\end{align}

Therefore, uniformly over all $\alpha\in \Acal_{n}^{\dagger}$,

\begin{align} \label{eq:lambda.2sum}
& \sum_{k=1}^n \left| \lambda_{k,n}(\alpha) - 1 \right|  \leq n\cdot L|\alpha/\alpha_0-1|^{\kappa} \leq Ln^{-1}, \nonumber \\
& \sum_{k=1}^n \left| \lambda_{k,n}(\alpha) - 1 \right|^2  \leq n\cdot L^2 |\alpha/\alpha_0-1|^{2\kappa} \leq L^2 n^{-3},
\end{align}

which has proved the first two inequalities in (<ref>).

For all $\alpha \in \Acal_{n}^{\dagger}$, for all sufficiently large $n$,

$$1-L|\alpha/\alpha_0-1|^{\kappa} \geq 1 - L n^{-2} \geq \frac{1}{2}. $$

Therefore, using (<ref>) in Lemma <ref> and (<ref>), we have that uniformly for all $\alpha \in \Acal_{n}^{\dagger}$, for all sufficiently large $n$,

\begin{align*}
& \min_{1\leq k\leq n} \lambda_{k,n}(\alpha) \geq 1/2, \quad \max_{1\leq k\leq n} \lambda_{k,n}(\alpha)^{-1}  \leq 2, \nonumber \\
& \sum_{k=1}^n \left| \lambda_{k,n}(\alpha)^{-1} - 1 \right| = \frac{\sum_{k=1}^n \left| \lambda_{k,n}(\alpha) - 1 \right|}{\min_{1\leq k\leq n} \lambda_{k,n}(\alpha)} \leq 2Ln^{-1}, \nonumber \\
& \sum_{k=1}^n \left| \lambda_{k,n}(\alpha)^{-1} - 1 \right|^2 = \frac{\sum_{k=1}^n \left| \lambda_{k,n}(\alpha) - 1 \right|}{\min_{1\leq k\leq n} \lambda_{k,n}(\alpha)^2} \leq 4L^2n^{-3}.
\end{align*}

which proves the last two inequalities in (<ref>).

§.§ Technical Lemmas for Evidence Lower Bound

Let $\widetilde Y_n = Y_n -  F_n  \beta_0$. Then $\widetilde Y_n \sim \Ncal(0,\theta_0 K_{\alpha_0,\nu}(S_n) + \tau_0 I_n)$, and $Y_n -  F_n  \beta = \widetilde Y_n -  F_n (\beta-\beta_0)$.

Suppose that Assumptions <ref> and <ref> hold. Let $\Acal_{n} = \big\{(\theta,\alpha,\tau,\beta)\in \RR_+^{3} \times \RR^p:$ $|\alpha/\alpha_0-1|\leq n^{-2/\kappa} \big\}$. There exists a large integer $N_3$ that only depends on $\nu,d,\theta_0,\alpha_0,L,r_0,\kappa$, such that for any $\epsilon_1,\epsilon_2\in (0,1/2)$, for all $n>N_3$, with probability at least $1-2\exp(-\log^2 n)$,

\begin{align} \label{eq:denom.lower12}
& \inf_{\Bcal_0(\epsilon_1,\epsilon_2)\cap \Acal_{n}} \Bigg\{ -\frac{1}{2} \left(\tr\left[\left\{\theta K_{\alpha,\nu}(S_n) + \tau I_n\right\}^{-1} \widetilde Y_n \widetilde Y_n^\T  \right] - \tr\left[\left\{\theta_0 K_{\alpha_0,\nu}(S_n) + \tau_0 I_n\right\}^{-1} \widetilde Y_n \widetilde Y_n^\T  \right] \right) \nonumber \\
&\qquad - \frac{1}{2} \log \frac{\dett\left\{\theta K_{\alpha,\nu}(S_n) +\tau I_n\right\}}{\dett\left\{\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right\}} \Bigg\} \nonumber \\
& \geq - \left(\frac{5}{2}n\epsilon_{1 \vee 2}  + 5Ln^{-1} \right),
\end{align}

where $\epsilon_{1\vee 2} = \max(\epsilon_1,\epsilon_2)$.

There are two terms inside the infimum in (<ref>). We provide lower bounds for each of them. For notational simplicity, we define

\begin{align} \label{eq:Omega}
\Omega_{\alpha} = \theta_0 K_{\alpha,\nu}(S_n) + \tau_0 I_n,\quad \Omega_0 = \theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n .
\end{align}

Lower bound for the first term in (<ref>):

In the set $\Bcal_0(\epsilon_1,\epsilon_2)$, $\theta>(1-\epsilon_{1 \vee 2})\theta_0$ and $\tau>(1-\epsilon_{1 \vee 2})\tau_0$. Since $K_{\alpha,\nu}(S_n)$ is positive definite, we have that $\theta K_{\alpha,\nu}(S_n) + \tau I_n \succeq (1-\epsilon_{1 \vee 2})\left\{\theta_0 K_{\alpha,\nu}(S_n) + \tau_0 I_n \right\}$, and hence $\left\{\theta K_{\alpha,\nu}(S_n) + \tau I_n\right\}^{-1} \preceq (1-\epsilon_{1 \vee 2})^{-1}\left\{\theta_0 K_{\alpha,\nu}(S_n) + \tau_0 I_n \right\}^{-1}$. Therefore, for $\widetilde Y_n = Y_n -  F_n  \beta_0 \sim \Ncal(0,\theta_0 K_{\alpha_0,\nu}(S_n) + \tau_0 I_n)$, by Lemma <ref> Part (ii), we have that

\begin{align} \label{eq:denom.expo1.1}
& \tr\left[\left\{\theta K_{\alpha,\nu}(S_n) + \tau I_n\right\}^{-1} \widetilde Y_n \widetilde Y_n^\T  \right] - \tr\left[\left\{\theta_0 K_{\alpha_0,\nu}(S_n) + \tau_0 I_n\right\}^{-1} \widetilde Y_n \widetilde Y_n^\T  \right] \nonumber \\
\leq{}& (1-\epsilon_{1 \vee 2})^{-1} \tr\left[\left\{\theta_0 K_{\alpha,\nu}(S_n) + \tau_0 I_n \right\}^{-1} \widetilde Y_n \widetilde Y_n^\T  \right] - \tr\left[\left\{\theta_0 K_{\alpha_0,\nu}(S_n) + \tau_0 I_n\right\}^{-1} \widetilde Y_n \widetilde Y_n^\T  \right] \nonumber \\
={}& (1-\epsilon_{1 \vee 2})^{-1} \left\{ \tr\left(\Omega_{\alpha}^{-1} \widetilde Y_n \widetilde Y_n^\T  \right) - \tr\left(\Omega_0^{-1} \widetilde Y_n \widetilde Y_n^\T  \right) \right\}  + \frac{\epsilon_{1 \vee 2}}{1-\epsilon_{1 \vee 2}} \tr\left(\Omega_0^{-1} \widetilde Y_n \widetilde Y_n^\T  \right).
\end{align}

We derive an upper bound for the first term inside the bracket in (<ref>). Using (<ref>) in Lemma <ref>, we have that $\theta_0 K_{\alpha,\nu}(S_n) =U_{\alpha}^{-\T} \Lambda_n(\alpha) U_{\alpha}^{-1}$ and $\theta_0 K_{\alpha_0,\nu}(S_n) =U_{\alpha}^{-\T} U_{\alpha}^{-1}$. Therefore, we repeatedly apply the Sherman-Morrison-Woodbury formula for inverse matrices to obtain that

\begin{align} \label{eq:denom.expo1.2}
& ~~\Omega_{\alpha}^{-1} - \Omega_0^{-1} \nonumber \\
={}& \left(U_{\alpha}^{-\T} \Lambda_n(\alpha) U_{\alpha}^{-1} + \tau_0 I_n\right)^{-1} - \Omega_0^{-1} \nonumber \\
={}& \left\{U_{\alpha}^{-\T} \left(\Lambda_n(\alpha)-I_n\right) U_{\alpha}^{-1} + \Omega_0 \right\}^{-1} - \Omega_0^{-1} \nonumber \\
={}& \Omega_0^{-1} - \Omega_0^{-1}U_{\alpha}^{-\T} \left[\left\{\Lambda_n(\alpha) - I_n\right\}^{-1} + U_{\alpha}^{-1} \Omega_0^{-1} U_{\alpha}^{-\T} \right]^{-1} U_{\alpha}^{-1} \Omega_0^{-1}  - \Omega_0^{-1} \nonumber \\
={}& \Omega_0^{-1}U_{\alpha}^{-\T} \left[\left\{I_n - \Lambda_n(\alpha) \right\}^{-1} - U_{\alpha}^{-1} \Omega_0^{-1} U_{\alpha}^{-\T} \right]^{-1} U_{\alpha}^{-1} \Omega_0^{-1}.
\end{align}

Next we focus on the matrix $T_{\alpha}\equiv \left[\left\{I_n - \Lambda_n(\alpha) \right\}^{-1} - U_{\alpha}^{-1} \Omega_0^{-1} U_{\alpha}^{-\T} \right]^{-1}$. Obviously $U_{\alpha}^{-1} \Omega_0^{-1} U_{\alpha}^{-\T}$ is positive definite. Using the definition of $\Omega_0$ in (<ref>), it follows from Lemma <ref> Part (i) that

\begin{align}\label{eq:UOU}
0_{n\times n} \prec U_{\alpha}^{-1} \Omega_0^{-1} U_{\alpha}^{-\T} & = U_{\alpha}^{-1} \left(U_{\alpha}^{-\T}U_{\alpha}^{-1} +\tau_0 I_n\right)^{-1} U_{\alpha}^{-\T} \nonumber \\
& \preceq  U_{\alpha}^{-1} \left(U_{\alpha}^{-\T} U_{\alpha}^{-1}\right)^{-1} U_{\alpha}^{-\T} = I_n.
\end{align}

Therefore, all the eigenvalues of $U_{\alpha}^{-1} \Omega_0^{-1} U_{\alpha}^{-\T}$ are between 0 and 1. Meanwhile, if $\alpha>\alpha_0$, then by Lemma <ref>, $K_{\alpha,\nu}(S_n)\preceq K_{\alpha_0,\nu}(S_n)$, and hence $\Omega_{\alpha}^{-1} \succeq \Omega_0^{-1}$. From (<ref>), this implies that $T_{\alpha}$ is positive definite when $\alpha> \alpha_0$. Similarly, $T_{\alpha}$ is negative definite when $\alpha< \alpha_0$. Let $T_{\alpha} = P_{\alpha}^{-1} \diag\{a_1,\ldots,a_n\} P_{\alpha}$ be the eigendecomposition of $T_{\alpha}$, where $P_{\alpha}$ is an $n\times n$ orthogonal matrix consisting of eigenvectors, and $a_1,\ldots,a_n$ are the eigenvalues of $T_{\alpha}$. Then $a_1^{-1},\ldots,a_n^{-1}$ are the eigenvalues of $T_{\alpha}^{-1}$. When $\alpha>\alpha_0$, since $T_{\alpha}^{-1}, \left\{I_n - \Lambda_n(\alpha) \right\}^{-1}, U_{\alpha}^{-1} \Omega_0^{-1} U_{\alpha}^{-\T}$ are all positive definite, by applying the Weyl's inequality to the relation $T_{\alpha}^{-1} = \left\{I_n - \Lambda_n(\alpha) \right\}^{-1}- U_{\alpha}^{-1} \Omega_0^{-1} U_{\alpha}^{-\T}$, we have that up to a permutation of $a_k$'s, $\{1-\lambda_{k,n}(\alpha)\}^{-1}-1 \leq a_k^{-1} \leq \{1-\lambda_{k,n}(\alpha)\}^{-1}$ for $k=1,\ldots,n$. Similarly, when $\alpha<\alpha_0$, $T_{\alpha}^{-1}$ is negative definite, and $\left\{I_n - \Lambda_n(\alpha) \right\}^{-1}$ is also negative definite since all $\lambda_{k,n}(\alpha)\geq 1$ by Lemma <ref>. We can apply the Weyl's inequality to the relation $(-T_{\alpha})^{-1} = \left\{\Lambda_n(\alpha) -I_n\right\}^{-1} + U_{\alpha}^{-1} \Omega_0^{-1} U_{\alpha}^{-\T}$ to obtain that $\{\lambda_{k,n}(\alpha)-1\}^{-1} \leq -a_k^{-1} \leq \{\lambda_{k,n}(\alpha)-1\}^{-1}+1$, or equivalently $\{1-\lambda_{k,n}(\alpha)\}^{-1}-1 \leq a_k^{-1} \leq \{1-\lambda_{k,n}(\alpha)\}^{-1} $ for $k=1,\ldots,n$.

On the other hand, from Lemma <ref>, we have that uniformly over all $\alpha\in \Acal_{n}$, for all sufficiently large $n$, $\max_{1\leq k\leq n}|\lambda_{k,n}(\alpha)-1|\leq Ln^{-2} < 1/2$, i.e., $\{1- \lambda_{k,n}(\alpha) \}^{-1} \in [-\infty,-2] \cup [2,+\infty]$ for all $k=1,\ldots,n$. Therefore, we have that for both the cases of  $\alpha>\alpha_0$ and $\alpha<\alpha_0$, the eigenvalues of $T_{\alpha}$ satisfy

\begin{align} \label{eq:lambda-1.bound}
& |a_k| \leq \max\left\{\left|\frac{1}{\{1-\lambda_{k,n}(\alpha)\}^{-1} }\right|,  \left|\frac{1}{\{1-\lambda_{k,n}(\alpha)\}^{-1}-1 }\right| \right\} \nonumber \\
&= \max\left\{|\lambda_{k,n}(\alpha)-1|, \left|\lambda_{k,n}(\alpha)^{-1} -1 \right| \right\} \leq |\lambda_{k,n}(\alpha)-1| + \left|\lambda_{k,n}(\alpha)^{-1} -1 \right|,
\end{align}

for $k=1,\ldots,n$. This also holds for the case $\alpha=\alpha_0$ since both sides of (<ref>) are equal to zero.

We define the $n\times n$ diagonal matrix $Q_n(\alpha)$:

\begin{align} \label{eq:Qn}
Q_n(\alpha) & = \diag\left\{|\lambda_{1,n}(\alpha)-1| + \left|\lambda_{1,n}(\alpha)^{-1} -1 \right|, \ldots, |\lambda_{n,n}(\alpha)-1| + \left|\lambda_{n,n}(\alpha)^{-1} -1 \right| \right\}.
\end{align}

Then by Lemma <ref> Part (i), (<ref>) and the eigendecomposition $T_{\alpha} = P_{\alpha}^{-1} \diag\{a_1,\ldots,a_n\} P_{\alpha}$ imply that

\begin{align}
T_{\alpha} &= \left[\left\{I_n - \Lambda_n(\alpha) \right\}^{-1} - U_{\alpha}^{-1} \Omega_0^{-1} U_{\alpha}^{-\T} \right]^{-1} \preceq P_{\alpha}^{-1}Q_n(\alpha)P_{\alpha}. \nonumber
\end{align}

This together with (<ref>) implies that for all $\alpha\in \Acal_{n}$,

\begin{align}
& \Omega_{\alpha}^{-1} - \Omega_0^{-1} \preceq \Omega_0^{-1}U_{\alpha}^{-\T}P_{\alpha}^{-1}Q_n(\alpha)P_{\alpha}U_{\alpha}^{-1} \Omega_0^{-1} .  \nonumber
\end{align}

Hence the first term in (<ref>) can be upper bounded by

\begin{align} \label{eq:denom.expo1.5}
& (1-\epsilon_{1 \vee 2})^{-1} \left[ \tr\left\{\left(\theta_0 K_{\alpha,\nu}(S_n) + \tau_0 I_n \right)^{-1} \widetilde Y_n \widetilde Y_n^\T  \right\} - \tr\left\{\left(\theta_0 K_{\alpha_0,\nu}(S_n) + \tau_0 I_n\right)^{-1} \widetilde Y_n \widetilde Y_n^\T  \right\} \right] \nonumber \\
\leq{}& (1-\epsilon_{1 \vee 2})^{-1} \widetilde Y_n^\T \Omega_0^{-1}U_{\alpha}^{-\T} P_{\alpha}^{-1}Q_n(\alpha)P_{\alpha} U_{\alpha}^{-1} \Omega_0^{-1} \widetilde Y_n.
\end{align}

Since $\Omega_0$ is positive definite, there exists an invertible matrix $V_0$ such that $V_0^\T \Omega_0 V_0 = I_n$ ($V_0$ depends on $\alpha_0$ and $\tau_0$) and hence $\Omega_0^{-1}=V_0V_0^\T$. Define the random vector $Z_n = V_0^{\T} \widetilde Y_n$. Since $\widetilde Y_n\sim \mathcal{N}(0,\Omega_0)$, we have that $Z_n \sim \mathcal{N}(0,I_n)$. Therefore, we can rewrite the right-hand side of (<ref>) as

\begin{align} \label{eq:denom.expo1.6}
& \widetilde Y_n^\T  \Omega_0^{-1}U_{\alpha}^{-\T} P_{\alpha}^{-1}Q_n(\alpha)P_{\alpha} \Omega_0^{-1} \widetilde Y_n
= Z_n^\T  V_0^\T  U_{\alpha}^{-\T} P_{\alpha}^{-1}Q_n(\alpha)P_{\alpha}  U_{\alpha}^{-1} V_0 Z_n.
\end{align}

Then we apply the Hanson-Wright inequality in Lemma <ref> to the quadratic form in (<ref>) to obtain that for any number $z>0$ and any $\alpha\in \Acal_{n}$,

\begin{align} \label{eq:denom.expo1.7}
& \pr \left\{Z_n^\T  \Sigma_{\alpha} Z_n \geq \tr(\Sigma_{\alpha}) + 2\left\{\tr(\Sigma_{\alpha}^2)z\right\}^{1/2} + 2\|\Sigma_{\alpha}\|_{\op}z \right\} \leq \ee^{-z},
\end{align}

where we define $\Sigma_{\alpha}=V_0^\T  U_{\alpha}^{-\T}  P_{\alpha}^{-1}Q_n(\alpha)P_{\alpha}  U_{\alpha}^{-1} V_0$.

We continue to find upper bounds for each term related to $\Sigma_{\alpha}$ in (<ref>), using the definition of $Q_n(\alpha)$ in (<ref>):

\begin{align}
\tr(\Sigma_{\alpha}) & = \tr\left\{V_0^\T  U_{\alpha}^{-\T}  P_{\alpha}^{-1}Q_n(\alpha)P_{\alpha}  U_{\alpha}^{-1} V_0 \right\} =  \tr\left\{ P_{\alpha}^{-1}Q_n(\alpha)P_{\alpha}  U_{\alpha}^{-1} V_0 V_0^\T  U_{\alpha}^{-\T} \right\} \nonumber \\
&= \tr\left\{ P_{\alpha}^{-1}Q_n(\alpha)P_{\alpha}  U_{\alpha}^{-1}\Omega_0^{-1} U_{\alpha}^{-\T} \right\} \stackrel{(i)}{\leq} \tr\left\{ P_{\alpha}^{-1}Q_n(\alpha)P_{\alpha} \right\} = \tr\left\{ Q_n(\alpha) \right\} \nonumber \\
&=  \sum_{k=1}^n |\lambda_1(\alpha)-1| + \sum_{k=1}^n \left|\lambda_1(\alpha)^{-1} -1 \right|. \label{eq:denom.expo1.8} \\
\tr(\Sigma_{\alpha}^2) & = \tr\left[\left\{V_0^\T  U_{\alpha}^{-\T} P_{\alpha}^{-1}Q_n(\alpha)P_{\alpha} U_{\alpha}^{-1} V_0\right\}^2 \right] \nonumber \\
&= \tr\left\{ V_0^\T  U_{\alpha}^{-\T} P_{\alpha}^{-1}Q_n(\alpha)P_{\alpha} U_{\alpha}^{-1} V_0 V_0^\T  U_{\alpha}^{-\T} P_{\alpha}^{-1}Q_n(\alpha)P_{\alpha} U_{\alpha}^{-1} V_0 \right\} \nonumber \\
&= \tr\left\{P_{\alpha}^{-1}Q_n(\alpha)P_{\alpha} U_{\alpha}^{-1} \Omega_0^{-1}  U_{\alpha}^{-\T} P_{\alpha}^{-1}Q_n(\alpha)P_{\alpha} U_{\alpha}^{-1} \Omega_0^{-1}  U_{\alpha}^{-\T} \right\} \nonumber \\
&\stackrel{(ii)}{\leq}  \tr\left\{P_{\alpha}^{-1}Q_n(\alpha)P_{\alpha} U_{\alpha}^{-1} \Omega_0^{-1}  U_{\alpha}^{-\T} P_{\alpha}^{-1}Q_n(\alpha)P_{\alpha} \right\} \nonumber \\
&=\tr\left\{P_{\alpha}^{-1}Q_n(\alpha)^2P_{\alpha} U_{\alpha}^{-1} \Omega_0^{-1}  U_{\alpha}^{-\T}  \right\} \stackrel{(iii)}{\leq} \tr\left\{P_{\alpha}^{-1}Q_n(\alpha)^2P_{\alpha}\right\} = \tr\left\{Q_n(\alpha)^2\right\} \nonumber \\
&= \sum_{k=1}^n \left\{\left|\lambda_{k,n}(\alpha) -1\right| + \left|\lambda_{k,n}(\alpha)^{-1}-1 \right| \right\}^2 \nonumber \\
&\stackrel{(iv)}{\leq} 2\sum_{k=1}^n \left| \lambda_{k,n}(\alpha) -1 \right|^2 + 2 \sum_{k=1}^n \left|\lambda_{k,n}(\alpha)^{-1}-1\right|^2 , \label{eq:denom.expo1.9} \\
\|\Sigma_{\alpha}\|_{\op}^2 & = \mathsf{s}_{\max} \left[ \left\{ V_0^\T  U_{\alpha}^{-\T} P_{\alpha}^{-1}Q_n(\alpha)P_{\alpha} U_{\alpha}^{-1} V_0\right\}^2 \right] \leq \tr\left[\left\{V_0^\T  U_{\alpha}^{-\T} P_{\alpha}^{-1}Q_n(\alpha)P_{\alpha} U_{\alpha}^{-1} V_0\right\}^2 \right]   \nonumber \\
&\leq 2\sum_{k=1}^n \left| \lambda_{k,n}(\alpha) -1 \right|^2 + 2 \sum_{k=1}^n \left|\lambda_{k,n}(\alpha)^{-1}-1\right|^2 , \label{eq:denom.expo1.10}
\end{align}

where (i), (ii) and (iii) follow from Lemma <ref> and (<ref>), (iv) follows from the inequality $(a+b)^2\leq 2(a^2+b^2)$ for any $a,b\in \RR$.

If we take $z=(\log n)^2$, using the upper bounds in Lemma <ref> (noticing that $\Acal_{n}^{\dagger}$ and $\Acal_{n}$ have the same condition on $\alpha$), we can obtain from (<ref>)-(<ref>) that for all sufficiently large $n$,

\begin{align} \label{eq:denom.expo1.11}
& \sup_{\alpha \in \Acal_{n}} \left\{ \tr(\Sigma_{\alpha}) + 2\left\{\tr(\Sigma_{\alpha}^2)\right\}^{1/2}\log n + 2\|\Sigma_{\alpha}\|_{\op} \log^2 n \right\}  \nonumber \\
& \leq Ln^{-1} + 2Ln^{-1} + 4 \left(L^2n^{-3} + 4L^2n^{-3} \right)^{1/2} \log n + 4 \left(L^2n^{-3} + 4L^2n^{-3} \right)^{1/2} \log^2 n \leq 4Ln^{-1}.
\end{align}

Because the matrix $K_{\alpha,\nu}(S_n)$ is non-increasing in $\alpha$ by Lemma <ref>, we can see from the definition of $\Omega_{\alpha}$ in (<ref>) that the first term in (<ref>) is an increasing function in $\alpha$. Let $\overline\alpha_n= \alpha_0 ( 1 + n^{-2/\kappa} ) $. Then the first term in (<ref>) attains its supremum at $\overline\alpha_n$. We combine (<ref>), (<ref>), (<ref>) and (<ref>) to obtain that

\begin{align} \label{eq:denom.expo1.12}
& \pr \Bigg\{ \sup_{\Acal_{n} }(1-\epsilon_{1 \vee 2})^{-1} \Big[ \tr\left\{\left(\theta_0 K_{\alpha,\nu}(S_n) + \tau_0 I_n \right)^{-1} \widetilde Y_n \widetilde Y_n^\T  \right\} \nonumber \\
&\quad - \tr\left\{\left(\theta_0 K_{\alpha_0,\nu}(S_n) + \tau_0 I_n\right)^{-1} \widetilde Y_n \widetilde Y_n^\T  \right\} \Big] \geq (1-\epsilon_{1 \vee 2})^{-1} 4Ln^{-1} \Bigg\} \nonumber \\
={}& \pr \Bigg\{ (1-\epsilon_{1 \vee 2})^{-1} \Big( \tr\left[\left\{\theta_0 K_{\overline\alpha_n,\nu}(S_n) + \tau_0 I_n \right\}^{-1} \widetilde Y_n \widetilde Y_n^\T  \right] \nonumber \\
&\quad - \tr\left\{\left(\theta_0 K_{\alpha_0,\nu}(S_n) + \tau_0 I_n\right)^{-1} \widetilde Y_n \widetilde Y_n^\T  \right\} \Big) \geq (1-\epsilon_{1 \vee 2})^{-1} 4Ln^{-1}  \Bigg\} \nonumber \\
\leq{}& \pr \Bigg\{  \widetilde Y_n^\T   \Omega_0^{-1}U_{\overline\alpha_n}^{-\T} P_{\overline\alpha_n}^{-1}Q_n(\overline\alpha_n)P_{\overline\alpha_n} U_{\overline\alpha_n}^{-1} \Omega_0^{-1} \widetilde Y_n \geq 4Ln^{-1}  \Bigg\} \nonumber \\
={}& \pr \Bigg\{ Z_n^\T \Sigma_{\overline\alpha_n} Z_n \geq 4Ln^{-1} \Bigg\} \nonumber \\
\leq{}& \pr \Bigg\{ Z_n^\T  \Sigma_{\overline\alpha_n} Z_n \geq  \tr(\Sigma_{\overline\alpha_n}) + 2\left\{\tr(\Sigma_{\overline\alpha_n}^2)\right\}^{1/2}\log n + 2\|\Sigma_{\overline\alpha_n}\|_{\op} \log^2 n \Bigg\} \nonumber \\
\leq{} & \exp(-\log^2 n).
\end{align}

For the second term on the right-hand side of (<ref>), with the relation $\tr\left(\Omega_0^{-1} \widetilde Y_n \widetilde Y_n^\T  \right) = Z_n^\T  Z_n$, we apply the Hanson-Wright inequality in Lemma <ref> with $\Sigma=I_n$ and $z=(\log n)^2$ to obtain that for all sufficiently large $n$,

\begin{align} \label{eq:denom.expo1.14}
&\quad~ \pr\left\{ \frac{\epsilon_{1 \vee 2}}{1-\epsilon_{1 \vee 2}}  \tr\left(\Omega_0^{-1} \widetilde Y_n \widetilde Y_n^\T  \right) \geq \frac{\epsilon_{1 \vee 2}}{1-\epsilon_{1 \vee 2}}  \cdot 2n \right\} \nonumber \\
&\leq \pr\left\{ Z_n^\T  Z_n \geq n + 2\sqrt{n} \log n + 2\log^2 n \right\} \leq \exp(-\log^2 n) .
\end{align}

Lower bound for the second term in (<ref>):

We first show a simple fact related to the Kullback-Leibler (KL) divergence. For two generic probability distributions $F_1,F_2$ with densities $f_1,f_2$, the KL divergence from $F_1$ to $F_2$ is defined as $\kl(F_1,F_2)=\int f_1\log (f_1/f_2)$. The KL divergence from $\mathcal{N}\left(0,\theta K_{\alpha,\nu}(S_n) +\tau I_n\right)$ to
$\mathcal{N}\left(0,\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right)$ is

\begin{align*}
&\kl\left(\mathcal{N}\left(0,\theta K_{\alpha,\nu}(S_n) +\tau I_n\right),\mathcal{N}\left(0,\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right)\right) \nonumber \\
={}& \frac{1}{2}\left[\log \frac{\dett\left(\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right)}{\dett\left(\theta K_{\alpha,\nu}(S_n) +\tau I_n\right)} + \tr\left\{\left(\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right)^{-1} \left(\theta K_{\alpha,\nu}(S_n) +\tau I_n\right) \right\} - n\right].
\end{align*}

Since the KL divergence is always nonnegative, this implies that

\begin{align} \label{eq:logdet.bound}
\log \frac{\dett\left\{\theta K_{\alpha,\nu}(S_n) +\tau I_n\right\}}{\dett\left\{\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right\}} \leq \tr\left[\left\{\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right\}^{-1} \left\{\theta K_{\alpha,\nu}(S_n) +\tau I_n\right\} \right] - n.
\end{align}

On the parameter set $\Bcal_0(\epsilon_1,\epsilon_2)$, $\theta\leq (1+\epsilon_1)\theta_0$ and $\tau \leq (1+\epsilon_2)\tau_0$. Therefore, we have $\theta K_{\alpha,\nu}(S_n) +\tau I_n\preceq (1+\epsilon_{1\vee 2}) \left(\theta_0 K_{\alpha,\nu}(S_n) +\tau_0 I_n\right)$.
By (<ref>) and Lemma <ref>, we have that

\begin{align} \label{eq:denom.expo2.1}
&\quad \log \frac{\dett\left\{\theta K_{\alpha,\nu}(S_n) +\tau I_n\right\}}{\dett\left\{\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right\}} \nonumber \\
& \leq \tr\left[\left\{\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right\}^{-1} \left\{\theta K_{\alpha,\nu}(S_n) +\tau I_n\right\} \right] - n \nonumber \\
&\leq \tr\left[\left\{\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right\}^{-1} \left\{(1+\epsilon_{1\vee 2})\theta_0 K_{\alpha,\nu}(S_n) + (1+\epsilon_{1\vee 2}) \tau_0 I_n\right\} \right] - n \nonumber \\
&= (1+\epsilon_{1\vee 2})\tr\left[\left\{\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right\}^{-1} \left\{\theta_0 K_{\alpha,\nu}(S_n) + \tau_0 I_n\right\} \right] - n \nonumber \\
&= (1+\epsilon_{1\vee 2})\theta_0 \tr\left[\left\{\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right\}^{-1} \left\{ K_{\alpha,\nu}(S_n) - K_{\alpha_0,\nu}(S_n)\right\} \right] \nonumber \\
&\quad +  (1+\epsilon_{1\vee 2})\tr\left[\left\{\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right\}^{-1} \left\{\theta_0 K_{\alpha_0,\nu}(S_n) + \tau_0 I_n\right\} \right]  - n \nonumber \\
&= (1+\epsilon_{1\vee 2})\theta_0 \tr\left[\left\{\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right\}^{-1} \left\{ K_{\alpha,\nu}(S_n) - K_{\alpha_0,\nu}(S_n)\right\} \right] + n \epsilon_{1\vee 2}.
\end{align}

We analyze the first term on the right-hand side of (<ref>). By Lemma <ref> and Lemma <ref>, for all $|\alpha/\alpha_0-1|\leq n^{-2/\kappa}$, for all sufficiently large $n$,

\begin{align*}
& \theta_0\left\{K_{\alpha,\nu}(S_n) - K_{\alpha_0,\nu}(S_n)\right\} = U_{\alpha}^{-\T}\left\{ \Lambda_n(\alpha)- I_n \right\} U_{\alpha}^{-1} \\
& \preceq \sup_{\Acal_{n}} \max_{1\leq k\leq n} | \lambda_{k,n}(\alpha)-1 | \cdot U_{\alpha}^{-\T} U_{\alpha}^{-1}
\preceq  (Ln^{-2}) U_{\alpha}^{-\T} U_{\alpha}^{-1}.
\end{align*}

\begin{align} \label{eq:denom.expo2.2}
& \theta_0\tr\left\{\left(\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right)^{-1} \left( K_{\alpha,\nu}(S_n) - K_{\alpha_0,\nu}(S_n)\right) \right\} \nonumber \\
\leq {}& \tr\left\{\left(\theta_0 K_{\alpha_0,\nu}(S_n) \right)^{-1} \left( \theta_0K_{\alpha,\nu}(S_n) - \theta_0 K_{\alpha_0,\nu}(S_n)\right) \right\} \nonumber \\
\leq{}& \tr\left\{ U_{\alpha} U_{\alpha}^\T  (Ln^{-2}) U_{\alpha}^{-\T} U_{\alpha}^{-1} \right\} \nonumber \\
={}& \tr\left( Ln^{-2} I_n \right) = Ln^{-1}.
\end{align}

We combine (<ref>) and (<ref>) to obtain that for all sufficiently large $n$,

\begin{align} \label{eq:denom.expo2.4}
& \quad \sup_{\Acal_{n}} \log \frac{\dett\left(\theta K_{\alpha,\nu}(S_n) +\tau I_n\right)}{\dett\left(\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right)} \leq (1+\epsilon_{1\vee 2}) Ln^{-1} + n \epsilon_{1\vee 2}.
\end{align}

Finally, we combine (<ref>), (<ref>), (<ref>) and (<ref>) to obtain that with probability for all sufficiently large $n$, at least $1-2\exp(-\log^2 n)$,

\begin{align}
&\quad \inf_{\Bcal_0(\epsilon_1,\epsilon_2)\cap \Acal_{n}} \Bigg\{ -\frac{1}{2} \Big(\tr\left[\left\{\theta K_{\alpha,\nu}(S_n) + \tau I_n\right\}^{-1} \widetilde Y_n \widetilde Y_n^\T  \right]  \nonumber \\
&\qquad - \tr\left[\left\{\theta_0 K_{\alpha_0,\nu}(S_n) + \tau_0 I_n\right\}^{-1} \widetilde Y_n \widetilde Y_n^\T  \right] \Big) - \frac{1}{2} \log \frac{\dett\left\{\theta K_{\alpha,\nu}(S_n) +\tau I_n\right\}}{\dett\left\{\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right\}} \Bigg\} \nonumber \\
&\geq  -\frac{1}{2}\left( \frac{4Ln^{-1}}{1-\epsilon_{1 \vee 2}} +  \frac{2n\epsilon_{1 \vee 2}}{1-\epsilon_{1 \vee 2}} + (1+\epsilon_{1\vee 2} ) Ln^{-1} +  n \epsilon_{1\vee 2} \right) \nonumber \\
&\geq - \left(\frac{5}{2}n\epsilon_{1 \vee 2}  + 5Ln^{-1} \right),  \nonumber
\end{align}

where the last inequality follows since $\epsilon_1,\epsilon_2\in (0,1/2)$.

Suppose that Assumptions <ref> and <ref> hold. Let $\Dcal_{n} = \big\{(\theta,\alpha,\tau,\beta)\in \RR_+^{3} \times \RR^p:$ $\|\beta-\beta_0\|\leq n^{-3} \big\}$. There exists a large integer $N_4$ that $\nu,d,\theta_0,\alpha_0,L,r_0,\kappa$, such that for any $\epsilon_1,\epsilon_2\in (0,1/2)$, for all $n>N_4$, with probability at least $1-\exp(-\log^2 n)$,

\begin{align} \label{eq:denom.lower34}
& \inf_{\Bcal_0(\epsilon_1,\epsilon_2)\cap \Dcal_{n}} \Big\{ \widetilde Y_n^\T  \left[\left\{\theta K_{\alpha,\nu}(S_n) +\tau I_n\right\}^{-1}-\left\{\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right\}^{-1}\right] F_n (\beta-\beta_0) \nonumber \\
&\quad - \frac{1}{2} (\beta-\beta_0)^\T   F_n ^\T  \left[\left\{\theta K_{\alpha,\nu}(S_n) +\tau I_n\right\}^{-1}-\left\{\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right\}^{-1}\right] F_n (\beta-\beta_0) \Big\} \nonumber \\
& \geq - 8C_{\ff}p^{1/2} \left\{\theta_0 K_{\alpha_0,\nu}(0) + \tau_0 \right\}^{1/2} \tau_0^{-1} n^{-1}- 2C_{\ff}^2 p \tau_0^{-1} n^{-4} .
\end{align}

We derive lower bounds for each of the two terms inside the infimum of (<ref>). On the event $\Dcal_{n}$, $\|\beta-\beta_0\|\leq n^{-3}$. By Assumption <ref>, we have that

\begin{align}\label{eq:F.bound1}
& \| F_n (\beta-\beta_0)\|  \leq \|F_n\|_F \|\beta-\beta_0\| \leq C_{\ff} p^{1/2} n^{1/2} \|\beta-\beta_0\| \leq C_{\ff} p^{1/2} n^{-2}.
\end{align}

On the set $\Bcal_0(\epsilon_1,\epsilon_2)$, since both matrices $\theta K_{\alpha,\nu}(S_n) +\tau I_n$ and $\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n$ are positive definite, we apply the triangle inequality for the matrix operator norm and obtain that

\begin{align}\label{eq:omega.diff1}
&\quad ~ \left\|\{\theta K_{\alpha,\nu}(S_n) +\tau I_n\}^{-1} - \{\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\}^{-1}\right\|_{\op} \nonumber \\
&\leq \left\|\{\theta K_{\alpha,\nu}(S_n) +\tau I_n\}^{-1}\right\|_{\op} + \left\|\{\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\}^{-1}\right\|_{\op} \nonumber \\
&\leq \frac{1}{\tau} + \frac{1}{\tau_0} \leq \frac{2}{(1-\epsilon_2)\tau_0}.
\end{align}

Since $\widetilde Y_n \sim \Ncal(0,\Omega_0)$, similar to the proof of Lemma <ref>, we let $V_0$ be an invertible matrix such that $V_0^\T  \Omega_0 V_0 = I_n$ and let $Z_n=V_0^\T  \widetilde Y_n$ such that $Z_n \sim \Ncal(0,I_n)$. Then by applying the Hanson-Wright inequality in Lemma <ref>, we have that for all sufficiently large $n$,

\begin{align} \label{eq:Ynorm}
& \Pr \left(\|V_0^\T  \widetilde Y_n\| \geq \sqrt{2n} \right)  = \Pr \left( \|Z_n\|^2 \geq 2n \right) \nonumber \\
&\leq \Pr \left( Z_n^\T  Z_n \geq n+ 2\sqrt{n} \log n + 2\log^2 n\right) \leq \exp(-\log^2 n).
\end{align}

Notice that since $V_0^\T  \Omega_0 V_0 = I_n$, we have that $V_0^{-\T}V_0^{-1} = \Omega_0$ and hence

\begin{align} \label{eq:omega0.op}
\|V_0^{-1}\|_{\op}^2 & \leq \|\Omega_0\|_{\op} \leq \left\|\theta_0 K_{\alpha_0,\nu}(S_n) + \tau_0 I_n \right\|_{\op} \nonumber \\
&\leq \left\|\theta_0 K_{\alpha_0,\nu}(S_n) \right\|_{\op} + \|\tau_0 I_n\|_{\op} \leq \theta_0 \tr\left\{K_{\alpha_0,\nu}(S_n)\right\} + \tau_0  \nonumber \\
&\leq n\theta_0 K_{\alpha_0,\nu}(0)  +\tau_0,
\end{align}

where $K_{\alpha_0,\nu}(0) \in \RR_+$ is a constant since $K_{\alpha_0,\nu}(\bfs-\bft)$ is a covariance function for $\bfs,\bft\in [0,1]^d$.

We combine (<ref>), (<ref>), (<ref>), and (<ref>) to conclude that with probability at least $1-\exp(-\log^2 n)$, on the event $\Bcal_0(\epsilon_1,\epsilon_2)\cap \Dcal_{n}$,

\begin{align}\label{eq:term3.lower1}
&~ \quad \widetilde Y_n^\T  \left[\left\{\theta K_{\alpha,\nu}(S_n) +\tau I_n\right\}^{-1}-\left\{\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right\}^{-1}\right] F_n (\beta-\beta_0) \nonumber \\
&\geq - \|\widetilde Y_n^\T  V_0\|\cdot \|V_0^{-1}\|_{\op} \cdot \left\|\left\{\theta K_{\alpha,\nu}(S_n) +\tau I_n\right\}^{-1}-\left\{\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right\}^{-1}\right\|_{\op} \cdot \| F_n (\beta-\beta_0)\| \nonumber \\
&\geq - \left[2n \cdot \left\{n\theta_0 K_{\alpha_0,\nu}(0) +\tau_0\right\}\right]^{1/2} \cdot \frac{2}{(1-\epsilon_2)\tau_0} \cdot C_{\ff} p^{1/2} n^{-2} \nonumber \\
&\geq - 8C_{\ff}p^{1/2} \left\{\theta_0 K_{\alpha_0,\nu}(0) + \tau_0 \right\}^{1/2} \tau_0^{-1} n^{-1},
\end{align}

for all sufficiently large $n$, where in the last inequality we used $\epsilon_2<1/2$.

We can also derive from (<ref>) and (<ref>) that on the event $\Bcal_0(\epsilon_1,\epsilon_2)\cap \Dcal_{n}$,

\begin{align}\label{eq:term4.lower1}
&~ \quad - \frac{1}{2} (\beta-\beta_0)^\T   F_n ^\T  \left[\left\{\theta K_{\alpha,\nu}(S_n) +\tau I_n\right\}^{-1}-\left\{\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right\}^{-1}\right] F_n (\beta-\beta_0) \nonumber \\
&\geq - \frac{1}{2} \left\|\left\{\theta K_{\alpha,\nu}(S_n) +\tau I_n\right\}^{-1}-\left\{\theta_0 K_{\alpha_0,\nu}(S_n) +\tau_0 I_n\right\}^{-1}\right\|_{\op} \cdot \| F_n (\beta-\beta_0)\|^2 \nonumber \\
&\geq - \frac{1}{2} \cdot \frac{2}{(1-\epsilon_2)\tau_0} \cdot \left(C_{\ff} p^{1/2} n^{-2}\right)^2 \nonumber \\
&\geq -\frac{2C_{\ff}^2 p}{\tau_0} n^{-4} ,
\end{align}

where in the last inequality we used $\epsilon_2<1/2$. The conclusion of Lemma <ref> follows from (<ref>) and (<ref>).

§.§ Other Technical Lemmas

The following results hold for any $n\times n$ symmetric matrices $A_1,A_2$.

(i) If $A_1\preceq A_2$, then for any $n\times n$ matrix $B$, $B^\T  A_1 B  \preceq B^\T  A_2 B$;

(ii) If $A_1\preceq A_2$, then for any positive semidefinite matrix $B$, $\tr(BA_1) \leq \tr(BA_2)$ and $\tr(A_1B)\leq \tr(A_2B)$.

(i) If $A_1\preceq A_2$, then $A_2-A_1$ is positive semidefinite, which means that $B^\T  (A_2-A_1) B $ is positive semidefinite. Therefore, $B^\T  A_1 B  \preceq B^\T  A_1 B  + B^\T  (A_2-A_1) B = B^\T  A_2 B$.

(ii) If $A_1\preceq A_2$, then $A_2-A_1$ is positive semidefinite, and

\begin{align*}
\tr(BA_2)-\tr(BA_1) &= \tr\left\{B(A_2-A_1)\right\} = \tr\left\{B^{1/2}(A_2-A_1)B^{1/2}\right\} \geq 0,\\
\tr(A_2B)-\tr(A_1B) &= \tr\left\{(A_2-A_1)B\right\} = \tr\left\{B^{1/2}(A_2-A_1)B^{1/2}\right\} \geq 0,
\end{align*}

where $B^{1/2}$ denotes the symmetric square root matrix of a positive semidefinite matrix $B$.

(Hanson-Wright Inequality, Proposition 1.1 in 19, Theorem 1.1 in 35) Let $Z_1,\ldots,Z_n$ be i.i.d. $\mathcal{N}(0,1)$ random variables and $Z=(Z_1,\ldots,Z_n)^\T $. Let $\Sigma$ be a positive definite matrix. Then for any $z>0$,

\begin{align*}
&\pr \left\{Z^\T  \Sigma Z \geq \tr(\Sigma) + 2\big\{\tr(\Sigma^2)z\big\}^{1/2} + 2\|\Sigma\|_{\op} z \right\} \leq \ee^{-z}.
\end{align*}

Furthermore, there exists a universal absolute constant $C_{\HW}>0$ that does not depend on $\Sigma$, such that for any $\epsilon>0$,

\begin{align*}
&\pr \left\{ \left|Z^\T  \Sigma Z - \tr(\Sigma)\right| \geq \epsilon \right\} \leq 2 \exp\left(-C_{\HW}\min\left\{\frac{\epsilon^2}{\|\Sigma\|_{F}^2}, \frac{\epsilon}{\|\Sigma\|_{\op}}\right\}\right).
\end{align*}

§ POSTERIOR INCONSISTENCY FOR $\ALPHA$ AND $\BETA$

Suppose that Assumptions <ref> and <ref> hold where the dimension of the domain $[0,1]^d$ is $d\in \{1,2,3\}$, and $\theta K_{\alpha,\nu}$ is the isotropic Matérn covariance function in Example <ref>. Suppose that we assign the uniform prior on the following set which covers the true parameter $(\theta_0,\alpha_0,\tau_0,\beta_0)$:

\begin{align*}
\Scal_{\Pi} & = \Big\{(\theta,\alpha,\tau,\beta): \theta \in [\theta_0/2,2\theta_0], \alpha \in [\alpha_0/2,2\alpha_0], \tau\in [\tau_0/2,2\tau_0],\\
& ~~~~ \beta_j\in [\beta_{0j}/2, 2\beta_{0j}],~ j=1,\ldots,p. \Big\}
\end{align*}

Then the posterior $\Pi(\cdot \mid Y_n, F_n )$ is inconsistent for $\alpha$ and $\beta$, i.e.,  there exist constants $\epsilon_0>0,\delta_0>0,\eta_0\in (0,1)$, such that for any $N\in \NN$, there exists $n>N$, such that

\begin{align*}
& \Pr\left(\Pi(|\alpha-\alpha_0|\geq \epsilon_0 |Y_n, F_n ) \geq \delta_0\right)\geq \eta_0, \\
\text{and } & \Pr\left(\Pi(\|\beta-\beta_0\|\geq \epsilon_0 |Y_n, F_n ) \geq \delta_0\right) \geq \eta_0.
\end{align*}

We prove the posterior inconsistency for $\alpha$ and $\beta$ by contradiction. Suppose that the conclusion does not hold, i.e.,  for any $\epsilon>0,\delta>0,\eta\in (0,1)$, there exists a $n_0(\epsilon,\delta,\eta)\in \NN$, such that for all $n>n_0(\epsilon,\delta,\eta)$, $\Pi(|\alpha-\alpha_0|\geq \epsilon |Y_n, F_n ) < \delta$ and $\Pi(\|\beta-\beta_0\|\geq \epsilon |Y_n, F_n ) < \delta$ with $\PP_{(\theta_0,\alpha_0,\tau_0,\beta_0)}$-probability at least $1-\eta$. For a sufficiently small $\epsilon>0$, the uniform prior density $\pi(\theta,\alpha,\tau,\beta)$ is a constant in an open neighborhood of $(\theta_0,\alpha_0,\tau_0,\beta_0)$, denoted by

$$\Ucal(\epsilon)= \left \{(\theta,\alpha,\tau,\beta): |\theta-\theta_0| < \epsilon, |\alpha-\alpha_0|<\epsilon, |\tau-\tau_0|<\epsilon,\|\beta-\beta_0\|_{\infty} <\epsilon/\sqrt{p} \right\},$$

where $\|\beta\|_{\infty}=\max_{1\leq j\leq p}|\beta_j|$. Under Assumption <ref>, the log-likelihood function $\Lcal_n(\theta,\alpha,\tau,\beta)$ is clearly finite for every $(\theta,\alpha,\tau,\beta)\in \Ucal(\epsilon)$. Therefore, the posterior distribution, provided exists, has a continuous density on $\Ucal(\epsilon)$. Now set $\delta=1/8$. Then with $\PP_{(\theta_0,\alpha_0,\tau_0,\beta_0)}$-probability at least $1-\eta$, the posterior median of $\alpha$, denoted by $\widetilde \alpha_n$, is in the interval $(\alpha_0-\epsilon,\alpha_0+\epsilon)$. It follows that for any $\epsilon>0,\eta>0$, for all $n>n_0(\epsilon,1/8,\eta)$, $|\widetilde\alpha_n-\alpha_0|<\epsilon$ with $\PP_{(\theta_0,\alpha_0,\tau_0,\beta_0)}$-probability at least $1-\eta$. Thus, $\widetilde\alpha_n$ is a consistent frequentist estimator for $\alpha_0$. Similarly, we can take the posterior median of each component $\beta_j$ ($j=1,\ldots,p$), denoted by $\widetilde \beta_{nj}$, and let $\widetilde\beta_n = (\widetilde \beta_{n1},\ldots,\widetilde \beta_{np})^\T $. Each $\widetilde \beta_{nj}$ is a consistent frequentist estimator for $\beta_{0j}$. As a result, $\widetilde\beta_n$ is a consistent frequentist estimator for $\beta_0$ given that $\big\|\widetilde\beta_n-\beta_0\big\|^2= \sum_{j=1}^p \big|\widetilde\beta_{nj}-\beta_{0j}\big|^2 < p\cdot (\epsilon^2/p) = \epsilon^2$. In this way, we have obtained two consistent frequentist estimators $\widetilde\alpha_n$ for $\alpha_0$ and $\widetilde\beta_n$ for $\beta_0$, both of which only depend on the data $(Y_n, F_n)$.