BOtied: Multi-objective Bayesian optimization with tied multivariate ranks

Bayesian optimization (BO) has demonstrated promise in a variety of scientific and industrial domains where the goal is to optimize an expensive black-box function using a limited number of potentially noisy function evaluations (Romero et al., 2013; Calandra et al., 2016; Kusne et al., 2020; Shields et al., 2021; Zuo et al., 2021; Bellamy et al., 2022; Park et al., 2022). In BO, we fit a probabilistic surrogate model on the available observations so far. Based on the model, the acquisition function determines the next candidate to evaluate by balancing exploration (evaluating highly uncertain candidates) with exploitation (evaluating designs believed to be optimal). Often, applications call for the joint optimization of $M{\rm\geq}2$ multiple, potentially competing objectives (Marler & Arora, 2004; Jain et al., 2017; Tagasovska et al., 2022). Unlike in single-objective settings, a single optimal solution may not exist and we must identify a set of solutions that represents the best compromises among the multiple objectives. The acquisition function in multi-objective Bayesian optimization (MOBO) navigates these trade-offs as it guides the optimization toward regions of interest.

Many bona fide MO acquisition functions without scalarization involve high-dimensional integrals and scale poorly with increasing $M$. Moreover, improvement-based acquisition functions including some variants of random scalarization are sensitive to non-informative monotonic transformations of the objectives. This is a pain point for many practical applications. For instance, in biochemistry, the dissociation constant $K_{D}$ is typically expressed in terms of its log transformation, the ${\rm p}K_{D}\equiv{\rm-}\log_{10}(K_{D})$. It would be desirable to work with acquisition functions that are invariant to the choice of such unit conversions.

To address these challenges, we propose BOtied¹¹1The name stems from non-dominated solutions considered to be ”tied.”, a novel acquisition function based on multivariate ranks. We show that BOtied has the desirable property of being invariant to relative rescaling or monotonic transformations of the objectives. While it maintains the multivariate structure of the objective space, its implementation has highly favorable time complexity and we report wall-clock time competitive with random scalarization.

In Figure 1, we present the intuition behind BOtied. Consider a minimization setup for $M{\rm=}2$ where we seek to identify solutions on the true Pareto front (solid red curves), unknown to us. Suppose we have many candidates, represented as circular posterior blobs in the objective space, where the predictive distributions have been output by our surrogate model. For simplicity, assume the uncertainties (sizes of blobs) are comparable among the candidates. How do we estimate the quality of our solutions (i.e., each candidate set’s proximity to the true Pareto front)? One quality indicator is the hypervolume (HV; Zitzler et al., 2003), defined as the size of the polytope bounded from above by a predefined reference point and dominated by the candidate set (shaded areas in the leftmost panel). When one candidate set dominates another, its HV is greater. We can visually confirm that the HV of the blue Pareto approximation is greater than that of the orange.

Next, let us adopt the related, but distinct, perspective of multivariate ranking (see, e.g., Ghosal & Sen, 2022). Define the random objective vector $Y=[f_{1}(X),\dotsc,f_{M}(X)]$ taking values in $\mathbb{R}^{M}$, which is the result of applying the objective function $f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{M}$ on the random design variable $X$ that takes values in $\mathbb{R}^{d}$. How do we compare two realizations of $Y$, say ${\bm{y}}$ and ${\bm{y}^{\prime}}$?

Ranking vectors is non-trivial, as there is no natural ordering in Euclidean spaces when $M{\rm\geq}2$. We propose to use the joint cumulative distribution function (CDF), defined as the probability of ${\bm{y}}$ being weakly dominated: $F_{Y}({\bm{y}})=\mathbb{P}(Y_{1}{\rm\leq}y_{1},\dotsc,Y_{M}{\rm\leq}y_{M})$, where ${\bm{y}}=[{y}_{1},\dotsc,{y}_{M}]\in\mathbb{R}^{M}$. The CDF formalizes the rank ordering of vectors as weak dominance in the joint minimization of $M$ objectives (Binois et al., 2015). Specifically, The CDF scores and their $\alpha$-level lines $L_{\alpha}^{F_{Y}}=\{{\bm{y}}:F_{Y}({\bm{y}})=\alpha\}$ are depicted in the middle panel of Figure 1 for multiple values of $\alpha\in[0,1]$. All candidates with equal multivariate rank, or ties, lie on the same level line, as shown in the rightmost panel.

