Kernel Packet: An Exact and Scalable Algorithm for Gaussian Process Regression with Matérn Correlations

Let $Y(\mathbf{x})$ denote a stationary GP prior on $\mathbb{R}^{d}$ with mean function $\mu(\mathbf{x})$, variance $\sigma^{2}$, and correlation function $K(\mathbf{x},\mathbf{x}^{\prime})$. The correlation function is also referred to as a “kernel function” in the language of applied math or machine learning (Rasmussen, 2006). When $d=1$, there are two types of popular correlation functions. The first type is the Matérn family (Stein, 1999):

$$K(x,x^{\prime})=\frac{2^{1-\nu}}{\Gamma(\nu)}\left(\sqrt{2\nu}\frac{\lvert x-x^{\prime}\rvert}{\omega}\right)^{\nu}K_{\nu}\left(\sqrt{2\nu}\frac{\lvert x-x^{\prime}\rvert}{\omega}\right),$$

(1)

for any $x,x^{\prime}\in\mathbb{R}$, where $\nu>0$ is the smoothness parameter, $\omega>0$ is the scale and $K_{\nu}$ is the modified Bessel function of the second kind. The smoothness parameter $\nu$ governs the smoothness of the GP $Y$ (Santner et al., 2003; Stein, 1999); the scale parameter $\omega$ determines the spread of the correlation (Rasmussen, 2006). Matérn correlation functions are widely used because of its great flexibility. The second type is the Gaussian family:

$$K(x,x^{\prime})=\exp\left(-\frac{\lvert x-x^{\prime}\rvert^{2}}{\omega}\right),$$

(2)

for any $x,x^{\prime}\in\mathbb{R}^{d}$. A Gaussian kernel function is the limit of a sequence of Matérn kernels with the smoothness parameter tending to infinity. The sample paths generated by GP with Gaussian correlation function are infinitely differentiable.

For multi-dimensional problems, a typical choice of the correlation structure is the separable or product correlation:

$$K(\mathbf{x},\mathbf{x}^{\prime})=\prod_{j=1}^{d}K_{j}(x_{j},x^{\prime}_{j}),$$

(3)

for any $\mathbf{x}=(x_{1},\ldots,x_{d})^{T},\mathbf{x}^{\prime}=(x^{\prime}_{1},\ldots,x^{\prime}_{d})^{T}$, where $K_{j}$ is a one-dimensional Matérn or Gaussian correlation function for each $j$. This assumption ensures that the GP lives in a tensor space, and is key to the “tensor-space” techniques, which reduces the multi-dimensional problems to one-dimensional ones.

Suppose that we have observed $\mathbf{Y}=\big{(}Y(\mathbf{x}_{1}),\cdots,Y(\mathbf{x}_{n})\big{)}^{T}$ on $n$ distinct points $\mathbf{X}=\{\mathbf{x}_{i}\}_{i=1}^{n}$. The aim of GP regression is to predict the output at an untried input $\mathbf{x}^{*}$ by computing the distribution of $Y(\mathbf{x}^{*})$ conditional on $\mathbf{Y}$, which is a normal distribution with the following conditional mean and variance (Santner et al., 2003; Banerjee et al., 2014):

	$\displaystyle\operatorname{\mathbb{E}}\left[Y(\mathbf{x}^{*})\big{|}\mathbf{Y}\right]=\mu(\mathbf{x}^{*})+K(\mathbf{x}^{*},\mathbf{X})\mathbf{K}^{-1}\left(\mathbf{Y}-\boldsymbol{\mu}\right),$		(4)
	$\displaystyle\operatorname{\mathrm{Var}}\left[Y(\mathbf{x}^{*})\big{|}\mathbf{Y}\right]=\sigma^{2}\bigg{(}K(\mathbf{x}^{*},\mathbf{x}^{*})-K(\mathbf{x}^{*},\mathbf{X})\mathbf{K}^{-1}K(\mathbf{X},\mathbf{x}^{*})\bigg{)},$		(5)

where $\sigma^{2}>0$ is the variance, $K(\mathbf{x}^{*},\mathbf{X})=\left(K(\mathbf{X},\mathbf{x}^{*})\right)^{T}=\left(K(\mathbf{x}^{*},\mathbf{x}_{1}),\cdots,K(\mathbf{x}^{*},\mathbf{x}_{n})\right)$, $\mathbf{K}=\left[K(\mathbf{x}_{i},\mathbf{x}_{s})\right]_{i,s=1}^{n}$ and $\boldsymbol{\mu}=\left(\mu(\mathbf{x}_{1}),\cdots,\mu(\mathbf{x}_{n})\right)^{T}$.

