Locally induced Gaussian processes for large-scale simulation experiments

Abstract

Gaussian processes (GPs) serve as flexible surrogates for complex surfaces, but buckle under the cubic cost of matrix decompositions with big training data sizes. Geospatial and machine learning communities suggest pseudo-inputs, or inducing points, as one strategy to obtain an approximation easing that computational burden. However, we show how placement of inducing points and their multitude can be thwarted by pathologies, especially in large-scale dynamic response surface modeling tasks. As remedy, we suggest porting the inducing point idea, which is usually applied globally, over to a more local context where selection is both easier and faster. In this way, our proposed methodology hybridizes global inducing point and data subset-based local GP approximation. A cascade of strategies for planning the selection of local inducing points is provided, and comparisons are drawn to related methodology with emphasis on computer surrogate modeling applications. We show that local inducing points extend their global and data subset component parts on the accuracy–computational efficiency frontier. Illustrative examples are provided on benchmark data and a large-scale real-simulation satellite drag interpolation problem.

Appendices

A IMSE and ALC overview

As mentioned in Sect. 2.3, variance-based sequential design criteria are better aligned with the goal of generating accurate GP predictions than using the likelihood. We consider variations on integrated mean-squared error (IMSE) over a domain $\mathcal {X}$, with smaller being better:

$$\begin{aligned} \mathrm {(IMSE)} \quad \quad I = \int _{\tilde{\mathbf {x}}\in \mathcal {X}}\sigma ^2(\tilde{\mathbf {x}}) \; \mathrm{d}\tilde{\mathbf {x}}. \end{aligned}$$

Choose $\sigma ^2(\cdot ) \equiv \sigma _N^2(\cdot )/\nu$ from Eq. (2), and I may be used to optimize the N coordinates of $\mathbf {X}_N$, or to choose the next ($N+1^\mathrm {st}$) one ($\tilde{\mathbf {x}}_{N+1})$ in a sequential setting.⁹ Closed-form expressions are available for rectangular $\mathcal {X}$ and common kernels (e.g., Ankenman et al. 2010; Anagnostopoulos and Gramacy 2013; Burnaev and Panov 2015; Leatherman et al. 2018). Analytic derivatives $\frac{\partial I}{\partial \tilde{\mathbf {x}}_{N+1}}$ facilitate numerical optimization (Binois et al. 2019; Gramacy 2020, Chapters 4 & 10). Approximations are common otherwise (Gramacy and Lee 2009; Gauthier and Pronzato 2014; Gorodetsky and Marzouk 2016; Pratola et al. 2017).

An analogue active learning heuristic from Cohn (1993), dubbed ALC, instead targets variance aggregated over a discrete reference set $\mathcal {X}$, originally for neural network surrogates:

$$\begin{aligned} \mathrm {(ALC)} \quad \quad \varDelta \sigma ^2 =\sum _{\tilde{\mathbf {x}}\in \mathcal {X}}\sigma ^2(\tilde{\mathbf {x}})-\sigma _{\mathrm {new}}^2(\tilde{\mathbf {x}}). \end{aligned}$$

Seo et al. (2000) ported ALC to GPs taking $\sigma ^2(\cdot ) = \sigma _N^2(\cdot )$ and $\sigma ^2_{\mathrm {new}}(\cdot ) \equiv \sigma ^2_{N+1}(\cdot )$. If discrete and volume-based $\mathcal {X}$ are similar, then $\varDelta \sigma ^2 \approx c - I$, where c is constant on $\mathbf {x}_{N+1}$. Discrete $\varDelta \sigma ^2$ via ALC is advantageous in transductive learning settings (Vapnik 2013), where $\mathcal {X}$ can be matched with a testing set. Otherwise, analytic I via IMSE may be preferred.

Against that backdrop, we propose employing ALC and IMSE to select inducing points ${\bar{\mathbf {X}}}_M$. To our knowledge, using such variance-based criteria is novel in the literature on the selection of inducing points. The criteria below are framed sequentially, for an $M+1^\mathrm {st}$ point given M collected already. Although we prefer this greedy approach—optimizing d coordinates one-at-a-time rather than Md all at once in a surface with many equivalent locally optimal configurations due to label-switching—either criteria is easily re-purposed for an all-at-once application. Under the diagonal-corrected Nyström approximation (3) and assuming coded $\mathcal {X} = [0,1]^d$,

