Unbiased Stochastic Optimization for Gaussian Processes on Finite Dimensional RKHS

TL;DR

This paper introduces exact stochastic inference algorithms for Gaussian Processes with kernels inducing finite or infinite dimensional RKHSs, improving convergence guarantees and performance under memory constraints.

Contribution

The authors develop novel algorithms for unbiased stochastic inference in GPs with finite-dimensional RKHS kernels, extending to infinite dimensions with approximate methods.

Findings

Outperforms existing methods in memory-limited scenarios

Achieves better convergence to true marginal likelihood

Effective for both finite and infinite dimensional RKHSs

Abstract

Current methods for stochastic hyperparameter learning in Gaussian Processes (GPs) rely on approximations, such as computing biased stochastic gradients or using inducing points in stochastic variational inference. However, when using such methods we are not guaranteed to converge to a stationary point of the true marginal likelihood. In this work, we propose algorithms for exact stochastic inference of GPs with kernels that induce a Reproducing Kernel Hilbert Space (RKHS) of moderate finite dimension. Our approach can also be extended to infinite dimensional RKHSs at the cost of forgoing exactness. Both for finite and infinite dimensional RKHSs, our method achieves better experimental results than existing methods when memory resources limit the feasible batch size and the possible number of inducing points.

Tables

Table 1. Table 1: Complexity analysis for one optimization iteration: C {\rm C} and M {\rm M} are the computational and the memory complexity in the computation of ∂ ϕ α ​ ( x ) ∂ α \frac{\partial\phi_{\alpha}\left(x\right)}{\partial\alpha} . b b is the batch size.

Algorithm	Computations	Storage
SVGP	$O\left(b^3\right)+O\left(b^{2}d\right)+bC$	$O\left(b^2\right)+O\left(bd\right)+bM$
BSGD	$O\left(b^3\right)+O\left(b^{2}d\right)+bC$	$O\left(b^2\right)+O\left(bd\right)+bM$
Minimax	$O\left(d^3\right)+O\left(b{d}^2\right)+bC$	$O\left(d^2\right)+O\left(bd\right)+bM$
CSGD	$O\left(d^3\right)+O\left(b{d}^2\right)+bC$	$O\left(d^2\right)+O\left(bd\right)+bM$

Table 2. Table 2: Negative log marginal likelihood with neural network and Linear kernel.

	method	MINIMAX(Ours)	SCGD(Ours)	SVGP	BSGD
batch size	name
32	bike	-1.482$\pm$0.249	-1.454$\pm$0.264	-0.895$\pm$0.296	-1.167$\pm$0.701
	elevators	-1.588$\pm$0.287	-1.530$\pm$0.195	-0.150$\pm$0.436	-1.125$\pm$0.017
	keggdirected	-0.745$\pm$0.076	-0.740$\pm$0.093	-0.151$\pm$1.057	-0.688$\pm$0.018
	keggundirected	-0.737$\pm$0.004	-0.739$\pm$0.011	-0.655$\pm$0.007	-0.716$\pm$0.025
	kin40k	-1.713$\pm$0.063	-1.760$\pm$0.090	0.792$\pm$0.881	-0.684$\pm$0.028
	pol	2.128$\pm$0.097	2.016$\pm$0.126	1.643$\pm$0.108	2.689$\pm$0.175
	protein	0.392$\pm$0.060	0.393$\pm$0.053	0.848$\pm$0.218	0.640$\pm$0.016
	slice	0.436$\pm$0.836	0.500$\pm$0.855	-0.045$\pm$0.127	0.459$\pm$0.089
	tamielectric	0.178$\pm$0.001	0.179$\pm$0.001	0.179$\pm$0.001	0.179$\pm$0.001
64	bike	-1.378$\pm$0.275	-1.366$\pm$0.259	-0.644$\pm$0.297	-1.282$\pm$0.617
	elevators	-1.523$\pm$0.276	-1.427$\pm$0.173	-0.152$\pm$0.440	-1.139$\pm$0.019
	keggdirected	-0.780$\pm$0.059	-0.748$\pm$0.105	-0.126$\pm$1.042	-0.726$\pm$0.020
	keggundirected	-0.747$\pm$0.004	-0.747$\pm$0.005	-0.670$\pm$0.005	-0.710$\pm$0.015
	kin40k	-1.654$\pm$0.045	-1.682$\pm$0.110	0.298$\pm$0.970	-0.952$\pm$0.059
	pol	2.206$\pm$0.299	2.171$\pm$0.250	1.610$\pm$0.128	2.413$\pm$0.069
	protein	0.325$\pm$0.013	0.381$\pm$0.060	0.659$\pm$0.104	0.618$\pm$0.021
	slice	0.471$\pm$1.059	0.086$\pm$0.110	-0.162$\pm$0.164	0.304$\pm$0.143
	tamielectric	0.178$\pm$0.001	0.178$\pm$0.001	0.179$\pm$0.001	0.179$\pm$0.001
128	bike	-1.340$\pm$0.039	-1.358$\pm$0.029	-0.832$\pm$0.231	-1.311$\pm$0.126
	elevators	-1.678$\pm$0.018	-1.399$\pm$0.022	0.044$\pm$0.003	-1.173$\pm$0.005
	keggdirected	-0.805$\pm$0.042	-0.778$\pm$0.091	-0.189$\pm$1.078	-0.763$\pm$0.009
	keggundirected	-0.725$\pm$0.022	-0.750$\pm$0.007	-0.679$\pm$0.007	-0.735$\pm$0.013
	kin40k	-1.590$\pm$0.024	-1.618$\pm$0.030	0.003$\pm$0.536	-1.224$\pm$0.038
	pol	2.215$\pm$0.156	2.150$\pm$0.140	1.583$\pm$0.014	2.223$\pm$0.037
	protein	0.289$\pm$0.008	0.327$\pm$0.010	0.709$\pm$0.069	0.576$\pm$0.003
	slice	0.019$\pm$0.043	0.044$\pm$0.068	-0.197$\pm$0.024	-0.038$\pm$0.119
	tamielectric	0.178$\pm$0.001	0.178$\pm$0.001	0.179$\pm$0.001	0.179$\pm$0.001

