A Tutorial on Bayesian Optimization of Expensive Cost Functions, with Application to Active User Modeling and Hierarchical Reinforcement Learning

Eric BrochuVlad M. CoraNando de Freitas

article2010arXiv2,701 citations

Presents a comprehensive guide to Bayesian optimization for expensive black-box functions, demonstrating how exploration-exploitation trade-offs solve practical problems in active user modeling and hierarchical reinforcement learning.

Listen

The article addresses the challenge of optimizing expensive, black-box objective functions that lack closed-form expressions, derivatives, or convexity guarantees. Such problems arise in machine learning when evaluations involve costly simulations, human feedback, or physical experiments, making traditional optimization methods inefficient or impractical.

The work sets out to provide a tutorial on Bayesian optimization as an efficient solution for these settings, along with two practical extensions for active user modeling via preferences and hierarchical reinforcement learning, supported by experiments and a discussion of advantages and drawbacks.

The authors review the core Bayesian optimization framework, which places a Gaussian process prior over the unknown objective, updates it with observations to form a posterior, and selects the next evaluation point by maximizing an acquisition function that trades off exploration and exploitation. They cover covariance functions, acquisition functions such as expected improvement and upper confidence bounds, and noise handling. Extensions include a probit model for inferring latent value functions from pairwise preferences and integration of Gaussian process active learning with the MAXQ hierarchy for control tasks. Experiments test the preference approach on a material design gallery task and the hierarchical approach on a simulated Vancouver taxi navigation domain.

Bayesian optimization requires far fewer evaluations than alternatives by leveraging prior knowledge and balancing exploration against exploitation. In the preference gallery study, selecting points via expected improvement cut the average number of user clicks needed to identify a target material by more than half compared with random or uncertainty sampling. In the taxi domain, active Gaussian process value learning within the hierarchy produced effective policies after roughly four hours of simulated driving time, outperforming standard MAXQ and recursive average-reward variants. The framework handles noisy observations through adjusted kernels and incumbent definitions, and acquisition functions can be tuned with a small exploration parameter that often works well at a fixed value.

These results show that Bayesian optimization can make previously intractable optimization tasks feasible in human-in-the-loop and simulation-heavy domains, lowering evaluation costs and enabling better policies or designs with limited feedback. The approach is especially valuable when each query is expensive or when the user must remain engaged.

Practitioners should initialize with a modest set of random or Latin-hypercube samples, then maximize the acquisition function with a deterministic global optimizer such as DIRECT. Portfolios of acquisition functions governed by a bandit strategy often outperform any single choice. Further gains are possible by extending the method to multi-step lookahead, batch selection, and automatic feature relevance determination.

The main limitations are sensitivity to the choice of kernel and hyperparameters, the myopic single-step nature of standard acquisition functions, and degraded performance in high dimensions without additional structure. Readers should therefore validate kernel assumptions on representative data and consider dimensionality reduction when the input space is large.

bayesian methods genetic selection hierarchical reinforcement learning Bayesian Optimization Gaussian Process Acquisition Function

arXiv: 1012.2599

Paper: Finite-time Analysis of the Multiarmed Bandit Problem, Peter Auer et al. (2002). Provides the foundational finite-time regret analysis and upper confidence bound (UCB) strategies for exploration-exploitation trade-offs that underpin Gaussian process bandit acquisition functions.
Paper: Using Confidence Bounds for Exploitation-Exploration Trade-offs, P. Auer (2003). Establishes the statistical confidence-bound framework for balancing exploration and exploitation in continuous and complex decision spaces.
Paper: Active Learning with Statistical Models, D. Cohn et al. (1996). Introduces variance-minimizing active learning principles using statistical models, motivating the active query selection criteria utilized in Bayesian optimization.
Paper: Between MDPs and semi-MDPs: A framework for temporal abstraction in reinforcement learning, Richard S. Sutton et al. (1999). Defines the core mathematical framework for temporal abstraction and semi-Markov decision processes foundational to the MAXQ hierarchical reinforcement learning methods used in the source.
Paper: Reinforcement Learning: A Survey, Leslie Pack Kaelbling et al. (1996). Surveys the core principles, exploration-exploitation dilemmas, and value-function estimation in reinforcement learning that the tutorial builds on.
Paper: Algorithms for Inverse Reinforcement Learning, Andrew Y. Ng et al. (2000). Introduces the foundational formulation for recovering latent reward functions from preferences and observed behavior.

Paper: Practical Bayesian Optimization of Machine Learning Algorithms, Jasper Snoek et al. (2012). Extends Gaussian process Bayesian optimization to automated machine learning hyperparameter tuning by incorporating fully Bayesian hyperparameter treatment, timing-aware acquisitions, and parallel evaluation.
Paper: A Tutorial on Bayesian Optimization, Peter I. Frazier (2018). Provides an updated, comprehensive tutorial extending standard GP Bayesian optimization to modern acquisition strategies like knowledge gradients and entropy search.
Paper: Algorithms for Hyper-Parameter Optimization, James Bergstra et al. (2011). Builds on sequential model-based optimization by developing the Tree-structured Parzen Estimator as an alternative surrogate to standard Gaussian processes for high-dimensional hyperparameter spaces.
Paper: Making a Science of Model Search: Hyperparameter Optimization in Hundreds of Dimensions for Vision Architectures, J. Bergstra et al. (2013). Applies sequential model-based Bayesian optimization strategies to large-scale, high-dimensional computer vision architecture searches.
Paper: Hyperband: A Novel Bandit-Based Approach to Hyperparameter Optimization, Lisha Li et al. (2016). Proposes a bandit-based multi-fidelity alternative to Bayesian optimization, comparing resource allocation tradeoffs directly against GP-based optimizers.
Paper: Optuna: A Next-generation Hyperparameter Optimization Framework, Takuya Akiba et al. (2019). Presents a modern hyperparameter optimization framework implementing scalable tree-structured Bayesian optimization and automated trial pruning.
Paper: Deep reinforcement learning from human preferences, Paul F. Christiano et al. (2017). Scales the concept of learning latent reward functions from pairwise human preferences to deep reinforcement learning in continuous control and gaming tasks.
Paper: On Hyperparameter Optimization of Machine Learning Algorithms: Theory and Practice, Li Yang et al. (2020). Offers a broad survey and benchmark evaluation contrasting Bayesian optimization with alternative black-box hyperparameter tuning paradigms.

