High-dimensional Bayesian optimization using low-dimensional feature spaces
Riccardo Moriconi Marc P. Deisenroth K. S. Sesh Kumar





High-dimensional Bayesian optimization using low-dimensional feature spaces        





Riccardo Moriconi Department of Computing 

              Imperial College London, UK

                         r.moriconi16@imperial.ac.uk 
            Department of Computer Science 

              University College London, UK 

                         m.deisenroth@ucl.ac.uk
              K. S. Sesh Kumar Data Science Institute 

              Imperial College London, UK

            s.karri@imperial.ac.uk






High-dimensional Bayesian optimization using low-dimensional feature spaces

    Riccardo Moriconi                
Marc Peter Deisenroth  
        K. S. Sesh Kumar 

    
========================================================================================




Bayesian optimization (BO) is a powerful approach for seeking the global optimum of expensive black-box functions and has proven successful for fine tuning hyper-parameters of machine learning models. However, BO is practically limited to optimizing 10–20 parameters. To scale BO to high dimensions, we usually make structural assumptions on the decomposition of the objective andor exploit the intrinsic lower dimensionality of the problem, e.g. by using linear projections. We could achieve a higher compression rate with nonlinear projections, but learning these nonlinear embeddings typically requires much data. This contradicts the BO objective of a relatively small evaluation budget. To address this challenge, we propose to learn a low-dimensional feature space jointly with (a) the response surface and (b) a reconstruction mapping. Our approach allows for optimization of BO's acquisition function in the lower-dimensional subspace, which significantly simplifies the optimization problem. We reconstruct the original parameter space from the lower-dimensional subspace for evaluating the black-box function. For meaningful exploration, we solve a constrained optimization problem[Implementation available at https://github.com/rm4216/BayesOpt].




§ INTRODUCTION

Bayesian optimization (BO) is a useful model-based approach to global optimization of black-box functions, which are expensive to evaluate <cit.>. This sample-efficient technique for optimization has been effective in experimental design of machine learning algorithms <cit.>, robotics applications <cit.> and medical therapies <cit.> for optimization of spinal-cord electro-stimulation. Despite its great success, BO is practically limited to optimizing 10–20 parameters. A large body of literature has been devoted to address scalability issues to elevate BO to high-dimensional optimization problems, such as discovery of chemical compounds <cit.> or automatic software configuration <cit.>.

The standard BO routine consists of two key steps: (i) estimating the black-box function from data through a probabilistic surrogate model, usually a Gaussian process (GP), referred to as the response surface; (ii) maximizing an acquisition function that trades off exploration and exploitation according to uncertainty and optimality of the response surface. 
As the dimensionality of the input space increases, these two steps become challenging. The sample complexity to ensure good coverage of inputs for learning the response surface is exponential in the number of dimensions <cit.>. With only a small evaluation budget, the learned response surface and the resulting acquisition function are characterized by vast flat regions interspersed with highly non-convex landscapes <cit.>. 
This renders the maximization of the acquisition in high dimensions inherently hard <cit.>.


High-dimensional optimization is often translated into low-dimensional problems, which are defined on subsets of variables <cit.>. These approaches apply a divide and conquer approach to decompose the problem into independent <cit.> and potentially dependent components <cit.>. 

However, high-dimensional data often possesses a lower intrinsic dimensionality, which can be exploited for optimization.
A feature mapping can then be used to map the original D-dimensional data onto a d≪ D-dimensional manifold. For example, in <cit.>, the authors used random linear mappings to reduce dimensionality of the optimization problem. Similar approaches, which use linear dimensionality reduction, drive exploration in BO to actively learn this linear embedding <cit.>. While these methods perform well in practice, they are restricted to linear subspaces of the original domain. With nonlinear embeddings, higher compression rates are possible. In our work, we focus on this nonlinear setting.



Using BO with nonlinear feature spaces was proposed in <cit.>. In <cit.>, a low-dimensional data representation is learned with variational autoencoders (VAEs) <cit.>. 

However, this approach requires both large amounts of data and learning the model offline without the possibility to update the learnt feature space during optimization. Nevertheless, in the specific application of automatic discovery of molecules, where libraries of existing compounds are readily available prior to optimization, this approach makes much sense. To accommodate fairly small evaluation budgets, in our work, we exploit a  probabilistic model based on GPs, which features superior data efficiency with respect to VAE-based approaches <cit.>.

VAE models <cit.> were used to propagate uncertainty of latent space representations through the response surface model with Gaussian process latent variable models <cit.>. However, in <cit.>, the latent space representation is not learned specifically for the regression task (learning the response surface). Gradient-based methods <cit.> have been used to learn a lower-dimensional Riemannian manifold for optimization and sampling.


Nonlinear embeddings also allow for modeling non-stationary objective functions. In this context, a hierarchical composition of GPs, referred to as deep GPs <cit.>, is especially useful when the response surface is characterized by abrupt changes or has constraints. An extensive investigation on the employment of deep GP models in BO is presented in <cit.>. In our work, we also exploit the idea of learning highly nonlinear functions through the composition of simpler functions <cit.>, but we focus on deterministic dimensionality reduction and optimization in feature space.











In this paper, we propose a BO algorithm for high-dimensional optimization, which learns a nonlinear feature mapping 𝐡:ℝ^D→ℝ^d to reduce the dimensionality of the inputs, and a reconstruction mapping 𝐠ℝ^d→ℝ^D based on GPs to evaluate the true objective function, jointly, see Figure <ref>. This allows us to optimize the acquisition function in a lower-dimensional feature space, so that the overall BO routine scales to high-dimensional problems that possess an intrinsic lower dimensionality. Finally, we use constrained maximization of the acquisition function in feature space to prevent meaningless reconstructions.







































