Bayesian Optimization for Policy Search via Online-Offline Experimentation

Benjamin Letham, Eytan Bakshy — Facebook AI Research, JMLR 2019 · arXiv:1904.01049

Relations

Builds on: Swersky et al. 2013 (Multi-task Bayesian Optimization, NIPS), Chai 2009 (MTGP generalization bound, NIPS), Letham et al. 2019 (Constrained BO with noisy EI, Bayesian Analysis) Concepts used: MTBO Prerequisites Curriculum (GP regression, ICM kernels, EI acquisition)

Table of Contents

• 1. The Value Model as a Policy Space
• 2. The Offline Simulator and Its Biases
  • 2.1 How the Naive Simulator Works
  • 2.2 Why Simulators Systematically Overestimate
• 3. Multi-Task Gaussian Process with ICM Kernel
  • 3.1 From Single-Task to Multi-Task GP
  • 3.2 The ICM Kernel
  • 3.3 Bias as a GP: The Key Structural Insight
  • 3.4 Hyperparameter Inference
• 4. The Bayesian Optimization Loop
  • 4.1 Noisy EI under Constraints
  • 4.2 Multi-Task Acquisition and the Interleaved Loop
  • 4.3 Simulator Non-Stationarity
• 5. Reading the Empirical Results
  • 5.1 LOO-CV Scatter Plots (Figs. 2–3)
  • 5.2 Learning Curves (Fig. 7)
  • 5.3 The Kernel Inference Ablation (Fig. 10)
• 6. The ρ² Bound on Multi-Task Generalization
  • 6.1 Chai (2009) Bound
  • 6.2 Bound vs. Empirical (Fig. 9)
• References

1. The Value Model as a Policy Space

A modern industrial recommender system (e.g., Facebook News Feed) relies on a cascade of ML models, each predicting a specific event probability: probability of click, share, comment, long view, etc. The value model is the function that combines these predictions into a single scalar score used to rank items:

\[s(i) = \sum_j w_j \cdot \hat{p}_j(i) + \text{conditional rules}(\mathbf{x}_i)\]

Here \(\hat{p}_j(i)\) is the event model prediction for item \(i\) on event type \(j\), and the weights \(w_j\) plus any conditional rules (e.g., “apply weight \(\alpha\) to videos with engagement rate \(> \theta\)”) constitute the policy \(\pi \in \mathbb{R}^d\).

The key consequence: outcomes of interest are nonlinear functions of the policy \(\pi\). Changing \(w_\text{comment}\) by 0.1 might increase comments by 2% but decrease long-views by 1%, with the trade-off depending on the current operating point. This nonlinear response surface is what Bayesian optimization is designed to navigate.

📐 Critically, the response surface can only be observed online — event models are trained on historical data and optimized for prediction, not for evaluating counterfactual policy changes. Only a live A/B test reveals how humans actually respond to a new ranking order. This is the bottleneck that the paper addresses.

2. The Offline Simulator and Its Biases

2.1 How the Naive Simulator Works

The simulator constructs a counterfactual estimate by replaying logged user sessions through the modified value model:

def simulate_policy(sessions, new_policy, event_models):
    outcomes = []
    for session in sessions:
        # Re-rank items with new policy weights
        scored_items = [(item, score(item, new_policy)) for item in session.candidates]
        ranked = sorted(scored_items, key=lambda x: -x[1])
        
        # Apply event models to selected items (e.g., top-k shown)
        shown = ranked[:K]
        predicted_events = sum(
            event_models[j].predict(item) for item in shown for j in OUTCOMES
        )
        outcomes.append(aggregate(predicted_events))
    return mean(outcomes)

This requires zero additional infrastructure — it reuses the existing event models and the ranking system itself. But it has two major sources of bias.

Source 1: Off-policy distribution shift. The event models \(\hat{p}_j\) were trained on data generated by the current policy \(\pi_0\). When the simulator evaluates \(\pi \neq \pi_0\), it shows items that the current policy would not have shown, and applies event models that were trained on a different item distribution. Predictions for out-of-distribution items are unreliable.

Source 2: Missing session dynamics. The simulator assumes independence across items and across sessions. In reality: - What a user does on item \(i\) at position \(k\) depends on items \(1, \ldots, k-1\) (cascade / attention effects) - Interactions in one session influence behavior in future sessions (habit formation, satiation)

2.2 Why Simulators Systematically Overestimate