Cover for A Tutorial on Bayesian Optimization of Expensive Cost Functions, with Application to Active User Modeling and Hierarchical Reinforcement Learning

Abstract

We present a tutorial on Bayesian optimization, a method of finding the maximum of expensive cost functions. Bayesian optimization employs the Bayesian technique of setting a prior over the objective function and combining it with evidence to get a posterior function. This permits a utility-based selection of the next observation to make on the objective function, which must take into account both exploration (sampling from areas of high uncertainty) and exploitation (sampling areas likely to offer improvement over the current best observation). We also present two detailed extensions of Bayesian optimization, with experiments---active user modelling with preferences, and hierarchical reinforcement learning---and a discussion of the pros and cons of Bayesian optimization based on our experiences.

Topics

bayesian methods genetic selection hierarchical reinforcement learning optimization penalty function

Keywords

Bayesian Optimization Gaussian Process Acquisition Function Expected Improvement Preference Gallery Hierarchical Reinforcement Learning Active Learning Kernel Hyperparameter Noise Model MAXQ

1 Introduction
1.1 An Introduction to Bayesian Optimization
1.2 Overview
2 The Bayesian Optimization Approach
2.1 Priors over functions
2.2 Choice of covariance functions
2.3 Acquisition Functions for Bayesian Optimization
2.3.1 Improvement-based acquisition functions
2.3.2 Exploration-exploitation trade-off
2.3.3 Confidence bound criteria
2.3.4 Maximizing the acquisition function
2.4 Noise
2.5 A brief history of Bayesian optimization
2.6 Kriging
2.7 Experimental design
2.8 Active learning
2.9 Applications
3 Bayesian Optimization for Preference Galleries
3.1 Probit model for binary observations
3.2 Application: Interactive Bayesian optimization for material design
3.2.1 User Study
4 Bayesian Optimization for Hierarchical Control
4.1 Hierarchical Reinforcement Learning
4.1.1 Semi-MDPs
4.1.2 Hierarchical Value Function Decomposition
4.2 Application: The Vancouver Taxi Domain
4.2.1 State Abstraction, Termination and Rewards
4.3 Bayesian Optimization for Hierarchical Policies
4.3.1 Active Policy Optimization
4.3.2 Active Value Function Learning
4.4 Simulations
5 Discussion and advice to practitioners
References

Knowls

Knowl 1 — Sequential Bayesian Optimization Framework
algorithm
Bayesian optimization is a sequential strategy for finding the global maximum of an expensive, black-box objective function $f: \mathcal{A} \to \mathbb{R}$ over a compact search space $\mathcal{A} \subset \mathbb{R}^d$ . The approach combines a probabilistic prior (typically a Gaussian process) with accumulated observations $\mathcal{D}_{1:t-1} = \{(x_i, y_i)\}_{i=1}^{t-1}$ to construct a posterior distribution over functions. An acquisition function $u(x \mid \mathcal{D}_{1:t-1})$ , which quantifies the expected utility of querying at $x$ by balancing exploration of uncertain regions and exploitation of high-value regions, is maximized to choose the next query point $x_t$ .
Input: Black-box objective function $f$ , search domain $\mathcal{A} \subset \mathbb{R}^d$ , GP prior over $f$ , acquisition function $u(\cdot \mid \mathcal{D})$
Initialize dataset $\mathcal{D}_{1:0} = \emptyset$ (or evaluate an initial set of sample points)
for $t = 1, 2, \dots$ do
    Find next query point by maximizing acquisition: $x_t = \arg\max_{x \in \mathcal{A}} u(x \mid \mathcal{D}_{1:t-1})$
    Sample the objective function: $y_t = f(x_t) + \varepsilon_t$ , with $\varepsilon_t \sim \mathcal{N}(0, \sigma_{\text{noise}}^2)$
    Augment dataset: $\mathcal{D}_{1:t} = \mathcal{D}_{1:t-1} \cup \{(x_t, y_t)\}$
    Update Gaussian process posterior surrogate using $\mathcal{D}_{1:t}$