§ BAYESIAN OPTIMIZATION









Bayesian optimization is a powerful tool for globally optimizing black-box functions that are expensive to evaluate <cit.>. In our setting, we consider the global minimization problem

    𝐱^* = 𝐱∈𝒳min  f_X(𝐱)

with input space 𝒳=[0,1]^D and objective function f_X𝒳→ℝ. We consider functions f_X that are costly to evaluate and for which we are allowed a small budget of evaluation queries to express our best guess of the optimum's location 𝐱^* in at most T_end iterations. We further assume we have access only to noisy evaluations of the objective y=f_X+ε, where ε∼𝒩(0, σ_n^2) is i.i.d. Gaussian measurement noise with variance σ_n^2. 

We restrict ourselves to the typical setting, where neither gradients nor convexity properties of f_X are available.












































 The main steps of a BO routine at iteration t involve (i) response surface learning, (ii) optimal input selection x⃗_t+1 and (iii) evaluation of the objective function f_X at x⃗_t+1. The first step trains a probabilistic surrogate model p(f_X), the response surface, which describes the black-box relationship between inputs x⃗ and observations y. In the (t+1)st iteration of BO, the optimal input selection step finds an input x⃗_t+1 that maximizes an acquisition function α(·), which describes the added value of input x⃗_t+1. The evaluation step returns a noisy observation of the true objective function f_X(x⃗_t+1)+ε at the selected location. These steps are summarized in lines <ref>, <ref> and <ref> of Algorithm <ref>, respectively. Having defined a probabilistic surrogate model for our objective function, which is usually modeled by a GP <cit.>, we can compute posterior predictions of objective function values at test locations. These posterior predictions are then fed to the acquisition function, which drives exploration during optimization. Posterior predictions of the GP are Gaussian distributed with mean μ and variance σ^2. Defining Z(x⃗):=(f_min - μ(x⃗))/σ(x⃗) and f_min:=x⃗∈X⃗_tmin f(x⃗), this  allows us to define three different acquisition functions to maximize:

 
    α(x⃗)    =Φ(Z(x⃗))  Probability of improvement (PI) <cit.>
    α(x⃗)    =σ(x⃗)Z(x⃗)Φ(Z(x⃗)) +σ(x⃗)ϕ(Z(x⃗))  Expected improvement (EI) <cit.>
    α(x⃗)    =-μ(x⃗) +β_tσ(x⃗) Upper confidence bound (UCB) <cit.>.


Here, ϕ and Φ denote the probability density function and the cumulative density function of the standard normal 𝒩(0,1), respectively. The parameter β_t controls the exploration exploitation trade-off. For a complete review on acquisition function the reader is referred to <cit.>. In high-dimensional settings (D>20), both the response surface learning and optimal input selection via optimization of the acquisition function are computationally challenging.



§ BAYESIAN OPTIMIZATION IN LOW-DIMENSIONAL FEATURE SPACES

In this section, we consider a setting, where the input space 𝒳 is high-dimensional and the objective function possesses an intrinsic lower dimensionality. In our work, we exploit the effective low dimensionality of the objective function for BO in a lower-dimensional feature space 𝒵⊂ℝ^d, where d≪ D. In particular, we express the true objective f_X:ℝ^D→ℝ as a composition of a feature mapping 𝐡:ℝ^D→ℝ^d and a  function f_Z𝒵→ℝ so that f_X = f_Z∘h⃗. The lower-dimensional feature space allows for both learning the response surface f_X and maximizing an acquisition function α with domain 𝒵, which yields optimizer z⃗_*. 





Since we cannot evaluate the true objective f_X directly at the low-dimensional features z⃗_*, we project z⃗_* back into the D-dimensional data space  𝒳 by means of a reconstruction mapping 𝐠𝒵→𝒳. We can think of this mapping as a decoder within an auto-encoder framework. We model both the composition f_X := f_Z∘𝐡 and the reconstruction with GPs <cit.>. Algorithm <ref> summarizes the main steps of this feature-space BO.

In the following, we detail the model (see Figure <ref>) for jointly learning the feature map 𝐡(·), the low-dimensional response surface in feature space f_Z, and the reconstruction mapping 𝐠(·).


















 §.§ Manifold Gaussian processes for response surface learning in feature space



We expect the response surface to predict the value of the black-box objective function f_X with calibrated uncertainty associated with each prediction. 
GPs are probabilistic models that allow for an analytic computation of posterior predictive function values within a Bayesian framework, and they are the standard model in BO for modeling the response surface. 

A GP is a distribution over functions f_Z∼𝒢𝒫(m(·), k(·, ·)) and is fully specified by a mean function m𝒵→ℝ, and a covariance functionkernel k𝒵×𝒵→ℝ. The kernel computes the covariance between pairs of function values as a function of the corresponding inputs, i.e. Cov(f_Z(𝐳), f_Z(𝐳')) = k(𝐳, 𝐳'), and thereby encodes regularity assumptions about f_Z, such as smoothness or periodicity. Common kernel choices in the BO literature include the squared exponential and Matérn  kernels <cit.>.




In our feature space optimization, we phrase lines <ref>–<ref> of Algorithm <ref> as a single learning problem. Therefore, we need a GP that learns useful representations z⃗ of inputs x⃗ for the regression task together with f_Z. A manifold GP (MGP) <cit.> addresses this issue by composing two mappings: The deterministic feature map 𝐡 with parameters θ_h and a GP f_Z∼𝒢𝒫(m, k) with kernel hyper-parameters θ_k. The GP  models the relationship between features z⃗ and function values y∈ℝ in observation space. The resulting composite model f_X := f_Z∘𝐡 is a GP so that f_X∼𝒢𝒫(m_M, k_M) with mean function given by m_M(𝐱) = m(𝐡(𝐱)) and the covariance function given by k_M(𝐱, 𝐱')  = k(𝐡(𝐱), 𝐡(𝐱')) respectively. Given high-dimensional training inputs 𝐗 and corresponding observations y⃗ of the objective function, we find model parameters {θ_h,θ_k} that maximize the marginal likelihood (evidence) {θ_h^*, θ_k^*}∈θ_h,θ_kmax  p(𝐲|𝐗, θ_h,θ_k).

