Safe Bayesian Optimization for the Control of High-Dimensional Embodied Systems
Authors: Yunyue Wei
Zeji Yi
Hongda Li
Saraswati Soedarmadji
Yanan Sui
Venue: CoRL 2024 (2024)
Keywords: Safe Bayesian Optimization
High-dimensional Embodied System

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.

Safe Bayesian Optimization for the Control of High-Dimensional Embodied Systems
Yunyue Wei, Zeji Yi, Hongda Li, Saraswati Soedarmadji, Yanan Sui
Tsinghua University
weiy veiyy200mails. ee00mmmlsstsiinnnaaa edu.
zejiy0andrew. cmu. edu 1hd21 heimmaisssssinnuuue edu.
chenxuying240mails tainggghsssssssiiuaaaaaaa tsinghua. edu, cn ysui@tsinghua.co edu.
Abstract: Learning to move a primary goal for animals and robots, where en- suring safety is often important when optimizing control policies the embodied systems. For complex tasks such as the control of human 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 the search process. In this paper, we propose High-dimensional Safe Bayesian Optimization with local optimistic exploration (HDSAFEBO), novel approach designed to handle nigh-dimensional sampling problems under probabilistic safety constraints. We introduce a local optimistic strategy to efficiently and safely optimize the objec- tive function, providing 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 main- taining 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.
Keywords: Safe Bayesian Optimization, High-dimensional Embodied System
1
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) an effective paradigm in optimizing black-box functions, and safe BO methods have been successfully used tune control parameters for various robotic systems [1, 2, 3]. Most existing safe BO methods use the Gaussian process (GP) model the underlying safety function, and discriminate safe regions with estimated function's 's lower confidence bound ensure safety with high probability [4]. Such conservative strategies are inefficient for objective optimization, and even infeasible in high-ding digh-dimensional input settings, such as human or humanoid system control.
motivating application of our work the control of musculoskeletal (tendon-driven) systems, where complex motions are driven by dozens 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
8th Conference Robot Learning (CoRL. 2024), Munich, Germany.

---

