Sample-efficient Bayesian Optimisation Using Known Invariances

Authors: theodore  brown
alexandru  cioba
ilija  bogunovic
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-02 15:06:20

Abstract: Bayesian optimisation (BO) is a powerful framework for global optimisation of costly functions, using predictions from Gaussian process models (GPs). In this work, we apply BO to functions that exhibit invariance to a known group of transformations.  We show that vanilla and constrained BO algorithms are inefficient when optimising such invariant objectives, and provide a method for incorporating group invariances into the kernel of the GP to produce invariance-aware algorithms that achieve significant improvements in sample efficiency. We derive a bound on the maximum information gain of these invariant kernels, and provide novel upper and lower bounds on the number of observations required for invariance-aware BO algorithms to achieve $\epsilon$-optimality. We demonstrate our method's improved performance on a range of synthetic invariant and quasi-invariant functions. We also apply our method in the case where only some of the invariance is incorporated into the kernel, and find that these kernels achieve similar gains in sample efficiency at significantly reduced computational cost. Finally, we use invariant BO to design a current drive system for a nuclear fusion reactor, finding a high-performance solution where non-invariant methods failed.

Main paper: 
 introduction: 
 Bayesian optimisation (BO) has been applied to many global optimisation tasks in science and engineering, such as improving the performance of web servers [30] , reducing the losses in particle accelerators [26] , or maximising the potency of a drug [1] . The most challenging real-world tasks tend to have large or continuous input domains, as well as target functions that are costly to evaluate.One particular example of interest is the design and operation of nuclear fusion reactors, which promise to deliver huge amounts of energy with zero greenhouse gas emissions. However, the development of high-performance fusion power plant scenarios is a difficult optimisation task, involving expensive simulations and physical experiments [39, 11] . For tasks like this, it is desirable to develop sample efficient algorithms that only require a small number of evaluations to find the optimal solution.Crucially, the physical world is characterized by underlying structures that can be leveraged to significantly improve sample efficiency. The concept of invariance, which describes properties that remain constant under various transformations, is one of the most fundamental structures in nature. A prime example of invariance is found in image classification tasks: a classifier should consistently identify a cat as a cat, regardless of its orientation, position, or size. The central question we answer is the extent to which the sample efficiency of BO can be improved by exploiting these invariances.We study a general approach for incorporating invariances into any BO algorithm. We represent invariant objective functions as elements of a reproducing kernel Hilbert space (RKHS), and produce invariance-aware versions of existing Gaussian process (GP) bandit optimisation algorithms by using the kernel of this RKHS. We apply these algorithms to synthetic and nuclear fusion optimisation tasks, and show that incorporating the underlying invariance into the kernel significantly improves the sample efficiency. We also investigate the behaviour of the algorithms when the full invariant kernel is not known, and provide an empirical study of the performance of invariance-aware BO on target functions that are 'quasi-invariant' (modelled as a sum of invariant and non-invariant components).Our main contribution is the derivation of novel upper and lower bounds on the regret achieved by invariance-aware BO algorithms.Invariance in deep learning. The role of invariance in deep learning has been extensively explored in the literature. Invariance, particularly to geometric transformations like rotation, scaling, and translation, is crucial for the robust performance of deep learning models. Many deep learning methods incorporate invariance via data augmentation, in which additional training examples are generated by transforming the input data [28, 29] . However, the cost of kernelised learning methods often scales poorly with the number of training points n -Gaussian process inference, for example, incurs O(n 3 ) time and O(n 2 ) memory cost [47] -which means that data augmentation is prohibitively expensive. In this work, we adopt an alternative approach by directly embedding invariance properties into the model's structure, which effectively reduces computational complexity [44, 15] .There is an extensive body of empirical BO literature that exploits known patterns present in the behaviour of a system. In [24] , BO is applied after a variety of physics-informed feature transformations to the state; in [48] the authors employ BO to handle optimization over crystal configurations with symmetries found using [40] . However, in these works, only empirical results are given and asymptotic rates on sample complexity are not estimated.Invariant kernel methods. The theory of group invariant kernels was introduced in [18, 27] . Using these ideas, [46] developed efficient variational methods for learning hyperparameters in invariant Gaussian processes. Invariant kernel ridge regression (KRR) for finite groups has been studied for translation-and permutation-invariant kernels [33, 4] . [38] quantify the gain in sample complexity for kernels invariant to the actions of Lie groups on arbitrary manifolds and tie the sample complexity bounds to the dimensions of the group. Motivated by the invariance properties of convolutional neural networks, [27, 33, 4] focus on kernelised regression with the Neural Tangent Kernel; we instead consider the kernelised bandit optimisation setting with the more general Matérn family of kernels. These methods can be readily extended to less strict notions of so-called quasi-invariance [18, 46, 4] .Matérn kernels on manifolds. The previous studies of invariant kernel algorithms rely on the spectral properties of kernels on the hypersphere [33, 4] ; the relevant spectral theory for Matérn kernels on Riemannian manifolds (including the hypersphere) was developed in [9] . These geometry-aware kernels were applied to optimisation tasks in [20] ; however, the structure they include is restricted to the geometry of the domain only, rather than invariance to a set of transformations.Regret bounds for Bayesian optimisation. The techniques for upper-bounding the regret of Bayesian optimisation algorithms using the kernel-dependent quantity known as maximum information gain were developed in [37] . Subsequent studies [43, 8, 14, 36, 19, 41] have focused on improving regret bounds in both cumulative and simple regret settings. [22] provided an upper bound for regret in BO using the convolutional NTK, which exhibits invariance only to the cyclic group. [23] explored graph neural network bandits, offering maximum information gain rates by embedding permutation invariance using additive kernels. In [25] , a permutation-invariant kernel is constructed in the same way as ours, and a regret analysis provided. Our work expands upon these approaches by considering more general (non-additive) kernels and broader transformation groups, demonstrating that the regret bounds from [22] can be shown to be a specialisation of ours. Lower bounds for kernelised bandits, including their robust variants ( [7, 6] ), have been examined in [35] and [12] . In this paper, we derive novel lower bounds on regret for the setting with invariant kernels.Our main contributions are:• We develop new upper bounds on sample complexity in Bayesian optimisation with invariant kernels, which explicitly quantify the gain in sample complexity due to the number of symmetries.We extend several results from the literature, as our statement applies to a wide class of groups. This is, to our knowledge, the first time such bounds are given in this degree of generality. • We show how results from the literature (e.g. [12] ) can be extended to obtain lower bounds on the sample complexity of invariant BO, using a novel construction and quantification of members of the RKHS of invariant functions. We extend the methods in the literature to the hyperspherical domain, giving a blueprint for similar constructions on other Riemannian manifolds. • We conduct experiments that support our theoretical findings. We demonstrate that incorporating invariance to some subgroup of transformations is sufficient to achieve good performance, a result of key importance when the full invariance is not known or is too expensive to compute.We also demonstrate the robustness of the invariant kernel method, in the case where the target function deviates from true invariance. • We apply these methods to solve a real-world design problem from nuclear fusion engineering that is known to exhibit invariance; our invariance-aware method finds better solutions in fewer samples compared to standard approaches. 
problem statement: 
 We begin by providing an overview of invariant function spaces. Consider a finite subgroup of isometries, G, on a compact manifold, X , embedded inEQUATIONWe provide some examples of invariant functions for different groups G in Figure 1 .In this paper, we study the optimisation of invariant functions that belong to a reproducing kernel Hilbert space (RKHS). The invariance of the RKHS elements translates to invariance properties of the kernel. For k : X × X → R, we say that k is simultaneously invariant to the action of G if, for all σ ∈ G and x, y ∈ X ,EQUATIONand is totally invariant to the action of G if, for all σ, τ ∈ G and x, y ∈ X ,EQUATIONas introduced in [18] . Note that both inner product kernels k(x, y) = κ (⟨x, y⟩) and stationary kernels k(x, y) = κ(∥x -y∥) satisfy simultaneous invariance, as G consists of isometries (which preserve distances and inner products).Next, we outline the relationship between the invariance of the kernel and the invariance of the functions in the associated RKHS. In particular, the following proposition (Proposition 1) identifies the space of invariant functions as a subspace of the RKHS of a simultaneously invariant kernel. Moreover, it turns out that this subspace is an RKHS in its own right, with a corresponding kernel that is totally invariant. For a proof of these properties, see Appendix A.1. Proposition 1 (RKHS of invariant functions). Let X and G be as defined above. Consider a positive definite and simultaneously invariant kernel k : X × X → R. The kernel defines a reproducing kernel Hilbert space H k , whose elements are functions f : X → R, and a corresponding Hilbert space norm ∥ • ∥ H k . Define the symmetrization operator S G :EQUATION)Figure 2 : UCB run with standard k and invariant kernel k G applied to an invariant objective; queried points in white. UCB with the invariant kernel requires significantly fewer samples to find the optimum.Then, S G is a self-adjoint projection operator, whose image Im[S G ] is the subspace of H k that contains G-invariant functions. Moreoever, Im[S G ] is an RKHS in its own right, H k G . For x, y ∈ X , the reproducing kernel ofEQUATIONWe consider a kernelised optimisation setup with bandit feedback, in which the goal is to find the maximum of an unknown function f : X → R that is invariant to a known group of transformations G.We model f as belonging to the invariant RKHS H k G from Proposition 1. Without further smoothness assumptions, this task is intractable; similarly to the common non-invariant setting, we make the regularity assumption that f has bounded RKHS norm, ∥f ∥ H k G < B. The learner is given k and G. In Section 3.2, we explore setting k in k G to common kernels, such as those from the Matérn family.At each round t, the learner selects a point x t ∈ X and receives a noisy function observationEQUATIONwhere η t is drawn independently from a σ n -sub-Gaussian noise distribution. The agent reports a point x * t and incurs a corresponding simple regret, defined asEQUATIONAfter T rounds, the cumulative regret is given byEQUATIONIn the next section, we quantify the impact of invariance on the number of function samples required to achieve a given regret, known as the sample complexity. 
