1 Introduction↩︎

In many sequential decision-making problems, we seek to improve performance while ensuring that each action is safe. We may define safety in absolute terms, such as keeping blood pressure below a threshold, but we may instead define it relative to an established baseline. For example, when determining the effectiveness of a new drug for cancer treatment, we may not want patients’ quality of life to be substantially worse than it would have been under the pre-existing standard of care. This type of safety constraint is challenging since the patients’ quality of life when taking the standard drug is not observed and can only be estimated.

We frame this problem using the framework of safe Bayesian optimization (SafeOpt),¹ which optimizes some objective function \(f(\mathbf{x})\), such as information gain, with respect to some constraint: \[\label{eq:constrained-opt} \max_{\mathbf{x}\in\mathcal{X}} f(\mathbf{x}) \text{ s.t. } q(\mathbf{x}) \ge 0.\tag{1}\] SafeOpt optimizes the objective within a set of safe states \(\mathcal{S}\) (where \(\mathcal{X}\subseteq\mathcal{S}\)) and expands \(\mathcal{S}\) as more information is gathered. We focus specifically on SafeOpt with online conformal prediction [2], which provides distribution-free safety guarantees. Conformal SafeOpt ensures that the rate that \(q(\mathbf{x}) \ge 0\) is violated over \(T\) optimization steps is at or below \(\alpha\), where \(\alpha\) is user-specified. See Sec. 3.1.

Potential safe \(\mathbf{x}\) values are typically determined based on the values of \(q\) that have been seen so far. However, we are interested in safety constraints that contain counterfactuals. For example, we may want the value of a variable \(Y\) under a treatment minus the counterfactual \(Y\) under a standard-of-care treatment to be at least some value \(-\omega\). Because \(q\) depends on a counterfactual that is not observed, it is difficult to provide typical SafeOpt-style safety guarantees. Simply making a point estimate of the counterfactual can result in a violate rate above \(\alpha\).

Therefore, we use split conformal prediction to estimate the outcome of the counterfactual \(\mathbf{x}\) under a standard of care [3]. Split conformal prediction creates interval estimates, rather than point estimates. There is a \(100(1-\epsilon)\)% chance that the true value of the counterfactual falls within the interval, where \(\epsilon\) is user-specified. \(\epsilon\) can be leveraged to prove the formal safety guarantee required by SafeOpt methods (proof sketch in Sec. 4.2 and full proof in Appx. 8).

We perform our primary experiments on two settings, both used in [2], [2]. The first uses a chemical reaction simulator [4], which gives us access to the true counterfactual outcomes. The second is the MovieLens dataset [5]. In both cases, we use the method in Sec. 4 to estimate the constraint values for different constraints with counterfactuals. Our baselines include standard (non-safe) Bayesian optimization and an oracle with full knowledge of the counterfactual. We find that our method optimizes the objective while meeting the safety requirement that is provided in Sec. 4.2.

In summary, our contributions include the following. We create a method that ensures that the actions taken by safe Bayesian optimization do not perform \(\omega\) worse than some known, safe policy for \(100(1-\alpha)\)% of samples, where \(\omega\) and \(\alpha\) are user-specified. In Sec. 4 and Appx. 10, we show that the method is robust to different kinds of covariate shifts when re-weighting the calibration set. We prove that the method meets safety guarantees under minimal assumptions, and we bound the conformal coverage gap under assumption violations in Appx. 9. Sec. 6 and Appx. 12 provide experimental evidence that the safety constraint is met for \(100(1-\alpha)\)% of samples and that the proposed method can adapt to different kinds of covariate shifts. We show this for constraints with one and two estimated terms. Finally, we provide a sensitivity analysis to examine how the method performs under different kinds and amounts of mis-specification.

1.1 Motivating example↩︎

Suppose you are running a clinical trial that tests a new drug for reducing systolic blood pressure. You give patients treatments one after the other for \(T\) timesteps. One of the many safety constraints in your study compares the reduction in blood pressure under the new treatment \(Y(\mathbf{x}_t)\) at time \(t\) with the reduction under the standard-of-care policy \(Y(\mathbf{x}_{soc,t})\). You want \(Y(\mathbf{x}_t)\) to be at least as good as \(Y(\mathbf{x}_{soc,t})\) with some tolerance level \(\omega\). Because you cannot ensure that the constraint will never be violated, you specify a small violation rate \(\alpha=0.01\). That is, the safety constraint can be violated \(1\%\) of the time over \(T\) timesteps.

You first gather data from the standard-of-care policy. The data includes your treatment, any covariates, and the blood pressure outcome under the treatment. You split the data into a training set \(\mathcal{D}_{train}\) and a calibration set \(\mathcal{D}_{cal}\). \(\mathcal{D}_{train}\) is used to get a point estimator for \(Y(\mathbf{x}_{soc})\), and \(\mathcal{D}_{cal}\) is used to create quantiles. You then start generating new treatments subject to the safety constraint using SafeOpt. For every new treatment that you generate, you estimate the upper bound of the counterfactual \(Y(\mathbf{x}_{soc})\) and use that to compute your safety constraint at each timestep \(q(\mathbf{x}_t)\).

Sec. 4 lays out the methodology in more detail. We also describe how to adapt the method to different scenarios by re-weighting the calibration data to look like it came from the standard-of-care policy with the test covariates observed during SafeOpt. For example, if the training and calibration data are observational, the calibration data can be re-weighted to look like it came from the standard-of-care policy (Scenario 1). Alternatively, suppose there is covariate shift. The training and calibration data are collected at Hospital A. SafeOpt is run at Hospital B, which serves an older population. By observing the ages of the people in the study, we can re-weight the calibration data (Scenario 2).

2 Related work↩︎

2.1 Safe Bayesian optimization under noisy constraint observations↩︎

Safe Bayesian optimization [1], [2], [6] imposes safety constraints on standard Bayesian optimization, as in Eq. 1 . Our work considers safety constraints with estimated components, implying that the true value of the safety constraint is never actually observed. [7] considers SafeOpt with constraint observations that may be corrupted by a known amount of Gaussian measurement noise and reformulates the safety constraint to be probabilistic. In addition to the deterministic setting laid out in Sec. 5, [2] considers the case when observations of the constraint have been corrupted by a known amount of measurement noise. They reformulate their probabilistic constraint to account for the noise. Other work in constrained Bayesian optimization reformulates hard constraints as probabilistic constraints under noisy constraint observations [8] or deals with the noise through better modeling of the constraint [9], [10]. Our work is distinct since we specifically consider safety constraints under estimation error for a variety of settings.

2.2 Constrained multi-armed bandits↩︎

Safe Bayesian optimization is similar to constrained bandits. For example, [11] and [12] learn a policy while taking actions subject to a safety constraint by maintaining a safe set, which is similar to SafeOpt. Other kinds of constraints include budget constraints [13]–[16] and expectation constraints [17]. In contrast to other work, our setting uses partially observed constraints.

3 Background↩︎