In GP regression, the mean function $\mu$ is usually parametrized as a linear form $\mu=\sum_{i=1}^{p}\beta_{i}f_{i}$ for some unknown coefficient vector $\boldsymbol{\beta}=\big{(}\beta_{1},\cdots,\beta_{p}\big{)}^{T}$ and known regression functions $f_{1},\cdots,f_{p}$. To improve the predictive performance of GP regression, the coefficient vector $\boldsymbol{\beta}$, variance $\sigma^{2}$ and scales $\boldsymbol{\omega}=\big{(}\omega_{1},\cdots,\omega_{d}\big{)}^{T}$ associated to each one-dimensional correlation function $k_{j}$ are usually estimated via maximum likelihood (Jones et al., 1998). The log-likelihood function given the data is:

$$L(\boldsymbol{\beta},\sigma^{2},\boldsymbol{\omega})=-\frac{1}{2}\left[n\log\sigma^{2}+\log\det(\mathbf{K})+\frac{1}{\sigma^{2}}\big{(}\mathbf{Y}-\mathbf{F}\boldsymbol{\beta}\big{)}^{T}\mathbf{K}^{-1}\big{(}\mathbf{Y}-\mathbf{F}\boldsymbol{\beta}\big{)}\right],$$

(6)

where $\det(\mathbf{K})$ denotes the determinant of the correlation matrix $\mathbf{K}$ and $\mathbf{F}$ is the $n\times p$ matrix whose $(i,s)^{\rm th}$ entry is $f_{s}(\mathbf{x}_{i})$. The Maximum Likelihood Estimator (MLE) is then defined as the maximimizer of the log-likelihood function: $(\hat{\boldsymbol{\beta}},\hat{\sigma}^{2},\hat{\boldsymbol{\omega}})=\operatorname*{argmax}_{\boldsymbol{\beta},\sigma^{2},\boldsymbol{\omega}}L(\boldsymbol{\beta},\sigma^{2},\boldsymbol{\omega})$.

In both GP regression and parameter estimation, the computation can become unstable or even intractable because it involves the pursuit of the inversion and the determinant of the correlation matrix $\mathbf{K}$. Each task takes $\mathcal{O}(n^{3})$ time if a direct method, such as the Cholesky decomposition, is applied, which is a fundamental computational challenge for GP regression.

When applicable, as a specialized algorithm, the proposed method is significantly superior to the existing alternatives. In this section, we compare the proposed method with a few popular existing approaches for large-scale GP regression. It is worth noting that, the fundamental mathematical theory for the proposed method differs from that of any of the existing methods. A summary of the comparisons is presented in Table 1.

In this section we introduce the theory that explains how we can find a compactly supported function in $\mathcal{K}$. Clearly, such a function must admit the representation $\phi(x)=\sum_{j=1}^{n}A_{j}K(x,x_{j}).$ Recall the requirement that the linear transform is sparse, which means that most of the coefficients $A_{j}$’s must be zero. This inspires the following definition.

At first sight, it seems to be too optimistic to expect the existence of KPs. But, surprisingly, these functions do exist for one-dimensional Matérn correlation functions with half-integer smoothness. We will show that if the smoothness parameter $\nu$ is a half integer, i.e., $\nu-1/2\in\mathbb{N}$, there is a KP of degree $k:=2\nu+2$ given any $k$ distinct input points.

For simplicity, we will use $k$ to parametrize the Matérn correlation, in other words, $\nu=(k-2)/2$ for $k=3,5,7,\ldots$. Let $\mathbf{a}=(a_{1},...,a_{k})^{T}$ be a vector with $a_{1}<\cdots<a_{k}$. The goal is to find coefficients $A_{j}$’s such that

$$\phi_{\mathbf{a}}(x):=\sum_{j=1}^{k}A_{j}K(x,a_{j})$$

(7)

is a KP. We will first find a necessary condition for $A_{j}$’s, and next we will prove that such a condition is also sufficient. We apply the Paley-Wiener theorem (see Lemma 15 in Appendix A and Stein and Shakarchi (2003)), which states that $\phi_{\mathbf{a}}(x)$ has a compact support only if the inverse Fourier transform of $\phi_{\mathbf{a}}$, denoted as $\tilde{\phi}_{\mathbf{a}}(x)$, can be extended to an entire function, i.e., a complex-valued function that is holomorphic on the whole complex plane. Let $i=\sqrt{-1}$. Our convention of inverse Fourier transform is $\tilde{f}(\xi)=(2\pi)^{-1/2}\int_{-\infty}^{\infty}f(x)e^{i\xi x}dx$. Direct calculations show

