Safe Bayesian Optimization for the Control of High-Dimensional Embodied Systems

In plain words

Tuning a controller means trying settings and watching what happens. Some settings can damage the system, and a method that tries many settings at once can waste trials or try unsafe ones. The paper proposes HdSafeBO, a way to tune controllers that have hundreds to thousands of settings while avoiding settings that are probably unsafe.

It has two main ideas. First, it compresses the settings into a smaller space, so the search works with fewer numbers. Second, it searches only a small neighbourhood around the best safe setting found so far. It accepts a candidate only when a cautious estimate says the candidate is probably safe. The neighbourhood grows after a run of good safe trials and shrinks after failures.

The authors test it on synthetic functions, on a simulated arm that holds a bottle, and on electrical stimulation of the spinal cord. The spinal test was run in simulation and in a clinical experiment with one patient.

TL;DR

• Safe black-box optimisation for controllers with up to thousands of parameters, under a probabilistic safety constraint §1 §4.
• Search runs in a latent space given by PCA or an isometric VAE (IRVAE) §4.3 app. B.1 app. B.2. A trust region around the best safe point is searched with Thompson sampling §4 alg. 1.
• Author-reported: 3.96 objective and 27.42 cumulative violation on 1000-D GP synthetic functions, against 3.07 and 36.59 for the best baseline (LineBO) tab. 1.
• Relevant if you tune controllers under safety limits and the parameter count is large. Read the Pitfalls first: the guarantee needs an unstated probability threshold, and the released code departs from the algorithm in its trust-region update.

Problem

Online tuning of a robot or neuromuscular controller is a sequential decision problem: maximise an unknown utility f(x) while an unknown safety function g(x) must stay safe during every trial §1 §3. The paper defines a point as safe when g(x) > 0 §3. Existing safe Bayesian optimisation methods, such as SafeOpt, are conservative: they become inefficient in high-dimensional settings, and the paper says they may struggle to expand the safe region once the parameter count exceeds about 10 §1 §2. Constrained BO methods handle high dimension but do not enforce safety during the search §1 §2. The paper targets inputs from dozens to thousands of variables §3.

Method

The method has four parts: Gaussian-process surrogates for f and g §3; an isometric embedding from the input space to a low-dimensional latent space §4.3; an optimistic safe set inside a trust region §4.1 §4.2; and a safety-aware update of the trust region §4.2 alg. 2.

1. Objective. Maximise f(x_t) subject to g(x_t) ≥ 0 for all t eq. 1. In Section 4 the requirement is relaxed to a probabilistic one, Pr(g(x_t) ≥ 0) ≥ α for every t §4.
2. Surrogates. Fit a GP to f and another to g, using the standard posterior mean and variance eq. 2. The upper confidence bound is u_t(x) = μ_{t-1}(x) + β_t σ_{t-1}(x) §3.
3. Embedding. An isometric encoder Π maps the data to the latent space Z, and Π⁻¹ maps latent points back to the input space alg. 1, line 4 §4.3. PCA is used for the synthetic and musculoskeletal tasks §5; an IRVAE is used for the neural-stimulation task §5.3 app. B.2.
4. Trust region. Around the best safe latent point x*, form a box L_t with side length l_t alg. 1, line 6 §4.2.
5. Optimistic safe set. Keep the latent points in L_t whose optimistic safety bound passes the threshold: μ_{t-1} + β_t σ_{t-1} ≥ h alg. 1, line 7. Algorithm 1 takes h as an input, but the line itself is written with 0 as the pass value alg. 1, line 7 alg. 1 input.
6. Choice of β_t. Choose β_t so that Φ(β_t) ≤ 1 − α, where Φ is the standard normal CDF. This makes the one-step probability bound hold for every x [§4.1 Prop. 4.1] eq. 3.
7. Acquisition. Thompson sampling: draw one sample of the objective posterior over the safe set and take its maximum alg. 1, line 8 §4. The sampling makes batch selection possible §4.
8. Evaluation and update. Decode z_t to x_t alg. 1, line 9, evaluate f and g alg. 1, lines 10-11, append the data alg. 1, line 12, then update the incumbent and trust region alg. 1, line 13 alg. 2.
9. Trust-region update (SafeUpdate). A trial is a success when it is safe and improves the incumbent. Any other trial is a failure §4.2. After enough successes the side length doubles, capped at l_max. After enough failures it halves, floored at l_min. When it reaches l_min the length restarts at l_0, and all past samples are kept §4.2 alg. 2.

