Surrogate modeling for Bayesian optimization beyond a single Gaussian process

In plain words

Bayesian optimization finds the best setting of something that is slow to test, such as a model's learning rate or layer sizes. It builds a guess of how good each setting is from the tests so far, then runs the next test where the guess looks most promising. The guess is usually a Gaussian process, and a Gaussian process needs a kernel: a rule for how smooth the guess is and how far its influence reaches. Choosing that kernel well is hard, and it is often done by hand.

This paper keeps several candidate kernels at once, each with its own guess of the function. Each candidate gets a weight that rises when it predicts the new results well. To choose the next test, the method picks one candidate at random in proportion to its weight, draws one plausible version of the function from that candidate, and tests the spot where that version is highest. Random features make drawing these functions cheap, and the same scheme runs several tests at once.

The authors also prove that the method's cumulative regret grows more slowly than the number of tests, so it keeps approaching the best setting. The proofs assume a smoothness condition on the functions.

TL;DR

EGP-TS is a Bayesian optimization method whose surrogate is a weighted mixture of Gaussian processes, one per kernel in a fixed dictionary §3. The next query is the maximiser of a function drawn from the mixture by Thompson sampling, with random features (50 per kernel) for scalable sampling §3.1, §5. The authors report lower simple regret than single-GP Thompson sampling, GP-EI and two bandit ensembles, on synthetic functions and robot pushing fig. 1–3. Use it when you cannot pick one kernel in advance and you need sublinear Bayesian regret. Expect to check the code before reusing it (see Paper vs code).

Problem

Bayesian optimization of an expensive black-box function needs a surrogate model and an acquisition rule, and both carry design choices §1, §2. With a single Gaussian process, the kernel form is usually preselected from domain knowledge, and the authors argue this is the design step they want to remove §1. Existing kernel-selection methods for Gaussian processes run in batch mode on many samples, so they do not suit the small-data setting of Bayesian optimization §1, §5. Parallel evaluation usually adds extra parameters or diversity rules to upper-confidence-bound or expected-improvement rules, while Thompson sampling needs none §1. The paper states the goal as a method that needs the least tuning, runs in parallel and has convergence guarantees §1.

Method

Ensemble prior §3. Take a dictionary of M kernels \kappa^1, \ldots, \kappa^M. Each GP m has prior weight w_0^m = \Pr(i=m). The function prior is the weighted sum of the component GP priors, with \sum_m w_0^m = 1 eq. 7. In the experiments, w_0^m = 1/M §5.

Posterior §3. The posterior is a mixture of the component posteriors, with per-component weight w_t^m \propto w_0^m\, p(\mathcal{D}_t \mid i=m), where p(\mathcal{D}_t \mid i=m) is the marginal likelihood of the data under GP m eq. 8, eq. 9.

Thompson sampling from the mixture §3. Sample a component m_t from the categorical distribution with weights \mathbf{w}_t, then sample a function from that component's posterior eq. 10. The paper uses random features for the function sample, so the second step becomes a Gaussian draw of the feature weights §3.1.

Random features §2.1, §3.1, §8. For a shift-invariant kernel with normalised spectral density \pi_{\bar\kappa}, draw D frequencies \mathbf{v}_j \sim \pi_{\bar\kappa} and set the feature map to \frac{1}{\sqrt{D}} times the stacked \sin(\mathbf{v}_j^\top \mathbf{x}) and \cos(\mathbf{v}_j^\top \mathbf{x}) terms eq. 5. The function is then linear in the 2D feature weights \boldsymbol\theta, with prior \mathcal{N}(\mathbf{0}, \sigma_\theta^2 \mathbf{I}_{2D}), where \sigma_\theta^2 is the kernel magnitude [eq. 5, §3.1]. The model is linear-Gaussian, so its posterior is Gaussian with closed-form updates §3.1, §8.

Query selection §3.1. Given the sampled component m_t and a weight vector \tilde{\boldsymbol\theta}_t \sim \mathcal{N}(\hat{\boldsymbol\theta}^{m_t}_t, \boldsymbol\Sigma^{m_t}_t), the next query maximises \boldsymbol\phi^{m_t}(\mathbf{x})^\top \tilde{\boldsymbol\theta}_t over \mathcal{X} eq. 12. The paper says a gradient-based solver can do this because the objective has a closed form; it does not say how to start the solver §3.1.

