qEHVI: Batch Multi-Objective Optimization

• qEHVI is a batch acquisition function that extends single-point EHVI, enabling simultaneous evaluations in multi-objective Bayesian optimization.
• It leverages Gaussian process surrogates with Monte Carlo and quasi-Monte Carlo integration for efficient gradient computation and acquisition maximization.
• The approach has demonstrated robust performance in real-world applications like high-throughput materials design and automated model selection.

The q-Expected Hypervolume Improvement (qEHVI) acquisition function is a central tool in parallel multi-objective Bayesian optimization. It quantifies the expected gain in hypervolume—measured with respect to a reference point—when evaluating a batch of $q$ candidate input locations, enabling efficient exploration and exploitation across multiple competing black-box objectives. qEHVI generalizes the single-point Expected Hypervolume Improvement (EHVI) criterion, extending its applicability to batch or parallel evaluation regimes. This approach underlies modern high-throughput experimental design, automated machine learning, and constrained optimization in both discrete and continuous domains. qEHVI admits highly efficient implementations using Gaussian process surrogates, Monte Carlo or quasi-Monte Carlo integration, and automatic differentiation, and supports both noisy and noise-free observed objectives.

1. Mathematical Definition and Probabilistic Foundation

Let $f:\mathbb{X} \rightarrow\mathbb{R}^M$ be a black-box vector-valued objective function, where $\mathbb{X}\subseteq\mathbb{R}^d$ and $M$ is the number of objectives. The current Pareto set $\mathcal{P}\subseteq\mathbb{R}^M$ is the set of nondominated vectors among observed outcomes, and $r\in\mathbb{R}^M$ is a user-specified reference point.

For a batch of $q$ candidates $X = \{x_1,\ldots,x_q\}$ and their unknown objectives $Y = \{f(x_1),\ldots,f(x_q)\}$, the hypervolume improvement (HVI) relative to $\mathcal{P}$ is

$$\mathrm{HVI}(Y; \mathcal{P}, r) = \mathrm{HV}(\mathcal{P} \cup Y; r) - \mathrm{HV}(\mathcal{P}; r),$$

where $\mathrm{HV}(\cdot; r)$ denotes the Lebesgue measure of the dominated region up to $r$.

Assuming a probabilistic model (commonly independent-output Gaussian processes) fitted to observed data, $f(X)$ has joint predictive density $p(Y| X, \mathcal{D})$. The q-Expected Hypervolume Improvement is then defined as

$$\mathrm{qEHVI}(X) = \mathbb{E}_{Y \sim p(\cdot| X, \mathcal{D})} \left[ \mathrm{HVI}(Y; \mathcal{P}, r) \right] = \int \mathrm{HVI}(Y; \mathcal{P}, r) \, p(Y | X, \mathcal{D})\, dY.$$

No general closed form exists for $q>1$ or $M>2$, so $\mathrm{qEHVI}(X)$ is approximated via (quasi-)Monte Carlo.

2. Computational Strategies and Efficient Gradient Computation

Efficient qEHVI computation critically leverages reparameterization for differentiability. Monte Carlo integration proceeds by sampling $S$ draws $Y^{(s)}$ from $p(Y| X, \mathcal{D})$: $$\mathrm{qEHVI}(X) \approx \frac{1}{S} \sum_{s=1}^S \mathrm{HVI}(Y^{(s)}; \mathcal{P}, r),$$ with $Y^{(s)} = \mu(X) + L(X)z^{(s)}$, where $z^{(s)} \sim \mathcal{N}(0, I)$ and $L(X)$ is the Cholesky factor of the predictive covariance $\Sigma(X)$. Inclusion–exclusion combinatorics decompose HVI across all non-empty subsets of the batch, scaling as $O(2^q)$ per sample.

Automatic differentiation is enabled by fixing base draws $z^{(s)}$ throughout acquisition optimization (sample average approximation), so gradients $\nabla_X \mathrm{qEHVI}(X)$ may be computed exactly via modern autodiff frameworks. This allows the use of first-order (L-BFGS-B, Adam) or quasi-second-order optimizers and supports efficient backpropagation through the entire acquisition process, including the Cholesky decomposition and batched kernel operations (Daulton et al., 2020, Daulton et al., 2021, Chen et al., 10 Dec 2025).

3. Extension to Noisy and Constrained Objective Settings

Noisy observation settings are addressed by integrating over the posterior of possible in-sample Pareto frontiers, leading to the NEHVI and qNEHVI acquisition functions: $$\alpha_{\mathrm{qNEHVI}}(X) = \int \mathrm{qEHVI}(X | \mathcal{P}_n(f)) \, p(f(X_n) | \mathcal{D}_n) \, df(X_n),$$ where $\mathcal{P}_n(f)$ is the Pareto frontier of a sample of the latent (denoised) outcomes at previous query points $X_n$. This double integration is again handled by Monte Carlo with outer samples over the past evaluations and inner samples conditional on them (Daulton et al., 2021).