This objective allows us to learn a low-dimensional embedding as a by-product of the supervised GP regression. 

Unsupervised dimensionality reduction usually solves an orthogonal task to that of learning a response surface. Algorithms, such as PCA or variational auto-encoders <cit.>, achieve compact data representations by optimizing objectives that are not necessarily useful in a supervised setting <cit.>.


The MGP, instead, leads to low-dimensional representations that are optimal (locally) for the regression task at hand.

We use a multi-layer feed-forward neural network with sigmoid activation functions as a feature map (encoder) h⃗, resulting in a feature space 𝒵=[0,1]^d. Neural networks as an explicit feature map within an MGP have already been applied successfully for modeling non-smooth responses in robot locomotion <cit.>. Deep networks have also proven useful for learning the orientation of images from high-dimensional images <cit.>. With a Gaussian likelihood, the MGP posterior predictive distribution at a test point 𝐱_⋆∈𝒳 is Gaussian distributed with mean and variance given by

    𝔼[f_X(𝐱_⋆)]    = m_M(𝐱_⋆) + k_M(𝐱_⋆, 𝐗)𝐊_My^-1(𝐲 - m_M(𝐗)) 
       = m(𝐳_⋆) + k(𝐳_⋆, 𝐙)𝐊_y^-1(𝐲 - m(𝐙)) 
    𝕍[f_X(x⃗_⋆)]   = k_M(𝐱_⋆, 𝐱_⋆) - k_M(𝐱_⋆,𝐗)𝐊_My^-1k_M(𝐗, 𝐱_⋆) 
       = k(𝐳_⋆, 𝐳_⋆) - k(𝐳_⋆,𝐙)𝐊_y^-1k(𝐙, 𝐳_⋆),

respectively, with  𝐳_⋆:=𝐡(𝐱_⋆) and 𝐙:=𝐡(𝐗). Moreover, k_M(𝐱_⋆, 𝐗) = k(𝐳_⋆, 𝐙)=[k(𝐳_⋆, 𝐳_i)]_i=1^N, where N is the size of the training dataset, 𝐊_My:=k_M(𝐗, 𝐗) + σ_n^2𝐈, 𝐊_y:=k(𝐙, 𝐙) + σ_n^2𝐈, k_M(𝐗, 𝐗):=k(𝐙, 𝐙), and m_M(𝐗):=m(𝐙)=[m(𝐳_i)]_i=1^N computes the prior mean function evaluated at the embedded training inputs Z⃗. 

Posterior predictions can be computed using both the feature and data space. Equations (<ref>)–(<ref>) appear in the definition of the acquisition functions in  (<ref>)–(<ref>) as mean μ(x⃗):=𝔼[f_X(𝐱)] and standard deviation σ(x⃗):=√(𝕍[f_X(x⃗)]) of the posterior predictions of the surrogate model.






The MGP defines a GP on 𝒳, but allows us to learn a response surface in the lower-dimensional feature space 𝒵. This is key for optimizing the acquisition function in a low-dimensional space 𝒵 instead of the original dataparameter space 𝒳. Thus far, we have detailed the feature-space BO procedure up to line <ref> in Algorithm <ref>. Once we found an optimizer z⃗_* of the acquisition function, we need to project it back into the original data space 𝒳 in order to evaluate the true objective f_X, whose domain is 𝒳. This can be done by means of a reconstruction mapping (decoder), which we detail in the following.
















 §.§ Input reconstruction with manifold multi-output Gaussian processes



In the following, we present the reconstruction part (decoder) of our feature space optimization model described in Figure <ref>.





We are interested in modeling the functional relationship between the feature space 𝒵 and the data space 𝒳 for step <ref> in Algorithm <ref>, which requires us to evaluate f_X. We therefore consider a vector-valued function 𝐠={g_i}_i=1^D, where each component g_i𝒵→𝒳_i maps vectors in feature space to the i-th coordinate of high-dimensional data, i.e. g_i(𝐳) = x̃^(i)∈𝒳_i. Multi-output GPs (MOGPs) <cit.> define a prior over vector-valued functions and explicitly allow for output correlations. An MOGP 𝒢𝒫(𝐦, 𝐊) is fully specified by a mean vector function 𝐦𝒵→ℝ^D and a positive, semi-definite matrix-valued covariance function 𝐊𝒵→ℝ^D× D, which computes the correlation between observations in the same output coordinate and cross-correlations between the D different outputs.


Here we consider the intrinsic coregionalization model (ICM) <cit.>, which structures the covariance matrix as a Kronecker product.
This model is particularly suitable for trading off number of model parameters and expressiveness of the vector valued function. In particular, the ICM facilitates information sharing across different tasks by adopting the same covariance function. It has been successfully adopted in robotics for learning inverse dynamics <cit.>. Hence, this model requires fewer parameters than the linear model of coregionalization <cit.> and allows for exploiting properties of the Kronecker product for efficient training and prediction. 