Table 3. Table 3: Negative log marginal likelihood with neural network and Gaussian kernel.

	method	MINIMAX(Ours)	SCGD(Ours)	SVGP	BSGD
batch size	name
16	bike	-1.013$\pm$0.082	-1.005$\pm$0.040	-1.038$\pm$0.029	0.105$\pm$0.107
	elevators	-1.094$\pm$0.018	-1.135$\pm$0.009	-0.915$\pm$0.012	-1.081$\pm$0.014
	keggdirected	-0.755$\pm$0.004	-0.754$\pm$0.007	-0.640$\pm$0.012	-0.566$\pm$0.012
	keggundirected	-0.714$\pm$0.000	-0.714$\pm$0.000	-0.641$\pm$0.000	-0.615$\pm$0.000
	kin40k	-1.559$\pm$0.015	-1.606$\pm$0.031	-0.697$\pm$0.027	-0.515$\pm$0.042
	pol	2.066$\pm$0.029	1.980$\pm$0.027	1.892$\pm$0.035	2.656$\pm$0.192
	protein	0.485$\pm$0.007	0.476$\pm$0.007	0.625$\pm$0.016	0.716$\pm$0.010
	slice	1.548$\pm$0.058	1.556$\pm$0.042	0.369$\pm$0.062	0.461$\pm$0.143
	tamielectric	0.179$\pm$0.001	0.179$\pm$0.001	0.179$\pm$0.001	0.179$\pm$0.001
32	bike	-0.929$\pm$0.023	-0.968$\pm$0.065	-1.098$\pm$0.030	0.083$\pm$0.055
	elevators	-1.042$\pm$0.023	-1.111$\pm$0.016	-0.926$\pm$0.007	-1.092$\pm$0.011
	keggdirected	-0.764$\pm$0.007	-0.767$\pm$0.007	-0.661$\pm$0.007	-0.650$\pm$0.012
	keggundirected	-0.735$\pm$0.000	-0.737$\pm$0.000	-0.656$\pm$0.000	-0.714$\pm$0.000
	kin40k	-1.585$\pm$0.018	-1.571$\pm$0.017	-0.796$\pm$0.031	-0.695$\pm$0.021
	pol	2.115$\pm$0.027	2.112$\pm$0.029	1.807$\pm$0.023	2.521$\pm$0.192
	protein	0.479$\pm$0.008	0.479$\pm$0.004	0.595$\pm$0.006	0.674$\pm$0.009
	slice	1.950$\pm$0.120	1.815$\pm$0.142	0.159$\pm$0.042	0.248$\pm$0.098
	tamielectric	0.179$\pm$0.001	0.179$\pm$0.001	0.179$\pm$0.001	0.179$\pm$0.001
64	bike	-0.912$\pm$0.035	-0.946$\pm$0.049	-1.153$\pm$0.011	-0.114$\pm$0.060
	elevators	-0.988$\pm$0.013	-1.056$\pm$0.018	-0.932$\pm$0.004	-1.111$\pm$0.007
	keggdirected	-0.762$\pm$0.006	-0.771$\pm$0.009	-0.672$\pm$0.005	-0.694$\pm$0.008
	keggundirected	-0.739$\pm$0.000	-0.743$\pm$0.000	-0.666$\pm$0.000	-0.691$\pm$0.000
	kin40k	-1.547$\pm$0.022	-1.539$\pm$0.039	-0.864$\pm$0.026	-0.947$\pm$0.040
	pol	2.433$\pm$0.081	2.399$\pm$0.059	1.754$\pm$0.023	2.237$\pm$0.030
	protein	0.479$\pm$0.007	0.474$\pm$0.007	0.585$\pm$0.011	0.644$\pm$0.011
	slice	1.964$\pm$0.846	1.899$\pm$0.764	0.092$\pm$0.036	-0.018$\pm$0.059
	tamielectric	0.179$\pm$0.001	0.179$\pm$0.001	0.179$\pm$0.001	0.179$\pm$0.001

Table 4. Table 4: RMSE with neural network and Linear kernel.

	method	MINIMAX(Ours)	SCGD(Ours)	SVGP	BSGD
batch size	name
32	bike	0.099$\pm$0.011	0.095$\pm$0.007	0.093$\pm$0.055	0.082$\pm$0.073
	elevators	0.111$\pm$0.008	0.158$\pm$0.114	0.217$\pm$0.071	0.090$\pm$0.003
	keggdirected	0.120$\pm$0.011	0.123$\pm$0.015	0.376$\pm$0.557	0.122$\pm$0.009
	keggundirected	0.120$\pm$0.006	0.121$\pm$0.008	0.119$\pm$0.003	0.119$\pm$0.003
	kin40k	0.061$\pm$0.012	0.056$\pm$0.011	0.651$\pm$0.310	0.126$\pm$0.010
	pol	2.280$\pm$0.130	2.328$\pm$0.123	2.593$\pm$0.174	3.402$\pm$0.152
	protein	0.426$\pm$0.024	0.428$\pm$0.023	0.550$\pm$0.079	0.483$\pm$0.006
	slice	0.642$\pm$0.327	0.697$\pm$0.282	0.591$\pm$0.119	0.992$\pm$0.339
	tamielectric	0.290$\pm$0.002	0.289$\pm$0.002	0.289$\pm$0.002	0.289$\pm$0.002