$$\begin{aligned} \begin{aligned} \text {ALC}_N^{(M+1)}&= \text {ALC}({\bar{\mathbf {x}}}_{M+1}; \mathbf {X}_N, \mathbf {Y}_N, \mathcal {X},{\bar{\mathbf {X}}}_M)\\&=c-\sum _{\tilde{\mathbf {x}}\in \mathcal {X}}\sigma _{M+1,N}^2(\tilde{\mathbf {x}}) , \quad \text{ and } \\ \text {IMSE}_N^{(M+1)}&= \text {IMSE}({\bar{\mathbf {x}}}_{M+1},\mathbf {X}_N, \mathbf {Y}_N, \mathcal {X},{\bar{\mathbf {X}}}_M)\\&=E-\text {tr}\Big \{\Big (\mathbf {K}_{M+1}^{-1}-\mathbf {Q}_{M+1}^{-1(N)}\Big )\mathbf {W}_{M+1}\Big \}, \end{aligned} \end{aligned}$$

(13)

where $E=\int _{\tilde{x}\in \mathcal {X}}k(\tilde{\mathbf {x}},\tilde{\mathbf {x}})\mathrm{d}\tilde{\mathbf {x}}$ and $\mathbf {W}_{M+1}$ is $(M+1)\times (M+1)$ via $w({\bar{\mathbf {x}}}_i,{\bar{\mathbf {x}}}_j)=\int _{\tilde{\mathbf {x}}\in \mathcal {X}}k({\bar{\mathbf {x}}}_i,\tilde{\mathbf {x}})k({\bar{\mathbf {x}}}_j,\tilde{\mathbf {x}})\mathrm{d}\tilde{\mathbf {x}}$ for $i,j\in \{1,\dots ,M+1\}$. This derivation is similar to the wIMSE calculations (10) and (11) following Binois et al. (2019).

Fig. 9

In both panels: $N=200$ training data points (black dots) and $M=19$ inducing points (blue dots), selecting the twentieth one (green) by two criteria: a variational lower bound of the log-likelihood; b ALC/IMSE. Yellow is higher/red lower. (Color figure online)

To explore inducing point optimization, consider Herbie’s tooth (Lee et al. 2011) described in Sect. 2.3. Figure 9 shows variational lower-bound of the log-likelihood (left) and ALC/IMSE surfaces (right) for ${\bar{\mathbf {x}}}_{20}$ given a modestly sized training data set $(\mathbf {X}_N,\mathbf {Y}_N)$ of size $N=200$. Similarities in the two surfaces are apparent. Many low/red areas coincide, but the optimizing locations (green dots), found via multi-start local optimization with identically fixed kernel hyperparameters, do not. Even after taking great care to humbly restrict searchers, e.g., from crossing ${\bar{\mathbf {X}}}_M$ locations, sometimes upward of 1000 evaluations were required to achieve convergence. Consequently, quadratic ALC/IMSE is faster.

1.1 A.1 Derivations of wIMSE and its gradient

For the predictive location $\mathbf {x}^*$, assign weight $k_\theta (\tilde{\mathbf {x}},\mathbf {x}^\star )$ and consider squared exponential kernel $k_\theta (\cdot , \cdot )$ with isotropic lengthscale (1). The following is based on predictive variance (8) and expectation of the quadratic form of a random vector (Binois et al. 2019, Section 3.1).