$$\tilde{\phi}_{\mathbf{a}}(x)\propto\left[\sum_{j=1}^{k}A_{j}\exp\{ia_{j}x\}\right](c^{2}+x^{2})^{-(k-1)/2},x\in\mathbb{R}.$$

Clearly, the analytic continuation of this function (up to a constant) is

$$\tilde{\phi}_{\mathbf{a}}(z)\propto\left[\sum_{j=1}^{k}A_{j}\exp\{ia_{j}z\}\right](c^{2}+z^{2})^{-(k-1)/2}\eqqcolon\gamma(z)(c^{2}+z^{2})^{-(k-1)/2},$$

and this function can be defined at any $z\in\mathbb{C}\setminus\{\pm ci\}$. Note that the function $(c^{2}+z^{2})^{-(k-1)/2}$ has poles at $z=\pm ci$, each with multiplicity $(k-1)/2$. According to Paley-Wiener theorem, we have to make $\tilde{\phi}_{\mathbf{a}}(z)$ an entire function, which implies that $\gamma(\pm ci)=0$, each with multiplicity $(k-1)/2$ as well. This condition leads to a set of equations¹¹1This statement is formalized as Lemma 20 in Appendix B.:

$\displaystyle\gamma^{(j)}(ci)=0,$

$\displaystyle\gamma^{(j)}(-ci)=0,$

for $j=0,\ldots,(k-3)/2$, where $\gamma^{(j)}$ denotes the $j$th derivative of $\gamma$. Clearly, there are $k-1$ equations, which can be rewritten as

$\displaystyle\sum_{j=1}^{k}A_{j}a_{j}^{l}\exp\{\delta ca_{j}\}=0,$

(8)

with $l=0,\ldots,(k-3)/2$ and $\delta=\pm 1$, which is a $(k-1)\times k$ linear system. All solutions to this system are real-valued vectors because all coefficients are real.

Next we study the property of the linear system (8) and the corresponding $\phi_{\mathbf{a}}$. Theorem 2 states that $\phi_{\mathbf{a}}$ can be uniquely determined by (8) up to a multiplicative constant.

Another important property of (8) is that its solution is not affected by a shift of $\mathbf{a}$. Define $\mathbf{a}+t=(a_{1}+t,\ldots,a_{n}+t)^{T}$.

Theorem 5 confirms that any non-zero $\phi_{\mathbf{a}}$ is indeed a KP.

In other words, we have the following Corollary 6.

Figure 1 illustrates that the linear combination of $5$ components $\{K(\cdot,a_{j})\}_{j=1}^{5}$ provides a compactly supported KP corresponding to Matérn-$3/2$ correlation function.

It is evident that KPs are highly non-trivial and precious. Their existence relies on the correlation function. Theorem 7 shows that many other correlation function do not admit any KP, and consequently, the proposed algorithm is not applicable to these correlations.

Theorem 8 shows that the KP constructed by (8) has the lowest degree.

Besides KPs, we need to introduce a set of functions to capture the “boundary effects” of Gaussian process regression. As before, let $\mathbf{a}=(a_{1},...,a_{s})^{T}$ be a vector with $a_{1}<\cdots<a_{s}$. We consider the functions

$$\phi_{\mathbf{a}}(x):=\sum_{j=1}^{s}A_{j}K(x,a_{j}),$$

(9)

with $(k+1)/2\leq s\leq k-1$ and a non-zero real vector $(A_{1},\ldots,A_{s})^{T}$. Then Theorem 8 suggests that $\phi_{\mathbf{a}}$ in (9) cannot have a compact support. Nevertheless, it is possible that the support of $\phi_{\mathbf{a}}$ is a half real line. In this case, we cal $\phi_{\mathbf{a}}$ a one-sided KP. Specifically, we call $\phi_{\mathbf{a}}$ a right-sided KP if $\operatorname{\mathrm{s}upp}\phi_{\mathbf{a}}=[a_{1},+\infty)$, and we call $\phi_{\mathbf{a}}$ a left-sided KP if $\operatorname{\mathrm{s}upp}\phi_{\mathbf{a}}=(-\infty,a_{s}]$.

First we consider right-sided KPs. We propose to identify $A_{j}$’s by solving