sample complexity bounds for invariance-aware bayesian optimisation: 
 The bandit problem outlined in the previous section can be tackled using Bayesian optimisation (BO) [17] . Intuitively, we expect BO algorithms that utilize the G-invariance of the target function to achieve improvements in sample complexity, as each observed x simultaneously provides information about additional transformed locations σ(x) for each σ ∈ G (see Figure 1 ).In this section, we quantify the performance improvements by deriving bounds on the sample complexity of BO with totally invariant kernels (Equation ( 3)). We begin with an outline of the BO framework.In Bayesian optimisation, predictions from a Gaussian process (GP) are used in conjunction with an acquisition function to select points to query [17] . After collecting a sequence of input-observation pairs, D = {(x 1 , y 1 ), . . . , (x t , y t )}, the learner models the invariant target function f ∈ H k G using a GP with a totally invariant kernel, GP(0, k G ). The posterior predictive distribution of the function value at a specified input location x will be Gaussian with mean and variance given byEQUATIONEQUATIONAlgorithm 1 Max. variance reduction [42] Input: k G , X 1: initialize D = {} 2: for t = 1, . . . , T do 3:Select x t = arg max x∈X σ t-1 (x) 4: Observe y t = f (x t ) + η t 5: Append (x t , y t ) to D 6:Update µ t , σ t from Eqs. 9, 10 7: end for 8: return xMVR opt = arg max x∈X µ T (x) Algorithm 2 Upper confidence bound [37] Input: k G , X 1: initialize D = {} 2: for t = 1, . . . , T do 3: whereSelect x t = arg max x∈X µ t-1 (x)+βσ t-1 (x)k G = [k G (x, x 1 ), . . . , k G (x, x t )] T ∈ R t×1 , K G is the kernel matrix with elements [K G ] i,j = k G (x i , x j ) ∀i, j ∈ 1, . . . , t, y = [y 1 , . . . , y t ]T , and λ is a regularising term. Note that as this GP has a totally G-invariant kernel, k G , it can be interpreted as a distribution over G-invariant functions.We consider two common acquisition functions. The tightest simple regret bounds in the literature are for the Maximum Variance Reduction algorithm (MVR), in which the learner iteratively queries the points with highest uncertainty before reporting the point with highest posterior mean [42] .Another algorithm that has been widely studied is Upper Confidence Bound (UCB), which balances exploration and exploitation by behaving optimistically in the face of uncertainty [37] . Pseudocode for MVR and UCB is provided in Algorithm 1 and Algorithm 2, respectively. Our analysis is framed in terms of MVR and simple regret, although in Section 4.1 we also investigate empirically the cumulative regret performance of UCB.A popular family of kernels used in BO is the Matérn kernel. On R d , the Matérn kernel with parameters l, ν is given byEQUATIONwhere K ν is the modified Bessel function of the second kind. The parameter ν controls the differentiability of the corresponding GP, which enables Matérn kernels to define a diverse range of function spaces. As the Matérn kernel k is simultaneously invariant (Equation ( 2)), we can construct a totally invariant k G and corresponding RKHS following Proposition 1. The bounds that we derive both hold for the Matérn family; our upper bound also extends to any kernel that exhibits polynomial eigendecay.In Figure 2 , we provide a visual example of the impact of using k G in Bayesian optimisation. We run UCB on a group-invariant objective with the standard and invariant kernels. With the non-invariant kernel k, the algorithm spends many samples querying suboptimal points with similar objective values. In contrast, using the invariant kernel k G allows the algorithm to identify suboptimal regions across the parameter space with only a handful of samples, resulting in significantly faster convergence.In this section, let X = S d-1 ⊂ R d and let G be a finite subgroup of the orthogonal group, O(d), equipped with its natural action on X . Let k be a simultaneously invariant kernel on S d-1 , and assume that the eigenvalues of the associated integral operator decay at a polynomial rateµ k = O(k -β * p ), for some constant β * p > 1.Many common kernels exhibit this property, such as the standard Matérn kernel on R d . For GPs defined intrinsically on Riemannian manifolds [45] , the corresponding Matérn kernel also exhibits polynomial eigendecay [9] . We impose the following assumption on the spectra of elements in the group G, which is in line with the assumptions and spectral decay rates in [4] . Assumption 1. All elements σ ∈ G satisfy the following spectral property:EQUATIONfor all λ ∈ Spec(σ), where m λ is the multiplicity of λ.This assumption excludes certain groups in low dimensions, but we expect it can be relaxed through tighter analysis. In fact, our experiments for d = 2, 3 show it is essential as a theoretical artefact only.Nevertheless, Equation ( 12) holds for all finite orthogonal groups in dimension d ≥ 10 (see Lemma 1) as well as for many groups of lower dimension; for example, groups for which m λ = 1, so that only distinct eigenvalues are present in the spectra. Under this assumption, we have the following upper bound on sample complexity. EQUATION)where Õ(•) hides logarithmic factors.Moreover, for the MVR algorithm with invariant kernel k G to achieve expected simple regret E[r T ] < ϵ, it must hold thatEQUATION)For the Matérn kernel, this corresponds toEQUATIONProof. See Appendix A.2.Our upper bound demonstrates that incorporating the group invariance into the kernel is guaranteed to result in sample efficiency improvements that increase with the size of the group, due to the 1 |G| factor. Our result for the MIG γ of invariant kernels is of independent interest, as it is applicable to a wide range of scenarios. Notably, the linear 1 |G| improvement in the information gain from Equation (13) resembles the linear 1 |G| improvement in generalisation error in the invariant kernel ridge regression setting [3] . Intuitively, this reflects how the MIG criterion is determined solely by the ability of the kernel to approximate the objective function, and is agnostic to the algorithm used in optimisation.In our proof, we begin with kernels on S d-1 that are simultaneously invariant w.r.t. the full orthogonal group O(d), and use Proposition 1 to construct kernels that totally invariant kernels w.r.t finite subgroups G ≤ O(d). Simultaneously O(d)-invariant kernels on S d-1 are in fact dot-product kernels, which allows us to use standard versions of Mercer's theorem in terms of known spectral decompositions. If an alternative decomposition is known for a particular k, then it is sufficient that k instead satisfies the weaker condition of simultaneous invariance w.r.t. G only.In the next theorem, we present a lower bound on sample complexity which nearly matches the upper bound in its dependence on |G|. Theorem 2 (Lower bound on sample complexity of invariance-aware BO with Matérn kernels). Let X be one of [0, 1] d or S d . Let k be the standard Matérn kernel on R d with smoothness ν, and G be a finite group of isometries acting on X . Let k G be the totally invariant kernel constructed from k and G according to Proposition 1. Consider noisy bandit optimisation of f ∈ H k G , f : X → R, where the observations are corrupted by σ-sub-Gaussian noise. Then, for sufficiently small ϵ B ,there exists a function family FX (ϵ, B) such that the number of samples required for an algorithm to achieve simple regret less than ϵ for all functions in FX with probability at least 1 -δ is bounded by:EQUATIONProof. See Appendix A.3.Our lower bound guarantees that algorithms that are run for a number of sample less than T fail to be ϵ-optimal with high probability. At fixed T , there is a trade-off between this probability and the order of the group. However, when combined with the upper bound our results imply that as the order of the group increases, the probability that an algorithm is ϵ-optimal after a fixed number of iterations also increases.We recover the previously known lower and upper bounds on sample complexity for non-invariant optimisation from [12] and [42] respectively by setting |G| = 1. We also note that there is a gap between our bounds, as the exponent of1 |G| in Equation ( 15) is1+ d-1 ν but is 1+ d-12ν in Equation ( 16). A likely reason for this gap is in the lower bound, which uses a packing argument for the support sets of functions in the RKHS. While in the absence of the group action, the lower bound and upper bound agree, it is not clear a priori how to achieve the tightest packing of invariant functions.To achieve the largest improvement in sample complexity, one would like to choose the largest possible G when constructing k G . However, as we will see in our experimental section, there is a computational tradeoff between the sample complexity and the computation of S G , as in the case of very large groups the sums over the group elements becomes a computational bottleneck. 
