Wrapped Gaussian Process Functional Regression Model for Batch Data on Riemannian Manifolds

Given a $d$-dimensional smooth differentiable manifold $\mathcal{M}$ and a tangent space $T_{p}\mathcal{M}$, for $p\in\mathcal{M}$, a Riemannian metric $g_{p}:T_{p}\mathcal{M}\times T_{p}\mathcal{M}\rightarrow\mathbb{R}$ on $\mathcal{M}$ is a family of positive definite inner products which vary smoothly with $p\in\mathcal{M}$. Equipped with this Riemannian metric, we call the pair $(\mathcal{M},g)$ a Riemannian manifold.

The tangent bundle of $\mathcal{M}$ is defined as a set $\mathcal{TM}=\cup_{p\in\mathcal{M}}(p\times T_{p}\mathcal{M})$, where $p$ is a point in $\mathcal{M}$ and $v\in T_{p}\mathcal{M}$ is a tangent vector at $p$. If $\gamma$ is a smooth curve in $\mathcal{M}$, then its length is defined as the integral of $\|\partial\gamma/\partial t\|$ where the norm is computed via the Riemannian metric at $\gamma(t)$. For any pair $(p,v)\in\mathcal{TM}$ with $\|v\|$ sufficiently small, there is a unique curve $\gamma:[0,1]\rightarrow\mathcal{M}$ of minimal length between two points $p\in\mathcal{M}$ and $q\in\mathcal{M}$, with initial conditions $\gamma(0)=p$, $\gamma^{\prime}(0)=v$ and $\gamma(1)=q$. Such a curve is called a geodesic. The exponential map $\text{Exp}(p,v):\mathcal{M}\times T_{p}\mathcal{M}\rightarrow\mathcal{M}$ at $p\in\mathcal{M}$ is defined as $\text{Exp}(p,v)=\gamma(1)$ for sufficiently small $v\in T_{p}\mathcal{M}$. Consider the ball of largest radius around the origin in $T_{p}\mathcal{M}$ on which $\text{Exp}(p,v)$ is defined and let $V(p)\subset\mathcal{M}$ denote the image of $\text{Exp}(p,\cdot)$ on this ball. Then the exponential map has an inverse on $V(p)$, called the Riemannian logarithm map at $p$, $\text{Log}(p,\cdot):V(p)\rightarrow T_{p}\mathcal{M}$. For any point $q\in V(p)$, the geodesic distance of $p\in\mathcal{M}$ and $q\in\mathcal{M}$ is given by $d_{\mathcal{M}}(p,q)=\|\text{Log}(p,q)\|$.

Figure 2 shows the concepts above. For example, $p$ and $q$ are two points on a Riemannian manifold $\mathcal{M}$, the tangent space at $p$ is denoted as $T_{p}(\mathcal{M})$ and $\gamma(t)$ refers to a geodesic between $p$ and $q$. The tangent vector $v$ is from $p$ to $q$ at $T_{p}(\mathcal{M})$ and $p$ moves towarding $v$ with distance $\|v\|$ is $q$, which represents $\text{Log}(p,q)$ and $\text{Exp}(p,v)$ respectively.

The exponential map and its inverse are analogs of vector space operations in the following way. Suppose $p_{\mathcal{M}}$ and $q_{\mathcal{M}}$ refer to two different points on $\mathcal{M}$ linked by a geodesic $\gamma$, and suppose the tangent vector from the former to the latter is denoted by $pq_{\mathcal{M}}=(\frac{d\gamma}{dt})_{p}$. In addition, suppose $p_{E}$ and $q_{E}$ refer to two Euclidean data points and the vector from the former to the latter is given by $pq_{E}$. Table 1 shows the relationship of operations between Riemannian manifolds and Euclidean spaces.

If $(\mathcal{M},d_{\mathcal{M}})$ is a complete metric space that every Cauchy sequence of points in $\mathcal{M}$ has a limit that is also in $\mathcal{M}$, then the Hopf-Rinow theorem shows that every pair of points $p,q\in\mathcal{M}$ is joined by at least one geodesic, a condition known as geodesic completeness, and, equivalently, $\text{Exp}(p,v)$ is defined for all $(p,v)\in\mathcal{TM}$. However, given two points on $\mathcal{M}$, the geodesic between the points may not be unique, even if $\mathcal{M}$ is geodesically complete. The injectivity radius at $p\in\mathcal{M}$, denoted $\textrm{inj}_{p}$, is defined to be the radius of the largest ball at the origin of $T_{p}\mathcal{M}$ on which $\text{Exp}(p,\cdot)$ is a diffeomorphism. It follows that the geodesic ball at $p$ of radius $\textrm{inj}_{p}$ is contained within $V(p)$ under which the Wrapped Gaussian process regression can be derived to a analytical form. In the following, we assume the data are on a Riemannian manifold with infinite injectivity radius, such as spheres and Kendall’s shape spaces.