Multivariate ranking via the CDF can be understood in relation to the associated probability density function (PDF), as the CDF is the integral of the PDF. Figure 2 shows the CDF and the associated PDF side by side for a bi-objective setting ($M{\rm=}2$). The right panel shows the PDF fit on 200 outcome samples (gray dots) via kernel density estimation, where the outcome samples were drawn from an elliptical Gaussian. The left panel shows the level lines of the corresponding CDF. The $\alpha$-level lines converge to the approximate Pareto front as $\alpha\rightarrow 0$. The lowermost level line ($\alpha{\rm\approx}0$) closely traces the convex shape of the true Pareto front shown as the solid red curve.

MO indicators and acquisition functions. A computationally attractive approach to MOBO scalarizes the objectives with random preference weights (Knowles, 2006; Paria et al., 2020) and applies a single-objective acquisition function. The distribution of the weights, however, may be insufficient to encourage exploration when there are many objectives with unknown scales.

Alternatively, we may preserve the MO structure by seeking improvement on a set-based performance metric, such as the HV indicator (Embrechts et al., 2003) or the R2 indicator (Deutz et al., 2019a, b). Improvement-based acquisition functions such as the expected hypervolume improvement (EHVI; Emmerich et al., 2011; Daulton et al., 2020, 2021) are sensitive to the rescaling of the objectives, which may carry drastically different natural units. In particular, computing the HV has time complexity that is super-polynomial in the number of objectives, because it entails computing the volume of an irregular polytope (Yang et al., 2019). Despite the efficiency improvement achieved by box decomposition algorithms (Dächert et al., 2017; Yang et al., 2019), HV computation remains slow when $M{\rm>}4$.

Another class of acquisition strategies is entropy search, which focuses on maximizing the information gain from the next observation (Villemonteix et al., 2009; Hennig & Schuler, 2012; Hoffman & Ghahramani, 2015; Shah & Ghahramani, 2015; Hernández-Lobato et al., 2016b; Belakaria et al., 2019; Tu et al., 2022). Entropy searches are commonly implemented in box decompositions as well, but are costly to evaluate without using more tractable bounds to serve as approximations.

Multivariate ranking. The scale-invariant properties of ranking makes it an attractive tool for optimization. Binois et al. (2015) relates the Pareto front to the extreme level line of the CDF, $F_{Y}$. Considering ranking dominance as an alternative to Pareto dominance, Kukkonen & Lampinen (2007) propose computing the ranks of individual objectives separately and combining them post-hoc with a simple aggregation function (min, max, average) to obtain the overall fitness value for a given candidate. Binois et al. (2020) explores the question of how to choose from the set of non-dominated solutions, which grows with $M$, and makes a game-theoretic argument for how to make the compromise. In particular, they define trade-offs in the copula space, which is the scale-invariant rank transformation of the original objective function. For single-objective BO, Picheny et al. (2019) propose Ordinal BO, which uses a Gaussian process (GP) surrogate that is only sensitive to the rankings of the inputs and objective values. Notably, this method can robustly handle ill-conditioned and multi-modal distributions in the objective function values, for which GP models are known to often fail. Similarly, Eriksson & Poloczek (2021) use the rank transformations to magnify values at the end of the observed ranges. To our knowledge, however, our work is the first to incorporate multivariate rankings enabled by the joint CDF into a MOBO algorithm. We explicitly account for the structure of the $M$-variate objective distribution in identifying the full Pareto front.

A more detailed overview and positioning of BOtied with respect to the MO literature can be found in Appendix A.

In MOBO, it is common to evaluate the quality of an approximate Pareto set ${\mathcal{X}}$ by computing its distance from the optimal Pareto set $\mathcal{X}^{*}$ in the objective space, defined by some distance metric $d:2^{\mathcal{Y}}\times 2^{\mathcal{Y}}\rightarrow\mathbb{R}$ where $2^{\mathcal{Y}}$ denotes the power set of the objective space $\mathcal{Y}$. HVI (Equation 3) is a popular metric, for instance. One advantage of HV is its sensitivity to any type of improvement; whenever an approximation set $A$ dominates another approximation set $B$, then the measure yields a strictly better quality value for the former (Zitzler et al., 2003). On the other hand, HV suffers from sensitivity to transformations of the objectives and scales super-polynomially with $M$, which hinders its practical value. An alternative approach is to use distance-based indicators (Miranda & Von Zuben, 2016; Shilton et al., 2018) that assign scores for the solutions based on a signed distance from each point to the approximate Pareto front, which is again computationally expensive.