$$\sum_{j=1}^{s}A_{j}a_{j}^{l}\exp\{-ca_{j}\}=0,\quad\sum_{j=1}^{s}A_{j}a_{j}^{r}\exp\{ca_{j}\}=0,$$

(10)

where $l=0,\ldots,(k-3)/2$ and the second term of (10) comprises auxiliary equations for the case $s\geq(k+3)/2$ with $r=0,\ldots,s-(k+3)/2$. Similar to (8), (10) is an $(s-1)\times s$ linear system.

The following theorems describes the properties of the linear system (10) and the corresponding $\phi_{\mathbf{a}}$. Specifically, Theorem 11 confirms that $\phi_{\mathbf{a}}$ is indeed a right-sided KP.

Left-sided KPs are constructed similarly by solving the following equations:

$$\sum_{j=1}^{s}A_{j}a_{j}^{l}\exp\{ca_{j}\}=0,\quad\sum_{j=1}^{s}A_{j}a_{j}^{r}\exp\{-ca_{j}\}=0,$$

(11)

where $l=0,\ldots,(k-3)/2$ and the second term comprises auxiliary equations for the case $s\geq(k+3)/2$ with $r=0,\ldots,s-(k+3)/2$. The properties of left-sided KPs are analogous to these stated in Theorems 9-11, for which we omit the statements.

Let $x_{1}<\cdots<x_{n}$ be the input data, and $K$ a Matérn correlation function with a half-integer smoothness. Suppose $n\geq k$. We can construct the following $n$ functions, as a subset of $\mathcal{K}$:

1. 1.

   $\phi_{1},\phi_{2},\ldots,\phi_{(k-1)/2}$, defined as left-sided KPs $\phi_{(x_{1},\ldots,x_{(k+1)/2})},\phi_{(x_{1},\ldots,x_{(k+1)/2+1})},\ldots,\phi_{(x_{1},\ldots,x_{k-1})}$,

2. 2.

   $\phi_{(k+1)/2},\phi_{(k+1)/2+1},\ldots,\phi_{n-(k-1)/2}$, defined as KPs $\phi_{(x_{1},\ldots,x_{k})},\phi_{(x_{2},\ldots,x_{k+1})},\ldots,\phi_{(x_{n-k+1},\ldots,x_{n})}$,

3. 3.

   $\phi_{n-(k-3)/2},\ldots,\phi_{n-1},\phi_{n}$, defined as right-sided KPs $\phi_{(x_{n-k+2},\ldots,x_{n})},\ldots,\phi_{(x_{n-(k-1)/2-1},\ldots,x_{n})},\phi_{(x_{n-(k-1)/2},\ldots,x_{n})}$.

Note that KPs and one-sided KPs given the input points cannot be uniquely defined. They are unique only up to a non-zero multiplicative factor. Here the choice of these factors are nonessential. The general theory and algorithms in this article will be valid for each specific choice. Now we present Theorem 13, which, together with the fact that the dimension of $\mathcal{K}$ is $n$, implies that $\{\phi_{j}\}_{j=1}^{n}$ forms a basis for $\mathcal{K}$, referred to as the KP basis.

Further, it is straightforward to check via Theorems 5 and 11 that, given any $x\in\mathbb{R}$, the vector $\boldsymbol{\phi}(x)=\left(\phi_{1}(x),\ldots,\phi_{n}(x)\right)^{T}$ has at most $k-1$ non-zero entries. As a result, we have constructed a basis for $\mathcal{K}$ satisfying the two sparse properties mentioned at the beginning of Section 2. Figure 2 illustrates a KP basis corresponding to Matérn-$3/2$ and Matérn-$5/2$ correlation function with input points $\mathbf{X}=\{0.1,0.2,\ldots,1\}$.

The theory in Section 2 shows that for one-dimensional problems, given ordered and distinct inputs $x_{1}<\cdots<x_{n}$, the correlation matrix $\mathbf{K}$ admits a sparse representation as

$$\mathbf{K}\mathbf{A}=\boldsymbol{\phi}(\mathbf{X}),$$

(12)

where both $\mathbf{A}$ and $\boldsymbol{\phi}(\mathbf{X})$ are banded matrices. In (12), the $(l,j)^{\rm th}$ entry of $\boldsymbol{\phi}(\mathbf{X})$ is $\phi_{j}(x_{l})$. In view of the compact supportedness of $\phi_{j}$, $\boldsymbol{\phi}(\mathbf{X})$ is a banded matrix with bandwidth $(k-3)/2$:

$$\boldsymbol{\phi}(\mathbf{X})=\begin{bmatrix}\ddots&&&&\\ \ddots&\phi_{j-\frac{k-3}{2}}(x_{j-2\frac{k-3}{2}})&&&\\ \ddots&\vdots&\ddots&&\\ &\phi_{j-\frac{k-3}{2}}(x_{j})&\cdots&\phi_{j+\frac{k-3}{2}}(x_{j})&\\ &&\ddots&\vdots&\ddots\\ &&&\phi_{j+\frac{k-3}{2}}(x_{j+2\frac{k-3}{2}})&\ddots\\ &&&&\ddots\end{bmatrix}.$$

The matrix of $\mathbf{A}$ consists of the coefficients to construct the KPs. In view of the sparse representation, $\mathbf{A}$ is a banded matrix with bandwidth $(k-1)/2$:

$$\mathbf{A}=\begin{bmatrix}\ddots&&&&\\ \ddots&A_{j-2\frac{k-1}{2},j-\frac{k-1}{2}}&&&\\ \ddots&\vdots&\ddots&&\\ &A_{j,j-\frac{k-1}{2}}&\cdots&A_{j,j+\frac{k-1}{2}}&\\ &&\ddots&\vdots&\ddots\\ &&&A_{j+2\frac{k-1}{2},j+\frac{k-1}{2}}&\ddots\\ &&&&\ddots\end{bmatrix}.$$

Computing $\mathbf{A}$ and $\boldsymbol{\phi}(\mathbf{X})$ takes $\mathcal{O}(k^{3}n)$ time, because in the construction of each $\phi_{j}$, at most $k$ kernel basis functions are needed and the time complexity for solving the coefficients $\{A_{w,j}:|w-j|\leq\frac{k-1}{2}\}$ satisfying equation (8), (10) or (11) is $\mathcal{O}(k^{3})$. The computational time $\mathcal{O}(k^{3}n)$ in this step will dominate that in the next step, which is $\mathcal{O}(k^{2}n)$. However, if the design points are equally spaces, the KP coefficients given by (8) will remain the same for each $k$ consecutive data points, so that we only need to compute these values once, and thus the computational time of this step is only $\mathcal{O}(k^{4})$. In this case, the computation time in the next step, i.e., $\mathcal{O}(k^{2}n)$, will be dominant, provided that $k\ll n$.

Now we solve the GP regression problem, by substituting the identity $K(\cdot,\mathbf{X})=\boldsymbol{\phi}(\cdot)\mathbf{A}^{-1}$ into (4) and (5) to obtain

	$\displaystyle\operatorname{\mathbb{E}}\left[Y(x^{*})\big{|}\mathbf{Y}\right]=\mu(x^{*})+\boldsymbol{\phi}^{T}(x^{*})\left[\boldsymbol{\phi}(\mathbf{X})\right]^{-1}\left(\mathbf{Y}-\boldsymbol{\mu}\right),$		(13)
	$\displaystyle\operatorname{\mathrm{Var}}\left[Y(x^{*})\big{|}\mathbf{Y}\right]=\sigma^{2}\bigg{(}K(x^{*},x^{*})-\boldsymbol{\phi}^{T}(x^{*})\left[\boldsymbol{\phi}(\mathbf{X})\right]^{-1}K(\mathbf{X},x^{*})\bigg{)}.$		(14)

The key to GP regression now becomes calculating the vector $\left[\boldsymbol{\phi}(\mathbf{X})\right]^{-1}{\mathbf{v}}$ with $\mathbf{v}=\mathbf{Y}-\boldsymbol{\mu}$ or $\mathbf{v}=K(\mathbf{X},x^{*})$. This is equivalent to solving the sparse banded linear system $\boldsymbol{\phi}(\mathbf{X})\mathbf{s}=\mathbf{v}$. There exists quite a few sparse linear solvers that can solve this linear system efficiently. For example, the algorithm based on the LU decomposition in Davis (2006) can be applied to solve for $\mathbf{s}$ in $\mathcal{O}(k^{2}n)$ time. MATLAB provides convenient and efficient builtin functions, such as mldivide or decomposition, to solve sparse banded linear system in this form.