Some examples of Riemannian manifolds are shown in Appendix, including explicit formulae for the exponential and logarithm maps in certain cases.

As described in [34], a general Gaussian process regression model is defined as

where $\boldsymbol{x}\in\mathbb{R}^{Q}$, $y\in\mathbb{R}$, $\epsilon\sim\mathcal{N}(0,\sigma^{2})$. The regression function $f$ is assumed to follow a Gaussian process prior, that is

where $\mu(\cdot)=E[f(\cdot)]$ is the mean function of Gaussian process, and $k(\cdot,\cdot;\boldsymbol{\theta})$ is the covariance function with hyper-parameter $\boldsymbol{\theta}$.

Given training data $\mathcal{D}=\{\boldsymbol{x}_{i},y_{i},i=1,\ldots,n\}$, we have a multivariate Gaussian distribution that is

where $\boldsymbol{X}$ refers to $(\boldsymbol{x}_{1},...,\boldsymbol{x}_{n})$, $\boldsymbol{\mu}(\boldsymbol{X})=(\mu(\boldsymbol{x}_{1}),...,\mu(\boldsymbol{x_{n}}))$ refers to a mean vector, and $K(\boldsymbol{X},\boldsymbol{X})$ refers to a $n\times n$ covariance matrix in which the entry of the $(i,j)$-th element is $k(\boldsymbol{x}_{i},\boldsymbol{x}_{j};\boldsymbol{\theta})$.

To obtain the predictive distribution at a new input $\boldsymbol{x}^{*}$, it is convenient to take advantage of the properties of Gaussian distribution:

where $K(\boldsymbol{x}^{*},\boldsymbol{X})$ is a $1\times n$ covariance matrix, $K(\boldsymbol{x}^{*},\boldsymbol{x}^{*})$ is scalar and their entries are similar to $K(\boldsymbol{X},\boldsymbol{X})$. Fortunately, the predictive distribution has an analytical form which is

For more details of Gaussian process regression, see [44].

The kernel function in the covariance function can be chosen from a parametric family, such as the radial basis function (RBF) kernel:

with the hyper-parameter $\boldsymbol{\theta}=\{\lambda_{1},\lambda_{2}\}$ chosen by maximizing the marginal likelihood. Specifically, the hyper-parameters can be estimated by conjugate gradient descent algorithm [32] or Markov Chain Monte Carlo [38].

In this subsection, we first review the definitions of the wrapped Gaussian distribution and wrapped Gaussian process on a Riemannian manifold $\mathcal{M}$, and then review wrapped Gaussian process regression (WGPR) which is introduced by [27].

Suppose $\mathcal{M}$ denotes a $d$-dimensional Riemannian manifold and $X$ is a random variable on $\mathcal{M}$. For $\mu\in\mathcal{M}$ and a symmetric positive definite matrix $K\in\mathbb{R}^{d\times d}$, we can define a wrapped Gaussian distribution as follows

where $\mu$ is called basepoint and $\boldsymbol{v}$ is a vector belongs to the tangent space at $\mu$ which follows a multi-variate Gaussian distribution with zero mean and covariance matrix $K$ (see Section 3.1 in [27]). The wrapped Gaussian distribution is formally denoted as

[27] describe a manifold-valued Gaussian process constructed by the wrapped Gaussian distribution, and the following material is based on their work. A collection $(f(X_{1}),...,f(X_{n}))$ of random points on a manifold $\mathcal{M}$ index over a set $\Omega$ is a wrapped Gaussian process, if every subcollection follows a jointly wrapped Gaussian distribution on $\mathcal{M}$. It is denoted as

where $\mu(\cdot)=E[f(\cdot)]$ is the mean function on a Riemannian manifold, and $k(\cdot,\cdot;\boldsymbol{\theta})$ is a covariance function with hyper-parameter $\boldsymbol{\theta}$ on a tangent space at $\boldsymbol{\mu}(\cdot)$.