In the following, the (weak) Pareto-dominance relation is used as a preference relation $\preccurlyeq$ on $\mathcal{Y}$ indicating that a solution ${\bm{y}^{\prime}}$ is at least as good as a solution ${\bm{y}}$ (denoted ${\bm{y}^{\prime}}\preccurlyeq{\bm{y}}$) iff $f_{i}({\bm{y}^{\prime}})\leq f_{i}({\bm{y}})\ \forall i\in[M]$. This relation can be canonically extended to sets of solutions where a set $A\subseteq X$ weakly dominates a set $B\subseteq X$ (denoted $A\preccurlyeq B$) iff $\forall{\bm{y}}\in B\ \exists{\bm{y}^{\prime}}\in A:{\bm{y}^{\prime}}\preccurlyeq{\bm{y}}$ (Zitzler et al., 2003). Given the preference relation, we consider the optimization goal of identifying a set of solutions that approximates the set of Pareto-optimal solutions and ideally this set is not strictly dominated by any other approximation set.

Since the generalized weak Pareto dominance relation defines only a partial order on $\mathcal{Y}$, there may be incomparable sets in $\mathcal{Y}$. Incomparability is a key challenge in search and performance assessment for multi-objective optimization and becomes more serious as $M$ increases (Fonseca et al., 2005). One way to circumvent this problem is to define a total order on $\mathcal{Y}$ which guarantees that any two objective vector sets are mutually comparable. To this end, quality indicators have been introduced that assign, in the simplest case, each approximation set a real number — that is, a (unary) indicator function $I:\mathcal{Y}\rightarrow\mathbb{R}$ (Zitzler et al., 2003). One important feature an indicator should have is Pareto compliance (Fonseca et al., 2005), which dictates that it must not contradict the order induced by the Pareto dominance relation.

In particular, this means that whenever $A\preccurlyeq B\land B\nsucceq A$, then the indicator value of A must not be worse than the indicator value of B. A stricter version of compliance would be to require that the indicator value of A is strictly better than the indicator value of B (if better means a lower indicator value):

Experimental setup. To empirically evaluate the sample efficiency of BOtied, we execute simulated BO rounds on a variety of problems. See Appendix F for more details about our setup. For all the experiments, the surrogate model was an independent GP with a Matern 5/2 ARD kernel. The GP hyperparameters were inferred via maximum a posteriori (MAP) estimation. The code that reproduces all of our experiments and plots is available at https://github.com/jiwoncpark/botied \faGithub.

Metrics. We use the HV indicator presented in section 4, a standard evaluation metric for MOBO, as well as our CDF indicator $I_{\rm CDF}$ on the noiseless function values. We rely on efficient algorithms for HV computation based on hyper-cell decomposition as described in (Fonseca et al., 2006; Ishibuchi et al., 2011) and implemented in BoTorch (Balandat et al., 2020).

BOtied implementation We implement two versions of BOtied that differ in the incorporation of predictive uncertainties in the CDF estimation. In one version (v1), we fit the CDF on all of the MC predictive posterior samples across all the candidates. This can sometimes result in poor CDF fit, particularly when uncertainties are large. The other version (v2) alleviates this issue by fitting the CDF on the posterior means of the candidates. The algorithms for both versions can be found in Algorithm 1, Appendix B. We optimize the BOtied acquisition values using the gradient-free CMA-ES algorithm (Hansen, 2006). The CDF estimation is detailed in Appendix E and BOtied optimization in Appendix F.

Baselines We compare BOtied with the noisy versions of popular acquisition functions. The baseline acquisition strategies are NEHVI (Daulton et al., 2020) described in Equation 4; noisy NParEGO (NParEGO; Knowles, 2006) which uses noisy EI on top of random augmented Chebyshev scalarization; predictive entropy search (PES; Hernández-Lobato et al., 2016a), maximum entropy search (MES; Belakaria et al., 2019), and joint entropy search (JES; Hvarfner et al., 2022) — the differences being the estimation of entropy in the inputs, objectives, or both, respectively; and random (Sobol) selection.

