A General $\mathcal{O}(n^2)$ Hyper-Parameter Optimization for Gaussian Process Regression with Cross-Validation and Non-linearly Constrained ADMM

TL;DR

This paper introduces two novel hyper-parameter optimization schemes for Gaussian process regression that achieve quadratic complexity using cross-validation and ADMM, outperforming traditional methods in accuracy and efficiency.

Contribution

It proposes general-purpose GP hyper-parameter training schemes replacing maximum likelihood with cross-validation and gradient descent with ADMM, reducing complexity to O(n^2).

Findings

Schemes demonstrate excellent convergence and accuracy.

Significant reduction in computational time.

Effective on both synthetic and real datasets.

Abstract

Hyper-parameter optimization remains as the core issue of Gaussian process (GP) for machine learning nowadays. The benchmark method using maximum likelihood (ML) estimation and gradient descent (GD) is impractical for processing big data due to its $O(n^3)$ complexity. Many sophisticated global or local approximation models, for instance, sparse GP, distributed GP, have been proposed to address such complexity issue. In this paper, we propose two novel and general-purpose GP hyper-parameter training schemes (GPCV-ADMM) by replacing ML with cross-validation (CV) as the fitting criterion and replacing GD with a non-linearly constrained alternating direction method of multipliers (ADMM) as the optimization method. The proposed schemes are of $O(n^2)$ complexity for any covariance matrix without special structure. We conduct various experiments based on both synthetic and real data sets, wherein the proposed schemes show excellent performance in terms of convergence, hyper-parameter estimation accuracy, and computational time in comparison with the traditional ML based routines given in the GPML toolbox.

Tables

Table 1. Table 1: Quantitative comparisons between GPCV-ADMM and GPML across nine synthetic data sets (combining three kernels and three data lengths). We recorded 1) the sample mean of hyper-parameter estimates and 2) the sample mean of the MSE averaged over 50 Monte Carlo simulations.

			Hyper-paramete Estimates (std)		Test MSE
Kernel	Setting	Size	GPCV-ADMM	GPML	GPCV-ADMM	GPML
SE	$l=0.5$	500	[0.46(0.054)]	[0.52(0.018)]	0.12	0.13
		1000	[0.50(0.044)]	[0.50(0.017)]	0.12	0.14
		2000	[0.50(0.017)]	[0.53(0.020)]	0.12	0.14
LP	$l=0.5$	500	[0.34(0.023),1.13(0.080)]	[0.55(0.052),1.06(0.090)]	0.13	0.36
	$p=1$	1000	[0.39(0.014),1.06(0.063)]	[0.52(0.058),1.16(0.130)]	0.17	0.44
		2000	[0.44(0.082),1.19(0.008)]	[0.53(0.013),1.02(0.019)]	0.26	0.28
SE+LP	$l_1=3$	500	[3.73(0.280),0.94(0.100),2.38(0.150)]	[3.62(0.630), 1.08(0.210),2.38(0.860)]	0.18	0.21
	$l_2=1$	1000	[3.74(0.250),0.95(0.110),2.14(0.090)]	[3.61(0.560),1.05(0.180),2.27(0.600)]	0.13	0.15
	$l_3=2$	2000	[3.94(0.075),0.99(0.120),2.10(0.080)]	[3.69(0.390),1.09(0.130),2.13(0.290)]	0.34	0.37

Equations

$$f(\mathbf{x})\sim \mathcal{G}\mathcal{P}(m(\mathbf{x}),k(\mathbf{x},{\mathbf{x}}^{\prime };{\mathbf{\theta }}_h)),$$

$${\mathbf{\theta }}_{ML}\triangleq \arg \underset{\mathbf{\theta }}{\min }l(\mathbf{\theta })={\mathbf{y}}^T{\mathbf{C}}^{-1}(\mathbf{\theta })\mathbf{y}+\log \det \left(\mathbf{C}(\mathbf{\theta })\right),$$

$${\theta }_i^{k+1}={\theta }_i^k-\mu \cdot \frac{\partial l(\mathbf{\theta })}{\partial {\theta }_i}{|}_{\mathbf{\theta }={\mathbf{\theta }}^k},\forall i=1,2,...,p,$$

$$\frac{\partial l(\mathbf{\theta })}{\partial {\theta }_i}=\text{tr}\left\{\left({\mathbf{C}}^{-1}(\mathbf{\theta })-\mathbf{\gamma }{\mathbf{\gamma }}^T\right)\frac{\partial \mathbf{C}(\mathbf{\theta })}{\partial {\theta }_i}\right\},$$

$$\bar{\mathbf{m}}$$

$$\bar{\mathbf{V}}$$

$$-\mathbf{K}({\mathbf{X}}_{\ast },\mathbf{X}){\left[\mathbf{K}(\mathbf{X},\mathbf{X})+{\sigma }_e^2{\mathbf{I}}_n\right]}^{-1}\mathbf{K}(\mathbf{X},{\mathbf{X}}_{\ast }).$$

$$\bar{\mathbf{m}}({\mathbf{X}}_V;\mathbf{\theta })\triangleq \mathbf{K}({\mathbf{X}}_V,{\mathbf{X}}_T;{\mathbf{\theta }}_h){\mathbf{z}}_T,$$

$${\mathbf{z}}_T\triangleq {\left[\mathbf{K}({\mathbf{X}}_T,{\mathbf{X}}_T;{\mathbf{\theta }}_h)+{\sigma }_e^2{\mathbf{I}}_n\right]}^{-1}{\mathbf{y}}_T.$$