⚠️ The paper makes a subtle empirical claim: simulators don’t just introduce noise — they systematically overestimate the magnitude of policy effects. The authors attribute this to attention and time elasticities (citing Simon 1971, 1996): humans have finite attention capacity, so showing “better” content does not produce proportionally more engagement. If the simulator predicts that switching to policy \(\pi\) increases comments by \(\Delta\), the true online effect is typically \(\alpha \cdot \Delta\) with \(\alpha < 1\).

The ICM model’s response to this: don’t correct the bias explicitly. Instead, model the bias as a function of policy parameters using the flexibility of the GP. The MTGP learns how the offline-to-online discrepancy varies across the design space.

3. Multi-Task Gaussian Process with ICM Kernel

3.1 From Single-Task to Multi-Task GP

Recall the GP posterior (kriging equations). Given observed data \(\{X, y\}\) with noise \(\eta^2\), the predictive distribution at \(x^*\) is

\[f(x^*) \mid x^*, X, y \;\sim\; \mathcal{N}\!\left(\mu(x^*, X, y),\; \sigma^2(x^*, X)\right)\]

\[\mu(x^*, X, y) = m(x^*) + K(x^*, X)\left[K(X,X) + \operatorname{diag}(\eta^2)\right]^{-1}(y - m(X)) \tag{kriging mean}\]

\[\sigma^2(x^*, X) = k(x^*, x^*) - K(x^*, X)\left[K(X,X) + \operatorname{diag}(\eta^2)\right]^{-1}K(X, x^*) \tag{kriging variance}\]

The Multi-Task GP extends this to \(D\) tasks by letting each task \(d\) have its own function \(f_d : \mathcal{X} \to \mathbb{R}\). Points in the combined dataset are indexed by \((d, x)\) pairs. We need a covariance function over the product space \(\{1,\ldots,D\} \times \mathcal{X}\).

3.2 The ICM Kernel

The Intrinsic Coregionalization Model (ICM) makes two assumptions:

Assumption 1 (Shared spatial structure). All tasks share the same base spatial covariance \(\kappa(x, x')\). In this paper, \(\kappa\) is the ARD RBF kernel:

\[\kappa(x, x') = \tau^2 \exp\!\left(-\frac{1}{2}\sum_{j=1}^m \frac{(x_j - x'_j)^2}{\ell_j^2}\right)\]

with per-dimension lengthscales \(\ell_j\) and output scale \(\tau^2\).

Assumption 2 (Separability). The cross-task covariance factorizes as a product of task covariance and spatial covariance:

\[\operatorname{Cov}[f_d(x), f_{d'}(x')] = B_{d,d'} \cdot \kappa(x, x') \tag{ICM}\]

where \(B \in \mathbb{R}^{D \times D}\) is a positive semi-definite task covariance matrix.

The ICM kernel has an equivalent characterization in terms of latent functions. For \(D = 2\) tasks, each task function is a linear combination of \(P\) independent latent GPs \(u_1, \ldots, u_P\) (each with covariance \(\kappa\)):

\[f_1(x) = \sum_{p=1}^P a_{1p}\, u_p(x), \qquad f_2(x) = \sum_{p=1}^P a_{2p}\, u_p(x)\]

The task covariance matrix is then \(B = AA^\top\) where \(A_{dp} = a_{dp}\). The rank of \(B\) equals \(P\), the number of latent functions. Using a low-rank decomposition \(B = \tilde{L}\tilde{L}^\top\) with \(\tilde{L} \in \mathbb{R}^{D \times P}\), \(P < D\), reduces the number of task covariance parameters from \(O(D^2)\) to \(O(DP)\) — important when \(D\) grows as new simulator batches are added.

3.3 Bias as a GP: The Key Structural Insight

💡 From the latent decomposition we get an elegant characterization of the simulator bias. Rearranging for \(P = 2\):

\[f_2(x) = f_1(x) + u'(x) \qquad \text{where } u'(x) = (a_{21} - a_{11})u_1(x) + (a_{22} - a_{12})u_2(x)\]