Guarantees, stated in the paper: Proposition 4.1 gives the one-step probability bound eq. 3. Theorem 4.2 bounds the expected cumulative safety violation E[V_T] by (1 − α)·sqrt(C₁ T β^c_T γ_T) eq. 4 app. A.2. Proposition 4.3 says that with a stationary kernel and an exact isometry, the probability bound transfers from the latent space to the original space §4.3 app. A.3.

The algorithm, condensed. Each line of the paper's box is cited, and the safety-aware trust-region update is in its own block:

# HdSafeBO main loop  [alg. 1]
# Input:  data (inputs x, objective f, safety g), threshold h,
#         encoder Pi, decoder Pi_inv, budget T
# Output: data extended by T evaluations

def hdsafebo(data, threshold, encoder, decoder, budget):
    # Initial incumbent: best observed point with g above h  [alg. 1, line 2]
    incumbent = best_safe_point(data, threshold)
    for step in range(budget):
        # Map the observed inputs into the latent space  [alg. 1, line 4]
        latents = encoder(data.inputs)
        # Fit one GP to f and one to g on the latent data  [alg. 1, line 5]
        f_gp, g_gp = fit_gaussian_processes(latents, data.f, data.g)
        # Trust region around the incumbent's latent code  [alg. 1, line 6]
        box = trust_region(encoder(incumbent.x), trust.length)
        # Keep candidates with optimistic safety above h  [alg. 1, line 7]
        safe = optimistic_safe_subset(box, g_gp, threshold, beta)
        # Thompson sample of f; take the best safe candidate  [alg. 1, line 8]
        latent_next = argmax(safe, sample_posterior(f_gp, safe))
        # Map the chosen latent point back to the input space  [alg. 1, line 9]
        x_next = decoder(latent_next)
        # Evaluate both objective and safety  [alg. 1, lines 10-11]
        f_next, g_next = evaluate(x_next)
        # Append the new observation  [alg. 1, line 12]
        data.append(x_next, f_next, g_next)
        # Update incumbent and trust region  [alg. 1, line 13] [alg. 2]
        incumbent, trust = safe_update(
            Sample(x_next, f_next, g_next), incumbent, trust, settings)

# SafeUpdate: incumbent and trust-region update  [alg. 2] [§4.2]
# Input:  new sample (x, f, g), incumbent, counters, length, settings
# Output: incumbent, counters, length

def safe_update(sample, incumbent, state, settings):
    # Success: safe (g > 0, as in the paper) and better than the incumbent
    if sample.g > 0 and sample.f > incumbent.f:
        incumbent = sample
        state.successes += 1
        state.failures = 0
    # Otherwise: a failure, whether unsafe or without improvement  [§4.2]
    else:
        state.successes = 0
        state.failures += 1
    # Expand after enough successes, capped at l_max  [alg. 2] [§4.2]
    if state.successes == settings.success_tolerance:
        state.length = min(2 * state.length, settings.length_max)
        state.successes = 0
    # Shrink after enough failures, floored at l_min  [alg. 2] [§4.2]
    elif state.failures == settings.failure_tolerance:
        state.length = max(state.length / 2, settings.length_min)
        state.failures = 0
        # At the floor, restart at l_0; all past samples are kept  [§4.2]
        if state.length == settings.length_min:
            state.length = settings.length_0
    return incumbent, state

Hyperparameters

