Paper deep dive
High-dimensional Multi-objective Bayesian Optimization with Learned Variable Interactions
Hongyan Wang, Jiayu Huang, Haotian Zheng, Xin Gao, Chi Ding, Ying Liu, Xia Wang, Qing Xu, Keqiang Li
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Multi-objective Bayesian optimization (MOBO) is effective in identifying the Pareto fronts for expensive black-box problems. However, most current MOBO approaches are limited to low-dimensional decision space due to its exponential sampling complexity. This paper presents decision variable interaction analysis-based MOBO, ViaMOBO, a generic framework for expensive multi-objective problems with high-dimensional decision space. The key idea of ViaMOBO is that it utilizes a variable interaction analysis model to determine whether the decision space can be completely or partially divided, and then performs local Bayesian optimization in the divided decision subspaces. Through the variable analysis model, it can be derived whether the objectives in black-box problems are separable, partially separable, or non-separable based on the potential independent or interdependent relationships among decision variables without any strong assumptions. We compare ViaMOBO with the state-of-the-art MOBO methods on both synthetic and real-world benchmarks. The experimental results demonstrate that ViaMOBO outperforms other related MOBO baselines in approximating the Pareto front of high-dimensional expensive multi-objective problems.
Tags
Links
- Source: https://arxiv.org/abs/2608.11713v1
- Canonical: https://arxiv.org/abs/2608.11713v1
Trouble viewing inline? Open PDF directly →
Full Text
67,244 characters extracted from source content.
Expand or collapse full text
High-dimensional Multi-objective Bayesian Optimization with Learned Variable Interactions 1st Hongyan Wang* Affiliation: School of Vehicle and Mobility Tsinghua University Beijing, China hongyan-wang@mail.tsinghua.edu.cn 2nd Jiayu Huang Affiliation: Faculty of Engineering The Hong Kong Polytechnic University Hong Kong, China jiayu1999.huang@connect.polyu.hk 3rd Haotian Zheng Affiliation: School of Vehicle and Mobility Tsinghua University Beijing, China zhenght21@mails.tsinghua.edu.cn 4th Xin Gao Affiliation: School of Vehicle and Mobility Tsinghua University Beijing, China gaoxin97@mail.tsinghua.edu.cn 5th Chi Ding Affiliation: School of Vehicle and Mobility Tsinghua University Beijing, China dc22@mails.tsinghua.edu.cn 6th Ying Liu Affiliation: School of Vehicle and Mobility Tsinghua University Beijing, China seuliuy@hotmail.com 7th Xia Wang Affiliation: School of Vehicle and Mobility Tsinghua University Beijing, China wang-xia@mail.tsinghua.edu.cn 8th Qing Xu Affiliation: School of Vehicle and Mobility Tsinghua University Beijing, China qingxu@tsinghua.edu.cn 9th Keqiang Li* Affiliation: School of Vehicle and Mobility Tsinghua University Beijing, China likq@tsinghua.edu.cn Abstract Multi-objective Bayesian optimization (MOBO) is effective in identifying the Pareto fronts for expensive black-box problems. However, most current MOBO approaches are limited to low-dimensional decision space due to its exponential sampling complexity. This paper presents decision variable interaction analysis-based MOBO, ViaMOBO, a generic framework for expensive multi-objective problems with high-dimensional decision space. The key idea of ViaMOBO is that it utilizes a variable interaction analysis model to determine whether the decision space can be completely or partially divided, and then performs local Bayesian optimization in the divided decision subspaces. Through the variable analysis model, it can be derived whether the objectives in black-box problems are separable, partially separable, or non-separable based on the potential independent or interdependent relationships among decision variables without any strong assumptions. We compare ViaMOBO with the state-of-the-art MOBO methods on both synthetic and real-world benchmarks. The experimental results demonstrate that ViaMOBO outperforms other related MOBO baselines in approximating the Pareto front of high-dimensional expensive multi-objective problems. Index Terms: Multi-objective Bayesian optimization (MOBO), Variable Interaction Analysis, High-dimensional Space, Expensive multi-objective Problems. I Introduction Multi-objective problems (MOPs) are often encountered in various experimental design problems, such as optimization in robotics [1], analog circuit design [2], aerospace engineering [3, 4], and hyper-parameter tuning in neural networks [5] etc. These problems are challenging to approximate the Pareto fronts (i.e., trade-offs among different objectives) since they are usually black-box problems with expensive evaluation costs, especially in high-dimensional space. Multi-objective Bayesian optimization (MOBO) ensures global convergence within limited computational cost, utilizing a probabilistic surrogate model (usually a Gaussian process surrogate) to approximate the original MOPs, and exploits an acquisition function to locate the next most promising points to be candidates for evaluation. Many researches have focused on improving the performance of MOBO on trading off the relationship among multiple objectives in black-box problems, mainly including EGO-based methods, hypervolume improvement (HVI)-based methods and predictive entropy-based methods, etc. EGO-based methods such as ParEGO [6], SMS-EGO [7] and MOEA/D-EGO [8], are extensions of EGO [9]. These EGO-based approaches ensures the global convergence for MOPs. HVI-based methods such as PEHVI [10], hypervolume scalarizariion approach [11], qEHVI [12], qNEHVI [13] mainly aim to improve the hypervolume calculation efficiency in solving MOPs. Entropy-based methods such as PESMO [14], PESMOC [15], MESMO [16], MESMOC [17], PFES [18] and MF-OSEMO [19] maximize the information gain of the Pareto optimum to select candidates for evaluation. Further MOBO studies mainly focus on uncertainty reduction in surrogate model PAL [20], improving the performance of scalarization [21, 11], solution diversity [22] and acquisition function [23]. Although MOBO approaches have made great advances in approximating the Pareto fronts of black-box multi-objective problems, most of them are limited to low-dimensional decision variable space due to the exponential sampling complexity with decision variable dimensions [24, 25]. The performances of these methods are compromised greatly over high-dimensional decision space. For high-dimensional cases, there are few research concerned MOBO, ReMO [26] and MORBO [27]. ReMO exploits random embedding [24, 25] to decompose the high-dimensional decision space into several subspaces with low dimensions to slove the black-box problems in low-dimensional space. MORBO uses local model with trust region [28] to perform local Bayesian optimization to search diverse solutions. However, the above methods manifest the following limitations. Firstly, ReMO is too hypothetical on the decision space. On one hand, it assumes that only few dimensions are effective while most dimensions do not affect the objectives. On the other hand, it assumes the same effective dimensions for all objectives in MOPs, which is obviously irrational for most optimization problems. Secondly, although without any strong assumption, MORBO aims at reducing the cubic time complexity with the number of observed data points, rather than emphasizing the high decision dimensionality and ignoring the sample efficiency of acquisition function. Nonparametric regression is challenging because of the curse of dimensionality [29], i.e., computational complexity depending exponentially in dimension. Besides, commonly used optimization heuristics for acquisition function also requires exponential complexity with dimension, leading to high cost and low optimization efficiency. In real-world high-dimensional expensive MOPs, the interaction relationship may exist among decision variables and could jointly influence the performance space. The fluctuation of one variable may cause the change of the other variables, and thus change the objective values. Take hyper-parameter tunning in neural networks as an example, the learning rate and dropout rate may influence the model performance simultaneously. However, current MOBO ignore considering this relationship, which results in severely decreased performance. According to whether there’s interaction relationship among decision variables, high-dimensional MOPs can be categoried into separable and non-separable problems. For separable problems, it is desirable to partition the decision space into disjoint subspaces by learning variable interaction relationship to reduce the sampling complexity and improve optimization efficiency, especially for high-dimensional problems. However, how to learn the interaction is a great challenge. For expensive high-dimensional separable MOPs, we present a MOBO framework, ViaMOBO, based on decision Variable Interaction Learning to alleviate the curse of dimensionality. ViaMOBO first determines whether a MOP is separable by learning the interaction among decision variables, and then divides the decision space into disjoint subspaces, with each subspace containing interdependent variables. To avoid real expensive function evaluation, ViaMOBO exploits a binary classifier to predict the magnitude relationship between two objective values of two decision variable vectors. In this way, ViaMOBO perform optimization for MOPs in low-dimensional decision sub-spaces, ensuring sample efficiency of acquisition function and avoiding the curse of dimensionality. Then, ViaMOBO selects a batch of candidates incorporating virtual derivate observations to avoid over-exploration of the search space as much as possible. Finally, trust region is used to avoid exponential complexity of observed points. Our contributions are as follows: • A general MOBO framework, ViaMOBO, is presented, in which an additive kernel structure is first introduced into multi-objective settings, and the sample efficiency of the acquisition function in MOBO is considered. • A decision variable interaction learning strategy is developed with the following properties: (a) the separability of high-dimensional MOPs is determined by learning the interaction relationships among decision variables; (b) expensive function evaluations are avoided through the use of a binary classifier for interaction learning. • The effectiveness of ViaMOBO is verified by comparisons with state-of-the-art MOBO methods on both synthetic and the real-world expensive problems. The experimental results demonstrate that ViaMOBO significantly outperforms the baselines. I Background I-A Preliminaries I-A1 Bayesian Optimization (BO) BO method uses a Gaussian process (GP) surrogate to model the objective, meanwhile using an acquisition function to recommend candidate solutions. Given a single-objective problem: ∗=argmin∈Xf(),X⊆ℝDx =arg _x∈ Xf(x),X ^D, f()f(x) is approximated with GP, i.e., f()∼(m(),k(,′))f(x) (m(x),k(x,x )), where m()m(x), k(,′)k(x,x ) are the mean and kernel function, respectively. Given the noisy observations 1:t=(xi,yi)∣yi=f(xi)+ϵi,i=1,…,tD_1:t=\(x_i,y_i)\ y_i=f(x_i)+ _i,i=1,…,t\, where ϵ∼(0,σnoise2)ε (0, _noise^2), the predictive model is: (t+1∣1:t,t+1)=(μt(t+1),σt2(t+1))P(y_t+1 _1:t,x_t+1)=N( _t(x_t+1), _t^2(x_t+1)). μt(t+1)=T(K+σnoise2I)−11:t _t(x_t+1)=k^T(K+ _noise^2I)^-1y_1:t and σt2(t+1)=(t+1,t+1)−T(K+σnoise2I)−1 _t^2(x_t+1)=k(x_t+1,x_t+1)-k^T(K+ _noise^2I)^-1k are the mean and variance function, respectively, and =[(t+1,1),(t+1,2),…,(t+1,t)]k=[k(x_t+1,x_1),k(x_t+1,x_2),…,k(x_t+1,x_t)]. Then, the acquisition function balancing exploitation (minimization of μt() _t(x)) and exploration (maximization of σt() _t(x)) is built to search for the most promising candidates. I-A2 Separable and Non-Separable Problems Definition I.1. Two decision variables xix_i and xjx_j are interacting (interdependent variables), if there exist x,a1,a2,b1x,a_1,a_2,b_1 and b2b_2 meeting [30, 31] f(x)|xi=a2,xj=b1<f(x)|xi=a1,xj=b1∧ f(x)|_x_i=a_2,x_j=b_1<f(x)|_x_i=a_1,x_j=b_1 (1) f(x)|xi=a2,xj=b2<f(x)|xi=a1,xj=b2, f(x)|_x_i=a_2,x_j=b_2<f(x)|_x_i=a_1,x_j=b_2, where f(x)|xi=a2,xj=b1≜f(x1,…,xi−1,a2,…,xj−1,b1,…,xD). f(x)|_x_i=a_2,x_j=b_1 f(x_1,…,x_i-1,a_2,…,x_j-1,b_1,…,x_D). (2) Definition I.2. A function f(x):ℝn→ℝf(x):R^n is partially separable if it is a sum of sub-functions f(x)=∑i=1Mfi(xi)f(x)= _i=1^Mf_i(x_i), and each one depends on a group of interdependent decision variables. The function f is called completely additively separable or fully separable function if all sub-functions are 1-D. f(x)f(x) is non-separable if it is neither partially separable nor fully separable. [32] Definition I.3. A MOP is a separable multi-objective problem if every objective is fully separable. Otherwise, F(x)F(x) is a non-separable MOP. [33] For convenience, hereafter, we refer to both partially and fully separable function as separable functions in this paper, since under both cases, the decision variables can be decomposed to groups. The global optimum of a separable problem can be obtained by considering each dimension (fully separable) or each group of dimensions (partially separable) of variable in turn, independently of one another [33]. For a MOP, this means that we can consider one dimension (fully separable) or a group of dimensions (partially separable) at one time to find at least some points on the Pareto front, which makes it easier for optimizer to approximate the Pareto front than non-separable MOPs. I-B High-dimensional Bayesian Optimization Various assumptions on the decision and objective/performance space have been made to extend BO to high (≥10≥ 10) dimensions [34]. Random embedding-based methods, REMBO [35], SI-BO [36], ALEBO [37], PCA-BO [38] assume that there are important and unimportant dimensions, and only consider the important dimensions of variables. Instead of directly optimizing in the original high-dimensional decision space, MS-UCB [39] performs acquisition function maximization over discretized low-dimensional subspaces embedded into the full space. By Assuming an additive structure for the function, ADD-GP-UCB [24] reduces the computation exponential complexity in dimension when using global optimization heuristics to maximize the acquisition function. P-GP-UCB [40]generalizes the additive assumption and the low-dimensional assumption to a projected-additive assumption to handle a broader class of functions. Batch-ADD-GP-UCB [25] extends ADD-GP-UCB to sample a batch of candidate solutions. G-Add-GP-UCB [41] extends GP-UCB by removing the restriction that the variable subsets must be mutually disjoint, thereby allowing additive structures with overlapping subsets. TuRBO [28] uses trust regions to construct a group of local models and searches global optimum across these local models. However, most of these methods make too strong hypotheses about the decision or objective spaces, and they are confined to high-dimensional single-objective problems. MORBO [27] improves TuRBO not only by selecting trust region center in a coordinate fashion to improve convergence but also chooses new candidates by collaboratively optimized a shared hypervolume-based global utility. However, the trust region-based TuRBO and MORBO targets to address the challenge in fitting a global surrogate model when the function is heterogeneous, while failing to overcome the problem that the search space increases exponentially with dimension, which make global solver for acquisition function. To address this, single-objective random embedding with high dimensions has been scaled to high-dimensional MOPs that assumes only few dimensions influence the objective values [26]. However, similar to other random embedding methods, ReMO assumes too ideal structure of the real-world black-box problems. Although other related approaches [42, 43] effectively address high-dimensional problems, the partitioning of the decision space involves elements of randomness and subjectivity. I The proposed ViaMOBO Algorithm 1 Framework of ViaMOBO 1: Input: An MOP F(x):ℝD→ℝMF(x):R^D ^M, maximum function evaluations maxEvalsmaxEvals, sample budget T 2: Output: Approximate Pareto front PF∗PF^* 3: Draw initial Sobol points XinitX_init and evaluate initial points XinitX_init with original MOP F(x)F(x) ⊳ Initialization 4: Generate training and testing data given D0D_0 5: for m=1,…,Mm=1,…,M do 6: train SVM model hi(⋅)h_i(·) for the iith objective fi(x)f_i(x) 7: Set x defined in Eqn.(3) as the decision vector that maximizes the iith objective fi(x)f_i(x) 8: Perturb x at the iith and jjth dimension with random value to obtain xi′,xj′x_i ,x_j and xij′x_ij , where i,j∈[1,D]∧i≠ji,j∈[1,D] i≠ j 9: Predict the magnitude relationship of the objective values of xi,xi′,xj′x_i,x_i ,x_j and xij′x_ij with the trained SVM model hi(⋅)h_i(·) 10: if [f(xi′),f(xj′)][f(x_i ),f(x_j )] dominates or is dominated by [f(x),f(xij′)][f(x),f(x_ij )] then 11: the iith decision variable interacts with the jjth decision variable 12: else 13: the iith decision variable does not interact with the jjth decision variable 14: end if 15: Learn all the interacting decision variables Ωi _i of fi(x)f_i(x) based on transitivity defined in Def.I.2 16: end for ⊳ Variable interaction learning 17: Learn the interdependent decision variables D=∪jNΩjD= _j^N _j of F(x)F(x) based on Def.I.3 to judge whether F(x)F(x) is separable 18: if F(x)F(x) is non-separable then 19: optimize F(x)F(x) with other high-dimensional MOBO approach (e.g., MORBO [27]) 20: else 21: while evals≤MaxEvalsevals do 22: Fit a local model for each objective fi(x)f_i(x), and obtain the additive structure of F(x)F(x) ⊳ Additive kernel structure 23: Select q candidates using the defined acquisition function (e.g., UCB, EI) ⊳ Batch candidate sampling 24: Evaluate candidates on F(x)F(x) and obtain new observations 25: end while 26: end if This section introduces ViaMOBO, a multi-objective BO method with decision variable interaction based on theoretical analysis. As shown in Algorithm 1, to improve the sampling efficiency of acquisition function, ViaMOBO learns whether a MOP is separable by analyzing the variable interaction or not. Unlike random embedding-based methods [24, 25], ViaMOBO makes no assumption on the structure of black-box problems, but deduces whether it is separable given prior observations or not. If the MOP is non-separable, we recommend to use other high-dimensional MOBO method, and we strongly recommend MORBO [27] here due to its high efficiency. For a separable MOP, ViaMOBO inferences the additive GPs and fits a local model of each objective in MOP in each lower dimensional decision sub-space. Then, multiple candidates are recommended by optimizing the acquisition function, and then evaluated by the true objective. I-A Re-definition of Separable MOP Def.I.2 gives the partially, fully separable and non-separable functions in single-objective cases. Hereafter, we refer to both partially and fully separable function as separable functions, since under both cases, the decision variables can be decomposed to groups. Given this, we redefine Def.I.3 for multi-objective cases as following. Definition I.1. A MOP is a separable multi-objective problem if every objective is separable (partially or fully separable). Otherwise, F(x)F(x) is a non-separable MOP if all the objective are non-separable. Definition I.2. If two decision variables xix_i and xjx_j are interacting, and xkx_k interacts with one of the two variables, then all these three variable are interacting. Above definition indicates that there’s transitivity among this characteristic and it can be extended to a case of infinite number of decision variable dimensions. This definition is necessary for MOPs due to the possible conflicting objectives. Definition I.3. Given a separable MOP F(x)=(f1(x),…,fM(x)):ℝD→ℝMF(x)=(f_1(x),…,f_M(x)):R^D ^M, and the sub-functions of fi(x)=∑j=1Mifj(x),x∈Ωif_i(x)= _j=1^M_if_j(x),x∈ _i where Ωi _i is the interacting variable partition of fi(x)f_i(x), then the interacting variable partition of F(x)F(x) is ∪i=1MΩi _i=1^M _i. I-B Variable Interaction Learning Based on Def.I.1, Def.I.1 and Def.I.2, we can find that if the interacting variables are obtained, then all the decision dimensions can be decomposed and then whether the objective is separable can be determined or not. Therefore, the key point is the variable interaction learning. For convenience, we re-write the decision vectors as below. x=(x1,…,xi−1,xi=a1,…,xj−1,xj=b1,…,xD) x=(x_1,…,x_i-1,x_i=a_1,…,x_j-1,x_j=b_1,…,x_D) (3) xi′=(x1,…,xi−1,xi=a2,…,xj−1,xj=b1,…,xD) x_i =(x_1,…,x_i-1,x_i=a_2,…,x_j-1,x_j=b_1,…,x_D) xj′=(x1,…,xi−1,xi=a1,…,xj−1,xj=b2,…,xD) x_j =(x_1,…,x_i-1,x_i=a_1,…,x_j-1,x_j=b_2,…,x_D) xij′=(x1,…,xi−1,xi=a2,…,xj−1,xj=b2,…,xD) x_ij =(x_1,…,x_i-1,x_i=a_2,…,x_j-1,x_j=b_2,…,x_D) Re-examining Def.I.1, it indicates that if there exists a decision vector x whose iith and jjth variable (a1a_1 and b1b_1) can be substituted with values a2a_2 and b2b_2 to generate a strong dominance relationship between [f(xi′),f(xj′)][f(x_i ),f(x_j )] and [f(x),f(xij′)][f(x),f(x_ij )] is established, then xix_i and xjx_j are interacting. However, the function evaluations of many real-world MOPs are so expensive that it is unacceptable to cost function evaluations on additional objective value calculation f(xi′),f(xj′),f(xij′)f(x_i ),f(x_j ),f(x_ij ). To avoid the expensive function calculation of xi′,xj′andxij′x_i ,x_j and\ x_ij , we use a binary classifier (e.g., SVM) to learn the magnitude relationship of function values instead of evaluating xi′,xj′andxij′x_i ,x_j and\ x_ij by using the true objective function. With initialized data points Xinit,F(xinit)|F(xinit)=(f1(x),f2(x),…,fM(x))\X_init,F(x_init)|F(x_init)=(f_1(x),f_2(x),…,f_M(x))\, we generate training and test dataset of each objective fi(x)f_i(x) in F(x)F(x) by concatenating two decision variables vectors x1,x2x_1,x_2 into one vector, denoted as x1⊕x2x_1 x_2. If fi(x1)<fi(x2)f_i(x_1)<f_i(x_2), then the label is set to positive while negative if otherwise. To make the training and test dataset balanced, we regard x1⊕x2x_1 x_2 and x2⊕x1x_2 x_1 as two different data. With the labeled samples, a binary classifier (e.g., supported Vector Machine (SVM)) is trained to learn the relationship of the objective values of two decision variable vectors. The training and test data ratio is 70%70\% and 30%30\%. In this paper, we adopt SVM because it can maintain satisfactory classification performance in sparsely sampled regions while offering both parametric and non-parametric flexibility. Specifically, to determine whether two dimensions in decision space of a MOP F(x):ℝD→ℝMF(x):R^D ^M are interacting or not, we first use Def.I.1 on the observed decision variable x that maximizes the iith objective fi(x)f_i(x) to separately learn the interacting variables of all M objectives in MOPs. When learning, we perturb x at the iith and jjth dimension with random value in problem bounds to obtain xi′,xj′,andxij′x_i ,x_j ,and\ x_ij . Then the trained SVM predicts the relationship f(xi′),f(xj′)f(x_i ),f_(x_j ) and f(xij′)f_(x_ij ) to compare whether there’s a dominance relationship between [f(xi′),f(xj′)][f(x_i ),f(x_j )] and [f(x),f(xij′)][f(x),f(x_ij )]. Then, Def.I.2 is used to learning all the possible interacting variables. Finally, we use Def.I.3 to deduce the interdependent variable partition D=∪jNΩjD= _j^N _j of F(x)F(x). I-C Additive Kernel for MOPs If the MOP F(x):ℝD→ℝMF(x):R^D ^M is separable, then F(x)=(∑i=1Nf1i(x),∑i=1Nf2i(x),…,∑i=1NfMi(x)), F(x)=( _i=1^Nf_1^i(x), _i=1^Nf_2^i(x),…, _i=1^Nf_M^i(x)), (4) where D=∪jNΩjD= _j^N _j and Ωi∩Ωj=∅,i,j∈[1,…,N],i≠j _i∩ _j= ,i,j∈[1,…,N],i≠ j. Similar to ADD-GP-UCB [24], we call Ωi _i disjoint sub-groups. For each sub-objective fi(x)f_i(x) in F(x)F(x), fi(x)f_i(x) is single-objective and additive, fi(x)=f(1)(x(1))+f(2)(x(2))+⋯+f(N)(x(N)), f_i(x)=f^(1)(x^(1))+f^(2)(x^(2))+…+f^(N)(x^(N)), (5) where x(j)∈Ωix^(j)∈ _i. Following single-objective ADD-GP-UCB [24], we assume f(j)∼GP(μi(j)(x),kj(j)(x(i),x(j)′))f^(j) GP( _i^(j)(x),k_j^(j)(x^(i),x^(j) )), then fi(x)∼GP(μi(x),ki(x,x′))f_i(x) GP( _i(x),k_i(x,x )) in noiseless case. The mean μi(x) _i(x) and kernel function ki(x,x′)k_i(x,x ) are μi(x) _i(x) =μ(1)(x(1))+⋯+μ(N)(x(N)) =μ^(1)(x^(1))+…+μ^(N)(x^(N)) (6) ki(x,x′) k_i(x,x ) =k(1)(x(1),x(1)′)+⋯+k(M)(x(M),x(M)′). =k^(1)(x^(1),x^(1) )+…+k^(M)(x^(M),x^(M) ). The sub-kernel k(j)(x(j),x(j)′)k^(j)(x^(j),x^(j) ) is defined in the decision sub-space Ωi _i. Then, we can inference the posterior distribution of the component ft+1(j)(x)f_t+1^(j)(x) given the observed data (Xt,Ft(x))|X=x1,…,xD,Ft(x)=(f1t(x),f2t(x),…,fMt(x))\(X_t,F_t(x))|X=\x_1,…,x_D\,F_t(x)=(f_1^t(x),f_2^t(x),…,f_M^t(x))\, is defined as the following. ft+1(j)(x)∼GP(μt+1(j)(x(j)),kt+1(j)(x(j),x(j)′)), f_t+1^(j)(x) GP(μ^(j)_t+1(x^(j)),k^(j)_t+1(x^(j),x^(j) )), (7) μt+1(j)(x(j))=k(j)(xt+1(j),X(j)k(X,X)−1), μ^(j)_t+1(x^(j))=k^(j)(x_t+1^(j),X^(j)k(X,X)^-1)f_i, kt+1(j)(x(j),x(j)′)=k(j)(xt+1(j),xt+1(j)′) k^(j)_t+1(x^(j),x^(j) )=k^(j)(x_t+1^(j),x_t+1^(j) ) −k(j)(xt+1(j),X(j))k(X,X)−1k(j)(X,x(j)) -k^(j)(x_t+1^(j),X^(j))k(X,X)^-1k^(j)(X,x^(j)) and thus ft+1(x)∼GP(μt+1(x),kt+1(x,x′))f_t+1(x) GP( _t+1(x),k_t+1(x,x )). More details about single-objective additive kernel can be found in ADD-GP-UCB [24]. In this way, given the observed data (Xt,Ft(x))|X=x1,…,xD,Ft(x)=(f1t(x),f2t(x),…,fMt(x))\(X_t,F_t(x))|X=\x_1,…,x_D\,F_t(x)=(f_1^t(x),f_2^t(x),…,f_M^t(x))\, we can obtain the posterior of F(x)F(x) as the following. F(x)∼(GP(μt+11(x),k1t+1(x,x′)), F(x) (GP(μ^t+1_1(x),k_1^t+1(x,x )), (8) GP(μ2t+1(x),k2t+1(x,x′)),…,GP(μMt+1(x),kMt+1(x,x′)). GP(μ^t+1_2(x),k_2^t+1(x,x )),…,GP(μ^t+1_M(x),k_M^t+1(x,x )). I-D Batch candidate sampling As discussed in I-A2, for a separable MOP, we can consider one dimension (fully separable) or a group of dimensions (partially separable) a time to find at least some points on the Pareto front. Therefore, after we learn the additive structure and inference the additive structure of posterior distribution of MOP, we can acquire candidate points by solving the acquisition functions in lower decision sub-spaces Ωi _i, such as EI [9], UCB [24] and EHVI [13]. Sampling with UCB. ADD-GP-UCB shows an additive UCB acquisition function based on the additive kernel, i.e., φt(x)=μt−1(x)+βt1/2∑j=1Nσt−1(j)(x(j)) _t(x)= _t-1(x)+ _t^1/2 _j=1^N _t-1^(j)(x^(j)), where μ(x)μ(x) and σ2(x)σ^2(x) (k(x,x′)k(x,x )) are defined in Eqn (6). The additive acquisition function can be maximized by maximizing each component in decision subspace Ωi _i, which is much easier than optimizing it in D-dimensional decision space. Based on single-objective ADD-GP-UCB, we can deduce the multi-objective UCB as the following. At(x)=(φ1t(x),φ2t(x),…,φ1M(x)), A_t(x)=( _1^t(x), _2^t(x),…, _1^M(x)), (9) where each φit(x),i∈[1,…,M] _i^t(x),i∈[1,…,M] is an additive UCB. To obtain a batch of candidates, we use qqUCB strategy [12] to scalarize all the M additive UCBs. Fig. 1: (Left) Illustration of all the related methods on synthetic problem with relative high 10D10D decision dimensions. (Medium) Illustration of all the related methods on synthetic problem with relative high 30D30D decision dimensions. (Right) Illustration of all the related methods on synthetic problem with higher 100D100D decision dimensions. Sampling with EI and EHVI.Although EI and EHVI can not be directly deduced as additive EI in multi-objective cases even based on our additive structure, we can acquire it in |Ωi|| _i|-dimensional decision space, which obviously improves the optimization efficiency for global optimization heuristics [24, 25]. We use qqEI and EHVI [13] (qqEHVI is out of memory in our experiments) as our sampling strategy in this paper. To draw large batch sizes q of candidates, we borrow the idea from MORBO[27] that uses Thompson sampling to obtain q posteriors from GP, and optimizes the acquisition function group by group. However, it is undesirable to sequentially maximize the first group and then the next group. Instead, we select the kkth candidate by maximizing the kkth group Ωk _k of decision space and all other groups we use values from the Pareto optimums from the previous observed points. After N iterations we cycle through the groups if q>Nq>N. N is the number of sub-groups of the decision space. To avoid over-exploration (boundary issue) [44] for high-dimensional BO methods when optimizing acquisition, we use the virtual derivative sign observations [45]. IV Experiments and Results IV-A Experimental Settings To validate the effectiveness of ViaMOBO, it is compared with NSGA-I11 1 https://github.com/anyoptimization/pymoo[46], ParEGO†footnotemark: [6], MOEA/D-EGO†footnotemark: [8], TSEMO22 2 https://github.com/yunshengtian/DGEMO[47], USeMO-EI†footnotemark: [23], DGEMO†footnotemark: [22], qqParEGO [12], qqqLogNEHVI[48], MORBO [27], and random search with Sobol—a quasi-random baseline [13].q-series methods, MORBO and ViaMOBO experiments are conducted using BoTorch 33 3 https://botorch.org/ [49]. Among the EHVI-family methods, we select qqLogNEHVI as a representative baseline due to its numerically stable logarithmic reformulation of the acquisition function. All methods are run for 10 replications, and the initial data points used to build GP are quasi-random for each replication. We use the same hypar-prameters for all the related methods. Hypervolume (HV) is used to evaluate the related performances. The acquisition functions used in this paper are EI [9], UCB [24] and EHVI [13]. All of them are Monte-Carlo acquisition functions [49]. The benchmarks used in this paper include synthetic functions and real-world benchmarks. Specifically, we first measure the performances of all methods on the synthetic function, separable DTLZ2 with D=10D=10, D=30D=30 and D≥100D≥ 100 with 33 objectives. Then, all the methods are evaluated on three real-world problems, including two aerodynamic shape-optimization[50] problems and trajectory planning[27]. Specifically, trajectory planning is with D=60D=60 and exhibits strong sequential coupling because neighboring control points jointly determine the resulting trajectory. The airfoil design benchmark optimizes 2020 and 4040 shape parameters with three objectives—drag, lift, and geometric regularity—using XFOIL under fixed aerodynamic conditions. The reference point for HV calculation of is min=[6,6,6]r_min=[6,6,6] and min=(0.25, 0.0, 0.5)r_ =(0.25,\;0.0,\;0.5) for 33-objective DTLZ2, airfoil design problem, respectively. All the related methods are initialized 200200 with Sobol points to construct the initial GP model. For all the synthetic problems, we use 20002000 maximum function evaluations, and 50 batch size for all the related methods. We find that these two methods even can not run when there is three objective. Therefore, we set the batch size to 55 for these two methods in this paper. For ViaMOBO, we use 70%70\% and 30%30\% of the initial 200200 points to generate our training and test datasets. We set the perturbation of decision values to 100 in the normalized standard bounds [0,1][0,1] of all the related problems to find if there is a2,b2a_2,b_2 and to determine whether two variables are interacting. We use a constant-mean independent GP with an RBF kernel is utilized, and its hyperparameters are estimated by maximizing the marginal log-likelihood. This GP model is consistently applied to all BO baseline methods. For all the related methods, we use Monte-Carlo (MC) acquisition functions [49], including EI, UCB and EHVI. Following MORBO [27], we use the same number of quasi-MC samples and Fourier basis functions. For qqParEGO and qqLogNEHVI, the acquisition function is optimized via L-BFGS-B using 20 random restarts. For MORBO and ViaMOBO, we use 4096 discrete points to optimize the related acquisition function for all problems. Details can be found in MORBO [27]. Besides, we use virtual deriviate information to overcome over-exploration in high-dimensional cases. TABLE I: Hypervolume performance and computational cost on the DTLZ2 problems with 10, 30, and 100 decision variables. The hypervolume results are reported as mean ± standard deviation. HV@500 and HV@1000 denote the hypervolume values obtained after 500 and 1,000 function evaluations, respectively. AUC-HV denotes the normalized area under the hypervolume convergence curve. Time per batch is calculated over 36 optimization batches with a batch size of q=50q=50. Relative runtime is normalized by the runtime of ViaMOBO at the corresponding dimensionality. All experiments are conducted on a server equipped with two Intel Xeon Platinum 8470Q CPUs (104 physical cores and 208 logical threads), 754 GiB of RAM, and an NVIDIA GeForce RTX 5090 GPU with 32 GB of memory. Dim. Algorithm HV@500 HV@1000 Final HV AUC-HV Runtime per seed Time per batch Relative runtime 10 Random 213.975 214.409 214.646±0.033214.646± 0.033 213.80 1.87 s – <0.001×<0.001× 10 NSGA-I 214.652 215.181 215.347±0.014215.347± 0.014 214.45 3.63 s – <0.001×<0.001× 10 ParEGO 214.777 214.839 214.851±0.065214.851± 0.065 214.59 17.92±2.0317.92± 2.03 h 1792 s 11.1×11.1× 10 MOEA/D-EGO 215.015 215.065 215.184±0.113215.184± 0.113 214.79 5.27±0.655.27± 0.65 h 527 s 3.3×3.3× 10 TSEMO 215.283 215.321 215.342±0.008215.342± 0.008 215.04 1.68±0.471.68± 0.47 h 168 s 1.0×1.0× 10 USeMO-EI 214.489 215.259 215.373±0.020215.373± 0.020 214.83 3.56±0.623.56± 0.62 h 356 s 2.2×2.2× 10 qParEGO 215.419 215.438 215.446±0.000215.446± 0.000 215.40 5.49±0.555.49± 0.55 h 549 s 3.4×3.4× 10 qLogNEHVI – – – – >48>48 h (timeout) – – 10 DGEMO 215.428215.428 215.445215.445 215.454±0.001215.454± 0.001 215.16215.16 9.13±0.559.13± 0.55 h 913 s 5.6×5.6× 10 MORBO 215.077 215.146 215.174±0.029215.174± 0.029 215.08 14.46±3.2514.46± 3.25 h 1446 s 8.9×8.9× 10 ViaMOBO 215.019 215.321 215.402±0.020215.402± 0.020 214.61 1.62±0.371.62± 0.37 h 162 s 1.0×1.0× 30 Random 202.885 204.536 205.837±0.207205.837± 0.207 202.68 1.95 s – <0.001×<0.001× 30 NSGA-I 207.729 212.909 214.674±0.217214.674± 0.217 209.31 3.45 s – <0.001×<0.001× 30 ParEGO 208.153 208.168 208.232±0.323208.232± 0.323 206.47 17.25±0.9917.25± 0.99 h 1725 s 10.0×10.0× 30 MOEA/D-EGO 203.874 210.534 212.524±0.654212.524± 0.654 207.66 5.57±0.135.57± 0.13 h 557 s 3.2×3.2× 30 TSEMO 205.112 207.070 210.108±0.455210.108± 0.455 206.22 3.69±0.433.69± 0.43 h 369 s 2.1×2.1× 30 USeMO-EI 203.997 206.015 208.959±0.957208.959± 0.957 205.29 4.77±0.134.77± 0.13 h 477 s 2.8×2.8× 30 qParEGO 212.747 214.760 215.396±0.023215.396± 0.023 213.95 6.42±0.266.42± 0.26 h 642 s 3.7×3.7× 30 qLogNEHVI – – – – >48>48 h (timeout) – – 30 DGEMO 214.659214.659 215.398215.398 215.439±0.003215.439± 0.003 212.52212.52 6.84±0.476.84± 0.47 h 684 s 4.0×4.0× 30 MORBO 212.052 213.511 213.704±0.235213.704± 0.235 212.05 5.42±0.665.42± 0.66 h 542 s 3.1×3.1× 30 ViaMOBO 212.253 214.571 215.184±0.217215.184± 0.217 210.97 1.73±0.621.73± 0.62 h 173 s 1.0×1.0× 100 Random 9.785 13.524 18.790±2.74118.790± 2.741 12.72 3.43 s – <0.001×<0.001× 100 NSGA-I 32.829 98.413 180.236±4.822180.236± 4.822 95.18 0.59 s – <0.001×<0.001× 100 ParEGO 7.167 7.167 37.925±21.96537.925± 21.965 13.97 >48>48 h (timeout) – >12.0×>12.0× 100 MOEA/D-EGO 10.786 12.248 58.452±8.14258.452± 8.142 21.66 18.38±1.0718.38± 1.07 h 1838 s 4.6×4.6× 100 TSEMO 30.010 58.495 94.433±4.48794.433± 4.487 55.22 26.74±1.0926.74± 1.09 h 2674 s 6.7×6.7× 100 USeMO-EI 7.230 8.334 15.660±7.13115.660± 7.131 9.33 25.10±0.3525.10± 0.35 h 2510 s 6.3×6.3× 100 qParEGO 6.806 20.123 143.954±36.732143.954± 36.732 50.80 8.69±0.878.69± 0.87 h 869 s 2.2×2.2× 100 qLogNEHVI 6.584 6.584 74.132±5.91674.132± 5.916 18.29 9.85±0.659.85± 0.65 h 985 s 2.5×2.5× 100 DGEMO 6.384 152.107 214.679±0.114214.679± 0.114 115.82 41.90±0.4041.90± 0.40 h 4190 s 10.5×10.5× 100 MORBO 51.898 151.351 202.670±2.099202.670± 2.099 127.110 2.836±1.3412.836± 1.341 h 283.6 s 0.71×0.71× 100 ViaMOBO 108.291108.291 186.882186.882 210.205±3.789210.205± 3.789 148.91148.91 3.99±2.163.99± 2.16 h 399 s 1.0×1.0× IV-B Results on Synthetic Benchmark Problems We first consider a separable synthetic problem, DTLZ2, a widely used test problem to study the performance of multi-objective methods. Strictly speaking, DTLZ2 is not separable, since optimizing one decision variable at a time will not identify all the global optimum. However, there are multiple global optimum and it is many-to-one, if optimizing along each coordinate direction is sufficient to attain a global optimum [33]. Therefore, it is usually regarded as a separable MOP. Figure 1 presents the HV convergence curves on the 10-, 30-, and 100-dimensional DTLZ2 problems. Table I shows HV@@500, HV@1000, final HV and their computational costs on DTLZ2. The first two metrics characterize early- and intermediate-stage convergence, while the final HV represents the performance after 2,000 function evaluations44 4 qLogNEHVI incurred substantial computational overhead due to Monte Carlo estimation, Pareto-set processing, and hypervolume partitioning for 10D and 3D 3-objective DTLZ2. Its runtime is not determined solely by the decision-space dimensionality; the larger non-dominated sets encountered in the low dimensional cases increased the partitioning cost. Under the given limit, only one 10D run completed successfully (HV@@500, HV@1000, AUC-HV and final HV are 215.272, 215.272, 215.299, 215.16 using 35.46h, respectively). Therefore, qLogNEHVI is reported only for the 100D problem.. Fig. 2: (Left) The average runtime per seed on a logarithmic scale from 10D to 100D. (Right) AUC-HV normalized by the best result at each dimensionality. TABLE I: Computational cost on the 20- and 40-dimensional Airfoil optimization problems. Runtime is reported as mean ± standard deviation. Time per batch is calculated using 36 optimization batches with q=50q=50. Relative runtime is normalized by the runtime of ViaMOBO at the corresponding dimensionality. 20D 40D Algorithm Runtime/seed Time/batch Relative Runtime/seed Time/batch Relative Random 0.077±0.0030.077± 0.003h – 0.097×0.097× 0.110±0.0090.110± 0.009h – 0.081×0.081× NSGA-I 0.082±0.0030.082± 0.003h – 0.103×0.103× 0.104±0.0040.104± 0.004h – 0.077×0.077× ParEGO 18.193±1.69118.193± 1.691h 1819.3s 22.76×22.76× >48>48h (timeout) – >35.47×>35.47× qParEGO 6.704±0.4136.704± 0.413h 670.4s 8.39×8.39× 8.532±0.4358.532± 0.435h 853.2s 6.31×6.31× MOEA/D-EGO 3.476±0.3963.476± 0.396h 347.6s 4.35×4.35× 6.394±1.0396.394± 1.039h 639.4s 4.72×4.72× TSEMO 3.760±0.4283.760± 0.428h 376.0s 4.70×4.70× 8.156±1.2868.156± 1.286h 815.6s 6.03×6.03× USeMO-EI 3.592±0.3823.592± 0.382h 359.2s 4.49×4.49× 7.867±0.6187.867± 0.618h 786.7s 5.81×5.81× DGEMO 9.045±1.6919.045± 1.691h 904.5s 11.32×11.32× 18.328±4.59018.328± 4.590h 1832.8s 13.54×13.54× MORBO 8.985±3.9738.985± 3.973h 898.5s 11.24×11.24× 9.371±5.1339.371± 5.133h 937.1s 6.93×6.93× ViaMOBO 0.799±0.1190.799± 0.119h 79.979.9s 1.00×1.00× 1.353±0.1741.353± 0.174h 135.3135.3s 1.00×1.00× From Figure 1 and Table I, we observe that on the 10-dimensional problem, DGEMO achieves the highest final HV, followed closely by qParEGO and ViaMOBO. The difference between DGEMO and ViaMOBO is only 0.024%0.024\%, whereas DGEMO and qParEGO require approximately 5.6×5.6× and 3.4×3.4× the runtime of ViaMOBO, respectively. Therefore, ViaMOBO provides a near-optimal Pareto-front approximation at substantially lower computational cost. On the 30-dimensional problem, DGEMO obtains the highest HV@500, HV@1000, final HV, and AUC-HV. qParEGO also achieves a competitive final HV only approximately 0.020%0.020\% below DGEMO, but requires approximately 3.7×3.7× the runtime of ViaMOBO. ViaMOBO is only 0.118%0.118\% below DGEMO, while requiring approximately one-quarter of its runtime. Other baselines require 2.12.1–10.0×10.0× the runtime of ViaMOBO while producing lower final HV. On the 100-dimensional problem, ViaMOBO exhibits a clear early- and intermediate-stage convergence advantage, achieving the highest HV@500 and HV@1000 values. It also obtains the highest AUC-HV, which is approximately 28.6%28.6\% higher than DGEMO’s. Although DGEMO ultimately achieves a final HV approximately 2.13%2.13\% higher than ViaMOBO, it requires approximately 10.5×10.5× the runtime. MORBO’s final HV and AUC-HV are approximately 3.58%3.58\% and 14.64%14.64\% lower than those of ViaMOBO, respectively. Therefore, although MORBO remains competitive in terms of final solution quality, ViaMOBO provides faster convergence and better overall performance under a limited evaluation budget. The observed MORBO runtime is 2.836±1.3412.836± 1.341 hours per seed. qParEGO requires approximately 2.2×2.2× the runtime of ViaMOBO but obtains a substantially lower final HV, indicating reduced effectiveness as the dimensionality increases. Figure 2 further demonstrates that ViaMOBO exhibits a relatively moderate increase in computational cost while maintaining competitive optimization quality as the dimensionality increases. Overall, ViaMOBO achieves near-optimal performance at low dimensionality and demonstrates an increasingly advantageous balance among convergence speed, final solution quality, and computational cost in high-dimensional decision spaces. IV-C Results on Real-world Problems This section considers real-world problems, including 20- and 40-dimensional airfoil benchmark optimization problem and a bi-objective trajectory planning problem with 6060 decision variables [27]. We compare the related methods using 2000 total function evaluations with batch size q=50q=50 at each iteration. Fig. 3: (Left) The results of log HV difference and (Right) HV–runtime of all methods on airfoil problem with 2020 decision variables. IV-C1 20D Airfoil shape-optimization problem As shown in Figure 3, TSEMO exhibits fast early-stage convergence on the 20-dimensional Airfoil problem, whereas qParEGO continues to improve and achieves the highest final HV. MORBO also demonstrates strong optimization performance. ViaMOBO converges relatively slowly during the early stage but improves consistently throughout the optimization process. Its final HV is approximately 2.16%2.16\% lower than that of TSEMO and 4.92%4.92\% lower than that of qParEGO, while remaining competitive with the other baselines. Table I and Figure 3 further illustrate the trade-off between optimization performance and computational cost. ViaMOBO requires only 0.799±0.1190.799± 0.119 hours per seed, making it the fastest surrogate-based method. Although TSEMO and qParEGO achieve higher final HV values, they require approximately 4.70×4.70× and 8.39×8.39× the runtime of ViaMOBO, respectively. MORBO requires 8.985±3.9738.985± 3.973 hours per completed seed, corresponding to 11.24×11.24× the runtime of ViaMOBO. DGEMO and ParEGO similarly require approximately 11.32×11.32× and 22.76×22.76× the runtime, respectively, without achieving better final performance than qParEGO or MORBO. Therefore, ViaMOBO provides a more favorable balance between solution quality and computational efficiency. Fig. 4: (Left) The results of log HV difference and (Right) HV–runtime of all methods on airfoil problem with 2020 decision variables. Overall, qParEGO achieves the highest final solution quality, while MORBO and TSEMO also exhibit strong optimization performance. In contrast, ViaMOBO attains competitive final performance using only approximately 21.3%21.3\%, 11.9%11.9\%, and 8.9%8.9\% of the runtimes required by TSEMO, qParEGO, and MORBO, respectively. These results demonstrate the suitability of ViaMOBO for high-dimensional practical optimization when computational resources are limited. IV-C2 40D Airfoil shape-optimization problem As shown in Figure 4 and Table I, MORBO achieves the highest final HV on the 40-dimensional Airfoil problem, but requires 9.371±5.1339.371± 5.133 hours on average, approximately 6.92×6.92× the runtime of ViaMOBO. In comparison, ViaMOBO attains a competitive final HV in only 1.353±0.1741.353± 0.174 hours, making it the most computationally efficient surrogate-based method among the completed runs. Although MORBO, qParEGO, DGEMO, and TSEMO improve the final HV over ViaMOBO by approximately 6.84%6.84\%, 4.92%4.92\%, 6.31%6.31\%, and 0.54%0.54\%, respectively, they require about 6.92×6.92×, 6.37×6.37×, 13.54×13.54×, and 6.03×6.03× more runtime. These results indicate that ViaMOBO provides a more favorable trade-off between Pareto-front quality and computational cost, particularly for high-dimensional expensive multi-objective optimization under limited computational budgets. MORBO’s computational cost is sensitive to its internal optimization trajectory. On the 100-dimensional problem, the trust regions remain relatively large, fewer invalid centers are encountered, and model fitting and candidate generation require only approximately 0.370.37–0.420.42 and 0.590.59–1.531.53 hours, respectively. By contrast, candidate generation alone requires approximately 8.588.58–12.0112.01 hours. Therefore, the shorter runtime observed in 100 dimensions reflects a less demanding search trajectory rather than improved scaling with dimensionality. This interpretation is further supported by its inferior optimization performance, with an AUC-HV of approximately 126.82126.82, compared with 148.91148.91 for ViaMOBO, together with a lower final HV. Accordingly, the reported MORBO runtime should be regarded as an empirical wall-clock measurement rather than evidence of reduced computational complexity at higher dimensions. A rigorous scalability comparison would require all dimensional settings to be evaluated under identical and isolated computational resources. Fig. 5: (Left) The results of HVs and (Right) log HV differences of all methods on bi-objective trajectory planning with 6060 decision variables. IV-C3 Trajectory planning problem Figure 5 shows the HV and log-HV difference of the evaluated algorithms on the 60-dimensional Rover trajectory-planning problem. The log -HV difference is log10(HV∗−HV) _10(HV^*-HV) where a lower value indicates better convergence toward the reference hypervolume. Figure 5, we can see that MORBO and NSGA-I achieve the best performance, respectively. ViaMOBO outperforms Random, ParEGO, MOEA/D-EGO, TSEMO, and USeMO-EI and remaining comparable to DGEMO, but falling behind MORBO and NSGA-I. MORBO exhibits both the fastest decrease and the lowest final value, followed by NSGA-I. Although ViaMOBO continuously reduces the hypervolume gap, its convergence becomes slower during the middle and later stages. These results suggest that the variable-grouping mechanism of ViaMOBO may omit important cross-group dependencies on strongly coupled sequential problems. Consequently, the Rover result illustrates an applicability boundary of ViaMOBO: it is better suited to high-dimensional problems with identifiable group structure or relatively weak inter-group coupling than to trajectory-optimization tasks with strong continuous coupling. IV-D Ablation study Finally, to study the effectiveness of ViaMOBO, we first show the model accuracy of our binary classifier SVM, and then study the characteristics of ViaMOBO with different acquisition functions, i.e., EI, UCB and EHVI. Fig. 6: The hypervolume performances of ViaMOBO with EI, UCB and EHVI on DTLZ2 with 1010 (Left) and 5050 (Right)decision variables. IV-D1 Accuracy of the binary classifier We record the mean of accuracy of 10 runs for each objective of 3-objective DTLZ2. Table I reports the variable-grouping accuracy of ViaMOBO on DTLZ2 with different dimensionalities. The overall grouping accuracy decreases from 90.68%90.68\% at D=10D=10 to 73.77%73.77\% at D=100D=100, indicating that identifying variable interactions becomes more challenging as the dimensionality increases. Nevertheless, ViaMOBO maintains an overall accuracy above 73%73\% even for the 100-dimensional problem, demonstrating that the grouping model remains effective in high-dimensional search spaces. The final-iteration accuracy follows a similar trend, decreasing from 90.23%90.23\% at D=10D=10 to 72.91%72.91\% at D=100D=100. The difference between the overall and final-iteration accuracies is small for all dimensionalities, suggesting that the learned grouping structures remain relatively stable throughout the optimization process. The larger standard deviation at D=100D=100 indicates increased variability across independent runs in higher-dimensional spaces. The D=50D=50 results should be interpreted with caution because only two independent runs are currently available. TABLE I: Variable-grouping model accuracy of ViaMOBO on DTLZ2. Dimension Overall accuracy (%) Final accuracy (%) 10 90.68±0.2690.68± 0.26 90.23±1.7890.23± 1.78 30 81.71±0.6381.71± 0.63 80.59±1.2680.59± 1.26 50 78.61±1.0478.61± 1.04 76.46±1.6776.46± 1.67 100 73.77±1.3273.77± 1.32 72.91±2.8072.91± 2.80 IV-D2 Performance of different acquisition functions This section study the HV and computation cost performances of different acquisition functions in our proposed ViaMOBO framework, including EI, UCB and EHVI. Figure 6 compares EI, UCB, and EHVI on the 2-objective DTLZ2 problem with D=10D=10 and D=50D=50 under identical experimental settings. It can bee seen that EHVI achieves better optimization performance because it directly maximizes the expected improvement in dominated hypervolume and is therefore naturally aligned with Pareto-front approximation. UCB, despite its computationally convenient additive form, provides a less direct measure of multi-objective improvement and consequently performs worse than EI and EHVI in this setting. This advantage, however, does not extend to the 3-objective problem. As shown in Figure 7, EHVI improves rapidly during the early stage but subsequently stagnates, whereas EI continues to improve and achieves the highest final HV. The deterioration of EHVI is associated with the rapidly increasing complexity of hypervolume partitioning and its complexity with the number of objectives. Figure 7 further shows that EHVI requires approximately 10.9×10.9× the runtime of EI, whereas EI and UCB require much less time. Moreover, several EHVI runs fail because of numerical instability during surrogate-model fitting. These results indicate that, although EHVI is attractive for low-dimensional problems with few objectives, EI and UCB provides a more favorable balance among convergence performance, computational efficiency, and numerical robustness for our ViaMOBO framework for high-dimensional real-world multi-objective optimization. V Conclusion This paper proposed ViaMOBO, a high-dimensional multi-objective Bayesian optimization framework that alleviates the curse of dimensionality by exploiting decision-variable interactions. Based on the redefined concepts of variable interaction and separability for multi-objective optimization, ViaMOBO employs a binary classifier to learn the underlying interaction structure from previously evaluated samples. By inferring interactions from the objective-value relationships between perturbed decision vectors, the proposed method avoids additional evaluations of expensive objective functions. When an additive structure is identified, the resulting variable groups are incorporated into a multi-objective additive kernel, effectively reducing the search dimensionality and facilitating acquisition-function optimization. Extensive experiments on high-dimensional synthetic benchmarks and real-world problems demonstrate that ViaMOBO can effectively exploit separable structures while achieving a favorable trade-off among solution quality, sample efficiency, and computational cost. Fig. 7: Comparison of HV (Left) and Wall time (Right) of three acquisition functions (EI, UCB and EHVI) under the same experimental settings. Nevertheless, ViaMOBO is primarily tailored to problems with separable or weakly coupled decision-variable structures, and its advantages may diminish for strongly coupled and non-separable problems. Future work will extend ViaMOBO toward a more general high-dimensional MOBO framework through more expressive interaction modeling, adaptive variable grouping, and enhanced capability for handling strongly coupled and non-separable expensive multi-objective problems. Acknowledgment This research was supported by the National Natural Science Foundation of China under Grant No. 52502380 and No. 52502508, and was also supported by the Postdoctoral Fellowship Program (Grade C) under Grant GZC20250905. References [1] M. Tesch, J. Schneider, and H. Choset, “Expensive multiobjective optimization for robotics,” in 2013 IEEE International Conference on Robotics and Automation, Karlsruhe, German, May 2013, p. 973–980. [2] W. Lyu, F. Yang, C. Yan, D. Zhou, and X. Zeng, “Batch Bayesian optimization via multi-objective acquisition ensemble for automated analog circuit design,” in Proceedings of the 35th International Conference on Machine Learning, ser. ICML’18, vol. 80. Stockholm, Sweden: PMLR, July 2018, p. 3306–3314. [3] Y. He and J. Cui, “Multi-Fidelity Surrogate Modeling via Hierarchical Kriging With Infinite-Width Bayesian Neural Network Correlation Function,” IEEE Access, vol. 13, p. 62 753–62 772, 2025. [4] “Multi-fidelity Kriging-assisted Preference-driven Multi-objective Optimization Method for Constrained Expensive Problems and Its Application,” Aerospace Science and Technology, p. 112258, 2026. [5] J. Snoek, H. Larochelle, and R. P. Adams, “Practical bayesian optimization of machine learning algorithms,” in Proceedings of the 25th International Conference on Neural Information Processing Systems, ser. NIPS’12, vol. 2. USA: MIT Press, 2012, p. 2951–2959. [6] J. Knowles, “Parego: a hybrid algorithm with on-line landscape approximation for expensive multiobjective optimization problems,” IEEE Transactions on Evolutionary Computation, vol. 10, no. 1, p. 50–66, Feb 2006. [7] W. Ponweiser, T. Wagner, D. Biermann, and M. Vincze, “Multiobjective optimization on a limited budget of evaluations using model-assisted s-metric selection,” in Parallel Problem Solving from Nature – PPSN X. Berlin, Heidelberg: Springer, 2008, p. 784–794. [8] Q. Zhang, W. Liu, E. Tsang, and B. Virginas, “Expensive multiobjective optimization by moea/d with gaussian process model,” IEEE Transactions on Evolutionary Computation, vol. 14, no. 3, p. 456–474, June 2010. [9] D. R. Jones, M. Schonlau, and W. J. Welch, “Efficient global optimization of expensive black-box functions,” Journal of Global Optimization, vol. 13, no. 4, p. 455–492, Dec 1998. [10] K. Yang and P. et al., “A multi-point mechanism of expected hypervolume improvement for parallel multi-objective bayesian global optimization,” in Proceedings of the Genetic and Evolutionary Computation Conference, ser. GECCO ’19. Prague, Czech Republic: ACM, July 2019, p. 656–663. [11] R. Zhang and D. Golovin, “Random hypervolume scalarizations for provable multi-objective black box optimization,” in Proceedings of the 37th International Conference on Machine Learning, ser. ICML’20, vol. 119. Virtual Event: PMLR, July 2020, p. 11 096–11 105. [12] S. Daulton, M. Balandat, and E. Bakshy, “Differentiable expected hypervolume improvement for parallel multi-objective bayesian optimization,” in Proceedings of the 33th International Conference on Neural Information Processing Systems, ser. NIPS’20. virtual: MIT Press, 2020. [13] M. Ranzato, A. Beygelzimer, Y. N. Dauphin, P. Liang, and J. W. Vaughan, Eds., Parallel Bayesian Optimization of Multiple Noisy Objectives with Expected Hypervolume Improvement, ser. NIPS’21, vol. 34. virtual: MIT Press, September 2021. [14] D. Hernández-Lobato, J. M. Hernández-Lobato, A. Shah, and R. P. Adams, “Predictive entropy search for multi-objective bayesian optimization,” in Proceedings of the 33rd International Conference on Machine Learning, ser. ICML’16, vol. 48. PMLR, 2016, p. 1492–1501. [15] E. C. Garrido-Merchán and D. Hernández-Lobato, “Predictive entropy search for multi-objective bayesian optimization with constraints,” Neurocomputing, vol. 361, p. 50–68, 2019. [16] S. Belakaria, A. Deshwal, and J. R. Doppa, “Max-value entropy search for multi-objective bayesian optimization,” in Proceedings of the 32th International Conference on Neural Information Processing Systems, ser. NIPS’19. BC, Canada: MIT Press, December 2019, p. 7823–7833. [17] D. Fernández-Sánchez, E. C. Garrido-Merchán, and D. Hernández-Lobato, “Max-value entropy search for multi-objective bayesian optimization with constraints,” CoRR, vol. abs/2011.01150, 2020. [18] S. Suzuki, S. Takeno, T. Tamura, K. Shitara, and M. Karasuyama, “Multi-objective bayesian optimization using pareto-frontier entropy,” in Proceedings of the 37th International Conference on Machine Learning, ser. ICML’20, vol. 119. Virtual Event: PMLR, July 2020, p. 9279–9288. [19] S. Belakaria, A. Deshwal, and J. R. Doppa, “Multi-fidelity multi-objective bayesian optimization: An output space entropy search approach,” in Proceedings of the 34th AAAI Conference on Artificial Intelligence, ser. AAAI’20, New York, USA, February 2020, p. 10 035–10 043. [20] M. Zuluaga, G. Sergent, A. Krause, and M. Püschel, “Active learning for multi-objective optimization,” in Proceedings of the 30th International Conference on Machine Learning, ser. ICML’13. GA, USA: PMLR, June 2013, p. 462–470. [21] B. Paria, K. Kandasamy, and B. Póczos, “A flexible framework for multi-objective bayesian optimization using random scalarizations,” in Proceedings of the 35th Conference on Uncertainty in Artificial Intelligence, ser. UAI’19, vol. 115. Tel Aviv, Israel: AUAI Press, July 2019, p. 766–776. [22] M. Konakovic Lukovic, Y. Tian, and W. Matusik, “Diversity-guided multi-objective bayesian optimization with batch evaluations,” in Proceedings of the 33th International Conference on Neural Information Processing Systems, ser. NIPS’20, vol. 33. virtual: MIT Press, September 2020. [23] S. Belakaria, A. Deshwal, N. K. Jayakodi, and J. R. Doppa, “Uncertainty-aware search framework for multi-objective bayesian optimization,” in Proceedings of the 34th AAAI Conference on Artificial Intelligence, ser. AAAI’20, New York, USA, February 2020, p. 10 044–10 052. [24] K. Kandasamy, J. Schneider, and B. Póczos, “High dimensional bayesian optimisation and bandits via additive models,” in Proceedings of the 32nd International Conference on Machine Learning, ser. ICML’15, vol. 37. PMLR, 2015, p. 295–304. [25] Z. Wang, C. Li, S. Jegelka, and P. Kohli, “Batched high-dimensional bayesian optimization via structural kernel learning,” in Proceedings of the 34th International Conference on Machine Learning, ser. ICML’17, vol. 70. Sydney, NSW, Australia: PMLR, August 2017, p. 3656–3664. [26] H. Qian and Y. Yu, “Solving high-dimensional multi-objective optimization problems with low effective dimensions,” in Proceedings of the 31st AAAI Conference on Artificial Intelligence, ser. AAAI’17. AAAI Press, February 2017, p. 875–881. [27] S. Daulton, D. Eriksson, M. Balandat, and E. Bakshy, “Multi-objective bayesian optimization over high-dimensional search spaces,” in Uncertainty in Artificial Intelligence. PMLR, 2022, p. 507–517. [28] D. Eriksson, M. Pearce, and J. R. G. et al., “Scalable global optimization via local bayesian optimization,” in Proceedings of the 32th International Conference on Neural Information Processing Systems, ser. NIPS’19. BC, Canada: MIT Press, December 2019, p. 5497–5508. [29] L. Györfi, M. Kohler, A. Krzyzak, and H. Walk, A Distribution-Free Theory of Nonparametric Regression, ser. Springer series in statistics. Springer, 2002. [30] W. Chen, T. Weise, Z. Yang, and K. Tang, “Large-scale global optimization using cooperative coevolution with variable interaction learning,” in Parallel Problem Solving from Nature - PPSN XI, vol. 6238. Kraków, Poland: Springer, September 2010, p. 300–309. [31] X. Ma, F. Liu, Y. Qi, X. Wang, L. Li, L. Jiao, M. Yin, and M. Gong, “A multiobjective evolutionary algorithm based on decision variable analyses for multiobjective optimization problems with large-scale variables,” IEEE Transactions on Evolutionary Computation, vol. 20, no. 2, p. 275–298, 2016. [32] Z. Bouzarkouna, A. Auger, and D. Y. Ding, “Local-meta-model CMA-ES for partially separable functions,” in 13th Annual Genetic and Evolutionary Computation Conference, ser. GECCO ’11. Dublin, Ireland: ACM, 2011, p. 869–876. [33] S. Huband, P. Hingston, L. Barone, and L. While, “A review of multiobjective test problems and a scalable test problem toolkit,” IEEE Transactions on Evolutionary Computation, vol. 10, no. 5, p. 477–506, Oct 2006. [34] C. Hvarfner, E. O. Hellsten, and L. Nardi, “Vanilla Bayesian optimization performs great in high dimensions,” in Proceedings of the 41st International Conference on Machine Learning, ser. Proceedings of Machine Learning Research, vol. 235. PMLR, 2024, p. 20 793–20 817. [35] Z. Wang, M. Zoghi, F. Hutter, D. Matheson, and N. De Freitas, “Bayesian optimization in high dimensions via random embeddings,” in Proceedings of the 23rd International Joint Conference on Artificial Intelligence, ser. IJCAI’13. Beijing, China: AAAI Press, August 2013, p. 1778–1784. [36] J. Djolonga, Krause, and A. et al., “High-dimensional gaussian process bandits,” in Proceedings of the 26th International Conference on Neural Information Processing Systems, ser. NIPS’13. USA: MIT Press, 2013, p. 1025–1033. [37] B. Letham, R. Calandra, A. Rai, and E. Bakshy, “Re-examining linear embeddings for high-dimensional bayesian optimization,” in Proceedings of the 33th International Conference on Neural Information Processing Systems, ser. NIPS’20. virtual: MIT Press, 2020, p. 1025–1033. [38] E. Raponi, H. Wang, M. Bujny, S. Boria, and C. Doerr, “High dimensional bayesian optimization assisted by principal component analysis,” in Parallel Problem Solving from Nature - PPSN XVI, ser. Lecture Notes in Computer Science, vol. 12269. Leiden, The Netherlands: Springer, September 2020, p. 169–183. [39] H. Tran-The, S. Gupta, S. Rana, and S. Venkatesh, “Trading convergence rate with computational budget in high dimensional bayesian optimization,” in Proceedings of the 34th AAAI Conference on Artificial Intelligence, ser. AAAI’20. New York, USA: AAAI Press, February 2020, p. 2425–2432. [40] C. Li, K. Kandasamy, B. Póczos, and J. G. Schneider, “High dimensional bayesian optimization via restricted projection pursuit models,” in Proceedings of the 19th International Conference on Artificial Intelligence and Statistics, ser. AISTAS’16, vol. 51. Cadiz, Spain: PMLR, May 2016, p. 884–892. [41] P. Rolland, J. Scarlett, I. Bogunovic, and V. Cevher, “High-dimensional bayesian optimization via additive models with overlapping groups,” in Proceedings of the 21st International Conference on Artificial Intelligence and Statistics, ser. AISTAS’18, vol. 84. Canary Islands, Spain: PMLR, April 2018, p. 298–307. [42] H. Wang, H. Xu, and Y. Yuan, “High-dimensional expensive multi-objective optimization via additive structure,” Intelligent Systems with Applications, vol. 14, p. 200062, 2022. [43] H. Wang, H. Xu, and Z. Zhang, “High-dimensional multi-objective bayesian optimization with block coordinate updates: Case studies in intelligent transportation system,” IEEE Transactions on Intelligent Transportation Systems, vol. 25, no. 1, p. 884–895, 2024. [44] C. Oh, E. Gavves, and M. Welling, “BOCK : Bayesian optimization with cylindrical kernels,” in Proceedings of the 35th International Conference on Machine Learning, ser. ICML’18, Stockholm, Sweden, July 2018, p. 3865–3874. [45] E. Siivola, A. Vehtari, and J. V. et al., “Correcting boundary over-exploration deficiencies in bayesian optimization with virtual derivative sign observations,” in 28th IEEE International Workshop on Machine Learning for Signal Processing, MLSP 2018, Aalborg, Denmark, September 17-20, 2018, 2018, p. 1–6. [46] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, “A fast and elitist multiobjective genetic algorithm: Nsga-i,” IEEE Transactions on Evolutionary Computation, vol. 6, no. 2, p. 182–197, April 2002. [47] E. Bradford, A. M. Schweidtmann, and A. Lapkin, “Efficient multiobjective optimization employing gaussian processes, spectral sampling and a genetic algorithm,” Journal of Global Optimization, vol. 71, p. 407–438, 2018. [48] S. Ament, S. Daulton, D. Eriksson, M. Balandat, and E. Bakshy, “Unexpected improvements to expected improvement for bayesian optimization,” in Advances in Neural Information Processing Systems, vol. 36, 2023. [49] B. K. Maximilian Balandat and D. R. J. et al., “Botorch: A framework for efficient monte-carlo bayesian optimization,” in Proceedings of the 33th International Conference on Neural Information Processing Systems, ser. NIPS’20. virtual: MIT Press, 2020. [50] C. M. Aye, K. Wansaseub, S. Kumar, G. G. Tejani, S. Bureerat, A. R. Yildiz, and N. Pholdee, “Airfoil shape optimisation using a multi-fidelity surrogate-assisted metaheuristic with a new multi-objective infill sampling technique,” Computer Modeling in Engineering & Sciences, vol. 137, no. 3, p. 2111–2128, 2023.