experiments: 
 We demonstrate the performance gains from incorporating known invariances in a number of synthetic and real-world optimisation tasks. Previous work has demonstrated the sample efficiency improvements from applying invariance-aware algorithms to regression tasks [3, 33, 38] ; our experiments complement the literature by observing the corresponding improvements in optimisation tasks. Code for our experiments and implementations of invariant kernels can be found on GitHub 1 .For each of our synthetic experiments, our goal is to find the maximum of a function drawn from an RKHS with a group-invariant isotropic Matérn-5 /2 kernel. The groups we consider act by permutation of the coordinate axes, and therefore include reflections and discrete rotations. For PermInv-2D (Figure 1a ) and PermInv-6D, the kernel of the GP is invariant to the full permutation group which acts on R 2 and R 6 respectively by permuting the coordinates. For CyclInv-3D it is invariant to the cyclic group acting on R 3 by permuting the axes cyclically (Figure 1b ). We provide the learner with the true kernel and observation noise variance to eliminate the effect of online hyperparameter learning. For more detail on our experiments, see Appendix B.In Figure 3 , we show that the algorithm's performance on invariant functions is improved by using an invariant kernel. In all cases, the performance of the invariance-aware algorithm significantly outperforms the baselines, converging to the optimum with many fewer samples. Notably, in the PermInv-6D task, the baseline algorithms fail to find the true optimum, whereas the permutationinvariant versions do so without difficulty. We also note that the performance improvement increases with increasing dimension and group size, as predicted by our regret bounds.For the PermInv tasks, we provide an additional baseline: constrained Bayesian optimisation (CBO, [16] ). In this case, we constrain the search space of the acquisition function optimiser to the fundamental region of the group's action. For arbitrary groups, finding an analytical formula for this region can be difficult, particularly when the dimension of the domain is high and the group is complicated; however, it is straightforward in the case of the full permutation group. We find that constrained BO performs worse than our method, and in high dimensions is identical to standard (non-invariant) BO. Intuitively, this matches our expectations, as in CBO the kernel is not aware of all of the structure that the RKHS elements possess -that is, while the domain of exploration is restricted, the function class is not. Observations in CBO contribute no information across the boundaries of the region, unlike in the invariant kernel case, and so the performance improvement of CBO is upper bounded by 1 |G| . For the PermInv-6D task, we also consider the case where the full invariance is not fully known to the learner; that is, the kernel is not totally invariant w.r.t. the full group, but only a subgroup. Our results show that incorporating any amount of the true invariance results in performance improvements. In this case, the true function is invariant to the full permutation group in 6D (720 elements). We evaluate Non-invariant kernels (blue) are outperformed by the full group invariant kernel (red) as well as partially specified (subgroup-invariant) kernels (green and yellow). For the permutation invariant function, the search space of standard BO can be constrained by the fundamental domain of the group (purple), but this performs worse than the invariant kernel. Mean ± s.d., 32 seeds. the performance of the full invariant kernel, the non-invariant kernel, and two subgroup invariant kernels. The subgroups we consider consist of 'block' permutations, which involve reordering the coordinate indices in chunks of length 2 and 3 (leading to groups with 6 and 2 elements respectively). For both algorithms, we see that all kernels that incorporate invariance outperform the non-invariant kernel, even when the difference in the group size considered is large.In real-world scenarios, it is possible that the underlying function only exhibits approximate symmetry. For example, in the fusion application from Section 4.3, the launchers may have slightly different properties and are no longer interchangeable. A 2D example of this 'quasi-invariance' is given in the first column of Figure 4 . Although the large-scale features (such as the layout of the peaks) remain approximately symmetric, invariance is not maintained in the detail (the values and shapes of the peaks).We model quasi-invariance by introducing an additive non-invariant disturbance to the target function. The function can then be considered as belonging to the RKHSH (1-ε)k G +εk ′ , where k G is a G- invariant kernel, k ′ isa non-invariant kernel, and ε sets the degree of deviation from invariance. In our experiments, we observe that performing BO with the fully invariant kernel still results in significant performance improvements over the non-invariant kernel when ε is small, achieving performance comparable to BO with the kernel of the function's RKHS. As this is a kind of model misspecification, we refer the reader to [5] for a theoretical discussion on the performance of BO with misspecified kernels.In this section, we use invariance-aware BO to efficiently solve a design task from fusion engineering.The leading concept for a controlled fusion power plant is the tokamak, a device that confines a hightemperature plasma using powerful magnetic fields. Finding steady-state configurations for tokamak plasma that are simultaneously robustly stable and high-performance is a challenging task. Moreover, evaluating the plasma operating point requires high-fidelity multi-physics simulations, which take several hours to evaluate one configuration due to the complexity and characteristic timescales of the physics involved [34] . Iterating on a design requires many such simulations, which has led to an increased focus on highly sample-efficient algorithms for tokamak design optimisation [13, 32, 10, 31] .We consider a scalarised version of the multi-objective optimisation task presented in [10] , in which the goal is to shape the output of a plasma current actuator to maximise various performance and ro- bustness criteria based on the so-called safety factor profile. Using a weighted sum, we reduce the vector objective to [0, 1], with component weights selected to represent the proposed SPR45 design of the STEP tokamak [39, 31] . As the objective components are in direct competition, the highest achievable scores are around 0.7 -0.8, while the lowest-scoring converged solutions achieve around 0.4 -0.5.In previous work, the actuator output was represented by arbitrary 1D functions (piecewise linear [31] and Bézier curves [10] ). However, in practice the actuator will have a finite number of 'launchers' that drive current in a localised region, and so will be unable to accurately recreate arbitrary profiles. In contrast to previous work, we directly optimise a parameterisation that reflects the actuator's real capabilities: a sum of 12 Gaussians, where each models a launcher targeted at a different point in the plasma cross-section (Figure 5a ). In this setting, the objective function is f : R 12 → R, where the input is the normalised radial coordinate of the 12 launchers. As all of the launchers have identical behaviour, the objective is invariant to permutations. However, the full invariant kernel is too costly to compute (|G| = 12! = 479 × 10 6 ), so we instead use a partially specified kernel (3-block invariant, |G| = 4! = 24).An additional challenge of this task is that large regions of parameter space are infeasible, corresponding to simulations that will not converge to a steady-state solution and will therefore be impossible to observe. Selecting an appropriate method of dealing with this situation can impact the optimiser's performance, and is an open area of research. In this work, we choose to assign a fixed low score to failed observations for simplicity. As MVR is purely exploratory it has no mechanism that discourages querying points from infeasible regions, leading to many failed observations. UCB is a more natural choice in this setting, as the β parameter introduces a balance between querying unknown regions (that could be infeasible) and exploiting known high-scoring regions.In Figure 5b , we see that the invariant version of UCB achieves significantly improved performance compared to the non-invariant version. In fact, the non-invariant algorithm totally fails to find the small region of high-performance solutions, and instead queries many suboptimal profiles that are similar under the action of permutations. The high sample efficiency of the invariant kernel method will enable tokamak designers to iterate on designs faster, with fewer calls to expensive simulations. 