$${\mathbf{\theta }}_{CV}=\arg \underset{\mathbf{\theta }}{\min }{\left|\left|{\mathbf{y}}_V-\bar{\mathbf{m}}({\mathbf{X}}_V;\mathbf{\theta })\right|\right|}_2^2.$$

$$\mathbf{C}(\mathbf{\theta }){\mathbf{z}}_T={\mathbf{y}}_T.$$

$$L_{\rho }\left(\mathbf{\theta },{\mathbf{z}}_T,\mathbf{\lambda }\right)$$

$$+{\mathbf{\lambda }}^T\left(\mathbf{C}(\mathbf{\theta }){\mathbf{z}}_T-{\mathbf{y}}_T\right)+\frac{\rho }{2}{\left|\left|\mathbf{C}(\mathbf{\theta }){\mathbf{z}}_T-{\mathbf{y}}_T\right|\right|}_2^2,$$

$${\mathbf{\theta }}^{\eta +1}$$

$${\mathbf{z}}_T^{\eta +1}$$

$${\mathbf{\lambda }}^{\eta +1}$$

$${\mathbf{\theta }}^{\eta +1}$$

$$\equiv {\mathbf{\theta }}^{\eta }-{\mu }_1\cdot {\nabla }_{\mathbf{\theta }}{\mathbf{g}}^{(\eta )}\left(\mathbf{\theta }\right){\mid }_{\mathbf{\theta }={\mathbf{\theta }}^{\eta }},$$

$$\arg \underset{{\mathbf{z}}_T}{\min }{\mathbf{g}}^{(\eta )}\left({\mathbf{z}}_T\right)=({\mathbf{b}}^{\eta }{)}^T{\mathbf{z}}_T+{\mathbf{z}}_T^T{\mathbf{\Sigma }}^{\eta }{\mathbf{z}}_T,$$

$${\mathbf{b}}^{\eta }$$

$${\mathbf{\Sigma }}^{\eta }$$

$${\mathbf{z}}_T^{\eta +1}={\mathbf{z}}_T^{\eta }+{\mu }_2\cdot {\mathbf{d}}_z^{\eta },$$

$${\mathbf{\theta }}_{CV}=\arg \underset{\mathbf{\theta }}{\min }\sum _{k=1}^K{\left|\left|{\mathbf{y}}_V^{(k)}-\bar{\mathbf{m}}({\mathbf{X}}_V^{(k)};\mathbf{\theta })\right|\right|}_2^2.$$

$${\mathbf{\theta }}_{CV}$$

$${\mathbf{\theta }}_1={\mathbf{\theta }}_2=,...,={\mathbf{\theta }}_K=\mathbf{z},$$

$${\mathbf{g}}^{(\eta )}\left(\mathbf{\theta }\right)$$

$$+({\mathbf{z}}_T^{\eta }{)}^T{\mathbf{K}}_{VT}({\mathbf{\theta }}_h{)}^T{\mathbf{K}}_{VT}({\mathbf{\theta }}_h){\mathbf{z}}_T^{\eta }$$

$$+{\mathbf{\lambda }}^T{\mathbf{K}}_{TT}({\mathbf{\theta }}_h){\mathbf{z}}_T^{\eta }$$

$$+\rho ({\sigma }_e^2{\mathbf{z}}_T^{\eta }-{\mathbf{y}}_T{)}^T{\mathbf{K}}_{TT}({\mathbf{\theta }}_h){\mathbf{z}}_T^{\eta }$$

$$+\frac{\rho }{2}({\mathbf{z}}_T^{\eta }{)}^T{\mathbf{K}}_{TT}({\mathbf{\theta }}_h){\mathbf{K}}_{TT}({\mathbf{\theta }}_h){\mathbf{z}}_T^{\eta }.$$

$$\frac{\partial {\mathbf{g}}^{(\eta )}\left(\mathbf{\theta }\right)}{\partial {\theta }_i}$$

$$+({\mathbf{z}}_T^{\eta }{)}^T\frac{\partial {\mathbf{K}}_{VT}^T({\mathbf{\theta }}_h)}{\partial {\theta }_i}{\mathbf{K}}_{VT}({\mathbf{\theta }}_h){\mathbf{z}}_T^{\eta }$$

$$+({\mathbf{z}}_T^{\eta }{)}^T{\mathbf{K}}_{VT}^T({\mathbf{\theta }}_h)\frac{\partial {\mathbf{K}}_{VT}({\mathbf{\theta }}_h)}{\partial {\theta }_i}{\mathbf{z}}_T^{\eta }$$

$$+{\mathbf{\lambda }}^T\frac{\partial {\mathbf{K}}_{TT}({\mathbf{\theta }}_h)}{\partial {\theta }_i}{\mathbf{z}}_T^{\eta }$$

$$+\rho ({\sigma }_e^2{\mathbf{z}}_T^{\eta }-{\mathbf{y}}_T{)}^T\frac{\partial {\mathbf{K}}_{TT}({\mathbf{\theta }}_h)}{\partial {\theta }_i}{\mathbf{z}}_T^{\eta }$$

$$+\frac{\rho }{2}({\mathbf{z}}_T^{\eta }{)}^T{\mathbf{K}}_{TT}({\mathbf{\theta }}_h)\frac{\partial {\mathbf{K}}_{TT}({\mathbf{\theta }}_h)}{\partial {\theta }_i}{\mathbf{z}}_T^{\eta }$$