It is worth noting that (13) can be executed in the following faster way when we need to evaluate $\operatorname{\mathbb{E}}\left[Y(x^{*})|\mathbf{Y}\right]$ for a many different $x^{*}$. First, we compute $\mathbf{s}:=\left[\boldsymbol{\phi}(\mathbf{X})\right]^{-1}\left(\mathbf{Y}-\boldsymbol{\mu}\right)$, which takes $\mathcal{O}(k^{2}n)$ time. Next we evaluate $\mu(x^{*})+\boldsymbol{\phi}^{T}(x^{*})\mathbf{s}$ for different $x^{*}$. As said before, $\boldsymbol{\phi}^{T}(x^{*})$ has at most $k-1$ non-zero entries; see Figure 2. If we know which $k-1$ entries are non-zero, the second step takes only $\mathcal{O}(k)$ time. To find the non-zero entry, a general approach is to use a binary search, which takes $\mathcal{O}(\log n)$ time. Sometime, these entries can be found within a constant time. For example, if the design points are equally spaced, there exist explicit expressions for the indices of the non-zero entries; if we need to predict for $x^{*}$ over a dense mesh (which is a typical task of surrogate modeling), we can use the indices of the non-zero entries for the previous point as an initial guess to find those for the current point.

Similar to the conditional inference, the log-likelihood function (6) can also be computed in $\mathcal{O}(k^{2}n)$ time. First, the log-determinant of $\mathbf{K}$ can be rewritten as $\log\det(\mathbf{K})=\log\det\left(\boldsymbol{\phi}(\mathbf{X})\right)-\log\det(\mathbf{A})$, according to identity (12): Because both $\mathbf{A}$ and $\boldsymbol{\phi}(\mathbf{X})$ are banded matrices, their determinants can be computed in $\mathcal{O}(k^{2}n)$ time by sequential methods (Kamgnia and Nguenang, 2014, section 4.1). Second, the same method for the conditional inference can be applied to compute $\left(\mathbf{Y}-\mathbf{F}\boldsymbol{\beta}\right)^{T}\mathbf{K}^{-1}\left(\mathbf{Y}-\mathbf{F}\boldsymbol{\beta}\right)$ in $\mathcal{O}(k^{2}n)$ time.

Suppose we observe data $\mathbf{Z}$, which is a noisy version of $\mathbf{Y}$. Specifically, $Z(\mathbf{x}_{i})=Y(\mathbf{x}_{i})+\varepsilon$, where $\varepsilon\sim\mathcal{N}(0,\sigma^{2}_{Y})$. In this case, the covariance of the observed noisy responses is $\text{Cov}\big{(}Z(\mathbf{x}_{i}),Z(\mathbf{x}_{j})\big{)}=\sigma^{2}K(\mathbf{x}_{i},\mathbf{x}_{j})+\sigma_{Y}^{2}\mathbb{I}(\mathbf{x}_{i}=\mathbf{x}_{j})$. In other words, the covariance matrix $\text{Cov}(\mathbf{Z},\mathbf{Z})$ is $\sigma^{2}\mathbf{K}+\sigma_{Y}^{2}\mathbb{I}$, where $\mathbb{I}$ is the identity matrix. The posterior predictor at a new point $\mathbf{x}^{*}$ is also normal distributed with the following conditional mean and variance:

	$\displaystyle\operatorname{\mathbb{E}}\left[Y(\mathbf{x}^{*})\big{|}\mathbf{Z}\right]=\mu(\mathbf{x}^{*})+K(\mathbf{x}^{*},\mathbf{X})\left[\mathbf{K}+\frac{\sigma_{Y}^{2}}{\sigma^{2}}\mathbb{I}\right]^{-1}\left(\mathbf{Z}-\boldsymbol{\mu}\right),$		(15)
	$\displaystyle\operatorname{\mathrm{Var}}\left[Y(\mathbf{x}^{*})\big{|}\mathbf{Z}\right]=\sigma^{2}\left(K(\mathbf{x}^{*},\mathbf{x}^{*})-K(\mathbf{x}^{*},\mathbf{X})\left[\mathbf{K}+\frac{\sigma_{Y}^{2}}{\sigma^{2}}\mathbb{I}\right]^{-1}K(\mathbf{X},\mathbf{x}^{*})\right),$		(16)

and the log-likelihood function given data $\mathbf{Z}$ is:

$$L(\boldsymbol{\beta},\sigma^{2},\boldsymbol{\omega})=-\frac{1}{2}\left[\log\det(\sigma^{2}\mathbf{K}+\sigma^{2}_{Y}\mathbb{I})+\big{(}\mathbf{Z}-\mathbf{F}\boldsymbol{\beta}\big{)}^{T}\big{[}\sigma^{2}\mathbf{K}+\sigma_{Y}^{2}\mathbb{I}\big{]}^{-1}\big{(}\mathbf{Z}-\mathbf{F}\boldsymbol{\beta}\big{)}\right].$$