In our reconstruction model, we  need to ensure that the output space is exactly 𝒳. If we start from a data space 𝒳 = [0,1]^D the reconstructions need to belong to this hypercube. This property is not guaranteed by the MOGP. In order to satisfy this constraint, we consider a strictly monotonic output squashing function Ψ, as introduced in the context of warped GPs <cit.>. This allows us to define a corresponding inverse transformation Ψ^-1 that is applied to the data in input to the model. The resulting output of the MOGP at test time is then squashed through the transformation Ψ. Since the reconstruction of the MOGP is a distribution p(x̃⃗̃_⋆|X⃗,X̃⃗̃,z⃗_⋆), we evaluate the expectation with respect to this distribution of the transformed outputs, i.e. x⃗_t+1=𝔼_p[Ψ(x̃⃗̃_⋆)|X⃗,X̃⃗̃,z⃗_⋆]. In this paper, we choose the Gaussian cumulative density function as a monotonic squashing function Ψ :=Φ for warping the outputs of our reconstruction model <cit.>. The motivation for this choice is twofold: the inverse mapping Ψ^-1 is defined as the Probit function, which is a well known function, and the expectation 𝔼_p[Ψ(x̃⃗̃_⋆)|X⃗,X̃⃗̃,z⃗_⋆]
 with respect to the distribution p(x̃⃗̃_⋆|X⃗,X̃⃗̃,z⃗_⋆) at test reconstructions can be derived analytically <cit.>.










Intrinsic coregionalization model. The ICM <cit.> applies a linear mapping to a set of latent functions. In particular, we consider a set of P latent functions u_i𝒵→ℝ, that are assumed to be sample paths, i.e. sample functions independently drawn from the same GP prior 𝒢𝒫(m_c, k_c). The ICM model expresses the vector-valued function as a linear combination of these sample functions
𝐠(𝐳)=𝐀𝐮(𝐳), where 𝐮(𝐳)∈ℝ^P is the collection of the P sample paths' evaluations at  𝐳, and 𝐀∈ℝ^D× P is the linear mapping that couples the independent vector and parameterizes the ICM model. As a result, 𝐠 is an MOGP 𝒢𝒫(𝐦, 𝐊) with mean function 𝐦=𝐀𝐦_c, where 𝐦_c=[m_c]_i=1^P is obtained by repeating the single-valued mean function m_c in a P-vector. The covariance function is  𝐊(𝐳, 𝐳') = 𝐀𝐀^T⊗ k_c(𝐳, 𝐳'), 



where k_c is the covariance function for the GP prior, ⊗ is the Kronecker product and the matrix 𝐀𝐀^T is denoted as the coregionalization matrix. Note that k_c may differ from the covariance function k used for the response surface f_Z.

Reconstruction Model. For the reconstruction task in line <ref> of Algorithm <ref>, we introduce the manifold MOGP with intrinsic coregionalization model (mMOGP), which shares the feature map 𝐡 with the MGP used for learning the response surface; see Section <ref>. Without loss of generality, we assume a prior zero-mean vector function for the mMOGP 𝒢𝒫(0, 𝐁⊗ k_MO), where k_MO(𝐱, 𝐱') = k_c(𝐡(𝐱), 𝐡(𝐱')) and the matrix 𝐁=𝐀𝐀^T. 

We can interpret this model as an auto-encoder, where the MGP 𝐠∘𝐡:𝒳→𝒵 plays the role of the encoder, and the MOGP the role of the decoder, mapping low-dimensional features back into data space.










 §.§ Joint training


The joint training of the MGP, which models the response surface, and the mMOGP, which is used for the reconstruction (see also  Figure <ref>), is performed by maximizing a rescaled version of the log-marginal likelihood 


    ℒ∝   -𝐲^T𝐊_y^-1𝐲-log |𝐊_y|
        -1D(𝐱_V^T𝐊_V^-1𝐱_V + log |𝐊_V|) + const.

Here, ℒ comprises terms from both the MGP and mMOGP models,
where 𝐊_y is defined in  (<ref>), and the covariance matrix of the mMOGP 𝐊_V=𝐊̅ + σ_n^2𝐈 is obtained by evaluating the Kronecker product 𝐊̅=𝐁⊗ k_c(𝐙, 𝐙) with the mMOGP kernel k_c. The vector 𝐱_V is a concatenation of the columns of the data 𝐗. 
The maximizers [θ_h^*,θ_k^*,θ_c^*] of the log-marginal likelihood are the parameters θ_h^* of the feature map 𝐡 (which is shared between the MGP and the mMOGP), the hyper-parameters θ_k^* of the kernel k and the hyper-parameters θ_c^* of k_c including the coregionalization matrix B⃗ for the mMOGP, respectively. The rescaling factor 1/D balances the contributions of the two log-marginal likelihood terms involved in  training. The dimensions of the matrix 𝐊_V are ND× ND which correspond to repeating the K⃗_y matrix D times in a block-diagonal fashion. This block diagonal would then have quadratic form equal to Dy⃗^TK⃗_y^-1y⃗ and log determinant equal to Dlog |𝐊_y|. Thus, an equivalent rescaling is to divide the reconstruction terms 𝐱_V^T𝐊_V^-1𝐱_V and log |𝐊_V| by D. Optimization of (<ref>) is performed via gradient-based methods <cit.>.


Modeling the black-box objective function f_X is orthogonal to the reconstruction problem. However, when training these tasks jointly, they have a regularization effect on the optimization of the parameters θ_h of the feature embedding in the sense that the mapping 𝐡 will not overfit to a single regression task: the parameters θ_h will give rise to a feature space embedding that is useful for both the modeling of the objective and the reconstruction of the original inputs.