Consider the wrapped Gaussian process regression model which relates the manifold-valued response $p\in\mathcal{M}$ to the vector-valued predictor $\boldsymbol{x}\in\mathbb{R}^{Q}$ though a link function $f:\mathbb{R}^{Q}\rightarrow\mathcal{M}$ which we assume to be a wrapped Gaussian process. Analogously to Gaussian process regression, the authors also derive inferences for wrapped Gaussian process regression in Bayesian framework. Specifically, given the training data $\mathcal{D}_{\mathcal{M}}=\{(\boldsymbol{x}_{i},p_{i})|\boldsymbol{x}_{i}\in\mathbb{R}^{Q},p_{i}\in\mathcal{M},i=1,\cdots,n\}$, the prior for $f(\cdot)$ is assumed to be a wrapped Gaussian process, which is $f(\cdot)\sim\mathcal{GP}_{\mathcal{M}}(\mu(\cdot),K(\cdot,\cdot;\boldsymbol{\theta}))$. Then the joint distribution between the training outputs $\boldsymbol{p}=(p_{1},\cdots,p_{n})$ and the test outputs ${p}^{*}$ at $\boldsymbol{x}^{*}$ is given by

where $\mu^{*}=\mu(\boldsymbol{x}^{*})$, $\boldsymbol{\mu}=(\mu(\boldsymbol{x}_{1}),...,\mu(\boldsymbol{x_{n}}))$, $K=K(\boldsymbol{x},\boldsymbol{x})$, $K^{*}=K(\boldsymbol{x},\boldsymbol{x}^{*})$, $K^{**}=K\left(\boldsymbol{x}^{*},\boldsymbol{x}^{*}\right)$ with $\boldsymbol{x}=(\boldsymbol{x}_{1},\cdots,\boldsymbol{x}_{n})$.

Therefore, by Theorem 1 in [27], the predictive distribution for new data given the training data can be obtained by

where

Similar to the comparison of operators between Riemannian manifolds and Euclidean spaces, Table 2 shows the analogs of Gaussian processes between these two spaces.

It is difficult to infer non-linear relationships non-parametrically between a functional response variable and multi-dimensional functional covariates even when both of them are observed in a space for real-valued functions [37]. Therefore, we assume that the mean structure $\mu_{m}(\cdot)$ depends on the scalar covariates $\boldsymbol{u}_{m}$ only. Specifically we assume

where $\mu_{*}(t)\in\mathcal{M}$ and $\boldsymbol{\beta}(t)\in T_{\mu_{*}(t)}\mathcal{M}$ play a similar role to the intercept term and slope term of the function-on-scalar linear regression model in Euclidean space. These are model parameters which we estimate. Furthermore, equation (19) can be seen as a functional version of geodesic regression proposed in [17], in which a geodesic regression model for manifold-valued response and scalar-valued predictor is proposed in the form $y=\text{Exp}(p,xv)$ with $(p,v)\in\mathcal{TM}$ where $x\in\mathbb{R}$ is the predictor variable.

Conditional on $\mu_{m}$, the $m$-th curve $y_{m}$ is modelled as

where $\tau_{m}(t)\in T_{\mu_{m}(t)}\mathcal{M}$ is a vector field along $\mu_{m}$. When the data are observed with random noise, [17] defines the geodesic model as $y=\text{Exp}(\text{Exp}(p,xv),\epsilon)$. Similarly, the proposed WGPFR model with error term can be defined as

where $\epsilon_{m}(t)$ are independent measurement errors following a multivariate normal distribution with mean zero and covariance matrix $\sigma^{2}\mathbf{I}$. For functional data, in addition to the measurement error term, the most important part we are interested is the underlying covariance structure contained in term $\tau_{m}(t)$, which is Equation (20). In what follows, we explain how $\tau_{m}(t)$ is modelled as depending on the functional covariate $\mbox{{\boldmath${x}$}${}_{m}$}(t)$.

We assume that the vector field $\tau_{m}(t)\in T_{\mu_{m}(t)}\mathcal{M}$ can be represented in local coordinates as a $d$-dimensional vector function $\tau_{m}(t)=(\tau_{m,1}(t),...,\tau_{m,d}(t))\in\mathbb{R}^{d}$. We further assume this vector-valued function follows a Gaussian process with mean zero and an underlying covariance structure $k_{m}(\cdot,\cdot;\boldsymbol{\theta})$, which is