(17)

When the input $\mathbf{x}$ is one dimensional, (15), (16), and (17) can be calculated in $\mathcal{O}(k^{2}n)$ as the noiseless case because the covariance matrix $\sigma^{2}\mathbf{K}+\sigma_{Y}^{2}\mathbb{I}$ admits the following factorization:

$$\sigma^{2}\mathbf{K}+\sigma_{Y}^{2}\mathbb{I}=\big{(}\sigma^{2}\boldsymbol{\phi}(\mathbf{X})+\sigma_{Y}^{2}\mathbf{A}\big{)}\mathbf{A}^{-1}.$$

(18)

By substituting (18) and the identity $K(\cdot,\mathbf{X})=\boldsymbol{\phi}(\cdot)\mathbf{A}^{-1}$ into (15), (16), and (17), we can obtain:

	$\displaystyle\operatorname{\mathbb{E}}\left[Y(x^{*})\big{|}\mathbf{Z}\right]=\mu(x^{*})+\boldsymbol{\phi}^{T}(x^{*})\left[\boldsymbol{\phi}(\mathbf{X})+\frac{\sigma_{Y}^{2}}{\sigma^{2}}\mathbf{A}\right]^{-1}\left(\mathbf{Z}-\boldsymbol{\mu}\right),$		(19)
	$\displaystyle\operatorname{\mathrm{Var}}\left[Y(x^{*})\big{|}\mathbf{Z}\right]=\sigma^{2}\bigg{(}K(x^{*},x^{*})-\boldsymbol{\phi}^{T}(x^{*})\left[\boldsymbol{\phi}(\mathbf{X})+\frac{\sigma_{Y}^{2}}{\sigma^{2}}\mathbf{A}\right]^{-1}K(\mathbf{X},x^{*})\bigg{)}.$		(20)

and

	$\displaystyle L(\boldsymbol{\beta},\sigma^{2},\boldsymbol{\omega})=$	$\displaystyle-\frac{1}{2}\bigg{[}\log\det\big{(}\sigma^{2}\boldsymbol{\phi}(\mathbf{X})+\sigma_{Y}^{2}\mathbf{A}\big{)}-\log\det(\mathbf{A})$
		$\displaystyle+\big{(}\mathbf{Z}-\mathbf{F}\boldsymbol{\beta}\big{)}^{T}\mathbf{A}\big{[}\sigma^{2}\boldsymbol{\phi}(\mathbf{X})+\sigma_{Y}^{2}\mathbf{A}\big{]}^{-1}\big{(}\mathbf{Z}-\mathbf{F}\boldsymbol{\beta}\big{)}\bigg{]}.$		(21)

We have shown that $\boldsymbol{\phi}(\mathbf{X})$ and $\mathbf{A}$ are banded matrices with bandwidth $(k-3)/2$ and $(k-1)/2$, respectively. Therefore, the matrix $\sigma^{2}\boldsymbol{\phi}(\mathbf{X})+\sigma_{Y}^{2}\mathbf{A}$ is also a banded matrix with bandwidth $(k-3)/2$. Time complexity for computing this sum is $\mathcal{O}(kn)$. We then can use the algorithms for banded matrices introduced in section 3.1 to compute (19), (20), and (21) in time complexity $\mathcal{O}(k^{2}n)$. Recall that the time complexities for computing $\boldsymbol{\phi}(\mathbf{X})$ and $\mathbf{A}$ are both $\mathcal{O}(k^{3}n)$. Therefore, in the noisy setting, the total time complexity for computing the posterior and MLE is still $\mathcal{O}(k^{3}n)$, which is the same as the noiseless case.

When data is noiseless, the exact algorithm proposed in Section 3.1 can be used to solve multi-dimensional problems if the input points are full or sparse grids.