$$+\frac{\rho }{2}({\mathbf{z}}_T^{\eta }{)}^T\frac{\partial {\mathbf{K}}_{TT}({\mathbf{\theta }}_h)}{\partial {\theta }_i}{\mathbf{K}}_{TT}({\mathbf{\theta }}_h){\mathbf{z}}_T^{\eta }.$$

$$k_{se}(x,x^{\prime })={\sigma }^{2}exp\left(-\frac{(x-x^{\prime }{)}^2}{2l^2}\right)$$

$$k$$

Full text

1 Introduction

Gaussian process (GP) models constitute a class of important Bayesian non-parametric models for machine learning and are tightly connected to several other salient models, such as support vector machines (SVM), single-layer Bayesian neural networks, regularized-least-squares, relevance vector machines, auto-regressive-moving-average (ARMA) and deep neural networks Neal (1996); MacKay (1998); Rasmussen and Williams (2006); Matthews et al. (2018). The idea behind GP models is to impose a Gaussian prior on the underlying function and then compute the predictive distribution over the function given the observed data. Due to their outstanding performance in function approximation with a natural uncertainty bound, GP models have been adopted in a plethora of applications, e.g., nonlinear system identification Frigola and Rasmussen (2013), state-space modeling and trajectory prediction Frigola et al. (2014), financial data modeling and prediction Han et al. (2016), to mention a few.

The predictive performance of GP regression models depends on the goodness of kernel selection and hyper-parameter estimation. There exist two classes of methods for GP hyper-parameter training. The first class of deterministic methods includes the maximum likelihood (ML) estimation based methods and the cross-validation (CV) based methods Rasmussen and Williams (2006); Krauth et al. (2017) among the others. The second class of stochastic methods includes for instance, the hybrid Monte-Carlo and Markov chain Monte-Carlo sampling methods Neal (1997); Hensman et al. (2015); Havasi et al. (2018). In this paper, we focus on the first class of deterministic methods.

The standard GP methods adopt the ML criterion to optimize the GP hyper-parameters via gradient descent (GD). However, the computational complexity scales cubically over the number of data points, i.e., $O(n^{3})$, which is deemed as the most prominent weakness of such methods. Low-complexity GP methods could be obtained with the ideas to 1) find a smaller subset (with $m\ll n$ data points) of the complete data set and construct a sparse representation of the original covariance matrix, yielding $\mathcal{O}(m^{2}n)$ complexity Quiñonero Candela and Rasmussen (2005); Titsias (2009); Wilson and Nickisch (2015); 2) adopt low-rank approximations of the covariance matrices, e.g., via hierarchical factorization of the covariance matrix into a product of block low-rank updates of the identity matrix, yielding $\mathcal{O}(n\log^{2}n)$ complexity Ambikasaran et al. (2016); 3) employ a number of $K$ local computing units, build GP models in smaller scales with a subset of data at each local computing unit, and fuse the hyper-parameter estimates via Bayesian committee machine type of approaches, yielding $\mathcal{O}(n^{3}/K^{3})$ complexity Deisenroth and Ng (2015). All these results rely on the use of ML criterion with different types of approximations. A more complete list of low-complexity methods were surveyed in Liu et al. (2018).

Cross-validation is yet another nice criterion for hyper-parameter training, demonstrated as an alternative of the ML criterion for model selection and GP hyper-parameter training in Rasmussen and Williams (2006). The idea of using CV for GP hyper-parameter training is perhaps more desirable than using the ML criterion due to the following two reasons. First, for Gaussian noises, the CV based method is essentially equivalent to the ML based method. Second, for general noises, the CV based method directly minimizes the mean-squared-error (MSE), which is widely used as the ultimate performance metric for regression tasks.

The main contributions of this paper are as follows. We propose two novel CV based schemes for GP hyper-parameter optimization. In order to eliminate the large-scale matrix inverse, an auxiliary variable is newly introduced and optimized jointly with the GP hyper-parameters using alternating direction method of multipliers (ADMM), where ADMM was first introduced in the mid-1970s and extended to handle a wide range of optimization problems in machine learning, control, and signal processing nowadays Eckstein and Bertsekas (1992); Bertsekas and Tsitsiklis (1997); Boyd et al. (2011); Hong et al. (2016). The proposed schemes demonstrate $\mathcal{O}(n^{2})$ computational complexity without making any approximation. Various experimental results show even better and faster GP hyper-parameter training results of the proposed schemes as compared to the GPML toolbox111Available on http://www.gaussianprocess.org/gpml/.

The remainder of the paper is organized as follows. In Section 2, we review the standard GP regression that applies ML and GD to train the hyper-parameters, yielding $\mathcal{O}(n^{3})$ complexity. In Section 3, we introduce two new schemes based on CV and ADMM, whose $\mathcal{O}(n^{2})$ complexity is validated in Section 4. Experimental results are given in Section 5. Finally, we conclude the paper in Section 6.

2 Standard Gaussian Process Regression

A Gaussian process is a collection of random variables, whose any finite subset follows a Gaussian distribution Rasmussen and Williams (2006). In the sequel, we consider the scalar output, real valued Gaussian processes that are completely specified by a mean function and a kernel function (a.k.a. covariance function) as:

[TABLE]

where $m(\boldsymbol{x})$ denotes the mean function and $k(\boldsymbol{x},\boldsymbol{x}^{\prime};\boldsymbol{\theta}_{h})$ represents the kernel function specified by a vector of hyper-parameters, $\boldsymbol{\theta}_{h}$.