In safe Bayesian optimization, we maximize an objective function \(f(\mathbf{x})\) while obeying a safety constraint, as in Eq. 1 . We assume that an optimal solution to Eq. 1 exists, and we assume that the initial set of safe candidate solutions \(\mathcal{S}_0\) is not empty. All decisions made by the system fall within the safe set, so \(\mathcal{S}_0\) needs to be nonempty: \(\mathcal{S}_0 \subseteq \{ \mathbf{x}\in\mathcal{X}: q(\mathbf{x}) \ge 0 \}\).

Both the safe set and the Bayesian posterior for the objective function are updated at each timestep \(t\). We can do this by creating credible intervals for \(f\) and \(q\), which use Bayesian priors to estimate the ranges of possible values for \(f\) and \(q\) [1], [2]. [2], [2] defines the intervals using Gaussian processes with mean \(\mu\) and standard deviation \(\sigma\) as follows: \[\label{eq:credible-f} \begin{align} \mathcal{I}_f(\mathbf{x}|\mathcal{O}_t) =& [ f_l(\mathbf{x}|\mathcal{O}_t), f_u(\mathbf{x}|\mathcal{O}_t) ] \\ =&[\mu_f(\mathbf{x}|\mathbf{X}_t, \mathbf{f}_t) - \beta_{t+1}\sigma_f(\mathbf{x}|\mathbf{X}_t, \mathbf{f}_t), \mu_f(\mathbf{x}|\mathbf{X}_t, \mathbf{f}_t) + \beta_{t+1}\sigma_f(\mathbf{x}|\mathbf{X}_t, \mathbf{f}_t)] \end{align}\tag{2}\] \[\label{eq:credible-q} \begin{align} \mathcal{I}_q(\mathbf{x}|\mathcal{O}_t) =& [ q_l(\mathbf{x}|\mathcal{O}_t), q_u(\mathbf{x}|\mathcal{O}_t) ] \\ =&[\mu_q(\mathbf{x}|\mathbf{X}_t, \mathbf{q}_t) - \beta_{t+1}\sigma_q(\mathbf{x}|\mathbf{X}_t, \mathbf{q}_t), \mu_q(\mathbf{x}|\mathbf{X}_t, \mathbf{q}_t) + \beta_{t+1}\sigma_q(\mathbf{x}|\mathbf{X}_t, \mathbf{q}_t)] \end{align}\tag{3}\] where \(\mathcal{O}_t=(\mathbf{X}_t, \mathbf{f}_t, \mathbf{q}_t)\) is the history up until time \(t\). \(\mathbf{X}_t\) is the decision history \((\mathbf{x}_1, ..., \mathbf{x}_t)\). \(\mathbf{f}_t\) is the history of \(f(\mathbf{x})\): \((f(\mathbf{x}_1), ..., f(\mathbf{x}_t))\). Similarly, \(\mathbf{q}_t\) is the history of \(q(\mathbf{x})\): \((q(\mathbf{x}_1), ..., q(\mathbf{x}_t))\). \(\beta_{t+1}\) controls the width of the interval and is defined in Sec. 3.1.

We then choose an \(\mathbf{x}_t\) that optimizes the objective function and/or expands the safe set by shrinking the credible interval of either \(f\) or \(q\). We do this by finding the action most decreases our uncertainty about \(f\) or \(q\) (whichever is more uncertain) while maximizing \(f\) or expanding the safe set. \(\mathcal{M}_{t+1}\) is the set of actions that maximizes \(f\), and \(\mathcal{G}_{t+1}\) is the set of actions that could expand the safe set: \[\label{eq:x95t1} \mathbf{x}_{t+1} = \underset{\mathbf{x}\in \mathcal{M}_{t+1}\cup\mathcal{G}_{t+1}}{\arg\max} \max\{ \sigma_f(\mathbf{x}|\mathcal{O}_t), \sigma_q(\mathbf{x}|\mathcal{O}_t) \}.\tag{4}\] The set of optimizers \(\mathcal{M}_{t+1}\) includes all safe solutions that could perform better than the best known pessimistic estimate: \[\label{eq:optimizers} \mathcal{M}_{t+1} = \{ \mathbf{x}\in\mathcal{S}_{t+1} | f_u(\mathbf{x}|\mathcal{O}_t) \ge \max_{\mathbf{x}'\in\mathcal{S}_{t+1}} f_l(\mathbf{x}'|\mathcal{O}_t) \}.\tag{5}\] \(\mathcal{G}_{t+1}\) is the set of solutions that could expand the safe set. That is, we define \(\mathcal{G}_{t+1}\) such that our optimistic prediction of the next safe set \(\tilde{\mathcal{S}}_{t+2}\) is larger than \(\mathcal{S}_{t+1}\). \(\tilde{\mathcal{S}}_{t+2}\) includes both the current history \(\mathcal{O}_t\) and the optimistic observation \((\mathbf{x}, q_u(\mathbf{x}|\mathcal{O}_t))\): \(\tilde{\mathcal{S}}_{t+2}(\mathbf{x}) = \mathcal{S}(\mathcal{O}_t \cup (\mathbf{x}, q_u(\mathbf{x}|\mathcal{O}_t))|\beta_{t+1} )\). \[\label{eq:expanders} \mathcal{G}_{t+1} = \{ \mathbf{x}\in\mathcal{S}_{t+1}: | \tilde{\mathcal{S}}_{t+2}(\mathbf{x}) \backslash \mathcal{S}_{t+1} |>0 \}\tag{6}\]

3.1 Safe optimization with conformal prediction↩︎

[2], [2] shows how to apply online conformal prediction [18], [19] to safe Bayesian optimization. In short, the parameter \(\beta_{t+1}\) in Eq. 2 and 3 is used to define the scale of the safe set so that it becomes more conservative when the safety constraint violation rate is too high. The violation rate is \[\label{eq:violation-rate} \text{violation-rate}(T) := \frac{1}{T}\sum^T_{t=1} \text{err}_t \le \alpha,\tag{7}\] where \(\alpha\) is the acceptable, user-defined violation rate, and the error at each \(t\) is defined as \[\text{err}_t = \indicator(q(\mathbf{x}_t) < 0).\] \(\beta\) is updated as follows: \[\label{eq:beta95t1} \beta_{t+1} = \varphi(\Delta\alpha_{t+1}),\tag{8}\] where \(\Delta\alpha_{t+1}\) is the excess violation rate, which measures the average number of errors above the allowed violation rate \(\alpha_{algo}\). \(\Delta\alpha_{t+1}\) is defined as follows and can be rewritten in terms of the violation rate: \[\begin{align} \begin{aligned} \Delta\alpha_{t+1} &= \Delta\alpha_t + \eta(\text{err}_t - \alpha_{algo}) =\Delta\alpha_1 + \eta\cdot(\sum\limits^t_{t'=1}\text{err}_{t'}-\alpha_{algo}\cdot t)\\ &=\Delta\alpha_1 + \eta\cdot t \cdot(\text{violation-rate}(t) - \alpha_{algo}), \end{aligned} \end{align}\] where \(\eta\) is an update rate hyperparameter; \(\alpha_{algo}\) is the per-step violation rate; and \(\alpha\) is the user-defined violation rate. \(\alpha_{algo}\) is defined as \[\label{eq:alpha-algo} \alpha_{algo} = \frac{1}{T-1}(T\alpha - 1 - \frac{1}{\eta} + \frac{\Delta\alpha_1}{\eta}).\tag{9}\] Finally, \(\varphi(\Delta\alpha_t)\) from Eq. 8 can be defined in a variety of ways, though the following method used in [2], [2] allows for fast adaptation of \(\beta\): \[\varphi(\Delta\alpha_t) = F^{-1}((\text{clip}(\Delta\alpha_t) + 1)/2),\] where \(\text{clip}(\cdot)\) is a function that clips the endpoints of the input between 0 and 1, and \(F^{-1}\) is the inverse cumulative distribution function of a standard Gaussian. \(\beta_t\) increases slowly for low values of \(\Delta\alpha_t\) and quickly as \(\Delta\alpha_t\) approaches 1. See Fig. 3 of [2], [2] for an illustration and [19], [19] for a discussion of adaptation functions.