64	bike	0.112$\pm$0.012	0.108$\pm$0.012	0.127$\pm$0.054	0.078$\pm$0.064
	elevators	0.111$\pm$0.013	0.114$\pm$0.023	0.216$\pm$0.071	0.094$\pm$0.010
	keggdirected	0.118$\pm$0.008	0.120$\pm$0.011	0.377$\pm$0.557	0.119$\pm$0.009
	keggundirected	0.119$\pm$0.006	0.121$\pm$0.008	0.121$\pm$0.008	0.123$\pm$0.006
	kin40k	0.065$\pm$0.015	0.060$\pm$0.010	0.443$\pm$0.338	0.097$\pm$0.009
	pol	2.449$\pm$0.160	2.263$\pm$0.151	2.415$\pm$0.115	3.067$\pm$0.145
	protein	0.430$\pm$0.007	0.426$\pm$0.021	0.494$\pm$0.046	0.473$\pm$0.003
	slice	0.479$\pm$0.066	0.615$\pm$0.288	0.554$\pm$0.160	0.680$\pm$0.173
	tamielectric	0.290$\pm$0.002	0.289$\pm$0.002	0.289$\pm$0.002	0.289$\pm$0.002
128	bike	0.125$\pm$0.005	0.123$\pm$0.011	0.111$\pm$0.024	0.056$\pm$0.004
	elevators	0.110$\pm$0.004	0.101$\pm$0.002	0.250$\pm$0.007	0.091$\pm$0.003
	keggdirected	0.116$\pm$0.008	0.118$\pm$0.011	0.370$\pm$0.561	0.123$\pm$0.018
	keggundirected	0.124$\pm$0.008	0.121$\pm$0.007	0.118$\pm$0.006	0.119$\pm$0.007
	kin40k	0.073$\pm$0.002	0.071$\pm$0.003	0.283$\pm$0.193	0.078$\pm$0.003
	pol	2.469$\pm$0.174	2.413$\pm$0.234	2.537$\pm$0.094	2.642$\pm$0.110
	protein	0.435$\pm$0.009	0.436$\pm$0.009	0.513$\pm$0.023	0.461$\pm$0.006
	slice	0.528$\pm$0.086	0.546$\pm$0.181	0.542$\pm$0.111	0.583$\pm$0.086
	tamielectric	0.290$\pm$0.002	0.289$\pm$0.002	0.289$\pm$0.002	0.289$\pm$0.002

Table 5. Table 5: RMSE with neural network and Gaussian kernel.

	method	MINIMAX(Ours)	SCGD(Ours)	SVGP	BSGD
batch size	name
16	bike	0.081$\pm$0.012	0.096$\pm$0.022	0.046$\pm$0.007	0.220$\pm$0.014
	elevators	0.096$\pm$0.003	0.098$\pm$0.004	0.089$\pm$0.002	0.089$\pm$0.003
	keggdirected	0.118$\pm$0.008	0.118$\pm$0.007	0.122$\pm$0.005	0.127$\pm$0.005
	keggundirected	0.131$\pm$0.000	0.128$\pm$0.000	0.116$\pm$0.000	0.138$\pm$0.000
	kin40k	0.039$\pm$0.002	0.040$\pm$0.002	0.110$\pm$0.007	0.147$\pm$0.007
	pol	2.212$\pm$0.146	2.200$\pm$0.104	2.674$\pm$0.081	3.774$\pm$0.440
	protein	0.387$\pm$0.009	0.388$\pm$0.009	0.471$\pm$0.008	0.507$\pm$0.012
	slice	0.519$\pm$0.149	0.496$\pm$0.110	0.562$\pm$0.067	0.661$\pm$0.090
	tamielectric	0.289$\pm$0.002	0.289$\pm$0.002	0.289$\pm$0.002	0.289$\pm$0.002
32	bike	0.093$\pm$0.010	0.095$\pm$0.012	0.047$\pm$0.010	0.235$\pm$0.024
	elevators	0.100$\pm$0.004	0.098$\pm$0.004	0.088$\pm$0.003	0.089$\pm$0.002
	keggdirected	0.118$\pm$0.008	0.118$\pm$0.007	0.121$\pm$0.007	0.125$\pm$0.008
	keggundirected	0.130$\pm$0.000	0.128$\pm$0.000	0.115$\pm$0.000	0.116$\pm$0.000
	kin40k	0.039$\pm$0.002	0.039$\pm$0.001	0.096$\pm$0.005	0.124$\pm$0.003
	pol	2.210$\pm$0.133	2.209$\pm$0.080	2.552$\pm$0.117	3.338$\pm$0.151
	protein	0.387$\pm$0.008	0.389$\pm$0.010	0.461$\pm$0.009	0.493$\pm$0.007
	slice	0.511$\pm$0.039	0.515$\pm$0.065	0.527$\pm$0.139	0.565$\pm$0.102
	tamielectric	0.289$\pm$0.002	0.289$\pm$0.002	0.289$\pm$0.002	0.289$\pm$0.001
64	bike	0.101$\pm$0.008	0.106$\pm$0.014	0.035$\pm$0.001	0.209$\pm$0.009
	elevators	0.089$\pm$0.003	0.098$\pm$0.007	0.088$\pm$0.003	0.088$\pm$0.003
	keggdirected	0.118$\pm$0.008	0.117$\pm$0.008	0.119$\pm$0.006	0.122$\pm$0.007
	keggundirected	0.128$\pm$0.000	0.127$\pm$0.000	0.113$\pm$0.000	0.117$\pm$0.000
	kin40k	0.041$\pm$0.001	0.041$\pm$0.002	0.088$\pm$0.005	0.096$\pm$0.003
	pol	2.348$\pm$0.077	2.353$\pm$0.088	2.512$\pm$0.081	3.020$\pm$0.071
	protein	0.388$\pm$0.008	0.389$\pm$0.005	0.456$\pm$0.005	0.481$\pm$0.004
	slice	0.616$\pm$0.163	0.546$\pm$0.129	0.507$\pm$0.116	0.512$\pm$0.120
	tamielectric	0.289$\pm$0.002	0.289$\pm$0.002	0.289$\pm$0.002	0.289$\pm$0.002

Equations

$$⟨A,B⟩=\mathbf{Tr}\left(A{B}^T\right).$$

$$\Vert A\Vert =\sqrt{⟨A,A⟩}.$$