The most widely used GP regression model is given by $y=f(\boldsymbol{x})+e$, where $y\in\mathbb{R}$ is a continuous valued, scalar output. The unknown function $f(\boldsymbol{x}):\mathbb{R}^{d}\mapsto\mathbb{R}$ is modeled via a zero mean Gaussian process for simplicity. The noise term $e$ is assumed to be Gaussian distributed with zero mean and variance $\sigma_{e}^{2}$, independent and identically distributed at different data points. The vector of unknown GP hyper-parameters is denoted by $\boldsymbol{\theta}\triangleq[\boldsymbol{\theta}_{h}^{T},\sigma^{2}_{e}]^{T}$, whose dimension is assumed to be $p$.

In the training phase, we are given a data set $\mathcal{D}\triangleq\{\boldsymbol{X},\boldsymbol{y}\}$, where $\boldsymbol{y}=[y_{1},y_{2},...,y_{n}]^{T}$ is the set of outputs and $\boldsymbol{X}=[\boldsymbol{x}_{1},\boldsymbol{x}_{2},...,\boldsymbol{x}_{n}]$ is the set of inputs. As the benchmark, we review the classic maximum likelihood (ML) based GP hyper-parameter optimization. Due to the Gaussian assumption on the noise term, the log-likelihood function could be obtained in closed form, and the GP hyper-parameters could be trained equivalently by minimizing the negative log-marginal likelihood function (ignoring the unrelated terms) as:

[TABLE]

where $\boldsymbol{C}(\boldsymbol{\theta})\triangleq\boldsymbol{K}(\boldsymbol{X},\boldsymbol{X};\boldsymbol{\theta}_{h})+\sigma_{e}^{2}\boldsymbol{I}_{n}$. The GD type of methods are most widely used for hyper-parameter optimization in the GP community. At the $(k+1)$-th iteration, each element of the GP hyper-parameters is updated as:

[TABLE]

where $\mu$ is a positive step size. The partial derivative of the $i$-th element in $\boldsymbol{\theta}$ can be derived in closed form as:

[TABLE]

where $\boldsymbol{\gamma}\triangleq\boldsymbol{C}^{-1}(\boldsymbol{\theta})\boldsymbol{y}$ is defined for notational brevity. At each iteration, $\boldsymbol{C}^{-1}(\boldsymbol{\theta})$ needs to be re-evaluated with the updated $\boldsymbol{\theta}=\boldsymbol{\theta}^{k}$, where multiplication of $n\times n$ matrices has to be performed for several times. Therefore, the total computational complexity scales as $\mathcal{O}(n^{3})$ per iteration, in general.

In the test phase, the trained GP models compute the predictive distribution of $\boldsymbol{y}_{*}=[y_{*,1},y_{*,2},...,y_{*,n_{*}}]^{T}$ for novel test inputs $\boldsymbol{X}_{*}=[\boldsymbol{x}_{*,1},\boldsymbol{x}_{*,2},...,\boldsymbol{x}_{*,n_{*}}]$. The test data set is denoted as $\mathcal{D}_{*}=\{\boldsymbol{y}_{*},\boldsymbol{X}_{*}\}$. According to the definition of Gaussian process and applying the canonical results of conditional Gaussian distribution, we could easily derive the predictive distribution to be $p(\boldsymbol{y}_{*}|\mathcal{D},\boldsymbol{X}_{*};\boldsymbol{\theta})\sim\mathcal{N}\left(\bar{\boldsymbol{m}},\bar{\boldsymbol{V}}\right)$, where the mean and variance are:

[TABLE]

Here, $\boldsymbol{K}(\boldsymbol{X},\boldsymbol{X})$ is an $n\times n$ matrix of covariances among the training inputs; $\boldsymbol{K}(\boldsymbol{X},\boldsymbol{X}_{*})$ is an $n\times n_{*}$ matrix of covariances between the training inputs and test inputs; $\boldsymbol{K}(\boldsymbol{X}_{*},\boldsymbol{X}_{*})$ is an $n_{*}\times n_{*}$ matrix of covariances among the test inputs.

3 Proposed Cross-Validation based GP Hyper-parameter Optimization Schemes

In this section, we introduce two new GP hyper-parameter optimization schemes by replacing ML with CV, and replacing GD with ADMM. The proposed schemes achieve $O(n^{2})$ computational complexity. We assume that the number of training data samples, $n$, is moderately large such that a single computing unit could handle computations in the order of $\mathcal{O}(n^{2})$, while $\mathcal{O}(n^{3})$ is beyond the processing limit.

3.1 Hold-out Cross-Validation based Scheme

We first focus on the training phase and consider the simple hold-out cross-validation (HOCV), which is widely used for model selection in machine learning. We divide the data set $\mathcal{D}$ into two non-overlapping subsets, namely the training set $\mathcal{D}_{T}=\{\boldsymbol{X}_{T},\boldsymbol{y}_{T}\}$ with $|\mathcal{D}_{T}|=n_{t}$ samples and the validation set $\mathcal{D}_{V}=\{\boldsymbol{X}_{V},\boldsymbol{y}_{V}\}$ with $|\mathcal{D}_{V}|=n_{v}$ samples. Note that, the ML based scheme instead use the whole $\mathcal{D}$ for training the GP hyper-parameters. We illustrate this difference with an example in Figure 1. In practice, we could pick $n_{v}$ samples randomly from $\mathcal{D}$ for validation and leave the rest for training. The overview of Gaussian process regression in Section 2 tells us that the predictive mean of the validation points in $\mathcal{D}_{V}$ given the training data set $\mathcal{D}_{T}$ is:

[TABLE]

where

[TABLE]

In the sequel, we use the short forms $\boldsymbol{K}_{VT}(\boldsymbol{\theta}_{h})$, $\boldsymbol{K}_{TT}(\boldsymbol{\theta}_{h})$ and $\boldsymbol{C}(\boldsymbol{\theta})$ to denote $\boldsymbol{K}(\boldsymbol{X}_{V},\boldsymbol{X}_{T};\boldsymbol{\theta}_{h})$, $\boldsymbol{K}(\boldsymbol{X}_{T},\boldsymbol{X}_{T};\boldsymbol{\theta}_{h})$, and $\boldsymbol{K}(\boldsymbol{X}_{T},\boldsymbol{X}_{T};\boldsymbol{\theta}_{h})+\sigma_{e}^{2}\boldsymbol{I}_{n}$ respectively.

We aim to find an optimal vector of the GP hyper-parameters that minimizes the differences between the validation outputs and their predictive means, which is formulated as:

[TABLE]

This optimization problem is non-convex in terms of $\boldsymbol{\theta}$ for most kernels. When applying GD for solving (9), the inverse of a possibly large $n\times n$ covariance matrix $\boldsymbol{C}(\boldsymbol{\theta})$ has to be computed at each iteration with $\mathcal{O}(n^{3})$ complexity in general, which forbids the practical use of Gaussian process regression for big data applications. In order to remedy this drawback, we let $\boldsymbol{z}_{T}$ be a vector-formed auxiliary variable and impose the following nonlinear equality constraint:

[TABLE]

In order to adopt ADMM for GP hyper-parameter optimization with $\mathcal{O}(n^{2})$ complexity, we formulate the augmented Lagrangian function as:

[TABLE]

where the regularization parameter $\rho\geq 0$ is pre-selected and $\boldsymbol{\lambda}$ is a vector of Lagrange multipliers. The complete method consists of a $\boldsymbol{\theta}$-minimization step, a $\boldsymbol{z}_{T}$-minimization step, and a closed form dual variable update step. Concretely, at the $(\eta+1)$-th iteration we have,

[TABLE]

We elaborate on the $\boldsymbol{\theta}$-minimization step first. Note that $L_{\rho}(\boldsymbol{\theta},\boldsymbol{z}_{T}^{\eta},\boldsymbol{\lambda}^{\eta})$ is often a non-convex function in terms of $\boldsymbol{\theta}$, which could be minimized as:

[TABLE]

where $\mu_{1}$ is a positive step size, $\boldsymbol{g}^{(\eta)}\left(\boldsymbol{\theta}\right)$ and its partial derivatives are introduced in the Appendix. When $\mu_{1}$ is chosen via line minimization or the Armijo’s rule, the sequence $\{\boldsymbol{\theta}^{\eta}\}$ is guaranteed to converge to a stationary point of subproblem (11a) according to Proposition 1.2.1 of Bertsekas (2016).

Second, we elaborate on the $\boldsymbol{z}_{T}$-minimization step. In subproblem (11b), it is easy to verify that $L_{\rho}(\boldsymbol{\theta}^{\eta+1},\boldsymbol{z}_{T},\boldsymbol{\lambda}^{\eta})$ is a quadratic function of $\boldsymbol{z}_{T}$ and minimizing it with respect to $\boldsymbol{z}_{T}$ is equivalent to

[TABLE]

where

[TABLE]

It is easy to verify that $\boldsymbol{\Sigma}^{\eta}$ is always positive definite. However, taking the derivative of $\boldsymbol{g}^{(\eta)}\left(\boldsymbol{z}_{T}\right)$ with respect to $\boldsymbol{z}_{T}$ and setting it equal to zero for closed form solution involves the inverse of $\boldsymbol{\Sigma}^{\eta}$ with $\mathcal{O}(n^{3})$ complexity. Thus, we instead solve the quadratic minimization problem in (13) numerically via the conjugate gradient method (CGM) Bertsekas (2016). Specifically, we update the auxiliary variable as

[TABLE]

where $\mu_{2}$ is a stepsize obtained via line minimization, the descent direction $\boldsymbol{d}^{\eta}_{z}$ is an outcome of Gram-Schmidt orthogonization process: $\boldsymbol{d}^{\eta}_{z}=-\boldsymbol{g}^{\eta}_{z}+\beta_{r}\boldsymbol{d}^{\eta-1}_{z}$ with $\boldsymbol{g}^{\eta}_{z}$ being short for $\nabla_{\boldsymbol{z}_{T}}g(\boldsymbol{z}^{\eta}_{T})=2\boldsymbol{\Sigma}^{\eta}\boldsymbol{z}^{\eta}_{T}+\boldsymbol{b}^{\eta}$ and $\beta_{r}=||\boldsymbol{g}^{\eta}_{z}||_{2}^{2}/||\boldsymbol{g}^{\eta-1}_{z}||_{2}^{2}$ , which starts from $\boldsymbol{d}^{0}_{z}=-\boldsymbol{g}^{0}_{z}=-\nabla_{\boldsymbol{z}_{T}}g(\boldsymbol{z}^{0}_{T})$ and terminates at an optimal solution after at most $n_{T}$ steps. Recall that we introduce an auxiliary variable $\boldsymbol{z}_{T}$ to eliminate the large matrix inverses, fortunately the efficient CGM exists for the solution with little computational effort.