Weight and posterior updates §3.1. After observing y_{t+1} at \mathbf{x}_{t+1}, each component's predictive mean is \hat y^m = \boldsymbol\phi^{m\top}\hat{\boldsymbol\theta}^m_t and its predictive variance is \boldsymbol\phi^{m\top}\boldsymbol\Sigma^m_t\boldsymbol\phi^m + \sigma_n^2 §3.1. The weights are multiplied by the Gaussian likelihood of y_{t+1} under each component and then normalised eq. 13. Each component's posterior is then updated with a rank-one step: the mean moves by the prediction error times \sigma^{-2}\boldsymbol\Sigma\boldsymbol\phi, and the covariance shrinks by \sigma^{-2}\boldsymbol\Sigma\boldsymbol\phi\boldsymbol\phi^\top\boldsymbol\Sigma eq. 14, eq. 15a, eq. 15b.

Initialization and re-initialization §3.1, §5, alg. 1. Start with t_0 random points, which is 10 in the experiments §5. For each component, fit the kernel hyperparameters by maximising the marginal likelihood, draw new frequencies, compute the batch posterior and the batch weight [§8 batch formulas, eq. 9]. The text says the hyperparameters and weights are refitted "every few iterations" from all acquired data §3.1; the experiments refit every 50 iterations for EGP-TS and every iteration for the single-GP baselines §5.

Batch and parallel queries §3, alg. 1, alg. 2. With K workers, the synchronous version draws K components and K weight vectors from the current mixture and posteriors, evaluates the K queries, and then applies the updates for all K pairs alg. 1. The asynchronous version (alg. 2 in the appendix) issues a new query whenever a worker finishes, and updates the posterior at that time §3.2, alg. 2. Theorems 2 and 3 give regret bounds for these two settings §4.

Algorithm 1 in pseudocode. The block below follows the paper's update rules for one query per round (K=1). The paper's own algorithm is written as pseudocode for K workers alg. 1; this version is ours.

# EGP-TS, one query per round (K = 1)  [alg. 1] [§3.1]
# gps[k]: features (map), theta (2D,), Sigma (2D, 2D), noise (scalar)
# weights: mixture weights over the M components, sum to 1  [eq. 9]

def egp_ts_round(gps, weights, objective):
    # Pick one component from the mixture weights  [eq. 10]
    m = sample_index(weights)
    # Draw feature weights from that component's posterior  [eq. 10]
    theta_tilde = sample_gaussian(gps[m].theta, gps[m].Sigma)
    # Query the maximiser of the sampled feature function  [eq. 12]
    x_next = argmax_over_domain(
        lambda x: dot(gps[m].features(x), theta_tilde)
    )
    y_next = objective(x_next)
    for k in range(M):
        # Predictive mean and variance of component k at x_next  [§3.1]
        phi = gps[k].features(x_next)                  # (2D,)
        mean_k = dot(phi, gps[k].theta)
        var_k = dot(phi, gps[k].Sigma @ phi) + gps[k].noise
        # Bayes update of the component weight  [eq. 13]
        weights[k] = weights[k] * gaussian_pdf(y_next, mean_k, var_k)
        # Rank-one posterior update of theta and Sigma  [eq. 15a, eq. 15b]
        gain = gps[k].Sigma @ phi / var_k              # (2D,)
        gps[k].theta = gps[k].theta + gain * (y_next - mean_k)
        gps[k].Sigma = gps[k].Sigma - var_k * outer(gain, gain)
    # Renormalise the weights so they sum to one  [eq. 13]
    weights = weights / sum(weights)
    return x_next, y_next

Re-initialization, which the paper runs every few iterations (every 50 in the experiments), rebuilds every component from all acquired pairs. The paper does not say whether the marginal likelihood in the weight is computed with the exact GP or with the feature model.

# Re-initialization from all acquired pairs  [§3.1] [alg. 1] [§5]
# prior_weights[k] = w_0^k = 1 / M  [§3, §5]