$$\begin{aligned}&\text {wIMSE}({\bar{\mathbf {x}}}_{m+1},\mathcal {X},{\bar{\mathbf {X}}}_m, \mathbf {X}_n,\theta , \mathbf {x}^\star ) \\&\quad = \int _{\tilde{\mathbf {x}}\in \mathcal {X}}k_\theta (\tilde{\mathbf {x}},\mathbf {x}^\star )\sigma ^2_{n,m+1}(\tilde{\mathbf {x}})\mathrm{d}\tilde{\mathbf {x}}\\&\quad = \int _{\tilde{\mathbf {x}}\in \mathcal {X}}k_\theta (\tilde{\mathbf {x}},\mathbf {x}^\star )\Big (k_\theta (\tilde{\mathbf {x}},\tilde{\mathbf {x}}) + \epsilon _K \\&\qquad - k_\theta (\tilde{\mathbf {x}}, {\bar{\mathbf {X}}}_{m+1})\Big [\mathbf {K}^{-1}_{m+1}-\mathbf {Q}^{-1(n)}_{m+1}\Big ]k_\theta (\tilde{\mathbf {x}}, {\bar{\mathbf {X}}}_{m+1})^\top \Big ) \mathrm{d}\tilde{\mathbf {x}}\\&= \prod _{k=1}^D \Big ((1+\epsilon _K)\int _{a_k}^{b_k}k_\theta (\tilde{\mathbf {x}}_k,\mathbf {x}^\star _k) \mathrm{d}\tilde{\mathbf {x}}_k \\&\qquad - \int _{a_k}^{b_k}k_\theta (\tilde{\mathbf {x}}_k, \mathbf {x}^\star _k)^{1/2}k_\theta (\tilde{\mathbf {x}}_k, {\bar{\mathbf {X}}}_{m+1,k})\Big [\mathbf {K}^{-1}_{m+1}-\mathbf {Q}^{-1(n)}_{m+1}\Big ]\\&\qquad \times k_\theta (\tilde{\mathbf {x}}_k, {\bar{\mathbf {X}}}_{m+1,k})^\top k_\theta (\tilde{\mathbf {x}}_k,\mathbf {x}^\star _k)^{1/2} \mathrm{d}\tilde{\mathbf {x}}_k\Big )\\&\quad =\frac{\sqrt{\theta \pi }(1+\epsilon _K)^D}{2}\prod _{k=1}^D\left( \text {erf}\left\{ \frac{\mathbf {x}^\star -a_k}{\sqrt{\theta }}\right\} -\text {erf}\left\{ \frac{\mathbf {x}^\star -b_k}{\sqrt{\theta }}\right\} \right) \\&\qquad -\text {tr}\Big \{(\mathbf {K}^{-1}_{m+1}-\mathbf {Q}^{-1}_{m+1})\mathbf {W}^\star _{m+1}\Big \}\\ \end{aligned}$$

where $\mathbf {W}^\star _{m+1}=\prod _{k=1}^D \mathbf {W}^\star _{m+1,k}$ in D dimensions. The entry in the $i^{\text {th}}$ row and $j^{\text {th}}$ column of $\mathbf {W}^\star _{m+1,k}$ is

$$\begin{aligned} w^{\star (i,j)}_{m+1,k}&=w^\star _{m+1,k}({\bar{\mathbf {x}}}_{i,k}, {\bar{\mathbf {x}}}_{j,k}) \\&= \int _{a_k}^{b_k}k_\theta (\tilde{\mathbf {x}}_k,\mathbf {x}^\star _k)k_\theta (\tilde{\mathbf {x}}_k,{\bar{\mathbf {x}}}_{i,k})k_\theta (\tilde{\mathbf {x}}_k,{\bar{\mathbf {x}}}_{j,k})\mathrm{d}\tilde{\mathbf {x}}_k \\&=\int _{a_k}^{b_k}\!\!\!\!\exp \left\{ -\frac{(\tilde{\mathbf {x}}_k-\mathbf {x}^\star _k)^2+(\tilde{\mathbf {x}}_k\!-\!{\bar{\mathbf {x}}}_{i,k})^2+(\tilde{\mathbf {x}}_k\!-\!{\bar{\mathbf {x}}}_{j,k})^2}{\theta }\right\} \mathrm{d}\tilde{\mathbf {x}}_k\\&=\sqrt{\frac{\pi \theta }{12}}\exp \Big \{\frac{2}{3\theta }({\bar{\mathbf {x}}}_{i,k}\mathbf {x}^\star _k + {\bar{\mathbf {x}}}_{j,k}\mathbf {x}^\star _k + {\bar{\mathbf {x}}}_{i,k}{\bar{\mathbf {x}}}_{j,k} \\&\quad - \mathbf {x}^{*2}_k - {\bar{\mathbf {x}}}_{i,k}^2 - {\bar{\mathbf {x}}}_{j,k}^2 )\Big \}\times \\&\quad \left( \text {erf}\left\{ \frac{\iota ^{(i,j)}_{k}-3a_k}{\sqrt{3\theta }}\right\} -\text {erf}\left\{ \frac{\iota ^{(i,j)}_{k}-3b_k}{\sqrt{3\theta }}\right\} \right) \end{aligned}$$