Constraints are accommodated by outcome-based feasibility weighting, multiplying each sample HVI by an indicator or smooth approximation of the feasibility condition, maintaining end-to-end differentiability (Daulton et al., 2020).

4. Numerical Pathologies and Robust Variants

Numerical pathologies are inherent in qEHVI, including vanishing acquisition values and gradients over most of the search space, especially when candidate means are far from the observed Pareto front. Piecewise-constant regions, caused by the max operator, result in zero gradients for large swathes of input space and complicate acquisition maximization (Ament et al., 2023).

Recent work introduces log-space smoothed approximations ("Log-qEHVI") employing softplus (fat-Softplus) to replace kinked [·]_+ operations and averaging in logarithmic space using numerically stable logsumexp routines. This regularization enables robust acquisition optimization and improves convergence rates, without altering theoretical optima (Ament et al., 2023).

5. Acquisition Maximization and Algorithmic Best Practices

Acquisition maximization poses a high-dimensional, non-convex optimization challenge. Application of gradient-based optimizers to the differentiable Monte Carlo approximation is standard for moderate $q$ and low to moderate search dimension $d$. However, due to multi-modality and vanishing gradients, metaheuristic strategies such as simulated annealing (SA) have been proposed as more robust alternatives—especially in high-dimensional or multi-modal landscapes. Empirical results show that SA achieves superior Pareto front coverage and hypervolume on challenging problems compared to SLSQP, particularly for large $q$ and complex objective spaces (Alvi et al., 12 Jan 2026).

The algorithmic loop is succinctly:

1. Fit GP surrogates to data.
2. For each batch $X$, compute $\mathrm{qEHVI}(X)$ via MC/QMC with cached box decomposition for HVI.
3. Maximize acquisition (gradient-based or SA).
4. Evaluate batch, update data, repeat.

6. Comparative Performance and Application Domains

qEHVI provides strong performance in both synthetic and real-world tasks, rapidly attaining near-optimal hypervolume in alloy design (1–3 objectives) and outperforming alternative acquisition functions (parEGO, qParEGO, TSEMO, etc.) in terms of both efficiency and final solution quality (Mamun et al., 2024, Daulton et al., 2021). NEHVI and qNEHVI demonstrably maintain robust Pareto set spacing and noise resilience, as shown on Branin–Currin, DTLZ2, adaptive bitrate control, vehicle design, and LLM Pareto set selection benchmarks (Daulton et al., 2021, Chen et al., 10 Dec 2025). The qEHVI framework is central to modern applications in high-throughput materials design, automated model selection, and multi-objective optimization for AI systems.

7. Computational Complexity and Implementation Details

The computational cost of qEHVI is dominated by the combination of Monte Carlo sample count $S$, Pareto front size $|\mathcal{P}|$, batch size $q$, and number of objectives $M$. The complexity per batch is generally $O(S \cdot 2^q \cdot M \cdot |\mathcal{P}|)$ using the inclusion–exclusion principle, but polynomial-time acceleration (e.g., cached box decomposition, QMC) ameliorates exponential scaling for moderate $q$ (Daulton et al., 2021, Daulton et al., 2020).

Efficient implementation is available in open-source frameworks such as BoTorch, leveraging CPU/GPU acceleration. Surrogate models are routinely GPs (independent outputs or joint), with training via MAP/maximum marginal likelihood, or SAASGP with HMC for high-dimensional search. qEHVI supports both continuous and discrete (batch) candidate selection, with adaptability to evolutionary or block-wise model selection pipelines (Mamun et al., 2024, Chen et al., 10 Dec 2025).

---

Summary Table: qEHVI Variants and Key Attributes

Variant	Supports Noise/Constraints?	Optimization Approach	Typical Complexity
qEHVI	No (deterministic only)	Gradient-based/SA	$O(S 2^q M |\mathcal{P}|)$
qNEHVI	Yes (noise, constraints)	Gradient-based/SA	Polynomial in $q$ (CBD)
Log-qEHVI	No/Yes (numerical fix)	Gradient-based	As above

---

The qEHVI framework, and its extensions to noise, constraints, and robust log-space variants, represent state-of-the-art methodology for batch multi-objective Bayesian optimization in both theory and empirical performance (Daulton et al., 2021, Ament et al., 2023, Daulton et al., 2020, Alvi et al., 12 Jan 2026, Chen et al., 10 Dec 2025, Mamun et al., 2024).