Synthetic datasets. We include synthetic test functions for direct evaluation of $f$. We focus on ones that support $M\geq 3$: DTLZ2 ($d{\rm=}6,M{\rm=}4$ and $d{\rm=}7,M{\rm=}6$; Deb & Gupta, 2005) and Penicillin ($d{\rm=}7,M{\rm=}3$; Liang & Lai, 2021), which simulates the penicillin yield, time to production, and undesired byproduct for various parameters of the production process. See Appendix F for more detail.

Real-world datasets. To emulate a multi-objective drug design setting, we postprocess the real-world dataset Caco2 (Wang et al., 2016) from the Therapeutics Data Commons database (Huang et al., 2021, 2022) to create Caco2+. The original Caco2 dataset consists of 906 drug molecules annotated with experimentally measured cell permeability, or the rates of passing through a human colon epithelial cancer cell line. Permeability is a key property in the absorption, distribution, metabolism, and excretion (ADME) profile of drugs. We augment the dataset with additional properties using RDKit (Landrum et al., 2023), including ClogP related to fat solubility and topological polar surface area (TPSA). Subsets of these properties (e.g., permeability and TPSA) are inversely correlated and thus compete with one another during optimization. In late-stage lead-molecule optimization, the trade-offs become more dramatic and as more properties are added (Sun et al., 2022). Demonstrating effective sampling of Pareto-optimal solutions in this setting is thus of great value. We represent each molecule as a concatenation of fingerprint and fragment feature vectors (Thawani et al., 2020).

We also include experiments over three datasets from the DDMOP benchmark (He et al., 2020). Differently from the synthetic test functions which have analytical solutions, each DDMOP dataset represents a complex objective function approximated by expensive numerical simulations. These datasets address cab car optimization, power system chip placement and neural network. Details and table results on each dataset can be found in subsection F.5. See Appendix F for more detail about these datasets.

Vine Copulas for MOBO in practice In Figure 4, we present a simple recipe for estimating CDFs with vine copulas, in three simple steps and fewer than five lines of code. We use the Caco2+ dataset ($M{\rm=}3$) as an example. First, the PIT transformation yields the uniform margins. We then choose a copula shape from parametric or non-parametric families, or we leave this undetermined and run model selection based on the Bayesian information criterion (BIC). In the case of Caco2+, we can use the domain knowledge that permeability and TPSA are negatively correlated and specify a Clayton copula. Once a vine copula has been fit, it is fully described by the trees (structure) and bivariate (pair) copula densities associated with each edge. This is all we need to evaluate the CDF on the data points in Step 3. Note that the darker shaded points corresponding to higher CDF scores indeed approach the Pareto front.

We introduce a new perspective on MOBO based on the multivariate CDF. Our proposed MOBO acquisition function, BOtied, is computed by fitting a multivariate CDF on the surrogate predictions and extracting the ranks associated with the CDF scores. It is computationally attractive, as the CDF can be efficiently fit with vine copulas even when $M$ is large. Moreover, it enables model-based estimation of the Pareto front. When domain knowledge about the distribution of the objective values is available, it can be injected into the specification of the CDF model family. We also propose a new Pareto-compliant indicator for measuring the quality of approximate Pareto fronts, the CDF indicator. The CDF indicator, equipped with desirable properties such as invariance to monotonic transformations of the objectives, promises to complement the popular HV indicator.

Our method is general and lends itself to a number of immediate extensions. First, whereas we have implemented gradient-free optimization of BOtied in this work, we can take advantage of gradient-based optimization for improved efficiency. In conjunction with a differentiable sorting algorithm (e.g., Cuturi et al., 2019; Blondel et al., 2020), the computation of our acquisition function can be made differentiable for many parametric copula families. Second, we can consider constrained or discrete extensions for broader applicability. Finally, as many applications carry noise in the input as well as the function of interest, accounting for input noise through the established connection between copulas and multivariate value-at-risk (MVaR) estimation will be of great practical interest.

This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.