$$f_{\alpha }\left(U\right)\sim \mathcal{N}\left(0,k_{\alpha }\left(U,U\right)\right).$$

$$p\left(\mathbf{y}|X,\alpha ,\sigma \right)=\mathcal{N}\left(\mathbf{y}|0,K\left(\alpha \right)+{\sigma }^{2}I\right)$$

$${\mathbf{y}}^T{\left(K\left(\alpha \right)+{\sigma }^{2}I\right)}^{-1}\mathbf{y}+\log \left|K\left(\alpha \right)+{\sigma }^{2}I\right|.$$

$$k_{\alpha }\left(x,x^{\prime }\right)={\varphi }_{\alpha }{\left(x\right)}^T{\varphi }_{\alpha }\left(x^{\prime }\right)$$

$${\mathbf{b}}^T{\left(V{V}^T+\lambda I\right)}^{-1}\mathbf{b}=\underset{\mathbf{w}}{\min }\frac{1}{\lambda }\Vert V\mathbf{w}-\mathbf{b}\Vert +{\Vert \mathbf{w}\Vert }^2$$

$$l\left(\theta \right)=\frac{1}{{\sigma }^2}{\Vert Z\left(\alpha \right)\mathbf{w}-\mathbf{y}\Vert }^2+{\Vert \mathbf{w}\Vert }^2+\log \left|F\left(\theta \right)\right|+\left(n-d\right)\log {\sigma }^2,$$

$$\theta$$

$$Z\left(\alpha \right)$$

$$F\left(\theta \right)$$

$${\nabla }_{\theta }\log \left|F\left(\theta \right)\right|=\sum _{i=1}^{n}G\left(\theta ;x_i\right)$$

$$g_i\left(\theta \right)$$

$$g\left(\theta \right)$$

$$F_i\left(\theta \right)$$

$$h\left(A\right)$$

$$F\left(\theta \right)$$

$$l\left(\theta \right)$$

$$\underset{\theta }{\min }g\left(\theta \right)+h\left(F\left(\theta \right)\right)$$

$$\underset{\theta ,A}{\min }g\left(\theta \right)+h\left(A\right)$$

$$\text{s.t. }A=F\left(\theta \right)$$

$$\underset{\zeta }{\min }{l}_{\mu }\left(\zeta \right)$$

$$l_{\mu }\left(\zeta \right)=g\left(\theta \right)+h\left(A\right)+\mu \frac{\Vert A-F\left(\theta \right)\Vert }{\Vert A\Vert }.$$

$$\frac{\Vert A-F\left(\theta \right)\Vert }{\Vert A\Vert }$$

$$\underset{\zeta }{\min }\underset{\Vert B\Vert \le 1}{\max }\Psi \left(\zeta ,B\right),$$

$$\Psi \left(\zeta ,B\right)=g\left(\theta \right)+h\left(A\right)+\mu \frac{⟨B,A-F\left(\theta \right)⟩}{\Vert A\Vert }$$

$$\Psi \left(\zeta ,B\right)=\sum _{i=1}^n\psi \left(\zeta ,B;x_i\right)$$

$$\psi \left(\zeta ,B;x_i,\mu \right)=g_i\left(\theta \right)+\frac{1}{n}h\left(A\right)+\mu \frac{⟨B,\frac{1}{n}A-F_i\left(\theta \right)⟩}{\Vert A\Vert }$$

$$1.$$

$$2.$$

$$1.$$

$$2.$$

$$u\left(\theta \right)$$

$$v\left(u_1,u_2\right)$$

$$\tilde{u}\left(\theta ;\omega \right)$$

$$\frac{\partial \log \left|A\right|}{\partial A}=A^{-1}$$

Full text

1 Introduction

Gaussian Processes (GPs) provide a powerful probabilistic framework which has been applied to a wide range of learning applications, such as multi-task learning (Alvarez et al., 2012; Liu et al., 2018a), active learning (Liu et al., 2018b), semi-supervised learning Jean et al. (2018), and reinforcement learning (Srinivas et al., 2010; Shahriari et al., 2015). These successes can be attributed to the natural way in which uncertainty is incorporated into the predictions via a Bayesian interpretation. However, hyperparameter learning scales poorly. For a general covariance function (kernel function in kernel method’s parlance) the computational cost grows as $O(n^{3})$ where $n$ denotes the number of samples, and the storage resources grow as $O(n^{2})$. As a consequence, approximations are required for any modern application that involves GPs and big data.

It is therefore unsurprising that much research effort has been devoted to approximate methods for learning GPs. For a comprehensive overview of this research, we refer the reader to the work of Liu et al. (2020). The most expensive operation in training GPs is the inversion of the $n\times n$ kernel matrix which is required during the maximization of the marginal likelihood and for the computation of the posterior of the responses. As a result, much research has been devoted to providing a cheap approximation of the kernel matrix.

Broadly speaking approximations of GPs can be split in to data dependent and data independent methods. Data dependent approximations usually come with valuable probabilistic interpretations. An important work in this line of research is the seminal work of Quinonero-Candela & Rasmussen (2005), which provides a unified view on previous work by using the concept of inducing points. Since then quite a few followup works have used inducing points, while exploiting the tool of variational inference as a theoretical platform (Titsias, 2009; Nguyen et al., 2014; Wilson & Nickisch, 2015; Zhao & Sun, 2016). This approach is closely related to the Nyström approximation of kernel matrices (Zhao & Sun, 2016).

In contrast, data independent methods rely on approximating the kernel function itself. Typically, it is approximated by an inner product between low dimensional vectors (Rahimi & Recht, 2007; Yang et al., 2015; Shustin & Avron, 2021).