limitations: 
 Although the number of terms involved in computing the totally invariant kernel only scales linearly with the group size (Equation ( 5)), this cost can still become a bottleneck for very large groups; for example, the size of the permutation group scales with the factorial of the dimension, leading to |G| = 24 for d = 4, but |G| = 479 × 10 6 for d = 12. In Figure 6 we compare the compute and memory costs of performing a marginal log likelihood fit of the GP hyperparameters for each of the kernels in the previous section, given 100 observations from the PermInv-6D function. The standard, 3-block invariant, and 2-block invariant kernel are faster than the fully invariant kernel by a factor of 2. As we observed that the 2and 3-block invariant kernels still achieve significantly improved sample efficiency over the standard kernel (Figure 3 ), we propose that subgroup approximations of the full invariance group can be used as a low-cost alternative when the group size is large, reducing computational cost while maintaining improved performance. Finally, we note that our invariant kernel method scales significantly more favourably than data augmentation, both in terms of memory and compute.Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes]Justification: We provide a running codebase with scripted experiment files. All will be made available online upon acceptance. We provide access to the relevant simulation environments where needed, otherwise, as our algorithms fundamentally deal with data acquisition, open access to data is not applicable. In particular, as the fusion experiment depends on proprietary code and models that belong to the United Kingdom Atomic Energy Authority, the full fusion experiment cannot be made public. We do, however, provide the optimisation script, with the calls to the proprietary models censored, so that reviewers can view our exact methodology.Guidelines:• The answer NA means that paper does not include experiments requiring code.• Please see the NeurIPS code and data submission guidelines (https://nips.cc/ public/guides/CodeSubmissionPolicy) for more details. • While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark). • The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: //nips.cc/public/guides/CodeSubmissionPolicy) for more details. • The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc. • The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only a subset of experiments are reproducible, they should state which ones are omitted from the script and why. • At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable). • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
conclusion: 
 We have developed new upper bounds on the sample complexity in Bayesian optimization with invariant kernels, highlighting the advantages gained from symmetries that hold for a broad range of practical groups. We also derived novel lower bounds by constructing and quantifying members of the RKHS of invariant functions, extending this approach to the hyperspherical domain and providing a framework for other Riemannian manifolds. Empirically, our experiments show that even partial invariance can effectively improve performance, which is crucial when full invariance is unavailable or computationally intensive. Additionally, we applied these methods to a real-world problem in nuclear fusion engineering, achieving superior outcomes with fewer samples compared to standard methods.Broader Impact. The incorporation of symmetries in Bayesian optimization can have a profound impact on science and engineering. By exploiting symmetrical properties in various applications, such as nuclear engineering, material and drug design, and robotic control, it can reduce the number of samples needed to achieve (near) optimal solutions, thus accelerating research and development processes.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes]Justification: All the special experimental settings needed to understand, appraise and judge the usefulness of our results have been provided, and connections to our theoretical setting have been drawn. Details needed for the implementation only are part of the code we supply; this includes exact hyperparameters and seeds used to generate each of the plots.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments? Answer: [Yes] Justification: In all our experiments we provide adequate statistics except where not applicable. This includes sample sizes for experimental runs, sources of randomness and the appropriate error bars where applicable. In the fusion experiment, providing estimates of the standard deviation is not meaningful due to the presence of failed simulations, as discussed in the text. Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments? Answer: [Yes] Justification: We make special mention of the hardware we use, the computational trade-offs in various applications of our method, with quantitative data summarized in a graph. Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?• Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: A.1 Proof of Proposition 1 (RKHS of invariant functions)The kernel k defines a reproducing kernel Hilbert space, H k , with inner product ⟨•, •⟩ H k and norm ∥f ∥ H k = ⟨f, f ⟩ H k for f ∈ H k .Recall from Equation ( 4) that S G (f ) = 1 |G| σ∈G f • σ. We begin by showing that S G is continuous and bounded. Observe that simultaneous invariance determines the following property of the representer elements kEQUATION)Now observe that:EQUATIONIf f = n i=1 k xi , we have, similarly:||f • σ|| 2 = || n i=1 k σ -1 (xi) || 2 = n i=1 n j=1 k(σ -1 (x i ), σ -1 (x j )) = n i=1 n j=1 k(x i , x j ) = ||f || 2 (19)Now, let f n → f be a Cauchy sequence, where each f n lies in the span of the k x elements. We haveEQUATIONTherefore,EQUATIONso S G is continuous, with ||S G || ≤ 1.We now verify that the image of S G , denoted as Im(S G ), contains precisely the G-invariant functions (Equation ( 1)). For f ∈ H k , σ ∈ G, we haveEQUATIONwhere Equation (23) follows by associativity of G, and Equation ( 24) follows by letting σ ′ := τ • σ and using the fact that G is closed under the• operator. Hence, S G (f ) is invariant to σ ∈ G (Equation (1)).Next, we show that the image of S G is a reproducing kernel Hilbert space with kernel k G (Equation ( 5)). First note that S G is a projection, as applying S G twice givesEQUATIONEQUATIONwhere Equation (27) follows by applying Equation (25) .Moreover, S G is self-adjoint. Let f ∈ H k be an element in the span of the representer elements k x , such that f (EQUATIONEQUATIONwhere Equation (30) follows by expanding f and recalling that for simultaneously invariant kernels 31) follows by using the reproducing property of k σ(x) , Equation (32) uses the decomposition of g • σ in terms of the reproducing element k xi , and Equation ( 33) uses the definition of f . Now, as before, let f n → f be a Cauchy sequence, where f n lies in the span of the k x elements. We havek x • σ = k(σ -1 (x), •), Equation (EQUATIONby continuity of the inner product. Hence, S G is self-adjoint for all f, g ∈ H k .As S G is a projection (Equation ( 28)) and is self-adjoint Im(S G ) is a closed subspace of H k . Now, denoting the evaluation functional L x : H k → R, L x (f ) = f (x), note that since Im(S G ) is closed, then if L x was bounded on H k , it remains bounded on Im(S G ). Hence Im(S G ) is itself a RKHS, which we'll denote as H ′ .We will now relate k H ′ , the kernel associated with H ′ , to k and G. Note that, for f ∈ H ′ ,EQUATIONSo S G (k x ) is the reproducing element for H ′ , giving usEQUATIONwhich recovers the formula for the kernel in Equation 5.Finally, note that ⟨S G (k y ), S G (k x )⟩ recovers the more general "double sum" formula for producing a totally invariant kernel out of a kernel which is not necessarily simultaneously invariant [18, 33] :EQUATIONIn this section we derive Theorem 1, the upper bound on regret for Bayesian optimisation algorithms that incorporate invariance. The main elements of this proof were introduced in [41] . Our derivation is similar to the one in section C.1 of [22] , but extends that work to more general groups and kernels.In Appendix A.2.1, we present a concise summary of already known bounds on the information gain of dot-product kernels with polynomial eigendecay on the hypersphere. In Appendix A.2.2, we will adapt this argument for information gain of totally invariant versions of these kernels. We then use the bounds on information gain to bound the regret in Appendix A.2.3.We begin with the following remark which we use throughout this section. Namely, with the assumptions of Theorem 1, a kernel k on S d-1 which is simultaneously invariant with respect to the whole of O(d) is a dot-product kernel. To see this, choose r ∈ [-1, 1] and pick x r , y r such that ⟨x r , y r ⟩ = r and define κ(r) = k(x r , y r ). Observe that κ is well defined, as for any other choice of (x ′ r , y ′ r ) there is a (non-unique) isometry M(x ′ r ,y ′ r )(xr,yr) such thatEQUATIONEQUATIONand therefore k(x r , y r ) = k(x ′ r , y ′ r ). Hence, throughout this section we can consider dot-product kernels only.Consider the hyperspherical domain S d-1 = {x ∈ R d s.t. ∥x∥ = 1}. The Mercer decomposition of a dot-product kernel on the hypersphere isEQUATIONwhere Y i,p is the i-th spherical harmonic polynomial of degree k and N (d, k) is the number of spherical harmonics of degree k on S d-1 , given byEQUATIONAsymptotically in terms of k we haveEQUATIONIn this representation, µ k can be interpreted as an eigenvalue that is repeated N (d, k) times, with corresponding eigenfunctions Y 1,k , . . . , Y N (d,k),k . We let ψ be the maximum absolute value of the eigenfunctions, such that|Y i,k (x)| ≤ ψ ∀x ∈ S d-1 .We project onto a subspace of eigenfunctions corresponding to the first D eigenvalues. Let K be the number of distinct values in the first D eigenvalues. It follows thatEQUATIONEQUATIONEQUATIONEQUATIONwhere Equation (46) follows by substituting Equation ( 44), and Equation (47) follows by bounding the sum with an integral, and c i , C i are constants that ensure the asymptotic equalities hold.Truncating the eigenvalues at D leaves a tail with massEQUATIONEQUATIONwhere Equation (50) follows by substituting Equation ( 44), Equation (51) follows if the distinct eigenvalues have polynomial decay rate, µ k = O(k -β * p );2 , Equation (52) follows by bounding the sum with an integral; Equation (53) follows because D = Θ(K d-1 ), and hence K = Θ(D1 d-1 ).Theorem 3 in [41] states thatEQUATIONwhere |k(x, x ′ )| ≤ k ∀x and τ is the observation noise variance. Substituting in Equation ( 53) givesEQUATIONWe have freedom in setting D. We choose it so that the first term dominates, by settingEQUATIONThis is satisfied byEQUATIONwhich gives:EQUATION)where we have used the fact that ⌈x⌉ ≤ x + 1. Hence,EQUATIONA.2.2 Information gain for symmetrised kernels on the hypersphere Our upper bound involves comparing the corresponding reduction in maximal information across the invariant vs the standard kernel. In this section we will use overbars to refer to properties of the symmetrised (invariant) kernel. We remind the reader that G ≤ O(d) is a finite subgroup of isometries of the sphere S d-1 . Our aim will be to provide bounds on the asymptotics of the eigenspace dimensions N (d,k) N (d,k) where N (d, k) is defined as in A.2.1 and N (d, k) is the corresponding eigenspace dimension for the symmetrized kernel k G .We will follow the exposition from [4] , where asymptotic bounds on γ(d, k) := N (d,k) N (d,k) are given with two minor differences: 1. our treatment will handle any finite group of isometries and 2. we pay the price for generalizing to arbitrary groups by obtaining less strict decay rates, which will not trouble us for our analysis.We will show the following: Lemma 1. With the notations described in this section, assume G ≤ O(d) for d ≥ 10 is a finite group. Then, if -I ∈ G:EQUATIONEQUATIONand, otherwise:EQUATION)where C is a constant that depends on k, d and the spectra of elements in G.Proof. Recall from [4] , the formula:EQUATIONwhere µ is the uniform measure on the sphere S d-1 , and P d,k is the corresponding Gegenbauer polynomial. To demystify the conditional statement in the Lemma, notice first that P d,k (⟨x, x⟩) = P d,k (1) = 1, and P d,k (⟨-x, x⟩) = P d,k (-1) = 1 if k is even and -1 if k is odd. For all other σ ∈ G, we'll follow the argument in [4] to show that γ(d, k, σ) → 0 with asymptotics at least O 1 k . We don't make significant contributions to the method in [4] , here, other than to remark on the fact that many of the arguments they present follow through for finite subgroups of O(d).First remark that for σ ∈ G ≤ O(d) there is a matrix representation A σ , which admits a canonical basis in which the matrix elements are:EQUATIONwhere each R λ is a 2 × 2 matrix of the form a b -b a with eigenvalues λ and λ.Moreover, notice that since A |G| σ = I, all eigenvalues are roots of unity λ = e 2πip q with q |G|. The the eigenvalues of such matrices are the same as the allowable eigenvalues of permutation matrices considered in [4] , from which Proposition 6 gives:EQUATIONwhere s = max λ∈Spec(Aσ) {m λ + 4 • 1(|λ| < 1)}. Now assume σ ̸ = ±I. We havem 1 ≤ d -2, m -1 ≤ d -2 and m λ ≤ d 2 if λ ̸ = ±1. Collectively this implies -d + s ≤ -1 if d ≥ 10.In the following assume that -I / ∈ G. The corresponding argument for -I ∈ G is similar. In view of Lemma 1EQUATIONWe take the same number of distinct eigenvalues for the invariant and noninvariant kernels, K.However, there will be different multiplicities for the eigenvalues, leading to a different number of eigenvalues D. The corresponding tail has massEQUATIONwhere Equation ( 69 /a n = λ > 0, then ∀ϵ > 0, ∃N such that ∀N ′ > N we have n<N ′ b n ≥ (1 -ϵ)λ n<N ′ a n .Proof. From the definition of the limit, noting that the 1 -ϵ gives slack to bound λ n<N a n above by ϵ N <n<N ′ b n .Using this proposition on the sequences N (d, k) and N (d, k) we have further:EQUATIONAs previously we'll choose D T as a function of T such that the logarithmic term dominates (notice that since the first term is increasing in D any choice of D T greater than this will also satisfy), i.e.D T log 1 + kT |G| τ (1 + ϵ) ≥ T τ ψ 2 C 4 D -β * p +d-1 d-1 T (78) or D T ≥   T τ ψ 2 C 4 log 1 + kT |G| τ (1+ϵ)   d-1 β * p (79)With this choice:γT ≤ 1 + ϵ |G| T τ ψ 2 C 4 d-1 β * p log 1 + kT |G| τ (1 + ϵ) β * p -d+1 β * p (80)Hence,EQUATIONUnder the assumption of sub-Gaussian noise, it is shown in [42, Remark 2] that the maximum variance reduction algorithm presented in Algorithm 1 incurs simple regret satisfyingEQUATION)with probability at least 1 -δ, where δ ∈ (0, 1).Note that we are able to use this result without modification, as our approach only incorporates G-invariance through the kernel; hence, any changes to the regret due to the G-invariance will only appear in the γ T term. Substituting γ T from Equation (80) into Equation (82) givesEQUATIONwhere we have absorbed factors that do not depend on T or G, and let α =β * p -d+1 2β * p .Our goal is for r T = O(ϵ). This can be achieved by settingEQUATIONwith appropriate implied constants. Substitution of T from Equation (85) into Equation (84) makes it evident that for this T we haveEQUATION)Substituting α into Equation (85) gives the final upper bound on sample complexity, which isEQUATION)For the Matérn kernel on the hypersphere, [9] showed thatEQUATIONand hence, the MVR algorithm with an invariant Matern GP defined on S d-1 will achieve simple regret ϵ with probability 1 -δ ifEQUATIONwhere Õ hides logarithmic factors in ϵ and |G|. In this section, we derive Theorem 2, the lower bound on regret for the invariant Matérn kernel under permutation groups. The main elements of this proof were introduced in [12] , which extended the method in [35] . This constructive method serves to demonstrate the tightness of the upper bound.As in [12] and [42] , our method involves building up a set of candidate optimization targets in H k G . However, unlike these papers, it is no longer clear in this setting that such constructions are optimal in terms of the resulting dependence of the lower bound on |G| (the current gap to the upper bound potentially measures this failure rate). Our methods will produce elements in S G (F) ⊆ H k G , for a particular set of functions F, and it's the object of further research if these constructions can be improved.In Figure 7a , we have an example of the types of functions we construct, by applying the operator S G to a simple bump function. In this example, G is a group of rotations along the vertical axis. In Figure 7b is an example of a function we cannot retrieve through our method. This function is similarly invariant (also to all other rotations along the axis). Our method relies on a packing argument for the supports of such functions and it is a priori unclear how the tightest packing can be obtained.We identify a bounded subset of the RKHS of invariant functions that contains functions with the property that each point in X is ϵ-optimal for at most one function. We then construct a set of noisy bandit problems using functions from this class, before applying [12, Lemma 1] to lower-bound the number of samples that a learner requires in order to distinguish between them, and hence the best-case regret for this configuration.The proof for our lower bound will be constructive. We will exhibit a family of functions for which there is a computable lower bound of samples needed to distinguish their maxima. The challenge comes from constructing these functions natively in the RKHS of invariant functions. The strategy we employ here is to construct functions by symmetrization of simple bump functions. For the benefit of the reader, we'll show an explicit version of this construction first, on the hypercube X = [0, 1] d equipped with the action of G ≤ S d and subsequently extend to general actions of subgroups of O(d) on the sphere S d-1 .Lemma 2 (Bounded support function construction [12, Lemma 4] ). Let h(x) = exp -11-∥x∥ 2 1{∥x∥ 2 < 1} be the d-dimensional bump function, and define g(x, w, ϵ) = 2ϵ h(0) h( x w ) for some w > 0 and ϵ > 0. Then, g satisfies the following properties:1. g(x, w, ϵ) = 0 for all x outside the ℓ 2 -ball of radius w centred at the origin; 2. g(x, w, ϵ) ∈ [0, 2ϵ] for all x, and g(0) = 2ϵ;3. ∥g∥ k ≤ c 1 2ϵ h(0) 1 w ν ∥h∥ k when k is the Matérn-ν kernel on R d , where c 1 is constant. In particular, we have ∥g∥ k ≤ B when w = 2ϵc1∥h∥ k h(0)B 1 ν .Let e 1 , . . . , e d be the canonical basis vectors in R d . Then, the setEQUATIONforms a lattice of points on R d with spacing α. Equivalently, L(α) defines a partitioning of R d into disjoint regions or 'lattice cells' {R(x 0 ) : x 0 ∈ L(α)}, where R(x 0 ) is a d-dimensional hypercube of side length α centred on x 0 .Define the function class F(w, ϵ) = {g(x -x 0 , w, ϵ) : x 0 ∈ L(w)}. By property (1) of Lemma 2, g(x -x 0 , w, ϵ) is zero outside of R(x 0 ). Hence, the support of functions f ∈ F(w, ϵ) are pairwise disjoint. In addition, as F(w, ϵ) consists only of translates of g(x, w, ϵ), we retain max f = 2ϵ and ∥f ∥ k ≤ B ∀f ∈ F(w, ϵ) for appropriate w from properties 2 and 3 of Lemma 2 respectively.The hypercube [0, 1] d contains N of these functions, whereEQUATIONNow, we construct a set of G-invariant bump functions with disjoint support on the hypercube [0, 1] d , centred on the lattice cells for a class of groups G equipped with an action which we define below.Lemma 3 (Set of invariant functions). Let G be a subset of S d , the permutation group on d elements.G acts on R d (and on any invariant subset X ) by permutation of the coordinate axes. Let F(w, ϵ) = {g(x -x 0 , w, ϵ) :x 0 ∈ L(w)}. Then, F(w, ϵ, G) = {S G (f ) : f ∈ F(|G| 1/ν w, |G|ϵ)},where S G is defined as in Proposition 1, satisfies the following properties:1. F(w, ϵ, G) is a set of G-invariant functions on R d with pairwise disjoint support; 2. max x fi (x) = 2|G|ϵ |G•xi| , where fi (x) = S G (g(x -x i , |G| 1/ν w, |G|ϵ)) ∈ F(w, ϵ, G); 3. ∥ f ∥ ≤ B ∀ f ∈ F(w, ϵ, G).Proof. First, we show that F(w, ϵ, G) is indeed a S G -invariant set of G-invariant functions with disjoint support. As all translates g(x -x 0 ) are symmetric around x 0 , it holds that g(σ(x) -x 0 ) = g(x -σ -1 (x 0 )). From Equation (90), it is evident that the lattice L(α) is invariant to permutations; that is, σ(x 0 ) ∈ L(w) for all x 0 ∈ L(w), σ ∈ G. Hence, F(w, ϵ) itself is invariant under S G , and a fortiori F(w, ϵ, G) = S G F(|G| 1/ν w, |G|ϵ) .Next, we show that functions in F(w, ϵ, G) have maximum value 2|G|ϵ |G•xi| . Let f i ∈ F(|G| 1/ν w, |G|ϵ) be the function centred on x i . From Lemma 2, property 2, f i has maximum value 2|G|ϵ. LetG xi = {σ : σ(x i ) = x i } denote the stabilizer group of G at x i .Since the lattice itself is invariant, the image G•x i of the centre will have sizeEQUATIONThe preimage of any element in G•x i consists of exactly |G xi | elements. Therefore,EQUATIONwhere σ H is a representative of the coset H and f i • σ H is the unique value of f i • σ on this coset, regardless of its representative.Note that for any cosets H and K, f i •σ H and f i •σ K have disjoint supports from the same arguments as before, and thus,EQUATIONFinally, we show that functions in F(w, ϵ, G) have norm bounded by B. Let f ∈ F(w, ϵ, G) andf ∈ F(|G| 1/ν w, |G|ϵ). From property (3) in Lemma 2, ∥f ∥ k ≤ c 1 2|G|ϵ h(0) 1 |G| 1/ν w ν ∥h∥ k = ∥g∥ k ≤ B,where w is the same as before. From Proposition 1, recall that S G is a projection; hence, ∥f∥ k ≤ B implies ∥S G (f )∥ k G ≤ B, and so ∥ f ∥ k G ≤ ∥f ∥ k ≤ B.We now restrict the function set to the unit hypercube. The number of functions from F(|G| 1/ν w, |G|ϵ) (the set of rescaled non-invariant functions) that fit in the unit hypercube is given by N ′ , whereEQUATION)where N is defined in Equation (91).Clearly, the number of orbits in F(|G| 1/ν w, |G|ϵ), i.e. the number of elements in F(w, ϵ, G) is equal to the number of orbits of G on L(w), by Lemma 3. We would like to count the number of orbits, or at least bound it from below. For specific groups G one can do this explicitly, however, in general it is impossible to give a count of the orbits, without further assumption. What we can, however say is that the number of orbits of size exactly |G| represents almost the total number of orbits when w is small.Proposition 3. Let M = | F(w, ϵ, G)| = | L(|G| 1/ν w) [0, 1] d /G|. Let K := ⌊ 1 |G| 1/ν w ⌋. Then M ≥ K d|G| , and this bound is asymptotically tight in the limit of w → 0. EQUATIONThis is easy to see asEQUATIONHence, we have that| F(w, ϵ, G)| approaches |F (|G| 1/ν w,|G|ϵ)|from above as the width w tends to 0:EQUATIONNote that M is also the number of functions from F(w, ϵ, G) that are supported in [0, 1] d . Then, we can rewrite the above inequality as:EQUATION)where N is the number of functions from F(w, ϵ) that fit in [0, 1] d , as defined in Equation (91). Equation (103) gives a direct relationship between the packing of invariant and non-invariant bump functions with bounded norm that fit in the unit hypercube.Note that the above analysis shows that this bound is tight for small enough w which depends only on d and ν since |G| < d!. If one is interested only in Equation ( 103), this is easy to obtain as M , the number of orbits is greater than N ′ divided by the size of the largest orbit.Note also that in the case of noiseless observations any optimisation algorithm will have to query at least M (resp. N ) points to be able to distinguish between functions in F(w, ϵ, G) (resp. F(w, ϵ)) defined on [0, 1] d ; hence, Equation (103) (resp. Equation ( 91)) can also be interpreted as a lower bound on the sample complexity of finding the maximum of these functions.The subsequent section of the proof mirrors the approach used to derive a simple regret lower bound in the standard (non-invariant) setting, as detailed in [12, Section 4.2] . The primary distinction lies in the construction of the representative functions that exhibit group invariance, as outlined in Equation ( 104) and Equation (105).Recall the definition of a region R(•) from the discussion following Equation (90). Define Ri as Ri = σ∈G R(σ(x i )), where x i ∈ L(|G| 1/ν w) is the center of the region. Note that these regions are self-contained under G, so that x ∈ Ri =⇒ σ(x) ∈ Ri , and are pairwise disjoint. Therefore, we define a set of M such regions { Ri } M i=1 , which partition our domain. The value of M is lower bounded in Equation (103).We replace B with B/3 in the construction of F(|G| 1/ν w, |G|ϵ) (see Lemma 3). 3 Let f i and f j (i ̸ = j) be shifted bump functions from the set F(|G| 1/ν w, |G|ϵ), centered at some fixed x i and x j respectively, where x i , x j ∈ L(|G| 1/ν w). Furthermore, these are selected such that the sizes of the orbits, |G • x i | and |G • x j |, are both equal to |G|. Then, we construct the pairEQUATIONEQUATION)The following holds ∥ f ∥ k G ≤ B/3 and ∥ f ′ ∥ k G ≤ B (from Lemma 3 and triangle inequality). We observe that f and f ′ coincide outside Rj i.e., f ′ -f is only non-zero in region Rj . Moreover, we have max x∈ Ri f (x) = 2ϵ and max x∈ Rj f ′ (x) = 4ϵ. It follows that only points from Rj and Ri can be ϵ-optimal for f ′ and f , respectively.In the subsequent steps, let P f [•] and E f [•] denote probabilities and expectations with respect to the random noise, given a suitably chosen underlying function f . Additionally, we define P f (y|x) as the conditional distribution N (f (x), σ 2 ), consistent with the Gaussian noise model.We fix T ∈ Z + , δ ∈ (0, 1/3), B > 0 and ϵ ∈ (0, 1/2), and assume the existence of an algorithm such that for any f ∈ F k G (B), it achieves an average simple regret r(x (T ) ) ≤ ϵ with a probability of at least 1 -δ. Here, x (T ) represents the final point reported by the algorithm. Note that an algorithm that achieves simple regret less than ϵ for all f ∈ F k G (B) automatically solves the decision problem of distinguishing between f and f ′ as above. If no such algorithm exists, the sample complexity necessary to achieve this bound on regret must be higher than T . Next, we establish a lower bound on the number of samples required for such an algorithm to differentiate between two bandit environments (with reward functions f and f ′ and the same Gaussian noise model) that share the same input domain but have distinct optimal regions.Let A denote the event that the returned point x (T ) falls within the region Ri . Assume that an algorithm achieves a simple regret of at most ϵ for both functions f and f ′ (both f , f ′ ∈ F k G (B) as established before) each with a probability of at least 1 -δ. Under our assumption on simple regret, this implies P f [A] ≥ 1 -δ and P f ′ [A] ≤ δ since only points within Ri can be ϵ-optimal under f , and only points within Rj can be ϵ-optimal under f ′ .We are now in position to apply [12, Lemma 1] . We apply this lemma using the previously defined f , f ′ , and A. 4 It therefore holds that, for δ ∈ (0, 1 3 ):EQUATIONwhere C m is the number of queried points in Rm up to time T .For the two bandit instances and any input x, the observation distributions are N ( f (x), σ 2 ) and N ( f ′ (x), σ 2 ). Utilising the standard result for the KL divergence between two Gaussian distributions with the same variance, we obtain:EQUATIONas f (x) = f ′ (x) for x / ∈ Rj , and for x ∈ Rj we have that f (x) = 0 and max x∈ Rj f ′ (x) = 4ϵ.Hence, Equation (106) becomesEQUATION)Let J = {j ∈ I|x j ∈ W(0)} be the set of indices j such that f j is centred around a point with no repeated indices such as in Proposition 3. Let ρ w = |J| |G|M be the fraction of all orbits of size exactly |G|. By Proposition 3, for any small δ w , there is a correspondingly small enough w such that ρ w ≥ 1 -δ w . Since the previous inequality holds for an arbitrary j, we sum over over j ∈ (J \ G{i}) /G to arrive at:EQUATION)By rearranging and noting thatEQUATIONSubstituting in M from Equation (103) givesEQUATIONand using N from Equation (91) gives the final boundEQUATIONwhere the implied constants can depend on d, l, and ν.A. 4 Lower bounds on different underlying spaces.In this section we'll investigate how to raise the lower bounds on sample complexity to embedded submanifolds X → R d . The full treatment below is valid only for the hypersphere, X = S d but many of the ingredients required are valid on other manifolds as well. When specific statements about the sphere are made in the treatment below, we will explicitly refer to S d instead of X .In the proof in A.3, we used the existence of a G-invariant lattice to underpin our bump function construction. Moreover, we've used the fact that the RKHS norm B := ∥g∥ k ∼ ϵ w ν . We will replicate the underlying function class construction for F(w, ϵ, G) here.We will recover the same property for bump functions on X but under a different notion of rescaling. Pick a point p ∈ U ⊂ X and consider (y 1 , ..., y m ) a system of coordinates on U so that x i (y 1 , ..., y m ) represent the standard coordinates on the ambient space R d . Let f p,w,ϵ : R d → R + be a bump function around p which is 0 outside B w (p; R d ) and supported in B w/2 (p; R d ) with height ϵ, as in section A.3. Let ϕ p,w,ϵ = f p,w,ϵ X .If X = S d , then the intersection B rw (p; X ) := B w (p; R d ) X is a geodesic neighbourhood with geodesic radius r w = arccos 2-w 2. For small enough w,EQUATION)by computing the Taylor expansion.We will lso need to compute the RKHS norm for ϕ p,w,ϵ for some suitably defined kernel. Let k ν : R d × R d be the standard Matérn kernel, and consider kν its restriction to X .Next, we recall the existence of a fundamental domain for the action of finite groups G on X . Following the definition of fundamental domain from [21] , we say F ⊂ X is a fundamental domain for the action of G if:• F is a domain. • G • F = X • gF F ̸ = ∅ =⇒ g = e,the identity.• For any compact K ⊂ X , the transporter set ( F |K) G is finite.Note that the last condition is trivial if G is finite. The construction of F is simple given that G acts by isometries of the underlying space R d and g X = ι * g R d .Indeed, let x be a point such that |Gx| = |G| (such a point trivially exists, as each member of G other than the identity fixes a subspace of dimension at most d -1) on the X = S d-1 . Then, |Gx| = |G|, and the Dirichlet domain:EQUATIONis a fundamental domain for the action (see [21] for proof). One can see this set of points as the interior of the points in X which are closest to x rather than any of its other translates under G. For example, for the action of S d on the positive orthant S d-1 , it is any of the simplices in the barycentric subdivision.Equipped with this definition, we construct the lattice of bump function centres on X as follows. First construct a 2r w -net L(r w ; F x ) on F x as the set of centres of an optimal geodesic ball packing of F x . Then defineEQUATIONWe will need to compare the cardinalities of L(r w ) and L(tr w ) for some scaling factor t > 0. Let Π(r w ) be the packing number of F ⊆ S d-1 by geodesic balls of radius r w .We have trivially thatEQUATIONwhere the first inequality comes trivially from noticing that a maximal r w packing is a 2r w covering. Moreover, since the exponential map is a diffeomorphism on small enough neighbourhoods (with derivative I at the origin), we have that:EQUATIONfor some 0 < K 1 < K 2 . Combining equations 117 and 118 gives:EQUATIONFollowing this discussion, we have the following lemma which is an extension to Lemma 3. Denote by w r the inverse mapping of w → r w , defined on small enough geodesic balls on S d-1 . We write F(w, ϵ) = {f p,w,ϵ : p ∈ L(r w )} for a set of functions on the ambient space and F X (r w , ϵ) = ϕ p,w,ϵ := f p,w,ϵ X : p ∈ L(r w ) for the restriction F(w, ϵ) X . Similarly, we write F(w, ϵ, G) = {S G (f ) : f ∈ F(|G| 1/ν w, |G|ϵ)} and FX (w, ϵ, G) = F(w, ϵ, G) X Lemma 4 (Invariant functions on the sphere). The sets FX (r w , ϵ, G) and F(w, ϵ, G) satisfy the following properties.1. F(w, ϵ, G) is a set of G-invariant functions on S d-1 with pairwise disjoint support; 2. max x S G (f p,w,ϵ )(x) = 2|G|ϵ |G•p| . 3. ∥ f ∥ ≤ B ∀ f ∈ FX (w, ϵ, G).Proof. The key ingredient for 1. is to note that each f p,w,ϵ (x) = g( x-p w ), and thus, each f is G-invariant around the point p, namely, for σ ∈ G, f p,w,ϵ • σ = f σ -1 (p),w,ϵ . Then, since L(r w ) is G-invariant by construction, 1. follows. For 2. the proof is identical to Lemma 3 2. For 3. notice that ∥ϕ p,tw,ϵ ∥ kν ≤ inf ϕ e ∥ϕ e p,tw,ϵ ∥ kν ≤ ∥f p,tw,ϵ ∥ kν = t ν ∥f p,w,ϵ ∥ kν (120)where the second inequality follows from the restriction and ϕ e is any extension of ϕ to the underlying space, since the norm of a restriction is less than that of all extensions.Let N = |F(w, ϵ)|, M = F(w, ϵ, G). Recall that these functions are defined in terms of the packing numbers of F , so that N = |G|Π(r w ) and M = Π(r |G| 1/ν w ). Pick δ w > 0. For all small enough w, we have that:EQUATIONwhere the first inequality comes from equation 114 and the fact that the packing number is a decreasing function of the radius, the second inequality comes from equation 119 and c 3 is a constant which depends on the ambient dimension and the curvature of the space X .Following the arguments of section A.3.2 we see, that by construction, L(r w ) admits only orbits of size |G|, end hence in Equation 110, the sum is over all elements in FX (r w , ϵ, G). Hence we recover the same lower bound on sample complexity:EQUATIONThis concludes the proof of Theorem 2.The key element is the original construction of the functions f p,w,ϵ (x), which are of the form g(x -p) for some bump function constructed around the origin which is itself symmetric with respect to the action of O(d) (and hence its subgroups).We use an isotropic Matérn-5 /2 kernel as the base kernel k, and compute k G according to Equation (5).To generate the objective functions, we first generate a large number n of sample pairs (x i , f (x i )) from a zero-mean GP prior with kernel lengthscale l = 0.12. We scale n with the dimensionality of the problem d, so that n = 64 for d = 2, n = 256 for d = 3, and n = 512 for d = 6, to mitigate sparsity. We then use BoTorch to fit the noise parameter of a noisy invariant GP with fixed lengthscale l = 0.12 to that data, and treat the mean of the trained GP as our objective function. In this way, we ensure that the objective function belongs to the desired RKHS.During the BO loop, the learner was provided the true lengthscale, magnitude, and noise variance of the trained kernel, to avoid the impact of online hyperparameter learning. For UCB, we use β = 2.0, except in the quasi-invariant case with ε = 0.05 where we found that increasing β was necessary to achieve good performance (β = 3.0). MVR has no algorithm hyperparameters.For the constrained BO baseline, we replaced BoTorch's default initialisation heuristic with rejection sampling in order to generate multiple random restarts that lie within the fundamental domain. We found that this improved the diversity of the initial conditions, and hence boosted the performance of the acquisition function optimiser.In all other cases BoTorch defaults were used.We (c) If the contribution is a new model (e.g., a large language model), then there should either be a way to access this model for reproducing the results or a way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset). (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility.In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.github.com/theo-brown/invariantkernels and github.com/theo-brown/bayesopt_with_invariancesMatérn kernels on S d-1 satisfy this property with β * p = 2ν + d -1[9].We note that this modification affects only the constant factors (also in the bound on M from Equation (103)), which are not important for our analysis.Following the approach in[12], we deterministically set the stopping time in Lemma 1 to T , corresponding to the fixed-length setting in which the time horizon is pre-specified. For an adaptation to scenarios where the algorithm may determine its stopping time, refer to Remark 1 in[12]. 