Lastly, the update of $\boldsymbol{\lambda}^{\eta+1}$ is conducted in light of (11c) after $\boldsymbol{\theta}^{\eta+1}$ and $\boldsymbol{z}_{T}^{\eta+1}$ are obtained. For the HOCV based scheme, we have the following theorem (for more details see Sections 4.2 and 5.2 of Bertsekas (2016)):

Theorem 1**.**

Local Convergence Property: When taking the initial guess $\boldsymbol{\lambda}^{0}$ close to the optimal Lagrange multiplier $\boldsymbol{\lambda}^{*}$ and taking $\rho$ large enough, solving the unconstrained minimization problem $L_{\rho}(\boldsymbol{\theta},\boldsymbol{z}_{T},\boldsymbol{\lambda})$ can yield points close to the local minimum $\boldsymbol{\theta}^{*}$ and $\boldsymbol{z}_{T}^{*}$ that satisfy the sufficient optimality conditions.

3.2 $K$-fold Cross Validation Based Scheme

In this subsection, we aim to design a $K$-fold CV based GP hyper-parameter optimization scheme, which is able to generate more robust result and exploit parallel computing. We let $\mathcal{D}_{V}^{(k)}=\{\boldsymbol{y}_{V}^{(k)},\boldsymbol{X}_{V}^{(k)}\}$ be the $k$-th partition of the complete data set $\mathcal{D}$ to be used for validation, and the corresponding training set $\mathcal{D}_{T}^{(k)}=\mathcal{D}\backslash\mathcal{D}_{V}^{(k)}$.

A Naive Scheme:

We train the GP hyper-parameters $\boldsymbol{\theta}_{CV}^{(k)}$, heuristically for every partition $k=1,2,...,K$, using the same routine given in the HOCV based scheme (see Algorithm 1), and average the results to obtain the final estimate $\boldsymbol{\theta}_{CV}=1/K\cdot\sum_{k=1}^{K}\boldsymbol{\theta}_{CV}^{(k)}$.

A Principled Scheme:

Alternatively, we formulate an optimization problem for the same purpose but with a sound rationale. For $K\geq 2$, the final estimate of the GP hyper-parameters is obtained as:

[TABLE]

To tackle (15), we introduce some local copies of $\boldsymbol{\theta}$ and solve the following linear equality-constrained optimization problem:

[TABLE]

where $l_{k}(\boldsymbol{\theta}_{k})\triangleq\left|\left|\boldsymbol{y}_{V}^{(k)}-\bar{\boldsymbol{m}}(\boldsymbol{X}_{V}^{(k)};\boldsymbol{\theta}_{k})\ \right|\right|_{2}^{2}$ is non-convex in terms of $\boldsymbol{\theta}_{k}$ for most kernels. The following theorem from Hong et al. (2016) is valuable and supportive:

Theorem 2**.**

When the following assumptions hold: 1) $l_{k}(\boldsymbol{\theta}_{k})$ satisfies the Lipschitz condition; 2) the augmented Lagrangian parameter, $\rho$, is chosen large enough; 3) the minimization problems are bounded from below and all ADMM subproblems are solved exactly, etc., it is guaranteed that any limit point of problem (16) is also a stationary solution.

4 Computational Complexity

We aim to verify the computational complexity of the proposed GP hyper-parameter optimization schemes in Section 3.

We start with the HOCV based scheme. Updating one particular element of the hyper-parameters (out of $p$ elements), say $\theta_{i}$, according to (12), mainly involves the computations of $\frac{\partial\boldsymbol{K}_{VT}(\boldsymbol{\theta}_{h})}{\partial\theta_{i}}\boldsymbol{z}_{T}^{\eta}$, and $\boldsymbol{K}_{VT}(\boldsymbol{\theta}_{h})\boldsymbol{z}_{T}^{\eta}$, $\frac{\partial\boldsymbol{K}_{TT}(\boldsymbol{\theta}_{h})}{\partial\theta_{i}}\boldsymbol{z}_{T}^{\eta}$, $\boldsymbol{K}_{TT}(\boldsymbol{\theta}_{h})\boldsymbol{z}_{T}^{\eta}$ and some cheap vector inner products. Therefore, the computational complexity scales as $\mathcal{O}(n_{v}\cdot n_{t}+n_{t}^{2})=\mathcal{O}(n\cdot n_{t})$ for this step. Similarly, updating the auxiliary parameter $\boldsymbol{z}_{T}$, according to (14), mainly involves the computations of $\boldsymbol{K}_{VT}(\boldsymbol{\theta}_{h}^{\eta+1})\boldsymbol{z}_{T}^{\eta}$, $\boldsymbol{C}(\boldsymbol{\theta}^{\eta+1})\boldsymbol{z}_{T}^{\eta}$ as well as some cheap vector inner products (for details see the Appendix), thus the computational complexity also scales as $\mathcal{O}(n_{v}\cdot n_{t}+n_{t}^{2})=\mathcal{O}(n\cdot n_{t})$. The third step involves only a closed form update, whose complexity scales as $\mathcal{O}(n_{t}^{2})$. Therefore, the overall computational complexity for running one complete ADMM iteration scales as $\mathcal{O}(p\cdot n\cdot n_{t})\approx\mathcal{O}(n^{2})$ for $p\ll n$, which is much lower than $\mathcal{O}(n^{3})$.