def reinitialize(gps, weights, prior_weights, data_x, data_y):
    for k in range(M):
        # Refit kernel hyperparameters by marginal likelihood  [alg. 1]
        gps[k].params = fit_marginal_likelihood(data_x, data_y, k)
        # Draw new frequencies from the spectral density  [§8]
        gps[k].features = draw_features(gps[k].params, D)
        # Batch posterior of the feature weights  [§8 batch formulas]
        gps[k].theta, gps[k].Sigma = batch_posterior(gps[k], data_x, data_y)
        # Batch weight, proportional to prior times evidence  [eq. 9]
        weights[k] = prior_weights[k] * evidence(gps[k], data_x, data_y)
    # Normalise the weights so they sum to one  [eq. 9]
    weights = weights / sum(weights)
    return gps, weights

Hyperparameters

Values stated in the paper and, where it differs or is silent, the value in the released notebook. Notebook line numbers are lines of the .ipynb JSON file.

Table 1: feasible ranges for the classifier hyperparameters tab. 1.

Results

Author-reported throughout. The main text gives no numeric regret or accuracy values; the comparisons are read from plots, and the figure files are not reproduced here. The regret results are theorems, not experiments.

Synthetic functions §5.1. Benchmarks are Ackley-5d, Zakharov, Drop-wave and Eggholder, with the sign flipped from the usual minimisation form. The candidate domains and the optimum values are in tab. 3. With the mixed dictionary, the text reports that EGP-TS gives substantially lower simple regret than GP-TS with each single kernel, and lower regret than BanditBO and EXP3BO fig. 1. On Zakharov, GP-EI is better than the GP-TS baselines, but EGP-TS is better than GP-EI §5.1. With the 11-RBF dictionary, the text reports that EGP-TS also beats the fully Bayesian GP-TS fig. 2.

Robot pushing §5.2. Two tasks, with 3 and 4 inputs, each averaged over 10 goal locations. With the mixed dictionary, EGP-TS beats GP-EI, GP-TS with a fixed kernel and the two ensembles fig. 3a. With the RBF dictionary it is also better fig. 3b. The text says EGP-TS beats fully Bayesian GP-TS in simple regret and runs much faster §5.2.

Hyperparameter tuning §5.3. Tuning a 2-layer FNN, an SVM and gradient boosting on Breast Cancer, Iris, Transportation and Wine, with validation accuracy as the objective. The text says EGP-TS is best in most cases, for both the mixed and RBF dictionaries fig. 4–9. The per-task numbers are not in the text.

EGP-EI §5.4. Using the same ensemble with expected improvement, the text reports EGP-EI beats GP-EI on three of four synthetic functions and on both robot tasks fig. 10, fig. 11.

Runtime §9.3, tab. 2. Average run times per method. The table gives no unit, and the text says EGP-TS needs slightly more time than single-GP TS.

Table 3: synthetic benchmarks tab. 3.

Theory §4, §7. The cumulative Bayesian regret over T evaluations has the bound below. It uses the maximum information gain \gamma_T \le \tilde O(T^c) with 0 \le c < 1 §4, Lemma 1, and Assumption 1 on sample-path smoothness §4.

\mathcal{BR}(T) \le c_1 \sqrt{M T^{c+1} \log T} + 2\sigma_n \sqrt{M T \log T} + c_2 §4, Theorem 1

with c_1 = (2+\sqrt{d})\,(2/\log(1+\sigma_n^{-2}))^{1/2} and c_2 = 9MB + \pi^2 d/6 + \sqrt{2\pi}\,M/12, where B is the expected supremum of a GP sample path §4, Lemma 3 in §7. The authors say the bound is sublinear for 0 \le c < 1, so the average regret per evaluation goes to zero §4.

For asynchronous workers the bound is \sqrt{\rho_K} times the sequential first term, where \rho_K \ge 1 comes from a posterior-variance ratio [§4, Theorem 2, Lemma 2]. For synchronous workers, the bound adds terms in K §4, Theorem 3. The bounds count evaluations, not wall-clock time §4.

Paper vs code

The released code is one notebook (EGP_BO_ro_PAMI.ipynb) and a helper file (push_world.py), both at the commit above. The notebook implements the sequential EGP-TS, its baselines and the robot-pushing objective. The paper and the code disagree in the places below. Notebook line numbers are JSON lines.