Decision: The reviewers discuss the incorporation of invariance under a group action into a standard UCB-based BO loop, and show that this leads to improvements both empirically and theoretically.  Coming out from the discussion period the reviewers were unanimously positive about the paper.  This seems like a clear case for acceptance.

Reviews: review: paper_summary: paper_summary: The paper targets the Bayesian optimization problem for a class of invariant functions, which is useful in many fields including machine learning and physics. Specifically, the paper proposes to incorporate the invariances into the kernel of the GP to produce invariance-aware algorithms, either fully or partially. The paper presents theoretical analysis regarding the lower bound on sample complexity when using invariant kernels in BO. Several experiments are shown to support the findings.
strength_weakness: strength_weakness: •	In Figure 3, the plots should be consistent, i.e., either use cumulative regret or simple regret. 
•	The experiments are only conducted on low-dim problems.
questions: questions: I don’t have many questions regarding the submission. I think my main question is on the performance of the algorithms for higher-dimensional problems. Will the methods work well in high-dim regime? And if not, are there anything we can do to make it work better in the high-dim setting?
limitations: limitations: I agree with the author’s limitations regarding the computational expense of fully invariant kernel on large group of transformation.

No negative societal impact needs to be addressed.
paper_summary: paper_summary: The paper introduces an Bayesian Optimisation method which is able to take into account known invariances of the objective function.  Specifically, it is assumed that the objective function remains invariant under a finite group action $G$.     