Approximating the inverse of the kernel matrix is enough for frequentist kernel methods such as Kernel Ridge Regression. However, in order to harness the full power of Bayesian kernel methods such as Gaussian Process Regression (GPR), we need to be also able to maximize the marginal likelihood, and for that an approximation of the inverse of the kernel matrix is not enough, due to an additional log-determinant term. This is especially true in cases where the covariance function depends on a large number of parameters, e.g., evaluating the covariance function involves a forward pass of a deep neural network (Wilson et al., 2016b; Calandra et al., 2016; Wilson et al., 2016a). In such cases, it is important to be able to use a stochastic mini-batches based optimization, since otherwise the cost of making a pass over the entire dataset typically proves too expensive. Moreover, it is well appreciated in the literature that when the model involves deep neural networks, e.g., when the covariance function is defined by a neural network, it is important to use stochastic gradients, since they enable more efficient optimization and better generalization (Goodfellow et al., 2016).

Currently, two main approaches are used for stochastic hyperparameter learning in GPs. The first relies on the inducing-points framework and employs stochastic variational inference (SVI) (Hoffman et al., 2013; Hensman et al., 2013; Hoang et al., 2015; Wilson et al., 2016a). The second approach is more direct: it is based on computing stochastic batches while ignoring the fact that they provide only biased estimates of gradients. Interestingly, a recent work shows that despite the bias, given a large enough batch size, the direct approach produces almost optimal models (in terms of the marginal likelihood) (Chen et al., 2020).

Both the SVI approach and the direct approach suffer from several disadvantages. For example, consider the case in which the covariance function is the inner product between two feature maps of a moderate dimension. This can be either because this decomposition serves as an approximation to another covariance function with an infinite dimensional * reproducing kernel Hilbert space* (RKHS), or because it has been defined this way to begin with, e.g., as an inner product between features that are created by passing the data through a neural network. Either way, using the SVI approach amounts to imposing an additional unnecessary approximation that comes from the need to use inducing points. As for the second approach, its applicability is highly dependent on our ability to use large batches. This can be a serious impediment in several cases, such as optimization on weak edge devices (Chen et al., 2016).

In this paper, we propose two stochastic optimization algorithms based on mini-batches for maximizing the marginal likelihood of GPs (i.e., learning the hyperparameters). The first algorithm is based on reframing the problem as a nonconvex-concave minimax problem. We then leverage recent advancements in the theory of solving such problems (Boţ & Böhm, 2020; Lin et al., 2020; Luo et al., 2020) to propose a concrete algorithm. The second algorithm is based on writing the loss function in compositional form, i.e., as the composition of a function and the expected value of another function. We then use the recently introduced Stochastic Compositional Gradient Descent method (Wang et al., 2017). Our novel algorithms are applicable for covariance functions connected with an RKHS of a moderate dimension, and guarantee convergence of the marginal likelihood to a local minimum for any batch size without need for further approximations. In the infinite dimensional case (e.g., Gaussian covariance function) one can use our method on top of a low rank approximation of the covariance function, e.g., using the random features method Rahimi & Recht (2007). Our experiments show that not only is our method superior to existing methods for stochastic optimization of the marginal likelihood in the finite dimensional case when the batches have a moderated size, it is also superior to the existing methods in the infinite dimensional case if the restriction on the batch size is more severe, although in the infinite dimensional case the models found by our method are no longer optimal.

Additional Related Work

Apart from the extensive literature on scaling GPs using approximations, several works are focused on exact inference using sophisticated distributed algorithms (Nguyen et al., 2019; Wang et al., 2019). In this context, the biased stochastic gradient proposed by Chen et al. (2020) can be considered as an economical method for exact inference, given that the covariance function and the system enable computation in large enough batches.

2 Preliminaries

2.1 Notations and Basic Definitions

For a function $f$:${\cal U}\to\mathbb{R}$ and $U=\left(u_{1},\dots,u_{m}\right)\in{\cal U}^{m}$, ${\mathbf{f}}=f\left(U\right)$ is a vector in $\mathbb{R}^{m}$ such that ${\mathbf{f}}_{i}=f\left(u_{i}\right)$. Similarly, for a binary function $k:{\cal U}\times{\cal U}\to\mathbb{R}$, $K=k\left(U,U\right)$ is a matrix in $\mathbb{R}^{m\times m}$ such that $K_{ij}=k\left(u_{i},u_{j}\right)$. Given a size $m$ vector ${\mathbf{b}}$ and ${\cal S}=\left(s_{1},\dots,s_{p}\right)\in\left\{1,\dots,m\right\}^{p}$, we use ${\mathbf{b}}_{{\cal S}}$ to denote the size $p$ vector such that the $i$’th coordinate of ${\mathbf{b}}_{{\cal S}}$ is equal to the $s_{i}$’th coordinate of ${\mathbf{b}}$. In a similar way if $C$ is an $m\times n$ matrix then $C_{{\cal S}}$ is a $p\times n$ matrix such that the $i$’th row of $C_{S}$ is equal to the $s_{i}$’th row of $C$. Finally, if ${\cal R}=\left(r_{1},\dots,r_{q}\right)\in\left\{1,\dots,n\right\}^{q}$ then $C_{{\cal SR}}$ is a $p\times q$ matrix such that the $\left(i,j\right)$ coordinate of $C_{{\cal SR}}$ is equal to the $\left(s_{i},r_{j}\right)$ coordinate of $C$.

For a square matrix $A$ we use the notation $\left|A\right|$ to denote the determinant of $A$. If ${\cal S}$ is a finite sequence or a set then $\left|{\cal S}\right|$ denotes its length or size. Whenever we use $\left\langle A,B\right\rangle$ where $A$ and $B$ are real matrices, $\left\langle\cdot,\cdot\right\rangle$ symbolizes the Frobenius inner product which is defined as

[TABLE]

If $A$ is a real matrix then $\left\|A\right\|$ is the Frobenius norm of $A$ so,

[TABLE]

For a vector ${\mathbf{v}}$, $\left\|{\mathbf{v}}\right\|$ is the Euclidean norm of ${\mathbf{v}}\in\mathbb{R}^{q}$.