• Refits are off by default. The paper refits hyperparameters every 50 iterations for EGP-TS §5. The notebook sets reinit_key = 0 [code: EGP_BO_ro_PAMI.ipynb#L217], and the refit block runs only when reinit_key == 1 EGP_BO_ro_PAMI.ipynb.
• Initial weights are uniform. Alg. 1 initializes the weights from the batch marginal likelihood alg. 1, eq. 9. The notebook sets all weights to 1/M and never updates them from the initial points [code: EGP_BO_ro_PAMI.ipynb#L553, #L567-L571].
• The first query repeats an initial point. The loop starts from the first initial point and evaluates it again, so that point enters the data twice [code: EGP_BO_ro_PAMI.ipynb#L554-L558, #L574].
• Refit drops the initial points. When reinit_key == 1, the refit uses X_list, which starts from that single repeated point, not the 10 initial pairs EGP_BO_ro_PAMI.ipynb. The paper says the refit uses all acquired data §3.1.
• Mixed-kernel refit resets the posterior. In the mixed-dictionary function, the refit block zeroes every component's mean and sets the covariance to the prior, but never conditions on any data again. The hyperparameters are fit to the initial points only EGP_BO_ro_PAMI.ipynb.
• Random-feature scale does not match the lengthscale. The RBF lengthscale is passed to scikit-learn as ker_var, and the same value is used to scale frequencies as 1/\sqrt{\ell} EGP_BO_ro_PAMI.ipynb. For an RBF kernel with lengthscale \ell the frequencies should have standard deviation 1/\ell. The ARD and Matérn generators repeat the 1/\sqrt{\cdot} scaling EGP_BO_ro_PAMI.ipynb. The paper gives no scaling rule beyond the spectral density §2.1, §8.
• Component lengthscales are fixed in EGP. For EGP components the notebook fixes each lengthscale to its dictionary value and fits only the amplitude EGP_BO_ro_PAMI.ipynb. The single-GP baselines fit both EGP_BO_ro_PAMI.ipynb. The paper describes the hyperparameters as fitted by marginal likelihood without this split §5.
• Noise prior differs. The fully Bayesian baseline uses noise \mathcal{U}[0.15, 0.25] in the code EGP_BO_ro_PAMI.ipynb and \mathcal{U}[0.1, 0.3] in the paper §5. The code also places a \mathcal{U}[0, 10^{-4}] prior on the constant mean, which the paper does not mention EGP_BO_ro_PAMI.ipynb.
• Robot 4-D ignores its fourth input. In the notebook, robot_pushing_4d overwrites its init_angle argument before use, so the push direction comes only from the first two inputs EGP_BO_ro_PAMI.ipynb. The paper says the fourth input is the push angle §5.2.
• One fixed goal, not ten. The notebook fixes the target location [code: EGP_BO_ro_PAMI.ipynb#L174, #L185]. The paper averages over ten random goals §5.2.
• Noise in the synthetic objective is broken. The notebook adds noise with y += random.normalvariate(0, tau) on a Python list [code: EGP_BO_ro_PAMI.ipynb#L178, #L189]. That line raises a TypeError as written, and if it ran it would use \tau = 0.2 as a standard deviation, while the GP uses it as a variance.
• The simulator is stochastic. Each push adds Gaussian noise with standard deviation 0.1 to the push velocity [code: push_world.py#L195, push_world.py#L210]. The objective is therefore noisy even when the tau noise is off.
• Prior-weight floor is off. The paper's floor of 10^{-4} appears only in commented-out text §5 comment. The notebook default is 0 [code: EGP_BO_ro_PAMI.ipynb#L54].
• Parts of the paper have no code here. The notebook has no code for the Ackley, Zakharov, Drop-wave or Eggholder benchmarks, no code for the FNN, SVM or gradient-boosting tuning, no data loaders for the four datasets, and no parallel or asynchronous loop. Results in §5.1, §5.3, §9.2 and alg. 2 cannot be reproduced from this repository.
• Weights update via density, not log density. The notebook multiplies the weights by the Gaussian density itself EGP_BO_ro_PAMI.ipynb. Matching the paper's eq. 13 in a numerically stable way needs logs.

Code: the repository is KonstantinosPol/EGP-BO, under MIT, with copyright held by Konstantinos D. Polyzos LICENSE. That is an author named on the paper, and the owner account's profile name is the same. The file push_world.py credits Ari Anders and Zi Wang in its header. The paper's footnote links a different repository (zi-w/Max-value-Entropy-Search), which we did not use. Search performed: the paper and the arXiv record have no code or project-page link. A GitHub user search for the author's name returned the owner account, whose repositories include EGP-BO, and that repository's licence names the author. The commit is the head of the default branch at clone time.

Pitfalls

• Equation numbers come from the .tex source. No PDF was compiled. Check equation numbers against the arXiv PDF before quoting them. Appendix equations are cited by section.
• Matérn frequencies are not specified. The paper does not say how to sample frequencies for the Matérn kernels. The code draws \mathbf{v} = \mathbf{z}/(\sqrt{\ell}\sqrt{g}) with g \sim \mathrm{Gamma}(\nu, 1/\nu) EGP_BO_ro_PAMI.ipynb, which carries the lengthscale mismatch above.
• Regret theory covers exact Gaussian processes. Theorems 1–3 bound the regret of the exact-kernel method. The source contains a commented-out remark that they ignore the random-feature approximation, and a bound for it is not in the paper §4 comment, §7. The random-feature method used in practice is not covered.
• Matérn smoothness is outside Assumption 1. The paper says Assumption 1 holds for Gaussian kernels and for Matérn kernels with \nu > 2 §4. The experiments use Matérn \nu = 3/2 and 5/2 §5, fig. 1.
• Constants are typeset with unbalanced parentheses. In the source, c_1 and c_3 read (2/\log(1+\sigma_n^{-2})^{1/2}. The intended form is (2/\log(1+\sigma_n^{-2}))^{1/2}, which matches the proof's step [§7.1, the bound on \mathcal{BR}_2].
• Proof cross-references are off. The main text points to "App. B" through "App. E" for the proofs and extra results. In this source, the proofs are §7.1–7.3, the random-feature section is §8, and the extra results are §9.
• Weight collapse. With the floor at 0, the mixture can put nearly all its weight on one component, and the others are then rarely sampled. The paper's floor appears only in commented-out text.
• Numerical risk in the rank-one update. The covariance update subtracts a rank-one term from a 2D \times 2D matrix with no Cholesky step or jitter EGP_BO_ro_PAMI.ipynb. For long runs the matrix can lose positive definiteness. The paper does not discuss this.
• Acquisition optimisation is inexact. Maximising the sampled feature function is non-convex, and the paper does not give starts or tolerances §3.1. The code uses 40 random starts with a box-constrained L-BFGS-B solver [code: EGP_BO_ro_PAMI.ipynb#L419-L430].
• Old pins in the setup. The notebook installs scipy==1.4.1 [code: EGP_BO_ro_PAMI.ipynb#L13], which may not run on current Python. It also needs Box2D and pygame for the robot simulator [code: EGP_BO_ro_PAMI.ipynb#L65-L66].
• Table 1 range looks odd. The gradient-boosting learning-rate range is printed as [0.1, 10] tab. 1. That is unusual for a learning rate; check the source before using it.
• Runtime units are missing. tab. 2 gives no unit, and neither does the text §9.3.
• Commented-out tables are not published. The source holds validation-accuracy tables in comments §5 comments. They are not in the compiled paper, so do not use their numbers.
• Figure-only results. The headline comparisons in §5.1–§5.4 are plots. Read them from the arXiv PDF, not from this card.

Lineage

• Builds on: the ensemble-of-GPs prior for supervised learning, with an online kernel selection method §1, §3; Thompson sampling for parallel Bayesian optimization §1, §2; random-feature approximation of shift-invariant kernels §2.1; Thompson sampling with a mixture prior in linear bandits, used for the regret proof §1, §4. No arXiv ids are given in the paper's bibliography, so none are linked here.
• Compared against: GP-EI §5; GP-TS with a fixed kernel and with ARD §5; fully Bayesian GP-TS §5; BanditBO and EXP3BO, adapted from categorical-input bandit methods §5.
• Superseded by: none in this wiki.