The approach is straightforward, one does standard Bayesian optimisation using a $G$-invariant kernel, following the approach of [Haasdonk et al, 2007] and others.    The authors then propose using standard BO with UCB or Maximum variance acquisition function. 
  The paper provides theoretical bounds on the regret for this invariance aware approach -- in the special case where the domain is a sphere.    The authors demonstrate the sample efficiency of the method empirically, on some synthetic experiments, and then a nice tokamak optimisation example.
strength_weakness: strength_weakness: Dealing with invariances in optimisation (Bayesian / black-box / otherwise) is certainly not a new problem, and there are typically well-established methods to deal with this, through symmetry breaking constraints etc.   This is particularly common in the context of computational chemistry and materials design properties where molecules have a wide range of symmetries.   One noteworthy example is Bayesian Optimization With Symmetry Relaxation Algorithm from [Zuo, Yunxing, et al. "Accelerating materials discovery with Bayesian optimization and graph deep learning." Materials Today 51 (2021): 126-135].     So while the approach is novel, it is far from unique, and I feel that the authors have not really engaged with existing methods at all.   The literature review seems to go from "Invariances in DL" -> "Invariant Kernels" -> "Kernels on Manifolds" (again missing a lot of important works there too) -> "regret bounds for BO". 

The sample complexity bounds are certainly novel, but for example the upper bound largely seems to follow Joan Bruna's existing work on this which tackles a very similar problem.    

