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

Authors: yunyue  wei
zeji  yi
hongda  li
saraswati  soedarmadji
yanan  sui
Venue: robot-learning_org_CoRL_2024_Conference
Published: 2024-06-07 13:52:10

Abstract: Learning to move is a primary goal for animals and robots, where ensuring safety is often important when optimizing control policies on the embodied systems. For complex tasks such as the control of human or humanoid control, the high-dimensional parameter space adds complexity to the safe optimization effort. Current safe exploration algorithms exhibit inefficiency and may even become infeasible with large high-dimensional input spaces. Furthermore, existing high-dimensional constrained optimization methods neglect safety in the search process. In this paper, we propose High-dimensional Safe Bayesian Optimization with local optimistic exploration (HdSafeBO), a novel approach designed to handle high-dimensional sampling problems under probabilistic safety constraints. We introduce a local optimistic strategy to efficiently and safely optimize the objective function, providing a probabilistic safety guarantee and a cumulative safety violation bound. Through the use of isometric embedding, HdSafeBO addresses problems ranging from a few hundred to several thousand dimensions while maintaining safety guarantees. To our knowledge, HdSafeBO is the first algorithm capable of optimizing the control of high-dimensional musculoskeletal systems with high safety probability. We also demonstrate the real-world applicability of HdSafeBO through its use in the safe online optimization of neural stimulation-induced human motion control.

Main paper: 
 introduction: 
 Some robotics applications require online optimization of control policies for performance while avoiding unsafe parameter tuning that could potentially damage the systems. These scenarios correspond to the problem of safe exploration, which involves the sequential optimization of an unknown objective function under the constraint of satisfying unknown safety conditions. Bayesian optimization (BO) is an effective paradigm in optimizing black-box functions, and safe BO methods have been successfully used to tune control parameters for various robotic systems [1, 2, 3] . Most existing safe BO methods use the Gaussian process (GP) to model the underlying safety function, and discriminate safe regions with estimated function's lower confidence bound to ensure safety with high probability [4] . Such conservative strategies are inefficient for objective optimization, and even infeasible in high-dimensional input settings, such as human or humanoid system control.A motivating application of our work is the control of musculoskeletal (tendon-driven) systems, where complex motions are driven by dozens to hundreds of muscle-tendon units rather than joints. Such overactuated embodied systems introduce additional control challenges within a large-scale parameter space. Under this high-dimensional input setting, efficiently optimizing the task function while maintaining safety remains a formidable challenge for existing safe optimization algorithms. Despite considerable efforts to leverage BO for solving high-dimensional constrained optimization problems, these methods often fail to incorporate safety considerations into the optimization procedure [5, 6, 7] . To our best knowledge, there are currently no methods that guarantee safety, or probabilistic safety during the high-dimensional optimization.In this paper, we introduce High-dimensional Safe Bayesian Optimization with local optimistic exploration (HDSAFEBO) to address probabilistic safety in high-dimensional sequential optimization problems, as shown in Figure 1 . To achieve efficient and safe optimization, we propose a local optimistic strategy with probabilistic safety guarantee and cumulative safety violation bound. By using isometric embedding for dimension reduction, we enable HDSAFEBO to handle even higher dimensional problems, ranging from a few hundred to several thousand dimensions, while maintaining the probabilistic safety guarantee. Our experimental results show that HDSAFEBO efficiently learns to control a musculoskeletal system with high safety probability -a task where all baseline methods fail. We also demonstrate the success of HDSAFEBO in real-world experiments to safely optimize the control of neural stimulation induced human motion. Our project page is at https://lnsgroup. cc/research/hdsafebo Figure 1 : Workflow of HDSAFEBO. Local optimistic safe optimization is employed to efficiently optimize the objective function while guaranteeing probabilistic safety. We utilize isometric embedding to reduce the problem dimension to deal with high-dimensional inputs.2 Related Work 
safe bayesian optimization: 
 The sequential decision-making problem with safety constraints has been extensively studied, varied by the definition of safety. To achieve full safety during exploratory sampling, algorithms have been proposed with theoretical guarantee in near-optimality and safety with high probability [4, 8, 1, 9, 10] . These methods have been applied in safe parameter tuning of quadrotor [11] , robot arm [2] and quadruped robots [3] . These safe optimization algorithms conservatively estimate and expand the safe region, leading to inefficient optimization performance.In contrast to a zero-tolerance approach to unsafe actions, an alternative approach allows for limited constraint violations within a predefined budget, trading safety for more efficient optimization [12, 13, 14] . A recent work CONFIG uses upper confidence bound to optimistically estimate the safe region, enjoying global optimal guarantee as unconstrained methods [15] .Another extreme case, called constrained Bayesian optimization, has also been used for finding safe controller parameters [16, 17, 18, 19, 20] . However, this line of methods aim only to find the best feasible solution, neglecting the safety during the optimization process. Constrained expected improvement (cEI) is a popular constrained BO algorithm that introduces feasibility constraints to acquisition function formulation [21, 16] .All of the aforementioned methods fall in the framework of Bayesian optimization, which is typically limited to low-dimensional problems [22, 23] . LineBO demonstrates success in optimizing problems with dimension up to 40 via searching over 1-dimensional subspace at each iteration [24, 25] . When parameter size exceeds 10, safe optimization algorithms might even struggle to expand the safe region due to sparse discretization of the input space.A large body of literature leverages dimension reduction to apply BO over a low-dimensional subspace. Several works use variable selection to identify and optimize important dimensions during optimization [28, 29, 30] . Another popular approach for reducing the search space is random linear embedding, which has been proven to contain the optimal solution with a certain probability relative to the objective function's effective dimension [31, 32, 33, 34, 35] . Many works also use autoencoder to learn a non-linear mapping between the original space and the latent space [36, 37, 38, 39, 40, 41, 6] .Another line of works utilizes local search to address the over-exploration issue in high-dimensional optimization [42] , achieving better empirical performance than global BO methods [43, 44, 45, 46] . This local optimization strategy is also able to be combined with dimension reduction [47] and extended to constrained setting (SCBO [7] ). Local BO methods have also been deployed for optimizing locomotion control in simulation [48, 30] .Besides BO methods, evolutionary algorithms such as CMA-ES are competitive to solve highdimensional problems [49] and can be extended to constrained setting [50, 51] . Although many works attempt to solve high-dimensional constrained optimization problems, to the best of our knowledge, there lacks work that addresses safety in high-dimensional sequential optimization with at least hundreds of variables. 