Since \(u_1, u_2\) are independent GPs with covariance \(\kappa\), their linear combination \(u'\) is also a GP with covariance

\[k_{u'}(x, x') = \left[(a_{21}-a_{11})^2 + (a_{22}-a_{12})^2\right] \kappa(x, x')\]

The ICM model implicitly assumes the offline-to-online bias is a GP with the same smoothness as the response surface itself. This is precisely the right assumption when the bias arises from aggregate behavioral dynamics: smooth, structured, but not zero.

3.4 Hyperparameter Inference

The MTGP has hyperparameters \(\Theta = \{B, \tau^2, \ell_1, \ldots, \ell_m\}\). Inference is by marginal likelihood maximization (or slice sampling for \(B\)):

\[\log p(y \mid X, \Theta) = -\frac{1}{2}y^\top (K_\Theta + \text{diag}(\eta^2))^{-1}y - \frac{1}{2}\log|K_\Theta + \text{diag}(\eta^2)| - \frac{n}{2}\log 2\pi\]

where \(K_\Theta\) is the full \(n \times n\) kernel matrix over all tasks and points. To ensure \(B \succ 0\), it is reparameterized as \(B = LL^\top\) with \(L\) lower-triangular; then \(L\) is optimized.

🔑 The critical advantage of the shared spatial kernel: the lengthscales \(\ell_j\) are inferred from all observations across tasks. Since the simulator provides \(n_S \gg n_T\) observations, the lengthscales are driven almost entirely by offline data. Online experiments can then focus on learning the task relationship \(B\) rather than kernel hyperparameters, which are notoriously hard to infer from \(O(20)\) points in 10–20 dimensions.

4. The Bayesian Optimization Loop

4.1 Noisy EI under Constraints

The optimization problem is: \[\max_{x \in \mathcal{X}} f(x) \quad \text{subject to } c_j(x) \geq 0, \; j = 1, \ldots, J\]

where \(f\) is the objective (e.g., meaningful interactions) and \(c_j\) are constraint metrics (e.g., ecosystem health). All are observed with noise.

The improvement over the current best feasible policy at a new point \(x\) is:

\[I(x, X) = \max\!\left(0,\; f(x) - \max_{x^* \in X: c(x^*) \geq 0} f(x^*)\right) \cdot \mathbf{1}_{c(x) \geq 0}\]

The Noisy Expected Improvement (NEI) acquisition function averages over both observation noise at \(x\) and noise at the existing best feasible:

\[\alpha_\text{NEI}(x \mid X, \mathcal{D}) = \mathbb{E}_{f(X), c(X)}\!\left[\mathbb{E}_{f(x)}\!\left[I(x,X) \mid f(X), c(X)\right] \Bigg\vert \mathcal{D}\right]\]

The inner expectation is analytic (standard EI formula). The outer expectation over the noisy current-best is handled with quasi-Monte Carlo integration.

4.2 Multi-Task Acquisition and the Interleaved Loop

In the multi-task setting, the acquisition function is modified so that only online observations count toward “best feasible”:

\[\alpha_\text{NEI}(x \mid X_T, \mathcal{D}_T \cup \mathcal{D}_S) \tag{MT-NEI}\]

Simulator observations improve the posterior over \(f(X_T)\) and \(c(X_T)\) but are never treated as candidate optimal solutions. This is principled: the goal is to find the best online policy.

The optimization loop interleaves offline and online batches:

# Algorithm 1: Online-Offline Bayesian Optimization
X_T, X_S = quasi_random_batch(n_T=20), quasi_random_batch(n_S=100)
D = run_online(X_T) + run_simulator(X_S)

for k in range(K):
    model = fit_MTGP(D)                         # fit ICM kernel to all data
    X_opt = optimize_acqf(MT_NEI, model, n=30)  # 30 optimized candidates
    
    D_S_k = run_simulator(X_opt)               # test candidates offline
    D = D ∪ D_S_k
    model = fit_MTGP(D)                        # re-fit with offline results
    
    X_T_k = thompson_sample(X_opt, model, n=8) # select best 8–10 for online
    D_T_k = run_online(X_T_k)
    D = D ∪ D_T_k

return best_online_policy(D)

4.3 Simulator Non-Stationarity

The simulator replays live sessions during a short window. Sessions vary by time of day, day of week, and content freshness — so two simulator batches run days apart on the same policy produce different outcomes (Figure 11 in the paper shows this clearly: linear shifts between batches).

The fix: treat each simulator batch as a separate task in the MTGP, using a low-rank ICM kernel to learn the batch-to-batch adjustment. Since the biases are empirically near rank-1 (linear shifts), a rank-1 \(\tilde{L}\) is usually sufficient. The online-offline covariance uses a full-rank (rank-2) \(B\).

5. Reading the Empirical Results

5.1 LOO-CV Scatter Plots (Figs. 2–3)

These are leave-one-out cross-validation scatter plots. Each point corresponds to a single policy that was: 1. Held out from the training set 2. Predicted by the GP/MTGP fit to the remaining observations 3. Actually tested online

The x-axis is the observed online outcome (ground truth); the y-axis is the model’s LOO prediction. A model with perfect predictions concentrates all points on the diagonal \(y = x\). A model that predicts the mean for every point produces a horizontal band around \(y = 0\) (since outcomes are standardized).

5.2 Learning Curves (Fig. 7)

Figure 7 shows mean squared prediction error (MSE) as a function of number of observations, estimated via 500-subsample cross-validation. There are two panels per column: - Top row: MSE vs. number of online observations \(n_T\), for fixed \(n_S \in \{0, 10, 30, 90\}\) - Bottom row: MSE vs. number of offline observations \(n_S\), for fixed \(n_T \in \{4, 10, 18\}\)

The three columns correspond to \(\rho^2 \approx \{0.95, 0.80, 0.09\}\).

📐 Reading the ρ² = 0.95 panel closely: at \(n_T = 2\) online observations, adding 90 offline observations reduces MSE from ~1.4 to ~0.2 — comparable to 18 online observations with no offline. So at high correlation, the offline-to-online information ratio is roughly 90:18 ≈ 5:1 in terms of equivalent samples. This ratio is determined by \(\rho^2\) and the shape of the learning curve, not by the simulator accuracy per se.

5.3 The Kernel Inference Ablation (Fig. 10)

This is the most surprising result in the paper. Figure 10 compares three models:

🔑 The offline kernel model (middle) nearly closes the gap between single-task GP and MTGP. At \(n_T = 18\) online observations, for \(\rho^2 = 0.8\) and \(\rho^2 = 0.09\), the performance difference between “GP + offline kernel” and MTGP is small — while the gap between “GP + online kernel” and MTGP is large.

Interpretation: In 10 dimensions with only 20 online observations, the MLE for ARD lengthscales is highly underdetermined. The optimizer may collapse to near-zero lengthscales (overfitting) or very long lengthscales (underfitting). The offline data, which provides 100 observations, gives a much better-conditioned MLE problem. Sharing the kernel across tasks means this well-estimated kernel gets used for online predictions too — essentially free regularization from the simulator.

6. The ρ² Bound on Multi-Task Generalization

6.1 Chai (2009) Bound

Define the single-task learning curve as the expected predictive variance with \(n\) i.i.d. training points:

\[\kappa(n) = \mathbb{E}_{X, x^*}\!\left[\sigma^2(x^*, X)\right] \tag{single-task}\]

For the MTGP with \(n_T\) primary (online) and \(n_S\) secondary (offline) observations, the multi-task learning curve is:

\[\mathcal{E}_T(\rho, n_T, n_S) = \mathbb{E}_{X_T, X_S, x^*}\!\left[\sigma_T^2(x^*, \rho, X_T, X_S)\right]\]

Under the normalized ICM with \(B = \begin{pmatrix}1 & \rho \\ \rho & 1\end{pmatrix}\), Chai (2009) proved:

\[\sigma_T^2(x^*, \rho, X_T, X_S) \geq \rho^2 \sigma_T^2(x^*, \rho=1, X_T, X_S) + (1-\rho^2)\sigma_T^2(x^*, \rho=0, X_T, X_S)\]

Taking expectations gives Corollary 2:

\[\boxed{\mathcal{E}_T(\rho, n_T, n_S) \geq \rho^2 \kappa(n_T + n_S) + (1-\rho^2)\kappa(n_T)}\]

The two extreme cases recover known single-task results: - \(\rho = 1\) (perfect task correlation): bound is \(\kappa(n_T + n_S)\) — equivalent to a single task with all observations combined - \(\rho = 0\) (independent tasks): bound is \(\kappa(n_T)\) — secondary task observations are ignored

The bound is a lower bound on the multi-task learning curve — i.e., an upper bound on how well the MTGP can do. The key practical property: the single-task curve \(\kappa(n)\) can be estimated from secondary task observations only (Figure 8). Combined with an estimate of \(\rho^2\) from past experiments, this lets you predict MTGP performance before collecting any primary observations.

6.2 Bound vs. Empirical (Fig. 9)

Figure 9 overlays the theoretical Chai bound (computed using offline-estimated \(\kappa\) and \(\rho^2\)) against the empirical MTGP learning curves. Key observations:

1. The bound is consistently optimistic (empirical is below the bound = better than predicted), especially at low \(n_T\). The paper attributes this to improved kernel inference with the MTGP — a mechanism not captured by Chai’s bound, which assumes fixed true kernel parameters.

2. The shape and \(\rho^2\)-ordering are well-predicted. The theory correctly predicts that the \(\rho^2 = 0.95\) curve descends much faster than the \(\rho^2 = 0.09\) curve as \(n_S\) grows.

3. Practical use: Since \(\kappa(n)\) can be estimated from offline data alone, and \(\rho^2\) can be estimated from past experiments, the bound gives a pre-deployment estimate of how many online vs. offline observations you need. This is a real operational tool for experiment sizing.

References