The paper gives values for β, the batch size, the synthetic-function setup and the neural-stimulation autoencoder. It does not give the probability threshold α, the trust-region constants, the Matérn smoothness or the number of candidate points. Rows marked "code only" are from the released repository, not from the paper. Where code and paper differ, the code value is marked.

Musculoskeletal reward, as the paper writes it app. B.5: the per-step reward is the sum of a pose term, a bonus, a grasp term, a survival term and an activation penalty, minus a penalty for leaving the range and a large drop penalty. The reward itself is computed by an external environment; see Paper vs code.

Neural stimulation selectivity index, simulation app. B.6: the index for muscle i is its normalised activation minus the mean of its neighbours' activations. The clinical index app. B.7 is the normalised peak-to-peak EMG divided by one plus the sum of the other target muscles' values.

Results

Author-reported. Setup: synthetic 1000-D GP functions, 500 evaluations, 100 runs tab. 1; musculoskeletal control with 50 runs fig. 2; neural stimulation in simulation with 10 runs fig. 3; one patient in the clinic §5.3.

Synthetic GP functions tab. 1. Objective higher is better. Violation lower is better.

The paper says SafeOpt failed to expand its safe region from the initial points and is not in the table §5.1. The baselines were run in the PCA subspace with d = 50 §5.1.

Embedding comparison, synthetic functions tab. 3. Columns marked L run in the latent space; O run in the original space.

Musculoskeletal control fig. 2. Results are only in the figure. The text says all baselines fail to improve the objective in the original 3,575-D space, so baselines were run in the PCA subspace; SCBO (L) and LineBO (L) never reached a positive reward and are omitted §5.2. The text reports that HdSafeBO samples fewer unsafe parameters than the baselines in the subspace §5.2. No table values are given.

Ablation, musculoskeletal task fig. 6. Without local search the algorithm over-explores. Without optimistic safety identification it makes more unsafe choices early on app. C.2. Figure only; no table values.

Confidence scalar β fig. 8. The paper runs β = 0, 2, 4, 8, 16. Safety and violation get slightly worse as β grows, and the objective is similar across the range app. C.4. Figure only.

Neural stimulation, simulation fig. 3. HdSafeBO gives the best objective, the highest safe-selection ratio and the lowest violation on the ST and GA selectivity tasks §5.3. Values are in the figure only.

Neural stimulation, one patient §5.3 app. B.7. The text gives a selectivity improvement for 7 of 8 target muscles. It lists eight values (left IL 0.112, left RF 0.143, left TA 0.097, left BF 0.380, right IL 0.00, right RF 0.266, right TA 0.216, right BF 0.141) without saying whether each is the baseline or the optimised value. Only three recommended configurations were rated unsafe, and they caused large leg movements but no injury §5.3.

Comparison with random-embedding BO tab. 5. Best objective found without safety constraints, neural stimulation task.

IRVAE versus conventional VAE tab. 4. Distance preservation on 10,000 simulated points, neural stimulation task.

Constrained digit generation tab. 6. Normalised objective (higher is better) and cumulative violation (lower is better), averaged over ten digit tasks. Budget 200 with 20 initial images.

The text says violations are 63% lower than the second-best method, SCBO app. C.6.

Run time per iteration, neural stimulation tab. 7. Mean ± 1 std over 80 iterations, seconds: HdSafeBO 11.75 ± 1.35; LineBO 9.88 ± 1.01; SCBO 11.99 ± 1.21; CONFIG 9.49 ± 1.05; cEI 20.38 ± 1.21; CMA-ES 9.59 ± 0.72 app. C.7. The authors say the run time is small next to a trial that takes one to two minutes app. C.7.

Paper vs code

The official repository is yunyuewei/HdSafeBO at commit 3c4eaf167c0a410f1a0f937119120b8b6f0a004d. It is linked from the project page that the paper cites in its introduction §1. Its README says it is "the code repository of the paper". The owner is the first author README.md. The repository has no LICENSE file, and the GitHub API reports no license, so we cite it by path and line only and quote nothing from it. The clone is at raw/code/2412.20350/.