The major computational bottleneck for evaluating the marginal likelihood comes from the term 𝐱_V^T𝐊_V^-1𝐱_V, which requires inverting an ND× ND covariance. We reduce the computational complexity of this operation to 𝒪(N^3) + 𝒪(D^3) by exploiting the properties of the Kronecker product, tensor algebra <cit.> and structured GPs <cit.> as shown in the following section. 



 §.§ Computationally efficient mMOGP



For the reconstruction mapping 𝐠, we use the posterior mean of the mMOGP with intrinsic coregionalization model and apply exact inference and training via rescaled marginal likelihood maximization. While the ICM enables modeling correlation between arbitrary pairs of dimensions, it also requires computing  a Kronecker product to evaluate the full covariance matrix of all outputs


    𝐊̅ = 𝐁⊗ k_c(𝐙, 𝐙),

where k_c(𝐙, 𝐙) is the covariance matrix obtained from the training inputs 𝐙 in feature space and the 𝐁 matrix is the coregionalization matrix of the ICM. Inverting the full covariance matrix 𝐊̅ requires 𝒪(N^3D^3) and easily becomes intractable in high-dimensional spaces even for small N. Storing this full covariance matrix  required 𝒪(N^2D^2) space and also becomes challenging in high dimensions. For an efficient implementation of the ICM, we exploit properties of the Kronecker product and apply results from structured GPs <cit.> that allow for efficient training and predictions in 𝒪(N^3) + 𝒪(D^3) time and 𝒪(ND) space. In particular, we are interested in the full covariance matrix under the assumption of a Gaussian likelihood for the multi-output observations, i.e. 𝐊̅ + σ_n^2𝐈. We first express the full covariance matrix in terms of its eigendecomposition, i.e. 𝐊̅ = 𝐐Λ𝐐^T. This allows expressing the inverse of the covariance from noisy targets as

    (𝐊̅ + σ_n^2𝐈)^-1 = 𝐐(Λ + σ_n^2𝐈)^-1𝐐^T,

where both Λ and σ_n^2𝐈 are diagonal and can be trivially inverted. However, the eigendecomposition of an ND× ND matrix would still be cubic in the product between the number of dimensions and the number of data. By the properties of the Kronecker product, we can express  the eigendecomposition itself with a Kronecker structure, i.e.

    ⊗_l=1^W 𝐊_l = ⊗_l=1^W 𝐐_l⊗_l=1^W Λ_l(⊗_l=1^W 𝐐_l)^T,

where each term of the Kronecker product on the left-hand side 𝐊_l∈ℝ^G_l× G_l has eigendecomposition 𝐊_l=𝐐_lΛ_l𝐐_l^T for l=1,...,W, where W is number of factors in the Kronecker product. In our ICM model W=2, because  the coregionalization matrix 𝐁 Kronecker multiplies k_c(Z⃗, Z⃗), the covariance matrix of the observations; see (<ref>). Thus, from  (<ref>)–(<ref>), we are allowed to invert the covariance from noisy targets by separately decomposing the covariance matrix k_c(𝐙,𝐙)=𝐐_kΛ_k𝐐_k^T and the coregionalization matrix 𝐁=𝐐_bΛ_b𝐐_b^T, which require 𝒪(N^3) and 𝒪(D^3) time, respectively; see line <ref> of Algorithm <ref>.



Storing this inverse matrix and multiplying it by a vector still requires 𝒪(N^2D^2) space and run time, respectively, so that this step becomes the main bottleneck for efficient mMOGP training and predictions computation. To address this issue we represent the expensive matrix-vector multiplication as a sequence of small matrix-tensor multiplications without computing the full Kronecker product  <cit.>. In particular, we are interested in efficiently evaluating

    𝐫 = (⊗_i=1^W 𝐊_i)𝐱.

We first represent the multiplication of a matrix with Kronecker structure by a vector as a tensor product. A tensor 𝐓_i_1,...,i_V can be interpreted as an extension of matrices to objects where elements are indexed using a set of V indices: i_1,...,i_V, where the number V is referred to as the order of the tensor.

With the definition of the Kronecker product we express the left-hand side of (<ref>) as a tensor

    [⊗_l=1^W 𝐊_l]_i, j   = [𝐊_1]_i_1, j_1· ... ·[𝐊_W]_i_W, j_W,
    
        1≤ i_l   ,j_l≤ G_l ,  
        1≤ i,j≤∏_l=1^WG_l.


The right-hand side of (<ref>) coincides with a tensor 𝐓^K_i_1,j_1,...,i_W,j_W, and a similar tensor-representation can be obtained for the ∏_l=1^WG_l-long vector 𝐱, i.e. 𝐓^X_j_W,...,j_1.
A tensor product between the tensors 𝐓^K_i_1,j_1,...,i_W,j_W and 𝐓^X_j_W,...,j_1 applies a contraction along the indices of the second tensor, i.e.

    ∑_j_1...∑_j_W 𝐓^K_i_1,j_1,...,i_W,j_W𝐓^X_j_W,...,j_1.

This tensor contraction can be expressed in terms of a sequence of tensor-transposed matrix-tensor products

    (⊗_l=1^W 𝐊_l)𝐱   = vec((𝐊_1⋯(𝐊_W𝐓^X)^⊤)^⊤)