So to summarise:
(+) I think this is nice work.
(+) I think its a novel approach.
(+) The theory is nice.
(-) However, it feels like a quite incremental contribution, which makes no attempt to engage with alternative approaches to the same problem.
(-) The theory, while nice, again seems a moderate modification of existing work.
questions: questions: 1. Could the authors help distinguish their work from other approaches in this field.   I strongly urge the authors to engage a bit more with other proposed methodology (some of which might not be very rigorous).  Understandably this is a fairly common issue, and many resolutions exist -- certainly not all as nice as this paper.

2. In particular, it would be good if the authors could identify strong advantages of this method over others, and maybe even provide other methods as a useful benchmark.

3. One possible advantage of this approach over others is situations where the invariance is not strongly constrained, e.g. quasi-invariances.  Could the authors demonstrate that this method still yields better sample complexity when the objective is nearly invariant?

I am v happy to adjust my score of the authors can provide a stronger case for the novelty and/or value-add of this work over other approaches.
limitations: limitations: The limitations addressed are around the complexity of computing the invariant kernel -- this has been addressed and some guidelines for mitigating this are discussed.
paper_summary: paper_summary: The authors proposed a new setting where the invariance is either known or partially known. They theoretically examined the upper and lower bounds of convergence rates concerning sample complexity. Their findings demonstrated both theoretical and empirical superior performance over the standard UCB approach in cases with known invariance. Additionally, they explored more practical settings with partially known invariance using synthetic tasks and a real-world example from tokamak optimization.
strength_weakness: strength_weakness: **Limited Applicability:**
I am not convinced that this new setting is widely applicable to real-world scenarios. Aside from the specific tokamak application, it is challenging to identify other examples where this setting would be useful. While regression tasks in previous work have broad applicability, is the same true for optimization? Can you provide more intuitive examples of invariance-aware BO? For instance, in classification tasks, such as translation invariance in cat images in the introduction section, the concept is clear. However, does this example apply to BO? Beyond the tokamak example, the authors mention material science, but it is difficult to envision practical situations in this domain where the search space has known invariance. For example, while crystal lattice invariance is known, is this relevant to the search space? This is merely one feature characterizing materials, and it is hard to imagine a scenario where we aim to maximize something within this symmetric lattice space. While band gap optimization can be the one (finding minimum and maximum of electronic energy), it is computationally simple and not computationally demanding.

