Paper deep dive
Adapting Nonstationary Multi-output Gaussian Processes to Bayesian Optimization
Zikai Xie
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) commonly relies on independent Gaussian processes (GPs) with stationary kernels, limiting its ability to represent nonstationary structure and share information between objectives. However, expressive nonstationary GPs do not necessarily make reliable BO decisions. We study this mismatch for the multi-output low-rank nonstationary (MO-LRN) GP: strong training fit can coexist with large off-design errors and optimistic acquisition predictions. We introduce MOLRN-BO, which combines a regularized shared-spectral surrogate with objective-specific residuals, prequential mean correction and tempered covariance scaling, and Pareto-local qLogEHVI optimization with periodic global search. Experiments on 12 deterministic bi-objective benchmarks show that MOLRN-BO substantially improves upon the original MO-LRN and achieves the best average problem ranks for final normalized hypervolume and normalized inverted generational distance among nine evaluated algorithms. It also achieves the strongest adverse-tail performance while remaining competitive with the leading baselines in anytime optimization. Ablation studies further show that the shared spectral construction improves off-design prediction, the local--global decision policy improves optimization performance, and hierarchical calibration reduces systematic candidate bias. These results demonstrate that nonstationary multi-output surrogates can deliver strong and robust MOBO performance when their structure and use are explicitly adapted to the demands of sequential optimization.
Tags
Links
- Source: https://arxiv.org/abs/2609.32464v1
- Canonical: https://arxiv.org/abs/2609.32464v1
Trouble viewing inline? Open PDF directly →
Full Text
79,695 characters extracted from source content.
Expand or collapse full text
Adapting Nonstationary Multi-output Gaussian Processes to Bayesian Optimization Zikai Xie †thanks: Corresponding author. Affiliation: State Key Laboratory of Precision and Intelligent Chemistry Affiliation: University of Science and Technology of China Affiliation: Hefei, Anhui 230000, China Email: zikaix@ustc.edu.cn Abstract Multi-objective Bayesian optimization (MOBO) commonly relies on independent Gaussian processes (GPs) with stationary kernels, limiting its ability to represent nonstationary structure and share information between objectives. However, expressive nonstationary GPs do not necessarily make reliable BO decisions. We study this mismatch for the multi-output low-rank nonstationary (MO-LRN) GP: strong training fit can coexist with large off-design errors and optimistic acquisition predictions. We introduce MOLRN-BO, which combines a regularized shared-spectral surrogate with objective-specific residuals, prequential mean correction and tempered covariance scaling, and Pareto-local qLogEHVI optimization with periodic global search. Experiments on 12 deterministic bi-objective benchmarks show that MOLRN-BO substantially improves upon the original MO-LRN and achieves the best average problem ranks for final normalized hypervolume and normalized inverted generational distance among nine evaluated algorithms. It also achieves the strongest adverse-tail performance while remaining competitive with the leading baselines in anytime optimization. Ablation studies further show that the shared spectral construction improves off-design prediction, the local–global decision policy improves optimization performance, and hierarchical calibration reduces systematic candidate bias. These results demonstrate that nonstationary multi-output surrogates can deliver strong and robust MOBO performance when their structure and use are explicitly adapted to the demands of sequential optimization. 1 Introduction Expensive scientific and engineering optimization often requires identifying good trade-offs among conflicting objectives with few black-box evaluations (Terayama et al., 2021; Shahriari et al., 2015; Luo et al., 2025). Bayesian optimization (BO) (Terayama et al., 2021) and its multi-objective variants (MOBO) (Laumanns and Ocenasek, 2002) address this problem by fitting a probabilistic surrogate and selecting experiments through an acquisition function. Its success depends on predictions at the points the acquisition function prefers, which may differ substantially from the observed training distribution. Methodological advances in MOBO have largely concentrated on acquisition design, including scalarization-, hypervolume-, and information-theoretic criteria (Knowles, 2006; Daulton et al., 2020; Daulton et al., 2021; Hernández-Lobato et al., 2016; Belakaria et al., 2019; Tu et al., 2022). By contrast, standard implementations typically retain independent GPs with stationary Matérn or RBF kernels, which neither explicitly share information across objectives nor represent input-dependent covariance structure. Recently, the multi-output low-rank nonstationary (MO-LRN) kernel offers a more expressive alternative through low-rank spectral factors (Xu et al., 2026). Yet, our preliminary experiments show that increasing its spectral capacity improved training fit while worsening off-design prediction. In addition, predictions at acquisition-selected candidates were biased and underdispersed. This regression–decision gap motivates our central question: how should MO-LRN be structured and deployed as a reliable surrogate for MOBO? To address these issues, we introduce MOLRN-BO, a decision-oriented nonstationary multi-output BO framework. MOLRN-BO stabilizes the surrogate through regularized shared LRN components with per-mixture coregionalization and objective-specific ARD-RBF residuals. It then constructs a decision posterior using strictly prequential bias estimates and bounded covariance rescaling, and optimizes qLogEHVI through Pareto-centered local search with periodic global rounds. Together, these components align surrogate fitting, predictive calibration, and acquisition optimization with the requirements of sequential multi-objective decision making. Our contributions are threefold. First, we identify and characterize the mismatch between training fit, off-design prediction, and acquisition decisions for MO-LRN and develop, to our knowledge, the first MOBO framework built on a harmonizable low-rank nonstationary multi-output GP. Second, evaluations on 12 bi-objective benchmarks show that MOLRN-BO achieves state-of-the-art final Pareto quality among the evaluated methods with the best average ranks in both normalized hypervolume and normalized IGD, while maintaining competitive anytime performance. Third, controlled ablations quantify the effects of spectral capacity and sharing, joint local–global search and mode-specific bias estimation, and hierarchical bias pooling on prediction and optimization. 2 Related Work Multi-objective surrogates and acquisitions. Research in MOBO has primarily focused on acquisition design. Scalarization-based methods, including ParEGO and random scalarization, reduce the vector-valued problem to a sequence of scalar subproblems (Knowles, 2006; Paria et al., 2020). Hypervolume-based methods, such as EHVI, qEHVI, qNEHVI, and their logarithmic variants, directly reward improvement in dominated hypervolume (Emmerich et al., 2006; Daulton et al., 2020; Daulton et al., 2021). Information-theoretic methods, including PESMO, MESMO, JES, and PF2ES, instead select evaluations that reduce uncertainty about the Pareto set or Pareto front (Hernández-Lobato et al., 2016; Belakaria et al., 2019; Tu et al., 2022). Surrogate design has received comparatively less attention: beyond the independent stationary GPs used by standard methods, a smaller line of work considers correlated multi-task or multi-output GPs and deep GP surrogates (Shah and Ghahramani, 2016; Chowdhury and Gopalan, 2021; Hebbal et al., 2023). However, nonstationary spectral multi-output surrogates remain largely unexplored in acquisition-driven MOBO. Nonstationary spectral models. Multi-output GPs place a joint prior over vector-valued functions through a matrix-valued kernel whose diagonal entries describe within-output covariance and whose off-diagonal entries describe cross-output dependence. Classical constructions such as the intrinsic coregionalization model and the linear model of coregionalization obtain valid cross-covariances by mixing latent scalar GPs (Alvarez et al., 2012). Spectral alternatives instead model these dependencies in the frequency domain: MOSM represents output-specific magnitudes, delays, and phases but remains stationary (Parra and Tobar, 2017), whereas MOHSM uses a harmonizable double-frequency representation to permit input-dependent auto- and cross-covariances (Altamirano and Tobar, 2022). MO-LRN further combines a generalized spectral–kernel duality with a low-rank factorization of flexible matrix-valued nonstationary spectral densities, reducing parameter growth in the number of outputs from quadratic to linear (Xu et al., 2026). Previous evaluations of MO-LRN have focused on passive prediction; we examine the additional demands that arise when it is deployed inside an acquisition-driven MOBO loop. Calibration and localized search. Acquisition functions depend on both posterior means and posterior uncertainties, making predictive calibration especially important under model misspecification and adaptive sampling. Previous work has used conformal prediction, posterior recalibration, and online residual correction to improve predictive coverage or quantile estimates in BO (Stanton et al., 2023; Deshpande et al., 2024). Most of these methods are developed for single-output BO or focus primarily on marginal predictive coverage, and calibration for correlated nonstationary surrogates under hypervolume-based acquisition remains comparatively unexplored. On the other hand, localized BO provides a complementary approach to stabilizing optimization in large or difficult search spaces. TuRBO maintains adaptive trust regions for single-objective BO, whereas MORBO coordinates multiple local regions to cover different parts of a multi-objective Pareto front (Eriksson et al., 2019; Daulton et al., 2022). Both approaches couple localization with local surrogate modeling and search. 3 Method Figure 1: Development diagnostics for the original MO-LRN surrogate on WFG4. (a) One-dimensional posterior slice at n=20n=20, obtained by varying x9x_9 while fixing the remaining coordinates at their training-data medians. Solid curves show posterior means, shaded regions show ±2± 2 posterior standard deviations, and dashed curves show the true objectives. (b) Strictly prequential standardized prediction errors, |y−μ|/σpred|y-μ|/ _pred, at acquisition-selected candidates over three trials. Predictions are evaluated before each candidate is added to the training set; the dashed line at 1.645 denotes the nominal two-sided 90% predictive limit. (c) Training, global, and Pareto-region RMSE together with top-HVI optimism for Q=2Q=2 and Q=4Q=4 MO-LRN models fitted to identical observations at n=110n=110. Bars denote means and points denote 15 fits from three data trajectories and five restarts. We maximize :⊂ℝD→ℝMf:X ^D ^M on a bounded box. After n evaluations, let n=(t,t)t=1nD_n=\(x_t,y_t)\_t=1^n and let nP_n be its nondominated objective vectors. Experiments use deterministic objectives and one candidate per round. Objective-wise standardization is recomputed on nD_n, with standard deviations ns_n. Raw posterior moments n() μ_n(x) and n() _n(x) below are on the original objective scale. From surrogate fit to decision reliability. Sequential model-based optimization generally relies on the premise that a more accurate surrogate provides a better basis for subsequent optimization decisions. However, when applying the original MO-LRN surrogate to WFG4, we observed that improved regression fit did not consistently translate into a more reliable posterior at the locations relevant to BO. Figure 1a provides a representative example: along a one-dimensional slice, MO-LRN with higher spectral capacity produces strongly distorted posterior means together with relatively narrow uncertainty bands, despite being fitted to the same observations. Compared to a stationary RBF kernel GP reference, this imbalance between posterior mean and uncertainty is particularly problematic for acquisition-based optimization, where confident extrapolation can directly affect candidate selection. Consistent with this observation, strictly prequential predictions at acquisition-selected candidates are severely underdispersed (Fig. 1b). Moreover, increasing the number of spectral mixtures from Q=2Q=2 to Q=4Q=4 reduces training RMSE while worsening both off-design prediction and HVI estimation (Fig. 1c). Additional pilot diagnostics revealed related symptoms, including spectral frequencies approaching their imposed bounds, weak overlap between output-specific spectra, and a tendency for unrestricted acquisition optimization to select distant or boundary-adjacent candidates (Appendix A). These observations suggest that the difficulty is not captured by regression accuracy alone. Instead, reliable use of MO-LRN in sequential optimization requires controlling spectral capacity, calibrating the posterior specifically at decision-relevant locations, and preventing acquisition optimization from overexploiting poorly supported regions. These considerations motivate the three corresponding components of MOLRN-BO developed below. Shared-spectral covariance. To share nonstationary structure across objectives without relying on the overlap of independently fitted output spectra, we use the LMC-style covariance Kij(,′)=∑q=1Q[q]ijkqsLRN(,′)+δijkiRBF(,′),q⪰0.K_ij(x,x )= _q=1^Q[B_q]_ijk_q^sLRN(x,x )+ _ijk_i^RBF(x,x ), _q 0. (1) Each scalar component kqsLRNk_q^sLRN retains the MO-LRN nonstationary construction with spectral parameters shared across objectives. The matrix qB_q represents component-specific cross-objective dependence, while an ARD-RBF residual kiRBFk_i^RBF captures smooth variation specific to objective i. We use Q=2Q=2 and bi-objective coregionalization matrices with unit diagonals. The coregionalization and residual parameterizations and the positive-semidefiniteness argument are given in Appendix B.1. Regularized fitting. To control spectral complexity during sequential refitting, we fit the exact GP to the standardized objectives ~n Y_n by minimizing θ^n∈argminθ∈Θ[−logp(~n∣n,θ)+ℛspectral(θ)+ℛother(θ)]. θ_n∈ _θ∈ [- p( Y_n _n,θ)+R_spectral(θ)+R_other(θ) ]. (2) The parameter domain Θ bounds frequency locations and spectral scales. The spectral penalty discourages rough functions, boundary-saturating frequencies, and mixture imbalance; the remaining penalties regularize coregionalization, residual, noise, and outputscale parameters. Each fit is warm-started from the preceding iteration, with objective standardization recomputed on nD_n. The full objective, parameter transformations, and hyperparameter settings are specified in Appendix B.1. Prequential mean correction. We correct systematic prediction errors using residuals recorded at previously selected candidates. Before evaluating tx_t, we store its raw GP moments, the current output scale t−1s_t-1, and the search mode St∈local,globalS_t∈\local,global\. After observing ty_t, we compute ~t=[t−t−1(t)]⊘t−1 r_t=[y_t- μ_t-1(x_t)] _t-1, where ⊘ denotes element-wise division. For the next search mode s, the correction is ^s,n=∑t∈n:St=swn,t~tclip+λs^0,n∑t∈n:St=swn,t+λs. b_s,n= _t _n:S_t=sw_n,t r_t^\,clip+ _s b_0,n _t _n:S_t=sw_n,t+ _s. (3) Here nW_n is a recent residual window, wn,tw_n,t discounts older observations, and clipping limits outlier influence. The shared estimate ^0,n b_0,n pools both modes and is shrunk toward zero; λs>0 _s>0 partially pools each mode toward this estimate when its own history is limited. Every residual used for the next decision comes from a prediction stored before the corresponding evaluation. Clipping, weighting, and pooling are detailed in Appendix B.2. Covariance rescaling and decision posterior. To address underdispersion without passing the full estimated scale correction to acquisition optimization, we distinguish predictive and acquisition scales. We estimate cpred,nc_pred,n using fixed-hyperparameter leave-one-scalar-out and spatial cross-validation, combined with Mahalanobis errors of past candidates’ stored raw predictions. For search mode s, the decision posterior is pacq(()∣n,Sn+1=s) p_acq(f(x) _n,S_n+1=s) =(n()+n⊙^s,n,cacq,nn()), =N\! ( μ_n(x)+s_n b_s,n,\;c_acq,n _n(x) ), (4) cacq,n c_acq,n =minccap,cpred,nγ. = \c_cap,c_pred,n^γ\. Here ⊙ denotes element-wise multiplication. Tempering and capping limit the covariance inflation exposed to qLogEHVI. For predictive diagnostics, we instead retain the raw mean and use cpred,nnc_pred,n _n. These posterior transformations leave the fitted GP unchanged, and scalar rescaling preserves its correlation structure. Cross-validation uses full-data hyperparameters; only the stored candidate errors are strictly prequential. The scale estimators and their settings are given in Appendix B.2. Pareto-local and periodic-global selection. We concentrate most acquisition searches around diverse observed Pareto points, while scheduling periodic global rounds to retain full-domain exploration. Let ℛn,k(rn)R_n,k(r_n) denote a box around the k-th selected Pareto center, with a common relative half-width rnr_n. The next candidate is obtained by n _n =⋃k=1Knℛn,k(rn),Sn+1=local,,Sn+1=global, = cases _k=1^K_nR_n,k(r_n),&S_n+1=local,\\ X,&S_n+1=global, cases (5) n+1 _n+1 ≈argmax∈nqLogEHVI(,pacq,n,HV), ≈ _x _nqLogEHVI(x;p_acq,P_n,r_HV), where HVr_HV is the acquisition reference point. On local rounds, farthest-point selection in normalized decision and objective coordinates determines the centers, and the restart budget is distributed approximately equally across their boxes. We select the best numerically optimized candidate. The common radius adapts to realized HVI on local rounds; global rounds leave the radius and local counters unchanged and occur more frequently early in optimization. The GP is always fitted to all observations: only acquisition optimization is localized. Center selection, restart allocation, scheduling, and radius updates are specified in Appendix B.3. Algorithm 1 MOLRN-BO 1: Initial data n0D_n_0, additional evaluation budget T, and acquisition reference point HVr_HV 2: Initialize empty calibration history, common local radius, and counters 3: for n=n0,…,n0+T−1n=n_0,…,n_0+T-1 do 4: Standardize objectives and fit the shared-spectral GP using Eqs. equation 1 and equation 2 5: Set u=n−n0+1u=n-n_0+1; choose Sn+1S_n+1 and construct nA_n using the local–global policy (Appendix B.3) 6: Estimate ^Sn+1,n b_S_n+1,n from past prequential residuals using Eq. equation 3 7: Estimate cpred,nc_pred,n from cross-validation on nD_n and stored raw candidate errors (Appendix B.2) 8: Construct pacqp_acq using Eq. equation 4 9: Obtain n+1x_n+1 by optimizing Eq. equation 5 10: Store n(n+1) μ_n(x_n+1), n(n+1) _n(x_n+1), ns_n, and Sn+1S_n+1 11: Evaluate n+1=(n+1)y_n+1=f(x_n+1) 12: Append the standardized residual and raw Mahalanobis error to the calibration history 13: n+1←n∪(n+1,n+1)D_n+1 _n∪\(x_n+1,y_n+1)\ 14: if Sn+1=localS_n+1=local then 15: Update the common radius and local counters using realized HVI (Appendix B.3) 16: end if 17: end for 18: return Nondominated observations in n0+TD_n_0+T 4 Experiments 4.1 Experimental settings Problems. We evaluate 12 deterministic bi-objective problems: Branin–Currin (D=2D=2) (Daulton et al., 2020); WFG1, WFG2, and WFG4–WFG9 (D=10D=10, k=4k=4) from the Walking Fish Group (WFG) suite (Huband et al., 2006); and DTLZ2, DTLZ4, and DTLZ7 (D=10D=10) from the DTLZ suite (Deb et al., 2002b). Branin–Currin tests sample-efficient modeling of differently shaped objective surfaces in a low-dimensional domain; the selected WFG problems introduce combinations of multimodality, deceptiveness, nonseparability, biased transformations, and varied Pareto-front geometries in higher-dimensional decision spaces; and the DTLZ problems assess front approximation under a regular concave geometry (DTLZ2), a biased decision-to-objective mapping (DTLZ4), and a disconnected front (DTLZ7). Comparators and budget. The eight baselines are independent BoTorch SingleTaskGP models (Balandat et al., 2020) with qLogEHVI or qLogNParEGO; lower-bound JES (Tu et al., 2022); MESMO (Belakaria et al., 2019); mobopt (Galuzio et al., 2020); original MO-LRN with Q=2Q=2 and qLogEHVI; NSGA-I (Deb et al., 2002a); and scrambled Sobol search. To ensure statistical reliability, each algorithm–task pair was evaluated over 20 trials using random seeds 100–119. Each problem–seed pair shares a 10-point scrambled Sobol initial design, followed by 100 sequential evaluations (batch size one). Runs were distributed across heterogeneous CPU/GPU systems, so no controlled runtime comparison is reported. Metrics and statistical analysis. All methods are rescored against fixed reference fronts under a common maximization convention, with the hypervolume reference point placed 10%10\% beyond the nadir in the worse direction along each objective. Normalized hypervolume (nHV; higher is better) is the hypervolume of the observed nondominated set divided by that of the fixed reference front; normalized inverted generational distance (nIGD; lower is better) uses ideal–nadir-normalized objectives. The four primary metrics are final nHV and nIGD and two anytime measures: normalized trapezoidal areas under the HV-regret and nIGD curves over sequential evaluation steps (both AUCs are lower-is-better). Secondary analyses assess run-to-run robustness using across-seed standard deviations and adverse-tail performance: the 10th percentile of nHV and the 90th percentiles of nIGD and both AUC measures. Additional IGD+, additive-ϵε, and convergence-curve results are reported in Appendix C. For each primary metric, we rank algorithms by their across-seed means within each problem and average ranks over the 12 problems. Friedman tests assess overall algorithm differences. Comparisons between MOLRN-BO and each baseline use Wilcoxon signed-rank tests: one-sided tests on the 12 paired problem-level means across problems, and two-sided tests on paired common-seed observations within each problem. Holm correction is applied separately for each metric, across baselines for cross-problem tests and across all problem–baseline comparisons for within-problem tests. 4.2 Experimental results Table 1: Mean problem ranks over 12 benchmarks using all valid completed trials. Lower rank is better. Bold and underlined entries denote the best and second-best method for each metric, respectively. Method Final nHV Final nIGD HV-regret AUC nIGD AUC MOLRN-BO 2.25 2.17 2.42 2.83 qLogNParEGO 2.42 3.00 2.17 2.75 qLogEHVI 3.42 4.33 3.33 3.83 mobopt 4.67 4.25 4.92 4.00 JES 5.50 4.83 5.50 4.75 Original MO-LRN 6.42 5.33 6.33 5.58 NSGA-I 5.54 6.50 6.04 7.00 MESMO 7.42 7.25 6.92 6.67 Sobol 7.38 7.33 7.38 7.58 Overall performance. As shown in Table 1, MOLRN-BO obtains the best average rank for both final Pareto-quality metrics: 2.252.25 for nHV and 2.172.17 for nIGD. qLogNParEGO ranks second at the final budget but is slightly better on both anytime AUC measures, placing the two methods in the leading group with different strengths. Algorithm effects are significant for all four metrics (Friedman p<1.7×10−6p<1.7× 10^-6). (a) Normalized hypervolume regret. (b) Normalized IGD. Figure 2: Per-problem optimization trajectories on the 12 benchmarks. Solid curves show across-trial means, and shaded regions indicate one sample standard deviation. Both metrics are lower-is-better. Figure 2 shows the per-problem HV-regret and nIGD trajectories. qLogNParEGO often progresses faster initially, whereas MOLRN-BO is competitive or leading at the final budget on many problems. This pattern may reflect the data requirements of the spectral surrogate and the limited early residual history. With only ten initial observations, identifying input-dependent spectral structure, mixture contributions, and cross-objective dependence may be less reliable than fitting a simpler stationary surrogate. Prequential corrections also have little residual history at this stage and are absent or strongly shrunk. As observations accumulate, both spectral fitting and mode-specific corrections can draw on more evidence. Comparison with the original MO-LRN. The direct comparison with the original MO-LRN evaluates the benefit of the complete MOLRN-BO framework. MOLRN-BO achieves better problem-level means on 11/1211/12 problems for final nHV and HV-regret AUC, and on all 12/1212/12 problems for final nIGD and nIGD AUC. The corresponding cross-problem tests remain significant after Holm correction (p≤0.00244p≤ 0.00244). The contribution of each individual module is examined separately in Sec. 4.3. Comparison with standard MOBO baselines. Compared with qLogEHVI, MOLRN-BO achieves better mean final nIGD on 10/1210/12 problems (cross-problem p=0.0637p=0.0637), including six statistically significant improvements; all significant differences favor MOLRN-BO. It also obtains better mean final nHV on 7/127/12 problems. The comparison with qLogNParEGO is closer: MOLRN-BO performs better on 6/126/12 problems in final nHV and 7/127/12 in final nIGD, whereas qLogNParEGO has slightly better anytime ranks. Although the aggregate differences between these two methods are not statistically significant, MOLRN-BO achieves slightly stronger final Pareto quality overall while remaining competitive throughout optimization. Robustness and limitations. MOLRN-BO has the smallest average trial-to-trial nIGD standard deviation (0.03580.0358) among the nine algorithms. It also obtains the best adverse-tail mean rank for all four primary outcomes: 2.082.08 for nHV-q10, 2.082.08 for nIGD-q90, 2.002.00 for HV-regret-AUC-q90, and 2.832.83 for nIGD-AUC-q90. These results are consistent with the design objective of reducing exposure to poor trials; module-level attribution is examined by the ablations in Sec. 4.3. However, the gains are not universal: MOLRN-BO is comparatively weak on WFG5, converges more slowly than the strongest baselines on WFG1, and exhibits slow early-to-middle nIGD improvement on DTLZ7. Detailed robustness statistics and per-problem curves are provided in Appendix C. 4.3 Ablation study We evaluate the principal surrogate and decision components under the experimental protocol described in Sec. 4.1, using matched initial designs and the same reference fronts and normalization. The surrogate diagnostic evaluates three models on a fixed set of 110 observations collected by MOLRN-BO: the original output-specific parameterization with Q=4Q=4, the same parameterization with Q=2Q=2, and the shared-spectrum Q=2Q=2 construction in Eq. equation 1. The RBF residual is disabled in all three variants. Consequently, the first transition isolates spectral capacity, while the second isolates spectral sharing and per-mixture coregionalization, conditional on a common observed design. The decision-policy ablation holds the shared-Q=2Q=2 and RBF surrogate, covariance calibration, and qLogEHVI fixed. Its control optimizes all restarts over the full domain and uses one bias estimate pooled across acquisition modes. The treatment instead applies the Pareto-local/periodic-global search policy and estimates independent biases for the local and global modes. This comparison measures the joint effect of search localization and mode-specific bias estimation. Finally, the bias-pooling ablation holds this local–global policy fixed and replaces the independent mode-specific estimates with the hierarchical estimator in Eq. equation 3. Table 2: Controlled diagnostics and ablations over 12 benchmarks. Values are means over available paired problem–seed observations, except candidate bias, reported as the median problem-level absolute bias. Panel (a) disables all RBF residuals. (a) Fixed-design surrogate Training RMSE Global RMSE Pareto-region RMSE Top-HVI optimism Output-specific, Q=4Q=4 0.018 1.154 1.101 1.213 Output-specific, Q=2Q=2 0.062 0.816 0.792 0.498 Shared-spectrum, Q=2Q=2 0.063 0.746 0.729 0.087 (b) Search / bias Final nHV ↑ Final nIGD ↓ HV-regret AUC ↓ nIGD AUC ↓ Full-domain / pooled 0.459 0.255 0.633 0.349 Local–global / independent 0.590 0.201 0.570 0.323 (c) Bias pooling Final nHV ↑ Final nIGD ↓ Median absolute candidate bias ↓ Independent 0.590 0.201 0.197 Hierarchical (MOLRN-BO) 0.582 0.207 0.079 Spectral structure. Reducing output-specific capacity from Q=4Q=4 to Q=2Q=2 increases training RMSE but lowers global and Pareto-region RMSE on 10/1210/12 problems and top-HVI optimism on 11/1211/12. At Q=2Q=2, the shared construction further reduces mean global and Pareto-region RMSE and lowers mean top-HVI optimism from 0.4980.498 to 0.0870.087, despite nearly unchanged training RMSE (Table 2a). These fixed-design diagnostics favor a lower-capacity shared surrogate for off-design prediction and acquisition estimation. Decision policy. The local–global/independent-bias bundle improves both final and anytime performance (Table 2b). Final nHV, final nIGD, HV-regret AUC, and nIGD AUC improve on 1111, 1111, 1111, and 1010 of the 12 problems, respectively. WFG1 is the principal exception for the final and hypervolume-based measures, and DTLZ4 for nIGD AUC. Bias pooling. Hierarchical pooling reduces candidate absolute bias on 11/1211/12 problems, with its problem-level median decreasing from 0.1970.197 to 0.0790.079. This reduction does not translate into consistent BO gains: mean final nHV decreases from 0.5900.590 to 0.5820.582, while nIGD increases from 0.2010.201 to 0.2070.207, and neither AUC improves consistently (Table 2c). The demonstrated benefit is therefore a substantial reduction in prediction bias at the cost of a slight decline in optimization performance. 5 Conclusion We introduced MOLRN-BO to address the mismatch between regression fit and acquisition-driven decisions when MO-LRN is used for multi-objective optimization. For the original surrogate, lower training error can coexist with larger off-design errors and overly optimistic acquisition predictions. The framework combines a regularized shared-spectral surrogate with objective-specific residuals, prequential mean correction, tempered covariance scaling, and Pareto-local search with periodic global rounds. Across 12 deterministic bi-objective benchmarks, it significantly outperforms the original MO-LRN on all four primary metrics in Holm-adjusted cross-problem tests and achieves the best average ranks for final nHV and nIGD, with second-best anytime ranks. Fixed-design diagnostics support using a lower-capacity spectral structure shared across objectives for off-design prediction, while sequential ablations show clear gains from combining local–global search with independent bias correction. Hierarchical pooling reduces systematic candidate bias while incurring only a small loss in BO performance. Together, these findings show that adapting both surrogate structure and its use in acquisition decisions can substantially improve sequential optimization performance compared with relying on training fit alone to select a surrogate. Despite these encouraging results, several limitations of this work should be acknowledged. The evaluation concerns small-budget, deterministic bi-objective optimization with one candidate per round. Larger budgets and additional outputs motivate sparse or variational inference, while noisy and batch settings require corresponding treatment of observation uncertainty and evaluation histories. Sensitivity to the fixed regularization settings, and the trade-off between optimization quality and computational overhead, remain to be quantified. AI use statement In this work, we used generative AI tools for assisting with translation, cleaning and reformatting datasets, supporting qualitative and thematic data analysis or interpreting results. We have not used generative AI tools for generating synthetic data sets or conceptual frameworks, designing or providing feedback on research methodology or experiments, and the rest of the required disclosure tasks are not applicable to this work. Additionally, we used generative AI tools for creating and editing software code, and editing the paper to improve readability. We have reviewed all AI-assisted work. We checked that the datasets reformatted by the LLM were consistent with the original datasets and that the statistical analysis tools provided by the LLM were correct. We take responsibility for the final content of this work, including text, claims or artifacts produced with the aid of generative AI. Ethics statement The authors have identified no ethical concerns associated with this work. Reproducibility statement The supplementary repository contains implementations of MOLRN-BO and all evaluated baselines, code for reproducing the ablation studies, and implementations of the 12 benchmark problems. Detailed formulations and hyperparameter settings for the shared-spectral surrogate, decision-posterior calibration, and local–global acquisition optimization are provided in Appendix sections B.1, B.2, and B.3, respectively. Section 4.1 specifies the common initial designs, random seeds, evaluation budgets, performance metrics, and statistical tests. The ablation configurations are described in Section 4.3, and additional performance and robustness analyses are reported in Appendix C. References Altamirano and Tobar (2022) M. Altamirano and F. Tobar Nonstationary multi-output gaussian processes via harmonizable spectral mixtures. In International Conference on Artificial Intelligence and Statistics, p. 3204–3218. Cited by: §2. Alvarez et al. (2012) M. A. Alvarez, L. Rosasco, and N. D. Lawrence Kernels for vector-valued functions: a review. Foundations and Trends® in Machine Learning 4 (3), p. 195–266. Cited by: §2. Balandat et al. (2020) M. Balandat, B. Karrer, D. Jiang, S. Daulton, B. Letham, A. G. Wilson, and E. Bakshy BoTorch: a framework for efficient monte-carlo bayesian optimization. Advances in neural information processing systems 33, p. 21524–21538. Cited by: §4.1. Belakaria et al. (2019) S. Belakaria, A. Deshwal, and J. R. Doppa Max-value entropy search for multi-objective bayesian optimization. Advances in neural information processing systems 32. Cited by: §1, §2, §4.1. Chowdhury and Gopalan (2021) S. R. Chowdhury and A. Gopalan No-regret algorithms for multi-task bayesian optimization. In International Conference on Artificial Intelligence and Statistics, p. 1873–1881. Cited by: §2. Daulton et al. (2020) S. Daulton, M. Balandat, and E. Bakshy Differentiable expected hypervolume improvement for parallel multi-objective bayesian optimization. Advances in neural information processing systems 33, p. 9851–9864. Cited by: §1, §2, §4.1. Daulton et al. (2021) S. Daulton, M. Balandat, and E. Bakshy Parallel bayesian optimization of multiple noisy objectives with expected hypervolume improvement. Advances in neural information processing systems 34, p. 2187–2200. Cited by: §1, §2. Daulton et al. (2022) S. Daulton, D. Eriksson, M. Balandat, and E. Bakshy Multi-objective bayesian optimization over high-dimensional search spaces. In Uncertainty in Artificial Intelligence, p. 507–517. Cited by: §2. Deb et al. (2002a) K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan A fast and elitist multiobjective genetic algorithm: nsga-i. IEEE transactions on evolutionary computation 6 (2), p. 182–197. Cited by: §4.1. Deb et al. (2002b) K. Deb, L. Thiele, M. Laumanns, and E. Zitzler Scalable multi-objective optimization test problems. In Proceedings of the 2002 congress on evolutionary computation. CEC’02 (Cat. No. 02TH8600), Vol. 1, p. 825–830. Cited by: §4.1. Deshpande et al. (2024) S. Deshpande, C. Marx, and V. Kuleshov Online calibrated and conformal prediction improves bayesian optimization. In International conference on artificial intelligence and statistics, p. 1450–1458. Cited by: §2. Emmerich et al. (2006) M. T. Emmerich, K. C. Giannakoglou, and B. Naujoks Single-and multiobjective evolutionary optimization assisted by gaussian random field metamodels. IEEE Transactions on Evolutionary Computation 10 (4), p. 421–439. Cited by: §2. Eriksson et al. (2019) D. Eriksson, M. Pearce, J. Gardner, R. D. Turner, and M. Poloczek Scalable global optimization via local bayesian optimization. Advances in neural information processing systems 32. Cited by: §2. Galuzio et al. (2020) P. P. Galuzio, E. H. de Vasconcelos Segundo, L. dos Santos Coelho, and V. C. Mariani MOBOpt—multi-objective bayesian optimization. SoftwareX 12, p. 100520. Cited by: §4.1. Hebbal et al. (2023) A. Hebbal, M. Balesdent, L. Brevault, N. Melab, and E. Talbi Deep gaussian process for multi-objective bayesian optimization: a. hebbal et al.. Optimization and Engineering 24 (3), p. 1809–1848. Cited by: §2. Hernández-Lobato et al. (2016) D. Hernández-Lobato, J. Hernandez-Lobato, A. Shah, and R. Adams Predictive entropy search for multi-objective bayesian optimization. In International conference on machine learning, p. 1492–1501. Cited by: §1, §2. Huband et al. (2006) 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 10 (5), p. 477–506. Cited by: §4.1. Knowles (2006) J. Knowles ParEGO: a hybrid algorithm with on-line landscape approximation for expensive multiobjective optimization problems. IEEE transactions on evolutionary computation 10 (1), p. 50–66. Cited by: §1, §2. Laumanns and Ocenasek (2002) M. Laumanns and J. Ocenasek Bayesian optimization algorithms for multi-objective optimization. In International Conference on Parallel Problem Solving from Nature, p. 298–307. Cited by: §1. Luo et al. (2025) M. Luo, Z. Xie, H. Li, B. Zhang, J. Cao, Y. Huang, H. Qu, Q. Zhu, L. Chen, J. Jiang, et al. Physics-informed, dual-objective optimization of high-entropy-alloy nanozymes by a robotic ai chemist. Matter 8 (4). Cited by: §1. Paria et al. (2020) B. Paria, K. Kandasamy, and B. Póczos A flexible framework for multi-objective bayesian optimization using random scalarizations. In Uncertainty in Artificial Intelligence, p. 766–776. Cited by: §2. Parra and Tobar (2017) G. Parra and F. Tobar Spectral mixture kernels for multi-output gaussian processes. Advances in neural information processing systems 30. Cited by: §2. Shah and Ghahramani (2016) A. Shah and Z. Ghahramani Pareto frontier learning with expensive correlated objectives. In International conference on machine learning, p. 1919–1927. Cited by: §2. Shahriari et al. (2015) B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, and N. De Freitas Taking the human out of the loop: a review of bayesian optimization. Proceedings of the IEEE 104 (1), p. 148–175. Cited by: §1. Stanton et al. (2023) S. Stanton, W. Maddox, and A. G. Wilson Bayesian optimization with conformal prediction sets. In International conference on artificial intelligence and statistics, p. 959–986. Cited by: §2. Terayama et al. (2021) K. Terayama, M. Sumita, R. Tamura, and K. Tsuda Black-box optimization for automated discovery. Accounts of Chemical Research 54 (6), p. 1334–1346. Cited by: §1. Tu et al. (2022) B. Tu, A. Gandy, N. Kantas, and B. Shafei Joint entropy search for multi-objective bayesian optimization. Advances in Neural Information Processing Systems 35, p. 9922–9938. Cited by: §1, §2, §4.1. Xu et al. (2026) Q. Xu, Z. Yang, Y. Li, M. M. Zhang, and P. Olmos Revisiting nonstationary kernel design for multi-output gaussian processes. In International Conference on Learning Representations, Vol. 2026, p. 97426–97458. Cited by: §B.1, §1, §2. Appendix A Preliminary diagnoses The original parameterization of MO-LRN assigns separate spectral factors to each output, with cross-output covariance induced by their spectral overlap. However, pilot sequential experiments revealed five BO-oriented problems when MO-LRN was used directly in MOBO. Three arose at the surrogate level: in-sample fit was poorly aligned with candidate-region accuracy, learned frequencies frequently approached their imposed bounds, and cross-output spectral overlap could become numerically negligible. A fourth concerned the decision posterior, which exhibited biased and underdispersed predictions at acquisition-selected points. The fifth arose during acquisition optimization, where unrestricted full-domain search favored distant and boundary-adjacent candidates. Table 3 summarizes these diagnostics on the ten-dimensional WFG4 development problem. These failure modes together motivate three corresponding changes in Section 3: a more structured and regularized surrogate, an error-calibrated posterior used only for decisions, and an acquisition search strategy that distributes local refinement across the Pareto set while retaining global exploration. Table 3: Preliminary BO-oriented diagnostics on the ten-dimensional WFG4 development problem. These descriptive pilot results motivate the design of MOLRN-BO; broader controlled evaluations are presented in Section 4.3. Failure mode Pilot experiment and sample Key result Capacity–decision mismatch WFG4 at n=110n=110; three data trajectories and five model fits per trajectory; Q=4Q=4 and Q=2Q=2 fitted to identical data Q=4Q=4 vs. Q=2Q=2: train RMSE 0.00220.0022 vs. 0.00820.0082; global RMSE 0.95180.9518 vs. 0.77300.7730; Pareto-region RMSE 0.93720.9372 vs. 0.85300.8530. Boundary-saturating spectra Original Q=4Q=4 MO-LRN on WFG4; 300 sequential fits with fmax=8f_ =8 At least one learned frequency was within 5%5\% of the imposed bound in 296/300296/300 fits (98.7%98.7\%); the mean per-fit maximum frequency was 7.967.96, and normalized roughness averaged 59.659.6. Vanishing spectral overlap Original Q=4Q=4 MO-LRN on ten-dimensional WFG4; spectral snapshots at n=20,50,80,n=20,50,80, and 110110 Fractions of mixtures with Rq,12<10−12R_q,12<10^-12 were 47%,95%,100%,47\%,95\%,100\%, and 100%100\%, respectively. Biased and overconfident predictions Original Q=4Q=4 MO-LRN on WFG4; 300 strictly prequential predictions at acquisition-selected points Normalized RMSE was 3.903.90, while the empirical coverage of nominal 90%90\% posterior intervals was only 3.7%3.7\%. Extrapolative acquisition search WFG4; three paired pilot trials and 300 acquisition decisions per policy Full-domain vs. local–global search: nearest-data distance 0.7410.741 vs. 0.4970.497; boundary-coordinate fraction 24.5%24.5\% vs. 16.7%16.7\%. Notes. The normalized spectral overlap Rq,12∈[0,1]R_q,12∈[0,1] measures the overlap between the two output-specific spectral factors of mixture q. Model fits in the first row share the same data. Appendix B Calculation Detail of MOLRN-BO We consider the maximization of an expensive vector-valued black-box function :⊂ℝD→ℝMf:X ^D ^M. After n evaluations, the available data are n=(t,t)t=1nD_n=\(x_t,y_t)\_t=1^n, where t=(t)+ϵty_t=f(x_t)+ ε_t. B.1 Regularized shared-spectral surrogate The first three diagnostics in Appendix A suggest that using MO-LRN reliably in sequential optimization requires controlling its spectral capacity while preventing cross-objective transfer from depending on the overlap of independently fitted output spectra. We therefore separate nonstationary input-space structure, cross-objective dependence, and objective-specific smooth variation within a single composite surrogate. At BO iteration n, each objective is standardized using the observations available before selecting the next candidate: y¯n,i=1n∑j=1nyj,i,sn,i=max[1n∑j=1n(yj,i−y¯n,i)2]1/2,εY,y~j,i=yj,i−y¯n,isn,i, y_n,i= 1n _j=1^ny_j,i, s_n,i= \! \ [ 1n _j=1^n(y_j,i- y_n,i)^2 ]^1/2, _Y \, y_j,i= y_j,i- y_n,is_n,i, (6) where εY>0 _Y>0 prevents division by zero. The GP is fitted to ~n Y_n, while its posterior moments are transformed back to the original objective scale before acquisition construction. The spectral locations, scales, mixture weights, coregionalization, and residual parameters are learned by regularized marginal-likelihood optimization. Parameter bounds, penalty strengths, and the calibration and search settings specified below are fixed hypoerparameter configurations. Composite covariance. Equation equation 1 combines shared nonstationary input covariance, a common cross-objective correlation matrix for all mixtures, and independent smooth residuals. The shared scalar components follow the real four-term MO-LRN construction (Xu et al., 2026), with spectral-factor parameters shared across outputs. For each mixture q and coordinate d, let μq,d(a) _q,d^(a) and σq,d(a)>0 _q,d^(a)>0, a∈1,2a∈\1,2\, denote the Gaussian spectral-factor means and standard deviations. The index a distinguishes the two frequency arguments, not the outputs: the parameters are shared across objectives, but those for a=1a=1 and a=2a=2 need not coincide. Frequency pairs are independent across coordinates, with a within-pair correlation ζq _q shared across coordinates within each mixture. Define Cq,d=[(σq,d(1))2ζqσq,d(1)σq,d(2)ζqσq,d(1)σq,d(2)(σq,d(2))2],Vq,d(0)=12Cq,d.C_q,d= bmatrix( _q,d^(1))^2& _q _q,d^(1) _q,d^(2)\\ _q _q,d^(1) _q,d^(2)&( _q,d^(2))^2 bmatrix, V_q,d^(0)= 12C_q,d. (7) The spectral correlation is parameterized as ζq=(1−ϵζ)[2sigmoid(ζ~q)−1],ϵζ=10−6, _q=(1- _ζ) [2sigmoid( ζ_q)-1 ], _ζ=10^-6, which ensures |ζq|<1| _q|<1 and Cq,d≻0C_q,d 0. It is distinct from the coregionalization parameter r and the relative residual scale ρi _i defined below. Let q,d=(μq,d(1),μq,d(2)) μ_q,d=( _q,d^(1), _q,d^(2)) T. The normalized product of the two identical Gaussian spectral factors has mean q,d μ_q,d and covariance Vq,d(0)=Cq,d/2V_q,d^(0)=C_q,d/2. For numerical stability, determinants in the inverse calculations are bounded below by τ=10−8τ=10^-8. Writing adjadj for the 2×22× 2 adjugate, we define the effective spectral moments by Pq,d P_q,d =adj(Cq,d)maxdetCq,d,τ, = adj(C_q,d) \ C_q,d,τ\, Tq,d T_q,d =2Pq,d,Vq,d=adj(Tq,d)maxdetTq,d,τ, =2P_q,d, V_q,d= adj(T_q,d) \ T_q,d,τ\, q,d ν_q,d =Vq,d(2Pq,d)q,d. =V_q,d(2P_q,d) μ_q,d. (8) The determinant floor is applied separately to each factor covariance and to the summed precision. Under spectral sharing, the two stabilized factor precisions are identical, yielding Tq,d=2Pq,dT_q,d=2P_q,d. If neither floor is active, these expressions reduce to Vq,d=Cq,d/2V_q,d=C_q,d/2 and q,d=q,d ν_q,d= μ_q,d. The parameters σq,d(a) _q,d^(a) describe the factors before multiplication; the fourth-moment regularization in Eq. equation 18 acts on these pre-product parameters. Write Vqab=diagd=1,…,D([Vq,d]ab),q(a)=(νq,d(a))d=1D,=−′.V_q^ab=diag_d=1,…,D\! ([V_q,d]_ab ), ν_q^(a)= ( _q,d^(a) )_d=1^D, δ=x-x . Let A>0A>0 be the shared-kernel outputscale and let πq>0 _q>0 be the normalized mixture weights in Eq. equation 16, with ∑qπq=1 _q _q=1. The shared scalar component is Fq(,′) F_q(x,x ) =exp[−12(Vq11−2Vq12′+′Vq22′)] = \! [- 12 (x TV_q^11x-2x TV_q^12x +x TV_q^22x ) ] ×cos((q(1))−(q(2))′), × \! (( ν_q^(1)) Tx-( ν_q^(2)) Tx ), Gq(a)() G_q^(a)( δ) =exp[−12Vqaa]cos((q(a))),a∈1,2, = \! [- 12 δ TV_q^a δ ] \! (( ν_q^(a)) T δ ), a∈\1,2\, kqsLRN(,′) k_q^sLRN(x,x ) =Aπqgq4[Fq(,′)+Fq(′,)+Gq(1)()+Gq(2)()], = A _qg_q4 [F_q(x,x )+F_q(x ,x)+G_q^(1)( δ)+G_q^(2)( δ) ], (9) where gq>0g_q>0 is the amplitude coefficient defined below. The frequency parameters are angular frequencies, so no additional 2π2π factor appears in the phase terms. With determinant floors inactive and in the zero-overlap-jitter limit, the amplitude coefficient is gq(0)=(2π)−D/2∏d=1Ddet(2Cq,d)−1/2.g_q^(0)=(2π)^-D/2 _d=1^D (2C_q,d)^-1/2. This coefficient is (2π)D/2(2π)^D/2 times the overlap integral of two identical, unit-normalized 2D2D-variate Gaussian spectral factors. For the stabilized amplitude, we use ϵo=10−6 _o=10^-6 and the same determinant floor τ=10−8τ=10^-8: Sq,d S_q,d =2Cq,d+ϵoI2, =2C_q,d+ _oI_2, Mq,d M_q,d =adj(Sq,d)maxdetSq,d,τ, = adj(S_q,d) \ S_q,d,τ\, gq g_q =(2π)−D/2∏d=1D[maxdet(Mq,d+ϵoI2),τ]1/2. =(2π)^-D/2 _d=1^D [ \ (M_q,d+ _oI_2),τ\ ]^1/2. (10) The overlap jitter regularizes both the covariance sum and its stabilized inverse, with determinant floors applied at both stages. These operations affect only gqg_q; the effective moments Vq,dV_q,d and q,d ν_q,d remain those defined in Eq. equation 8, without overlap jitter. The determinant floors in the effective-moment calculation can, however, alter the covariance and mean entering the exponential and phase terms in Eq. equation 9. When the amplitude determinant floors are inactive, Eq. equation 10 becomes gq=(2π)−D/2∏d=1Ddet[(2Cq,d+ϵoI2)−1+ϵoI2]1/2,g_q=(2π)^-D/2 _d=1^D \! [(2C_q,d+ _oI_2)^-1+ _oI_2 ]^1/2, which converges to gq(0)g_q^(0) as ϵo→0 _o→ 0, provided the floors remain inactive. For the bi-objective problems considered here, all mixtures use the same unit-diagonal coregionalization matrix =[1r1],r=0.995tanh(r~).B= bmatrix1&r\\ r&1 bmatrix, r=0.995 ( r). (11) Since |r|<1|r|<1, B is positive definite. With αq=Aπq _q= A _q, the spectral-factor product contributes AπqA _q once through Eq. equation 9. The unit-diagonal matrix B controls cross-objective correlation. Each mixture’s variance contribution also depends on gqg_q; at the origin, kqsLRN(,)=Aπqgqk_q^sLRN(0,0)=A _qg_q. The positive semidefiniteness of the scalar component follows from a feature representation. The stabilized covariance Vq,dV_q,d is positive definite: Cq,d≻0C_q,d 0, and both adjugate divisions in Eq. equation 8 have positive denominators. Let the frequency pairs be distributed independently across coordinates as (ωq,d(1)ωq,d(2))∼2(q,d,Vq,d), pmatrix _q,d^(1)\\ _q,d^(2) pmatrix _2\! ( ν_q,d,V_q,d ), and write q(a)=(ωq,d(a))d=1D ω_q^(a)=( _q,d^(a))_d=1^D. Define ϕq()=12(cos((q(1)))+cos((q(2)))sin((q(1)))+sin((q(2)))). φ_q(x)= 12 pmatrix (( ω_q^(1)) Tx)+ (( ω_q^(2)) Tx)\\ (( ω_q^(1)) Tx)+ (( ω_q^(2)) Tx) pmatrix. Expanding the expected inner product gives kqsLRN(,′)=Aπqgq[ϕq()ϕq(′)].k_q^sLRN(x,x )=A _qg_q\,E\! [ φ_q(x) T φ_q(x ) ]. (12) Hence, kqsLRNk_q^sLRN is positive semidefinite for every admissible set of spectral parameters and every positive scalar amplitude gqg_q. The residual kernel for objective i is an ARD-RBF covariance, kiRBF(,′)=aiexp[−12∑d=1D(xd−xd′)2ℓi,d2],ai>0,ℓi,d>0.k_i^RBF(x,x )=a_i \! [- 12 _d=1^D (x_d-x_d )^2 _i,d^2 ], a_i>0, _i,d>0. (13) In the implementation, ai=ρivishareda_i= _iv_i^shared, where visharedv_i^shared is a reference marginal variance of the shared kernel and ρi _i is a learned relative residual scale. This prevents the RBF component from becoming numerically negligible merely because the marginal scale of the LRN kernel changes with the input dimension. For any finite input set, the contribution of component q has Gram matrix q⊗qsLRNB_q _q^sLRN, which is positive semidefinite because both factors are positive semidefinite. The residual kernels form a block-diagonal positive-semidefinite matrix. Equation equation 1, being a sum of these terms, therefore defines a valid matrix-valued covariance kernel. Spectral constraints and regularization. To prevent unstable high-frequency fits, the frequency locations and spectral scales are represented by bounded transformations: μq,d(a) _q,d^(a) =μmaxtanh(μ~q,d(a)),a∈1,2, = _ \! ( μ_q,d^(a) ), a∈\1,2\, (14) σq,d(a) _q,d^(a) =σmin+(σmax−σmin)sigmoid(σ~q,d(a)). = _ +( _ - _ )sigmoid\! ( σ_q,d^(a) ). (15) The normalized mixture weights are πq=πmin+(1−Qπmin)exp(π~q)∑q′=1Qexp(π~q′),∑q=1Qπq=1. _q= _ +(1-Q _ ) ( π_q) _q =1^Q ( π_q ), _q=1^Q _q=1. (16) The lower bound πmin _ prevents a component from disappearing completely during optimization. For a Gaussian spectral coordinate, define its fourth moment as m4(μ,σ)=μ4+6μ2σ2+3σ4.m_4(μ,σ)=μ^4+6μ^2σ^2+3σ^4. (17) We penalize the aggregate normalized fourth moment, ℛrough=λroughLspec4∑q=1Q∑d=1D∑a=12m4(μq,d(a),σq,d(a)),Lspec=maxμmax,σmax.R_rough= _roughL_spec^4 _q=1^Q _d=1^D _a=1^2m_4\! ( _q,d^(a), _q,d^(a) ), L_spec= \ _ , _ \. (18) This discourages spectral solutions that imply unnecessarily rough input-space functions. In addition to the bounded transformation in Eq. equation 14, we penalize frequencies near the limits of their allowed range: ℛbdry=λbdry2∑a=12meanq,d[−logmaxεbdry,1−(μq,d(a)μmax)2].R_bdry= _bdry2 _a=1^2mean_q,d [- \! \ _bdry,1- ( _q,d^(a) _ )^2 \ ]. (19) Mixture collapse is additionally discouraged by ℛmix=−λmix∑q=1Qlog(Qπq),R_mix=- _mix _q=1^Q (Q _q), (20) which is minimized at uniform mixture weights. The unconstrained coregionalization parameters receive a zero-centered Gaussian penalty, ℛB=12σB2∑q=1Qr~q2.R_B= 12 _B^2 _q=1^Q r_q^2. (21) For positive parameters, we use the generic log-scale penalty (,0,σz)=12‖log−log0σz‖22.P(z;z_0, _z)= 12 \| - _0 _z \|_2^2. (22) The residual, observation-noise, and shared-outputscale penalties are then ℛRBF _RBF =(,ρ0,σρ)+(ℓ,ℓ0,σℓ), =P( ρ; _01, _ρ)+P( ; _01, _ ), (23) ℛnoise _noise =(,η0,ση),ℛout=(A,A0,σA), =P( η; _01, _η), _out=P(A;A_0, _A), (24) where η denotes the objective-wise noise levels and A is the shared-kernel output scale. The composite exact GP is fitted by minimizing n(θ)=−logp(~n∣n,θ)+ℛrough+ℛbdry+ℛmix+ℛB+ℛRBF+ℛnoise+ℛout.J_n(θ)=- p( Y_n _n,θ)+R_rough+R_bdry+R_mix+R_B+R_RBF+R_noise+R_out. (25) The solution from the preceding BO iteration is used to initialize the next fit, while the objective standardization in Eq. equation 6 is recomputed from the currently available observations. In all experiments, Q=2Q=2, μmax=4 _ =4, (σmin,σmax)=(0.02,8)( _ , _ )=(0.02,8), πmin=0.05 _ =0.05, λrough=1 _rough=1, λmix=1 _mix=1, λbdry=0.5 _bdry=0.5, εbdry=10−6 _bdry=10^-6, and σB=1.5 _B=1.5. The RBF relative scale is initialized at ρ0=0.05 _0=0.05, with bounds [10−5,10][10^-5,10] and log-prior width σρ=0.75 _ρ=0.75. Its ARD lengthscales are initialized at ℓ0=0.5 _0=0.5, with bounds [0.02,10][0.02,10] and log-prior width σℓ=1.5 _ =1.5. B.2 Mean correction and covariance-scale adjustment The surrogate modifications above regularize the fitted covariance structure, but they do not ensure accurate uncertainty at the points preferred by the acquisition function. The fourth diagnostic in Appendix A shows systematic mean error and posterior underdispersion at acquisition-selected points. We address these errors through mean correction and covariance rescaling of the fitted GP posterior for acquisition decisions. Suppose the candidate selected at BO iteration t is tx_t. Before evaluating it, we record the raw GP predictive moments t=GP,t−1(t),t=GP,t−1(t), μ_t= μ_GP,t-1(x_t), _t= _GP,t-1(x_t), (26) together with the current output scale t−1s_t-1 and search mode StS_t. Only after observing ty_t do we append the resulting error to the calibration history. Thus, every correction used to select tx_t depends only on errors observed before iteration t. Hierarchical mean-bias correction. The standardized prequential residual is ~t=(t−t)⊘t−1, r_t=(y_t- μ_t) _t-1, (27) where ⊘ denotes element-wise division. Let nW_n contain the most recent WbW_b residuals available before iteration n+1n+1. For each objective, we compute the component-wise median and MAD, n=mediant∈n~t,n=mediant∈n|~t−n|,n=max1.4826n,τmin.m_n=median_t _n r_t, _n=median_t _n | r_t-m_n |, τ_n= \1.4826d_n, _ 1\. (28) The residuals are clipped component-wise according to ~tclip=clip(~t,n−κbn,n+κbn). r_t^\,clip=clip\! ( r_t,\,m_n- _b τ_n,\,m_n+ _b τ_n ). (29) Recent observations receive exponentially decaying weights wn,t=2−an,t/Hb,w_n,t=2^-a_n,t/H_b, (30) where an,t=0a_n,t=0 for the newest residual in nW_n, 11 for the second newest, and so forth, and HbH_b is the half-life. A shared bias estimate is first computed from all search modes: ^0,n=∑t∈nwn,t~tclipλ0+∑t∈nwn,t. b_0,n= _t _nw_n,t r_t^\,clip _0+ _t _nw_n,t. (31) For the mode s∈local,globals∈\local,global\ planned for the next iteration, the mode-specific estimate is partially pooled toward the shared estimate: ^s,n=∑t∈nSt=swn,t~tclip+λs^0,n∑t∈nSt=swn,t+λs. b_s,n= _ subarrayct _n\\ S_t=s subarrayw_n,t r_t^\,clip+ _s b_0,n _ subarrayct _n\\ S_t=s subarrayw_n,t+ _s. (32) When few global evaluations are available, this estimate remains close to the shared residual trend. As mode-specific evidence accumulates, the corresponding local or global correction receives greater weight. Fixed-hyperparameter spatial cross-validation. We next estimate covariance-scale error under spatial extrapolation. Let nK_n be the predictive training covariance, including observation noise, evaluated at the fitted hyperparameters, and let n=n−1,n=nvec(~n).P_n=K_n^-1, α_n=P_nvec( Y_n). (33) The leave-one-scalar-out Mahalanobis sum can be obtained directly from the precision matrix: HLOO=∑a=1nMαn,a2Pn,aa.H_LOO= _a=1^nM _n,a^2P_n,a. (34) To form spatial folds, we center the observed inputs and compute their first principal direction, 1∈argmax‖2=1n,cn,c.v_1∈ _\|v\|_2=1v TX_n,c TX_n,cv. (35) The observations are sorted by their projection n,c1X_n,cv_1 and divided into G contiguous folds. Let gG_g contain the stacked scalar-output indices associated with fold g. With the GP hyperparameters fixed at their full-data values, the conditional residual and covariance for this fold are g=n,gg−1n,g,g=n,gg−1.e_g=P_n,G_gG_g^-1 α_n,G_g, _g=P_n,G_gG_g^-1. (36) Consequently, Hsp=∑g=1Ggg−1g=∑g=1Gn,gn,gg−1n,g.H_sp= _g=1^Ge_g TV_g^-1e_g= _g=1^G α_n,G_g TP_n,G_gG_g^-1 α_n,G_g. (37) No hyperparameters are refitted for individual folds. Because the fitted hyperparameters depend on the held-out observations, these residuals provide internal scale diagnostics rather than an independent assessment of predictive coverage. For a Mahalanobis sum H containing NHN_H scalar residuals, define the shrunk covariance-scale estimate c(H,NH)=clip[cmin,cmaxCV](λCVc0+HλCV+NH).c(H,N_H)=clip_[c_ ,c_ ^CV] ( _CVc_0+H _CV+N_H ). (38) The LOO and spatial estimates are combined geometrically: cLOO=c(HLOO,nM),csp=c(Hsp,nM),cCV=cLOO 1−ωcspω.c_LOO=c(H_LOO,nM), c_sp=c(H_sp,nM), c_CV=c_LOO^\,1-ωc_sp^\,ω. (39) Acquisition-selected uncertainty error. For each previously evaluated acquisition-selected point, we calculate the raw normalized Mahalanobis error ht=1M(t−t)t−1(t−t).h_t= 1M(y_t- μ_t) T _t^-1(y_t- μ_t). (40) Under a correctly scaled M-variate predictive distribution, [ht]E[h_t] is approximately one. Let nU_n contain the most recent WcW_c valid values. We first define ut=logclip[1,cmaxsel](ht),ucap=minlogcmaxsel,mediant∈nut+1.4826κcMADt∈n(ut),u_t= _[1,c_ ^sel](h_t), u_cap= \! \ c_ ^sel,median_t _nu_t+1.4826 _cMAD_t _n(u_t) \, (41) and winsorize only the upper tail, utclip=minut,ucap.u_t^clip= \u_t,u_cap\. (42) Using weights vn,t=2−an,t/Hcv_n,t=2^-a_n,t/H_c, the selected-point scale is csel=clip[1,cmaxsel][exp(∑t∈nvn,tutclipλc+∑t∈nvn,t)].c_sel=clip_[1,c_ ^sel] [ \! ( _t _nv_n,tu_t^clip _c+ _t _nv_n,t ) ]. (43) The full predictive covariance scale is cpred=mincmaxsel,max1,cCV,csel.c_pred= \! \c_ ^sel, \1,c_CV,c_sel\ \. (44) Using cpredc_pred directly in qLogEHVI can overemphasize uncertain regions. The covariance scale exposed to the acquisition function is therefore tempered and capped: cacq=minccap,cpredγ.c_acq= \! \c_cap,c_pred^γ \. (45) This gives two posterior views: ppred(()∣n) p_pred(f(x) _n) =(GP(),cpredGP()), =N\! ( μ_GP(x),c_pred _GP(x) ), (46) pacq(()∣n,Sn+1=s) p_acq(f(x) _n,S_n+1=s) =(GP()+n⊙^s,n,cacqGP()). =N\! ( μ_GP(x)+s_n b_s,n,c_acq _GP(x) ). (47) The first is used for predictive diagnostics, whereas the second is the decision posterior supplied to qLogEHVI. We use Wb=30W_b=30, Hb=8H_b=8, λ0=5 _0=5, λs=3 _s=3, κb=2.5 _b=2.5, and τmin=0.10 _ =0.10 for mean correction. Spatial calibration uses G=5G=5, c0=1c_0=1, λCV=20 _CV=20, ω=0.75ω=0.75, cmin=0.25c_ =0.25, and cmaxCV=1000c_ ^CV=1000. Selected-point calibration uses Wc=12W_c=12, Hc=4H_c=4, λc=2 _c=2, κc=2.5 _c=2.5, and cmaxsel=64c_ ^sel=64. Finally, γ=1/4γ=1/4 and ccap=3c_cap=3. B.3 Pareto-local and periodic-global qLogEHVI The fifth diagnostic in Appendix A concerns the numerical optimization of the acquisition function. Full-domain multistart optimization can send several restarts toward similar distant or boundary-adjacent regions, where the nonstationary surrogate is weakly supported. We retain standard qLogEHVI but organize most of its restarts around diverse observed Pareto points. Periodic global rounds preserve exploration of the complete domain. Importantly, the GP itself remains global and is fitted to all observations; localization is applied only to acquisition optimization. Let ℐn=i∈1,…,n:∄j such that j dominates iI_n= \i∈\1,…,n\: j such that y_j dominates y_i \ (48) index the currently observed nondominated points. For each decision coordinate and objective, we apply min–max normalization over this set: x~i,d x_i,d =xi,d−minj∈ℐnxj,dmaxmaxj∈ℐnxj,d−minj∈ℐnxj,d,εz, = x_i,d- _j _nx_j,d \ _j _nx_j,d- _j _nx_j,d, _z\, (49) y~i,m y_i,m =yi,m−minj∈ℐnyj,mmaxmaxj∈ℐnyj,m−minj∈ℐnyj,m,εz. = y_i,m- _j _ny_j,m \ _j _ny_j,m- _j _ny_j,m, _z\. (50) The center-selection feature is i=[~i~i],i∈ℐn.z_i= bmatrix x_i\\ y_i bmatrix, i _n. (51) Let B denote the total acquisition-restart budget and KmaxK_ the maximum number of regions. The number of centers is Kn=minKmax,B,|ℐn|.K_n= \K_ ,B,|I_n|\. (52) If |ℐn|≤Kn|I_n|≤ K_n, all observed nondominated points are retained. Otherwise, the first center is the observed extreme maximizing the first objective, i1∈argmaxi∈ℐnyi,1,i_1∈ _i _ny_i,1, (53) and subsequent centers are selected by farthest-point traversal: ik∈argmaxi∈ℐn∖i1,…,ik−1minℓ<k1D+M‖i−iℓ‖22,k=2,…,Kn.i_k∈ _i _n \i_1,…,i_k-1\ _ <k 1D+M \|z_i-z_i_ \|_2^2, k=2,…,K_n. (54) Around center ikx_i_k, we define the axis-aligned region ℛn,k(rn)=[max¯,ik−rn(¯−¯),min¯,ik+rn(¯−¯)].R_n,k(r_n)= [ \ x,x_i_k-r_n( x- x)\,\; \ x,x_i_k+r_n( x- x)\ ]. (55) Thus, rnr_n is a half-width expressed as a fraction of the full range of each decision coordinate. The B restarts are allocated as evenly as possible: Bn,k=⌊BKn⌋+k≤BmodKn,∑k=1KnBn,k=B.B_n,k= BK_n +I \k≤ B K_n \, _k=1^K_nB_n,k=B. (56) For region k, the raw local initialization pool has size Nn,kraw=max16Bn,k,⌈NrawBn,kB⌉.N_n,k^raw= \! \16B_n,k, N_raw B_n,kB \. (57) Its points are sampled around the center and clipped to the region: n,k,j=Πℛn,k(rn)[ik+σn(¯−¯)⊙j],j∼(,),σn=minσ0,rn/2.u_n,k,j= _R_n,k(r_n) [x_i_k+ _n( x- x) ξ_j ], ξ_j (0,I), _n= \ _0,r_n/2\. (58) The Bn,kB_n,k pool points with the largest acquisition values initialize bounded gradient optimization inside ℛn,k(rn)R_n,k(r_n). Let u=1,…,Tu=1,…,T index the BO evaluations after an initial design of size n0n_0. Before iteration u, n=n0+u−1n=n_0+u-1 evaluations have been completed. The search mode for the next candidate n+1x_n+1 is Sn+1=global,u≤40andumod5=0,global,u>40andumod10=0,local,otherwise.S_n+1= casesglobal,&u≤ 40\ and\ u 5=0,\\ global,&u>40\ and\ u 10=0,\\ local,&otherwise. cases (59) Hence, global rounds constitute 20%20\% of the first 40 post-initial-design iterations and 10%10\% thereafter. For each mode s∈local,globals∈\local,global\, define the acquisition function using the corresponding decision posterior in Eq. equation 4: αn(,s)=qLogEHVI(;pacq(()∣n,Sn+1=s)). _n(x;s)=qLogEHVI (x;p_acq(f(x) _n,S_n+1=s) ). (60) On a local round, let n,k,b⋆x_n,k,b denote the result of restart b optimized within region k. We select the best region–restart pair and return its candidate: (k⋆,b⋆) (k ,b ) ∈argmax1≤k≤Kn1≤b≤Bn,kαn(n,k,b⋆,local), ∈ _ subarrayc1≤ k≤ K_n\\ 1≤ b≤ B_n,k subarray _n(x_n,k,b ;local), (61) n+1 _n+1 =n,k⋆,b⋆. =x_n,k ,b . On a global round, a scrambled Sobol pool of size Nglobalraw=maxNraw,16BN_global^raw= \N_raw,16B\ (62) is generated over X. Its B highest-valued points initialize full-domain optimization. Let n,b⋆∈x_n,b denote the result of restart b. We select the best restart and return its candidate: b⋆ b ∈argmax1≤b≤Bαn(n,b⋆,global), ∈ _1≤ b≤ B _n(x_n,b ;global), (63) n+1 _n+1 =n,b⋆. =x_n,b . Adaptive local radius. After evaluating a local candidate n+1x_n+1, define the realized hypervolume improvement ΔHVn+1=HV(n+1)−HV(n). _n+1=HV(D_n+1)-HV(D_n). (64) A local round is considered successful if ΔHVn+1>ϵHV _n+1> _HV. On local rounds, the common radius is updated as rn+1=minrmax,ηincrn,after three consecutive local successes,maxrmin,ηdecrn,after five consecutive local failures,rn,otherwise.r_n+1= cases \r_ , _incr_n\,&after three consecutive local successes,\\ \r_ , _decr_n\,&after five consecutive local failures,\\ r_n,&otherwise. cases (65) The corresponding counter is reset whenever the other outcome occurs and after a radius expansion or contraction. Global rounds satisfy rn+1=rnr_n+1=r_n and leave both local counters unchanged. We use B=8B=8, Nraw=128N_raw=128, Kmax=6K_ =6, rn0=0.20r_n_0=0.20, rmin=0.10r_ =0.10, rmax=0.25r_ =0.25, σ0=0.10 _0=0.10, ϵHV=10−4 _HV=10^-4, ηinc=1.2 _inc=1.2, and ηdec=0.8 _dec=0.8. Computational cost. With n complete observations of M objectives, dense exact inference uses an nM×nMnM× nM covariance matrix. A Cholesky factorization requires ((nM)3)O((nM)^3) time and ((nM)2)O((nM)^2) storage, with repeated factorizations during hyperparameter optimization. Fixed-hyperparameter calibration uses the fitted covariance and precision submatrices without refitting the GP for each fold. Sequential fits are warm-started, and local regions share a fixed acquisition-restart budget. Localization changes the acquisition domain but not the GP training set. The reported benchmarks contain at most 110 observations of two objectives, corresponding to 220 scalar observations. Appendix C Additional performance and robustness analyses Run-to-run robustness and alternative Pareto-quality metrics. We first complement the mean-performance results with adverse-tail statistics computed across seeds within each problem. For nHV, the adverse tail is its 10th percentile, whereas for nIGD and the two AUC measures it is the 90th percentile. Table 4 reports the corresponding mean problem ranks. MOLRN-BO ranks first on all four tail measures and also has the smallest mean across-seed standard deviation of final nIGD among the nine algorithms (0.03580.0358). These results indicate that its main robustness benefit is a reduced incidence of poorly performing trials, rather than a uniform reduction in variance on every problem and metric. We additionally evaluate two distance-based indicators in the ideal–nadir normalized objective space. For an approximation set A and reference front R, normalized IGD+ is defined as nIGD+(A,R)=1|R|∑∈Rmin∈A‖max(−,)‖2,nIGD^+(A,R)= 1|R| _r∈ R _a∈ A \| (r-a,0) \|_2, (66) and the normalized additive ϵε-indicator is Iϵ+(A,R)=max∈Rmin∈Amaxj(rj−aj).I_ε+(A,R)= _r∈ R _a∈ A _j(r_j-a_j). (67) Both indicators are lower-is-better. MOLRN-BO obtains the best mean problem rank for nIGD+ (2.252.25) and lies in the top three on 11/1211/12 problems. It outperforms the original MO-LRN on all 12 problems under this metric. For the additive ϵε-indicator, qLogNParEGO ranks first (2.832.83) and MOLRN-BO ranks second (3.173.17). MOLRN-BO is better than qLogNParEGO on only 5/125/12 problems, although it remains better than the original MO-LRN on 10/1210/12. This difference is consistent with the indicators’ emphases: nIGD+ averages directional shortfalls over the reference front, whereas Iϵ+I_ε+ is controlled by the most poorly covered reference point. Thus, MOLRN-BO provides strong average front coverage, while qLogNParEGO is slightly stronger in worst-location coverage. Table 4: Mean problem ranks for adverse-tail and additional Pareto-quality metrics over 12 benchmarks using all valid trials. Lower is better. Tail quantiles are computed across seeds within each problem. Bold and underlined entries denote the best and second-best ranks, respectively. Method nHV-q10 nIGD-q90 HV-regret-AUC-q90 nIGD-AUC-q90 Final nIGD+ Final Iϵ+I_ε+ MOLRN-BO 2.08 2.08 2.00 2.83 2.25 3.17 qLogNParEGO 2.83 3.50 3.00 3.33 2.50 2.83 qLogEHVI 4.08 4.83 3.58 3.92 4.00 4.67 mobopt 5.12 4.33 5.54 4.67 4.58 3.67 JES 5.21 4.58 5.21 4.25 4.75 4.42 Original MO-LRN 6.12 5.08 6.29 5.75 5.83 5.33 NSGA-I 5.62 6.50 6.29 7.17 6.33 7.50 MESMO 7.21 7.00 7.12 6.50 7.00 5.75 Sobol 6.71 7.08 5.96 6.58 7.75 7.67 Convergence behavior. Figure 2 reports the nHV and nIGD trajectories over the complete evaluation budget. MOLRN-BO’s advantage is generally more pronounced late in optimization. Across the nine algorithms, its mean nIGD problem rank improves from 3.333.33 after 25 post-initial evaluations to 2.672.67, 2.422.42, and 2.172.17 after 50, 75, and 100 evaluations, respectively. Its nHV rank improves more gradually from 2.582.58 at 25 evaluations to 2.252.25 at the final budget. This late improvement explains why MOLRN-BO achieves the best final nHV and nIGD ranks while qLogNParEGO remains slightly better on the two anytime AUC measures. The trajectories also reveal meaningful problem dependence. MOLRN-BO is comparatively slow on WFG1, performs weakly in both final nHV and anytime performance on WFG5, and exhibits an extended nIGD plateau on the disconnected DTLZ7 front. In contrast, its final coverage improves substantially on many of the remaining ten-dimensional problems and is consistently better than that of the original MO-LRN. The curves therefore support a conclusion of strong final quality and reduced poor-trial risk, but not uniformly faster convergence over the entire benchmark suite.