For any closed convex set $\Omega\subseteq\mathbb{R}^{q}$ and ${\mathbf{v}}\in\mathbb{R}^{q}$, $\operatorname{proj}_{\Omega}\left({\mathbf{v}}\right)=\arg\min_{{\mathbf{v}}^{\prime}\in\Omega}\left\|v^{\prime}-v\right\|$ is the Euclidean projection of ${\mathbf{v}}$ on $\Omega$.

2.2 Hyperparameter Learning in Gaussian Process Regression

Let ${\cal X}$ be a feature space and let $k_{\alpha}:{\cal X}\times{\cal X}\to\mathbb{R}$ be a positive definite covariance function parameterized by hyperparameters $\alpha\in\mathbb{R}^{m}$. For any $\alpha\in\mathbb{R}^{m}$, let $f_{\alpha}$ be a random function on ${\cal X}$ distributed as a zero mean GP whose covariance is $k_{\alpha}$, that is, for any $j\in\mathbb{N}$, $U\in{\cal X}^{j}$:

[TABLE]

Let $X\in{\cal X}^{n}$, ${\mathbf{y}}\in\mathbb{R}^{n}$ be a training set, where $n$ is the number of training samples. In Gaussian Process Regression (GPR), it is assumed that ${\mathbf{y}}$ is a sample of the random vector $f_{\bar{\alpha}}\left(X\right)+\bar{\sigma}^{2}\mathbf{\epsilon}$ where $\mathbf{\epsilon}\sim{\cal\mathcal{N}}\left(0,I_{n}\right)$, and $\bar{\alpha},\bar{\sigma}^{2}$ are the true hyperparameters of the model.

In GPR, hyperparameters are learned by solving a maximum likelihood type II problem, i.e., maximizing the marginal likelihood.

[TABLE]

with respect to $\alpha$ and $\sigma$, where $K\left(\alpha\right)=k_{\alpha}\left(X,X\right)$. The maximization of the marginal likelihood $p\left({\mathbf{y}}|X,\alpha,\sigma\right)$ is equivalent to the minimization of

[TABLE]

See Rasmussen (2003) for details. In this work we further assume that the covariance function has the following form:

[TABLE]

for some feature map $\phi_{\alpha}:{\cal X}\to\mathbb{R}^{d}$. It can be shown that for all $\lambda>0$, $V\in\mathbb{R}^{n\times d}$ and ${\mathbf{b}}\in\mathbb{R}^{d}$ we have that

[TABLE]

(see Appendix B). As a result, maximizing $p\left({\mathbf{y}}|X,\alpha,\sigma\right)$ is equivalent to minimizing

[TABLE]

where

[TABLE]

The goal of this work is to propose algorithms for minimizing $l\left(\theta\right)$ using a stochastic gradient method based on mini-batches. This would have been straightforward if we could write

[TABLE]

for some function $G(\cdot;\cdot)$. However, there is no obvious decomposition of this form.

3 Stochastic Optimization for Gaussian Processes

This section contains our main contribution: two novel stochastic mini-batched based algorithms for minimizing $l(\theta)$. Before detailing our approaches, we make few additional notations:

[TABLE]

So we can write,

[TABLE]

3.1 A Minimax Approach

Our first step is to replace the problem

[TABLE]

with the equivalent problem

[TABLE]

The next step is to replace the last problem with a parameterized problem, such that the hard constraint is replaced with a penalty term, and the penalty term is driven to infinity. To do so, let us first define the optimization problems

[TABLE]

where $\zeta=\left(\theta,A\right)$ and

[TABLE]

Now, suppose that Problem (3.1) admits an optimal solution. Let the sequences $\{\mu_{k}\}$ and $\{\zeta_{k}\}$ satisfy $\mu_{k}\to\infty$ and, for each $k$, let $\zeta_{k}$ minimize $l_{\mu_{k}}(\zeta)$. If $\zeta^{\star}$ is an accumulation point of $\{\zeta_{k}\}$, then $\zeta^{\star}$ is an optimal solution of Problem (3.1) (Ruszczynski, 2006, Theorem 6.6).

The reason we use

[TABLE]

instead of just $\left\|A-F\left(\theta\right)\right\|$ is to avoid finding critical points around $F\left(\theta\right)=0$.

Unfortunately, the term $\left\|A-F\left(\theta\right)\right\|$ exhibits the same issues as $h\left(F\left(\theta\right)\right)=\log\det F(\theta)$ and does not allow straightforward unbiased mini-batch based stochastic gradients.

In order to handle this, we use the fact that for an element $v$ in Euclidean space we have that $\left\|v\right\|=\max_{\left\|u\right\|\leq 1}\left\langle u,v\right\rangle$. Thus, minimizing $l_{\mu}\left(\zeta\right)$ is equivalent to

[TABLE]

where

[TABLE]

Now, we can write

[TABLE]

where

[TABLE]

Thus, Problem (3.2) naturally admits unbiased stochastic gradients based on mini-batches taken separately for the minimization and the maximization parts.

For better stability, we restrict $A$ to be positive semidefinite such that $A\succeq\sigma^{2}I$, and restrict $\sigma\geq\sigma_{\min}$, where $\sigma_{m}in$ is a non-learned hyperparameter. To further restrict the search space, we constrain $\zeta\preceq\zeta_{\max}$, where $\zeta_{\max}$ is another non-learned hyperparameter. Together, the constraints on $\zeta$ can be represented by $\zeta\in\Omega_{1}$, where $\Omega_{1}$ is a convex set.

Let $\Omega_{2}=\left\{B\in\mathbb{R}^{d\times d}\mid\left\|B\right\|\leq 1\right\}$. The algorithm we propose is based on the following update rule for solving Problem (3.2):

[TABLE]

where ${\cal S}_{t}$ and $\bar{{\cal S}}_{t}$ are sets of $s$ indices randomly chosen from $\left\{1,\dots,n\right\}$ independently from all previous iterations. Recent work by Boţ & Böhm (2020) shows that in case of a minimax problem where maximization is of a concave function over a convex constraint111The setup in Boţ & Böhm (2020) is more general., then with a few additional mild conditions on the objective function, an algorithm based on the above update rule converges to a stationary point in $O\left(\epsilon^{-8}\right)$ iterations.