Specifically, given two observation times $t_{i}$ and $t_{mj}$ from the $m$-th batch of data, we assume that the covariance of the two random vectors $\tau_{m}(t_{i})$ and $\tau_{m}(t_{mj})$ a function of the observed covariate values $\mbox{{\boldmath${x}$}${}_{m}$}(t_{i})$ and $\mbox{{\boldmath${x}$}${}_{m}$}(t_{mj})$.

Under the framework of the wrapped Gaussian process in [27], we can represent the above Gaussian process as the wrapped Gaussian process for manifold functional data. That is,

where $\mu_{m}(t)\in\mathcal{M}$ and $k_{m}(\cdot,\cdot;\boldsymbol{\theta})\in T_{\mu_{m}(t)}\mathcal{M}$ are defined as above.

To summarize, the WGPFR model studies the non-linear relationship between the manifold-valued functional response and mixed (functional and vector) Euclidean-valued predictors in two parts, as shown in Figure 3. The first part is the mean structure which is defined through a functional regression to model the manifold-valued functional mean with the vector-valued predictor. The second part is called the covariance structure, and for which the inference is based on a Bayesian method with the prior assumed to be Gaussian process.

Compared to the geodesic regression in [17], the WGPFR additionally incorporates the information of dependence between data modelled by covariance structure which is the main focus of interest for functional data analysis. Compared to the WGPR model in [27], this paper further models the mean part which relates the functional mean on the manifold to the scalar predictor on Euclidean space, while the mean part in [27] is not related to scalar predictors. Note that the scalar predictor plays an important role in the prediction of the model, especially for the batch data we consider in this paper. The proposed WGPFR model is more complicated but much flexible than existing models and induces a more complex inference procedure which will be illustrated in detail in the next subsection.

Since the mean structure model (19) can be seen as a functional version of the geodesic regression model [17], an intuitive inference method consists of optimizing $\mu_{*}(t)$ and $\boldsymbol{\beta}(t)$ simultaneously for each fixed $t$. However, this would ignore information about the correlation of the function $\mu_{*}$ and $\boldsymbol{\beta}(t)$ at different times. On the other hand, if the number of observed times $t_{i}$ is large, and we account fully for correlations between all different times, then the computational cost will be very high, unless we use some efficient methods, such as [12]. Due to the dependent structure, we cannot conduct optimization at each fixed $t$ separately, otherwise the correlation information across different t will be lost. It can also be seen that inference of the functional mean structure model is indeed challenging.

In this paper we borrow an idea of [11] to estimate $\mu_{*}(t)$ and $\boldsymbol{\beta}(t)$. Note that $\mu_{*}(t)$ plays the role of an intercept in the linear regression in Euclidean space, which is usually defined as the mean of the response. Similarly, for functional data $y_{m}(t),m=1,...,M$ on a Riemannian manifold, we can define the intrinsic population Fréchet mean function under the assumption of existence and uniqueness:

where $d_{\mathcal{M}}(\cdot,\cdot)$ denotes geodesic distance. When $t$ is fixed, $\mu_{0}(t)$ refers to a point on a Riemannian manifold; when $t$ is variable, $\mu_{0}(t)$ refers to a curve on a Riemannian manifold. Additionally, since each $y_{m}(t)\in\mathcal{M}$ is continuous, $\mu_{0}(t)$ is also continuous (see equation (1) in [11]).

The intrinsic Fréchet mean function $\mu_{0}(t)$ can be a good choice for $\mu_{*}(t)$ since it is a natural mean of manifold-valued functional data and it is reasonable to assume $y_{m}(t),m=1,...,M$ are confined to stay within the radius of injectivity at $\mu_{0}(t)$ for all $t\in\mathcal{T}$. Then the inverse exponential map can be well-defined. We estimate $\mu_{*}(t)$ at each time $t$ via the sample Fréchet mean of the points $y_{m}(t)$. This requires that the times $t_{i}$ are the same across all curves $m=1,\ldots,M$. Standard gradient descent methods [29] can be used to obtain the Fréchet sample mean at each time $t=t_{i}$ to obtain estimates $\hat{\mu}_{0}(t_{i})$ of $\mu_{*}(t_{i})$, $i=1,\ldots,N$. Denote the logarithm process data as