• Trust-region success and failure rules. The paper says a success is a safe sample that improves the incumbent, and any unsafe or non-improving sample is a failure §4.2 alg. 2. The code counts a batch as a success only if every sample in the batch is safe and the best of them improves on the incumbent hdsafebo.py. A safe batch with no improvement counts as neither success nor failure hdsafebo.py.
• Success never fires in the released code. The improvement test uses a relative margin, best + 1e-3·|best|, and best_value starts at −∞ [code: algorithm/hdsafebo.py#L21] hdsafebo.py. Then the margin is −∞ + ∞, which is NaN, and a comparison with NaN is false. best_value is never set anywhere else in the repository hdsafebo.py. So as written, the success counter never increments and the trust region never expands. It only shrinks on failures and restarts at 0.8 when it falls below the minimum. I did not run this code, because executing Python was not approved in this session. The NaN reasoning is from IEEE arithmetic. TuRBO in the same repository seeds best_value from the data [code: algorithm/turbo.py#L98], and HdSafeBO does not.
• Trust-region constants. The paper does not give them app. B.3. The code uses initial length 0.8, maximum 1.6 and minimum 0.5⁷ [code: algorithm/hdsafebo.py#L14-L16]. It sets the success tolerance to 10, with a comment that the original paper uses 3 hdsafebo.py. The failure tolerance is ⌈max(4/B, d/B)⌉ for batch size B and dimension d [code: algorithm/hdsafebo.py#L25-L27].
• Restart value. The paper resets to l_0 §4.2. The code hard-codes 0.8 and sets a global_samp flag that only changes log messages [code: algorithm/hdsafebo.py#L203-L209].
• Safety bound scale. The paper's bound is μ + β·σ §3 alg. 1, line 7. The code uses confidence_region() and sets var = (upper − lower)/2, then ucb = mean + safe_beta · var hdsafebo.py. If gpytorch's confidence_region() returns mean ± 2σ, which is my recollection of the library and not checked here, the effective multiplier is 2β = 4 standard deviations. Check the installed gpytorch version before relying on either reading.
• Incumbent when nothing is safe. Algorithm 1 initialises x* with the best point satisfying g > 0 alg. 1, line 2. When no safe point exists, the code uses the best observed point regardless of safety [code: algorithm/hdsafebo.py#L109-L112].
• Fallback to unrestricted sampling. When the trust region is too small to hold a safe candidate, the code halves it until its length reaches the minimum and then sets the safety mask to all True hdsafebo.py. Thompson sampling then picks from every candidate in that tiny box, including unsafe ones. The flag is named global_samp, but the box is not global. The paper does not describe this case.
• Candidate count. Each iteration draws 5,000 scrambled Sobol points inside the trust region hdsafebo.py [code: algorithm/hdsafebo.py#L132-L134]. The paper does not state this number.
• Threshold convention. The paper's SafeUpdate uses g > 0 alg. 2. The candidate filter uses the threshold h alg. 1, line 7. The code uses >= for the batch check against the threshold hdsafebo.py. Experiments use h = 0.05, −0.75 and −4.4 in different units app. B.4–B.6.
• GP fitting. The paper says kernel parameters are refit each iteration app. B.3. The code uses a Matérn-2.5 ARD kernel with length scales in [0.005, 4.0], noise in [1e-8, 1e-3] and standardised objective values [code: algorithm/base_optimizer.py#L61-L66]. The safety GP is not standardised [code: algorithm/base_optimizer.py#L131]. The noise bounds and length-scale bounds are code-only choices.
• Synthetic function construction. The paper says functions are samples from GPs with random linear projections and d_e randomly selected variables [app. B.4]. The code samples each new query from the GP conditioned on earlier queries, and reuses the stored value for exact repeats [code: task/gpfun/gp_fun.py#L72-L93]. The effective variables are drawn with replacement, so duplicates are possible gpfun_task.py. The PCA is fitted to the 200 initial points gpfun_task.py. The noise argument of 0.01 is stored and never used [code: task/gpfun_task.py#L24] gpfun_task.py.
• Musculoskeletal safety and latent bounds. The paper defines safety as the landing speed of the bottle, with threshold 4.4 app. B.5. The code computes the value as the negative end-of-episode speed, clipped at −10 [code: task/muscle/task_env_mjhand.py#L208-L209], with threshold −4.4 [code: task/muscle_task.py#L20]. The paper does not mention the clip. The code clips latent weights to ±3 task_env_mjhand.py [code: task/muscle/task_env_mjhand.py#L176] muscle_task.py. The paper gives no bound. The task environment is loaded from an external package (msmodel_gym/HandFloat-v1) task_env_mjhand.py, which is not in this repository. The reward formula therefore cannot be checked here. The muscle dimensions do match the paper: 55 × 65 = 3,575 and 5 × 65 = 325 optimization_muscle.py optimization_muscle.py. The driver runs 50 rounds of 1,000 evaluations with batch 10 [code: optimization_muscle.py#L27-L28] optimization_muscle.py.
• Neural-stimulation oracle. The paper says the nearest unevaluated point is found by distance on the electric-field map app. B.6. The evaluation path finds the nearest pole vector in the 17-D input space [code: task/scs_task.py#L282-L287]. Only the latent-to-input decoder uses map distance scs_task.py. A zero current returns utility 0 and safety 1 without a lookup scs_task.py. The paper says nothing about this case.
• Neural-stimulation targets. The paper reports two targets, ST and GA §5.3. The driver runs only SI_idx = 4 optimization_scs.py. The code does not say which muscle index 4 is.
• Autoencoder. The paper's table lists 4 layers of 256 units and latent dimension 16, trained for 300 epochs with Adam at 1e-4 tab. 2 app. B.2. The released config matches these [code: ae_model/SCS/DIM_16/irvae_scs/scs_irvae_z16.yml#L35] [code: ae_model/SCS/DIM_16/irvae_scs/scs_irvae_z16.yml#L54] [code: ae_model/SCS/DIM_16/irvae_scs/scs_irvae_z16.yml#L64] [code: ae_model/SCS/DIM_16/irvae_scs/scs_irvae_z16.yml#L67] [code: ae_model/SCS/DIM_16/irvae_scs/scs_irvae_z16.yml#L69]. The isometry loss has weight 1000 scs_irvae_z16.yml and uses eta = 0.2 ae.py. The paper states neither. I did not trace which model class irvae_scs maps to, so the loss citation is for the first VAE class in that file.
• Digit task. The code matches the paper on 20 initial images, budget 200, latent dimension 6, batch 1 and softmax temperature 200 [code: task/digit_task.py#L39] optimization_digital.py [code: optimization_digital.py#L28] optimization_digital.py [code: task/digit_task.py#L158]. The code uses a safety threshold of 0.5 [code: task/digit_task.py#L31], which the paper does not state. The objective is the pixel sum digit_task.py, but I found no normalisation by the best MNIST point, which the paper describes tab. 6.
• Synthetic table versus table. Table 1 and Table 3 should share the same runs. For LineBO, CONFIG, cEI and CMA-ES they match the O columns. For SCBO they do not: Table 1 gives 2.97 and 39.05, while Table 3 gives 2.95 and 38.47 in the O column and 3.15 and 38.51 in the L column tab. 1 tab. 3.
• Clinical experiment. No clinical code or data is in the repository. The README covers only simulation and says the spinal data must be downloaded separately README.md.
• Wrong driver in the README. The spinal-cord section runs optimization_muscle.py README.md. The spinal driver is optimization_scs.py optimization_scs.py.
• Musculoskeletal model not released. The paper says the model will be released and is available on request app. B. It is not in the repository.

Pitfalls

• The β = 2 setting and the probability threshold conflict. Proposition 4.1 needs Φ(β) ≤ 1 − α. With β = 2, the normal tail 1 − Φ(2) is about 0.023, so the bound holds only for α ≤ 0.023 §4.1 eq. 3. The paper never states α, so the probability level of the experiments is unknown. Pick α, then check Φ(β) ≤ 1 − α, and say which α you used.
• Theorem 4.2 uses |X|. The constant β^c_T = 2 log(|X| T² π²/(6δ)) needs the cardinality of the domain eq. 4 app. A.2. The domain is continuous, and the proof does not discuss discretisation. The bound is therefore not directly applicable as written.
• The proof switches between β_t and β^c_t. One step uses the confidence scalar from Proposition 4.1, β_t app. A.2 (eq. "beta"), and the next uses the Srinivas-style bound with β^c_t app. A.2 (eq. "cb"). The paper does not reconcile the two. The expectation bound is also stated "with probability at least 1 − δ" eq. 4, which mixes an expected and a high-probability claim.
• Isometry is only approximate. Proposition 4.3 needs exact distance preservation app. A.3. The IRVAE reaches a linear distance correlation of 0.9729 tab. 4, so the guarantee does not strictly hold for the neural task. The paper says the algorithm still works in practice app. C.3. A PCA decoder is an affine isometry onto its span, but it reaches only that span. Policies outside the principal subspace are not reachable.
• Trust-region constants decide the behaviour. The paper gives none of τ_s, τ_f, l_0, l_max or l_min app. B.3. The released defaults differ from the algorithm, and the success rule never fires (see Paper vs code). Expect the trust region to shrink and restart in the released code. Do not assume the paper's expansion behaviour.
• Safety sign and units. Section 3 defines a safe point as g > 0. The experiments use three different thresholds: 0.05 for the spinal task with g = 1 − max activation app. B.6, −0.75 for the synthetic task app. B.4, and −4.4 for the musculoskeletal task with negative landing speed app. B.5. Algorithm 1 takes h, while Algorithm 2 uses 0. Choose one convention and apply it everywhere.
• The guarantee is per step and assumes the model is right. The paper's limitations say real systems may break the regularity assumptions, that a poor embedding can make the reduced subspace unsafe, and that an unexplored problem can cause more unsafe choices because the method is optimistic about safety §6, Limitations paragraph. The paper advises checking the embedding's distance preservation before online use.
• Cubic cost of the GP. The paper says the number of evaluations is limited by the cubic cost of the GP, which rules out very long horizons §6.
• Clinical numbers conflict. Section 5.3 says 504 additional trials followed 132 initial configurations. Appendix B.7 says 441 of 504 trials were recommended by the algorithm, and that the first 132 trials were the initial data. Whether the 504 includes the 132 is not clear.
• Neural-stimulation improvements are ambiguous. The eight values in §5.3 are not labelled as baseline or result. Do not read them as improvements without checking the source.
• Clinical data cannot be reproduced. The experiment involved one patient, with therapist-supervised trials and safety ratings from the patient and the therapists app. B.7. Only the setup is described. No raw data is released.
• Synthetic functions. Because effective variables can repeat in the code, the true effective dimension may be below 40 for some runs [code: task/gpfun_task.py#L43].
• Figure-only results. The musculoskeletal, ablation, β-sweep and neural-stimulation results are in figures only fig. 2 fig. 3 fig. 6 fig. 8. Their values are not in this card.

Lineage

• Builds on: local trust-region BO, cited as Eriksson et al. 2019 §4.2, and SCBO, cited as Eriksson et al. 2021 §2. The safe BO baseline SafeOpt is cited as Sui et al. 2015 §2. The GP-UCB confidence bound is cited as Srinivas et al. 2009 app. A.2. The IRVAE is cited as Yonghyeon et al. 2021 §5. None of these has an arXiv id in this card's sources.
• Compared against: LineBO, SCBO, CONFIG, cEI and CMA-ES §5.1; HesBO and BAxUS in the neural-stimulation task app. C.5; SafeOpt, which failed to run §5.1.
• Superseded: not stated in this card's sources.