3.2 Stochastic Compositional Gradient Descent Approach

Consider a loss function $l:\Theta\to\mathbb{R}$ of the form $l\left(\theta\right)=v\left(u\left(\theta\right)\right)$ where $u:\Theta\to\mathbb{R}^{p}$ and $v:\mathbb{R}^{p}\to\mathbb{R}$ are differentiable functions, and assume that $u\left(\theta\right)=\operatorname{\mathbb{E}}_{\omega}\tilde{u}\left(\theta;\omega\right)$ for a differentiable function $\tilde{u}\left(\theta;\omega\right)$ that depends on a random variable $\omega$.* Stochastic Compositional Gradient Descent* (SCGD) (Wang et al., 2017) is an intuitive algorithm that alternates between two steps: updating the solution $\theta_{t}$ by a stochastic gradient iteration, and estimating $u\left(\theta_{t}\right)$ using an iterative weighted average of past values. More precisely, the update rule of of SCGD is given by:

[TABLE]

where $\omega_{1},\omega_{2},\dots$ are samples from $\omega$ in an i.i.d. manner, and $a_{0},a_{1},\dots,b_{0},b_{1},\dots$ degrees of freedom in the algorithm. Under few additional standard conditions on $u,v,\tilde{u}$, Wang et al. (2017) showed a convergence rate of ${\cal O}\left(\epsilon^{-4}\right)$ if we choose $a_{t}=t^{-\frac{3}{4}}$ and $b_{t}=t^{-\frac{1}{2}}$, see (Wang et al., 2017, Theorem 8).

In order to use SCGD for our loss, we define:

[TABLE]

where $\omega={\cal S}$ is a random set of indices chosen from $\left\{1,\dots,n\right\}.$ With that, given the fact that

[TABLE]

we can write an explicit update rule for our SCGD-based algorithm:

[TABLE]

4 Comparison to Existing Methods

In this work we suggest two novel methods for stochastic optimization of the marginal likelihood. We recognize two main competing methods that also optimize the marginal likelihood stochastically. The simplest among them is what we henceforth refer to as Biased Stochastic Gradient Decent (BSGD) (Chen et al., 2020). The idea in BSGD is to take the gradient of the marginal likelihood using data only from the current batch, ignoring the fact that this produces only a biased estimate of the full gradient. On the other hand Scalable Variational Gaussian Process (SVGP) of Hensman et al. (2015) is a sophisticated approach that approximates the original inference problem using inducing points and stochastic variational inference.

Unlike our algorithms, how well each of the competing methods approximates the solution of the true problem depends on memory consumption. In the case of BSGD, the bias in stochastic gradients shrinks as batch size increases. On the other hand, inference quality SVGP crucially depends on the number of inducing points which need to be processed forward and backward by a neural network at each iteration.

4.1 Complexity Analysis

Let ${\rm C}$ and ${\rm M}$ be the computational and the memory complexity in the computation of $\frac{\partial\phi_{\alpha}\left(x\right)}{\partial\alpha}$, and let $b$ be the batch size. For SVGP, assume that the number of the inducing points is in the same order of $b$.

Number of inducing points vs. batch size.

For the sake if comparison, We choose the number of inducing points to be roughly the same as the mini-batch size because each inducing point requires a forward and backward pass through the network, storing activations and gradients just like a data sample. Thus, both sets consume memory in a similar way during each iteration, and matching their sizes keeps the per-iteration memory cost balanced and comparable.

The computational complexity of computing one optimization iteration of both BSGD and SVGP is $O\left(b^{3}\right)+O\left(b^{2}d\right)+b{\rm C}$ while both our algorithm take $O\left(d^{3}\right)+O\left(bd^{2}\right)+b{\rm C}$ per iteration. For storage, BSGD and SVGP take $O\left(b^{2}\right)+O\left(bd\right)+b{\rm M}$ while ours takes $O\left(d^{2}\right)+O\left(bd\right)+b{\rm M}$.

Typically ${\rm M}$ is relatively large as it is the storage used by AD to compute the gradient of a neural network. In small devices that means one might have to use small batches. Our approach aims not to harm the exactness of the algorithm in this case. Naturally, to use our algorithms effectively, the feature dimension $d=\dim\!\phi_{\alpha}(x)$ must remain moderate, since memory scales as $O(d^{2})$ and computation as $O(d^{3})$.

4.2 Relation to BSGD

Our approach is close to the BSGD approach of Chen et al. (2020). BSGD is based on the following update formula:

[TABLE]

where ${\cal S}_{t}$ are an independent random batch of indices and

[TABLE]

In our case, since $k_{\alpha}\left(x\right)=\phi_{\alpha}\left(x\right)^{T}\phi_{\alpha}\left(x\right)$, we can write $l\left(\theta;{\cal S}\right)$ as

[TABLE]

This leads to the update formula

[TABLE]

Writing the update formula of the BSGD algorithm this way emphasizes the fact that the only difference between the BSGD and the SCGD algorithms lays in the way in which $F\left(\theta_{t}\right)$ is approximated: using $\tilde{F}_{t}=\sum_{i\in{\cal S}}F_{i}\left(\theta_{t}\right)$ in BSGD, and $\tilde{F}_{t+1}=\left(1-b_{t}\right)\tilde{F}_{t}+b_{t}\frac{n}{\left|{\cal S}\right|}\sum_{i\in{\cal S}}F_{i}\left(\theta_{t}\right)$ in SCGD. Intuitively, the exponential smoothing that occurs SCGD should provide additional numerical robustness beyond the theoretical advantage of converging to a stationary point without an error which depends on the batch size.

We can also see that in BSGD the scale of $\tilde{F}_{t}$ is incorrect because it misses the multiplication of $\sum_{i\in{\cal S}}F_{i}\left(\theta_{t}\right)$ by $\frac{n}{\left|{\cal S}\right|}$.