4 Methodology↩︎

Rather than creating constraints using the observed values of an outcome (or post-treatment covariate), we may wish to create constraints using counterfactuals. Specifically, we may want to ensure that the output of Bayesian optimization is not too much worse than a safe, standard-of-care policy \(\pi_{soc}\) at every timestep \(t\). In our setting, \(\mathbf{x}\) is observed after every timestep. However, the outcome of the standard-of-care policy is not observed since that policy was not applied.

To perform valid conformal prediction, we need to specify a probability of failure \(\epsilon\) in addition to \(\alpha\). For example, we might want to ensure that a constraint of the form \(Y(\mathbf{x}_{t}) - Y(\mathbf{x}_{soc,t})\) is at least \(-\omega\) for 95% of the samples (\(\alpha=0.05\)). However, if \(\epsilon=0.01\), then there is a \(1\%\) chance of an uncaught constraint violation. Therefore, we need to adjust \(\alpha\) to be slightly lower to account for this potential violations. Suppose that \(\mathbf{x}_t\) is the output of Bayesian optimization at time \(t\) and \(\mathbf{x}_{soc,t}\) is chosen by the standard of care policy \(\pi_{soc}\).

4.1 SafeOpt with counterfactual policy constraints (SafeOpt-CPC)↩︎

\(Y(\mathbf{x}_{t}) - Y(\mathbf{x}_{soc,t}) + \omega\) is not measurable due to the counterfactual \(Y(\mathbf{x}_{soc,t})\), so we need to replace \(Y(\mathbf{x}_{soc,t})\) with an estimate when it is not observed. To maintain our safety guarantee (Eq. 7 ), we cannot simply plug in a point prediction for \(Y(\mathbf{x}_{soc,t})\) without risking constraint violations. Rather, we will create a conformal prediction for the counterfactual.

When \(\mathbf{x}_{soc,t}\) is chosen by Bayesian optimization, the safety constraint is met trivially: \[\label{eq:q} q(\mathbf{x}_t) = \begin{cases} Y(\mathbf{x}_t) - \hat{U}^{soc}_{t}(\mathbf{x}_{soc,t}) + \omega & \text{if } \mathbf{x}_t \neq \mathbf{x}_{soc,t} \\ \omega & \text{if } \mathbf{x}_t = \mathbf{x}_{soc,t} \end{cases}\tag{10}\] where \[\label{eq:ideal-conformal-guarantee} \mathbb{P}(Y(\mathbf{x}_{soc,t}) \leq \hat{U}^{soc}_{t}(\mathbf{x}_{soc,t})) \ge 1 - \epsilon.\tag{11}\] We use the upper bound \(\hat{U}^{soc}_{t}(\mathbf{x}_{soc,t})\) in Eq. 10 since in the worst case, the outcome under the standard of care is higher than the actual outcome. Because the outcome under the counterfactual policy is not observed at every step, we cannot apply online conformal prediction based just on the error signal as in [18], [18] and [19], [19]. Rather, we will use offline, split conformal prediction to find \(\hat{U}^{soc}_{t}(\mathbf{x}_{soc,t})\).

In split conformal prediction, there is a training set \(\mathcal{D}_{train}\) that is used to train a point predictor and a calibration set \(\mathcal{D}_{cal}\) of size \(n\) that is used to create quantiles. The point prediction and quantiles can then be used to create an interval estimate for a new data point. In our setting, we will have some pre-collected data that we will split into a training set \(\mathcal{D}_{train}\) and calibration set \(\mathcal{D}_{cal}\), and we will observe our test data sequentially.

To create prediction intervals, we take the quantile of the scores for the calibration set and our new test point. Scores measure how unusual our new test point is compared to our training data. Because we are using an upper one-sided conformal prediction, we define the score \(S\) (Eq. 12 ) as the difference between the true \(Y(\mathbf{x}))\) and the predicted \(\hat{Y}(\mathbf{x})\) without taking the absolute value. \(\hat{Y}(\mathbf{x})\) is a prediction from the point predictor trained on \(\mathcal{D}_{train}\). \[\label{eq:counterfactual-score} S(\mathbf{x}, Y(\mathbf{x})) = Y(\mathbf{x}) - \hat{Y}(\mathbf{x}).\tag{12}\] The scores for the calibration data are expected to be exchangeable with the score of the new test point we are predicting. If they are not exchangeable, we can use weighting to restore exchangeability. In our case the training and calibration data will ideally have been collected by our standard-of-care policy since we are predicting the outcome under the standard of care. If not, we will use weighting.

\(\hat{U}^{soc}_{t}(\mathbf{x}_{soc,t})\) can be defined as \[\label{eq:upper-bound-counterfactual} \begin{align} &\hat{U}^{soc}_{t}(\mathbf{x}_{soc,t}) = \hat{Y}(\mathbf{x}_{soc,t}) + Q^{soc}(1-\epsilon; \sum\limits_{i=1}^n \tilde{w}_i \delta_{S_i^{soc}} + \tilde{w}_{t} \delta_{+\infty}) \end{align}\tag{13}\] where \(\delta_S\) is the point mass at score \(S\); \(\delta_{+\infty}\) is a point mass at a high score corresponding to the new data point at timestep \(t\) of SafeOpt; \(Q^{soc}\) denotes the quantile; and \(\tilde{w}\) are normalized weights [20]: \[\label{eq:norm95w} \tilde{w}_i = w_i/(\sum\limits_{j=1}^n w_j + w_{t}) \text{ and } \tilde{w}_{t} = w_{t}/(\sum\limits_{j=1}^n w_j + w_{t}).\tag{14}\] Note that if the calibration data was collected by the standard-of-care policy, and there is no distribution shift between the calibration set and the new test point, then the weights in Eq. 14 are all set to 1. This is equivalent to unweighted conformal prediction.