**Feasibility of Listing Possible Invariances:**
It is hard to imagine users being able to list possible invariances. This raises concerns about the motivation behind the work.

**Lack of Simpler Baselines:**
The work lacks comparisons with simpler baselines, such as periodic kernels or constrained BO. For example, in Figure 2, if the search space is known to be symmetric for 10 cycles, why not constrain the search space to 1/10 as in typical constrained BO approaches? Finding all peaks in Figure 2 seems wasteful since we know they are repetitive. The same applies to permutation invariance; why not algorithmically reject repetitive candidates?

**Clarity in Experimental Section:**
The experimental section is difficult to follow. The tasks in Figure 3 are unclear, making it hard for readers to reproduce the results.

**Minor Points**
- Numerous typos (e.g., lines 102-103, 116-117, repeating the RKHS acronym definition).
- Limited baseline comparisons (e.g., permutation kernel, periodic kernel, and constrained BO).
questions: questions: The questions are detailed in the weakness section above.
limitations: limitations: Limitations are well discussed.
paper_summary: paper_summary: This paper provides sample complexity bounds for Bayesian optimization under settings with invariant kernels. Invariant kernels are able to model functions which are invariant under transformation families. Such models allow us to carry out Bayesian optimization efficiently due to having to being able to obtain far much information with a single query than if we used standard kernels (as we obtain information over the whole group of invariances with a single observations, while standard BayesOpt will carry out redundant sampling). The theoretical effect of invariances is investigated in this paper.

A simultaneously invariant kernel satisfies the property that $k(\sigma(x), \sigma(y)) = k(x, y) \forall \sigma \in G$. Where $G$ is the group of invariance. The paper focuses on finite subgroups of isometries, which maintain distances and dot products and therefore any stationary kernel will be simultaneously invariant for any group in question. The paper shows in Proposition 1 how a simultaneously invariant kernel can be used to build an RKHS with the invariances translated to the function space as well. We can then use this space of functions with their corresponding kernel to model the black-box functions for Bayesian Optimization. 

The analysis is carried out for two standard BayesOpt algorithms: (a) maximum variance reduction, which selects as the next query the point with maximal variance, and (b) the standard upper confidence bound acquisition function. The first main result is shown in Theorem 1, where the authors obtain an upper bound on the maximal information gain of any bandit algorithm, and of particular importance is a factor of $O(1 / |G|)$. They then show this translates to a $O(1 / |G|)^{(2v + d - 1) / 2v}$ dependency on the sample complexity of the maximum variance reduction algorithm under a Matern kernel. A lower bound is then shown, with complexity $O(1 / |G|)^{(2v + d - 1) / v}$. This shows the upper bound is relatively tight, and that algorithms are likely to fail when run for anything lower than the sample complexity.

Two main experiments are ran:

1. The first experiment looks at the optimization of Matern function samples under different invariance groups. This results in significant speed-ups for optimization in all three benchmarks. In the PermInv6D benchmark, further analysis is carried out where only certain subgroups of the invariances are included in the known invariance group; these result in the expected outcome and the regret improves as more invariances are added.

2. A real-world example is then introduced, where a simulated nuclear reactor configuration is optimized. A subgroup of the invariances is used for computational efficiency, and the results show a very significant improvement from using invariant kernels.
strength_weakness: strength_weakness: - The modeling comes from previous literature, and the bandit strategies followed are standard. The paper's contributions are strictly the theoretical additions, and the convincing, yet somewhat short, empirical study.

- From my understanding, sample complexity upper bounds are only given for the maximum variance reduction algorithm, which could be argued is less important for many optimization applications than UCB.

- Perhaps more applications of the algorithm could be showcased or discussed.
questions: questions: - How do you do the batching of the acquisition function in the nuclear fusion experiment? 

- Is there anything that can be said about the sample complexity of the UCB acquisition function? From L195 I am under the impression that the upper bound only holds for MVR.

- Is there more intuition into why the bounds for maximal information gain are independent of the bandit strategy chosen?
limitations: limitations: Limitations are addressed in Section 5.

Comments: title: Rebuttal summary and acknowledgements
Response
Request for response
comment: Thank you for the response. The response addresses my concerns so I keep my current score.
To all reviewers and area chairs,

As the rebuttal period it nearing its end, we would like to thank all the reviewers again for taking the time to read our manuscript, to provide constructive criticism and raise many important points that have helped improve both our understanding and the write-up. 

We would like to doubly thank reviewers for engaging in the rebuttal discussion and for providing _unanimous approval for our paper_. We will improve the manuscript according to the discussion below to ensure that it passes the high quality standards expected from a NeurIPS publication.
Thank you for the author for their effort and clear rebuttal. My concerns are adequately addressed, so I raise the score.
I thank the authors for their rebuttal and for their suggested updates.   I feel that I understand a bit better the authors contributions on the theoretical aspects of this work, which are the main contribution of this paper, in particular how they distinguish with other works in this area.    I think this is broadly a well written paper, but I do agree with the other reviewers that perhaps some other applications should have been showcased to demonstrate the feasibility of identifying (partial or otherwise) invariance.   I will update my score.
We thank all the reviewers again for their time taken to review our work.

**As the author-reviewer discussion period is nearing its end, we would appreciate it if the reviewers would respond to our rebuttals.**

If there any further questions, we are more than happy to clarify and discuss further.
Thanks for the response, the clarifications, and for the additional experiments. 

> ...but believe they make sufficient advancement into the field of BO and learning with symmetries, that they present a benchmark of interest to the community

I do agree with this, which is why I recommended acceptance. However, as the contributions are mainly theoretical, and I am not well acquainted with the theoretical side of BO, it is difficult for me to measure the strengths of the paper (which is the reason for my low confidence). 

Regarding the difficulty to come up with bounds for the UCB algorithm, it remains a weakness of the paper even if a difficult task.

Based on these, I will stand by my original review; I liked the paper and I recommend acceptance, however, with low confidence.