where $\iota ^{(i,j)}_{k}=\mathbf {x}^\star _{k}+{\bar{\mathbf {x}}}_{i,k}+{\bar{\mathbf {x}}}_{j,k}$. ${\bar{\mathbf {x}}}_{i,k},{\bar{\mathbf {x}}}_{j,k}$ are entries from the $i^{\text {th}}$ and $j^{\text {th}}$ rows and $k^{\text {th}}$ column of ${\bar{\mathbf {X}}}_{m+1}$ ($i, j \in \{1,\dots ,m+1\}$) and $\mathbf {x}^\star _k$ is the $k^{\text {th}}$ coordinate of $\mathbf {x}^\star$.

The gradient of weighted integrated mean-squared error with respect to the $k^{\text {th}}$ dimension of ${\bar{\mathbf {x}}}_{m+1}$ is:

$$\begin{aligned}&\frac{\partial \text {wIMSE}({\bar{\mathbf {x}}}_{m+1},\mathcal {X},{\bar{\mathbf {X}}}_m, \mathbf {X}_n,\theta , \mathbf {x}^\star )}{\partial {\bar{\mathbf {x}}}_{m+1,k}} \\&\quad = - \text {tr}\left\{ \left( \frac{\partial \mathbf {K}^{-1}_{m+1}}{\partial {\bar{\mathbf {x}}}_{m+1,k}}-\frac{\partial \mathbf {Q}^{-1(n)}_{m+1}}{\partial {\bar{\mathbf {x}}}_{m+1,k}}\right) \mathbf {W}^\star _{m+1}\right\} \\&\qquad - \text {tr}\left\{ \left( \mathbf {K}^{-1}_{m+1}-\mathbf {Q}^{-1(n)}_{m+1}\right) \frac{\partial \mathbf {W}^\star _{m+1}}{\partial {\bar{\mathbf {x}}}_{m+1,k}}\right\} \\&\quad = \text {tr}\Big \{\Big (\mathbf {K}^{-1}_{m+1}\frac{\partial \mathbf {K}_{m+1}}{\partial {\bar{\mathbf {x}}}_{m+1,k}}\mathbf {K}^{-1}_{m+1}\\&\qquad -\mathbf {Q}^{-1(n)}_{m+1}\frac{\partial \mathbf {Q}^{(n)}_{m+1}}{\partial {\bar{\mathbf {x}}}_{m+1,k}}\mathbf {Q}^{-1(n)}_{m+1}\Big )\mathbf {W}^\star _{m+1}\Big \}\\&\qquad - \text {tr}\left\{ \left( \mathbf {K}^{-1}_{m+1}-\mathbf {Q}^{-1(n)}_{m+1}\right) \frac{\partial \mathbf {W}^\star _{m+1}}{\partial {\bar{\mathbf {x}}}_{m+1,k}}\right\} . \end{aligned}$$

In the matrix $\frac{\partial \mathbf {W}^\star _{m+1}}{d{\bar{\mathbf {x}}}_{m+1,k}}$, all entries are zero except the row/column that corresponds to the row of ${\bar{\mathbf {X}}}_{m+1}$ that contains ${\bar{\mathbf {x}}}_{m+1}$, which we place in the last $m+1^{\text {st}}$ row. For the nonzero entries in $\frac{\partial \mathbf {W}^\star _{m+1}}{\partial {\bar{\mathbf {x}}}_{m+1,k}}$, we re-express them as