If our calibration data did not come from the standard-of-care policy (Scenario 1, Sec. 1.1), we will use inverse propensity scoring. Let \(\pi_{cal}(\mathbf{x}_i)\) be the probability of \(\mathbf{x}_i\) under the policy that collected the calibration data. Then, propensity scoring uses weights of the form: \[\label{eq:w95i} w_i = \pi_{soc}( \mathbf{x}_i)/\pi_{cal}( \mathbf{x}_i)\tag{15}\] Sometimes, there is a covariate shift between the calibration and test data, but no labeled data can be collected to create a new calibration set (Scenario 2, Sec. 1.1). If unlabeled data can be collected, it can be leveraged to learn weights [20]: \[\label{eq:w95i-covariate-shift} w_i =\mathbb{P}_{test}( \mathbf{x}_i) / \mathbb{P}_{cal}( \mathbf{x}_i),\tag{16}\] where \(\mathbb{P}_{test}( \mathbf{x}_i)\) is the probability of \(\mathbf{x}_i\) under the test distribution, and \(\mathbb{P}_{cal}(\mathbf{x}_i)\) is the probability of \(\mathbf{x}_i\) under the calibration distribution.

In Appx. 10 we discuss how non-stationarity affects the weights \(w_i\). We specifically examine two common covariate shift scenarios (changepoints and drift) and offer principled algorithms for both. When a changepoint in the covariates is detected, we can pause SafeOpt and re-calculate the weights before resuming. If drift is detected, we can create time-decaying weights. To avoid the calibration set becoming out-of-date, we take a standard-of-care action every other timestep. See Appx. 10 for details, algorithms, and experiments.

4.1.1 Extensions to different \(q(\mathbf{x})\) constraints↩︎