A full grid is defined as the Cartesian product of one dimensional point sets: $\mathbf{X}^{\rm FG}=\times_{j=1}^{d}\mathbf{X}^{(j)}$ where each $\mathbf{X}^{(j)}$ denotes any one-dimensional point set. Assuming a separable correlation function (3) comprising $d$ one-dimensional Matérn correlation functions with half-integer smoothness, and inputs on a full grid $\mathbf{X}^{\rm FG}$, the covariance vector $K(\mathbf{x}^{*},\mathbf{X}^{\rm FG})$ and covariance matrix $\mathbf{K}$ decompose into Kronecker products of matrices over each input dimension (Saatçi, 2012; Wilson, 2014):

	$\displaystyle K(\mathbf{x}^{*},\mathbf{X}^{\rm FG})=\bigotimes_{k=1}^{d}K_{j}(x^{*}_{j},\mathbf{X}^{(j)})=\bigotimes_{j=1}^{d}\boldsymbol{\phi}^{T}_{j}(x_{j}^{*})\mathbf{A}^{-1}_{j}=\bigg{(}\bigotimes_{j=1}^{d}\boldsymbol{\phi}^{T}_{j}(x_{j}^{*})\bigg{)}\bigg{(}\bigotimes_{j=1}^{d}\mathbf{A}^{-1}_{j}\bigg{)}$		(22)
	$\displaystyle\mathbf{K}=\bigotimes_{k=1}^{d}K_{j}(\mathbf{X}^{(j)},\mathbf{X}^{(j)})=\bigotimes_{j=1}^{d}\boldsymbol{\phi}_{j}(\mathbf{X}^{(j)})\mathbf{A}^{-1}_{j}=\bigg{(}\bigotimes_{j=1}^{d}\boldsymbol{\phi}_{j}(\mathbf{X}^{(j)})\bigg{)}\bigg{(}\bigotimes_{j=1}^{d}\mathbf{A}^{-1}_{j}\bigg{)}.$		(23)

When we compute the vector $K(\mathbf{x}^{*},\mathbf{X})\mathbf{K}^{-1}$, the matrix $\bigotimes_{j=1}^{d}\mathbf{A}^{-1}_{j}$ is cancelled as the one dimensional case. Therefore, (4) and (5) can be expressed as

	$\displaystyle\operatorname{\mathbb{E}}\left[Y(\mathbf{x}^{*})\big{|}\mathbf{Y}\right]=\mu(\mathbf{x}^{*})+\left(\bigotimes_{j=1}^{d}\boldsymbol{\phi}^{T}_{j}(x_{j}^{*})\right)\left(\bigotimes_{j=1}^{d}\left[\boldsymbol{\phi}_{j}(\mathbf{X}^{(j)})\right]^{-1}\right)\left(\mathbf{Y}-\boldsymbol{\mu}\right)$		(24)
	$\displaystyle\operatorname{\mathrm{Var}}\left[Y(\mathbf{x}^{*})\big{|}\mathbf{Y}\right]=\sigma^{2}\bigg{(}K(\mathbf{x}^{*},\mathbf{x}^{*})-\prod_{j=1}^{d}\boldsymbol{\phi}^{T}_{j}(x_{j}^{*})\left[\boldsymbol{\phi}_{j}(\mathbf{X}^{(j)})\right]^{-1}K_{j}(\mathbf{X}^{(j)},x_{j}^{*})\bigg{)}$		(25)

and the log-likelihood function (6) becomes

	$\displaystyle L(\boldsymbol{\beta},\sigma^{2},\boldsymbol{\omega})$	$\displaystyle=-\frac{1}{2}\bigg{[}n\log\sigma^{2}+\sum_{j=1}^{d}\frac{n}{n_{j}}\left(\log\det\boldsymbol{\phi}_{j}(\mathbf{X}^{(j)})-\log\det\mathbf{A}_{j}\right)$
		$\displaystyle+\frac{1}{\sigma^{2}}\big{(}\mathbf{Y}-\mathbf{F}\boldsymbol{\beta}\big{)}^{T}\left(\bigotimes_{j=1}^{d}\mathbf{A}_{j}\right)\left(\bigotimes_{j=1}^{d}\left[\boldsymbol{\phi}_{j}(\mathbf{X}^{(j)})\right]^{-1}\right)\big{(}\mathbf{Y}-\mathbf{F}\boldsymbol{\beta}\big{)}\bigg{]},$		(26)

where $\boldsymbol{\phi}_{j}$’s are the KPs associated to correlation function $K_{j}$ and point set $\mathbf{X}^{(j)}$, $\mathbf{A}_{j}$ is the coefficient matrix for constructing $\boldsymbol{\phi}_{j}$ defined in (12), and $n=\prod_{j=1}^{d}n_{j}$ is the size of $\mathbf{X}^{\rm FG}$, $n_{j}$ is the size of $\mathbf{X}^{(j)}$. We can also note that entries of vector $\bigotimes_{j=1}^{d}\boldsymbol{\phi}_{j}(\cdot)$ are products of one-dimensional KPs. Therefore, similar to the one-dimensional case, $\bigotimes_{j=1}^{d}\boldsymbol{\phi}_{j}(\cdot)$ is a vector of compactly supported functions.