5 Experimental Result

In our experiments we consider a covariance function of the form

[TABLE]

where $g_{w}$ is a neural network comprised of two fully connected layers, both with output dimension of 128 and ReLU activation function, and $k^{\prime}_{u}$ which is either the linear kernel

[TABLE]

or the Gaussian kernel parametrized with two hyper parameter, length scale $u_{1}$ and magnitude $u_{2}$, that is

[TABLE]

Setup:

The experiment is designed such that we can see the influence of the batch size on the result. We used nine regression datasets from the UCI repository (Asuncion & Newman, 2007), all with number of samples above 14,000 and less the 60,000. We tested the algorithms using different batch sizes: for the linear kernel we examined batch sizes of 32, 64, 128, 256 and 512, and for the Gaussian kernel we examined batch sizes of 16, 32, 64, 128 and 256. For SVGP, we set the number of inducing point equal to batch size. For our algorithms, labeled MINIMAX (Subsection 3.1) and SCGD (Subsection 3.2), we approximated the Gaussian kernel $k^{\prime}_{u}$ using the Random Fourier Features method with random features $\varphi$ of dimension 1000 such that

[TABLE]

We ran each test on five different splits of 90% train 10% test. We used AdaDelta for all methods and for each combination of dataset, split and method, and used grid search in order to select the learning rate that achieves minimal marginal likelihood. For SCGD, we fixed $b_{t}=0.9$. Note that this is very close to $1$ which means that SCGD becomes quite similar to BSGD and much of the improvement comes just from the correct scaling of $\tilde{F}_{t}$ . For MINIMAX we fixed $\mu_{t}=1.0$, since we found it enough for achieving good results despite that in theory it should be increased in an outer loop. We ran each algorithm for 100 epochs and used the hyperparameters from the iteration in which the marginal likelihood achieved its minimal value.

Result with linear kernel - the exact case:

It seems that the fact that our algorithm does an exact stochastic optimization brings a significant improvement over existing methods in the optimization of the marginal likelihood. As expected, we can see that this advantage is more significant when the batch size is smaller (see Table 2). However, lower negative log marginal likelihood does not always translate to lower MSE on the test set as we can see in Table 4. We saw that sometimes by using a suboptimal learning rate that does not achieve the minimal loss our algorithms can achieve better results in terms of test RMSE. However, since our work is focused on optimization, we do not use procedures such as early stopping or cross validation which could potentially improve the result from the test RMSE perspective.

Result with Gaussian kernel:

Since here for our algorithms, MINIMAX and SCGD, we use an approximated Gaussian kernel based random Fourier features, we are no longer performing an exact optimization in this case. However, we can see in Table 2 that although we use an approximated kernel eventually when the restrictions on batch size are high our methods do a better job than the existing methods in the optimization of the marginal likelihood. The advantage of our methods in the optimization is also reflected in the test error (see Table 5).

We see that in both cases, the finite dimensional RKHS and the infinite dimensional RKHS, unlike the existing inference algorithms, there is no degradation in the result of our algorithms with the decreasing of the batch size. This property can be vital for inference on weak edge devices where the memory restrictions limit the possible batch size.

6 Conclusion

In many cases the covariance function of a GP is defined as an inner product between features of a finite and moderate dimension. In this case, the problem of minimizing the negative-log-marginal-likelihood has a shape of a standard ridge regression problem with a non standard regularization term in the form of the log-determinant of the covariance matrix of the representations plus $\sigma^{2}I$. In this work we developed two techniques that enable solving this problem with stochastic mini-batches which unlike the exiting methods does not depend on large batches in order to be exact. When the inference involves forward and backward passes of a feed forward neural net this property is of great importance and can be an enabler of such inference architectures on weak edge devices.

Limitations:

We remark that in comparison to BSGD, both MINIMAX and SCGD are more complex and require tuning of additional hyperparameters in order to achieve the minimal negative log marginal likelihood. In addition, the optimization advantage is not fully reflected in the test error, so cross-validation is still needed to select the best algorithm.

Appendix A Complexity Analysis

We analyze the computational complexity for each one of the algorithms. We ignore the complexity of for the computation of $\phi_{\alpha}\left(x_{i}\right)$ and $\frac{\phi_{\alpha}\left(x_{i}\right)}{\partial\alpha}$ as they are the same for all the algorithms. We let $d$ be the dimension of $\phi_{\alpha}\left(x_{i}\right)$ and $b$ be the batch size.

A.1 Minimax

Computational Complexity Forward

[TABLE]

Computational Complexity Backward

[TABLE]

Memory Complexity Forward

[TABLE]

Memory Complexity Backward

[TABLE]

A.2 SCGD

Computational Complexity Forward

[TABLE]

Computational Complexity Backward

[TABLE]

Memory Complexity Forward

[TABLE]

Memory Complexity Backward

[TABLE]

A.3 BSGD

[TABLE]

Computational complexity: $O\left(b^{3}\right)+O\left(b^{2}d\right)$.

Memory complexity: $O\left(b^{2}\right)+O\left(bd\right)$

Appendix B Missing Proofs

Theorem B.1**.**

For all $\lambda>0$, $V\in\mathbb{R}^{n\times d}$ and ${\mathbf{b}}\in\mathbb{R}^{d}$ we have that ${\mathbf{b}}^{T}\left(VV^{T}+\lambda I\right)^{-1}{\mathbf{b}}=\min_{{\mathbf{w}}}\frac{1}{\lambda}\left\|V{\mathbf{w}}-{\mathbf{b}}\right\|+\left\|{\mathbf{w}}\right\|^{2}$

Proof.

The minimizer of $\frac{1}{\lambda}\left\|V{\mathbf{w}}-{\mathbf{b}}\right\|^{2}+\left\|{\mathbf{w}}\right\|^{2}$ is given by $\hat{{\mathbf{w}}}=\left(V^{T}V+\lambda I\right)^{-1}V^{T}{\mathbf{b}}$. With this we can calculate

[TABLE]

∎