Eq. 10 assumes that higher values of \(Y\) are better, so we want \(Y(\mathbf{x}_t)\) to be higher than \(Y(\mathbf{x}_{soc,t})\). If we instead want to minimize \(Y\), then we can reformulate the constraint as: \[\label{eq:q-lower} q(\mathbf{x}_t) = \begin{cases} \hat{L}^{soc}_{t}(\mathbf{x}_t) - Y(\mathbf{x}_t) + \omega & \text{if } \mathbf{x}_t \neq \mathbf{x}_{soc,t} \\ \omega & \text{if } \mathbf{x}_t = \mathbf{x}_{soc,t} \end{cases}\tag{17}\] where \(\hat{L}^{soc}_t\) is defined as \[\label{eq:lower-bound-counterfactual} \begin{align} &\hat{L}^{soc}_{t}(\mathbf{x}_t) = \hat{Y}(\mathbf{x}_{soc,t}) - Q^{soc}(1-\epsilon; \sum\limits_{i=1}^n \tilde{w}_i \delta_{S_i^{soc}} + \tilde{w}_{t} \delta_{+\infty}) \end{align}\tag{18}\] \[\label{eq:conformal-guarantee-lower} \mathbb{P}(Y(\mathbf{x}_{soc,t}) \geq \hat{L}^{soc}_{t}(\mathbf{x}_t)) \ge 1 - \epsilon\tag{19}\] The corresponding score function would then be \[\label{eq:score-lower} S(\mathbf{x}, Y(\mathbf{x})) = \hat{Y}(\mathbf{x}) - Y(\mathbf{x}).\tag{20}\] Second, we may wish to define \(q(\mathbf{x})\) with more than one estimated component. Suppose we would like \(Y(\mathbf{x}_t)\) to be higher than the values of \(Y\) under two different standards of care: \(Y(\mathbf{x}_t) - \max(Y(\mathbf{x}_{soc1}), Y(\mathbf{x}_{soc2})) + \omega\). Then, we could define the constraint function as \[\label{eq:q-max} q(\mathbf{x}_t) = \begin{cases} Y(\mathbf{x}_t) - \max(\hat{U}^{soc1}_{t}(\mathbf{x}_t), \hat{U}^{soc2}_{t}(\mathbf{x}_t)) + \omega & \text{if } \mathbf{x}_t \neq \mathbf{x}_{soc1} \neq \mathbf{x}_{soc2} \\ \omega & \text{if } \mathbf{x}_t = \mathbf{x}_{soc1} = \mathbf{x}_{soc2} \end{cases}\tag{21}\] Note that both \(\hat{U}^{soc1}_{t}(\mathbf{x}_t)\) and \(\hat{U}^{soc2}_{t}(\mathbf{x}_t)\) in Eq. 21 have conformal guarantees of the form in Eq. 11 . We can combine them into one conformal guarantee using a union bound: \[\label{eq:conformal-guarantee-multiple} \mathbb{P}(Y(\mathbf{x}_{soc,1}) \leq \hat{U}^{soc,1}_{t}(\mathbf{x}_t) \text{ and } ... \text{ and } Y(\mathbf{x}_{soc,m}) \leq \hat{U}^{soc,m}_{t}(\mathbf{x}_t)) \ge 1 - \epsilon_1 - ... - \epsilon_m.\tag{22}\] In this case \(m=2\). We can then combine \(\epsilon_1 + ... + \epsilon_m = \epsilon'\).

4.1.2 Adapting \(\alpha\)↩︎

Because we are nesting conformal predictions, we will need to adjust the target violation rate \(\alpha\) since the violation rate will be slightly higher than expected due to values of \(Y(\mathbf{x}_{soc,i})\) that fall above \(\hat{U}^{soc,i}\). \(100\epsilon'\)% of the values are expected to cause a potential violation of the safety constraint. (We will deal with the random variation that could cause more violations in Sec. 4.2 and Appx. 8.1.) Therefore, we define a practical maximum violation rate \(\alpha'\): \[\label{eq:alpha-prime} \alpha' + \epsilon(1-\alpha') \le \alpha \equiv \alpha' \le (\alpha - \epsilon)/(1 - \epsilon).\tag{23}\] After computing \(\alpha'\), Eq. 9 can be rewritten as \(\alpha_{algo} = \frac{1}{T-1}(T\alpha' - 1 - \frac{1}{\eta} + \frac{\Delta\alpha_1}{\eta}).\)

As an example, if we want to ensure that the safety constraint is met for 95% of samples and \(\epsilon = 0.01\), then \((0.05 - 0.01)/(1-0.01) \approx 0.04\). In this case, we would target a violation rate of \(\alpha'\approx0.04\) since the stricter target rate would account for the uncaught violations.

Alg. 1 adapts conformal SafeOpt to use counterfactual policy constraints (SafeOpt-CPC).

4.2 Safety proof sketch↩︎

The method in Sec. 4.1 assumes weighted exchangeability ([20], [20], Definition 1), which means that we assume that our weights are correctly specified for the test set. Weighted exchangeability also implies that when non-stationarity is present, it can be immediately detected and the weights can be recalculated. If weighted exchangeability is violated, there is a coverage gap bounded by Eq. 24 or Eq. 25 in Appx. 9.

Thm. 1 states the safety constraint from Eq. 7 is met by the proposed SafeOpt method laid out in Sec. 4. Because we used conformal prediction, we can leverage \(\epsilon\) to guarantee safety under weighted exchangeability. See Appx. 8.1 for the full proof.

Proof sketch. We are performing split conformal prediction to estimate \(Y_{soc}\). If the training and calibration sets come from the standard-of-care policy, then the calibration set is exchangeable with the new test point. If not, we can use weighted split conformal prediction, as described in Sec. 4.1. Therefore, \(\mathbb{P}(Y(\mathbf{x}_{soc,t}) \leq \hat{U}^{soc}_{t}(\mathbf{x}_{soc,t})) \ge 1 - \epsilon\) holds for each conformal estimate. If there is more than one conformal estimate, we can use a union bound to create one combined conformal guarantee, as in Eq. 22 . In particular, we can combine the \(\epsilon\) values to create \(\epsilon_1 + ... + \epsilon_m = \epsilon'\).

SafeOpt with online conformal prediction stays below the violation rate \(\alpha\), as proven by Thm. 2 of [2], [2]. It is expected that about \(100\epsilon'\%\) of the time, Eq. 22 is not met. We define an “excess" as follows: \(excess = \indicator[Y(\mathbf{x}_{soc,1}) > \hat{U}^{soc,1}_{t}(\mathbf{x}_t) \text{ or } ... \text{ or } Y(\mathbf{x}_{soc,m} > \hat{U}^{soc,m}_{t}(\mathbf{x}_t))]\). This means that if there are \(T\)=100 timesteps, \(100\epsilon'\) of the timesteps with no perceived violations are expected to actually contain a violation. We set \(\alpha'\) in Eq. 23 so that the violation rate is by definition less than or equal to \(\alpha\) if the total number of excesses is less than or equal to \(T\epsilon'\). Finally, there is a small probability that the total number of excesses is more than \(T \epsilon'\). We can bound the total number of excesses using the cumulative distribution function of a binomial:
\(\mathbb{P}(\sum^T_t excess_t \le T \epsilon') = \sum^{T\cdot \epsilon'}_{k=0} \binom{T}{k} (\epsilon')^k (1-\epsilon')^{T-k}\). \(\qed\)

The purpose of Algs. 4 and 5 is to calculate better weights under non-stationarity. Because Thm. 1 assumes weighted exchangeability, Algs. 4 and 5 also satisfy the safety constraint. See Corollary 1.

5 Experiments↩︎

We evaluate the proposed method on data generated by a chemical reaction simulator [4] and the MovieLens dataset [5]. The advantage of the chemical reaction simulator is that we have full access to the counterfactuals. For MovieLens, we assume that the counterfactual rating is the same as the user’s rating for the standard-of-care recommendation.

For both the chemical reaction simulator and the MovieLens dataset, we compare SafeOpt-CPC to Bayesian optimization without constraints, as well as to SafeOpt with full knowledge of the counterfactual (an oracle). We consider both the case where the calibration and training sets use the standard of care policy, as well as the case where the calibration and training sets use a different policy. That policy is random in the case of the simulator and is the data collection policy in the case of MovieLens. We examine constraints of the form in Eq. 10 and Eq. 21 for both datasets.

For the chemical reaction experiments, \(f(\mathbf{x})\) is the percent yield, while \(Y\) is the percent selectivity minus the percent yield (see [4]). The standard-of-care policy randomly recommends actions within the range Temperature\(=[100.0,120.0]\) and pH\(=[0.7,0.9]\). The second standard of care policy \(\pi_{soc2}\) (when relevant) recommends Temperature\(=110.0\) and pH=\(0.8\).

For the MovieLens experiment, suppose we recommend movies to specific users chosen from a pool. Each movie has a popularity score of \(Z\), and each user gives each movie a rating \(Y\). We wish to recommend less popular movies while receiving sufficiently high ratings. Therefore, \(f(\mathbf{x})=-Z(\mathbf{x})\), and the safety constraint \(q(\mathbf{x})\) is of the form Eq. 10 or Eq. 21 . The standard-of-care policy randomly recommends one of the top five most popular movies in the training set. The second standard of care policy \(\pi_{soc2}\) (when relevant) recommends the top movie. The pool of test users consists of the five users that (1) have rated all the movies possibly recommended by the standard of care policy and (2) have the most rated movies in common. The pool of movies to recommend consists of the movies these users have in common. This ensures that we have a valid rating for the counterfactual and that the popularity score does not influence which movies are available to recommend at each step. To add realism, the training and calibration sets consist of movies from the pool and users not from the pool, so there is a covariate shift between calibration and testing. We re-weight the calibration set to account for this shift.

Finally, we randomly generate synthetic datasets, as described in Appx. 11.3. The purpose of the synthetic data is to generate a variety of \(f(\mathbf{x})\) and \(q(\mathbf{x})\) functions in case there are any idiosyncrasies in the MovieLens or chemical reaction settings. In short we generate two variables \(Z\) and \(Y\) using different Gaussian processes. \(Y\) is affected by a context variable. \(f(\mathbf{x}) = Z(\mathbf{x})\) and \(q(\mathbf{x})\) is defined in Eq. 10 . We run a sensitivity analysis examining the effects of the quality of \(\hat{Y}(\mathbf{x}_t)\), the amount of noise, and weight misspecification on the results.

See Appx. 11 for further experimental details and Appx. 10 for experiments in nonstationary settings.

6 Results↩︎

As shown in the bar plots in Fig. 2, SafeOpt-CPC identifies potential violations below the rate \(\alpha\). Of these potential violations caught by the algorithm, only a few end up being true violations. The oracle baseline gets the true violation rate closer to the rate \(\alpha\) since it is not as conservative. Standard Bayesian optimization with no safety constraints has a potential and true violation rate far above the allowed violation rate \(\alpha\).

The \(\omega\) tolerance graphs in Fig. 2 show that as the tolerance term increases, SafeOpt-CPC falls back to the standard-of-care less often. \(\omega\) is the tolerance term in Eq. 10 and 21 . Note that for the MovieLens experiments, we can see some sharp drops near the boundary of each star rating. The maximum \(f(\mathbf{x})\) value found increases as the tolerance term increases.

Fig. 2 and Appx. 12.1 show that these trends are the same when the training and calibration set are collected from a different policy. The top row of Fig. 2 shows the results when the training and calibration sets are collected from the standard-of-care policy. The second row shows the results when the sets are from the data collection policy. In both cases, \(q(\mathbf{x})\) is defined as in Eq. 10 .

Fig. 2 and Appx. 12.2 show that the above trends hold for the constraints defined in Eq. 10 and Eq. 21 . The third row of Fig. 2 corresponds to Eq. 10 , and the fourth row corresponds to Eq. 21 .

Overall, these results support the claim that SafeOpt-CPC stays below the allowed violation rate under different data collection policies for training and calibration and under different numbers of estimated components in \(q(\mathbf{x})\).

In Appx. 12.5 we perform a sensitivity analysis and show that the quantiles in Eq. 13 become larger under noise in the data and under high-bias and high-variance estimates of \(Y\). We also show that the coverage gap increases with the weights’ inaccuracy for covariate shift and non-stationarity. Finally, we show that coverage holds under non-stationarity with well-specified weights.

7 Discussion and conclusion↩︎

In our results, we found that SafeOpt-CPC consistently keeps violations below the specified rate, despite the estimated components in the constraint. We also find that the number of potential or perceived violations is generally higher than the number of true violations. As expected, the number of standard of care fallback steps decreases as the tolerance increases, and the optimized value of \(f(\mathbf{x})\) increases as the tolerance increases.

Recall that this work focuses on safety constraints that include estimates of counterfactuals. That said, the split conformal prediction method laid out in Sec. 4 can be trivially extended to any constraint that requires an estimated quantity, rather than an observed one.

A limitation of our work is that if the quantiles are too large and/or the tolerance term is very low, the method may revert back to a given standard of care policy, which can be overly conservative. Further, if there are too many estimated quantities in the constraint, then \(\alpha'\) can become very low, allowing for very few violations. For example, if \(\alpha=0.05\) and there are four estimated quantities each at \(\epsilon=0.01\), then \(\alpha'\) decreases to only \(0.01\).

This work targets applications with strict safety constraints, which is important for high-risk applications of machine learning. We would like to emphasize that controlling the violation rate alone is not sufficient for most high-risk applications, such as medical applications. Though we discuss only the violation rate, this work can be easily extended to other constraint functions. [19], [19] discusses constraints functions other than violation rate, as well as how to control for multiple constraints simultaneously. We leave the extension of these methods to Bayesian optimization to future work.

8 Proofs↩︎

8.1 Proof of Theorem 1↩︎

For this proof, we assume weighted exchangeability, defined in Definition 1 of [20]. Note that if the calibration and test sets are exchangeable, then this assumption is not required.

Second, we will need to prove that Eq. 22 has valid conformal coverage:

\(P(Y(\mathbf{x}_{soc,1}) \le \hat{U}^{soc,1}_{t}(\mathbf{x}_t) \text{ and } ... \text{ and } Y(\mathbf{x}_{soc,m}) \leq \hat{U}^{soc,m}_{t}(\mathbf{x}_t)) \\= 1 - P(Y(\mathbf{x}_{soc,1}) > \hat{U}^{soc,1}_{t}(\mathbf{x}_t) \text{ or } ... \text{ or } Y(\mathbf{x}_{soc,m}) > \hat{U}^{soc,m}_{t}(\mathbf{x}_t))\)

By the union bound,

\(P(Y(\mathbf{x}_{soc,1}) > \hat{U}^{soc,1}_{t}(\mathbf{x}_t) \text{ or } ... \text{ or } Y(\mathbf{x}_{soc,m}) > \hat{U}^{soc,m}_{t}(\mathbf{x}_t)) \\ \le P(Y(\mathbf{x}_{soc,1}) > \hat{U}^{soc,1}_{t}(\mathbf{x}_t)) + ... + P(Y(\mathbf{x}_{soc,m}) > \hat{U}^{soc,m}_{t}(\mathbf{x}_t))\)

By definition,

\(P(Y(\mathbf{x}_{soc,1}) > \hat{U}^{soc,1}_{t}(\mathbf{x}_t)) \leq \epsilon_1\), ..., \(P(Y(\mathbf{x}_{soc,m}) > \hat{U}^{soc,m}_{t}(\mathbf{x}_t)) \leq \epsilon_m\)

Therefore,

\(P(Y(\mathbf{x}_{soc,1}) \le \hat{U}^{soc,1}_{t}(\mathbf{x}_t) \text{ and } ... \text{ and } Y(\mathbf{x}_{soc,m}) \leq \hat{U}^{soc,m}_{t}(\mathbf{x}_t)) \\ \ge 1 - (\epsilon_1 + ... + \epsilon_m) = 1 - \epsilon_1 - ... - \epsilon_m\)

0◻

Theorem 1 Under weighted exchangeability, Alg. 1 satisfies the safety constraint in Eq. 7 under the constraint function \(q(\mathbf{x})\) for user-defined violation rates \(\alpha \in (0,1]\) and \(\epsilon \in (0,1]\) with probability \(\sum^{T\cdot \epsilon'}_{k=0} \binom{T}{k} (\epsilon')^k (1-\epsilon')^{T-k}\).

We are performing split conformal prediction to estimate \(Y_{soc}\). If the training and calibration sets come from the standard-of-care policy, then the calibration set is exchangeable with the new test point.

If not, we can use weighted split conformal prediction, as described in Sec. 4.1. See Sec. 2.2 of [20], [20], which shows that weighted split conformal prediction satisfies Eq. 13 if the weights are well-specified. [21], [21], Thm. 2 proves the coverage gap in Eq. 24 and extends Thm. 2 to data-dependent weights in Sec. 4.5 (Eq. 25 ).

Therefore, \(P(Y(\mathbf{x}_{soc,t}) \leq \hat{U}^{soc}_{t}(\mathbf{x}_{soc,t})) \ge 1 - \epsilon\) holds for each conformal estimate. If there is more than one conformal estimate, we can use a union bound to create one combined conformal guarantee, as proven in Lemma 1. In particular, we can combine the \(\epsilon\) values to create \(\epsilon_1 + ... + \epsilon_m = \epsilon'\).

SafeOpt with online conformal prediction stays below the violation rate \(\alpha\), as proven by Thm. 2 of [2], [2]. It is expected that about \(100\epsilon'\%\) of the time, Eq. 22 is not met. We define an “excess" as follows: \(excess = \indicator[Y(\mathbf{x}_{soc,1}) > \hat{U}^{soc,1}_{t}(\mathbf{x}_t) \text{ or } ... \text{ or } Y(\mathbf{x}_{soc,m} > \hat{U}^{soc,m}_{t}(\mathbf{x}_t))]\). This means that if there are \(T\)=100 timesteps, \(100\epsilon'\) of the timesteps with no perceived violations are expected to actually contain a violation. We set \(\alpha'\) in Eq. 23 so that the violation rate is by definition less than or equal to \(\alpha\) if the total number of excesses is less than or equal to \(T\epsilon'\).

Finally, there is a small probability that the total number of excesses is more than \(T \epsilon'\). We can bound the total number of excesses using the cumulative distribution function of a binomial.

\(\mathbb{P}(\sum^T_t excess_t \le T \epsilon') = \sum^{T\cdot \epsilon'}_{k=0} \binom{T}{k} (\epsilon')^k (1-\epsilon')^{T-k}\)

In this case, we assume that each excess is drawn from a binomial distribution: \(excess\sim Binomial(T, \epsilon')\). We believe that this assumption is reasonable given that the probability of an excess is \(\epsilon'\) at each timestep and given no dependence of the \(\mathbf{x}_{soc}\) values across time. \(\qed\)

The purpose of Algs. 4 and 5 is to calculate better weights under non-stationarity. Because Thm. 1 assumes weighted exchangeability, Algs. 4 and 5 also satisfy the safety constraint.

9 Coverage gaps and analysis↩︎

When the weighted exchangeability assumption is not met, we have a coverage gap: \[P(Y(\mathbf{x}_{soc,t}) \leq \hat{U}^{soc}_{t}(\mathbf{x}_{soc,t})) \ge 1 - \epsilon_t + gap_t\] The gap is measured in terms of total variation distance (TV). If the weights are not data-dependent, then the gap can be bounded by \[\label{eq:gap-not-dependent} gap \le \sum\limits_{i=1}^n \tilde{w_i} \cdot TV(R(\mathcal{D}_{cal+1}), R(\mathcal{D}_{cal+1}^i))\tag{24}\] \(\mathcal{D}_{cal+1}\) is the calibration set with the new test point appended to the end: \((\mathcal{D}_{cal}, \mathbf{x}_t)\). \(\mathcal{D}_{cal+1}^i\) is \(\mathcal{D}_{cal+1}\) with point \(i\) swapped for point \(n+1\). \(R(\mathcal{D}_{cal+1})\) is the vector of residuals for \(\mathcal{D}_{cal+1}\). Note that the residuals and the score are not necessarily the same. If the weights are data-dependent, \[\label{eq:gap-data-dependent} gap \le \mathbb{E}[ \sum\limits_{i=1}^n \tilde{w_i} \cdot TV(R(\mathcal{D}_{cal+1}), R(\mathcal{D}_{cal+1}^i) | w_1,...,w_n) ]\tag{25}\] Data dependent weights are computed from the observed data, rather than being fixed in advance or known from the true distribution. The weight of one point can depend on the weights of the other points, so we condition on the weights and take the expectation [21].² In Eq. 15 , \(P_{cal}\) is usually estimated from the sample, so the weights would be data-dependent. Alternatively, if a practitioner accounts for covariate drift by deciding the decay rate for the weights in advance, then the weights would not be data-dependent.

Note that we cannot calculate the coverage gap in practice because \(Y\) is not observed until after an action is taken. Further, the gap is very expensive to calculate because we need to measure the TV distance for every new test point swapped with every point \(i\) in the calibration set. In theory, you could estimate the gap based on the unlabeled data used to calculate the weights, though you would need to assume no shifts in \(P(Y | \mathbf{x})\) or \(P(\mathbf{x} | Y)\). \[\label{eq:gap-guarantee} \begin{align} &\epsilon_t' = \epsilon + gap_t \\ &P(Y(\mathbf{x}_{soc,t}) \leq \hat{U}^{soc}_{t}(\mathbf{x}_{soc,t})) \ge 1 - \epsilon_t' \end{align}\tag{26}\] The overall chance that \(Y(\mathbf{x}_{soc})\) falls outside of the interval between timesteps 1 and \(T_0\) is therefore \(\epsilon' = (\sum_t^{T_0} \epsilon_t') / T_0\). \(\epsilon'\) could then be substituted for \(\epsilon\) in Eq. 23 .

10 Non-stationarity extension↩︎

SafeOpt with online conformal prediction [2] does not assume stationarity, and we would like to preserve this property as much as possible for safety constraints with estimated components. Not accounting for non-stationarity when it is present will result in a larger coverage gap (see Sec. 9).

Unfortunately, we need to know or be able to detect the kind of non-stationarity in our application. Detection is often possible for non-stationarity in the covariates, so we limit our discussion to covariate shift. We outline two common scenarios, changepoints and drift, below and provide algorithms to account for each scenario. See [21], [21] for a discussion of nonexchangeability, including nonexchangeability from non-stationarity, such as changepoints and drift.

10.1 Motivating example↩︎

Scenario 3: Sudden shift. We begin running SafeOpt at Hospital A, but we move to Hospital B halfway through the study. By observing the ages of the people in the study, we can re-weight the calibration data for the second half of the study.

Scenario 4: Gradual drift. Our training and calibration sets were collected at Hospital A, and the patients recruited and scheduled for the study earlier were from Hospital A. The patients recruited and scheduled later tended to be from Hospital B. We look at when patients are scheduled to receive treatment via SafeOpt and calculate the rate of covariate drift. We then weight each point in the calibration set such that points collected further in the past are less important. To avoid the entire calibration set being far in the past by the end of the study, we alternate between running SafeOpt and running the standard-of-care policy. The standard-of-care step is added to our calibration set.

10.2 Algorithmic description↩︎

If covariate shift is detected at a changepoint (Scenario 3), we can observe unlabeled data and recalculate the weights as in Eq. 16 . After the weights have been calculated, we can create conformal predictions for \(Y_{soc}\) again and SafeOpt continues (Alg. 4).

If there is gradual covariate drift (Scenario 4), we can observe the rate of drift before starting Bayesian optimization so that the weights are set appropriately. Each point in the calibration set is weighted such that points collected further in the past are less important. For example, we could learn a monotone function \(g\) such that \[\label{eq:drift} w_i = \exp(-g(t-t_{i})) \cdot P_{test}( \mathbf{x}_i)/P_{cal}( \mathbf{x}_i)\tag{27}\] where \(t_{i}\) is the time of data collection for point \(i\), and \(t\) is the current timestep. To ensure that the calibration set does not become out of date, we alternate between taking a step of Bayesian optimization and taking an action from \(\pi_{soc}\). This \(\pi_{soc}\) step is added to our calibration set (Alg. 5).

11 Additional experimental details↩︎

Part of our code is based on [6], which uses an MIT License. The MovieLens 100k dataset [5] is available at https://grouplens.org/datasets/movielens/100k/ with a special license listed here: https://files.grouplens.org/datasets/movielens/ml-100k-README.txt. The chemical reaction simulator [4] is available at https://github.com/VlachosGroup/Fructose-HMF-Model. It also has an MIT License.

Each experiment was run on one CPU and used about 300 megabytes. Running 10 seeds took between 4 and 20 minutes, depending on the experiment. For example, experiments with no weighting ran in a few minutes, while experiments with different policies between calibration and test took about 20 minutes.

For the MovieLens 100k dataset, we also wanted to include variation in the data, rather than running the algorithm on the same dataset every time. Therefore, we used bootstrap sampling on the user level to get a slightly different dataset for every random seed.

The results for every dataset or simulator were run on 10 random seeds each. Ninety-five percent confidence intervals were computed by multiplying the standard error of the mean by a critical t-value from scipy.stats.t.ppf().

The hyperparameter \(\eta\) was set to 0.2 for all experiments. In preliminary tests, we found that this was an appropriate value for 100 timesteps. We chose \(\alpha=0.1\) since that value was used in many of the [2] experiments, and it was a high enough value to accommodate a couple of estimated components. \(\epsilon=0.01\) since error rates in conformal prediction are usually set to some small value, and \(\epsilon\) needed to be small enough that \(\alpha'\) was positive.

Because the weights and scores are calculated and stored before the looping through each new test point, not much computational overhead is added compared to standard SafeOpt with online conformal prediction. To calculate the quantiles for each new test point, we need to iterate through at most \(n+1\) data points, where \(n\) is the size of the calibration set.

All code can be found at [Github link redacted].

11.1 MovieLens experiments↩︎

We use matrix factorization on the training set to describe each movie as a set of movie features and each user as a set of user features. User-item pairs from the test set are not included when creating the factorization since they would not have been observed. Each movie is a set of 20 features, and each observed user is a set of 20 features.

We initialize the safe set with one decision \(\mathbf{x}_0\) from the standard-of-care policy and its corresponding observations of \(f(\mathbf{x}_0)\) and \(q(\mathbf{x}_0)\).

In contextual Gaussian processes, the context is often appended to the decision. For the MovieLens experiments, the context is appended to each decision \(\mathbf{x}\). The context includes the standard of care action, followed by the user features. Each decision is a set of item features.

The kernels are factored linear kernels that treat the standard of care action context and the user context differently: \[\begin{align} &k(item, item', item_{soc}, item_{soc}', user, user') = \\& k(item, item') \cdot k(item_{soc}, item_{soc}') + k(item, item') \cdot k(user, user') \end{align}\] where \(k(a, a') = a^\top a'\).

11.2 Chemical reaction experiments↩︎

As in the MovieLens experiments, we initialize the safe set with a decision and observations from the standard-of-care policy. The context only includes the standard-of-care action in this case. The kernel is as follows: \[\begin{align} &k(item, item', item_{soc}, item_{soc}') = k(item, item') \cdot k(item_{soc}, item_{soc}') \end{align}\] where \(k(a, a') = a^\top a'\).

11.3 Synthetic experiments↩︎

To create the synthetic functions \(Z\) and \(Y\), we sample functions from two Gaussian processes. The process for \(Z\) has an RBF kernel with lengthscale 1.0 and variance 1.0. The process for \(Y\) has an RBF kernel variance 1.0. There is a single context variable \(c\) that affects \(Y\). The lengthscale is 0.1 for \(item\), the choice that SafeOpt will make, and 1.0 for \(c\). \[\begin{align} Z& \sim GP(0, RBF(1.0, 1.0)) \\ Y& \sim GP(0, RBF(1.2, (1.0, 0.1))) + 5\cdot c \end{align}\] We seed the Gaussian processes with support points between 0 and 1 for \(item\) and 0 and 5 for \(c\); sample a prior; and then solve for the weights that reconstruct the sampled GP function from the support points. See synthetic_data.py for details.

We initialize the safe set with a decision and observations from the standard-of-care.

The kernel we use for \(q\) is similar to the MovieLens kernel, except \(c\) is our 1D context variable: \[\begin{align} &k(item, item', item_{soc}, item_{soc}', c, c') = \\& k(item, item') \cdot k(item_{soc}, item_{soc}') + k(item, item') \cdot k(c, c') \end{align}\] where \(k(a, a') = a^\top a'\).

12 Additional results↩︎

Note that \(\alpha=0.1\) and \(\epsilon=0.01\) for all plots in this section. We create confidence bounds based on 10 random seeds.

12.1 Additional MovieLens experiments↩︎

Fig. 2 shows the results when the training and calibration data come from the standard-of-care policy or the data collection policy when \(q(\mathbf{x})\) is defined as in Eq. 10 . Fig. 6 shows the full results for both Eq. 10 and Eq. 21 under the standard-of-care policy, and Fig. 7 shows the full results for both Eq. 10 and Eq. 21 under the data collection policy. Note that some results are repeated from Fig. 2. The results do not differ much from the results in the main text.

12.2 Additional chemical reaction experiments↩︎

Fig. 2 shows the chemical reaction results under a random policy for training and calibration when \(q(\mathbf{x})\) is defined as in Eq. 10 and Eq. 21 . Fig. 8 shows the results when the training and calibration data have been collected from the standard-of-care policy. The results do not differ much from the results in the main text.

12.3 Synthetic data results↩︎

We did not include the synthetic data experiments in the main text due to space constraints, but it is useful to include them in order to understand the later experiments in the non-stationarity experiments section and sensitivity analysis section. Because the non-stationarity experiments and sensitivity analysis use \(q(\mathbf{x})\) from Eq. 10 and use training and calibration data from the standard-of-care policy, Fig. 9 also uses these conditions. Note that for each seed, we randomly generate a new \(Z\) and \(Y\).

We observe the same general trends in the synthetic results as we do in the MovieLens and reaction results. The primary difference is that the synthetic results quickly reach near-zero standard-of-care fallback steps and quickly reach the maximum \(f\) value.

12.4 Non-stationarity results↩︎

We add drift to the synthetic data’s covariates so that the covariates transition from (0, 1) to (4, 5) at a rate proportional to \(\gamma^{\Delta t}\), where \(\gamma\) is a parameter plotted in Fig. 10.

The left-hand plot of Fig. 10 shows the violation rate when we correctly re-weight for drift according to \(\gamma^{\Delta t}\). Recent observations are weighted so that they are more relevant, and past observations are less relevant. As described in Sec. 10, we alternate between Bayesian optimization and taking a standard-of-care action for the calibration set. In the right-hand plot, we ignore drift altogether and run Alg. 1 as normal.

In Fig. 11 we examine the effects of changepoints. Before the changepoint, the covariates are between 0 and 1, while after the changepoint, they are between 4 and 5. We plot the violation rate versus the time of the changepoint. When we correctly re-weight for the changepoint, the violation rate stays below \(\alpha\). The violation rate is slightly higher at lower timesteps because we only observe 100 points of new calibration data after the changepoint, so we have better estimates before the changepoint since we have more old calibration data.

In the right-hand plot, we ignore the changepoint. In the worst case, the changepoint occurs early. Because we are not re-weighting, conformal coverage does not hold.

12.5 Sensitivity analysis↩︎

When we have noisy data or we have a high-bias or high-variance estimator, conformal coverage still holds. Therefore, the actual and perceived violations are still below \(\alpha\). However, the number of standard-of-care fallback steps increases as the estimator worsens or as noise increases (Fig. 12). For these plots, \(\omega=0.1\) since the number of fallback steps is low under default conditions, as shown in Fig. 9.

When the weights are mis-specified, conformal coverage does not hold. In Fig. 13 we plot the violation rate versus the percent overlap between the training and calibration and the testing covariates. One hundred percent overlap means that both sets have covariates between 0 and 1. Zero percent overlap means that the training and calibration data have covariates between 0 and 1, while the testing data has covariates between 1 and 2. For this plot, \(\omega=0.001\) since small changes to the weights have a large impact at low \(\omega\) values.

Note that different \(\alpha\) values will change the maximum number of fallback steps and allowed violation rate threshold. High \(\omega\) values result in fewer fallback steps in general and therefore less sensitive to noise or worse \(\hat{Y}\) estimators, while low \(\omega\) values result in more fallback steps.

References↩︎

---

1. SafeOpt sometimes refers to the specific method proposed in [1], [1] and sometimes refers to safe Bayesian optimization in general. In this work, it refers to the latter.↩︎

2. Because we are doing split conformal prediction, we ignore tags.↩︎