For the $K$-fold CV based scheme, each computing unit updates a local copy of the global variable, $\boldsymbol{\theta}_{k}$, incurring $\mathcal{O}(n^{2})$ complexity according to the above analysis. The overall complexity remains low for practical $K$ and $p$ values.

5 Experiments

In this section, we aim to evaluate our proposed schemes, termed as GPCV-ADMM, using both synthetic and real data sets. Specifically, we choose the naive two-fold CV based scheme, which is more practical to implement. As the benchmark, we choose the state-of-the-art ML based scheme given in the GPML toolbox.

Simulation Platform:

Our GPCV-ADMM is implemented in R (version 3.5.2), and compared with the GPML toolbox executed in MATLAB 2018b. All the experiments were conducted on a MacBook Pro with 2.2 GHz Intel Core i7.

Algorithmic Setup:

For both fairness and clarity of comparisons, all the hyper-parameters are initialized with same values for both GPCV-ADMM and GPML when testing on the synthetic data sets under a fixed kernel configuration. However, random restarts are recommended for initialization in practice, and was also adopted in our experiment for the real $CO_{2}$ concentration data set. The auxiliary variable $\boldsymbol{z}_{T}$ of GPCV-ADMM is initialized according to (10) with a perturbed $\boldsymbol{\theta}^{0}$, and the dual variable $\boldsymbol{\lambda}$ initialized to be a vector of all ones. The regularization parameter is pre-selected to be $\rho=5$. The error tolerance for ADMM is set to be $\epsilon=10^{-2}$ and the maximum number of iterations is set to be $maxItr=100$. The source code is available online222 Available on https://github.com/eveneveno/GPCV-ADMM-2019.

5.1 Data Sets

Synthetic Data:

The synthetic data sets used in our experiments were generated from: 1) a squared exponential (SE) kernel with $\boldsymbol{\theta}=l$, 2) a local periodic (LP) kernel with $\boldsymbol{\theta}=[l,p]^{T}$, and 3) a composite kernel by adding up a SE and LP kernel (short as SE+LP in the sequel) with $\boldsymbol{\theta}=[l_{1},l_{2},p]^{T}$. Explicit forms of these selected kernels and their interpretations are given in the Appendix. For each selected kernel, we generated three kinds of data sets with sample sizes $n=500,1000,2000$, respectively, for the primary purpose of verifying the $O(n^{2})$ complexity of GPCV-ADMM as compared to the $O(n^{3})$ complexity of GPML. The synthetic data inputs lie in the range of $[0,10]$ for $n=500$, while in the range of $[0,20]$ for $n=1000,2000$. For a fixed kernel configuration and a fixed sample size, we ran 50 independent Monte Carlo trials to evaluate the goodness of the hyper-parameter estimates obtained by the two schemes. Note that the two schemes use the full data sets differently for training the GP hyper-parameters, as is shown in Figure 1 for clarity.

Real Atmospheric $CO_{2}$ Concentration Data:

This data set consists of the monthly average atmospheric $CO_{2}$ concentrations at the Mauna Loa Observatory, Hawaii, from 1958 to 2015, and is widely used in the GP community333Available on http://scrippsco2.ucsd.edu/data/atmospheric_co2. In Rasmussen and Williams (2006), the authors proposed to use a composite kernel that consists of four parts: 1) an SE kernel modeling the rising trend, 2) an LP kernel representing the seasonal repetition, 3) a rational quadratic kernel for small irregularities, and 4) and another SE kernel for the noise term. Here, we adopt an SE+LP kernel to capture the main features of the data. The explicit expression of the adopted kernel is given in the Appendix.

5.2 Result Analysis

Conclusion:

The following experimental results confirm that GPCV-ADMM is able to achieve comparable (even better) hyper-parameter estimation performance compared to GPML with much reduced computation complexity.

Estimation Performance:

Table 1 gives the quantitative comparisons between GPCV-ADMM and GPML across all synthetic data sets, where the underlying true hyper-parameter values are given as references. The results show that GPCV-ADMM hyper-parameter estimates are fairly close to the GPML estimates, and moreover they are both close to the true values. It is not surprising that GPCV-ADMM tends to generate lower test MSE ($MSE_{test}=\frac{1}{N}\left|\left|\boldsymbol{y}_{*}-\bar{\boldsymbol{m}}(\boldsymbol{X}_{*};\boldsymbol{\theta}^{*})\ \right|\right|_{2}^{2}$), simply because GPCV-ADMM is designed intentionally to minimize MSE and thus better reveals the predictive performance compared with the ML based scheme. According to the Monte Carlo simulation results, GPCV-ADMM estimation is more robust, demonstrating smaller sample standard deviation of both the hyper-parameter estimates and the test MSE.

In the above tests with synthetic data sets, we implicitly assumed that the underlying kernel function was precisely known, i.e., there was no model mismatch. In the following test performed on the real atmospheric $CO_{2}$ data, recorded from 1958 to 2015, GPCV-ADMM demonstrated outstanding training and prediction performance as shown in Figure 4, despite of model mismatch. We use the $CO_{2}$ data ranging from 1958 to 2008 for training. The hyper-parameter estimates obtained by the GPML in Rasmussen and Williams (2006) with several random restarts are $l_{1}=67$ years, $l_{2}=90$ years, and $l_{3}=1.3$ respectively. Though the adopted kernel function in our experiment is only a portion of the original one in Rasmussen and Williams (2006), it could be viewed as a more general function approximation with the ignorance of some small irregularities and noises, thus these reference values are still considered to be good enough. For GPCV-ADMM, we tried a few random restarts and picked the one with the lowest training MSE as suggested in Rasmussen and Williams (2006). We obtained $l_{1}=27$ years, $l_{2}=51$ years, and $l_{3}=1.26$. To evaluate the predictive performance of the two competing schemes, we predicted the real $CO_{2}$ concentrations from 2009 to 2015. The standardized test MSE obtained by GPML is 1.408, while our GPCV-ADMM gives a better prediction with a lower test MSE of 1.307.