with 𝐊_l𝐓^X = ∑_k=1^G_l [𝐊_l]_i_1, k𝐓^X_k, j_2,...,j_W. The function vec(·) returns the vectorized form of a matrix by stacking its columns vertically. The tensor transposition ⊤ applies a cyclic permutation to the order of the indices in a tensor. As a result, the right-hand side in (<ref>) allows us to evaluate the expensive matrix-vector product without computing and storing the Kronecker product. Algorithm <ref> shows the main steps of the efficient matrix inversion and matrix-vector multiplication for matrices that feature a Kronecker structure. The matrix vector multiplication subroutine is expressed as a sequence of tensor-transpose matrix-tensor products.





















































































§ CONSTRAINED ACQUISITION



We defined a joint probabilistic model for the response surface learning and the input reconstruction tasks; see lines <ref>–<ref> and <ref> of Algorithm <ref>, respectively. We are now concerned with the maximization of the acquisition function in feature space; see line <ref> of Algorithm <ref>. We aim at maximizing the acquisition function in a low-dimensional feature space of the original dataparameter space 𝒳. However, one problem that arises with the mMOGP decoding is that locations in feature space, which are too far away from data, will be mapped back to the mMOGP prior. 

Since the acquisition function is a key driver of exploration in BO, this is a problem. We address this limitation by introducing a constraint based on the Lipschitz continuity of the mMOGP posterior. This will ensure that candidates z⃗_*∈𝒵 selected in feature space will not collapse to the origin 0∈ℝ^D if the reconstruction is defined as x̃⃗̃_*=μ(z⃗_*), where μ is the posterior mean of the mMOGP.



We want to leverage information from observed data for the multi-output mapping and exploit it when optimizing the acquisition function in feature space. This can be achieved by upper-bounding the Euclidean distance


    dist(𝐳, 𝐙_t) = 1≤ i≤ N_tmin 𝐳_i-𝐳_2