maintaining safety remains 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 handle even higher dimensional problems, ranging from few hundred several thousand dimensions, while maintaining the probabilistic safety guarantee. Our experimental results show that HDSAFEBO efficiently learns to control 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::/t// //1ngroon. cc/ h/researrchhhdsafeeeeee.
GP Regression
Local Region Identification
Optimistic Safe identification
Safety
UCB
Isometric Encoder
W
Z;
Isometric Decoder
X
D
Objective
Acq Fun
DJ=
Figure 1: Workflow of HDSAFEBO. Local optimistic safe optimization is employed efficiently optimize the objective function while guaranteeing probabilistic safety. utilize isometric embedding to reduce problem dimension to deal with high-dimensional inputs.
Related Work
2.1 Safe Bayesian Optimization
The sequential decision-maie 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, 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 zero-tolerance approach unsafe actions, alternative approach allows for limited constraint violations within predefined budget, trading safety for more efficient optimization[12, [12, 13, 14]. A recent work CONFIG uses upper confidence bound to optimistically estimate the safe region, enjoying global optimal guarantee 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 methods aim only find the best feasible solution, neglecting the safety during the optimization process. Constrained expected improvement (cEl) is popular constrained BO algorithm that introduces feasibility constraints to acquisition function formulation [21, 16].
All the aforementioned methods fall the framework of Bayesian optimization, which typically limited 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 struggle to expand the safe region due to sparse discretization of the input space.

---

2.2 High-dimension Bayesian Optimization
Recent research efforts have delved into the utilization of Bayesian optimization in high-dimensional problems [26, 27]. In the main paper, our focus lies on dimension reduction- based and local BO methods. Further discussion about other relevant works can be found in Appendix D.
A large body of literature leverages dimension reduction to apply BO over low-dimensional subspace. Several works use variable selection to identify and optimize important dimensions during optimization [28, 29, 30J. 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 's effective dimension[31, 32, 33, 34, 35]. Many works also use autoencoder to learn 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 CMA-ES competitive to solve high- dimensional problems [49] and can be extended to constrained setting[50, 51]. Although many works attempt to solve high- -igh--imensional constrained optimization problems, to the best of our knowledge, there lacks work that addresses safety high-dimensional sequential optimization with at least hundreds variables.
3 Problem Formulation
We aim to optimize an unknown objective function f X + R by sequentially sampling points $1,.... 3y E X. We can also get observations of safety measurement from another unknown function
X' - IR. We define a point æ safe when g(x) > 0. We can formally write our optimization
problem as follows:
max f(x) subject to g(x,) N 0, 1 m
(1)
We consider the above problem formulation widely exists in the application of robotic control, where a represents a parameterized controller, and f(x), g(x) are the utility and safety measurement of the controller single experimental trial. In our target applications such as the control of embodied systems, the input space X' e F. high-dimensional. encompassing dozens thousands of variables.
In Bayesian optimization, Gaussian process is usually used as the surrogate model to learn the unknown functions. Taking as an example, for samples at points X, = [x,"....,""" we have noise-perturbed observations y, = uu!"". The GP posterior over is also Gaussian with mean (x)), covariance (xx x²) and variance (x x²) under kernel function k:
(=(x) = (xx) (Ki
k (z. x') = k<xx,, x') - (x) (K, + 02)) -
(xx)) (xx,
where (x)) = [R<(0; (n)aa.))::)) ...... a(),,,))) m) is the covariance between æ and sampled points, K; the covariance of sampled positions: [xx x x')) Similarly we can use GP derive posterior of f. Using the posterior of GP, confidence bound (UCB) of the function estimation, where 3r is 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:
define define (x)) #4-1 (x) + 8, (æ) as the upper
Assumption 3.1. I and f are samples of two Gaussian processes defined by the kernels k(-,.). k(.,)) respectively. The observations y (æ,) f(x)) ni, where ~ (1 n ~ (0,
perturbed by i..d. Gaussian noise: 3d (xx)) g(x.) ni.
3

---

Algorithm HDSAFEBO
Algorithm 2 SAFEUPDATE
Input Sample set X, GP priors GPI, GPs, safety threshold h, acquisition function A, initial dataset Do trust region length lo, isometric encoder II and decoder II
Input Current
sample Vt; vf current best sample x", ,y", current counters ca, cur- rent trust region length 4, success tol- erance Tex failure tolerance T,, initial trust region length l0, trust region limits mmax: mmin
1: Cs. Cy 4 0
2: x", ,, f- argmax (m.y)) (e=p/))DDorP())000" (x)
1: if yi (ax) 0 and yf > y then
3: for 1= == to do
4: D f 11(
2: x'. + 2t, ye cs t c + +ico t-
5: Update GPI, GP using D
3: else
6:
Le {= E 2 [1(x())) 4 S 2 = 11((=')))) ) +4}
4: Ca + 0g t cf +1
5: end if
7:
S. - Iz' E Lt Mt-1 (2) + 1,0)-1(2))
6: if C, = Ts then
8:
21 - argniax (A(z))
7: 1, x (nin(222-11 lmxx) ); Cs. Cy + 0
9:
Kt + (z))
8: else if ey = Ty then
10:
V t f (x1) +h
9:
l t- max 2411: mmin ); Cg. Cf +-
yi + g (1)) + n4
10:
if mmin then
11:
le f- lo
12:
Dj + D,-1 U+ Vi v
12: end if
13:
x¹, N¹ CC. 4
13: else
SAFEUPDATE(22: vf x², ,
Ca, 44,
14:
l. f- 4-11
14: end for
15: end if
4 High-dimensional Safe Bayesian Optimization
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 probability above predefined threshold a, that Pr(g(x:) (x;) 2 0) 2 a, Vt < T. In this sense, we introduce HDSAFEBO, an innovative algorithm designed for ensuring probabilistic safety while optimizing in high-dimensional space. The workflow of this algorithm is illustrate in Figure 1 and Algorithm 1.
HDSAFEBO first uses the isometric encoder II to reduce the problem dimension, converting the original dataset from X' to the dinneimenioonal latent space Z E Rd (Line 4). Isometry means the embedded subspace is able to preserve the distance of the original space according to the corresponding metric dx, dz, i.e. (x x²) = d2 (II(z), II(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 define 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 II-1 (line 8-9). Thompson sampling (TS) [52] was selected as the acquisition function due to 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 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) = æ, which 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.

---

4.1 Optimistic safety identification
To mitigate inefficiency concerns, HDSAFEBO distinguishes the safe region through the Gaussian process upper confidence bound the safety function. This optimistic strategy allows for optimization over distinct safe regions and viable in high-dimensional problems. In many real-world applications, the search space is inherently constrained by domain priors, where certain proportion of decisions safe even under random search. Therefore, contend that optimistic safety identification is more practical strategy. Here we derive the appropriate choice of the scalar 3r ensure step-wise probabilistic safety.
Proposition 4.1. Let Assumptions 3.1 holds for the latent safety function g, and set Be satisfying (3,) S 1 - a. Then at every time step t:
Pr(g(x) 2 p(- + 310(-1(x))) 2 a. Vx E X.
(3)
where (--) is the cumulative distribution function (CDF) of the standard normal distribution (0, 1).
Our use Be here different from the common-used setting [53], which would make uy > g(x)) with 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 more conservative strategy. utilizing the lower confidence bound to identify safety. We also derive the upper bound of cumulative safety violation Vi- = ='== max (0, -(x))) of our optimistic strategy, where the proper choice of 3e contributes to the reduction of safety violations:
Theorem 4.2. Let assumptions 3.1 holds for safety function g. Define 31 == 2log X and C = 8/ log(1 + o -)). With probability at least 1 5, the sample points of Algorithm at time step satisfy
B[V+] (1-
where YT is the maximum information gain for g over X.
4.2 Local search via trust region
In addition to optimistic safety identification, trust region method is employed to dynamically pinpoint local search regions, which has demonstrated impressive empirical performance over high- dimensional 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 sample dataset. We design safety sensitive strategy to update the trust region state, illustrated in Algorithm 2. Specifically, sampling round is considered "successful" if it finds a better reward while maintaining comprehensive safety (line 1-2). Conversely, it is labeled "failure" if any unsafe points found if there is no discernible improvement (line 3-5). The side length is decreased for dillrres-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 4 to its initial length and retains all previous samples, ensuring a different, safer trajectory sampling than the initial instance (line 10-12).
for successes and
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.
4.3 Safe optimization with isometric embedding
In addition identity mapping, when applying safe optimization over the subspace of other embed- dings, pertinent question arises regarding whether safety guarantees persist in the original space.
5

---

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 dx (a, x') = |æ - æ'), dz( =²) = 1z - z'1, and up (z) = 14-1(2)
- is the upper confidence bound estimated from GP over the embeded subspace Z. Suppose the GP kernel is stationary and (z( dII (2), -I-II-1-60 (2')) At every time step t, if the point
satisfies that Pr(ut (z) 2 0) 2 then we have (r(g(II (z)) 2 0) 2 a
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 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 apply HDSAFEBO to safely optimize the control of high-dimenssonal 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.
5.1 Synthetic Function Optimization
To evaluate the algorithm performance under full assumption satisfaction, we sample both objective and safety functions from Gaussian processes, with an effective dimension d. 40, which is much lower than the input dimension = 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 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.
Table 1: Algorithm performance on GP synthetic functions. The total sample budget 500 including 200 initial points. We show the averaged performance over 100 independent runs.
Metric
HDSAFEBO
LineBO
SCBO
CONFIG
cEI
CMAES
Objective (t)
3.96 + 0.15
3.07 + 0.04
2.97 1.076
2.91 : 0.05
2.9 0.03
2.7 + 0.04 0.78 - 0.01 38.65 = 1.89
Safety (t)
0.81 + 0.02
0.78 1 0.0
0.77 0.0
0.77 1 0.0
0.77 0.0
Violation (L)
27.42 + 4.01
36.59 4: 0.82
39.05 + 0.92
38.96 + 0.96
39.61 0.58
5.2
Optimization for the Control Musculoskeletal System
We establish upper limb control task utilizing musculoskeletal system [55]. The objective is to optimize the activities of 55 hand-related muscles to rotate and hold a bottle in vertical position

---

HdSateB0
CONFIG (L)
dEl
UMAES
1030
1.00
50
Dbjective
830
0.96
Violstin
630
0.90
400
0.85
20
0.80
0.75
200
400
800
800
200
400
600
800
Evalustion number
Evaluction number
Evaluction number
(a)
(b)
Figure 2: Optimization for the control of a museuloskeletal system. (a) Task illustration. (b) Optimization
performance averaged 50 independent runs. Arrows indicate the better direction.
(Figure 2 (a)). As Mania al. [56], we formulate the original reinforcement learning (RL) problem sampling problem, where the algorithms need to optimize linear policy: R |a|x ol 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 very high- dimensional 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 musele activation and use PCA to build the muscle synergies of performing the task, reducing action dimension from 55 5 (3575 325 for policy dimension). While the search space significantly reduced, the remaining optimization problem is still -igennimeniioal..
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 Figure (b). SCBO (L.) and LineBO (L.) are omitted from the figure as they fail 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.
#55s150
LiviO
SCBO
CONPIG
CMANES
Osedive
0.28
1.0
Moaton
0.24
08-
08:
0
200
400
600
500
200
400
800
800
200
400
600
000
Objesive
0.22
10-
0.20
08
06
Stimaation Neural Shimutttin induced Human Mtttion
200
400
600
800
200
400
600
803
200 400
600 100
Parameter
Muscla Activation
Evalusson number
Eviluation number
iuollatoonummer
(a)
Figure 3: Optimization for the control of neural stimulation induced human motion. (a) Task illustration. II, RF, TA, BF, ST, GA different group of target musele the lower limb. (b) Optimization performance the control of neuromuseular model for semitendinosus (ST), and gastrocnemius (GA), averaged over 10 independent runs. Arrows indicate the better direction.
5.3 Optimization for the Control of Neural Stimulation Induced Human Motion
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 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.
7

---

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 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 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 x 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 preseervngg checks in Appendix C.3.
Simulation over neuromuscular model. In simulation, we 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 gastroenemius (GA) for their importance during human walking. The simulation results are depicted Figure (b). Despite the subspace space dimensions being comparable to the original space, HDSAFEBO achieves safe exploration by optimizing 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, total of 504 additional trials were conducted with the patient over 1-month period. We observed selectivity improvement of 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 worl-world high-dimensional embodied systems online.
Limitations. While we provide the probabilistic safety guarantee for HDSAFEBO, real-world ap- plications 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. important to pre-check the eestance-preserving quality of this subspace before conduct- ing 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 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.
completely unexplored problem may cause
8

---

Acknowledgments
This project is funded by STI 2030-Major Projects 20222D0209400 and Tsinghua DuShi fund.
References
[1] Berkenkamp. Schoellig, and A. Krause. Safe controller optimization for quadrotors with gaussian processes. In 2016 IEEE international conference on robolics and automation (ICRA), pages 491-496. IEEE, 2016.
[2] B. Sukhija, M. Turchetta, D. Lindner, A. Krause, S. Trimpe, and D. Baumann. Scalable safe exploration for global optimization of dynamical systems. arXiv preprint arXiv:2201.09563 2022.
[3] D. Widmer, D. Kang, B. Sukhija, J. Hübotter, A. Krause, and S. Coros. Tuning legged locomotion controllers via safe bayesian optimization. Conference on Robot Learning, pages 2444 2464. PMLR, 2023.
[4] Y. Sui, A. Gotovos, J. Burdick, and A. Krause. Safe exploration for optimization with gaussian processes. In International conference on machine learning, pages 997-1005, PMLR, 2015.
[5] R.-R. Griffiths and J. M. Hernández-Looato. Constrained bayesian optimization for automatic chemical design using variational autoencoders. Chemical science, 11(2):577-586, 2020.
[6] Notin, J. M. Hermández-Lobato, and Y. Gal. Improving black-box optimization in latent space using decoder uncertainty. Advances in Neural Information Processing Systems, 34: 802-814, 2021.
[7] D. Eriksson and M. Poloczek. Scalable constrained bayesian optimization, In International Conference Artificial Intelligence and Statistics, pages 730-738. PMLR, 2021.
[8] Y. Sui, V. Zhuang, Burdick, and Y. Yue. Stagewise safe bayesian optimization with gaussian processes. In International conference on machine learning, pages 4781-4789 PMIR, 2018.
[9] M. Turchetta, F. Berkenkamp, and A. Krause. Safe exploration for interactive machine learning. Advances in Neural Information Processing Systems, 32, 2019.
[10] D. Baumann, A. Marco, M. Turchetta, and S. Trimpe. Gosafe: Globally optimal safe robot learning. In 2021 IEEE International Conference on Robotics and Automation (ICRA), pages 4452-4458. IEEE, 2021.
[11] F. Berkenkamp, A. Krause, and A. P. Schoellig. Bayesian optimization with safety constraints: safe and automatic parameter tuning in robotics. Machine Learning, pages 1-35, 2021.
[12] X. Zhou and B. Ji. On kernelized multi- armee bandits with constraints. Advances in Neural Information Processing Systems, 35:14-26 2022.
[13] C. Lu and J. A. Paulson. No-regret bayesian optimization with unknown equality and inequality constraints using exact penalty functions. Pane- 55(7):895-902, 2022.
[14] H. Guo, Z. Qi, and X. Liu. Rectified pessinistic-opiiiistic learning for stochastic continuum- armed bandit with constraints. In Learning for Dynamics and Control Conference, pages 1344. PMLR, 2023.
[15] W. Xu, Y. Jiang, B. Svetozarevic, and C. Jones. Constrained efficient global optimization of expensive black-box functions. In International Conference on Machine Learning, pages 38485-38498. PMLR, 2023.
[16] M. A. Gelbart, J. Snoek, and R. P. Adams. Bayesian optimization with unknown constraints. 2014.
arXiv preprint

---

[17] J. M. Hermández-Lobato, M. A. Gelbart, R. P. Adams, M. W. Hoffman, and Z. Ghahramani. A general framework for constrained bayesian optimization using information-han search. 2016.
[18] A. Marco, D. Baumann, M. Khadiv, P. Hennig. L. Righetti, and S. Trimee. Robot learning with crash constraints. IEEE Robotics and Automation Letters, ((2):1439-144)) 1446, 2021.
[19] R. Gardner, M. J. Kusner, Z. E. Xu, K. Q. Weinberger, and P. Cunningham. Bayesian optimization with inequality constraints, In ICML, volume 2014, pages 937-945, 2014.
[20] A. Marco, A. von Rohr, D. Baumann, J. M. Heráández-Lobato, and S. Trimpe. Excursion search for constrained bayesian optimization under a limited budget of failures. arXiv preprint arXiv:2005.07443, 2020.
[21] M. Schonlau, W. Welch, and D. R. Jones. Global versus local search in constrained optimiza- tion of computer models. Lecture motes-monograph series, pages 11-25, 1998.
[22] P. I. Frazier. A tutorial on bayesian optimization. arXiv preprint arXiv: 180.0020 2018.
[23] B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, and N. De Freitas. Taking the human out of the loop: A review of bayesian optimization. Proceedings of the IEEE, 104(1):148-177 2015.
[24] J. Kirschner, M. Mutny, N. Hiller, R. Ischebeck, and A. Krause. Adaptive and safe bayesian optimization in high dimensions via one-dimensional subspaces. International Conference on Machine Learning, pages 3429 -438.. PMLR, 2019.
[25] S. De Blasi, A. Neifer, and A. Gepperth. Multi-pronged safe bayesian optimization for high dimensions. In 2021 IEEE International Conference Systems, Man, and Cybernetics (SMC), pages 1976.1973. IEEE, 2021.
[26] R. Turner, D. Eriksson, M. McCourt, J. Kiili, E. Laaksonen, Z. Xu, and I. Guyon. Bayesian optimization is superior to random search for machine learning hyperparameter tuning: Analysis of the black-box optimization challenge 2020. In NeurlPS 2020 Competition and Demonstration Track, pages 3 26. PMLR, 2021.
[27] M. Binois and N. Wycoff. survey on high-dimensional gaussian process modeling with appli- cation to bayesian optimization. ACM Transactions Evolutionary Learning and Optimization, 2(2):1-26, 2022.
[28] B. Chen, R. Castro, and A. Krause. Joint optimization and variable selection of high-dimensional gaussian processes. International Conference on Machine Learning, 2012.
[29] M. Zhang, H. Li, and S. Su. High dimensional bayesian optimization via supervised dimension reduction. arXiv preprint XXii::997000955 2019.
[30] L. Song, K. Xue, X. Huang, and C. Qian. Monte carlo tree search based variable selection for high dimensional bayesian optimization. Advances in Neural Information Processing Systems, (422848-22800) 2022.
[31] Z. Wang, F. Hutter, M. Zoghi, D. Matheson, and N. de Feitas. Bayesian optimization in a billion dimensions via random embeddings. Journal of Artificial Intelligence Research, 55:361-387 2016.
[32] A. Nayebi, A. Munteanu, and M. Poloczek. framework for bayesian optimization in embedded subspaces. In International Conference on Machine Learning, pages 4752-4761. PMLR, 2019.
[33] C. Li, S. Gupta, S. Rana, V. Nguyen, S. Venkatesh, and A. Shilton. High dimensional bayesian optimization using dropout. arXiv preprint arXiv:1802.05-000 2018.
10

---

[34] B. Letham, R. Calandra, A. Rai, and E. Bakshy, Re-examining linear embeddings for high- dimensional bayesian optimization. Advances in neural information processing systems, 33: 1546-1558, 2020.
[35] L.. Papenmeier, L. Nardi, and M. Poloczek. Increasing the scope you learn: Adaptive bayesian optimization in nested subspaces. Advances in Neural Information Processing Systems, 35: 11586-1 1601, 2022.
[36] R. Gómez-Bombarelle N. Wei, D. Duvenaud, J. M. Hernández-Lobato, B. Sánchez-Lengelingg D. Sheberla, J. Aguilera-Iparraguiree T. D. Hirzel, R. P. Adams, and A. Aspuru-Guzik. Auto- matic chemical design using a data-driven continuous representation of molecules. ACS central science, 4(2):268-22) 276, 2018.
[37] R. Moriconi, M. P. Deisenroth, and K. Sesh Kumar. High-dimensional bayesian optimization using fow-dimensions feature spaces, Machine Learning. 109(9):1925-194) 2020.
[38] A. Tripp, E. Daxberger, and M. Hernández-Lobato. Sample-efficient optimization in the latent space of deep generative models via weighted retraining. Advances in Neural Information Processing Systems, 33:1 1 1259-11272 2020.
[39] A. Deshwal and Doppa. Combining latent space and structured kernels for bayesian opti- mization over combinatorial spaces. Advances in Neural Information Processing Systems, 34: 8 3185-8200, 2021.
[40] A. Grosnit, R. Tutunov, M. Maraval, R.-R. Griffiths, I. Cowen-Rivers, Yang. Zhu, W. Lyu, Z. Chen, Wang, al. High-dimensional bayesian optimisation with variational autoencoders and deep metric learning. arXiv preprint arXiv:2106.03609 2021.
[41] E. Siivola, A. Paleyes, J. González, and Vehtari. Good practices for bayesian optimization of high dimensional structured spaces. Applied Al Letters, 2(2):e2,, 2021.
[42] C. Oh, F. Gavves, and M. Welling. Bock: Bayesian optimization with cylindrical kernels. In International Conference on Machine Learning, pages 3868-3877. PMLR, 2018.
[43] D. Eriksson, M. Pearce, J. Gardner, R. D. Turner, and M. Poloczek. Scalable global optimization via local bayesian optimization. Advances in neural information processing systems, 32, 2019.
[44] S. Müller, A. von Rohr, and S. Trimpe. Local policy search with bayesian optimization. Advances in Neural Information Processing Systems, 34:20708-20722 2021.
[45] Q Nguyen, K. Wu, Gardner, and R. Garnett. Local bayesian optimization via maximizing probability of descent. Advances in Neural Information Processing Systems, 15:13190-13202. 2022.
[46] K. Wu, K. Kim, R. Garnett, and R. Gardner. The behavior and convergence of local bayesian optimization. arXiv preprint arXiv:2305.15571 2023.
[47] N. Maus, H. Jones, J. Moore, M. J. Kusner, J. Bradshaw, and J. Gardner. Local latent space bayesian optimization over structured inputs. Advances in Neural Information Processing Systems, 35:34505-34511 34518, 2022.
[48] L. Wang, R. Fonseca, and Y. Tian. Learning search space partition for black-box optimization using monte carlo tree search. Advances in Neural Information Processing Systems, 33:19511 19522, 2020.
[49] N. Hansen. The ema evolution strategy: a comparing review. Towards a new evolutionary computation: Advances in the estimation of distribution algorithms, pages 75-102, 2006.
[50] O. Kramer. A review of constraint-handlngg techniques for evolution strategies. Applied Computational Intelligence and Soft Computing, 2010:1-19, 2010.

---

[51] D. V. Armold and N. Hansen. A (1+ e)-em-ess for constrained optimisation. In Proceedings of the 14th annual conference on Genetic and evolutionary computation, pages 297-304, 2012.
[52] K. Kandasamy, A. Krishnamurthy, J. Schneider, and B. Póczos. Parallelised bayesian opti- misation via thompson sampling. In International Conference on Artificial Intelligence and Statistics, pages 133-142. PMLR, 2018.
[53] N. Srinivas, A. Krause, S. M. Kakade, and M. Seeger. Gaussian process optimization in the bandit setting: No regret and experimental design. arXiv preprint arXiv:0912.3995, 2009,
[54] L.. Yonghyeon, S. Yoon, M. Son, and F. C. Park. Regularized autoencoders for isometric representation learning. In International Conference on Learning Representations, 2021.
[55] C. Zuo, K. He, J. Shao, and Y. Sui. Self model for embodied intelligence: Modeling full- body human musculoskeletal system and locomotion control with hierarchical low-dimensional representation. arXiv preprint arXiv:23/2.05477, 2023.
[56] H. Mania, A. Guy, and B. Recht. Simple random search provides a competitive approach to reinforcement learning. arXiv preprint arXiv:/803.07055 2018.
[57] O. Aydin, X. Zhang, S. Nuethong, G. J. Pagan-Diaz, R. Bashir, M. Gazzola, and M. T. A. Saif. Neuromuscular actuation of biohybrid motile bots. Proceedings of the National Academy of Sciences, 1 116(40);1844111998)
2019.
[58] S. Kurumaya, K. Suzumori, H. Nabae, and S. Wakimoto. Musculoskeletal lower-limb robot driven by multifilament muscles. Robomech Jourmal, 3:1-15, 2016.
[59] S. R. Chowdhury and A. Gopalan. On kernelized multi-armed bandits. In International Conference on Machine Learning, pages 844-853. PMLR, 2017.
[60] D. P. Kingma and Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
[61] M. Balandat, B. Karrer, D. Jiang. S. Daulton, B. Letham, A. G. Wilson, and E. Bakshy, Botorch: A framework for efficient monte-carlo bayesian optimization. Advances in neural information processing systems, 33:21524-21533 2020.
[62] E. Todorov, T. Erez, and Y. Tassa. Mujoco: physics engine for model-based control. In 2012 IEEE/RSJ international conference on intelligent robots and systems, pages 5026-5033, IEEE, 2012.
[63] T. Haarnoja, A. Zhou, P. Abbeel, and S. Levine. Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with stochastic actor. In International conference on machine learning, pages 1861-1870. PMLR, 2018.
[64] A. Thomson. Fifth annual report of the committee of collective investigation of the anatomical society of great britain and ireland for the year 1893-94. Journal of Anatomy, 29 1):35660, 1894.
[65] R. E. Mccotter. Regarding the length and extent of the human medulla spinalis. The Anatomical Record, 10(9):559-566 1916.
[66] M. W. Zhou, W. T. Wang, H. S. Huang, G. Y. Zhu, Y. P. Chen, and C. M. Zhou. Microsurgical anatomy of lumbosacral nerve rootlets for highly selective rhizotomy in chronic spinal cord injury. The Anatomical Record: Advances in Integrative Anatomy and Evolutionary Biology, 293, 2010.
[67] A. Rowald, S. Komi, R. Demesmaeker, E. Baaklini, S. D. Hermandez-Charpak, E. Paoles, H. Montanaro, A. Cassara, F, Becce, B. Lloyd, al. Activity-dependent spinal cord neuromod- ulation rapidly restores trunk and leg motor functions after complete paralysis. Nature medicine, 28(2):260-271, 2022.
12

---

[68] T. Kameyama, Y. Hashizume, and G. Sobue. Morphologic features of the normal human cadaveric spinal cord. Spine, 21 1): 1215, 1996.
[69] J. Ladenbauer, K. Minassian, U. S. Hofstoetter, M. R. Dimitrijevic, and F. Rattay. Stimulation of the human lumbar spinal cord with implanted and surface electrodes: A computer simulation study. IEEE Transactions on Neural Systems Rehabilitation Engineering A Publication of the IEEE Engineering in Medicine Biology Society, 18(6):637, 2010.
[70] W. Sharrard. The segmental innervation of the lower limb muscles in man. Annals of the Royal College of Surgeons of England, 35(2):106, 1964.
[71] C. M. Schirmer, L. Shils, J. E. Arle, R. Cosgrove, and S. Magge. Heuristic map of myotomal innervation in humans using direct intraoperative nerve root stimulation. J Neurosurg Spine, 15(1)664-70, 2011.
Analysis of models for external stimulation of axons. IEEE transactions on [72] Rattay and Frank.
bio-medical engineering, 33(10)9974, 1986.
[73] C. R. Butson and C. C. McIntyre. Role of electrode design on the volume of tissue activated during deep brain stimulation. Journal of Neural Engineering, 3:1-8, 2006. ISSN 1741-2560. doi: 1088//111 04008/7771125500//0000.0
[74] M. W. Seeger, C. K. Williams, and N. D. Lawrence. Fast forward selection to speed up sparse gaussian process regression. In International Workshop Artificial Intelligence and Statistics, pages
PMLR, 2003.
[75] E. Snelson and Z. Ghahramani. Sparse gaussian processes using pseudo-inputs. Advances in neural information processing systems, 18, 2005.
[76] Hensman, N. Fusi, and N. D. Lawrence. Gaussian processes for big data. arXiv preprint arXiv: 13300666355 2013.
[77] K. Kandasamy. J. Schneider, and B. Póczos. High dimensional bayesian optimisation and handits via additive models. In International conference on machine learning, pages 295-304.
PMLR, 2015.
[78] J. Gardner, C. Guo, K. Weinberger, R. Garnett, and R. Grosse. Discovering and exploiting additive structure for bayesian optimization. In Artificial Intelligence and Statistics, pages 1311-1319. PMLR, 2017.
[79] Z. Wang, C. Gehring, Kohli, and Jegelka. Batched large-scale bayesian optimization in nigh-dimensional spaces. In International Conference Artificial Intelligence and Statistics, pages 745-754. PMLR, 2018.
[80] M. Mutny and A. Krause. Efficient high dimensional bayesian optimization with additivity and quadrature fourier features. Advances in Neural Information Processing Systems, 31, 2018.
[81] D. Eriksson and M. Jankowiak. High-dimensional bayesian optimization with sparse axis- aligned subspaces. In Uncertainty Artificial Intelligence, pages 493-503. PMLR, 2021.
13

---

A Theoretical Analysis
A.1 Proof of Proposition 4.1
Proof. Fix t 2 and 2. Conditioned on y.-1 (vf..... ...,8/1)) (z..... - istic, and g(z) ~ N ( ( 0- Fo0r~N(00, 1), Pr(r c) = ((c).
are determin-
Therefore, by applying (g(z) (1-1 and c = 3t, the statement holds.
A.2 Proof of Theorem 4.2
We first introduce the results in Chowdhury and Gopalan [59] which prove the selection of confidence scalar to contain the function with high probability, and the upper bound of cumulative standard deviations.
Lemma A.1. (Lemma 5.1 in Srinivas al. [53]) Let assumptions 3.1 and ?? hold for g. For any 5 E (0, 1), with probability at least 1 5, the following holds for all x E X and 1 I T, EN,
[g(x) pt-1(x:) (8
(5)
where B% 200g(1X()?>>2//)
Then we can bound instantaneous safety violation of HDSAFEBO.
Lemma A.2.
Let assumptions 3. hold for safety function g- With probability at least 1 6, the sample points of Algorithm 1 for all time steps 1 S S TT satisfy
E Ele = E[mmxx(O, - ((x,))) - 2(1 - 1 (z.))
(6)
Proof. We denote [.] max(-;0) 0). Then we have
Ejer] = [[[-g([x)]]]]
(7)
[[(- - u, (x+) + [
(8)
S B[[-g(x;) ) + u (x; + (x))+++
(9)
[[[-9(x;) + (x,))
(10)
+ + [[[-9(x;)
(xx)] ((g(x;)))))))))))))))){())))11) u.(x:)) + (--9(x))) (zo))) (xx))]
(11)
= (1 - a) a){-)(x;, + (x(x))++
(12)
< (1 a)| - (ax)) (89) 0 (xx:) ( (35)
(. (2:,)
(13)
2(1 (3))) 0 (
(14)
where the inequality (9) follows by the fact that [a + b/+ < [a] [b] Ya, b R. the equality (10) is derived from the fact that ux,)) 0, the equality (12) is derived from Proposition 4.1, and the inequality (13) is derived from Lemma A.l.
Lemma A.3. (Theorem I in Srinivas et al. [53].) Let "..". ar be the points selected by the algorithms. With C₁ == 8/ log(1 -2)
T
Σ
2 5 V
(15)
Finally bound the summation of instantancous safety violations.
14

---

Proof.
E V
E
(16)
1=1
E3 [1/1
(17)
T
Σ
2(1 - ()))
(18)
t=1
T
< 2(1 a)
Σ
(19)
1=1
(1-a)VCiTB9T)
(20)
where the inequality (18) follows by Lemma A.2, the inequality (19) follows by the monotonicity of B8, and the inequality (20) follows by Lemma A.3.
A.3 Proof of Proposition 4.3
Proof. Define k* h² are the kernel function of GP with same hyperparameter over X' and Z respectively. For stationary kernel, the kernel function value between and x' is only depend on l x'). therefore we can write the k* k² as a function of metric dx,, dz:
(x, = xxx (dx(x, x²)), (z, = h2 (d2(2z
(21)
We define x- as the kernel function of GP X, and u," the estimated upper confidence bound obtain equivalent GP estimation
of GP over X. Utilizing the property of isometric embedding, between original space and embeded subspace:
(d2 k (dx (x) (z))).
(22)
Therefore the estimated upper confidence bound also equivalent. For every sample z from HDSAFEBO, have
Pr(g(11 (x)) 2 0) 2 Pr(u, (II (II 0) -r (z) 2 0)
(23)
(24)
a,
(25)
the statement holds.
B Experimental Details
The full implementation of our experiments can be found on our project page: https://Insgroup cc/r bbasearchhhdsafee Our musculoskeletal model will be released soon. In the meantime, the model can be accessed for research purposes upon request (ysui ssuiittsinhuua.du.cn))
B.1 PCA training
We employ PCA using the scikit-leamn library with default parameter settings.
'https://scikitt learn. http:///ww//eeeeeee///////////ta wtt.://eww...eecmmomiiio........ html
15

---

B.2 Autoencoder training
We employ IRVAE on neural stimulation task and digital generation task by directly using the paper's original repository". We use MLP as the VAE module for all tasks, and list the model detail in Table We train all models for 300 epochs using Adam[60] optimizer with a learning rate of 0.0001.
Task
Layer number
Hiddien number
Latent dimension
Neural stimulation
4
256
16
Digital generation
4
256
6
Table 2: Autoencoder model detail.
B.3 Algorithm Implementation
For implementation of HDSAFEBO, We use BoTorch as the GP inference component [61]. We also use BoTorrh to replicate I LineBO, SCBO, CONFIG, cEI and cEI-Prob. For LineBO, we choose to implement the random line embedding version shown the main paper. We use the package pycma² 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 Tx; Ty, Ty, lu, hmmx, Imin 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 3 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 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.
B.4 Synthetic Function
We sample objective and safety function from separate Gaussian process with Matérn kernel and length scale as 0.05, which is implemented using GPyTorch'. We set the safety threshold to -0.75. For each independent run, latent space, and randomly select de variables the function effective dimension.
a random linear projection Trand RDxd create d-dimensional
B.5 Musculoskeletal Model Control
We use 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 bottle. At the beginming of the episode, the bottle is initially positioned horizontally in the hand. The initial task involves first rotating the bottle achieve and maintain vertical orientation. At each time step, the reward from the environment is computed follows:
P' = Posse + Pbonus 10 Penalty + 7 Tgrasp + 2 # Tsurvive Factivation
100 # Tarop
(26)
where Tpose is the difference between the bottle and vertical orientation, computed by Euler angle. Tbouus is the reward when the difference falls below a predefined threshold. Tpenalty is positive when the bottle position is out of from the predefined range. r rggasp is the distance between the centroid of
hnttpe://github.c
Caacaaub..oc//MM-BBE /CMMmS/ "nteps://gpytorch.tt
16

---

the bottle and palm joints. Tsurvive is the reward for not dropping the bottle during the current time step. Factivation is the penalty for large muscle activations. 'roo is the penalty when the bottle drop from hand. The overall simulation based on Mujoco [62].
When the height of the bottle is below 0.4m, consider it be dropped from the hand and the episode concludes. We record the speed of episode ending the landing speed of the bottle. We use 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.
B.6 Modeling of the human neuromuscular system
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 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 museles of bilateral lower limbs were studied: iliopsoas (II.), 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:
Meneighbar
SI:
14
(27)
Maeighhor
where mueighbor represents the number of muscles whose motor neuron pools adjacent to the ith muscle's. The selectivity index ranges from to 1, where -1 represents the maximum of activation of all undesired muscles with complete absence of activation of the targeted muscle, O indicates that all muscles are activated at the same level, and means the targeted muscle is activated at the greatest extent while no undesired muscles are activated. And μ4 is the normalized activation of the ith muscle and is defined follows in the simulation:
JfS f(z, v, 2)drdydz 14
(28)
SJJ,, Idrdydz
if 2) Althreshold if AF y, 2) < AFthrestold P2 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 axon([72, 73]).
f1. (0.
fxx y, v,,
(29)
We use the spinal model to traverse all stimulation parameters with 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 g(x) 1 maz;(u;). 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.
17

---

We convert electrode parameters from 17d vector to 2d electric field image using simplified computation. In concrete, we map contact combinations to the spatial position 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 Figure 5.
We the constructed SCS dataset as the function oracle. For given 2d map, set its function value as the function value of the nearest unevaluated point in the dataset measured by the electrical field map.
1031
104
697
-11
-32
11
04
-3
04
SI distriution IL
distribution of RF
distribution of
4400
203
83
252
1301
204
477
412
61
6.2
02 63
0.2
6.1
00
6.1
05
00
04
SI distribution BF
SI distribttio of ST
distriutioon of GA
Figure 4: Distribution of SI for six muscle groups of different configurations used in SCS simulation experiment
27-
53
Figure 5: 2d electrical map computation.
B.7 Clinical experiment of neuromuscular system control
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 the patient while the evoked muscle activities were recorded using EMG. These data (132 trials) were used initialize HDSAFEBO. Except for the first 132 trials the initial data, 441 out of 504 trials recommended by HDSAFEBO. The other trials were conducted by the therapist.
18

---

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:
14
Sh
(30)
where 14 represents the normalized peak-to-peat to-peak value of the evoked EMG for the ith muscle and m the total number of the target muscles.
During the optimization, We restrict the number of cathodes and anodes to only evaluate configura- tions consistent with clinical priors. We used the Electromyographhy (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 onto the stimulator, which would deliver electrical stimulation to the patient for 800 ms a frequency of 10 Hz. Peak-to-peak values averaged and normalized obtain selectivity indices of different muscles after stimu- lation. The calculated feedback and queried safety score were used to update the optimizer and it would recommend new parameter. The tasks were refined sequentially and all the history data were reused when optimizing a new task.
C Additional Experiment
C.1 Baselines performance over embeded subspace
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 high- dimensional for them.
Table 3: Algorithm performance 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.
LinetO
SLBO
CONPIO
Objentee (1)
3,56 & 0.15
51G 013
3,00 000
3.1.1...8 0.06
290 077
23 0.04
29 6.65
280 001
24 0.01
2.T 0.0 0.27 1.0.1
2/00 0.0.
Saley (1)
00
077.180
0.00
Vallation
0.58
19:5
C.2 Ablation on algorithm components
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 over-explore, leading degraded optimization and safety performance in this dighedimmnsional problem. Without optimistic safe identification, the algorithm makes more unsafe selections during the early stages of optimization. Combining these two components enables safe and efficient optimization procedure.
HdSafeBO
HdSafeBO w/o local search
HdSafeBO w/o optimistic safe
1000
1.00
50
800
0.95
40
Objective
Violation
600
Safety
0.90
30
400
0.85
20
200
0.80
10
0.75
200
400
600
800
200
400
600
800
200
400
600
800
Evaluation number
Evaluation number
Evaluation number
Figure Ablation on components HdSafeBO. Optimization performance averaged over 50 independent runs.
19

---

C.3 Distance Preserving of IRVAE
In the implementation of HDSAFEBO, explore the use of IRVAE to learn mapping with distance- preserving property. We randomly sample 10,000 datapoints from SCS dataset. Using conventional VAE and IRVAE trained from same dataset, we compute the point-wise distance as well as GP estimation difference between the original and latent space, as shown Table The results
Using the 218,000 synthetic training data in human neural stimulation control, plotted the pair- wise distance of points in both original input space and latent space learns by IRVAE (Figure 7). We can observe that the learned embedding space exhibits approximately scaled isometry.
Model
Lincar correlation of distance (t)
GP mean estimation difference (3)
GP variance estimation difference ())
IRVAF
0.9729
1.2936
2.9953
conventional
0.8923
2.1839
4.5608
Table 4: Distance preserving comparison between IRVAE and conventional VAE
L0O2 pornpo dranne
prgeel ditance
Figure Pair-wise Distance training data points neural stimulation induced human motion control task.
C.4 Ablation on Confidence Bound Scalar
Here HDSAFEBO with 0, 16 musculoskeletal model control task, and shown the results in Figure 8. The safety and violation metric may become slightly worse 3 increases. We observe the algorithm performance is similar under a wide range choice of /3 on this high-dimensional task.
tebrO
bebr
253
-000
10
ND D8
(a)
(b)
(c)
Figure 8: Ablation study on confidence bound scalar J3.
C.5 Comparison with random-embedding BO
We additionally run HesBO and BAxUS on the musculoskeletal model control task and the neural stimulation task with no safety constraints in simulation. Due to the algorithmic mechanism of HesBO and BAxUS, we cannot directly use the same initial point as HDSAFEBO. Therefore we randomly sample initial points from their corresponding latent space. In HesBO, set the same latent dimension number as in HDSAFEBO. Table 5 shows the best objective function values found by algorithms (shown mean + std).
20

---

Table 5: Best objective function values found by algorithms.
Algorithm
SCS-ST
SCS-GA
HDSAFEBO
0.26 + 0.02
0.22 + 0.01
HesBO
0.23 + 0.0
0.19 0.02
BAXUS
0.24 + 0.03
0.19 0.05
We observe HDSAFEBO still outperforms HesBO and BAxUS across all tasks, even when optimizing under safety constraint. We think using IRVAE enables utilizing the pre-collected unlabeled data to learn better representation than random projection.
C.6 Hand-writing Digital Generation
In hand-writing digital generation task, the goal is to generate images of target digit thick as possible, while keeping the image valid for the required number. Using this task we can test the algorithm performances when the latent dimension is low. We trained fully-connected IRVAE with a latent dimension as the sum of image pixel intensities and the safety function as the prediction probability of target number. Since the CNN prediction is very sharp, we wrap the CNN output via softmax layer with temperature as 200. We set the sample budget 200 including with 20 images of target digit the initial data.
and use a two-layer CNN model as the predictor. We set the objective function
We summarize the averaging performance Table 6. We observed HDSAFEBO outperforms all baselines in terms of optimization performance and safety violation. Note that while HDSAFEBO effi- ciently finds highest objective, its safety violations is 63% less than the second best method original SCBO.
Table Experiment results of constrained hand-writing digital generation. We evaluate algorithm performance generating digital from to in terms of best found feasible objective value (higher is better) and cumulative safety violation (lower is better). Objective values are normalized by best feasible point in the MNIST dataset. The results are shown mean performance + one standard deviation across ten tasks.
Metric
HDSAFEBO
SCBO
CONFIG
cEI
CMAES
Objective
1.14
0.12
1.08 0.18
0.77
0.26
0.77 0.26
0.68
0.12
Violation
18.99 + 10.82
52.01
20.28
72.09 17.85
72.04 + 17.82
72.01 + 7.06
C.7 Run time performance
We HdSafeBO and baseline algorithms the human neural stimulation control problem and recorded the average run-time of each iteration in Table 7. Most baseline algorithms have processing time of around 10 seconds.
In our problem setting, we consider the run-time difference to be marginal compared to the actual experiment time for each trial. For instance, in our real experiments, applying recommended stimu- lation parameters typically takes 1-2 minutes. Therefore, we compared optimization performances based on the number of evaluations.
Method
HdSafeBO
LineBO
SCBO
CONFIG
cEI
CMAES
Iteration time
11.75 + 1.35s
9.88 1.01s
11.99
1.21s
9.49 + L.05s 20.38 1 1.21s
9.59 1: 0.72s
Table 7: Run-time per iteration in neural induced human motion control task. Results show mean 1. 1 std
averaged over 80 iterations.
21

---

D Additional Related Work
Here we additionally discuss more works about high-dimensional Bayesian optimization besides dimension-reduucion based BO and local BO.
Due to the inversion of kernel function matrix, the complexity of GP inference scales exponentially with the sample number, limiting the search budget of high-dimensional problems. Sparse GP or
variational GP is used to achieve scalable sampling over the high-dimensional space [74, 75, 76].
Another line of work assume the addictive structure the objective functions, and decomposes the function to solve the low-dimensional sub- troblrrobem decentrally(77, 78, 79, 80].
To overcome over-exploration issue over the high-dimensional space, several works also propose to shape kernel prior to sample points near the search center [42, 81].
22


Reviews:
[{"review_id": "WaNUqQursc", "reviewer": "Reviewer_7rTb", "summary": "", "questions": "", "limitations": "", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": 3, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "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.", "summary_of_recommendation": "I think the paper should be accepted at the conference if there is no significant algorithmic overlap with prior work.", "technical_quality": 4, "clarity_of_presentation": 5, "potential_impact": 3, "robotics_focus": 3, "quality_of_limitations": 3, "questions_for_rebuttal": "Sea weaknesses and comments section above", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "5UsmoxFJOu", "reviewer": "Reviewer_KPXx", "summary": "", "questions": "", "limitations": "", "rating": 2, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": 2, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "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.", "summary_of_recommendation": "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.", "technical_quality": 2, "clarity_of_presentation": 2, "potential_impact": 2, "robotics_focus": 4, "quality_of_limitations": 2, "questions_for_rebuttal": "Please see review", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "mF5IJd1AZG", "reviewer": "Reviewer_kViv", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": 3, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "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).", "summary_of_recommendation": "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.", "technical_quality": 3, "clarity_of_presentation": 2, "potential_impact": 3, "robotics_focus": 3, "quality_of_limitations": 2, "questions_for_rebuttal": "-\tIn 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. \n-\tAdditionally, 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.  \n-\tWhat 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?\n-\tRegarding 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?\n-\tIn 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?\n-\tHow 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?\n-\tConcerning 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. \n\n**Small comments**\n\n-\tProblem 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. \n-\tL. 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. \n-\tAlgorithm 1: line 13, please check type setting again, arrow is weird\n-\tL. 227: what is the observation space? (num of muscles + num of joints?) \n-\tL. 266 should be semitendinosus (s missing)", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: 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.
Decision: Accept