problem formulation: 
 We aim to optimize an unknown objective function f : X → R by sequentially sampling points x 1 , . . . , x t ∈ X . We can also get observations of safety measurement from another unknown function g : X → R. We define a point x is safe when g(x) > 0. We can formally write our optimization problem as follows:EQUATIONWe consider the above problem formulation widely exists in the application of robotic control, where x represents a parameterized controller, and f (x), g(x) are the utility and safety measurement of the controller in a single experimental trial. In our target applications such as the control of embodied systems, the input space X ∈ R D is high-dimensional, encompassing dozens to thousands of variables.In Bayesian optimization, Gaussian process is usually used as the surrogate model to learn the unknown functions. Taking g as an example, for samples at points X t = [x 1 ...x t ] T , we have noise-perturbed observations y t = [y g 1 ...y g t ] T . The GP posterior over g is also Gaussian with mean µ t (x), covariance k t (x, x ′ ) and variance σ 2 t (x, x ′ ) under kernel function k:EQUATIONwhere k t (x) = [k(x 1 , x), ..., k(x t , x)] is the covariance between x and sampled points, K t is the covariance of sampled positions: [k(x, x ′ )] x,x ′ ∈Xt . Similarly we can use GP to derive posterior of f . Using the posterior of GP, we can define define u t (x) := µ t-1 (x) + β t σ t-1 (x) as the upper confidence bound (UCB) of the function estimation, where β t is a scalar which can be properly set to contain g(x) with desired probability. we make the following regularity assumptions that are commonly used in the field of Bayesian optimization: Assumption 3.1. g and f are samples of two Gaussian processes defined by the kernels k(•, •), k ′ (•, •) respectively. The observations are perturbed by i.i.d. Gaussian noise:y g (x t ) = g(x t ) + n t , y f (x t ) = f (x t ) + n ′ t where n t ∼ N (0, σ 2 ), n ′ t ∼ N (0, σ ′2 ) Algorithm 1 HDSAFEBOInput Sample set X , GP priors GP f , GP g , safety threshold h, acquisition function A, initial dataset D X 0 , trust region length l 0 , isometric encoder Π and decoder Π -1 1: c s , c f ← 0 2: x * , y * ← argmax {x,y f }∈D0,y g (x)>0 y f (x) 3: for t = 1 to . . . do 4:D Z t-1 ← Π(D X t-1 ) 5:Update GP f , GP g using D Zt-16:L t ← {z ∈ Z | Π(x * ) -l t ≤ z ≤ Π(x * ) + l t } 7: S t ← {z ′ ∈ L t | µ t-1 (z) + β t σ t-1 (z) ≥ 0} 8:z t ← argmax z∈St (A(z))9:x t ← Π -1 (z t )10: x * , y * ← x t , y f t ; c s ← c s +1; c f ← 0 3: else 4: In high-dimensional space, existing safe optimization methods are too conservative to efficiently optimize the objective function, and would be infeasible due to sparse discretization. Therefore we aim to improve sample efficiency by slightly relaxing the safety constraint to a probabilistic version. The probabilistic safety means each sample point is safe with a probability above predefined threshold α, that is Pr(g(x t ) ≥ 0) ≥ α, ∀t ≤ T . In this sense, we introduce HDSAFEBO, an innovative algorithm designed for ensuring probabilistic safety while optimizing in a high-dimensional space. The workflow of this algorithm is illustrate in Figure 1 and Algorithm 1.y f t ← f (x t ) +c s ← 0; c f ← c f + 1 5: end if 6: if c s = τ s then 7: l t ← min(2l t-1 , l max ); c s , c f ← 0 8: else if c f = τ f then 9: l t ← max(HDSAFEBO first uses the isometric encoder Π to reduce the problem dimension, converting the original dataset from X to the low-dimensional latent space Z ∈ R d (Line 4). Isometry means the embedded subspace is able to preserve the distance of the original space according to the corresponding metric dX , d Z , i.e. d X (x, x ′ ) = d Z (Π(x), Π(x ′ )).Leveraging the historical record of function observations, we compute the posterior and confidence interval of both the objective and safety functions through distinct Gaussian processes. (Line 5). Then we define a local region to search over (Line 6), and identify the safe space within the local region using GP upper confidence bound (Line 7). We optimize the acquisition function over the safe space, and project the recommended latent point back to original input space using the decoder Π -1 (line 8-9). Thompson sampling (TS) [52] was selected as the acquisition function A due to its compatibility with the discrete nature of our safety estimation and search space, and its innate ability for batch optimization by sampling the GP posterior-an appropriate choice for high-dimensional tasks that support parallel evaluations. Finally the history data is updated via evaluating new inputs. (Line 10-12), and the local region parameters are updated based on the sample results (Line 13).In conjunction with isometric embedding, we highlight two important components to improve the efficiency and safety: Optimistic Safety Identification and Local Search via Trust Region. For clear illustration, we first introduce algorithm details under identity mapping I(x) = x, which is a special case of isometric embedding. In this way, the embedded subspace is equivalent to the original space X , and inherits all the assumptions defined in Section 3. 
optimistic safety identification: 
 To mitigate inefficiency concerns, HDSAFEBO distinguishes the safe region through the Gaussian process upper confidence bound of the safety function. This optimistic strategy allows for optimization over distinct safe regions and is viable in high-dimensional problems. In many real-world applications, the search space is inherently constrained by domain priors, where a certain proportion of decisions are safe even under random search. Therefore, we contend that optimistic safety identification is a more practical strategy. Here we derive the appropriate choice of the scalar β t to ensure step-wise probabilistic safety. Proposition 4.1. Let Assumptions 3.1 holds for the latent safety function g, and set β t satisfying Φ(β t ) ≤ 1 -α. Then at every time step t:EQUATIONwhere Φ(•) is the cumulative distribution function (CDF) of the standard normal distribution N (0, 1).Our use of β t here is different from the common-used setting [53] , which would make u t (x t ) > g(x) with a high probability and violate the probabilistic requirement. Our results indicate that the upper confidence bound is permissible when the probability threshold is less than 0.5, representing the maximum acceptable level of safety violation. Otherwise, the algorithm reverts to a more conservative strategy, utilizing the lower confidence bound to identify safety. We also derive the upper bound of cumulative safety violation V T = T t=1 max(0, -g(x t )) of our optimistic strategy, where the proper choice of β t contributes to the reduction of safety violations: EQUATIONwhere γ T is the maximum information gain for g over X .In addition to optimistic safety identification, a trust region method is employed to dynamically pinpoint local search regions, which has demonstrated impressive empirical performance over highdimensional problems [43, 48, 7] . At each optimization round, a local search space is defined as a hypercube trust region around the current best safe point in the sample dataset. We design a safety sensitive strategy to update the trust region state, as illustrated in Algorithm 2. Specifically, a sampling round is considered "successful" if it finds a better reward while maintaining comprehensive safety (line 1-2). Conversely, it is labeled a "failure" if any unsafe points are found or if there is no discernible improvement (line 3-5). The side length is adjusted-increased for successes and decreased for failures-upon reaching a preset threshold (line 6-9). Unlike conventional local BO methods that discard all data and restart when the side length reaches its minimum, our approach resets l t to its initial length and retains all previous samples, ensuring a different, safer trajectory sampling than the initial instance (line 10-12).We also provide the theoretical implications of our local search strategy in reducing the safety violations. The maximum information gain γ T is positively correlated with the size of search space X . During the optimization, the actual search region of HDSAFEBO is restricted by the adaptive trust region, resulting in a lower maximum information gain compared to entire input space. Therefore the safety violation bound in Theorem 4.2 can be further reduced compared to global search.In addition to identity mapping, when applying safe optimization over the subspace of other embeddings, a pertinent question arises regarding whether safety guarantees persist in the original space.Here we demonstrate that, through the use of isometric embedding, the probabilistic safety guarantee in the embeded subspace can still be fulfilled in the original space when using stationary kernels.Proposition 4.3. Define d X (x, x ′ ) = |x -x ′ |, d Z (z, z ′ ) = |z -z ′ |, and u Z t (z) := µ t-1 (z) + β t σ t-1 (z)is the upper confidence bound estimated from GP over the embeded subspace Z. Suppose the GP kernel is stationary andd Z (z, z ′ ) = d X (Π -1 (z), Π -1 (z ′ )). At every time step t, if the point z satisfies that Pr(u Z t (z) ≥ 0) ≥ α, then we have Pr(g(Π -1 (z)) ≥ 0) ≥ αThe choice of the isometry embedding can be some linear mapppings, such as principal component analysis (PCA) when the number of principal components exceeds the effective dimension of the problem. Recent research efforts have also explored the utilization of deep neural networks to learn a subspace with both good representation and approximated distance-preserving property. We will show in experiment that HDSAFEBO is still able to maintain high safe probability when using approximated isometric embeddings. 
experiments: 
 In this section, we first utilize synthetic functions to evaluate the algorithm performance. Then we apply HDSAFEBO to safely optimize the control of high-dimensional musculoskeletal system. Finally we demonstrate the potential of HDSAFEBO in optimizing neural stimulation induced human motion control through both simulation and real experiments. Additional experiments are presented in Appendix C.Evaluation metrics. We assess the performance of the algorithm according to three metrics: the best feasible objective function value (Objective), the safe decision ratio of all samples (Safety), and the cumulative safety violation (Violation). The presented plots and tables display the means along with one standard error.Isometric embedding. We utilize PCA in synthetic function and musculoskeletal system control tasks. In the neural stimulation induced human motion control task, we explore the use of isometric regularized variational autoencoder (IRVAE [54] ) as the dimension reduction component.To evaluate the algorithm performance under full assumption satisfaction, we sample both objective and safety functions from Gaussian processes, with an effective dimension d e = 40, which is much lower than the input dimension D = 1000. We contrast the optimization performance with the following competitive algorithms target for high-dimensional constrained optimization: LineBO [24] , SCBO [7] , CONFIG [15] , cEI [21] and CMA-ES [49] . We attempted to run SafeOpt [4] , but it failed to expand the safe region from the initial points. We set PCA subspace dimension as d = 50.We show the optimization result in Table 1 . HDSAFEBO achieves significantly better optimization performance, higher safety decision ratio, and lower cumulative safety violation compared against all baselines. We establish an upper limb control task utilizing a musculoskeletal system [55] . The objective is to optimize the activities of 55 hand-related muscles to rotate and hold a bottle in vertical position (Figure 2 (a) ). As in Mania et al. [56] , we formulate the original reinforcement learning (RL) problem as a sampling problem, where the algorithms need to optimize a linear policy: x ∈ R |a|×|o| , where |a| = 55 and |o| = 65 are the dimensions of the action space (number of muscles) and the observation space, respectively. The policy to be optimized has D = 3575 parameters, making it a very highdimensional task. The objective function is set as the accumulated reward from the environment, and the safety function is defined as the landing speed of the bottle (see Appendix B.5). We collected muscle activation and use PCA to build the muscle synergies of performing the task, reducing action dimension from 55 to 5 (3575 to 325 for policy dimension). While the search space is significantly reduced, the remaining optimization problem is still high-dimensional.We find that all baselines fail to improve the objective when optimizing over the original parameter space. Therefore we conducted baseline runs over the subspace derived from PCA (denoted with (L)), where dimension reduction facilitated effective optimization. The optimization results are presented in Figure 2 (b). SCBO (L) and LineBO (L) are omitted from the figure as they fail to attain a positive reward. We observe even within the PCA subspace, all the shown baselines demonstrated lower efficiency and sampled more unsafe parameters compared to HDSAFEBO. Our proposed method stands as the pioneering algorithm to achieve efficient and safe optimization over high-dimensional musculoskeletal system control. We also conducted an ablation study on the components of HDSAFEBO in Appendix C.2, where both optimistic safe identification and local search were shown to contribute to a safer optimization process. Apart from direct activation of muscle-tendon units, locomotion of musculoskeletal systems can also be governed by muscle synergies arising from spinal cord stimulation, as observed in the central nervous system of vertebrates and certain neuromuscular robotic designs [57, 58] . While the neuromuscular system holds the potential to achieve robust control performance with higher energy efficiency compared to direct muscle actuation, the mapping from stimulation inputs to motion becomes less straightforward under complex neural systems. In this section, we showcase the success of HDSAFEBO in optimizing the control of intricate neural stimulation-actuated human motion through simulation and real experiments.As shown in Figure 3 (a), we want to improve the control of lower limb muscles via an electrode array implanted in the human spinal cord. By setting different parameters, the electrical stimulation delivered by the electrode array can induce patient's muscle activities, allowing us to control the lower limb movements of the patient. The stimulation parameter space consists of 16-contact spatial configuration (discrete) and current intensity (continuous). Our optimization goal is the selectivity index of target muscle groups (see Appendix B.6). A higher selectivity index indicates better control over the target muscle group and less influence over non-target muscles. To enhance the physical interpretability of discrete input representation, we transformed 17-dimensional vectors into an electric field map with 52 × 14 pixels using simplified computation. Subsequently, we generate unlabeled synthetic data with clinical priors and train an IRVAE with a 16-dimensional latent space to embed and reconstruct the electric field map, with distance-preserving checks in Appendix C.3.Simulation over neuromuscular model. In simulation, we use a human neuromuscular model as the function oracle, which is capable of computing the evoked electric field around the spinal cord given certain stimulation parameter, and inferring the lower limb muscle activation (see Appendix B.6). The maximum induced muscle activation is used as the safety function to avoid causing harm to the patient during optimization. We choose to optimize the selectivity of emitendinosus (ST) and gastrocnemius (GA) for their importance during human walking. The simulation results are depicted in Figure 3 (b) . Despite the subspace space dimensions being comparable to the original space, HDSAFEBO achieves safe exploration by optimizing on the continuous manifold, which obtains the best control performance while maintaining the highest safety selection ratio and lowest cumulative safety violation compared to other algorithms.Real-world experiments on paraplegic patient. We further applied HDSAFEBO to improve the motor function of a paraplegic patient with the same electrode array implanted. Starting with 132 initial configurations from clinical prior, a total of 504 additional trials were conducted with the patient over a 1-month period. We observed selectivity improvement of 7 out of 8 target muscles compared to the baseline (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). During the whole experimental procedure, only three configurations recommended by HDSAFEBO were rated as unsafe, which evoked large lower limb movements but no physical damage or pain. Our real-world experiment underscores the practicality of HDSAFEBO in safely optimize the control of complex neuromuscular systems. 
conclusion: 
 We develop HDSAFEBO for optimizing of the control over high-dimensional embodied systems under safety constraints. Our proposed method employs a local optimistic safe strategy to optimize the objective function and expand the safe region, with probabilistic safety guarantee and cumulative safety violation bound. HDSAFEBO is able to optimize high-dimensional input ranging from a few dozen to several thousand variables with safety guarantee. The algorithm can efficiently optimize the control of high-dimensional human musculoskeletal systems with high safety probability, and successfully optimize human motion control via neural stimulation in real clinical experiments. HD-SAFEBO has great potential to safely optimize the control of real-world high-dimensional embodied systems online.Limitations. While we provide the probabilistic safety guarantee for HDSAFEBO, real-world applications may fail to fully satisfy the theoretical assumptions regarding function regularity. The imperfection of the trained embedding could lead to unsafe behavior when optimizing over a reduced subspace. It is important to pre-check the distance-preserving quality of this subspace before conducting online optimization, and improvements can be made by synthesizing additional unsupervised data based on domain knowledge. The parameter space in our real human experiment has been restricted by domain prior. Directly applying HDSAFEBO to a completely unexplored problem may cause more unsafe decisions due to its optimism to safety. The total evaluation number of HDSAFEBO is constrained by cubic complexity of Gaussian process. The computational complexity prevents the use of HDSAFEBO in very long time horizon optimization problems.Proof.EQUATIONEQUATIONEQUATIONEQUATIONEQUATIONwhere the inequality (18) follows by Lemma A.2, the inequality (19) follows by the monotonicity of β c t , and the inequality (20) follows by Lemma A.3. 
 
 Appendix: Proof. Define k X , k Z are the kernel function of GP with same hyperparameter over X and Z respectively. For stationary kernel, the kernel function value between x and x ′ is only depend on |x -x ′ |, therefore we can write the k X , k Z as a function of metric d X , d Z :EQUATIONWe define k X as the kernel function of GP over X , and u X t as the estimated upper confidence bound of GP over X . Utilizing the property of isometric embedding, we can obtain equivalent GP estimation between original space and embeded subspace:EQUATIONTherefore the estimated upper confidence bound is also equivalent. For every sample z from HDSAFEBO, we haveEQUATIONEQUATIONthe statement holds.The full implementation of our experiments can be found on our project page: https://lnsgroup. cc/research/hdsafebo. Our musculoskeletal model will be released soon. In the meantime, the model can be accessed for research purposes upon request (ysui@tsinghua.edu.cn).We employ PCA using the scikit-learn library1 with default parameter settings.We employ IRVAE on neural stimulation task and digital generation task by directly using the paper's original repository 2 . We use MLP as the VAE module for all tasks, and list the model detail in Table 2 . We train all models for 300 epochs using Adam [60] optimizer with a learning rate of 0.0001. For implementation of HDSAFEBO, We use BoTorch as the GP inference component [61] . We also use BoTorch to replicate LineBO, SCBO, CONFIG, cEI and cEI-Prob. For LineBO, we choose to implement the random line embedding version shown in the main paper. We use the package pycma 3to run CMA-ES on benchmarks.All GP-based methods uses matérn kernel and fits kernel parameters after each iteration. During the experiment, we set τ s , τ f , l 0 , l max , l min for HDSAFEBO and SCBO (the default setting of SCBO). We use Thompson sampling for LineBO and CONFIG, and use inherent acquisition function for other BO baselines. Confidence scalar β is set as 2 for HDSAFEBO, LineBO and CONFIG across all experiments. We set the latent optimization bound as the upper bound and lower bound of training points in the latent space. Other baseline parameters are set to default values as in the original implementation.During the experiment, we set the sample size to 10 for all tasks in the main paper, and sample one point each iteration in constrained digital generation task.We sample objective and safety function from separate Gaussian process with Matérn kernel and length scale as 0.05, which is implemented using GPyTorch4 . We set the safety threshold to -0.75.For each independent run, we use a random linear projection Π rand ∈ R D×d to create d-dimensional latent space, and randomly select d e variables as the function effective dimension.We use a full-body musculoskeletal model which actuates locomotion by controlling muscle activation.Here we only control the right hand part (below elbow), and fix other joints, leading to 55 muscles and 28 joints. The primary task is to control hand muscles to maintain a steady vertical grip on a bottle. At the beginning of the episode, the bottle is initially positioned horizontally in the hand. The initial task involves first rotating the bottle to achieve and maintain a vertical orientation. At each time step, the reward from the environment is computed as follows: r = r pose + r bonus -10 * r penalty + r grasp + 2 * r survive -r activation -100 * r drop (26) where r pose is the difference between the bottle and vertical orientation, computed by Euler angle. r bonus is the reward when the difference falls below a predefined threshold. r penalty is positive when the bottle position is out of from the predefined range. r grasp is the distance between the centroid of the bottle and palm joints. r survive is the reward for not dropping the bottle during the current time step. r activation is the penalty for large muscle activations. r drop is the penalty when the bottle drop from hand. The overall simulation is based on Mujoco [62] .When the height of the bottle is below 0.4m, we consider it to be dropped from the hand and the episode concludes. We record the speed of episode ending as the landing speed of the bottle. We use a safety threshold of 4.4, which is the average landing speed when randomly sampling the environment. We train a Soft Actor-Critic agent [63] for 150k time steps under different observation setting, and rollout for 1000 episode to build a muscle activation dataset with 107, 439 datapoints.The performance video of HDSAFEBO and other baselines (in the subspace from PCA) is shown in the supplementary file. HDSAFEBO is capable of quickly rotating the bottle and hold vertically and steadily, while other algorithms either hold the bottle non-vertically, or learn drop the bottle safely to avoid penalty of wrong orientation.We developed an average model of the human spinal cord based on anatomical statistics(model paper under review, [64, 65, 66, 67, 68] ). The model contains gray matter, white matter, nerve roots, cerebrospinal fluid (CSF), and dura mater of T12-S2 segments of the spinal cord which are related with the motor control of lower limbs. The specific conductivity values of the modeled tissues were set refer to [69] . Electric fields induced by different stimulation parameters were derived using finite element method (FEM). To calculate the stimulation effects for different muscles, we redistricted the cord model according to reported results of the segmental innervation for lower limb muscles( [70, 71] ). Six groups of muscles of bilateral lower limbs were studied: iliopsoas (IL), vastus lateralis and rectus femoris (VL&RF), tibialis anterior (TA), biceps femoris muscle and gluteus maximus (BF&GM), semitendinosus (ST), and gastrocnemius (GA).To evaluate the selectivity of stimulation for certain muscle, we used a selectivity index (SI) to characterize the distribution of the electric field. The selectivity index for the ith muscle was defined as follows:EQUATIONwhere m neighbor represents the number of muscles whose motor neuron pools are adjacent to the ith muscle's. The selectivity index ranges from -1 to 1, where -1 represents the maximum of activation of all undesired muscles with a complete absence of activation of the targeted muscle, 0 indicates that all muscles are activated at the same level, and 1 means the targeted muscle is activated at the greatest extent while no undesired muscles are activated. And µ i is the normalized activation of the ith muscle and is defined as follows in the simulation:EQUATIONEQUATION)Ω i is the segmental volume of the ith muscle in the cord. AF is the activating function, defined as the second spatial derivative of extracellular voltage along an axon( [72, 73] ).We use the spinal model to traverse all stimulation parameters with 1 cathode with anodes no more than 3, and 2 cathodes with anodes no more than 2, leading to a spinal cord stimulation (SCS) dataset with 218,000 stimulation parameters and predicted muscle activation. We compute the objective function using 27, and compute the safety function as g(x) = 1 -max i (µ i ). The selectivity index distribution is shown in Figure 4 . We set the safe threshold as 0.05, with nearly half of the traversed parameters are safe.We convert electrode parameters from 17d vector to a 2d electric field image using simplified computation. In concrete, we map contact combinations to the spatial position in the electrode, linearly compute the diffusion of electrical field from each cathode and anode, and multiply the map by current intensity. One example of generated electric field image is shown in Figure 5 .We use the constructed SCS dataset as the function oracle. For a given 2d map, we set its function value as the function value of the nearest unevaluated point in the dataset measured by the electrical field map. We employ HDSAFEBO in the treatment of spinal cord stimulation to find more selective stimulation parameters for different muscles. The clinical experiments received approval from the Institutional Review Board (IRB) of the hospital. All the trials were conducted under the supervision of therapists.The patient was seated in the wheelchair in a comfortable way and was told to relax. During the first period, typical parameters which were usually used in the therapy (e.g. bipolar stimulation) were delivered to the patient while the evoked muscle activities were recorded using EMG. These data (132 trials) were used to initialize HDSAFEBO. Except for the first 132 trials as the initial data, 441 out of 504 trials are recommended by HDSAFEBO. The other trials were conducted by the therapist.We focus on 8 groups of muscles: iliopsoas (IL), rectus femoris (RF), tibialis anterior (TA), and biceps femoris (BF) for both sides. The clinical selectivity index was defined as following: (30) where µ i represents the normalized peak-to-peak value of the evoked EMG for the ith muscle and m is the total number of the target muscles.SI i = µ i 1 + m j̸ =i µ jDuring the optimization, We restrict the number of cathodes and anodes to only evaluate configurations consistent with clinical priors. We used the Electromyography (EMG) to compute the selectivity index and queried safety scores from the patient and the therapists. For each trial, our algorithm recommended the parameter based on the history data and configured it onto the stimulator, which would deliver electrical stimulation to the patient for 800 ms at a frequency of 10 Hz. Peak-to-peak values were averaged and normalized to obtain selectivity indices of different muscles after stimulation. The calculated feedback and queried safety score were used to update the optimizer and it would recommend a new parameter. The tasks were refined sequentially and all the history data were reused when optimizing a new task.We run baselines on the same embedded space as HdSafeBO in optimizing GP synthetic functions (Table 3 ). We observe that, overall, the optimization and safety performances slightly improved. However, HdSafeBO still outperforms these baselines, as the reduced space remains too highdimensional for them.Table 3 : Algorithm performance on GP synthetic functions. We show the averaged performance of 500 evaluations over 100 independent runs. "L" indicates baselines optimize over the embeded latent space, "O" indicates baselines optimize over the original input space. In the musculoskeletal system control task, we conducted an ablation study of two components in HdSafeBO: local search and optimistic safe identification (Figure 6 ). We observed that without local search, the algorithm tends to over-explore, leading to degraded optimization and safety performance in this high-dimensional problem. Without optimistic safe identification, the algorithm makes more unsafe selections during the early stages of optimization. Combining these two components enables a safe and efficient optimization procedure. 8th Conference on Robot Learning (CoRL 2024), Munich, Germany.https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.htmlhttps://github.com/Gabe-YHLee/IRVAE-publichttps://github.com/CMA-ES/pycmahttps://gpytorch.ai/ 


Decision: In general, the reviewers are positive about this paper. Most reviewers like the simplicity and effectiveness of the method.
There are still some points to clarify. Specifically, the authors should try to highlight the novelty of the approach and frame it better in the literature. On top of that, I would strongly advise the authors to add more BO baselines to the comparison.
Another minor issue is the proper handling of the limitation section, which should make clear when the safety properties hold and how the safety problem is relaxed, particularly when the encoder is not perfectly trained.

=====

After the rebuttal, most of the reviewers were convinced by the author's response. The additional ablations to assess the impact of algorithmic design choices were particularly appreciated.
Please, make sure to open-source your code to enhance the reproducibility of this work, as discussed with the reviewers.

Reviews: review: title: title: A good submission (although I cannot judge novelty)
paper_summary: paper_summary: The paper presents a Bayesian optimisation algorithm that targets optimisation in high-dimensional spaces while providing probabilistic safety guarantees. The authors use a combination of isometric dimensionality reduction, trust-region approaches, and confidence bounds to realise the algorithm. The proposed combination of these techniques consistently outperforms alternative methods in experiments.
main_review: main_review: Strengths:

The paper is well written and easy to follow (with a working knowledge of Bayesian optimization).
The algorithmic ingredients are simple yet seem to work effectively in combination.
Although no robot experiments are presented, the authors performed a real-world study on a paraplegic patient.
For all experiments, the authors provide either explanations of the objective functions in the appendix and/or provide videos of the optimised behaviour.
I appreciate the connection between the trust-region size and the maximum information gain which in turn influences the cumulative safety violation. Although no formal proof, this argument helps understanding the safer exploration of HDSafeBO.

Weaknesses:

I do not see major weaknesses in general with the paper. The method seems to be sound and working. 

Apart from experimental improvements which I'll mention next, I simply want to state that I am not an expert in the related literature of this domain and hence cannot judge whether similar approaches have already been developed. I am only mentioning this because the individual algorithmic ingredients are relatively simple which is why I would be surprised if nobody else had tried this before. If indeed the authors are the first to try this out, I congratulate them for creating this simple yet effective recipe for BO.

In terms of experiments, I would have the following suggestions for improvement:
 * It would be interesting to see the baseline method performances in Section 5.1 and 5.3 if they were granted the same embedding as HDSafeBO. In Section 5.2, the authors already granted access to the embedding since no baseline method could otherwise deal with the search space otherwise.
 * Given the theoretical impact of the trust-region on the safety of the algorithm, I would be highly interested in an ablation of the method without the trust-region.
 * Furthermore, I would appreciate an experiment that investigates the algorithm's performance w.r.t. to the quality of the isometric embedding. While the authors note that they did not observe problems with imperfect embeddings, it would be interesting to see just how imperfect they can be.

Comments:

 * Is the reward function for the musculoskeletal model control in the appendix wrong? It seems like the signs are flipped since dropping the bottle leads to a large positive reward.
 * In Algorithm 2, shouldn't it be $\min(2 l_{t-1}, l_{max})$ in line 7 and $max(0.5 l_{t-1}, l_{min})$ in line 9?
questions: questions: Sea weaknesses and comments section above
limitations: limitations: 3
review_summary: review_summary: I think the paper should be accepted at the conference if there is no significant algorithmic overlap with prior work.
title: title: Review
paper_summary: paper_summary: This paper proposes a method for safe high-dimensional control through Bayesian optimization (BO). The typical scalability problems of BO are mitigated through various dimensionality reduction techniques. Optimistic safety identification and local search are used to improve performance. The method is evaluated on two high-dimensional control tasks, outperforming baselines.
main_review: main_review: Strengths:
- The paper tackles the important problem of high-dimensional safety-critical control.

Concerns:
- It’s not 100% clear what the paper’s technical contribution is over existing methods. For instance, several papers in the related work (like [36-41]) also use dimensionality reduction techniques, e.g., autoencoders, and then run BO in the latent space. To my understanding, the primary novelty is in using optimism in sample selection (to provide some safety), and in integrating local search. Thus, testing two ablations where HDSafeBO is run without the optimistic sample selection/the local search procedure would help make the contribution more clear.
- In all practical situations, autoencoder training is not exact, e.g., encoding and decoding a high-dimensional sample does not exactly return the original sample. As such, \Pi^{-1} is not exact, and thus to my understanding Proposition 4.3 does not hold for this practical setting. It should be made clear in the text, or in the limitations section, that the safety guarantees do not hold in practical situations where the decoder is not trained to be a perfect inverse of the encoder.
- How well upper-confidence sample selection works seems entirely dependent on the size/shape of the feasible set of the specific problem. One can imagine that for highly-constrained problems, upper confidence would lead to frequent safety violations. Some discussion in the paper with regard to when upper-confidence sample selection can be expected to work would be helpful.
- Line 156-158: the chosen 0.5 threshold is not motivated. Why is 0.5 the maximum acceptable level of violation?
- Proposition 4.1 reasons about how each individual point is safe with some probability, and then a guarantee on the cumulative safety violation is constructed using the individual safety violation probabilities per trial using I believe Boole’s inequality. However, wouldn’t one want to directly compute the joint probability of the T trials violating the safety constraints, as this would be less conservative? This is done, for instance, in [A].
- Line 134-135: “the discrete nature of our safety estimation and search space” -> this is unclear to me - aren’t the experiments run for continuous control problems? For instance, Sec. 5.2 is optimizing over continuous parameters in a linear control policy (line 226).
- Why is Assumption 3.1 not also required for the function f?
- \Pi is never defined formally in the text (only mentioned in the algorithm block)

[A] Chou, Wang, Berenson. Gaussian Process Constraint Learning for Scalable Chance-Constrained Motion Planning from Demonstrations. RA-L 2022.
questions: questions: Please see review
limitations: limitations: 2
review_summary: review_summary: I currently recommend a weak reject due to a number of unclear points and claims in the writing, and an unclear technical contribution. I'm open to raising my score if the concerns raised in the review are addressed.
title: title: Safe Bayes Opt Corl 2024
paper_summary: paper_summary: This paper presents a novel approach (HDSAFEBO) to safely optimize the control of high-dimensional systems at the example of musculoskeletal systems. The results indicate that the HDSAFEBO outperforms existing approaches in terms of efficiency while still having a high safety probability. Two examples are showcased: the control of an upper limb (55 muscles) rotating a bottle in simulation and a neural stimulation of a lower limb (both in simulation and with a paraplegic patient).
main_review: main_review: Strengths:
- The authors introduce a novel method that allows to optimize high-dimensional control inputs to robotics systems with high safety probability, specifically targeting musculoskeletal systems.
- The inclusion of both simulated and real-world examples, such as the control of an upper limb with 55 muscles and neural stimulation of a lower limb in a paraplegic patient, enhances the practical relevance and impact of the proposed method.


Weaknesses:


- Only one derivative-free algorithm that is not a BO method (CMA-ES) has been used for the comparison. 
- Clarification Needed: In Figure 2b, there is an inconsistency regarding whether the comparison with HdSafeBO was conducted over the subspace (L) derived from PCA. This needs clarification.
- Safety Constraints Definition: The specific safety constraints should be more clearly defined and justified. 
- The focus on Bayesian optimization methods could be better contextualized within the broader scope of optimization techniques.
Please also see questions for rebuttal below
questions: questions: -	In Fig. 2b, the performance of the optimization is shown. It is unclear whether the comparison with HdSafeBO was also conducted over the subspace (L) derived from PCA. The text (line 232) mentions this, but the figure does not have the subspace denotation L? Please clarify this inconsistency. 
-	Additionally, if HdSafeBO was not run over the subspace derived from PCA, it would be valuable to see this comparison to understand he influence of dimension reduction.  
-	What are the specific safety constraints in your study? Landing speed of the bottle/ speed of bottle at the end of the episode? In the 2nd case: Is the maximum induced muscle activation the constraint? If so, isn't this an input parameter to the model? Would it not be more meaningful to characterize maximum joint angle velocities or joint violations as measures for pain and wear/tear of the system?
-	Regarding the fairness of comparison, can you provide an estimation of the run-time comparison between the algorithms? In Fig. 2 and 3 you plot the (commonly used) performance over the evaluation number. However, in real robotics applications, the actual time to optimize is crucial. If I can run twice as many optimizations with LineBO compared to HdSafeBO in the same time, the difference observed in Figure 3 might be less pronounced. Could you address this aspect?
-	In Chapter 2.2, the title claims "high-dimensional BO," but the last paragraph also discusses CMA-ES as another method. This seems inconsistent. Additionally, the focus of the comparison is heavily on Bayesian optimization methods, while the introduction generally states that BO has advantages for optimizing black-box functions, among which many other methods exist. Could you clarify the reasoning behind this focus?
-	How does your method differ / outperform GOSAFEOPT that has been also used to safely tune legged locomotion controllers for high-dimensional dynamical systems using safe BO?
-	Concerning code and model availability, Reference 55 refers to a model that is not openly available. Please ensure that both the model and the code are provided as open-source upon acceptance. 

**Small comments**

-	Problem formulation: first you use x_n as index (l. 96), later x_t (eq. 1, l. 105) to denote the sampling points -> is this on purpose for spatial and temporal resolution or should it always be x_t? Please check carefully. 
-	L. 108: covariance between z and sampled points -> what is z? I assume you refer to the latent space but this is not clear here, please clarify. 
-	Algorithm 1: line 13, please check type setting again, arrow is weird
-	L. 227: what is the observation space? (num of muscles + num of joints?) 
-	L. 266 should be semitendinosus (s missing)
limitations: limitations: 2
review_summary: review_summary: This paper is technically sound, and the question of safe optimization in high-dimensions addresses a core robotics challenge. However, the comparison focuses almost solely on BO methods and neglecs other high-dimensional optimization methods.