in feature space
 between the optimization variable 𝐳 and the embedded training data 𝐙_t=[𝐳_1,...,𝐳_N_t]. Here, N_t is the number of data points available at BO iteration t. The desired upper bound is obtained by exploiting the Lipschitz continuity property of the multi-output posterior mean for which 


    |[μ(𝐳)]_i - [μ(𝐳')]_i|≤ L 𝐳-𝐳'.

Here, L denotes the Lipschitz constant of the posterior mean μ of the mMOGP. For common kernels, such as Matérn_52 and squared exponential, the posterior mean is Lipschitz continuous. The upper bound

    dist(𝐳, 𝐙_t)≤μ_max(𝐳^*)/L

allows us to specify how far from the data we can move in feature space without falling back to the prior on all coordinates of the reconstruction. Here 𝐳^* minimizes the distance in (<ref>), while the numerator on the right-hand side is the component-wise maximum of μ(𝐳^*).




We estimate the Lipschitz constant as the maximum norm of the Jacobian of the posterior mean of the mMOGP <cit.>

    L = 𝐳∈𝒵max ∇_z⃗μ(𝐳).

This maximization returns a valid Lipschitz constant <cit.> for the multi-output mapping for any choice of norm in  (<ref>). The Jacobian of the posterior mean is represented by a D× d matrix and we adopt the max norm ∇_z⃗μ(𝐳)_∞=max |μ_i,j^'| for i=1,...,D and j=1,...,d. 




Lower values of valid Lipschitz constants L allow for exploration in larger regions of the feature space that still satisfy the nonlinear constraint in  (<ref>).



§ EXPERIMENTS




We report results on a set of high-dimensional benchmark functions that possess an intrinsic low dimensionality. In particular, we (i) assess the benefits of adopting a model structure as presented in Figure <ref>; (ii) analyze the benefits of the constrained optimization of the acquisition function. Our purpose is to compare empirical performances across (a) different characterizations of the feature spaces, e.g. linearnonlinear subspaces; (b) different properties of the objective function, e.g. additivitynon additivity; (c) a real problem set. 











Approaches. We compare our approach (MGPC-BO) with the random embeddings optimization (REMBO) <cit.>, which performs BO on a random linear subspace of the inputs. Additional baselines include additive models (ADD-BO) <cit.>, which assumes an additive structure (across dimensions) of the objective f_X, and one recently proposed VAE-based model (VAE-BO) <cit.> that learns a feature space with deep networks offline. We also include a version of our model presented in Figure <ref> (HMGPC-BO) that uses a hierarchical ICM for the input reconstruction mapping 𝐠. The hierarchical ICM partitions the data space into low-dimensional disjoint subsets, i.e. {𝒳_i}_i=1^Q, 𝒳_i⊂ℝ^3, and assumes independence between reconstructions of different subsets, i.e. 𝐱̃^(i)⊥𝐱̃^(j), where 𝐱̃^(i)∈𝒳_i,
𝐱̃^(j)∈𝒳_j for i ≠ j. Moreover, the baselines MGP-BO and HMGP-BO correspond to same modeling as in MGPC-BO and HMGPC-BO, respectively, but without applying the nonlinear constraint in  (<ref>). We also compare with a different parametrization of the covariance function of the decoder 𝐠. The baseline DMGP-BO and DMGPC-BO define a single kernel k_c for the reconstruction task while HMGP-BO and HMGPC-BO define different kernels {k_c^i}_i=1^Q, one for each subset of the partitioning. Here, DMGPC-BO and DMGP-BO denote the baseline with and without Lipschitz regularization, respectively. For all the approaches we specify the dimensionality d_fs of a feature space where the optimization is performed. Note that this value may differ from the intrinsic dimensionality d of the objective functions i.e. d_fs≠ d.



Acquisition functions. We evaluate the performances of all baselines across common acquisition functions: EI <cit.>, UCB <cit.> and PI <cit.>, which are also given  in (<ref>)–(<ref>). The motivation in selecting the above acquisition functions is that we wish to explore performances of our BO approach on a range of different decision strategies: aggressive exploitation (PI), aggressive exploration (UCB) and one-time-step optimal selection (EI). 

In our experiments we set the β_t parameter of UCB in (<ref>) to √(3). Moreover, we do not have access to the true f_min required in  (<ref>) and (<ref>). Therefore, we compute the improvement based acquisitions (EI and PI) using y_min:=miny⃗_t, which is the best noisy observation obtained up to iteration t.

The maximization of the acquisition function is identical for all baselines: we first perform a random search step with 5,000 samples drawn uniformly at random and select the best 100 locations to apply gradient-based optimization from these starting locations. For box-constrained acquisition optimization we use L-BFGS-B  <cit.>. For constrained acquisition optimization with nonlinear constraints we use a trust-region interior point method <cit.>.

Model parameters. In our experiments, we select the Matérn_5/2 kernel as the covariance function for the GPs in each baseline. For the neural network employed in the encoder, the architecture was a single hidden layer with 20 units, and as the activation function we use the sigmoid activation 1/(1+exp(-x)).



Experiment setup. Each BO progression curve shows the mean and standard error of the immediate logarithmic regret log_10|f(𝐱_best(t))-f_min|, where f_min is the true minimum of f_X and 𝐱_best(t)∈min_i=1:t f_X(𝐱_i). Mean and standard error are computed over 20 experiments with different random initializations. All optimization experiments start with a budget of 10 data points and perform a total of 300 iterations. The noise variance is σ^2_n=10^-4.  








 §.§ Linear feature space



We consider benchmark functions that are defined in a d=10-dimensional space. We map their input space to a D=60-dimensional space using an orthogonal matrix R⃗^d× D so that the overall objective is f_X(𝐱)=f(z⃗) = f(𝐑𝐱). 




  §.§.§ Additive objective



We minimize the Rosenbrock benchmark function 

    f(𝐳) = ∑_i=1^d-1[100(z_i+1-z_i^2)^2 + (z_i-1)^2]

in a d_fs=10-dimensional feature space. Figure <ref> shows that HMGPC-BO baseline descends quickly to relatively low regret in the early stages of optimization and recovers better regret at termination than the unconstrained baseline HMGP-BO. 

The VAE-BO baseline improves quickly but lacks exploration due to an insufficiently expressive reconstruction mapping from feature space to data space. We highlight that the VAE-BO model was trained on a budget of 500 inputs-observations pairs prior to starting the BO experiments. This additional budget, however, still does not allow the VAE-BO to compare well with baselines that learn a feature mapping during optimization. REMBO shows a competitive descent for two main reasons: the fact that the baseline conforms with the linear embedding assumption that characterizes the objective function and the employment of an orthonormal linear mapping which is supposed to improve performances and conforms to structural assumption about the linear embedding R⃗. The ADD-BO baseline suffers from the coupling effects of the linear dimensionality reduction 𝐑. Overall, Figure <ref> highlights the fast learning of feature space representations that are effective for optimization with MGPC-BO, HMGPC-BO and DMGPC-BO baselines.












  §.§.§ Non-additive objective




    

Here, we optimize the product of sines with intrinsic dimensionality d=10 

    f(𝐳) = 10sin(z_1) ∏_i=1^dsin(z_i)

and compare results when the additivity assumption is not satisfied. Figure <ref> shows the regret curves obtained optimizing the objective on a d_fs=10-dimensional feature space. Solid lines describe the Lipschitz-regularized baselines MGPC-BO, HMGPC-BO and DMGPC-BO (with nonlinear constraint), while dashed lines are baselines that apply box-constrained maximization of the acquisition in feature space. The HMGP-BO, MGP-BO and DMGP-BO regrets flatten early in both improvement-based acquisition functions (EI and PI) since these acquisition functions highlight locations in feature space that are too far away from the training data. In this setting, the decoder 𝐠 returns the same high-dimensional reconstruction, which prevents BO from exploring. The constrained maximization of the acquisition is beneficial for all our models. We also note that the REMBO baseline conforms to the intrinsic linear low-dimensionality assumption described in Section <ref> and is  the most competitive baseline especially for UCB acquisition. However, the linear reconstruction mapping applied by REMBO also suffers from non-injectivity, and this slows down exploration in the high-dimensional space. The linear projection  deteriorates performances of the additive model. ADD-BO assumes independence between axis-aligned projections of the high-dimensional space, while the linear mapping 𝐑 couples all subsets of dimensions. This linear mapping, therefore, penalizes optimization with independent additive components. The VAE-BO approach requires much larger amounts of data to learn a meaningful reconstruction mapping than available in our experiment. Thus, most locations in feature space are mapped to similar reconstructions. This explains the flat curve observed on all VAE-BO progressions with different acquisitions. 




 §.§ Nonlinear feature space with non-additive objective


We consider the product of sines functions and apply a nonlinear dimensionality reduction. We define a single-layer neural network mapping to elevate the dimensionality of the objective to D=60, i.e. f_X(𝐱)=f(γ(𝐑𝐱)). Here γ is the sigmoid activation function. We select a dimensionality of the feature space as in previous sections d_fs=10 which is equal to the intrinsic dimensionality of the objective function d=10. Figure <ref> shows the progression of the regret over 300 BO iterations. We can observe consistent improvements of MGPC-BO, HMGPC-BO and DMGPC-BO with respect to VAE-BO which also assumes a nonlinear embedding for the objective.  The performance of MGPC-BO, HMGPC-BO and DMGPC-BO also retain better regret at termination than with box-constrained acquisition maximization (MGP-BO, HMGP-BO, DMGP-BO). Here we apply a significance testing with the Wilcoxon signed-rank test <cit.> at termination of the optimization between the best performing of our baselines, namely DMGPC-BO and the best competitive baseline that is ADD-BO. We observe a significance of at least 0.014% for all acquisition functions (largest p-value p=0.00014 for UCB acquisition) meaning that our best baseline DMGPC-BO is highly significantly different than the ADD-BO baseline and attains better regret than ADD-BO at termination of the optimization. 


Overall, we observe that the constrained maximization of the acquisition function is beneficial for the proposed model. The advantages with respect to ADD-BO, REMBO and VAE-BO baselines are more evident with the product of sines objective with nonlinear embedding while with the Rosenbrock we retain no worse regret.





















 §.§ Sensitivity analysis on real data


Here we apply a sensitivity analysis with respect to the dimensionality of the feature space d_fs on a D=12-dimensional real problem. We consider the Thomson problem of finding the lowest potential configuration of a set of electrons on a sphere <cit.>. This is a central problem in physics and chemistry for identifying a structure with respect to atomic locations <cit.>. The potential of a set of n_p electrons on a unit sphere is given by the objective

    ∑_i=1^n_p-1∑_j=i+1^n_p((x_i-x_j)^2+(y_i-y_j)^2+(z_i-z_j)^2)^-1/2.

This is a constrained minimization problem with constraints x_i^2+y_i^2+z_i^2=1 for i=1,...,n_p, which means that all electrons must lie on a unit sphere. We represent the variables of the problem as spherical coordinates with unit radius. This allows us defining two variables per point with a total number of 2n_p (azimuthal and polar angles) parameters to optimize within box constraints. For optimization, we select n_p=6, which results in a D=12-dimensional problem and we optimize it on low-dimensional feature spaces of dimensionalities d_fs=6, 4, 3, 2 to observe the effect of this hyper-parameter in the optimization. Figure <ref> shows a comparison of our approaches with ADD-BO, REMBO and VAE-BO baselines on a single acquisition function PI. Overall, we observe a deterioration of performances with diminishing dimensionality of the feature space. The regret clearly increases for our baselines when we select d_fs=2 meaning that, with a high compression rate, the probabilistic model for MGPC-BO, DGPC-BO and HMGPC-BO learns less useful features for optimization. We observe the most competitive baseline to be ADD-BO, which decomposes the 12-dimensional problems into D/d_fs sub-problems with dimensionality d_fs. Another competitive baseline is REMBO. We apply a significance test and compare our nonlinearly constrained baseline MGPC-BO with the most competitive baseline (ADD-BO or REMBO) for each plot of Figure <ref> at termination of the optimization. We select the Wilcoxon signed-rank test <cit.>, which does not assume that the difference between the sample populations is Gaussian. 
For feature space dimensionality d_fs=2 we do not observe values significantly different since the p-value is p=0.135. This is due to the deterioration of performances at d_fs=2. For d_fs=3 we observe a more significant difference between MGPC-BO and ADD-BO with p-value p≤ 0.002. With hyper-parameter values d_fs≥ 4 we observe significantly different baselines with significance at 0.6% (difference between MGPC-BO and ADD-BO for d_fs=4 with p-value p≤ 0.003 and between MGPC-BO and REMBO for d_fs=6 with p-value p≤ 0.006).
Overall, we observe our constrained baselines to perform better than ADD-BO and REMBO and to reach the lowest value in notably less BO iterations.



 §.§ Run-time complexity



The computational complexity of MGPC-BO is 𝒪(D^3+N^3) due to the eigen-decomposition of both the coregionalization (D^3) and kernel matrix (N^3). The baseline HMGPC-BO scales with 𝒪(d_out^3Q+N^3Q) with Q being the number of independent subsets of dimensions, i.e. Q=D/d_out, with d_out being a small constant value (d_out=3). This baseline achieves faster computations when having small number of data points N, for large number of data points and large number of dimensions (both tending to infinity) the MGPC-BO results more efficient. The baseline DMGPC-BO instead has complexity 𝒪(d_out^3Q+N^3), which is  faster than the MGPC-BO. MGP-BO, DMGP-BO and HMGP-BO have the same complexity of MGPC-BO, DMGPC-BO and HMGPC-BO, respectively. The remaining baselines have all computational complexity 𝒪(N^3) due to the matrix inversion of the covariance matrix for GP training which is used in ADD-BO, REMBO and VAE-BO. Our baseline has an additional overhead of at least a linear term d_out^3Q, which implies slower training times for our probabilistic model. This is a reasonable trade off for improved optimization performances and better data efficiency in our reconstruction model.




§ CONCLUSION


We proposed a framework for efficient Bayesian optimization of intrinsically low-dimensional black-box functions based on nonlinear embeddings. In our model, a manifold GP learns useful low-dimensional feature representations of high-dimensional data by jointly learning the response surface and a reconstruction mapping. 

Our approach allows for optimizing acquisition functions in a low-dimensional feature space. 
Since exploration in feature space (driven by the acquisition function) does not necessarily mean exploration in the high-dimensional parameter space, we introduce a nonlinear constraint based on Lipschitz continuity of predictions of the reconstruction mapping, which encourages exploration in the vicinity of the training data and mitigates un-identifiability issues in data space, which hinder optimization. 




§ ACKNOWLEDGMENTS

We thank James T. Wilson for valuable feedback on early drafts of the manuscript. This work has been supported by the EPSRC Centre for Doctoral Training in High Performance Embedded and Distributed Systems (HiPEDS, Grant EP/L016796/1) and the Data Science Institute, Imperial College London.




spmpsci      

   