$$\begin{aligned} \frac{\partial w^\star _{m+1}({\bar{\mathbf {x}}}_i,{\bar{\mathbf {x}}}_{m+1})}{\partial {\bar{\mathbf {x}}}_{m+1,k}}=\frac{\partial w^{\star (i,m+1)}_{m+1}}{\partial {\bar{\mathbf {x}}}_{m+1,k}}\prod _{k'=1,k'\ne k}^D w^{\star (i,m+1)}_{m+1,k'} \end{aligned}$$

where

$$\begin{aligned}&\frac{\partial w^{\star (i,m+1)}_{m+1}}{\partial {\bar{\mathbf {x}}}_{m+1,k}}\\&\quad =\sqrt{\frac{\pi \theta }{12}}\text {exp}\Bigg \{\frac{2}{3\theta }\Big ({\bar{\mathbf {x}}}_{i,k}\mathbf {x}^\star _{k} + {\bar{\mathbf {x}}}_{m+1,k}\mathbf {x}^\star _{k} + {\bar{\mathbf {x}}}_{i,k'}{\bar{\mathbf {x}}}_{m+1,k} \\&\qquad - \mathbf {x}^{*2}_{k} - {\bar{\mathbf {x}}}_{i,k}^2 - {\bar{\mathbf {x}}}_{m+1,k}^2 \Big )\Bigg \}\\&\qquad \times \Bigg [\frac{2}{3\theta }(\mathbf {x}_{k}^\star -2{\bar{\mathbf {x}}}_{m+1,k}-{\bar{\mathbf {x}}}_{i,k})\\&\qquad \times \left( \text {erf}\left\{ \frac{\iota ^{(i,m+1)}_{k}-3a_{k}}{\sqrt{3\theta }}\right\} -\text {erf}\left\{ \frac{\iota ^{(i,m+1)}_{k}-3b_{k}}{\sqrt{3\theta }}\right\} \right) \\&\qquad +\frac{2}{\sqrt{3\pi \theta }}\Bigg (\text {exp}\Big \{-\frac{(\iota ^{(i,m+1)}_{k}-3a_{k})^2}{3\theta }\Big \}\\&\qquad -\text {exp}\Big \{-\frac{(\iota ^{(i,m+1)}_{k}-3b_{k})^2}{3\theta }\Big \}\Bigg )\Bigg ]. \end{aligned}$$

Working with $\mathbf {K}_{m+1}$ and $\mathbf {Q}_{m+1}^{(n)}$ is cubic in m, yet even that is overkill. Thrifty evaluation of Eqs. (10–12) lies in construction of $\mathbf {Q}_{m+1}^{(n)}$ which is equivalent to $\mathbf {K}_{m+1} + \mathbf {k}^\top _{n,m+1}\varvec{\Gamma }_{n,m+1}$. Evaluating $\mathbf {k}^\top _{n,m+1}\varvec{\Gamma }_{n,m+1}$ requires $2n-1$ products for each of $(m+1)^2$ entries, incurring costs in $\mathcal {O}(m^2n)$ flops. Assuming $n\gg m$, this dominates the $\mathcal {O}(m^3)$ cost of decomposition.

More time can be saved through partitioned inverse (Barnett 1979) sequential updates to $\mathbf {K}^{-1}_{m+1}$ after the new ${\bar{\mathbf {x}}}_{m+1}$ is chosen, porting LAGPs frugal updates to the LIGP context. Writing $\mathbf {K}_{m+1}$ as an m-submatrix with new $m+1^\mathrm {st}$ column gives

$$\begin{aligned} \begin{aligned} \mathbf {K}_{m+1}&=\begin{bmatrix} \mathbf {K}_m &{}\quad \mathbf {k}_m({\bar{\mathbf {x}}}_{m+1}) \\ \mathbf {k}_m({\bar{\mathbf {x}}}_{m+1})^\top &{}\quad k_\theta ({\bar{\mathbf {x}}}_{m+1},{\bar{\mathbf {x}}}_{m+1}) \end{bmatrix} \quad \text{ so } \text{ that } \\ \mathbf {K}_{m+1}^{-1}&= \begin{bmatrix} \mathbf {K}_m^{-1}+\rho \varvec{\eta } \varvec{\eta }^\top &{}\quad \varvec{\eta } \\ \varvec{\eta }^\top &{}\quad \rho ^{-1}&{} \end{bmatrix} \end{aligned} \end{aligned}$$

(14)

using $\rho = k_\theta ({\bar{\mathbf {x}}}_{m+1},{\bar{\mathbf {x}}}_{m+1}) -\mathbf {k}_m^\top ({\bar{\mathbf {x}}}_{m+1})\mathbf {K}_m^{-1}\mathbf {k}_m({\bar{\mathbf {x}}}_{m+1})$ and m-length column vector $\varvec{\eta }=-\rho ^{-1}\mathbf {K}_m^{-1}\mathbf {k}_m({\bar{\mathbf {x}}}_{m+1})$. Updating $\mathbf {K}_{m+1}^{-1}$ requires calculation of $\rho$, $\varvec{\eta }$, and $\varvec{\eta } \varvec{\eta }^\top$, each of which is in $\mathcal {O}(m^2)$. Thus, we reduce the computational complexity of $\mathbf {K}_{m+1}^{-1}$ from $\mathcal {O}\left( m^3\right)$ to $\mathcal {O}\left( m^2\right)$. Similar partitioning provides sequential updates to $\varOmega _n^{(m+1)}$, a diagonal matrix:

$$\begin{aligned} \varOmega _n^{(m+1)}&=\text {Diag}\!\left( \mathbf {K}_{n}+\epsilon _K\mathbb {I}_n-\mathbf {k}_{n,m+1}\mathbf {K}_{m+1}^{-1} \mathbf {k}_{n,m+1}^\top \!\right) \nonumber \\&=\varOmega _n^{(m)}-\rho ^{-1}\text {Diag}\left\{ \zeta \zeta ^\top \right\} \end{aligned}$$

(15)

where $\zeta =\mathbf {k}_{nm}\mathbf {K}_m^{-1}\mathbf {k}_m({\bar{\mathbf {x}}}_{m+1})-\mathbf {k}_n({\bar{\mathbf {x}}}_{m+1})$. Updates of $\varOmega _n^{(m+1)}$ without partitioning, driven by matrix–vector product(s) $\mathbf {k}_{n,m+1}\mathbf {K}_{m+1}^{-1} \mathbf {k}_{n,m+1}^\top$ involve $m^2n$ flops. Using (15) reduces that to $\mathcal {O}(mn)$.

Unlike in Eq. (14), $\mathbf {Q}_{m}^{(n)}$ cannot be trivially augmented to construct $\mathbf {Q}_{m+1}^{(n)}$ due to the presence of $\varOmega _n^{(m)}$ which is also embedded in $\mathbf {Q}_{m}^{(n)}$. Yet there are some time savings to be found in the partitioned inverse

$$\begin{aligned} \begin{aligned} \mathbf {Q}_{m+1}^{(n)}&=\begin{bmatrix} \mathbf {Q}_{m*}^{(n)} &{} \varvec{\gamma }({\bar{\mathbf {x}}}_{m+1}) \\ \varvec{\gamma }({\bar{\mathbf {x}}}_{m+1})^\top &{} \psi ({\bar{\mathbf {x}}}_{m+1}) \end{bmatrix}\\ \mathbf {Q}_{m+1}^{-1(n)}&=\begin{bmatrix} \mathbf {Q}_{m*}^{-1(n)}+\upsilon \varvec{\xi } \varvec{\xi }^\top &{}\quad \varvec{\xi }\\ \varvec{\xi }^\top &{}\quad \upsilon ^{-1} \end{bmatrix} \end{aligned} \end{aligned}$$

(16)

with $\mathbf {Q}_{m*}^{(n)}= \mathbf {K}_m + \mathbf {k}_{nm}^\top \varOmega _n^{(m+1)-1} \mathbf {k}_{nm}$ built via updated values of $\varOmega _n^{(m+1)}$, $\varvec{\gamma }({\bar{\mathbf {x}}}_{m+1}) = \mathbf {k}_m({\bar{\mathbf {x}}}_{m+1})+\mathbf {k}^\top _{nm}\varOmega _n^{-1(m+1)} \mathbf {k}_n({\bar{\mathbf {x}}}_{m+1})$, $\psi ({\bar{\mathbf {x}}}_{m+1}) = k_\theta ({\bar{\mathbf {x}}}_{m+1},{\bar{\mathbf {x}}}_{m+1})+k_n({\bar{\mathbf {x}}}_{m+1})^\top \varOmega _n^{-1(m+1)}k_n({\bar{\mathbf {x}}}_{m+1})$, $\upsilon = \psi ({\bar{\mathbf {x}}}_{m+1})-\varvec{\gamma }({\bar{\mathbf {x}}}_{m+1})^\top \mathbf {Q}_{m*}^{-1(n)}\varvec{\gamma }({\bar{\mathbf {x}}}_{m+1})$ and $\xi = -\upsilon ^{-1}\mathbf {Q}_{m*}^{-1(n)}\varvec{\gamma }({\bar{\mathbf {x}}}_{m+1})$. Similar to $\mathbf {Q}_{m}^{(n)}$, calculating $\mathbf {Q}_{m*}^{(n)}$ requires in flops in $\mathcal {O}(m^2n)$. Consequently, the entire scheme can be managed in $\mathcal {O}(m^2n)$.

B Determining neighborhood size

Little attention is paid in the literature to the choosing the number of (global) inducing points (Seeger et al. 2003; Titsias 2009b; Azzimonti et al. 2016) relative to problem size (N, d), except on computational grounds—smaller M is better. The same is true for local neighborhood size n in LAGP. Although there is evidence that the laGP default of $n=50$ is too small (Gramacy 2016), especially with larger input dimension d, cubically growing expense in n limits the efficacy of larger n in practice. With local inducing points this is mitigated through cubic-in-m proxies, allowing larger local neighborhoods, thus implying more latitude to explore/choose good (m, n) combinations.

Toward that end, we considered a coarse grid of (m, n) and predictive RMSEs on Herbie’s tooth ($d=2$) and borehole $(d=8)$ toy problems. Setup details are identical to descriptions in Sects. 3.2 and 5.2, respectively, and we used the qNorm ($\varPhi ^{-1}$) template throughout. An LHS testing set of size $N'=1000$ was used to generate the response surfaces of RMSEs reported in Fig. 10. These are shown in log space for a more visually appealing color scheme, and were obtained after GP smoothing to remove any artifacts from random testing. Grid elements where $m>n$ were omitted from the simulation on the grounds that there are no run-time benefits to those choices.

Fig. 10

$\log (\mathrm {RMSE})$ over inducing points m and neighborhood n: Herbie’s tooth (top) and borehole (bottom). (Color figure online)

Observe that both surfaces are fairly flat across a wide swath of m, excepting quick ascent (decrease in accuracy) for smaller numbers of inducing points in the left panel. The situation is similar for n. Best settings are apparently input-dimension dependent. Numbers of inducing points as low as $m=10$ seems sufficient in 2d (top panel), whereas $m=80$ is needed in 8d. For borehole, it appears that larger neighborhoods n are better, perhaps because the response surface is very smooth and the likelihood prefers long lengthscales (Gramacy 2016). A setting like $n=150$ seems to offer good results without being too large. The situation is different for Herbie’s tooth. Here, larger n has deleterious effects. Its non-stationary nature demands reactivity which is proffered by smaller local neighborhood. A setting of $n=100$ looks good.

These are just two problems, and it is clearly not reasonable to grid-out (m, n) space for all future applications. But nevertheless, we have found that these rules of thumb port well to our empirical work in Sect. 5. Our satdrag example ($d=8$) and classic $d=21$ benchmark work well with the settings found for borehole, for example. Some ideas for automating the choice of (m, n) are discussed in Sect. 6.