where $V_{m}(t)$ is a $d$-dimensional tangent vector of $y_{m}(t)$ at $\mu_{0}(t)$ for a fixed $t$. Since we assume $\mathcal{M}$ is a $d$ dimensional closed Riemannian submanifold of a Euclidean space, it is reasonable to identify the tangent space $T_{\mu_{0}(t)}\mathcal{M}$ as a subspace of the Euclidean space $\mathbb{R}^{d_{0}}$ with $d_{0}\geq d$. Then $V_{m}(t)$ can be seen as $\mathbb{R}^{d_{0}}$ valued square integrable functions, which can be represented by a set of basis functions. Similarly, the slope function $\boldsymbol{\beta}(t)$ can also be represented by this set of basis functions. In this way, we can not only easily estimate the intercept term $\mu_{0}(t)$ but also the slope term $\boldsymbol{\beta}(t)$ can be obtained by least square methods on $\mathbb{R}^{d_{0}}$. This reduces the computation complexity compared with the optimization method in [17]. Next we illustrate the estimation procedure in detail.

For the estimation of the intrinsic Fréchet mean $\mu_{0}(t)$, the sample Fréchet mean is calculated by minimizing the following loss function for each $t\in\mathcal{T}$:

It is natural to use the gradient descent algorithm introduced by [29] to estimate the sample Fréchet mean for each fixed $t$. Using [14], the gradient is given by

An implementation of the gradient descent algorithm to solve equation (25) is given in Algorithm 1. In practice, we choose the step size $l=\frac{1}{N}$ where $N$ is the number of samples and $\epsilon=10^{-8}$ which is relatively small.

For the estimation of $\boldsymbol{\beta}(t)=\big{(}\boldsymbol{\beta}_{1}(t),\cdots,\boldsymbol{\beta}_{p}(t)\big{)}$, the first step is to find a reasonable basis to represent $\boldsymbol{\beta}(t)$ in the tangent space. Under the framework of [11], suppose we have an arbitrary system of $K$ orthonormal basis functions

where $\mathbb{H}=\{\nu:\mathcal{T}\rightarrow\mathbb{R}^{d},\int_{\mathcal{T}}\nu(t)^{T}\nu(t)dt<\infty\}$ is a $L^{2}$ Hilbert space of $\mathbb{R}^{d}$ with inner product $\langle\nu,u\rangle=\int_{\mathcal{T}}\nu(t)^{T}u(t)dt\}$ and norm $\|v\|=\langle v,v\rangle^{1/2}$, $\delta_{kl}=1$ if $k=l$ and 0 otherwise. Note that the values of each $\phi_{k}(t)$ at each time $t\in\mathcal{T}$ is restricted to the tangent space $T_{\mu_{0}(t)}\mathcal{M}$, which are identified with $\mathbb{R}^{d_{0}}$. Define the $K$-dimensional linear subspace of $\mathbb{H}$:

Then the slope functions $\boldsymbol{\beta}_{j}(t)\in T_{\mu_{0}(t)}\mathcal{M}$, $j=1,\cdots,p$ can be approximated by a linear span on $\mathcal{M}_{K}\left(\Phi_{K}\right)$ with expansion coefficients $b_{jk}$:

Then the estimation of $\boldsymbol{\beta}_{j}(t)$ can be transformed to the estimation of $b_{jk}$, $k=1,2,\cdots,K$, $j=1,2,\cdots,p$.

Defining $V_{m}(t)=\text{Log}(\mu_{0}(t),y_{m}(t))$, $\hat{V}_{m}(t)=\text{Log}(\hat{\mu}_{0}(t),y_{m}(t))$ and $W_{m}(t)=\text{Log}(\mu_{0}(t),\mu_{m}(t))$, then the multiple linear function-on-scalar regression model becomes

where $e_{m}(t)=V_{m}(t)-W_{m}(t)$ which are assumed to be independent of the covariate $\boldsymbol{u}_{m}$ for $m=1,\cdots,M$.

Since $V_{m}(t)$ can not be obtained directly, we use $\hat{V}_{m}(t)$ as the response of the multiple linear regression model, so that the parameters of the regression model are estimated by minimizing the following loss function with observed data:

The $b_{jk}$ coefficients are calculated by standard least squares methods.