end for
The acquisition optimization step is solved using global derivative-free optimizers such as the DIRECT (Dividing Rectangles) algorithm, multi-start local optimizers, or Monte Carlo sampling.
Knowl 2 — Gaussian Process Predictive Distribution for Objective Functions
model/method
Given an unknown objective function $f(x)$ modeled by a zero-mean Gaussian process prior $f(x) \sim \mathcal{GP}(0, k(x, x'))$ and a set of $t$ observations $\mathcal{D}_{1:t} = \{(x_i, y_i)\}_{i=1}^t$ with additive Gaussian observation noise $y_i = f(x_i) + \varepsilon_i$ , where $\varepsilon_i \overset{\text{iid}}{\sim} \mathcal{N}(0, \sigma_{\text{noise}}^2)$ , the joint distribution of observed targets $\mathbf{y}_{1:t} = [y_1, \dots, y_t]^T$ and the latent value $f_{t+1} = f(x_{t+1})$ at a candidate point $x_{t+1} \in \mathcal{A}$ is Gaussian:

$\begin{bmatrix} \mathbf{y}_{1:t} \\ f_{t+1} \end{bmatrix} \sim \mathcal{N}\left( \mathbf{0}, \begin{bmatrix} \mathbf{K} + \sigma_{\text{noise}}^2 \mathbf{I} & \mathbf{k} \\ \mathbf{k}^T & k(x_{t+1}, x_{t+1}) \end{bmatrix} \right)$

where $\mathbf{K} \in \mathbb{R}^{t \times t}$ has entries $\mathbf{K}_{i,j} = k(x_i, x_j)$ , $\mathbf{I}$ is the $t \times t$ identity matrix, and $\mathbf{k} = [k(x_{t+1}, x_1), k(x_{t+1}, x_2), \dots, k(x_{t+1}, x_t)]^T \in \mathbb{R}^t$ .

The posterior predictive distribution is Gaussian $P(f_{t+1} \mid \mathcal{D}_{1:t}, x_{t+1}) = \mathcal{N}(\mu_t(x_{t+1}), \sigma_t^2(x_{t+1}))$ with sufficient statistics:

$\mu_t(x_{t+1}) = \mathbf{k}^T [\mathbf{K} + \sigma_{\text{noise}}^2 \mathbf{I}]^{-1} \mathbf{y}_{1:t}$

$\sigma_t^2(x_{t+1}) = k(x_{t+1}, x_{t+1}) - \mathbf{k}^T [\mathbf{K} + \sigma_{\text{noise}}^2 \mathbf{I}]^{-1} \mathbf{k}$

In noise-free settings ( $\sigma_{\text{noise}}^2 = 0$ ), observations equal the function values $\mathbf{f}_{1:t}$ , and the term $\mathbf{K} + \sigma_{\text{noise}}^2 \mathbf{I}$ reduces directly to $\mathbf{K}$ . In noisy settings, the incumbent used in acquisition criteria is defined as the evaluated location with the highest posterior mean $\mu^+ = \arg\max_{x_i \in \{x_1, \dots, x_t\}} \mu_t(x_i)$ to prevent optimizing against noisy overestimates.
Knowl 3 — Expected Improvement Acquisition Function with Exploration-Exploitation Parameter
equation
Let $f(x^+) = \max_{i=1,\dots,t} f(x_i)$ be the incumbent maximum observation in noise-free optimization (or $\mu^+$ in noisy settings). The improvement metric at a candidate query point $x$ is $I(x) = \max\{0, f_{t+1}(x) - f(x^+) - \xi\}$ , where $\xi \ge 0$ is a parameter controlling exploration-exploitation trade-off. Evaluating the expectation of $I(x)$ under the Gaussian process predictive distribution $\mathcal{N}(\mu(x), \sigma^2(x))$ yields the closed-form Expected Improvement (EI) function:

$\mathrm{EI}(x) = \begin{cases} (\mu(x) - f(x^+) - \xi) \Phi(Z) + \sigma(x) \phi(Z) & \text{if } \sigma(x) > 0 \\ 0 & \text{if } \sigma(x) = 0 \end{cases}$

where

$Z = \begin{cases} \dfrac{\mu(x) - f(x^+) - \xi}{\sigma(x)} & \text{if } \sigma(x) > 0 \\ 0 & \text{if } \sigma(x) = 0 \end{cases}$

$\Phi(\cdot)$ is the standard normal cumulative distribution function (CDF), and $\phi(\cdot)$ is the standard normal probability density function (PDF). Setting $\xi = 0.01$ (scaled by signal variance if necessary) provides robust empirical performance across diverse global optimization tasks, whereas cooling schedules for $\xi$ generally do not improve long-horizon optimization.
Knowl 4 — Preference-Based Bayesian Optimization via Probit Gaussian Process Regression
model/method
In interactive applications where absolute numerical evaluations of an objective function $f(x)$ are difficult or biased for human users, preferences over pairs of candidates are collected instead. Given $M$ observed pairwise comparisons $\mathcal{D} = \{r_i \succ c_i\}_{i=1}^M$ over $N$ distinct items $x_{1:N}$ , where $r_i \succ c_i$ denotes that item $r_i$ is preferred to $c_i$ , user valuations are modeled via a latent function $v(x) = f(x) + \varepsilon$ with independent Gaussian noise $\varepsilon \sim \mathcal{N}(0, \sigma_{\text{noise}}^2)$ and a GP prior $f \sim \mathcal{GP}(0, k(x, x'))$ .

The observation likelihood follows the Thurstone-Mosteller binomial-probit model:

$P(r_i \succ c_i \mid f(r_i), f(c_i)) = P(v(r_i) > v(c_i)) = \Phi(Z_i), \quad Z_i = \frac{f(r_i) - f(c_i)}{\sqrt{2}\sigma_{\text{noise}}}$

where $\Phi(\cdot)$ is the standard normal CDF. The posterior $P(\mathbf{f} \mid \mathcal{D}) \propto P(\mathbf{f}) \prod_{i=1}^M \Phi(Z_i)$ is approximated around the maximum a posteriori (MAP) estimate $\mathbf{f}_{\text{MAP}}$ using Laplace approximation: $P(\mathbf{f} \mid \mathcal{D}) \approx \mathcal{N}(\mathbf{K}\mathbf{b}, (\mathbf{K}^{-1} + \mathbf{C})^{-1})$ , where $\mathbf{f}_{\text{MAP}}$ is obtained via Newton-Raphson iterations $\mathbf{f}^{\text{new}} = \mathbf{f}^{\text{old}} - \mathbf{H}^{-1}\mathbf{g}$ with gradient $\mathbf{g} = -\mathbf{K}^{-1}\mathbf{f} + \mathbf{b}$ and Hessian $\mathbf{H} = \mathbf{K}^{-1} + \mathbf{C}$ .

The vector $\mathbf{b} \in \mathbb{R}^N$ and matrix $\mathbf{C} \in \mathbb{R}^{N \times N}$ have elements:

$b_j = \frac{1}{\sqrt{2}\sigma_{\text{noise}}} \sum_{i=1}^M \frac{\phi(Z_i)}{\Phi(Z_i)} h_i(x_j)$

$\mathbf{C}_{m,n} = \frac{1}{2\sigma_{\text{noise}}^2} \sum_{i=1}^M h_i(x_m) h_i(x_n) \left[ \frac{\phi(Z_i)}{\Phi^2(Z_i)} + \frac{\phi^2(Z_i)}{\Phi(Z_i)} Z_i \right]$

where $h_i(x_j) = 1$ if $x_j = r_i$ , $-1$ if $x_j = c_i$ , and $0$ otherwise, and $\phi(\cdot)$ is the standard normal PDF. Convolving with the GP conditional distribution yields the predictive distribution at candidate point $x_{t+1}$ :

$P(f_{t+1} \mid \mathcal{D}, x_{t+1}) \sim \mathcal{N}\left(\mathbf{k}^T \mathbf{K}^{-1} \mathbf{f}_{\text{MAP}}, \; k(x_{t+1}, x_{t+1}) - \mathbf{k}^T (\mathbf{K} + \mathbf{C}^{-1})^{-1} \mathbf{k}\right)$
Knowl 5 — Gaussian Process Upper Confidence Bound Acquisition Function with Regret Guarantees
theoretical result
When Bayesian optimization is formulated in the multi-armed bandit framework over a compact search space $\mathcal{A} \subset \mathbb{R}^d$ , the instantaneous regret is $r(x_t) = f(x^*) - f(x_t)$ , where $x^* = \arg\max_{x \in \mathcal{A}} f(x)$ . The Gaussian Process Upper Confidence Bound (GP-UCB) acquisition function is defined as:

$\mathrm{GP\text{-}UCB}(x) = \mu(x) + \sqrt{\nu \tau_t} \sigma(x)$

where $\mu(x)$ is the GP posterior mean, $\sigma(x)$ is the GP posterior standard deviation, $\nu > 0$ is a tuning parameter, and $\tau_t$ is an iteration-dependent exploration schedule.

Under regular smoothness assumptions on the GP covariance kernel, setting $\nu = 1$ and $\tau_t = 2 \log\left( \frac{t^{d/2 + 2}\pi^2}{3\delta} \right)$ guarantees that the algorithm is no-regret with probability at least $1 - \delta$ for $\delta \in (0, 1)$ , satisfying $\lim_{T \to \infty} \frac{R_T}{T} = 0$ , where cumulative regret is $R_T = \sum_{t=1}^T [f(x^*) - f(x_t)]$ . This provides a theoretical lower bound on the convergence rate to the global maximum.
Knowl 6 — Probability of Improvement Acquisition Function
equation
The Probability of Improvement (PI) acquisition function, also known as Maximum Probability of Improvement (MPI) or Kushner's P-algorithm, selects the candidate point $x \in \mathcal{A}$ that maximizes the probability that the objective function value improves by at least a target margin $\xi \ge 0$ over the current incumbent value $f(x^+)$ , where $x^+ = \arg\max_{x_i \in \{x_1, \dots, x_t\}} f(x_i)$ :

$\mathrm{PI}(x) = P(f(x) \ge f(x^+) + \xi) = \Phi\left( \frac{\mu(x) - f(x^+) - \xi}{\sigma(x)} \right)$

where $\mu(x)$ and $\sigma(x)$ are the posterior mean and standard deviation predicted at $x$ by the Gaussian process surrogate, and $\Phi(\cdot)$ denotes the cumulative distribution function of the standard normal distribution.

When $\xi = 0$ , the criterion is purely exploitative and tends to query points with high certainty of infinitesimal improvements over points offering potentially larger gains with greater uncertainty. When $\xi$ is chosen too large, the search becomes excessively global and fails to refine promising local optima.

Knowl 7 — User Study Comparison of Sampling Strategies in BRDF Preference Gallery

data/table

A user study with 5 human subjects was conducted on an interactive material design task using the MERL database (restricted to a representative subset of 38 measured materials). The objective was to locate a randomly chosen target Bidirectional Reflectance Distribution Function (BRDF) material through iterative pairwise comparisons. In each iteration, users were shown a pair containing the incumbent $\arg\max_{x_i} \mu(x_i)$ alongside a candidate point selected by one of three strategies: uniform random sampling, maximum uncertainty ( $\arg\max_x \sigma(x)$ ), or maximum Expected Improvement ( $\arg\max_x \mathrm{EI}(x)$ ).

The results demonstrate that Expected Improvement required less than half the user interactions (8.56 clicks on average) compared to random sampling (18.40) and maximum uncertainty (17.87). Selecting queries based purely on maximum predictive variance ( $\arg\max_\sigma$ ) did not significantly improve over random selection because it preferentially selected boundary points of the parameter domain, providing low utility for identifying the target.

Knowl 8 — Active Path Value Function Learning Algorithm with Gaussian Process Exploration
algorithm
Active Path Learning integrates Gaussian process regression and Expected Improvement into a MAXQ task hierarchy to learn the value function over continuous map coordinates for a navigation subtask, minimizing the number of world samples needed to find optimal multi-step paths.
function NavigateTaskLearner(NavigateTask $i$ , State $s$ )
    Let $trajectory = ()$
    Let $intersections = ()$
    Let $visits = 0$
    while $Terminated_i(s)$ is false do
        Choose adjacent waypoint $WP$ via active GP Expected Improvement (or $\epsilon$ -greedy)
         $childSeq = Follow(WP, s)$
        Append $childSeq$ to $trajectory$
        Observe resulting state $s'$
         $N_{steps} = \text{length}(childSeq)$
         $R = \sum_{j=1}^{N_{steps}} \gamma^{N_{steps}-j} r_j$
         $V'_s = V(TaxiLoc_{s'}, Target_i)$
         $V_s = V(TaxiLoc_s, Target_i)$
        if $Terminated_i(s')$ is true then
             $V_s \leftarrow (1 - \alpha) V_s + \alpha R$
            for $j = 1$ to $\text{length}(intersections)$ do
                 $\{s_{prev}, N_{prev}, R_{prev}\} = intersections[j]$
                 $R \leftarrow R_{prev} + \gamma^{N_{prev}} R$
                 $V'_{prev} \leftarrow V(TaxiLoc_{s_{prev}}, Target_i)$
                 $V'_{prev} \leftarrow (1 - \alpha) V'_{prev} + \alpha R$
            end for
        else
            Append $\{s, N_{steps}, R\}$ to $intersections$
             $visits(TaxiLoc_s) \leftarrow visits(TaxiLoc_s) + 1$
             $penalty \leftarrow V_s \times visits(TaxiLoc_s)$
             $V_s \leftarrow (1 - \alpha) V_s + \alpha (penalty + R + \gamma^{N_{steps}} V'_s)$
        end if
         $s \leftarrow s'$
    end while
    return $trajectory$
end function
The 4D value function $V(x_C, y_C, x_T, y_T)$ is parameterized by continuous 2D coordinates of the current intersection and target destination. The GP response surface interpolates values across spatial coordinates, while the visit penalty prevents infinite looping during exploration before the destination is reached.
Knowl 9 — Active Policy Optimization with Gaussian Processes in Hierarchical RL
algorithm
In hierarchical reinforcement learning, low-level continuous controllers (such as steering and throttle trajectory tracking) require finding a vector of policy parameters $x \in \mathbb{R}^d$ that maximizes the trajectory return $V(x)$ . Because trajectory evaluations through full simulation episodes are computationally expensive, Gaussian process Bayesian optimization is applied directly to the policy parameter space.
Input: Trajectory policy evaluator $V(\cdot)$ , parameter dimension $d = 15$ , initial sample size $N_0 = 30$
Generate initial parameter configurations $\{x_1, \dots, x_{N_0}\}$ using Latin hypercube sampling
Evaluate initial configurations: $V_i = V(x_i)$ for $i = 1, \dots, N_0$
Initialize dataset $\mathcal{D}_{1:N_0} = \{(x_i, V_i)\}_{i=1}^{N_0}$
$N \leftarrow N_0$
while stopping criterion is not met do
    Fit GP surrogate on $\mathcal{D}_{1:N}$
    Update Expected Improvement function $\mathrm{EI}(x)$ over $\mathcal{D}_{1:N}$
    Select candidate policy parameters: $x_{N+1} = \arg\max_x \mathrm{EI}(x)$
    Evaluate candidate: $V_{N+1} = V(x_{N+1})$ via simulation rollouts
    Augment dataset: $\mathcal{D}_{1:N+1} = \mathcal{D}_{1:N} \cup \{(x_{N+1}, V_{N+1})\}$
     $N \leftarrow N + 1$
end while
return Optimal policy parameters $x^* = \arg\max_{x_i} V_i$
In vehicle trajectory-following experiments in the TORCS simulator with a 15-parameter neural network controller, the policy evaluation averages returns across 10 distinct starting angles along a trajectory requiring u-turns, acceleration, and deceleration. Using Expected Improvement, a high-performance driving policy is acquired in approximately 50 parameter samples beyond the initial 30 Latin hypercube samples.
Knowl 10 — Empirical Convergence Speedup of Active Value Learning in the Vancouver Taxi Domain
empirical result
In the Vancouver Taxi domain simulated in TORCS (comprising 61 road intersection waypoints, 22 potential pickup/drop-off destinations, and roughly 28 km of navigable roads with full 3D vehicle physics), hierarchical learners were evaluated across 10 independent experimental runs. The methods compared were:
1. Standard Recursive Average Reward (RAR),
2. Standard MAXQ with discounted rewards,
3. MAXQ with $\epsilon$ -greedy path exploration ( $\epsilon = 0.1$ ),
4. MAXQ with Gaussian Process Active Value Learning ( $V_{TM}$ GP with kernel bandwidth $k = 0.01$ and annealing exploration parameter $\epsilon = 0.2$ ).
The MAXQ $V_{TM}$ GP learner converged to an accumulated reward of over 2000 per episode within approximately $4 \times 10^4$ reward samples (corresponding to roughly 4 hours of real-time driving at 60 km/h). In contrast, standard MAXQ and RAR required between $4 \times 10^5$ and $5 \times 10^5$ reward samples to reach similar accumulated reward levels, demonstrating an order-of-magnitude reduction in sample complexity due to GP-based spatial value interpolation and active exploration.
Knowl 11 — Practical Limitations and Scalability Bottlenecks of Bayesian Optimization
limitation
The practical application of Bayesian optimization is constrained by several fundamental limitations:
1. Prior and Kernel Sensitivity: Optimization performance depends strongly on the GP prior and kernel choice (e.g., length-scale hyperparameters). Without sufficient initial data or accurate domain knowledge, poor kernel specifications lead to excessive unproductive exploration or premature entrapment in local extrema.
2. Curse of Dimensionality: As input dimensionality $d$ increases, sample requirements to cover the space grow exponentially, hyperparameter estimation becomes difficult, and optimizing the acquisition function itself becomes an expensive multi-modal global optimization problem.
3. Myopic Horizon: Standard acquisition functions (PI, EI, GP-UCB) are myopic, selecting only the immediate one-step-ahead sample ( $t+1$ ) rather than planning across a multi-step observation budget.
4. Sequential Querying Constraint: Standard formulations select a single point per iteration, limiting efficiency in batched or parallel experimental setups unless explicit batch extensions are used.
5. Stationarity and Noise Assumptions: Typical formulations rely on stationary covariance structures and homoskedastic Gaussian noise, requiring complex extensions to handle time-varying dynamics, heteroskedasticity, or non-Gaussian observations.

Coverage note — Detailed historical literature reviews regarding classical geostatistical kriging formulas, Wiener process global optimization, and DACE experimental design taxonomy were deliberately omitted as they represent background and contextual survey material rather than the primary methodological, theoretical, and experimental contributions.

References

1.D. Andre. Programmable Reinforcement Learning Agents. PhD thesis, University of California at Berkley, 2003.
2.C. Audet, J. Jr, Dennis, D. W. Moore, A. Booker, and P. D. Frank. Surrogate-model-based method for constrained optimization. In AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, 2000.
3.J. Azimi, A. Fern, and X. Z. Fern. Batch Bayesian optimization via simulation matching. In Advances in Neural Information Processing Systems 24, 2011.
4.A. G. Barto and S. Mahadevan. Recent advances in hierarchical reinforcement learning. Discrete Event Dynamic Systems, 13(1-2):41–77, 2003.
5.T. Bartz-Beielstein, C. Lasarczyk, and M. Preuss. Sequential parameter optimization. In Proc. CEC-05, 2005.
6.J. Baxter and P. L. Bartlett. Infinite-horizon policy-gradient estimation. Journal of Artificial Intelligence Research, 15:319–350, 2001.
7.D. P. Bertsekas and J. N. Tsitsiklis. Neuro-Dynamic Programming. Athena Scientific, 1996.
8.B. Betr`o. Bayesian methods in global optimization. J. Global Optimization, 1:1–14, 1991.
9.P. Boyle. Gaussian Processes for Regression and Optimisation. PhD thesis, Victoria University of Wellington, Wellington, New Zealand, 2007.
10.E. Brochu, N. de Freitas, and A. Ghosh. Active preference learning with discrete choice data. In Advances in Neural Information Processing Systems 21, 2007.
11.E. Brochu, A. Ghosh, and N. de Freitas. Preference galleries for material design. In ACM SIGGRAPH 2007 Posters, page 105, 2007.
12.E. Brochu, T. Brochu, and N. de Freitas. A Bayesian interactive optimization approach to procedural animation design. In Eurographics/ ACM SIGGRAPH Symposium on Computer Animation, 2010.
13.E. Brochu, M. Hoffman, and N. de Freitas. Hedging strategies for Bayesian optimization. eprint arXiv:1009.5419, arXiv.org, September 2010.
14.D. Busby. Hierarchical adaptive experimental design for Gaussian process emulators. Reliability Engineering and System Safety, 94(7):1183–1193, July 2009.
15.W. Chu and Z. Ghahramani. Extensions of Gaussian processes for ranking: semi-supervised and active learning. In Learning to Rank workshop at NIPS-18, 2005.
16.W. Chu and Z. Ghahramani. Preference learning with Gaussian processes. In Proc. 22nd International Conf. on Machine Learning, 2005.
17.V. M. Cora. Model-based active learning in hierarchical policies. Master's thesis, University of British Columbia, Vancouver, Canada, April 2008.
18.D. D. Cox and S. John. A statistical method for global optimization. In Proc. IEEE Conference on Systems, Man and Cybernetics, volume 2, pages 1241–1246, 1992.
19.D. D. Cox and S. John. SDO: A statistical method for global optimization. In M. N. Alexandrov and M. Y. Hussaini, editors, Multidisciplinary Design Optimization: State of the Art, pages 315–329. SIAM, 1997.
20.T. G. Dietterich. Hierarchical reinforcement learning with the maxq value function decomposition. Journal of Artificial Intelligence Research, 13:227–303, 2000.
21.P. J. Diggle and P. J. Ribeiro. Model-Based Geostatistics. Springer Series in Statistics. Springer, 2007.
22.P. J. Diggle, J. A. Tawn, and R. A. Moyeed. Model-based geostatistics. Journal of the Royal Statistical Society: Series C (Applied Statistics), 47(3):299–350, 1998.
23.J. F. Elder, IV. Global $\mathbb{R}^d$ optimization when probes are expensive: The GROPE algorithm. In Proc. IEEE International Conference on Systems, Man and Cybernetics, 1992.
24.Á. Élő. The Rating of Chess Players: Past and Present. Arco Publishing, New York, 1978.
25.M. Frean and P. Boyle. Using Gaussian processes to optimize expensive functions. In W. Wobcke and M. Zhang, editors, AI 2008: Advances in Artificial Intelligence, volume 5360 of Lecture Notes in Computer Science, pages 258–267. Springer Berlin / Heidelberg, 2008.
26.R. Garnett, M. Osborne, S. Reece, A. Rogers, and S. Roberts. Sequential Bayesian prediction in the presence of changepoints and faults. The Computer Journal, 2010.
27.R. Garnett, M. Osborne, and S. Roberts. Bayesian optimization for sensor set selection. In Proceedings of the 9th ACM/IEEE International Conference on Information Processing in Sensor Networks, pages 209–219. ACM, 2010.
28.M. G. Genton. Classes of kernels for machine learning: A statistics perspective. Journal of Machine Learning Research, 2:299–312, 2001.
29.M. Ghavamzadeh. Hierarchical Reinforcement Learning in Continuous State and Multi-agent Environments. PhD thesis, University of Massachusetts Amherst, 2005.
30.D. Ginsbourger, R. Le Riche, and L. Carraro. A Multipoints Criterion for Deterministic Parallel Global Optimization based on Gaussian Processes. 2008.
31.P. W. Goldberg, C. K. I. Williamsn, and C. M. Bishop. Regression with input-dependent noise: A Gaussian process treatment. In Advances in Neural Information Processing Systems 10, 1998.
32.R. Herbrich and T. Graepel. Trueskill: A Bayesian skill rating system. Technical Report MSR-TR-2006-80, Microsoft Research, June 2006.
33.G. Hinton and R. Salakhutdinov. Reducing the dimensionality of data with neural networks. Science, 313(5786):504 – 507, 2006.
34.C. Holmes and L. Held. Bayesian auxiliary variable models for binary and multinomial regression. Bayesian Analysis, 1(1):145–168, 2006.
35.D. Huang, T. T. Allen, W. I. Notz, and N. Zheng. Global optimization of stochastic black-box systems via sequential Kriging meta-models. J. Global Optimization, 34(3):441–466, March 2006.
36.F. Hutter, H. H. Hoos, K. Leyton-Brown, and K. P. Murphy. An experimental investigation of model-based parameter optimisation: SPO and beyond. In Proc. GECCO’09, 2009.
37.F. Hutter. Automating the Configuration of Algorithms for Solving Hard Computational Problems. PhD thesis, University of British Columbia, Vancouver, Canada, August 2009.
38.D. R. Jones, C. D. Perttunen, and B. E. Stuckman. Lipschitzian optimization without the Lipschitz constant. J. Optimization Theory and Apps, 79(1):157–181, 1993.
39.D. R. Jones, M. Schonlau, and W. J. Welch. Efficient global optimization of expensive black-box functions. J. Global Optimization, 13(4):455–492, 1998.
40.D. R. Jones. A taxonomy of global optimization methods based on response surfaces. J. Global Optimization, 21:345–383, 2001.
41.D. Kahneman and A. Tversky. Prospect theory: an analysis of decision making under risk. Econometrica, 47:263–291, 1979.
42.A. Kapoor, K. Grauman, R. Urtasun, and T. Sarrell. Active learning with Gaussian processes for object categorization. In Proc. International Conference on Computer Vision (ICCV), 2007.
43.M. Kendall. Rank Correlation Methods. Griffin Ltd, 1975.
44.D. C. Kingsley. Preference uncertainty, preference refinement and paired comparison choice experiments. Working Paper 06-06, Center for Economic Analysis, University of Colorado at Boulder, 2006. Dept. of Economics, University of Colorado.
45.A. Krause, A. Singh, and C. Guestrin. Near-optimal sensor placements in Gaussian processes: Theory, efficient algorithms and empirical studies. J. Machine Learning Research, 9:235–284, 2008.
46.D. G. Krige. A statistical approach to some basic mine valuation problems on the Witwatersrand. J. the Chemical, Metallurgical and Mining Soc. of South Africa, 52(6), 1951.
47.H. J. Kushner and G. G. Yin. Stochastic Approximation Algorithms and Applications. Springer-Verlag, 1997.
48.H. J. Kushner. A new method of locating the maximum of an arbitrary multipeak curve in the presence of noise. J. Basic Engineering, 86:97–106, 1964.
49.D. Lewis and W. Gale. A sequential algorithm for training text classifiers. In Proc. ACM SIGIR Conference on Research and Development in Information Retreival, 1994.
50.L. Liberti and N. Maculan, editors. Global Optimization: From Theory to Implementation. Springer Nonconvex Optimization and Its Applications. Springer, 2006.
51.D. Lizotte, T. Wang, M. Bowling, and D. Schuurmans. Automatic gait optimization with Gaussian process regression. In IJCAI, 2007.
52.D. Lizotte. Practical Bayesian Optimization. PhD thesis, University of Alberta, Edmonton, Alberta, Canada, 2008.
53.M. Locatelli. Bayesian algorithms for one-dimensional global optimization. J. Global Optimization, 1997.
54.D. J. C. Mackay. A practical bayesian framework for backpropagation networks. Neural Computation, 4(3):448–472, 1992.
55.B. Marthi, D. Latham, S. Russell, and C. Guestrin. Concurrent hierarchical reinforcement learning. In Proceedings of the 19th International Joint Conference on Artificial Intelligence, 2005.
56.R. Martinez–Cantin, N. de Freitas, and J. Castellanos. Analysis of particle methods for simultaneous robot localization and mapping and a new algorithm: Marginal-SLAM. In Proc. IEEE International Conference on Robots and Automation, 2006.
57.R. Martinez–Cantin, N. de Freitas, A. Doucet, and J. A. Castellanos. Active policy learning for robot planning and exploration under uncertainty. Robotics: Science and Systems (RSS), 2007.
58.R. Martinez–Cantin, N. de Freitas, E. Brochu, J. Castellanos, and A. Doucet. A Bayesian exploration-exploitation approach for optimal online sensing and planning with a visually guided mobile robot. Autonomous Robots, 27(2):93–103, 2009.
59.B. Matérn. Spatial Variation. Springer-Verlag, 2nd (1986) edition, 1960.
60.G. Matheron. The theory of regionalized variables and its applications. Cahier du Centre de Morphologie Mathematique, Ecoles des Mines, 1971.
61.D. McFadden. Econometric models for probabilistic choice among products. Journal of Business, 53(3):13–29, 1980.
62.D. McFadden. Economic choices. The American Economic Review, 91:351–378, 2001.
63.J. Močkus, V. Tiesis, and A. Žilinskas. Toward Global Optimization, volume 2, chapter The Application of Bayesian Methods for Seeking the Extremum, pages 117–128. Elsevier, 1978.
64.J. Močkus. The Bayesian approach to global optimization. In R. Drenick and F. Kozin, editors, System Modeling and Optimization, volume 38, pages 473–481. Springer Berlin / Heidelberg, 1982.
65.J. Močkus. Application of Bayesian approach to numerical methods of global and stochastic optimization. J. Global Optimization, 4(4):347 – 365, 1994.
66.M. Mongeau, H. Karsenty, V. Rouzé, and J.-B. Hiriart-Urruty. Comparison of public-domain software for black-box global optimization. Technical Report LAO 98-01, Universite Paul Sabatier, Toulouse, France, 1998.
67.F. Mosteller. Remarks on the method of paired comparisons: I. the least squares solution assuming equal standard deviations and equal correlations. Psychometrika, 16:3–9, 1951.
68.R. Murray-Smith and A. Girard. Gaussian process priors with ARMA noise models. In Irish Signals and Systems Conference, 2001.
69.R. Murray-Smith and D. Sbarbaro. Nonlinear adaptive control using non-parametric Gaussian process prior models. In 15th IFAC World Congress on Automatic Control. Citeseer, 2002.
70.A. Y. Ng and M. I. Jordan. Pegasus: A policy search method for large MDPs and POMDPs. In Uncertainty in Artificial Intelligence (UAI2000), 2000.
71.A. O’Hagan. On curve fitting and optimal design for regression. Journal of the Royal Statistical Society B, 40:1–42, 1978.
72.M. Osborne, R. Garnett, and S. Roberts. Active data selection for sensor networks with faults and changepoints. In IEEE International Conference on Advanced Information Networking and Applications, 2010.
73.M. Osborne. Bayesian Gaussian Processes for Sequential Prediction, Optimization and Quadrature. PhD thesis, University of Oxford, 2010.
74.R. E. Parr. Hierarchical control and learning for markov decision processes. PhD thesis, 1998. Chair-Stuart Russell.
75.J. W. Payne, J. R. Bettman, and E. J. Johnson. The Adaptive Decision Maker. Cambridge University Press, 1993.
76.G. Poyiadjis, A. Doucet, and S. S. Singh. Particle methods for optimal filter derivative: Application to parameter estimation. In IEEE ICASSP, pages 925–928, 2005.
77.W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipes: The Art of Scientific Computing. Cambridge University Press, 3rd edition, 2007.
78.C. E. Rasmussen and C. K. I. Williams. Gaussian Processes for Machine Learning. MIT Press, Cambridge, Massachusetts, 2006.
79.J. Sacks, W. J. Welch, T. J. Welch, and H. P. Wynn. Design and analysis of computer experiments. Statistical Science, 4(4):409–423, 1989.
80.T. J. Santner, B. Williams, and W. Notz. The Design and Analysis of Computer Experiments. Springer, 2003.
81.M. J. Sasena. Flexibility and Efficiency Enhancement for Constrained Global Design Optimization with Kriging Approximations. PhD thesis, University of Michigan, 2002.
82.M. Schonlau. Computer Experiments and Global Optimization. PhD thesis, University of Waterloo, Waterloo, Ontario, Canada, 1997.
83.B. Settles. Active learning literature survey. Computer Science Technical Report 1648, University of Wisconsin-Madison, January 2010.
84.S. Siegel and N. J. Castellan. Nonparametric Statistics for the Behavioral Sciences. McGraw-Hill, 1988.
85.N. Srinivas, A. Krause, S. M. Kakade, and M. Seeger. Gaussian process optimization in the bandit setting: No regret and experimental design. In Proc. International Conference on Machine Learning (ICML), 2010.
86.M. L. Stein. Interpolation of Spatial Data: Some Theory for Kriging. Springer Series in Statistics. Springer, 1999.
87.H. Stern. A continuum of paired comparison models. Biometrika, 77:265–273, 1990.
88.S. Streltsov and P. Vakili. A non-myopic utility function for statistical global optimization algorithms. J. Global Optimization, 14:283–298, 1999.
89.B. Stuckman. A global search method for optimizing nonlinear systems. IEEE Transactions on Systems, Man and Cybernetics, 18(6):965–977, 1988.
90.R. S. Sutton and A. G. Barto. Reinforcement Learning: An Introduction. MIT Press, 1998.
91.R. S. Sutton, D. Precup, and S. P. Singh. Between MDPs and semi-MDPs: A framework for temporal abstraction in reinforcement learning. Artificial Intelligence, 112(1-2):181–211, 1999.
92.L. Thurstone. A law of comparative judgement. Psychological Review, 34:273–286, 1927.
93.A. Törn and A. Žilinskas. Global Optimization. Springer-Verlag, 1989.
94.A. Tversky and D. Kahneman. Advances in prospect theory: Cumulative representation of uncertainty. J. Risk and Uncertainty, 5:297–323, 1992.
95.E. Vasquez and J. Bect. On the convergence of the expected improvement algorithm. Technical Report arXiv:0712.3744v2, arXiv.org, Feb 2008.
96.B. J. Williams, T. J. Santner, and W. I. Notz. Sequential design of computer experiments to minimize integrated response functions. Statistica Sinica, 10:1133–1152, 2000.
97.B. Wymann, C. Dimitrakakis, and C. Alexopoulos. The open racing car simulator (http://torcs.sourceforge.net/), 2009.
98.L. Younes. Parameter estimation for imperfectly observed Gibbsian fields. Prob. Theory and Rel. fields, 82:625–645, 1989.
99.A. Zhigljavsky and A. Žilinskas. Stochastic Global Optimization. Springer Optimization and Its Applications. Springer, 2008.
100.A. Žilinskas and J. Žilinskas. Global optimization based on a statistical model and simplical partitioning. Computers and Mathematics with Applications, 44:957–967, 2002.
101.A. Žilinskas. Lecture Notes in Control and Information Sciences, chapter On the Use of Statistical Models of Multimodal Functions for the Construction of Optimization Algorithms. Number 23. Springer-Verlag, 1980.

Citation

MLA

Brochu, E., et al. “A Tutorial on Bayesian Optimization of Expensive Cost Functions, with Application to Active User Modeling and Hierarchical Reinforcement Learning”. arXiv, 2010, https://doi.org/10.48550/arxiv.1012.2599.

APA

Brochu, E., Cora, V. M., & de Freitas, N. (2010). A Tutorial on Bayesian Optimization of Expensive Cost Functions, with Application to Active User Modeling and Hierarchical Reinforcement Learning. arXiv. https://doi.org/10.48550/arxiv.1012.2599

Chicago

Brochu, E., V. M. Cora, and N. de Freitas. 2010. “A Tutorial on Bayesian Optimization of Expensive Cost Functions, with Application to Active User Modeling and Hierarchical Reinforcement Learning”. Preprint, ArXiv. https://doi.org/10.48550/arxiv.1012.2599.

Harvard

Brochu, E., Cora, V.M. and de Freitas, N. (2010) “A Tutorial on Bayesian Optimization of Expensive Cost Functions, with Application to Active User Modeling and Hierarchical Reinforcement Learning”. arXiv. Available at: https://doi.org/10.48550/arxiv.1012.2599.

Vancouver

1. Brochu E, Cora VM, de Freitas N (2010) A Tutorial on Bayesian Optimization of Expensive Cost Functions, with Application to Active User Modeling and Hierarchical Reinforcement Learning. https://doi.org/10.48550/arxiv.1012.2599

BibTeX

@misc{https://doi.org/10.48550/arxiv.1012.2599,
  doi = {10.48550/ARXIV.1012.2599},
  url = {https://arxiv.org/abs/1012.2599},
  author = {Brochu, Eric and Cora, Vlad M. and de Freitas, Nando},
  keywords = {Machine Learning (cs.LG), FOS: Computer and information sciences, FOS: Computer and information sciences, G.1.6; G.3; I.2.6},
  title = {A Tutorial on Bayesian Optimization of Expensive Cost Functions, with Application to Active User Modeling and Hierarchical Reinforcement Learning},
  publisher = {arXiv},
  year = {2010},
  copyright = {arXiv.org perpetual, non-exclusive license}
}

Metadata:DOI registry

Access the Paper

This paper is available from its original source. Click below to access the PDF.

Open PDF

License: Authors

Algorithm	Trials	Number of Selections $n$ (mean $\pm$ std)
random	50	18.40 $\pm$ 7.87
$\arg\max_\sigma$	50	17.87 $\pm$ 8.60
$\arg\max_{\mathrm{EI}}$	50	8.56 $\pm$ 5.23

A Tutorial on Bayesian Optimization of Expensive Cost Functions, with Application to Active User Modeling and Hierarchical Reinforcement Learning

Abstract

Topics

Keywords

Table of Contents

Knowls

Knowl 1 — Sequential Bayesian Optimization Framework

Knowl 2 — Gaussian Process Predictive Distribution for Objective Functions

Knowl 3 — Expected Improvement Acquisition Function with Exploration-Exploitation Parameter

Knowl 4 — Preference-Based Bayesian Optimization via Probit Gaussian Process Regression

Knowl 5 — Gaussian Process Upper Confidence Bound Acquisition Function with Regret Guarantees

Knowl 6 — Probability of Improvement Acquisition Function

Knowl 7 — User Study Comparison of Sampling Strategies in BRDF Preference Gallery

Knowl 8 — Active Path Value Function Learning Algorithm with Gaussian Process Exploration

Knowl 9 — Active Policy Optimization with Gaussian Processes in Hierarchical RL

Knowl 10 — Empirical Convergence Speedup of Active Value Learning in the Vancouver Taxi Domain

Knowl 11 — Practical Limitations and Scalability Bottlenecks of Bayesian Optimization

References

Citation

BibTeX

Access the Paper