Safe Bayesian Optimization with Counterfactual Policies

# Safe Bayesian Optimization with Counterfactual Policies
Source: [https://arxiv.org/html/2607.05620](https://arxiv.org/html/2607.05620)
Katherine Avery College of Computer Science University of Massachusetts Amherst Amherst, MA 01002 kavery@cs\.umass\.edu &Bruno Castro da Silva College of Computer Science University of Massachusetts Amherst Amherst, MA 01002 bsilva@cs\.umass\.edu&David Jensen College of Computer Science University of Massachusetts Amherst Amherst, MA 01002 jensen@cs\.umass\.edu

###### Abstract

In many decision\-making settings, new interventions are acceptable only if they do not reduce outcomes below some established threshold\. For example, in clinical medicine, new treatments are often acceptable only if they do not worsen outcomes relative to an established standard of care\. Safe Bayesian optimization maximizes an objective subject to safety constraints\. In the setting that we consider here, safety is defined relative to a known baseline policy whose outcomes are counterfactual and therefore unobserved\. Thus, the counterfactual outcomes of the baseline policy must be estimated and those \(uncertain\) estimates must be used to safely optimize the objective\. We address this estimation problem by using conformal prediction to construct valid uncertainty intervals for counterfactual baseline outcomes, and we show how these intervals can be integrated into safe Bayesian optimization to ensure that constraint violations occur at or below a user\-specified rate\. We also show how to adapt these conformal estimates to different kinds of covariate shift\. We provide a safety proof, experimental evidence, and a sensitivity analysis\.

## 1Introduction

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\),111SafeOpt sometimes refers to the specific method proposed in\([2015](https://arxiv.org/html/2607.05620#bib.bib7)\)and sometimes refers to safe Bayesian optimization in general\. In this work, it refers to the latter\.which optimizes some objective functionf​\(𝐱\)f\(\\mathbf\{x\}\), such as information gain, with respect to some constraint:

max𝐱∈𝒳⁡f​\(𝐱\)​s\.t\.​q​\(𝐱\)≥0\.\\max\_\{\\mathbf\{x\}\\in\\mathcal\{X\}\}f\(\\mathbf\{x\}\)\\text\{ s\.t\. \}q\(\\mathbf\{x\}\)\\geq 0\.\(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\(Zhanget al\.,[2024](https://arxiv.org/html/2607.05620#bib.bib13)\), which provides distribution\-free safety guarantees\. Conformal SafeOpt ensures that the rate thatq​\(𝐱\)≥0q\(\\mathbf\{x\}\)\\geq 0is violated overTToptimization steps is at or belowα\\alpha, whereα\\alphais user\-specified\. See Sec\.[3\.1](https://arxiv.org/html/2607.05620#S3.SS1)\.

Potential safe𝐱\\mathbf\{x\}values are typically determined based on the values ofqqthat have been seen so far\. However, we are interested in safety constraints that contain counterfactuals\. For example, we may want the value of a variableYYunder a treatment minus the counterfactualYYunder a standard\-of\-care treatment to be at least some value−ω\-\\omega\. Becauseqqdepends 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\(Lei and Candès,[2021](https://arxiv.org/html/2607.05620#bib.bib10)\)\. Split conformal prediction creates interval estimates, rather than point estimates\. There is a100​\(1−ϵ\)100\(1\-\\epsilon\)% chance that the true value of the counterfactual falls within the interval, whereϵ\\epsilonis user\-specified\.ϵ\\epsiloncan be leveraged to prove the formal safety guarantee required by SafeOpt methods \(proof sketch in Sec\.[4\.2](https://arxiv.org/html/2607.05620#S4.SS2)and full proof in Appx\.[A](https://arxiv.org/html/2607.05620#A1)\)\.

We perform our primary experiments on two settings, both used in[Zhanget al\.](https://arxiv.org/html/2607.05620#bib.bib13)\([2024](https://arxiv.org/html/2607.05620#bib.bib13)\)\. The first uses a chemical reaction simulator\(Kanget al\.,[2019](https://arxiv.org/html/2607.05620#bib.bib18)\), which gives us access to the true counterfactual outcomes\. The second is the MovieLens dataset\(Harper and Konstan,[2015](https://arxiv.org/html/2607.05620#bib.bib17)\)\. In both cases, we use the method in Sec\.[4](https://arxiv.org/html/2607.05620#S4)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](https://arxiv.org/html/2607.05620#S4.SS2)\.

In summary, our contributions include the following\. We create a method that ensures that the actions taken by safe Bayesian optimization do not performω\\omegaworse than some known, safe policy for100​\(1−α\)100\(1\-\\alpha\)% of samples, whereω\\omegaandα\\alphaare user\-specified\. In Sec\.[4](https://arxiv.org/html/2607.05620#S4)and Appx\.[C](https://arxiv.org/html/2607.05620#A3), 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\.[B](https://arxiv.org/html/2607.05620#A2)\. Sec\.[6](https://arxiv.org/html/2607.05620#S6)and Appx\.[E](https://arxiv.org/html/2607.05620#A5)provide experimental evidence that the safety constraint is met for100​\(1−α\)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\.1Motivating 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 forTTtimesteps\. One of the many safety constraints in your study compares the reduction in blood pressure under the new treatmentY​\(𝐱t\)Y\(\\mathbf\{x\}\_\{t\}\)at timettwith the reduction under the standard\-of\-care policyY​\(𝐱s​o​c,t\)Y\(\\mathbf\{x\}\_\{soc,t\}\)\. You wantY​\(𝐱t\)Y\(\\mathbf\{x\}\_\{t\}\)to be at least as good asY​\(𝐱s​o​c,t\)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α=0\.01\\alpha=0\.01\. That is, the safety constraint can be violated1%1\\%of the time overTTtimesteps\.

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𝒟t​r​a​i​n\\mathcal\{D\}\_\{train\}and a calibration set𝒟c​a​l\\mathcal\{D\}\_\{cal\}\.𝒟t​r​a​i​n\\mathcal\{D\}\_\{train\}is used to get a point estimator forY​\(𝐱s​o​c\)Y\(\\mathbf\{x\}\_\{soc\}\), and𝒟c​a​l\\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 counterfactualY​\(𝐱s​o​c\)Y\(\\mathbf\{x\}\_\{soc\}\)and use that to compute your safety constraint at each timestepq​\(𝐱t\)q\(\\mathbf\{x\}\_\{t\}\)\.

Sec\.[4](https://arxiv.org/html/2607.05620#S4)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\)\.

## 2Related work

### 2\.1Safe Bayesian optimization under noisy constraint observations

Safe Bayesian optimization\(Suiet al\.,[2015](https://arxiv.org/html/2607.05620#bib.bib7); Berkenkampet al\.,[2023](https://arxiv.org/html/2607.05620#bib.bib3); Zhanget al\.,[2024](https://arxiv.org/html/2607.05620#bib.bib13)\)imposes safety constraints on standard Bayesian optimization, as in Eq\.[1](https://arxiv.org/html/2607.05620#S1.E1)\. Our work considers safety constraints with estimated components, implying that the true value of the safety constraint is never actually observed\.Bergmann and Graichen \([2020](https://arxiv.org/html/2607.05620#bib.bib19)\)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,Zhanget al\.\([2024](https://arxiv.org/html/2607.05620#bib.bib13)\)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\(Gelbartet al\.,[2014](https://arxiv.org/html/2607.05620#bib.bib21)\)or deals with the noise through better modeling of the constraint\(Lethamet al\.,[2019](https://arxiv.org/html/2607.05620#bib.bib20); Wang and Li,[2024](https://arxiv.org/html/2607.05620#bib.bib22)\)\. Our work is distinct since we specifically consider safety constraints under estimation error for a variety of settings\.

### 2\.2Constrained multi\-armed bandits

Safe Bayesian optimization is similar to constrained bandits\. For example,Amaniet al\.\([2019](https://arxiv.org/html/2607.05620#bib.bib23)\)andMoradipariet al\.\([2021](https://arxiv.org/html/2607.05620#bib.bib26)\)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\(Zhou and Tomlin,[2018](https://arxiv.org/html/2607.05620#bib.bib25); Agrawal and Devanur,[2016](https://arxiv.org/html/2607.05620#bib.bib28),[2014](https://arxiv.org/html/2607.05620#bib.bib29); Wuet al\.,[2015](https://arxiv.org/html/2607.05620#bib.bib30)\)and expectation constraints\(Pacchianoet al\.,[2021](https://arxiv.org/html/2607.05620#bib.bib27)\)\. In contrast to other work, our setting uses partially observed constraints\.

## 3Background

In safe Bayesian optimization, we maximize an objective functionf​\(𝐱\)f\(\\mathbf\{x\}\)while obeying a safety constraint, as in Eq\.[1](https://arxiv.org/html/2607.05620#S1.E1)\. We assume that an optimal solution to Eq\.[1](https://arxiv.org/html/2607.05620#S1.E1)exists, and we assume that the initial set of safe candidate solutions𝒮0\\mathcal\{S\}\_\{0\}is not empty\. All decisions made by the system fall within the safe set, so𝒮0\\mathcal\{S\}\_\{0\}needs to be nonempty:𝒮0⊆\{𝐱∈𝒳:q​\(𝐱\)≥0\}\\mathcal\{S\}\_\{0\}\\subseteq\\\{\\mathbf\{x\}\\in\\mathcal\{X\}:q\(\\mathbf\{x\}\)\\geq 0\\\}\.

Both the safe set and the Bayesian posterior for the objective function are updated at each timesteptt\. We can do this by creating credible intervals forffandqq, which use Bayesian priors to estimate the ranges of possible values forffandqq\(Zhanget al\.,[2024](https://arxiv.org/html/2607.05620#bib.bib13); Suiet al\.,[2015](https://arxiv.org/html/2607.05620#bib.bib7)\)\.[Zhanget al\.](https://arxiv.org/html/2607.05620#bib.bib13)\([2024](https://arxiv.org/html/2607.05620#bib.bib13)\)defines the intervals using Gaussian processes with meanμ\\muand standard deviationσ\\sigmaas follows:

ℐf​\(𝐱\|𝒪t\)=\[fl​\(𝐱\|𝒪t\),fu​\(𝐱\|𝒪t\)\]=\[μf​\(𝐱\|𝐗t,𝐟t\)−βt\+1​σf​\(𝐱\|𝐗t,𝐟t\),μf​\(𝐱\|𝐗t,𝐟t\)\+βt\+1​σf​\(𝐱\|𝐗t,𝐟t\)\]\\begin\{split\}\\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\{split\}\(2\)ℐq​\(𝐱\|𝒪t\)=\[ql​\(𝐱\|𝒪t\),qu​\(𝐱\|𝒪t\)\]=\[μq​\(𝐱\|𝐗t,𝐪t\)−βt\+1​σq​\(𝐱\|𝐗t,𝐪t\),μq​\(𝐱\|𝐗t,𝐪t\)\+βt\+1​σq​\(𝐱\|𝐗t,𝐪t\)\]\\begin\{split\}\\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\{split\}\(3\)where𝒪t=\(𝐗t,𝐟t,𝐪t\)\\mathcal\{O\}\_\{t\}=\(\\mathbf\{X\}\_\{t\},\\mathbf\{f\}\_\{t\},\\mathbf\{q\}\_\{t\}\)is the history up until timett\.𝐗t\\mathbf\{X\}\_\{t\}is the decision history\(𝐱1,…,𝐱t\)\(\\mathbf\{x\}\_\{1\},\.\.\.,\\mathbf\{x\}\_\{t\}\)\.𝐟t\\mathbf\{f\}\_\{t\}is the history off​\(𝐱\)f\(\\mathbf\{x\}\):\(f​\(𝐱1\),…,f​\(𝐱t\)\)\(f\(\\mathbf\{x\}\_\{1\}\),\.\.\.,f\(\\mathbf\{x\}\_\{t\}\)\)\. Similarly,𝐪t\\mathbf\{q\}\_\{t\}is the history ofq​\(𝐱\)q\(\\mathbf\{x\}\):\(q​\(𝐱1\),…,q​\(𝐱t\)\)\(q\(\\mathbf\{x\}\_\{1\}\),\.\.\.,q\(\\mathbf\{x\}\_\{t\}\)\)\.βt\+1\\beta\_\{t\+1\}controls the width of the interval and is defined in Sec\.[3\.1](https://arxiv.org/html/2607.05620#S3.SS1)\.

We then choose an𝐱t\\mathbf\{x\}\_\{t\}that optimizes the objective function and/or expands the safe set by shrinking the credible interval of eitherfforqq\. We do this by finding the action most decreases our uncertainty aboutfforqq\(whichever is more uncertain\) while maximizingffor expanding the safe set\.ℳt\+1\\mathcal\{M\}\_\{t\+1\}is the set of actions that maximizesff, and𝒢t\+1\\mathcal\{G\}\_\{t\+1\}is the set of actions that could expand the safe set:

𝐱t\+1=arg⁡max𝐱∈ℳt\+1∪𝒢t\+1​max⁡\{σf​\(𝐱\|𝒪t\),σq​\(𝐱\|𝒪t\)\}\.\\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\}\)\\\}\.\(4\)The set of optimizersℳt\+1\\mathcal\{M\}\_\{t\+1\}includes all safe solutions that could perform better than the best known pessimistic estimate:

ℳt\+1=\{𝐱∈𝒮t\+1\|fu​\(𝐱\|𝒪t\)≥max𝐱′∈𝒮t\+1⁡fl​\(𝐱′\|𝒪t\)\}\.\\mathcal\{M\}\_\{t\+1\}=\\\{\\mathbf\{x\}\\in\\mathcal\{S\}\_\{t\+1\}\|f\_\{u\}\(\\mathbf\{x\}\|\\mathcal\{O\}\_\{t\}\)\\geq\\max\_\{\\mathbf\{x\}^\{\\prime\}\\in\\mathcal\{S\}\_\{t\+1\}\}f\_\{l\}\(\\mathbf\{x\}^\{\\prime\}\|\\mathcal\{O\}\_\{t\}\)\\\}\.\(5\)𝒢t\+1\\mathcal\{G\}\_\{t\+1\}is the set of solutions that could expand the safe set\. That is, we define𝒢t\+1\\mathcal\{G\}\_\{t\+1\}such that our optimistic prediction of the next safe set𝒮~t\+2\\tilde\{\\mathcal\{S\}\}\_\{t\+2\}is larger than𝒮t\+1\\mathcal\{S\}\_\{t\+1\}\.𝒮~t\+2\\tilde\{\\mathcal\{S\}\}\_\{t\+2\}includes both the current history𝒪t\\mathcal\{O\}\_\{t\}and the optimistic observation\(𝐱,qu​\(𝐱\|𝒪t\)\)\(\\mathbf\{x\},q\_\{u\}\(\\mathbf\{x\}\|\\mathcal\{O\}\_\{t\}\)\):𝒮~t\+2​\(𝐱\)=𝒮​\(𝒪t∪\(𝐱,qu​\(𝐱\|𝒪t\)\)\|βt\+1\)\\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\}\)\.

𝒢t\+1=\{𝐱∈𝒮t\+1:\|𝒮~t\+2​\(𝐱\)\\𝒮t\+1\|\>0\}\\mathcal\{G\}\_\{t\+1\}=\\\{\\mathbf\{x\}\\in\\mathcal\{S\}\_\{t\+1\}:\|\\tilde\{\\mathcal\{S\}\}\_\{t\+2\}\(\\mathbf\{x\}\)\\backslash\\mathcal\{S\}\_\{t\+1\}\|\>0\\\}\(6\)
### 3\.1Safe optimization with conformal prediction

[Zhanget al\.](https://arxiv.org/html/2607.05620#bib.bib13)\([2024](https://arxiv.org/html/2607.05620#bib.bib13)\)shows how to apply online conformal prediction\(Gibbs and Candes,[2021](https://arxiv.org/html/2607.05620#bib.bib11); Feldmanet al\.,[2023](https://arxiv.org/html/2607.05620#bib.bib5)\)to safe Bayesian optimization\. In short, the parameterβt\+1\\beta\_\{t\+1\}in Eq\.[2](https://arxiv.org/html/2607.05620#S3.E2)and[3](https://arxiv.org/html/2607.05620#S3.E3)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

violation\-rate​\(T\):=1T​∑t=1Terrt≤α,\\text\{violation\-rate\}\(T\):=\\frac\{1\}\{T\}\\sum^\{T\}\_\{t=1\}\\text\{err\}\_\{t\}\\leq\\alpha,\(7\)whereα\\alphais the acceptable, user\-defined violation rate, and the error at eachttis defined as

errt=1​\(q​\(𝐱t\)<0\)\.\\text\{err\}\_\{t\}=\\text\{1\}\(q\(\\mathbf\{x\}\_\{t\}\)<0\)\.\(8\)β\\betais updated as follows:

βt\+1=φ​\(Δ​αt\+1\),\\beta\_\{t\+1\}=\\varphi\(\\Delta\\alpha\_\{t\+1\}\),\(9\)whereΔ​αt\+1\\Delta\\alpha\_\{t\+1\}is the excess violation rate, which measures the average number of errors above the allowed violation rateαa​l​g​o\\alpha\_\{algo\}\.Δ​αt\+1\\Delta\\alpha\_\{t\+1\}is defined as follows and can be rewritten in terms of the violation rate:

Δ​αt\+1=Δ​αt\+η​\(errt−αa​l​g​o\)=Δ​α1\+η⋅\(∑t′=1terrt′−αa​l​g​o⋅t\)=Δ​α1\+η⋅t⋅\(violation\-rate​\(t\)−αa​l​g​o\),\\displaystyle\\begin\{split\}\\Delta\\alpha\_\{t\+1\}&=\\Delta\\alpha\_\{t\}\+\\eta\(\\text\{err\}\_\{t\}\-\\alpha\_\{algo\}\)=\\Delta\\alpha\_\{1\}\+\\eta\\cdot\(\\sum\\limits^\{t\}\_\{t^\{\\prime\}=1\}\\text\{err\}\_\{t^\{\\prime\}\}\-\\alpha\_\{algo\}\\cdot t\)\\\\ &=\\Delta\\alpha\_\{1\}\+\\eta\\cdot t\\cdot\(\\text\{violation\-rate\}\(t\)\-\\alpha\_\{algo\}\),\\end\{split\}\(10\)whereη\\etais an update rate hyperparameter;αa​l​g​o\\alpha\_\{algo\}is the per\-step violation rate; andα\\alphais the user\-defined violation rate\.αa​l​g​o\\alpha\_\{algo\}is defined as

αa​l​g​o=1T−1​\(T​α−1−1η\+Δ​α1η\)\.\\alpha\_\{algo\}=\\frac\{1\}\{T\-1\}\(T\\alpha\-1\-\\frac\{1\}\{\\eta\}\+\\frac\{\\Delta\\alpha\_\{1\}\}\{\\eta\}\)\.\(11\)Finally,φ​\(Δ​αt\)\\varphi\(\\Delta\\alpha\_\{t\}\)from Eq\.[9](https://arxiv.org/html/2607.05620#S3.E9)can be defined in a variety of ways, though the following method used in[Zhanget al\.](https://arxiv.org/html/2607.05620#bib.bib13)\([2024](https://arxiv.org/html/2607.05620#bib.bib13)\)allows for fast adaptation ofβ\\beta:

φ​\(Δ​αt\)=F−1​\(\(clip​\(Δ​αt\)\+1\)/2\),\\varphi\(\\Delta\\alpha\_\{t\}\)=F^\{\-1\}\(\(\\text\{clip\}\(\\Delta\\alpha\_\{t\}\)\+1\)/2\),\(12\)whereclip​\(⋅\)\\text\{clip\}\(\\cdot\)is a function that clips the endpoints of the input between 0 and 1, andF−1F^\{\-1\}is the inverse cumulative distribution function of a standard Gaussian\.βt\\beta\_\{t\}increases slowly for low values ofΔ​αt\\Delta\\alpha\_\{t\}and quickly asΔ​αt\\Delta\\alpha\_\{t\}approaches 1\. See Fig\. 3 of[Zhanget al\.](https://arxiv.org/html/2607.05620#bib.bib13)\([2024](https://arxiv.org/html/2607.05620#bib.bib13)\)for an illustration and[Feldmanet al\.](https://arxiv.org/html/2607.05620#bib.bib5)\([2023](https://arxiv.org/html/2607.05620#bib.bib5)\)for a discussion of adaptation functions\.

Algorithm 1SafeOpt with Counterfactual Policy Constraints \(SafeOpt\-CPC\)Input:data

𝒟\\mathcal\{D\}, prior models for

ffand

qq, initial safe set

𝒮0\\mathcal\{S\}\_\{0\}, initial decision

𝐱0\\mathbf\{x\}\_\{0\}, total iterations

TT, total observation iterations

T0T\_\{0\}, target rate

α\\alpha, target rate

ϵ\\epsilon, update rate

η\>0\\eta\>0
Split data

𝒟\\mathcal\{D\}into

𝒟t​r​a​i​n\\mathcal\{D\}\_\{train\}and

𝒟c​a​l\\mathcal\{D\}\_\{cal\}
Train estimator

Y^\\hat\{Y\}on

𝒟t​r​a​i​n\\mathcal\{D\}\_\{train\}
ifthere may be a covariate shift between the calibration and test datathen

Observe test data for

T0T\_\{0\}timesteps

endif

for

i=0​to​l​e​n​g​t​h​\(𝒟c​a​l\)i=0\\text\{ to \}length\(\\mathcal\{D\}\_\{cal\}\)do

Compute weights

wiw\_\{i\}according to Eq\.[18](https://arxiv.org/html/2607.05620#S4.E18)or Eq\.[19](https://arxiv.org/html/2607.05620#S4.E19)

endfor

Calculate

∑j=1nwj\\sum\_\{j=1\}^\{n\}w\_\{j\}
Calculate the scores on

𝒟c​a​l\\mathcal\{D\}\_\{cal\}as in Eq\.[15](https://arxiv.org/html/2607.05620#S4.E15)or Eq\.[23](https://arxiv.org/html/2607.05620#S4.E23)

Calculate the violation rate

α′\\alpha^\{\\prime\}according to Eq\.[26](https://arxiv.org/html/2607.05620#S4.E26)

for

t=0t=0to

TTdo

Observe

f​\(𝐱t\)f\(\\mathbf\{x\}\_\{t\}\)and

Y​\(𝐱t\)Y\(\\mathbf\{x\}\_\{t\}\)
if

𝐱t\\mathbf\{x\}\_\{t\}is the same as the standard\(s\) of care or

𝐱t\\mathbf\{x\}\_\{t\}is the decision under some safe fallback policy \(usually

πs​o​c\\pi\_\{soc\}\)then

q​\(𝐱𝐭\)=ωq\(\\mathbf\{x\_\{t\}\}\)=\\omega

else

Compute weight

wtw\_\{t\}for current

𝐱t\\mathbf\{x\}\_\{t\}
Compute the normalized weights

w~i\\tilde\{w\}\_\{i\}and

w~t\\tilde\{w\}\_\{t\}as in Eq\.[17](https://arxiv.org/html/2607.05620#S4.E17)

Form the conformal upper or lower bounds according to Eq\.[16](https://arxiv.org/html/2607.05620#S4.E16)or Eq\.[21](https://arxiv.org/html/2607.05620#S4.E21)

Form

q​\(𝐱𝐭\)q\(\\mathbf\{x\_\{t\}\}\)according to Eq\.[13](https://arxiv.org/html/2607.05620#S4.E13),[20](https://arxiv.org/html/2607.05620#S4.E20), or[24](https://arxiv.org/html/2607.05620#S4.E24)

endif

𝐱t\+1←\\mathbf\{x\}\_\{t\+1\}\\leftarrowSafeOpt\-Helper\(\) \{Continue online conformal SafeOpt \(Alg\.[2](https://arxiv.org/html/2607.05620#alg2)\)\}

endfor

## 4Methodology

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πs​o​c\\pi\_\{soc\}at every timesteptt\. 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ϵ\\epsilonin addition toα\\alpha\. For example, we might want to ensure that a constraint of the formY​\(𝐱t\)−Y​\(𝐱s​o​c,t\)Y\(\\mathbf\{x\}\_\{t\}\)\-Y\(\\mathbf\{x\}\_\{soc,t\}\)is at least−ω\-\\omegafor 95% of the samples \(α=0\.05\\alpha=0\.05\)\. However, ifϵ=0\.01\\epsilon=0\.01, then there is a1%1\\%chance of an uncaught constraint violation\. Therefore, we need to adjustα\\alphato be slightly lower to account for this potential violations\. Suppose that𝐱t\\mathbf\{x\}\_\{t\}is the output of Bayesian optimization at timettand𝐱s​o​c,t\\mathbf\{x\}\_\{soc,t\}is chosen by the standard of care policyπs​o​c\\pi\_\{soc\}\.

### 4\.1SafeOpt with counterfactual policy constraints \(SafeOpt\-CPC\)

Y​\(𝐱t\)−Y​\(𝐱s​o​c,t\)\+ωY\(\\mathbf\{x\}\_\{t\}\)\-Y\(\\mathbf\{x\}\_\{soc,t\}\)\+\\omegais not measurable due to the counterfactualY​\(𝐱s​o​c,t\)Y\(\\mathbf\{x\}\_\{soc,t\}\), so we need to replaceY​\(𝐱s​o​c,t\)Y\(\\mathbf\{x\}\_\{soc,t\}\)with an estimate when it is not observed\. To maintain our safety guarantee \(Eq\.[7](https://arxiv.org/html/2607.05620#S3.E7)\), we cannot simply plug in a point prediction forY​\(𝐱s​o​c,t\)Y\(\\mathbf\{x\}\_\{soc,t\}\)without risking constraint violations\. Rather, we will create a conformal prediction for the counterfactual\.

When𝐱s​o​c,t\\mathbf\{x\}\_\{soc,t\}is chosen by Bayesian optimization, the safety constraint is met trivially:

q​\(𝐱t\)=\{Y​\(𝐱t\)−U^ts​o​c​\(𝐱s​o​c,t\)\+ωif​𝐱t≠𝐱s​o​c,tωif​𝐱t=𝐱s​o​c,tq\(\\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\}\(13\)where

ℙ​\(Y​\(𝐱s​o​c,t\)≤U^ts​o​c​\(𝐱s​o​c,t\)\)≥1−ϵ\.\\mathbb\{P\}\(Y\(\\mathbf\{x\}\_\{soc,t\}\)\\leq\\hat\{U\}^\{soc\}\_\{t\}\(\\mathbf\{x\}\_\{soc,t\}\)\)\\geq 1\-\\epsilon\.\(14\)We use the upper boundU^ts​o​c​\(𝐱s​o​c,t\)\\hat\{U\}^\{soc\}\_\{t\}\(\\mathbf\{x\}\_\{soc,t\}\)in Eq\.[13](https://arxiv.org/html/2607.05620#S4.E13)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[Gibbs and Candes](https://arxiv.org/html/2607.05620#bib.bib11)\([2021](https://arxiv.org/html/2607.05620#bib.bib11)\)and[Feldmanet al\.](https://arxiv.org/html/2607.05620#bib.bib5)\([2023](https://arxiv.org/html/2607.05620#bib.bib5)\)\. Rather, we will use offline, split conformal prediction to findU^ts​o​c​\(𝐱s​o​c,t\)\\hat\{U\}^\{soc\}\_\{t\}\(\\mathbf\{x\}\_\{soc,t\}\)\.

In split conformal prediction, there is a training set𝒟t​r​a​i​n\\mathcal\{D\}\_\{train\}that is used to train a point predictor and a calibration set𝒟c​a​l\\mathcal\{D\}\_\{cal\}of sizennthat 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𝒟t​r​a​i​n\\mathcal\{D\}\_\{train\}and calibration set𝒟c​a​l\\mathcal\{D\}\_\{cal\}, and we will observe our test data sequentially\.

To create prediction intervals, we take the quantile of thescoresfor 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 scoreSS\(Eq\.[15](https://arxiv.org/html/2607.05620#S4.E15)\) as the difference between the trueY\(𝐱\)\)Y\(\\mathbf\{x\}\)\)and the predictedY^​\(𝐱\)\\hat\{Y\}\(\\mathbf\{x\}\)without taking the absolute value\.Y^​\(𝐱\)\\hat\{Y\}\(\\mathbf\{x\}\)is a prediction from the point predictor trained on𝒟t​r​a​i​n\\mathcal\{D\}\_\{train\}\.

S​\(𝐱,Y​\(𝐱\)\)=Y​\(𝐱\)−Y^​\(𝐱\)\.S\(\\mathbf\{x\},Y\(\\mathbf\{x\}\)\)=Y\(\\mathbf\{x\}\)\-\\hat\{Y\}\(\\mathbf\{x\}\)\.\(15\)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\.

U^ts​o​c​\(𝐱s​o​c,t\)\\hat\{U\}^\{soc\}\_\{t\}\(\\mathbf\{x\}\_\{soc,t\}\)can be defined as

U^ts​o​c​\(𝐱s​o​c,t\)=Y^​\(𝐱s​o​c,t\)\+Qs​o​c​\(1−ϵ;∑i=1nw~i​δSis​o​c\+w~t​δ\+∞\)\\begin\{split\}&\\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\{split\}\(16\)whereδS\\delta\_\{S\}is the point mass at scoreSS;δ\+∞\\delta\_\{\+\\infty\}is a point mass at a high score corresponding to the new data point at timestepttof SafeOpt;Qs​o​cQ^\{soc\}denotes the quantile; andw~\\tilde\{w\}are normalized weights\(Tibshiraniet al\.,[2019](https://arxiv.org/html/2607.05620#bib.bib16)\):

w~i=wi/\(∑j=1nwj\+wt\)​and​w~t=wt/\(∑j=1nwj\+wt\)\.\\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\}\)\.\(17\)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\.[17](https://arxiv.org/html/2607.05620#S4.E17)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](https://arxiv.org/html/2607.05620#S1.SS1)\), we will use inverse propensity scoring\. Letπc​a​l​\(𝐱i\)\\pi\_\{cal\}\(\\mathbf\{x\}\_\{i\}\)be the probability of𝐱i\\mathbf\{x\}\_\{i\}under the policy that collected the calibration data\. Then, propensity scoring uses weights of the form:

wi=πs​o​c​\(𝐱i\)/πc​a​l​\(𝐱i\)w\_\{i\}=\\pi\_\{soc\}\(\\mathbf\{x\}\_\{i\}\)/\\pi\_\{cal\}\(\\mathbf\{x\}\_\{i\}\)\(18\)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](https://arxiv.org/html/2607.05620#S1.SS1)\)\. If unlabeled data can be collected, it can be leveraged to learn weights\(Tibshiraniet al\.,[2019](https://arxiv.org/html/2607.05620#bib.bib16)\):

wi=ℙt​e​s​t​\(𝐱i\)/ℙc​a​l​\(𝐱i\),w\_\{i\}=\\mathbb\{P\}\_\{test\}\(\\mathbf\{x\}\_\{i\}\)/\\mathbb\{P\}\_\{cal\}\(\\mathbf\{x\}\_\{i\}\),\(19\)whereℙt​e​s​t​\(𝐱i\)\\mathbb\{P\}\_\{test\}\(\\mathbf\{x\}\_\{i\}\)is the probability of𝐱i\\mathbf\{x\}\_\{i\}under the test distribution, andℙc​a​l​\(𝐱i\)\\mathbb\{P\}\_\{cal\}\(\\mathbf\{x\}\_\{i\}\)is the probability of𝐱i\\mathbf\{x\}\_\{i\}under the calibration distribution\.

In Appx\.[C](https://arxiv.org/html/2607.05620#A3)we discuss how non\-stationarity affects the weightswiw\_\{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\.[C](https://arxiv.org/html/2607.05620#A3)for details, algorithms, and experiments\.

#### 4\.1\.1Extensions to differentq​\(𝐱\)q\(\\mathbf\{x\}\)constraints

Eq\.[13](https://arxiv.org/html/2607.05620#S4.E13)assumes that higher values ofYYare better, so we wantY​\(𝐱t\)Y\(\\mathbf\{x\}\_\{t\}\)to be higher thanY​\(𝐱s​o​c,t\)Y\(\\mathbf\{x\}\_\{soc,t\}\)\. If we instead want to minimizeYY, then we can reformulate the constraint as:

q​\(𝐱t\)=\{L^ts​o​c​\(𝐱t\)−Y​\(𝐱t\)\+ωif​𝐱t≠𝐱s​o​c,tωif​𝐱t=𝐱s​o​c,tq\(\\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\}\(20\)whereL^ts​o​c\\hat\{L\}^\{soc\}\_\{t\}is defined as

L^ts​o​c​\(𝐱t\)=Y^​\(𝐱s​o​c,t\)−Qs​o​c​\(1−ϵ;∑i=1nw~i​δSis​o​c\+w~t​δ\+∞\)\\begin\{split\}&\\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\{split\}\(21\)ℙ​\(Y​\(𝐱s​o​c,t\)≥L^ts​o​c​\(𝐱t\)\)≥1−ϵ\\mathbb\{P\}\(Y\(\\mathbf\{x\}\_\{soc,t\}\)\\geq\\hat\{L\}^\{soc\}\_\{t\}\(\\mathbf\{x\}\_\{t\}\)\)\\geq 1\-\\epsilon\(22\)The corresponding score function would then be

S​\(𝐱,Y​\(𝐱\)\)=Y^​\(𝐱\)−Y​\(𝐱\)\.S\(\\mathbf\{x\},Y\(\\mathbf\{x\}\)\)=\\hat\{Y\}\(\\mathbf\{x\}\)\-Y\(\\mathbf\{x\}\)\.\(23\)Second, we may wish to defineq​\(𝐱\)q\(\\mathbf\{x\}\)with more than one estimated component\. Suppose we would likeY​\(𝐱t\)Y\(\\mathbf\{x\}\_\{t\}\)to be higher than the values ofYYunder two different standards of care:Y​\(𝐱t\)−max⁡\(Y​\(𝐱s​o​c​1\),Y​\(𝐱s​o​c​2\)\)\+ωY\(\\mathbf\{x\}\_\{t\}\)\-\\max\(Y\(\\mathbf\{x\}\_\{soc1\}\),Y\(\\mathbf\{x\}\_\{soc2\}\)\)\+\\omega\. Then, we could define the constraint function as

q​\(𝐱t\)=\{Y​\(𝐱t\)−max⁡\(U^ts​o​c​1​\(𝐱t\),U^ts​o​c​2​\(𝐱t\)\)\+ωif​𝐱t≠𝐱s​o​c​1≠𝐱s​o​c​2ωif​𝐱t=𝐱s​o​c​1=𝐱s​o​c​2q\(\\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\}\(24\)Note that bothU^ts​o​c​1​\(𝐱t\)\\hat\{U\}^\{soc1\}\_\{t\}\(\\mathbf\{x\}\_\{t\}\)andU^ts​o​c​2​\(𝐱t\)\\hat\{U\}^\{soc2\}\_\{t\}\(\\mathbf\{x\}\_\{t\}\)in Eq\.[24](https://arxiv.org/html/2607.05620#S4.E24)have conformal guarantees of the form in Eq\.[14](https://arxiv.org/html/2607.05620#S4.E14)\. We can combine them into one conformal guarantee using a union bound:

ℙ​\(Y​\(𝐱s​o​c,1\)≤U^ts​o​c,1​\(𝐱t\)​and​…​and​Y​\(𝐱s​o​c,m\)≤U^ts​o​c,m​\(𝐱t\)\)≥1−ϵ1−…−ϵm\.\\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\}\)\)\\geq 1\-\\epsilon\_\{1\}\-\.\.\.\-\\epsilon\_\{m\}\.\(25\)In this casem=2m=2\. We can then combineϵ1\+…\+ϵm=ϵ′\\epsilon\_\{1\}\+\.\.\.\+\\epsilon\_\{m\}=\\epsilon^\{\\prime\}\.

#### 4\.1\.2Adaptingα\\alpha

Because we are nesting conformal predictions, we will need to adjust the target violation rateα\\alphasince the violation rate will be slightly higher than expected due to values ofY​\(𝐱s​o​c,i\)Y\(\\mathbf\{x\}\_\{soc,i\}\)that fall aboveU^s​o​c,i\\hat\{U\}^\{soc,i\}\.100​ϵ′100\\epsilon^\{\\prime\}% 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](https://arxiv.org/html/2607.05620#S4.SS2)and Appx\.[A\.1](https://arxiv.org/html/2607.05620#A1.SS1)\.\) Therefore, we define a practical maximum violation rateα′\\alpha^\{\\prime\}:

α′\+ϵ​\(1−α′\)≤α≡α′≤\(α−ϵ\)/\(1−ϵ\)\.\\alpha^\{\\prime\}\+\\epsilon\(1\-\\alpha^\{\\prime\}\)\\leq\\alpha\\equiv\\alpha^\{\\prime\}\\leq\(\\alpha\-\\epsilon\)/\(1\-\\epsilon\)\.\(26\)After computingα′\\alpha^\{\\prime\}, Eq\.[11](https://arxiv.org/html/2607.05620#S3.E11)can be rewritten asαa​l​g​o=1T−1​\(T​α′−1−1η\+Δ​α1η\)\.\\alpha\_\{algo\}=\\frac\{1\}\{T\-1\}\(T\\alpha^\{\\prime\}\-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ϵ=0\.01\\epsilon=0\.01, then\(0\.05−0\.01\)/\(1−0\.01\)≈0\.04\(0\.05\-0\.01\)/\(1\-0\.01\)\\approx 0\.04\. In this case, we would target a violation rate ofα′≈0\.04\\alpha^\{\\prime\}\\approx 0\.04since the stricter target rate would account for the uncaught violations\.

Alg\.[1](https://arxiv.org/html/2607.05620#alg1)adapts conformal SafeOpt to use counterfactual policy constraints \(SafeOpt\-CPC\)\.

### 4\.2Safety proof sketch

The method in Sec\.[4\.1](https://arxiv.org/html/2607.05620#S4.SS1)assumesweighted exchangeability\([Tibshiraniet al\.](https://arxiv.org/html/2607.05620#bib.bib16)\([2019](https://arxiv.org/html/2607.05620#bib.bib16)\), 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\.[28](https://arxiv.org/html/2607.05620#A2.E28)or Eq\.[29](https://arxiv.org/html/2607.05620#A2.E29)in Appx\.[B](https://arxiv.org/html/2607.05620#A2)\.

Thm\.[4\.1](https://arxiv.org/html/2607.05620#S4.Thmtheorem1)states the safety constraint from Eq\.[7](https://arxiv.org/html/2607.05620#S3.E7)is met by the proposed SafeOpt method laid out in Sec\.[4](https://arxiv.org/html/2607.05620#S4)\. Because we used conformal prediction, we can leverageϵ\\epsilonto guarantee safety under weighted exchangeability\. See Appx\.[A\.1](https://arxiv.org/html/2607.05620#A1.SS1)for the full proof\.

###### Theorem 4\.1\.

Under weighted exchangeability, Alg\.[1](https://arxiv.org/html/2607.05620#alg1)satisfies the safety constraint in Eq\.[7](https://arxiv.org/html/2607.05620#S3.E7)under the constraint functionq​\(𝐱\)q\(\\mathbf\{x\}\)for user\-defined violation ratesα∈\(0,1\]\\alpha\\in\(0,1\]andϵ∈\(0,1\]\\epsilon\\in\(0,1\]with probability∑k=0T⋅ϵ′\(Tk\)​\(ϵ′\)k​\(1−ϵ′\)T−k\\sum^\{T\\cdot\\epsilon^\{\\prime\}\}\_\{k=0\}\\binom\{T\}\{k\}\(\\epsilon^\{\\prime\}\)^\{k\}\(1\-\\epsilon^\{\\prime\}\)^\{T\-k\}\.

Proof sketch\.We are performing split conformal prediction to estimateYs​o​cY\_\{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](https://arxiv.org/html/2607.05620#S4.SS1)\. Therefore,ℙ​\(Y​\(𝐱s​o​c,t\)≤U^ts​o​c​\(𝐱s​o​c,t\)\)≥1−ϵ\\mathbb\{P\}\(Y\(\\mathbf\{x\}\_\{soc,t\}\)\\leq\\hat\{U\}^\{soc\}\_\{t\}\(\\mathbf\{x\}\_\{soc,t\}\)\)\\geq 1\-\\epsilonholds 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\.[25](https://arxiv.org/html/2607.05620#S4.E25)\. In particular, we can combine theϵ\\epsilonvalues to createϵ1\+…\+ϵm=ϵ′\\epsilon\_\{1\}\+\.\.\.\+\\epsilon\_\{m\}=\\epsilon^\{\\prime\}\.

SafeOpt with online conformal prediction stays below the violation rateα\\alpha, as proven by Thm\. 2 of[Zhanget al\.](https://arxiv.org/html/2607.05620#bib.bib13)\([2024](https://arxiv.org/html/2607.05620#bib.bib13)\)\. It is expected that about100​ϵ′%100\\epsilon^\{\\prime\}\\%of the time, Eq\.[25](https://arxiv.org/html/2607.05620#S4.E25)is not met\. We define an “excess" as follows:e​x​c​e​s​s=1​\[Y​\(𝐱s​o​c,1\)\>U^ts​o​c,1​\(𝐱t\)​or​…​or​Y​\(𝐱s​o​c,m\>U^ts​o​c,m​\(𝐱t\)\)\]excess=\\text\{1\}\[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 areTT=100 timesteps,100​ϵ′100\\epsilon^\{\\prime\}of the timesteps with no perceived violations are expected to actually contain a violation\. We setα′\\alpha^\{\\prime\}in Eq\.[26](https://arxiv.org/html/2607.05620#S4.E26)so that the violation rate is by definition less than or equal toα\\alphaif the total number of excesses is less than or equal toT​ϵ′T\\epsilon^\{\\prime\}\. Finally, there is a small probability that the total number of excesses is more thanT​ϵ′T\\epsilon^\{\\prime\}\. We can bound the total number of excesses using the cumulative distribution function of a binomial: ℙ​\(∑tTe​x​c​e​s​st≤T​ϵ′\)=∑k=0T⋅ϵ′\(Tk\)​\(ϵ′\)k​\(1−ϵ′\)T−k\\mathbb\{P\}\(\\sum^\{T\}\_\{t\}excess\_\{t\}\\leq T\\epsilon^\{\\prime\}\)=\\sum^\{T\\cdot\\epsilon^\{\\prime\}\}\_\{k=0\}\\binom\{T\}\{k\}\(\\epsilon^\{\\prime\}\)^\{k\}\(1\-\\epsilon^\{\\prime\}\)^\{T\-k\}\.∎

The purpose of Algs\.[3](https://arxiv.org/html/2607.05620#alg3)and[4](https://arxiv.org/html/2607.05620#alg4)is to calculate better weights under non\-stationarity\. Because Thm\.[4\.1](https://arxiv.org/html/2607.05620#S4.Thmtheorem1)assumes weighted exchangeability, Algs\.[3](https://arxiv.org/html/2607.05620#alg3)and[4](https://arxiv.org/html/2607.05620#alg4)also satisfy the safety constraint\. See Corollary[A\.2](https://arxiv.org/html/2607.05620#A1.Thmtheorem2)\.

## 5Experiments

We evaluate the proposed method on data generated by a chemical reaction simulator\(Kanget al\.,[2019](https://arxiv.org/html/2607.05620#bib.bib18)\)and the MovieLens dataset\(Harper and Konstan,[2015](https://arxiv.org/html/2607.05620#bib.bib17)\)\. 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\.[13](https://arxiv.org/html/2607.05620#S4.E13)and Eq\.[24](https://arxiv.org/html/2607.05620#S4.E24)for both datasets\.

For the chemical reaction experiments,f​\(𝐱\)f\(\\mathbf\{x\}\)is the percent yield, whileYYis the percent selectivity minus the percent yield \(seeKanget al\.\([2019](https://arxiv.org/html/2607.05620#bib.bib18)\)\)\. The standard\-of\-care policy randomly recommends actions within the range Temperature=\[100\.0,120\.0\]=\[100\.0,120\.0\]and pH=\[0\.7,0\.9\]=\[0\.7,0\.9\]\. The second standard of care policyπs​o​c​2\\pi\_\{soc2\}\(when relevant\) recommends Temperature=110\.0=110\.0and pH=0\.80\.8\.

For the MovieLens experiment, suppose we recommend movies to specific users chosen from a pool\. Each movie has a popularity score ofZZ, and each user gives each movie a ratingYY\. We wish to recommend less popular movies while receiving sufficiently high ratings\. Therefore,f​\(𝐱\)=−Z​\(𝐱\)f\(\\mathbf\{x\}\)=\-Z\(\\mathbf\{x\}\), and the safety constraintq​\(𝐱\)q\(\\mathbf\{x\}\)is of the form Eq\.[13](https://arxiv.org/html/2607.05620#S4.E13)or Eq\.[24](https://arxiv.org/html/2607.05620#S4.E24)\. 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πs​o​c​2\\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\.[D\.3](https://arxiv.org/html/2607.05620#A4.SS3)\. The purpose of the synthetic data is to generate a variety off​\(𝐱\)f\(\\mathbf\{x\}\)andq​\(𝐱\)q\(\\mathbf\{x\}\)functions in case there are any idiosyncrasies in the MovieLens or chemical reaction settings\. In short we generate two variablesZZandYYusing different Gaussian processes\.YYis affected by a context variable\.f​\(𝐱\)=Z​\(𝐱\)f\(\\mathbf\{x\}\)=Z\(\\mathbf\{x\}\)andq​\(𝐱\)q\(\\mathbf\{x\}\)is defined in Eq\.[13](https://arxiv.org/html/2607.05620#S4.E13)\. We run a sensitivity analysis examining the effects of the quality ofY^​\(𝐱t\)\\hat\{Y\}\(\\mathbf\{x\}\_\{t\}\), the amount of noise, and weight misspecification on the results\.

See Appx\.[D](https://arxiv.org/html/2607.05620#A4)for further experimental details and Appx\.[C](https://arxiv.org/html/2607.05620#A3)for experiments in nonstationary settings\.

## 6Results

As shown in the bar plots in Fig\.[1](https://arxiv.org/html/2607.05620#S6.F1), 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α\\alphasince 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ω\\omegatolerance graphs in Fig\.[1](https://arxiv.org/html/2607.05620#S6.F1)show that as the tolerance term increases, SafeOpt\-CPC falls back to the standard\-of\-care less often\.ω\\omegais the tolerance term in Eq\.[13](https://arxiv.org/html/2607.05620#S4.E13)and[24](https://arxiv.org/html/2607.05620#S4.E24)\. Note that for the MovieLens experiments, we can see some sharp drops near the boundary of each star rating\. The maximumf​\(𝐱\)f\(\\mathbf\{x\}\)value found increases as the tolerance term increases\.

Fig\.[1](https://arxiv.org/html/2607.05620#S6.F1)and Appx\.[E\.1](https://arxiv.org/html/2607.05620#A5.SS1)show that these trends are the same when the training and calibration set are collected from a different policy\. The top row of Fig\.[1](https://arxiv.org/html/2607.05620#S6.F1)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​\(𝐱\)q\(\\mathbf\{x\}\)is defined as in Eq\.[13](https://arxiv.org/html/2607.05620#S4.E13)\.

Fig\.[1](https://arxiv.org/html/2607.05620#S6.F1)and Appx\.[E\.2](https://arxiv.org/html/2607.05620#A5.SS2)show that the above trends hold for the constraints defined in Eq\.[13](https://arxiv.org/html/2607.05620#S4.E13)and Eq\.[24](https://arxiv.org/html/2607.05620#S4.E24)\. The third row of Fig\.[1](https://arxiv.org/html/2607.05620#S6.F1)corresponds to Eq\.[13](https://arxiv.org/html/2607.05620#S4.E13), and the fourth row corresponds to Eq\.[24](https://arxiv.org/html/2607.05620#S4.E24)\.

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 inq​\(𝐱\)q\(\\mathbf\{x\}\)\.

In Appx\.[E\.5](https://arxiv.org/html/2607.05620#A5.SS5)we perform a sensitivity analysis and show that the quantiles in Eq\.[16](https://arxiv.org/html/2607.05620#S4.E16)become larger under noise in the data and under high\-bias and high\-variance estimates ofYY\. 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\.

![Refer to caption](https://arxiv.org/html/2607.05620v1/x1.png)Figure 1:MovieLens and chemical reaction results\.\(Left\) The number of violations and perceived violations for the oracle, standard BO, and SafeOpt\-CPC \(ω=0\.01\\omega=0\.01\)\. \(Middle\) The number of timesteps that SafeOpt\-CPC reverts to the standard\-of\-care policy for different values ofω\\omega\(tolerance levels\)\. \(Right\) The maximum safef​\(𝐱\)f\(\\mathbf\{x\}\)value for different values ofω\\omega\. Note that the ratings for the MovieLens results are normalized between 0 and 1, so the range ofω\\omegais also \[0,1\]\.
## 7Discussion 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 off​\(𝐱\)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](https://arxiv.org/html/2607.05620#S4)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^\{\\prime\}can become very low, allowing for very few violations\. For example, ifα=0\.05\\alpha=0\.05and there are four estimated quantities each atϵ=0\.01\\epsilon=0\.01, thenα′\\alpha^\{\\prime\}decreases to only0\.010\.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\.[Feldmanet al\.](https://arxiv.org/html/2607.05620#bib.bib5)\([2023](https://arxiv.org/html/2607.05620#bib.bib5)\)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\.

## References

- S\. Agrawal and N\. R\. Devanur \(2014\)Bandits with concave rewards and convex knapsacks\.InProceedings of the fifteenth ACM conference on Economics and computation,pp\. 989–1006\.Cited by:[§2\.2](https://arxiv.org/html/2607.05620#S2.SS2.p1.1)\.
- S\. Agrawal and N\. Devanur \(2016\)Linear contextual bandits with knapsacks\.Advances in neural information processing systems29\.Cited by:[§2\.2](https://arxiv.org/html/2607.05620#S2.SS2.p1.1)\.
- S\. Amani, M\. Alizadeh, and C\. Thrampoulidis \(2019\)Linear stochastic bandits under safety constraints\.Advances in Neural Information Processing Systems32\.Cited by:[§2\.2](https://arxiv.org/html/2607.05620#S2.SS2.p1.1)\.
- R\. F\. Barber, E\. J\. Candes, A\. Ramdas, and R\. J\. Tibshirani \(2023\)Conformal prediction beyond exchangeability\.The Annals of Statistics51\(2\),pp\. 816–845\.Cited by:[§A\.1](https://arxiv.org/html/2607.05620#A1.SS1.p13.1),[Appendix B](https://arxiv.org/html/2607.05620#A2.p1.9),[Appendix C](https://arxiv.org/html/2607.05620#A3.p2.1)\.
- D\. Bergmann and K\. Graichen \(2020\)Safe bayesian optimization under unknown constraints\.In2020 59th IEEE Conference on Decision and Control \(CDC\),pp\. 3592–3597\.Cited by:[§2\.1](https://arxiv.org/html/2607.05620#S2.SS1.p1.1)\.
- F\. Berkenkamp, A\. Krause, and A\. P\. Schoellig \(2023\)Bayesian optimization with safety constraints: safe and automatic parameter tuning in robotics\.Machine Learning112\(10\),pp\. 3713–3747\.Cited by:[Appendix D](https://arxiv.org/html/2607.05620#A4.p1.1),[§2\.1](https://arxiv.org/html/2607.05620#S2.SS1.p1.1)\.
- S\. Feldman, L\. Ringel, S\. Bates, and Y\. Romano \(2023\)Achieving risk control in online learning settings\.Transactions on Machine Learning Research2023\.Cited by:[§3\.1](https://arxiv.org/html/2607.05620#S3.SS1.p1.1),[§3\.1](https://arxiv.org/html/2607.05620#S3.SS1.p1.18),[§4\.1](https://arxiv.org/html/2607.05620#S4.SS1.p2.3),[§7](https://arxiv.org/html/2607.05620#S7.p4.1)\.
- M\. A\. Gelbart, J\. Snoek, and R\. P\. Adams \(2014\)Bayesian optimization with unknown constraints\.arXiv preprint arXiv:1403\.5607\.Cited by:[§2\.1](https://arxiv.org/html/2607.05620#S2.SS1.p1.1)\.
- I\. Gibbs and E\. Candes \(2021\)Adaptive conformal inference under distribution shift\.Advances in Neural Information Processing Systems34,pp\. 1660–1672\.Cited by:[§3\.1](https://arxiv.org/html/2607.05620#S3.SS1.p1.1),[§4\.1](https://arxiv.org/html/2607.05620#S4.SS1.p2.3)\.
- F\. M\. Harper and J\. A\. Konstan \(2015\)The movielens datasets: history and context\.ACM Trans\. Interact\. Intell\. Syst\.5\(4\)\.External Links:ISSN 2160\-6455Cited by:[Appendix D](https://arxiv.org/html/2607.05620#A4.p1.1),[§1](https://arxiv.org/html/2607.05620#S1.p5.1),[§5](https://arxiv.org/html/2607.05620#S5.p1.1)\.
- P\. Kang, C\. Shang, and Z\. Liu \(2019\)Glucose to 5\-hydroxymethylfurfural: origin of site\-selectivity resolved by machine learning based reaction sampling\.Journal of the American Chemical Society141\(51\),pp\. 20525–20536\.Cited by:[Appendix D](https://arxiv.org/html/2607.05620#A4.p1.1),[§1](https://arxiv.org/html/2607.05620#S1.p5.1),[§5](https://arxiv.org/html/2607.05620#S5.p1.1),[§5](https://arxiv.org/html/2607.05620#S5.p3.7)\.
- L\. Lei and E\. J\. Candès \(2021\)Conformal inference of counterfactuals and individual treatment effects\.Journal of the Royal Statistical Society Series B: Statistical Methodology83\(5\),pp\. 911–938\.Cited by:[§1](https://arxiv.org/html/2607.05620#S1.p4.4)\.
- B\. Letham, B\. Karrer, G\. Ottoni, and E\. Bakshy \(2019\)Constrained bayesian optimization with noisy experiments\.Cited by:[§2\.1](https://arxiv.org/html/2607.05620#S2.SS1.p1.1)\.
- A\. Moradipari, S\. Amani, M\. Alizadeh, and C\. Thrampoulidis \(2021\)Safe linear thompson sampling with side information\.IEEE Transactions on Signal Processing69,pp\. 3755–3767\.Cited by:[§2\.2](https://arxiv.org/html/2607.05620#S2.SS2.p1.1)\.
- A\. Pacchiano, M\. Ghavamzadeh, P\. Bartlett, and H\. Jiang \(2021\)Stochastic bandits with linear constraints\.InProceedings of The 24th International Conference on Artificial Intelligence and Statistics,A\. Banerjee and K\. Fukumizu \(Eds\.\),Proceedings of Machine Learning Research, Vol\.130,pp\. 2827–2835\.Cited by:[§2\.2](https://arxiv.org/html/2607.05620#S2.SS2.p1.1)\.
- Y\. Sui, A\. Gotovos, J\. Burdick, and A\. Krause \(2015\)Safe exploration for optimization with gaussian processes\.InInternational Conference on Machine Learning,pp\. 997–1005\.Cited by:[§2\.1](https://arxiv.org/html/2607.05620#S2.SS1.p1.1),[§3](https://arxiv.org/html/2607.05620#S3.p2.7),[footnote 1](https://arxiv.org/html/2607.05620#footnote1)\.
- R\. J\. Tibshirani, R\. Foygel Barber, E\. Candes, and A\. Ramdas \(2019\)Conformal prediction under covariate shift\.Advances in neural information processing systems32\.Cited by:[§A\.1](https://arxiv.org/html/2607.05620#A1.SS1.p1.1),[§A\.1](https://arxiv.org/html/2607.05620#A1.SS1.p13.1),[§4\.1](https://arxiv.org/html/2607.05620#S4.SS1.p5.7),[§4\.1](https://arxiv.org/html/2607.05620#S4.SS1.p6.7),[§4\.2](https://arxiv.org/html/2607.05620#S4.SS2.p1.1)\.
- S\. Wang and K\. Li \(2024\)Constrained bayesian optimization under partial observations: balanced improvements and provable convergence\.InProceedings of the AAAI Conference on Artificial Intelligence,Vol\.38,pp\. 15607–15615\.Cited by:[§2\.1](https://arxiv.org/html/2607.05620#S2.SS1.p1.1)\.
- H\. Wu, R\. Srikant, X\. Liu, and C\. Jiang \(2015\)Algorithms with logarithmic or sublinear regret for constrained contextual bandits\.Advances in Neural Information Processing Systems28\.Cited by:[§2\.2](https://arxiv.org/html/2607.05620#S2.SS2.p1.1)\.
- Y\. Zhang, S\. Park, and O\. Simeone \(2024\)Bayesian optimization with formal safety guarantees via online conformal prediction\.IEEE Journal of Selected Topics in Signal Processing19\(1\),pp\. 45–59\.Cited by:[§A\.1](https://arxiv.org/html/2607.05620#A1.SS1.p15.8),[Appendix C](https://arxiv.org/html/2607.05620#A3.p1.1),[Appendix D](https://arxiv.org/html/2607.05620#A4.p5.5),[§1](https://arxiv.org/html/2607.05620#S1.p2.8),[§1](https://arxiv.org/html/2607.05620#S1.p5.1),[§2\.1](https://arxiv.org/html/2607.05620#S2.SS1.p1.1),[§3\.1](https://arxiv.org/html/2607.05620#S3.SS1.p1.1),[§3\.1](https://arxiv.org/html/2607.05620#S3.SS1.p1.13),[§3\.1](https://arxiv.org/html/2607.05620#S3.SS1.p1.18),[§3](https://arxiv.org/html/2607.05620#S3.p2.7),[§4\.2](https://arxiv.org/html/2607.05620#S4.SS2.p4.11)\.
- D\. Zhou and C\. Tomlin \(2018\)Budget\-constrained multi\-armed bandits with multiple plays\.InProceedings of the AAAI Conference on Artificial Intelligence,Vol\.32\.Cited by:[§2\.2](https://arxiv.org/html/2607.05620#S2.SS2.p1.1)\.

Algorithm 2SafeOpt\-HelperInput:

𝐱t\\mathbf\{x\}\_\{t\},

q​\(𝐱t\)q\(\\mathbf\{x\}\_\{t\}\),

f​\(𝐱t\)f\(\\mathbf\{x\}\_\{t\}\),

TT,

tt,

α\\alpha,

ϵ\\epsilon,

η\\eta,

𝒪t−1\\mathcal\{O\}\_\{t\-1\},

Δ​αt\\Delta\\alpha\_\{t\}, current models for

ffand

qq, current safe set

𝒮\\mathcal\{S\}
Calculate

errt=1​\(q​\(𝐱𝐭\)<0\)\\text\{err\}\_\{t\}=\\text\{1\}\(q\(\\mathbf\{x\_\{t\}\}\)<0\)
if

t==0t==0then

𝒪t=\{𝐱t,f​\(𝐱t\),q​\(𝐱t\)\}\\mathcal\{O\}\_\{t\}=\\\{\\mathbf\{x\}\_\{t\},f\(\\mathbf\{x\}\_\{t\}\),q\(\\mathbf\{x\}\_\{t\}\)\\\}

else

Update

𝒪t=𝒪t−1∪\{𝐱t,f​\(𝐱t\),q​\(𝐱t\)\}\\mathcal\{O\}\_\{t\}=\\mathcal\{O\}\_\{t\-1\}\\cup\\\{\\mathbf\{x\}\_\{t\},f\(\\mathbf\{x\}\_\{t\}\),q\(\\mathbf\{x\}\_\{t\}\)\\\}
endif

Update the models \(e\.g\., Gaussian processes\) for

ffand

qqwith

𝒪t\\mathcal\{O\}\_\{t\}
Update

Δ​αt\+1=Δ​αt\+η​\(errt−αa​l​g​o\)\\Delta\\alpha\_\{t\+1\}=\\Delta\\alpha\_\{t\}\+\\eta\(\\text\{err\}\_\{t\}\-\\alpha\_\{algo\}\)
Update

βt\+1=φ​\(Δ​αt\+1\)\\beta\_\{t\+1\}=\\varphi\(\\Delta\\alpha\_\{t\+1\}\)
Create credible intervals

ℐf​\(𝐱\|𝒪t\)\\mathcal\{I\}\_\{f\}\(\\mathbf\{x\}\|\\mathcal\{O\}\_\{t\}\)and

ℐq​\(𝐱\|𝒪t\)\\mathcal\{I\}\_\{q\}\(\\mathbf\{x\}\|\\mathcal\{O\}\_\{t\}\)
Update safe set

𝒮t\+1\\mathcal\{S\}\_\{t\+1\}
Create set of optimizers

ℳt\+1\\mathcal\{M\}\_\{t\+1\}and expanders

𝒢t\+1\\mathcal\{G\}\_\{t\+1\}
if

t==Tt==Tthen

𝐱t\+1=arg⁡max𝐱∈St\+1⁡fl​\(𝐱\|𝒪T\)\\mathbf\{x\}\_\{t\+1\}=\\arg\\max\_\{\\mathbf\{x\}\\in S\_\{t\+1\}\}f\_\{l\}\(\\mathbf\{x\}\|\\mathcal\{O\}\_\{T\}\)

else

Form decision

𝐱t\+1\\mathbf\{x\}\_\{t\+1\}as in Eq\.[4](https://arxiv.org/html/2607.05620#S3.E4)

endif

## Appendix AProofs

### A\.1Proof of Theorem[4\.1](https://arxiv.org/html/2607.05620#S4.Thmtheorem1)

For this proof, we assumeweighted exchangeability, defined in Definition 1 ofTibshiraniet al\.\[[2019](https://arxiv.org/html/2607.05620#bib.bib16)\]\. Note that if the calibration and test sets are exchangeable, then this assumption is not required\.

Second, we will need to prove that Eq\.[25](https://arxiv.org/html/2607.05620#S4.E25)has valid conformal coverage:

###### Lemma A\.1\.

Given that we havemmdifferent conformal guarantees of the form P​\(Y​\(𝐱s​o​c,i\)≤U^ts​o​c,i​\(𝐱t\)\)≥1−ϵiP\(Y\(\\mathbf\{x\}\_\{soc,i\}\)\\leq\\hat\{U\}^\{soc,i\}\_\{t\}\(\\mathbf\{x\}\_\{t\}\)\)\\geq 1\-\\epsilon\_\{i\}, then P​\(Y​\(𝐱s​o​c,1\)≤U^ts​o​c,1​\(𝐱t\)​and​…​and​Y​\(𝐱s​o​c,m\)≤U^ts​o​c,m​\(𝐱t\)\)≥1−ϵ1−…−ϵmP\(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\}\)\)\\geq 1\-\\epsilon\_\{1\}\-\.\.\.\-\\epsilon\_\{m\}\(Eq\.[25](https://arxiv.org/html/2607.05620#S4.E25)\)\.

P​\(Y​\(𝐱s​o​c,1\)≤U^ts​o​c,1​\(𝐱t\)​and​…​and​Y​\(𝐱s​o​c,m\)≤U^ts​o​c,m​\(𝐱t\)\)=1−P​\(Y​\(𝐱s​o​c,1\)\>U^ts​o​c,1​\(𝐱t\)​or​…​or​Y​\(𝐱s​o​c,m\)\>U^ts​o​c,m​\(𝐱t\)\)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\}\)\)\\\\ =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​\(𝐱s​o​c,1\)\>U^ts​o​c,1​\(𝐱t\)​or​…​or​Y​\(𝐱s​o​c,m\)\>U^ts​o​c,m​\(𝐱t\)\)≤P​\(Y​\(𝐱s​o​c,1\)\>U^ts​o​c,1​\(𝐱t\)\)\+…\+P​\(Y​\(𝐱s​o​c,m\)\>U^ts​o​c,m​\(𝐱t\)\)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\}\)\)\\\\ \\leq 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​\(𝐱s​o​c,1\)\>U^ts​o​c,1​\(𝐱t\)\)≤ϵ1P\(Y\(\\mathbf\{x\}\_\{soc,1\}\)\>\\hat\{U\}^\{soc,1\}\_\{t\}\(\\mathbf\{x\}\_\{t\}\)\)\\leq\\epsilon\_\{1\}, …,P​\(Y​\(𝐱s​o​c,m\)\>U^ts​o​c,m​\(𝐱t\)\)≤ϵmP\(Y\(\\mathbf\{x\}\_\{soc,m\}\)\>\\hat\{U\}^\{soc,m\}\_\{t\}\(\\mathbf\{x\}\_\{t\}\)\)\\leq\\epsilon\_\{m\}

Therefore,

P​\(Y​\(𝐱s​o​c,1\)≤U^ts​o​c,1​\(𝐱t\)​and​…​and​Y​\(𝐱s​o​c,m\)≤U^ts​o​c,m​\(𝐱t\)\)≥1−\(ϵ1\+…\+ϵm\)=1−ϵ1−…−ϵmP\(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\}\)\)\\\\ \\geq 1\-\(\\epsilon\_\{1\}\+\.\.\.\+\\epsilon\_\{m\}\)=1\-\\epsilon\_\{1\}\-\.\.\.\-\\epsilon\_\{m\}

∎

Theorem[4\.1](https://arxiv.org/html/2607.05620#S4.Thmtheorem1)Under weighted exchangeability, Alg\.[1](https://arxiv.org/html/2607.05620#alg1)satisfies the safety constraint in Eq\.[7](https://arxiv.org/html/2607.05620#S3.E7)under the constraint functionq​\(𝐱\)q\(\\mathbf\{x\}\)for user\-defined violation ratesα∈\(0,1\]\\alpha\\in\(0,1\]andϵ∈\(0,1\]\\epsilon\\in\(0,1\]with probability∑k=0T⋅ϵ′\(Tk\)​\(ϵ′\)k​\(1−ϵ′\)T−k\\sum^\{T\\cdot\\epsilon^\{\\prime\}\}\_\{k=0\}\\binom\{T\}\{k\}\(\\epsilon^\{\\prime\}\)^\{k\}\(1\-\\epsilon^\{\\prime\}\)^\{T\-k\}\.

We are performing split conformal prediction to estimateYs​o​cY\_\{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](https://arxiv.org/html/2607.05620#S4.SS1)\. See Sec\. 2\.2 of[Tibshiraniet al\.](https://arxiv.org/html/2607.05620#bib.bib16)\[[2019](https://arxiv.org/html/2607.05620#bib.bib16)\], which shows that weighted split conformal prediction satisfies Eq\.[16](https://arxiv.org/html/2607.05620#S4.E16)if the weights are well\-specified\.[Barberet al\.](https://arxiv.org/html/2607.05620#bib.bib15)\[[2023](https://arxiv.org/html/2607.05620#bib.bib15)\], Thm\. 2 proves the coverage gap in Eq\.[28](https://arxiv.org/html/2607.05620#A2.E28)and extends Thm\. 2 to data\-dependent weights in Sec\. 4\.5 \(Eq\.[29](https://arxiv.org/html/2607.05620#A2.E29)\)\.

Therefore,P​\(Y​\(𝐱s​o​c,t\)≤U^ts​o​c​\(𝐱s​o​c,t\)\)≥1−ϵP\(Y\(\\mathbf\{x\}\_\{soc,t\}\)\\leq\\hat\{U\}^\{soc\}\_\{t\}\(\\mathbf\{x\}\_\{soc,t\}\)\)\\geq 1\-\\epsilonholds 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[A\.1](https://arxiv.org/html/2607.05620#A1.Thmtheorem1)\. In particular, we can combine theϵ\\epsilonvalues to createϵ1\+…\+ϵm=ϵ′\\epsilon\_\{1\}\+\.\.\.\+\\epsilon\_\{m\}=\\epsilon^\{\\prime\}\.

SafeOpt with online conformal prediction stays below the violation rateα\\alpha, as proven by Thm\. 2 of[Zhanget al\.](https://arxiv.org/html/2607.05620#bib.bib13)\[[2024](https://arxiv.org/html/2607.05620#bib.bib13)\]\. It is expected that about100​ϵ′%100\\epsilon^\{\\prime\}\\%of the time, Eq\.[25](https://arxiv.org/html/2607.05620#S4.E25)is not met\. We define an “excess" as follows:e​x​c​e​s​s=1​\[Y​\(𝐱s​o​c,1\)\>U^ts​o​c,1​\(𝐱t\)​or​…​or​Y​\(𝐱s​o​c,m\>U^ts​o​c,m​\(𝐱t\)\)\]excess=\\text\{1\}\[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 areTT=100 timesteps,100​ϵ′100\\epsilon^\{\\prime\}of the timesteps with no perceived violations are expected to actually contain a violation\. We setα′\\alpha^\{\\prime\}in Eq\.[26](https://arxiv.org/html/2607.05620#S4.E26)so that the violation rate is by definition less than or equal toα\\alphaif the total number of excesses is less than or equal toT​ϵ′T\\epsilon^\{\\prime\}\.

Finally, there is a small probability that the total number of excesses is more thanT​ϵ′T\\epsilon^\{\\prime\}\. We can bound the total number of excesses using the cumulative distribution function of a binomial\.

ℙ​\(∑tTe​x​c​e​s​st≤T​ϵ′\)=∑k=0T⋅ϵ′\(Tk\)​\(ϵ′\)k​\(1−ϵ′\)T−k\\mathbb\{P\}\(\\sum^\{T\}\_\{t\}excess\_\{t\}\\leq T\\epsilon^\{\\prime\}\)=\\sum^\{T\\cdot\\epsilon^\{\\prime\}\}\_\{k=0\}\\binom\{T\}\{k\}\(\\epsilon^\{\\prime\}\)^\{k\}\(1\-\\epsilon^\{\\prime\}\)^\{T\-k\}

In this case, we assume that each excess is drawn from a binomial distribution:e​x​c​e​s​s∼B​i​n​o​m​i​a​l​\(T,ϵ′\)excess\\sim Binomial\(T,\\epsilon^\{\\prime\}\)\. We believe that this assumption is reasonable given that the probability of an excess isϵ′\\epsilon^\{\\prime\}at each timestep and given no dependence of the𝐱s​o​c\\mathbf\{x\}\_\{soc\}values across time\.∎

The purpose of Algs\.[3](https://arxiv.org/html/2607.05620#alg3)and[4](https://arxiv.org/html/2607.05620#alg4)is to calculate better weights under non\-stationarity\. Because Thm\.[4\.1](https://arxiv.org/html/2607.05620#S4.Thmtheorem1)assumes weighted exchangeability, Algs\.[3](https://arxiv.org/html/2607.05620#alg3)and[4](https://arxiv.org/html/2607.05620#alg4)also satisfy the safety constraint\.

###### Corollary A\.2\.

Under weighted exchangeability, Algs\.[3](https://arxiv.org/html/2607.05620#alg3)and[4](https://arxiv.org/html/2607.05620#alg4)satisfy the safety constraint in Eq\.[7](https://arxiv.org/html/2607.05620#S3.E7)under the constraint functionq​\(𝐱\)q\(\\mathbf\{x\}\)for user\-defined violation ratesα∈\(0,1\]\\alpha\\in\(0,1\]andϵ∈\(0,1\]\\epsilon\\in\(0,1\]with probability∑k=0T⋅ϵ′\(Tk\)​\(ϵ′\)k​\(1−ϵ′\)T−k\\sum^\{T\\cdot\\epsilon^\{\\prime\}\}\_\{k=0\}\\binom\{T\}\{k\}\(\\epsilon^\{\\prime\}\)^\{k\}\(1\-\\epsilon^\{\\prime\}\)^\{T\-k\}\.

## Appendix BCoverage gaps and analysis

When the weighted exchangeability assumption is not met, we have a coverage gap:

P​\(Y​\(𝐱s​o​c,t\)≤U^ts​o​c​\(𝐱s​o​c,t\)\)≥1−ϵt\+g​a​ptP\(Y\(\\mathbf\{x\}\_\{soc,t\}\)\\leq\\hat\{U\}^\{soc\}\_\{t\}\(\\mathbf\{x\}\_\{soc,t\}\)\)\\geq 1\-\\epsilon\_\{t\}\+gap\_\{t\}\(27\)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

g​a​p≤∑i=1nwi~⋅T​V​\(R​\(𝒟c​a​l\+1\),R​\(𝒟c​a​l\+1i\)\)gap\\leq\\sum\\limits\_\{i=1\}^\{n\}\\tilde\{w\_\{i\}\}\\cdot TV\(R\(\\mathcal\{D\}\_\{cal\+1\}\),R\(\\mathcal\{D\}\_\{cal\+1\}^\{i\}\)\)\(28\)𝒟c​a​l\+1\\mathcal\{D\}\_\{cal\+1\}is the calibration set with the new test point appended to the end:\(𝒟c​a​l,𝐱t\)\(\\mathcal\{D\}\_\{cal\},\\mathbf\{x\}\_\{t\}\)\.𝒟c​a​l\+1i\\mathcal\{D\}\_\{cal\+1\}^\{i\}is𝒟c​a​l\+1\\mathcal\{D\}\_\{cal\+1\}with pointiiswapped for pointn\+1n\+1\.R​\(𝒟c​a​l\+1\)R\(\\mathcal\{D\}\_\{cal\+1\}\)is the vector of residuals for𝒟c​a​l\+1\\mathcal\{D\}\_\{cal\+1\}\. Note that the residuals and the score are not necessarily the same\. If the weights are data\-dependent,

g​a​p≤𝔼​\[∑i=1nwi~⋅T​V​\(R​\(𝒟c​a​l\+1\),R​\(𝒟c​a​l\+1i\)\|w1,…,wn\)\]gap\\leq\\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\}\)\]\(29\)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\[Barberet al\.,[2023](https://arxiv.org/html/2607.05620#bib.bib15)\]\.222Because we are doing split conformal prediction, we ignore tags\.In Eq\.[18](https://arxiv.org/html/2607.05620#S4.E18),Pc​a​lP\_\{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 becauseYYis 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 pointiiin 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 inP​\(Y\|𝐱\)P\(Y\|\\mathbf\{x\}\)orP​\(𝐱\|Y\)P\(\\mathbf\{x\}\|Y\)\.

ϵt′=ϵ\+g​a​ptP​\(Y​\(𝐱s​o​c,t\)≤U^ts​o​c​\(𝐱s​o​c,t\)\)≥1−ϵt′\\begin\{split\}&\\epsilon\_\{t\}^\{\\prime\}=\\epsilon\+gap\_\{t\}\\\\ &P\(Y\(\\mathbf\{x\}\_\{soc,t\}\)\\leq\\hat\{U\}^\{soc\}\_\{t\}\(\\mathbf\{x\}\_\{soc,t\}\)\)\\geq 1\-\\epsilon\_\{t\}^\{\\prime\}\\end\{split\}\(30\)The overall chance thatY​\(𝐱s​o​c\)Y\(\\mathbf\{x\}\_\{soc\}\)falls outside of the interval between timesteps 1 andT0T\_\{0\}is thereforeϵ′=\(∑tT0ϵt′\)/T0\\epsilon^\{\\prime\}=\(\\sum\_\{t\}^\{T\_\{0\}\}\\epsilon\_\{t\}^\{\\prime\}\)/T\_\{0\}\.ϵ′\\epsilon^\{\\prime\}could then be substituted forϵ\\epsilonin Eq\.[26](https://arxiv.org/html/2607.05620#S4.E26)\.

Algorithm 3SafeOpt with counterfactual policy constraints \(changepoint\)Input:data

𝒟\\mathcal\{D\}, prior models for

ffand

qq, initial safe set

𝒮0\\mathcal\{S\}\_\{0\}, initial decision

𝐱0\\mathbf\{x\}\_\{0\}, total iterations

TT, observation iterations

T0T\_\{0\}, target overall violation rate

α\\alpha, target counterfactual policy violation rate

ϵ\\epsilon, update rate

η\>0\\eta\>0
Split data

𝒟\\mathcal\{D\}into

𝒟t​r​a​i​n\\mathcal\{D\}\_\{train\}and

𝒟c​a​l\\mathcal\{D\}\_\{cal\}
n←l​e​n​g​t​h​\(𝒟c​a​l\)n\\leftarrow length\(\\mathcal\{D\}\_\{cal\}\)

Train estimator

Y^\\hat\{Y\}on

𝒟t​r​a​i​n\\mathcal\{D\}\_\{train\}
ifthere may be a covariate shift between the calibration and test datathen

Observe test data for

T0T\_\{0\}timesteps

endif

fori=0 to ndo

Form weights

wiw\_\{i\}according to Eq\.[18](https://arxiv.org/html/2607.05620#S4.E18)or Eq\.[19](https://arxiv.org/html/2607.05620#S4.E19)

endfor

Calculate

∑j=1nwj\\sum\_\{j=1\}^\{n\}w\_\{j\}
Calculate the scores on

𝒟c​a​l\\mathcal\{D\}\_\{cal\}as in Eq\.[15](https://arxiv.org/html/2607.05620#S4.E15)or Eq\.[23](https://arxiv.org/html/2607.05620#S4.E23)

Calculate violation rate

α′\\alpha^\{\\prime\}according to Eq\.[26](https://arxiv.org/html/2607.05620#S4.E26)

for

t=0t=0to

TTdo

Observe

f​\(𝐱t\)f\(\\mathbf\{x\}\_\{t\}\)and

Y​\(𝐱𝐭\)Y\(\\mathbf\{x\_\{t\}\}\)
if

𝐱t\\mathbf\{x\}\_\{t\}is the same as the standard\(s\) of care or

𝐱t\\mathbf\{x\}\_\{t\}is the decision under the safe fallback policythen

q​\(𝐱𝐭\)=ωq\(\\mathbf\{x\_\{t\}\}\)=\\omega

else

Form weight

wtw\_\{t\}for current

𝐱t\\mathbf\{x\}\_\{t\}
Compute the normalized weight

w~i\\tilde\{w\}\_\{i\}as in Eq\.[17](https://arxiv.org/html/2607.05620#S4.E17)

Compute the normalized weight

w~t\\tilde\{w\}\_\{t\}as in Eq\.[17](https://arxiv.org/html/2607.05620#S4.E17)

Form the conformal upper or lower bounds according to Eq\.[16](https://arxiv.org/html/2607.05620#S4.E16)or Eq\.[21](https://arxiv.org/html/2607.05620#S4.E21)

Form

q​\(𝐱𝐭\)q\(\\mathbf\{x\_\{t\}\}\)according to Eq\.[13](https://arxiv.org/html/2607.05620#S4.E13),[20](https://arxiv.org/html/2607.05620#S4.E20), or[24](https://arxiv.org/html/2607.05620#S4.E24)

endif

ifa changepoint is detectedthen

Observe test data for

T0T\_\{0\}timesteps

fori=0 to ndo

Form weights

wiw\_\{i\}according to Eq\.[19](https://arxiv.org/html/2607.05620#S4.E19)

endfor

Calculate

∑j=1nwj\\sum\_\{j=1\}^\{n\}w\_\{j\}
endif

𝐱t\+1←\\mathbf\{x\}\_\{t\+1\}\\leftarrowSafeOpt\-Helper\(\)

endfor

Algorithm 4SafeOpt with counterfactual policy constraints \(drift\)Input:data

𝒟\\mathcal\{D\}, prior models for

ffand

qq, initial safe set

𝒮0\\mathcal\{S\}\_\{0\}, initial decision

𝐱0\\mathbf\{x\}\_\{0\}, total iterations

TT, observation iterations

T0T\_\{0\}, target overall violation rate

α\\alpha, target counterfactual policy violation rate

ϵ\\epsilon, update rate

η\>0\\eta\>0
Split data

𝒟\\mathcal\{D\}into

𝒟t​r​a​i​n\\mathcal\{D\}\_\{train\}and

𝒟c​a​l\\mathcal\{D\}\_\{cal\}
n←l​e​n​g​t​h​\(𝒟c​a​l\)n\\leftarrow length\(\\mathcal\{D\}\_\{cal\}\)

Train estimator

Y^\\hat\{Y\}on

𝒟t​r​a​i​n\\mathcal\{D\}\_\{train\}
ifthere may be a covariate shift between the calibration and test datathen

Observe test data for

T0T\_\{0\}timesteps

endif

fori=0 to ndo

Form time\-decaying weights

wiw\_\{i\}\(for example, according to Eq\.[31](https://arxiv.org/html/2607.05620#A3.E31)\)

endfor

Calculate

∑j=1nwj\\sum\_\{j=1\}^\{n\}w\_\{j\}
Calculate the scores on

𝒟c​a​l\\mathcal\{D\}\_\{cal\}as in Eq\.[15](https://arxiv.org/html/2607.05620#S4.E15)or Eq\.[23](https://arxiv.org/html/2607.05620#S4.E23)

Calculate violation rate

α′\\alpha^\{\\prime\}according to Eq\.[26](https://arxiv.org/html/2607.05620#S4.E26)

for

t=0t=0to

TTdo

if

t%2==0t\\%2==0then

Observe

f​\(𝐱t\)f\(\\mathbf\{x\}\_\{t\}\)and

Y​\(𝐱𝐭\)Y\(\\mathbf\{x\_\{t\}\}\)
if

𝐱t\\mathbf\{x\}\_\{t\}is the same as the standard\(s\) of care or

𝐱t\\mathbf\{x\}\_\{t\}is the decision under the safe fallback policythen

q​\(𝐱𝐭\)=ωq\(\\mathbf\{x\_\{t\}\}\)=\\omega

else

Form weight

wtw\_\{t\}for current

𝐱t\\mathbf\{x\}\_\{t\}
Compute the normalized weight

w~i\\tilde\{w\}\_\{i\}as in Eq\.[17](https://arxiv.org/html/2607.05620#S4.E17)

Compute the normalized weight

w~t\\tilde\{w\}\_\{t\}as in Eq\.[17](https://arxiv.org/html/2607.05620#S4.E17)

Form the conformal upper or lower bounds according to Eq\.[16](https://arxiv.org/html/2607.05620#S4.E16)or Eq\.[21](https://arxiv.org/html/2607.05620#S4.E21)

Form

q​\(𝐱𝐭\)q\(\\mathbf\{x\_\{t\}\}\)according to Eq\.[13](https://arxiv.org/html/2607.05620#S4.E13),[20](https://arxiv.org/html/2607.05620#S4.E20), or[24](https://arxiv.org/html/2607.05620#S4.E24)

endif

𝐱t\+1←πs​o​c\\mathbf\{x\}\_\{t\+1\}\\leftarrow\\pi\_\{soc\}

Add

𝐱t\+1\\mathbf\{x\}\_\{t\+1\}to

𝒟c​a​l\\mathcal\{D\}\_\{cal\}
else

𝐱t\+1←\\mathbf\{x\}\_\{t\+1\}\\leftarrowSafeOpt\-Helper\(\)

endif

endfor

## Appendix CNon\-stationarity extension

SafeOpt with online conformal prediction\[Zhanget al\.,[2024](https://arxiv.org/html/2607.05620#bib.bib13)\]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\.[B](https://arxiv.org/html/2607.05620#A2)\)\.

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[Barberet al\.](https://arxiv.org/html/2607.05620#bib.bib15)\[[2023](https://arxiv.org/html/2607.05620#bib.bib15)\]for a discussion of nonexchangeability, including nonexchangeability from non\-stationarity, such as changepoints and drift\.

### C\.1Motivating 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\.

### C\.2Algorithmic description

If covariate shift is detected at a changepoint \(Scenario 3\), we can observe unlabeled data and recalculate the weights as in Eq\.[19](https://arxiv.org/html/2607.05620#S4.E19)\. After the weights have been calculated, we can create conformal predictions forYs​o​cY\_\{soc\}again and SafeOpt continues \(Alg\.[3](https://arxiv.org/html/2607.05620#alg3)\)\.

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 functionggsuch that

wi=exp⁡\(−g​\(t−ti\)\)⋅Pt​e​s​t​\(𝐱i\)/Pc​a​l​\(𝐱i\)w\_\{i\}=\\exp\(\-g\(t\-t\_\{i\}\)\)\\cdot P\_\{test\}\(\\mathbf\{x\}\_\{i\}\)/P\_\{cal\}\(\\mathbf\{x\}\_\{i\}\)\(31\)wheretit\_\{i\}is the time of data collection for pointii, andttis 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πs​o​c\\pi\_\{soc\}\. Thisπs​o​c\\pi\_\{soc\}step is added to our calibration set \(Alg\.[4](https://arxiv.org/html/2607.05620#alg4)\)\.

## Appendix DAdditional experimental details

Part of our code is based onBerkenkampet al\.\[[2023](https://arxiv.org/html/2607.05620#bib.bib3)\], which uses an MIT License\. The MovieLens 100k dataset\[Harper and Konstan,[2015](https://arxiv.org/html/2607.05620#bib.bib17)\]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\[Kanget al\.,[2019](https://arxiv.org/html/2607.05620#bib.bib18)\]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η\\etawas set to 0\.2 for all experiments\. In preliminary tests, we found that this was an appropriate value for 100 timesteps\. We choseα=0\.1\\alpha=0\.1since that value was used in many of theZhanget al\.\[[2024](https://arxiv.org/html/2607.05620#bib.bib13)\]experiments, and it was a high enough value to accommodate a couple of estimated components\.ϵ=0\.01\\epsilon=0\.01since error rates in conformal prediction are usually set to some small value, andϵ\\epsilonneeded to be small enough thatα′\\alpha^\{\\prime\}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 mostn\+1n\+1data points, wherennis the size of the calibration set\.

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

### D\.1MovieLens 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𝐱0\\mathbf\{x\}\_\{0\}from the standard\-of\-care policy and its corresponding observations off​\(𝐱0\)f\(\\mathbf\{x\}\_\{0\}\)andq​\(𝐱0\)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:

k​\(i​t​e​m,i​t​e​m′,i​t​e​ms​o​c,i​t​e​ms​o​c′,u​s​e​r,u​s​e​r′\)=k​\(i​t​e​m,i​t​e​m′\)⋅k​\(i​t​e​ms​o​c,i​t​e​ms​o​c′\)\+k​\(i​t​e​m,i​t​e​m′\)⋅k​\(u​s​e​r,u​s​e​r′\)\\begin\{split\}&k\(item,item^\{\\prime\},item\_\{soc\},item\_\{soc\}^\{\\prime\},user,user^\{\\prime\}\)=\\\\ &k\(item,item^\{\\prime\}\)\\cdot k\(item\_\{soc\},item\_\{soc\}^\{\\prime\}\)\+k\(item,item^\{\\prime\}\)\\cdot k\(user,user^\{\\prime\}\)\\end\{split\}\(32\)wherek​\(a,a′\)=a⊤​a′k\(a,a^\{\\prime\}\)=a^\{\\top\}a^\{\\prime\}\.

### D\.2Chemical 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:

k​\(i​t​e​m,i​t​e​m′,i​t​e​ms​o​c,i​t​e​ms​o​c′\)=k​\(i​t​e​m,i​t​e​m′\)⋅k​\(i​t​e​ms​o​c,i​t​e​ms​o​c′\)\\begin\{split\}&k\(item,item^\{\\prime\},item\_\{soc\},item\_\{soc\}^\{\\prime\}\)=k\(item,item^\{\\prime\}\)\\cdot k\(item\_\{soc\},item\_\{soc\}^\{\\prime\}\)\\end\{split\}\(33\)wherek​\(a,a′\)=a⊤​a′k\(a,a^\{\\prime\}\)=a^\{\\top\}a^\{\\prime\}\.

### D\.3Synthetic experiments

To create the synthetic functionsZZandYY, we sample functions from two Gaussian processes\. The process forZZhas an RBF kernel with lengthscale 1\.0 and variance 1\.0\. The process forYYhas an RBF kernel variance 1\.0\. There is a single context variableccthat affectsYY\. The lengthscale is 0\.1 fori​t​e​mitem, the choice that SafeOpt will make, and 1\.0 forcc\.

Z∼G​P​\(0,R​B​F​\(1\.0,1\.0\)\)Y∼G​P​\(0,R​B​F​\(1\.2,\(1\.0,0\.1\)\)\)\+5⋅c\\begin\{split\}Z&\\sim GP\(0,RBF\(1\.0,1\.0\)\)\\\\ Y&\\sim GP\(0,RBF\(1\.2,\(1\.0,0\.1\)\)\)\+5\\cdot c\\end\{split\}\(34\)We seed the Gaussian processes with support points between 0 and 1 fori​t​e​mitemand 0 and 5 forcc; 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 forqqis similar to the MovieLens kernel, exceptccis our 1D context variable:

k​\(i​t​e​m,i​t​e​m′,i​t​e​ms​o​c,i​t​e​ms​o​c′,c,c′\)=k​\(i​t​e​m,i​t​e​m′\)⋅k​\(i​t​e​ms​o​c,i​t​e​ms​o​c′\)\+k​\(i​t​e​m,i​t​e​m′\)⋅k​\(c,c′\)\\begin\{split\}&k\(item,item^\{\\prime\},item\_\{soc\},item\_\{soc\}^\{\\prime\},c,c^\{\\prime\}\)=\\\\ &k\(item,item^\{\\prime\}\)\\cdot k\(item\_\{soc\},item\_\{soc\}^\{\\prime\}\)\+k\(item,item^\{\\prime\}\)\\cdot k\(c,c^\{\\prime\}\)\\end\{split\}\(35\)wherek​\(a,a′\)=a⊤​a′k\(a,a^\{\\prime\}\)=a^\{\\top\}a^\{\\prime\}\.

## Appendix EAdditional results

Note thatα=0\.1\\alpha=0\.1andϵ=0\.01\\epsilon=0\.01for all plots in this section\. We create confidence bounds based on 10 random seeds\.

### E\.1Additional MovieLens experiments

Fig\.[1](https://arxiv.org/html/2607.05620#S6.F1)shows the results when the training and calibration data come from the standard\-of\-care policy or the data collection policy whenq​\(𝐱\)q\(\\mathbf\{x\}\)is defined as in Eq\.[13](https://arxiv.org/html/2607.05620#S4.E13)\. Fig\.[2](https://arxiv.org/html/2607.05620#A5.F2)shows the full results for both Eq\.[13](https://arxiv.org/html/2607.05620#S4.E13)and Eq\.[24](https://arxiv.org/html/2607.05620#S4.E24)under the standard\-of\-care policy, and Fig\.[3](https://arxiv.org/html/2607.05620#A5.F3)shows the full results for both Eq\.[13](https://arxiv.org/html/2607.05620#S4.E13)and Eq\.[24](https://arxiv.org/html/2607.05620#S4.E24)under the data collection policy\. Note that some results are repeated from Fig\.[1](https://arxiv.org/html/2607.05620#S6.F1)\. The results do not differ much from the results in the main text\.

![Refer to caption](https://arxiv.org/html/2607.05620v1/x2.png)Figure 2:MovieLens results with standard\-of\-care policy training and calibration\.![Refer to caption](https://arxiv.org/html/2607.05620v1/x3.png)Figure 3:MovieLens results with data collection policy training and calibration\.
### E\.2Additional chemical reaction experiments

Fig\.[1](https://arxiv.org/html/2607.05620#S6.F1)shows the chemical reaction results under a random policy for training and calibration whenq​\(𝐱\)q\(\\mathbf\{x\}\)is defined as in Eq\.[13](https://arxiv.org/html/2607.05620#S4.E13)and Eq\.[24](https://arxiv.org/html/2607.05620#S4.E24)\. Fig\.[4](https://arxiv.org/html/2607.05620#A5.F4)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\.

![Refer to caption](https://arxiv.org/html/2607.05620v1/x4.png)Figure 4:Reaction results with standard\-of\-care policy training and calibration\.
### E\.3Synthetic 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 useq​\(𝐱\)q\(\\mathbf\{x\}\)from Eq\.[13](https://arxiv.org/html/2607.05620#S4.E13)and use training and calibration data from the standard\-of\-care policy, Fig\.[5](https://arxiv.org/html/2607.05620#A5.F5)also uses these conditions\. Note that for each seed, we randomly generate a newZZandYY\.

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 maximumffvalue\.

![Refer to caption](https://arxiv.org/html/2607.05620v1/x5.png)Figure 5:Synthetic results\.
### E\.4Non\-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γΔ​t\\gamma^\{\\Delta t\}, whereγ\\gammais a parameter plotted in Fig\.[6](https://arxiv.org/html/2607.05620#A5.F6)\.

The left\-hand plot of Fig\.[6](https://arxiv.org/html/2607.05620#A5.F6)shows the violation rate when we correctly re\-weight for drift according toγΔ​t\\gamma^\{\\Delta t\}\. Recent observations are weighted so that they are more relevant, and past observations are less relevant\. As described in Sec\.[C](https://arxiv.org/html/2607.05620#A3), 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](https://arxiv.org/html/2607.05620#alg1)as normal\.

![Refer to caption](https://arxiv.org/html/2607.05620v1/x6.png)Figure 6:Drift results\.In Fig\.[7](https://arxiv.org/html/2607.05620#A5.F7)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\.

![Refer to caption](https://arxiv.org/html/2607.05620v1/x7.png)Figure 7:Changepoint results\.
### E\.5Sensitivity 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\.[8](https://arxiv.org/html/2607.05620#A5.F8)\)\. For these plots,ω=0\.1\\omega=0\.1since the number of fallback steps is low under default conditions, as shown in Fig\.[5](https://arxiv.org/html/2607.05620#A5.F5)\.

![Refer to caption](https://arxiv.org/html/2607.05620v1/x8.png)Figure 8:Sensitivity analysis: estimator and variability\.When the weights are mis\-specified, conformal coverage does not hold\. In Fig\.[9](https://arxiv.org/html/2607.05620#A5.F9)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,ω=0\.001\\omega=0\.001since small changes to the weights have a large impact at lowω\\omegavalues\.

![Refer to caption](https://arxiv.org/html/2607.05620v1/x9.png)Figure 9:Sensitivity analysis: weights\.Note that differentα\\alphavalues will change the maximum number of fallback steps and allowed violation rate threshold\. Highω\\omegavalues result in fewer fallback steps in general and therefore less sensitive to noise or worseY^\\hat\{Y\}estimators, while lowω\\omegavalues result in more fallback steps\.