To be more specific, let $\hat{\textbf{V}}=(\hat{V}_{im})_{i=1,\cdots,nd_{0};m=1,\cdots,M}$ denote a $nd_{0}\times M$ matrix such that the $m$-th column denotes the $m$-th batch, and the $i$-th row denotes the vectored element $\bigg{(}\hat{V}_{m1}(t_{1}),\cdots,\hat{V}_{m1}(t_{n}),\cdots,\hat{V}_{md_{0}}(t_{1}),\cdots,\hat{V}_{md_{0}}(t_{n})\bigg{)}$. Similarly, let $\boldsymbol{\Phi}=(\phi_{ik})_{i=1,\cdots,nd_{0};k=1,\cdots,K}$ denote a $nd_{0}\times K$ matrix; let $U$ denote a $M\times p$ matrix for which each row $\textbf{u}_{m}$ indicates a data point; and let $B$ denote the $p\times K$ coefficient matrix with element $b_{jk}$ for $j=1,\cdots,p$, $k=1,\cdots,K$. Then the above loss function can be rewritten as

The least squares estimate of $B$ is obtained by

Consistency of the estimator of the sample Fréchet mean was proved in [11]. The following theorem establishes consistency of the estimator of the regression coefficients; the proof is left to the appendix.

It follows that the functional slope coefficients can be estimated by $\hat{\boldsymbol{\beta}}(t)=\big{(}\hat{\boldsymbol{\beta}}_{1}(t),\cdots,\hat{\boldsymbol{\beta}}_{p}(t)\big{)}$ where

for $j=1,\cdots,p$. The estimated mean structure for the $m$-th curve is given by

The function $\hat{\mu}_{m}(t)$ is expected to approximate $\mu_{m}(t)$. However we omit the rigorous proof of the consistency of $\hat{\mu}_{m}(t)$ on the manifold due to some technique difficulties.

From the section above, we obtain an estimate of the mean structure $\mu_{m}(t)$ which is a continuous function on $\mathcal{M}$. In this section, we will focus on inference of the covariance structure $\tau_{m}(t)\in T_{\mu_{m}(t)}\mathcal{M}$. We assume that the tangent spaces $T_{\mu_{m}(t)}\mathcal{M}$ for $t\in\mathcal{T}$ can be identified with $\mathbb{R}^{d_{0}}$ via some smooth local basis of orthonormal vector fields along $\mu_{m}(t)$. We have mentioned above that the covariance structure can be related to another functional covariate $x_{m}(t)$, which is a $\mathbb{R}^{q}$-valued function. If we denote $\tau_{m}(t)=\tau_{m}(x_{m}(t))=(\tau_{m1}(t),\cdots,\tau_{md_{0}}(t))$, then the correlation of different components in $\tau_{m}(t)$ could be estimated via a cross-covariance function model, such as the convolved Gaussian process [6]. However, if there are $n$ observations $\tau_{md}(t_{i})$, $i=1,\cdots,n$, then the size of a covariance matrix in a Gaussian process of $\tau_{md}(t)$, $d=1,\cdots,d_{0}$ is $n\times n$ while the size of a cross-covariance matrix in a convolved Gaussian process of $\tau_{m}(t)$ is $n^{d_{0}}\times n^{d_{0}}$, which is computationally expensive. As a result, we will consider different dimensions independently. Specifically, we assume

where $K_{md}(\cdot,\cdot;\boldsymbol{\theta}_{md})$ denotes a covariance kernel depending on hyper-parameter $\boldsymbol{\theta}_{md}$.

Given the value of $\tau_{md}(t_{i})$, $i=1,\cdots,n$, denote $\boldsymbol{\tau}_{md}=(\tau_{md}(t_{m1}),...,\tau_{md}(t_{mn}))$, then for any new input $t_{m}^{*}$, the conditional distribution of $\tau_{md}(t_{m}^{*})$ is

where $\boldsymbol{k}_{md}^{*}=(K_{md}(t_{m1},t_{m}^{*}),...,K_{md}(t_{mn},t_{m}^{*}))$ and $k_{md}^{**}=K_{md}(t_{m}^{*},t_{m}^{*})$.

In Gaussian process regression, the hyper-parameters for the $d$-th dimension, $\boldsymbol{\theta}_{md}$, can be estimated by maximizing the sum of marginal likelihood functions for each batch, i.e. maximizing

Under some regularity conditions, the estimator $\hat{\theta}_{md}$ converges to $\theta_{md}$ [9]. With the estimated hyper-parameters, the estimated covariance structure of the $d$-th coordinate for the $m$-th batch is given by

where $\boldsymbol{k}_{md}(t)=(K_{md}(t,\boldsymbol{x}_{m}(t_{m1})),\cdots,K_{md}(t,\boldsymbol{x}_{m}(t_{mn})))$. Under some mild conditions, the above estimator of the covariance structure is information consistent, as the following theorem states.