Convergence Property:

GPCV-ADMM has its merit with rapid convergence to a moderately good estimation within a few iterations. The convergence curves of the selected parameters as well as the ADMM objective value are shown in Figure 2 for one particular Monte Carlo trial. It is noticed that the number of iterations needed to converge is not influenced as the data size increases. In our experiments, we manually set the maximum number of iterations to be 100 such that both the CT and the risk of over-fitting could be well reduced.

Computational Complexity:

We aim to verify that GPCV-ADMM has only $O(n^{2})$ complexity while GPML has $O(n^{3})$ complexity. Since the total number of iterations needed is not determined by increased sample size for both schemes, we only need to verify the complexity for each iteration. To this end, we fixed the kernel configuration and varied the data size from $n=500$ to 1000 and 2000, to see whether a quadratic increase in the computational time (CT) would be witnessed. Here, we treat the data size $n=500$ as the baseline, and compute the scaling factor (defined as $\frac{CT_{A}/CT_{B}}{n_{A}/n_{B}}$, where the subscript $B$ stands for the baseline data sets with $n=500$ and $A$ stands for data sets with larger sample sizes $n_{A}=1000,2000$).

The scaling factor across data sets generated from three kernels are shown in Figure 3. It is clear that the scaling factor of GPML is consistently larger than that of GPCV-ADMM. It is also noticed that a quadratic increase in the CT of GPCV-ADMM and a cubic increase in the CT of GPML would become more apparent as data size increases.

5.3 Implementation Details

The practical implementation of GPCV-ADMM requires special attentions to the following aspects.

Initialization:

A good starting point for both the hyper-parameters $\boldsymbol{\theta}$ and the auxiliary variable $\boldsymbol{z}_{T}$, will lead to faster and smoother convergence of GPCV-ADMM as observed in Figure 2. Random restarts could be adopted to alleviate the adverse impact of bad initializations.

Numerical Search:

We follow (12) to update the GP hyper-parameters numerically. Coordinate descent Bertsekas (2016) is adopted when $\boldsymbol{\theta}$ has more than one element. New GD type of methods such as the ADAM algorithm Kingma and Ba (2015) and other variants could be used for faster and more stable numerical search.

Choice of the regularization parameter $\rho$:

The magnitude of $\rho$ controls both the descent speed and the convexity of the ADMM objective function. A large $\rho$ endows a strong convexity of the ADMM objective function, yet often requiring more iterations to converge. A smaller $\rho$ endows faster descent, but the training procedure may get stuck at a bad local minimum more easily. When a suitable $\rho$ value is difficult to determine, one possible remedy, as suggested in Hong et al. (2016), is to use a different and smaller $\rho^{\prime}$ in (11c) for updating the dual variable.

6 Conclusion

In this paper, we proposed two general CV based GP hyper-parameter optimization schemes suitable for big data applications. By introducing a nonlinear equality constraint to avoid large-scale matrix inverse, the resulting GPCV-ADMM scheme was proven to reduce the $O(n^{3})$ computational complexity of the state-of-the-art GPML scheme considerably to $O(n^{2})$. Unlike the existing low-complexity GP methods, GPCV-ADMM does not make any sophisticated approximations, and it reduces the gap between the training- and test performance, and most favorably, is extremely easy to implement. Especially, the $K$-fold CV based scheme has the potential to exploit the multi-core processing in modern computing platforms, and is robust to over-fitting. Various experimental results validated the performance of the proposed scheme, which in some cases outperforms its GPML counterpart with smaller sample standard deviation of the hyper-parameter estimates, lower test MSE, and most importantly, a significant reduction in the computational complexity.

7 Appendix (supplementary)

7.1 Partial Derivatives and Gradient

The expression of $\boldsymbol{g}^{(\eta)}\left(\boldsymbol{\theta}\right)$ used for updating the GP hyper-parameters, $\boldsymbol{\theta}$, in (12) is obtained as:

[TABLE]

For each element of $\boldsymbol{\theta}$ (denoted as $\theta_{i}$), its partial derivative is computed as:

[TABLE]

7.2 Explicit Form of Kernel Functions

The expressions for the selected kernels that we use for the synthetic data are listed below.

• •

Squared Exponential (SE) Kernel

SE kernel is usually regarded as the default kernel for GP models, due to its great universality as well as many good properties. The length scale $l$ in an SE kernel specifies the width of the kernel and thereby determines the smoothness of the regression function.

[TABLE]

• •

**Locally Periodic (LP) Kernel

**Periodicity is another important pattern that people always get interested, especially in modeling time series data. Most of the real data do not repeat themselves exactly. Therefore combining a local kernel together with a periodic kernel into a locally periodic kernel, is considered to allow the shape of the repeating patterns to evolve over time:

[TABLE]

• •

**Composite SE + LP Kernel

**One good thing about using kernel function is its flexibility in combining various kernel components, which allows multiplications and/or additions over different kernels to capture different features of the data. In our experiments, we added up one SE kernel and one LP kernel to model local periodicity with trend.

[TABLE]