arXiv CSopen access

Beyond Expected Information Gain: Stable Bayesian Optimal Experimental Design with Integral Probability Metrics and Plug-and-Play Extensions

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

machine learning, deep learning, neural networks

Beyond EIG: Stable BOED with IPMs

arXiv:2604.21849v1 [stat.ML] 23 Apr 2026

Di Wu∗

[email protected]

Department of Mathematics University of Maryland, College Park MD 20742, USA

Ling Liang∗

[email protected]

Department of Mathematics The University of Tennessee, Knoxville TN 37916, USA

Haizhao Yang†

[email protected]

Department of Mathematics and Department of Computer Science University of Maryland, College Park MD 20742, USA

Abstract Bayesian Optimal Experimental Design (BOED) provides a rigorous framework for decisionmaking tasks in which data acquisition is often the critical bottleneck, especially in resourceconstrained settings. Traditionally, BOED typically selects designs by maximizing expected information gain (EIG), commonly defined through the Kullback-Leibler (KL) divergence. However, classical evaluation of EIG often involves challenging nested expectations, and even advanced variational methods leave the underlying log-density-ratio objective unchanged. As a result, support mismatch, tail underestimation, and rare-event sensitivity remain intrinsic concerns for KL-based BOED. To address these fundamental bottlenecks, we introduce an IPM-based BOED framework that replaces density-based divergences with integral probability metrics (IPMs), including the Wasserstein distance, Maximum Mean Discrepancy, and Energy Distance, resulting in a highly flexible plug-and-play BOED framework. We establish theoretical guarantees showing that IPM-based utilities provide stronger geometry-aware stability under surrogate-model error and prior misspecification than classical EIG-based utilities. We also validate the proposed framework empirically, demonstrating that IPM-based designs yield highly concentrated credible sets. Furthermore, by extending the same sample-based BOED template in a plug-and-play manner to geometry-aware discrepancies beyond the IPM class, illustrated by a neural optimal transport estimator, we achieve accurate optimal designs in high-dimensional settings where conventional nested Monte Carlo estimators and advanced variational methods fail. Keywords: Bayesian Optimal Experimental Design, Integral Probability Metrics, Neural Optimal Transport

∗. Equal contribution. †. Corresponding author.

Wu, Liang and Yang

1 Introduction Optimal Experimental Design (OED) provides a principled framework for selecting experiments to learn unknown model parameters as efficiently as possible (Steinberg and Hunter, 1984). It has been widely used across science and engineering, including parameter inference in scientific computing, such as partial differential equations (Huan and Marzouk, 2013; Hellmuth et al., 2025), and inverse problems (Haber et al., 2008; Alexanderian et al., 2016; Ruthotto et al., 2018; Alexanderian, 2021; Helin et al., 2022; Wu et al., 2023; Liang and Yang, 2024; Chang et al., 2025; Chen et al., 2026). In many of these applications, data acquisition is a major bottleneck because experiments can be expensive, simulations may require substantial computational resources, and only a limited number of measurements can be collected. This challenge arises in a wide range of settings including materials discovery (Lookman et al., 2019; Chen et al., 2021), calibration of expensive computational simulators (Sacks et al., 1989; Dunbar et al., 2022; Sun et al., 2024), and data selection and training in modern machine learning (Chen et al., 2015; Yang et al., 2016; Fiez et al., 2019; Camilleri et al., 2021; Wagenmaker and Jamieson, 2022; Bhatt et al., 2024; Hu et al., 2025; Zhu et al., 2025; Gao et al., 2026). In such resource-constrained regimes, one must carefully decide which experiments to perform in order to maximize the information gained about the underlying system. In this context, the use of prior knowledge plays an essential role. Bayesian Optimal Experimental Design (BOED) extends OED by explicitly incorporating prior uncertainty about the unknown parameters (Lindley, 1956). Let x ∈ X denote the unknown model parameters, Ξ the design domain, ξ ∈ Ξ an experimental design, and y ∈ Y the corresponding experimental outcome with Y being the space of all possible outcomes. Given a prior distribution p(x), classical BOED seeks a design ξ ∗ that maximizes the expected information gain (EIG), commonly defined through the Kullback–Leibler (KL) divergence from the posterior to the prior (Shannon, 1948): Z Z p(x | y, ξ) ∗ ξ := argmax EIG(ξ) := log p(x | y, ξ) p(y | ξ) dx dy . p(x) ξ∈Ξ Y X Despite its conceptual appeal, the classical EIG formulation raises significant challenges. From a computational perspective, evaluating EIG requires estimating nested expectations, which is typically expensive and difficult to do accurately (Rainforth et al., 2018), and standard Monte Carlo methods often suffer from high variance and slow convergence. From a modeling perspective, the main issue is structural rather than estimator-specific. Even when advanced estimators (Foster et al., 2019) reduce the computational burden of EIG evaluation, the utility remains a log-density-ratio functional. Consequently, it is inherently sensitive to tail underestimation, near-zero likelihood regions, and rare but highly informative observations. Moreover, modern BOED pipelines often replace continuous priors or posteriors with empirical approximations, and the KL divergence is not always a natural perturbation metric between such exact and approximate measures because the required absolute continuity may fail. To address these limitations, we propose a geometry-aware reformulation of BOED based on Integral Probability Metrics (IPMs) (Müller, 1997). IPMs, including the Wasserstein distance (Dudley, 2018), Maximum Mean Discrepancy (MMD) (Gretton et al., 2009; Berlinet and Thomas-Agnan, 2011; Gretton et al., 2012), and Energy Distance (Székely 2

Beyond EIG: Stable BOED with IPMs

et al., 2004; Székely and Rizzo, 2005), quantify the discrepancy between two distributions through their action on a prescribed class of test functions. We introduce an IPM-based Bayesian Optimal Experimental Design framework by defining the design utility as Z γF p(x), p(x | y, ξ) p(y | ξ)dy, max UF (ξ) := ξ∈Ξ

where γF is the IPM induced by a function class F. This formulation is naturally flexible such that the same sample-based BOED template can also be extended in a plug-and-play manner to other discrepancies based on domain knowledge. This IPM-based framework offers both theoretical and computational advantages over the traditional EIG-based formulation. On the theoretical side, we conduct a systematic stability analysis to characterize how surrogate-model and prior errors affect BOED utilities. In particular, we show that bounded-kernel MMD, as a bounded IPM, admits universal stability bounds controlled directly by the surrogate likelihood error or the prior perturbation, without requiring strict tail assumptions. We further analyze the unbounded IPMs considered in this paper, namely Energy Distance and the 1-Wasserstein distance, and clarify the tail and smoothness conditions required for stability. By contrast, the classical KL divergence does not admit comparable robustness guarantees because of its dependence on density ratios. On the computational side, IPMs enable stable sample-based utility estimation without explicit density-ratio evaluation, leading to substantial efficiency gains. Moreover, in our experiments, IPM-based utilities exhibit broader near-optimal design regions and more stable optimization landscapes than KL-based criteria. More generally, BOED can be built from any discrepancy between the prior and posterior. In this work, our theoretical development focuses on IPM utilities because they admit a unified stability analysis. At the same time, the resulting sample-based BOED pipeline is modular, which allows us to incorporate geometry-aware discrepancies beyond the IPM class at the computational level. We illustrate this extension by incorporating a neural optimal transport estimator (Amos et al., 2017) in high-dimensional settings. Contributions.

The main contributions of this work are summarized as follows:

• A stable BOED framework based on Integral Probability Metrics. We develop a BOED formulation in which the design utility is defined through an Integral Probability Metric between the prior and posterior distributions. This provides a geometry-aware alternative to classical KL-based expected information gain and avoids direct dependence on density ratios. Within this framework, we establish utility perturbation bounds with respect to both surrogate likelihood error and prior misspecification, and we make explicit the different assumptions required for bounded and unbounded IPMs. • Reliable sample-based utility estimation and improved optimization behavior. We show empirically that IPM-based utilities admit stable sample-based estimation without explicit density estimation or log-ratio evaluation. In our experiments, this leads to substantial computational advantages over KL-based criteria, while also producing broader near-optimal design regions and more reliable optimization landscapes. 3

Wu, Liang and Yang

• Plug-and-play extensions beyond the IPM class. Beyond the IPM theory developed in this paper, we show that the same sample-based BOED pipeline is modular enough to incorporate other geometry-aware discrepancies. As a representative example, we integrate a neural optimal transport estimator for the squared 2-Wasserstein distance in high-dimensional settings where nested Monte Carlo approaches and variational methods become computationally prohibitive or unreliable, illustrating that the proposed computational framework extends naturally beyond the IPM class. The rest of this paper is organized as follows. Section 2 reviews the relevant literature. Section 3 provides the formal definition of IPMs, together with numerical results that motivate our approach. In Section 4, we establish the theoretical guarantees for the proposed framework and provide a detailed discussion. Section 5 presents a series of experiments that validate the robustness and demonstrate the scalability and flexibility of the plug-andplay framework. Finally, Section 6 concludes the paper with a summary and directions for future research. All formal proofs and supplementary mathematical details are deferred to the appendix.

2 Related Work Classical Bayesian Optimal Experimental Design and Expected Information Gain. Classical BOED relies on an information-theoretic strategy, where Shannon entropy (Shannon, 1948) forms the basis for the expected information gain (EIG) (Lindley, 1956; Bernardo, 1979; Shewry and Wynn, 1987; Sebastiani and Wynn, 2000). Expected Fisher information has also been studied as an alternative optimality criterion (Walker, 2016; Overstall, 2022; Prangle et al., 2023). Despite its elegant theoretical foundation, estimating the EIG is difficult in practice because it involves nested expectations, which can be computationally expensive and statistically unstable. To overcome this computational bottleneck, various estimation strategies have been developed, including nested Laplace approximations (Long et al., 2013; Beck et al., 2018) and nested Monte Carlo estimators (Myung et al., 2013; Rainforth et al., 2018). To further improve efficiency, recent variational approaches utilize function approximation to relax the computational burden imposed by strict nesting (Barber and Agakov, 2004; Huan and Marzouk, 2013; Kleinegesse and Gutmann, 2019; Foster et al., 2019, 2021; Jin et al., 2024; Ohn and Lin, 2024). For broader discussions of the development of BOED, we refer the reader to several surveys (Chaloner and Verdinelli, 1995; Ryan et al., 2016; Rainforth et al., 2024; Huan et al., 2024). Our goal is not to claim that KL-based BOED is always computationally infeasible. Rather, we stress that estimator-level improvements do not change the underlying stability properties of the KL objective. In contrast, our work changes the utility itself and studies stability at the level of probability measures. Optimal Transport and Wasserstein Information Criteria. While KL divergence is widely used in BOED, it is fragile when dealing with non-overlapping supports (Arjovsky et al., 2017). To address this, recent concurrent works (Helin et al., 2025; Kerrigan et al., 2025) have explored using optimal transport (Peyré and Cuturi, 2019; Cang et al., 2022; Wu et al., 2025) and the Wasserstein distance (Villani et al., 2009; Liu et al., 2021; Song et al., 2026) as measures between probability distributions, which is often more natural 4

Beyond EIG: Stable BOED with IPMs

when supports are poorly aligned or when empirical particle approximations are used. This line of work provides an important step toward geometry-aware BOED. Our work is closely related in spirit, but differs in scope: rather than focusing on a single transport criterion, we develop a unified BOED view based on Integral Probability Metrics, which allows us to compare several geometry-aware discrepancies within one theoretical framework. IPMs and Neural Discrepancy Estimators. Integral Probability Metrics (IPMs) (Müller, 1997) have gained significant attention in machine learning as robust, geometryaware alternatives to density-based divergences. Recent studies (Amos et al., 2017; Korotin et al., 2019; Makkuva et al., 2020; Birrell et al., 2022) also explore learning-based discrepancy estimation. These methods provide scalable sample-based tools for comparing probability measures in high-dimensional settings by parameterizing transport potentials, maps, or dual objectives with neural networks. Our use of OT-ICNN in the high-dimensional experiments should be understood in this computational sense. The theory developed in this paper focuses on IPM-based BOED utilities, while the sample-based BOED pipeline is modular enough to accommodate geometry-aware discrepancies beyond the IPM class.

3 Integral Probability Metrics for BOED To construct a stable utility function for BOED, we replace the standard Kullback-Leibler (KL) divergence with Integral Probability Metrics (IPMs). In this section, we formally define the IPM framework, introduce the three IPMs that form the theoretical core of the paper, and illustrate the fundamental fragility of KL divergence in rare-event regimes. 3.1 Integral Probability Metrics Integral probability metrics are a family of distances between two probability distributions defined by the maximum difference in their expectations over a specified class of test functions. Let P and Q be probability measures on the measurable parameter space X , and let F be a set of measurable real-valued functions f : X → R. The IPM induced by the function class F is defined as: γF (P, Q) = sup EX∼P [f (X)] − EY ∼Q [f (Y )] . f ∈F

By restricting the complexity of the test function class F , IPMs offer a geometry-aware measure of distance that does not require the distributions to have overlapping support or well-behaved density ratios. In this work, we focus on three representative IPMs characterized by distinct function classes: 1. 1-Wasserstein Distance (W1 ): By the Kantorovich–Rubinstein duality theory (Hanin, 1992), the 1-Wasserstein distance is defined as W1 (P, Q) = sup

EX∼P [f (X)] − EY ∼Q [f (Y )]

∥f ∥L ≤1

where ∥ · ∥L denotes the Lipschitz constant. Thus, W1 is the IPM induced by the class of 1-Lipschitz functions; 5

Wu, Liang and Yang

2. Maximum Mean Discrepancy (MMD): Let Hk be a reproducing kernel Hilbert space (RKHS) (Smola et al., 2007) with a kernel k. The corresponding MMD is MMDk (P, Q) =

EX∼P [f (X)] − EY ∼Q [f (Y )] .

sup ∥f ∥Hk ≤1

Hence, MMD is the IPM induced by the unit ball of the RKHS Hk ; 3. Energy Distance (ED): Let ρ : X × X → R be a semimetric of negative type. The associated energy distance is 1/2 EDρ (P, Q) = 2Eρ(X, Y ) − Eρ(X, X ′ ) − Eρ(Y, Y ′ ) , where X, X ′ ∼ P and Y, Y ′ ∼ Q are independent, respectively. Asshown by Sejdinovic et al. (2013), if we set kx0 (x, y) = 21 ρ(x, x0 ) + ρ(y, x0 ) − ρ(x, y) for a fixed x0 ∈ X , then kx0 is a positive definite kernel and √ EDρ (P, Q) = 2 MMDkx0 (P, Q). Therefore, energy distance is equivalent to an IPM induced by the unit ball of the RKHS associated with kx0 . While Energy Distance (ED) can be viewed as a special case of MMD, the two exhibit significant differences in both theory and practice. As we will discuss in Section 4, the standard configurations including MMD with bounded kernels and ED with Euclidean distance lead to distinct theoretical guarantees. Practically, MMD offers greater flexibility through various kernel choices, whereas ED is often more robust. Note that the high-dimensional experiments in Section 5.3 should be viewed as plug-and-play extensions beyond the IPM class and are not covered by the theory. 3.2 Motivating Example: The Fragility of KL Divergence To demonstrate where KL-based measures typically fail in practice, we compare these sample-based estimators using a controlled rare-event contamination setup. Consider two distributions that differ only by a low-probability component: Q ∼ N (0, 1),

P = (1 − ϵ)N (0, 1) + ϵN (10, 0.12 ),

with ϵ = 10−2 .

Here, the difference between P and Q comes down to tail events, which rarely appear consistently in small sample batches. The bottom row of Figure 1 shows that all discrepancies respond to the presence of contaminated samples, but they do so in qualitatively different ways. The KL estimate is almost locked to the exact contamination count, producing an approximately linear response and hence large trial-to-trial variance. By contrast, the IPM estimates are influenced by the contaminated samples through pairwise geometry or transport cost, which leads to a smoother and less count-dominated dependence. In particular, ED and W1 show a clear but gradual increase with ncontam , while MMD exhibits a noisier finite-sample response due to kernel interactions. The key point is therefore not that IPMs are insensitive to rare events, but that they avoid the near-deterministic count domination exhibited by the KL estimate. 6

Beyond EIG: Stable BOED with IPMs

Figure 1: Rare-event contamination experiment with ϵ = 0.01. Top row: Distributions of sample-based distance estimates across trials as the sample size n increases. Bottom row: How each estimator responds to the actual contamination count ncontam in the P batch for a fixed sample size n = 1000. The KL divergence is strongly sensitive to ncontam (showing an almost linear dependence), indicating extreme variance driven by rare samples. In contrast, the empirical W1 shows a stable two-sample noise floor that predictably decreases with n, while ED and MMD vary smoothly and robustly with ncontam .

This phenomenon is relevant to BOED because many designs are evaluated under observations for which most samples are weakly informative, while a small subset of rare observations can induce much larger posterior shifts. In such settings, estimators built on ratio-type quantities may become highly unstable unless additional structure or variancereduction machinery is introduced. This motivates our subsequent focus on utility formulations whose perturbation behavior can be controlled directly at the measure level.

4 Stability Analysis of IPMs Motivated by Helin et al. (2025), we ask the following question regarding utility stability: If one slightly perturbs the mathematical model (i.e., the likelihood) or the initial assumptions (i.e., the prior), how much does the design utility change? As shown in Section 3, the severe instability of the KL divergence in rare-event regimes motivates the paradigm shift toward IPMs. Our theoretical framework extends the recent work of Helin et al. (2025), who established the use of the Wasserstein distance as a valid, geometrically meaningful OED criterion, while we adopt a broader IPM framework. A key advantage of these metrics is that they avoid density ratios entirely, thereby circumventing the tail sensitivity that limits standard EIG. Within this unified IPM framework, our theory 7

Wu, Liang and Yang

highlights a fundamental theoretical trade-off between geometric richness and mathematical regularity. For MMD, we primarily consider the most commonly-used setting, translation-invariant radial kernels (Tolstikhin et al., 2017), such as Gaussian, Mixture of Gaussians, Inverse Multiquadrics, and Matérn kernels. For these kernels, the test functions are strictly bounded. By the reproducing property, we have sup ∥f ∥Hk ≤1

|f (x)| =

|⟨f, k(x, ·)⟩| ≤

sup ∥f ∥Hk ≤1

∥f ∥Hk

k(x, x) = O(1),

where the final equality follows directly from transition invariance. Consequently, MMD acts as a bounded IPM in the following analysis. In contrast, ED is generally unbounded. For the standard Euclidean energy distance, the growth rate of the test functions is |f (x)| ≤

sup ∥f ∥Hk

p ρ(x, x0 ) = O(∥x∥1/2 ).

≤1

The 1-Wasserstein distance exhibits even faster growth, with an envelope of O(∥x∥). To guarantee utility stability when using unbounded IPMs like ED and Wasserstein, stricter mathematical conditions are required: for instance, the prior exhibits sub-Gaussian tails, and the forward model is globally Lipschitz continuous. 4.1 Mathematical Setup and Notation We begin by formalizing our mathematical setting. Let X and Y denote the parameter and observation spaces, respectively. Because we consider a fixed experimental design throughout this section, we streamline our notation by leaving the design condition implicit within both the forward model and the resulting likelihood. For a given parameter x ∈ X , we define the forward model G : X → Y as a mapping to the observation space. We denote the class of all such potential forward mappings as G, where G ∈ G. The model generates data y ∈ Y according to y = G(x) + η, where η represents observational noise, thereby inducing the likelihood p(y|x). With a slight abuse of notation, for aRgiven prior distribution µ(dx), the evidence (i.e., marginal likelihood) is given by p(y) = X p(y|x)µ(dx), yielding the posterior distribution µy (dx) = p(y|x) p(y) µ(dx). Building on the IPM framework from Section 3, we define the expected utility of an experimental design based on the discrepancy between the prior and posterior. To guarantee that the IPM γF behaves as a true distance metric—specifically, satisfying symmetry such that γF (P, Q) = γF (Q, P ), the measurable function class F must be symmetric (i.e., f ∈ F ⇒ −f ∈ F ). This property holds natively for many standard function classes, including those defining the Wasserstein distance, MMD, and energy distance. The expected IPM utility for a given experimental design is therefore formulated as: Z

UF = Ep(y) [γF (µ, µ )] = Y

γF (µ, µy )p(y)dy.

Beyond EIG: Stable BOED with IPMs

4.2 Likelihood Stability In computationally intensive scientific applications, the true forward model G(x) is frequently replaced by a cheaper surrogate, such as a neural network. Because these surrogates are only approximations, they inevitably introduce localized errors into the likelihood function (Liu and Zhao, 2023). A robust experimental design framework must be resilient to these inaccuracies, meaning its utility degrades smoothly. As illustrated in our motivating example, the traditional EIG fails this criterion since a minor surrogate error in a low-probability region can cause the underlying KL divergence to explode logarithmically. To rigorously demonstrate how IPMs resolve this fragility and maintain stability under surrogate modeling errors, we shall first quantify the “topological severity” of the metric’s test functions. We capture this via a pointwise growth profile. Definition 4.1 (Growth Profile) For an IPM defined by a symmetric function class F , its pointwise growth profile is defined as ωF (x) := supf ∈F |f (x)|. The expected growth under R the prior µ is denoted as ω̄ = X ωF (x)µ(dx). Essentially, the growth profile characterizes how severely the metric penalizes large deviations. As discussed in the beginning of this section, for bounded IPMs such as MMD with standard radial kernels, the growth profile is strictly bounded: ωF (x) ≤ M < ∞. For unbounded metrics, we have ωF (x) = O(∥x∥1/2 ) for energy distance with standard Euclidean distance, and ωF (x) = O(∥x∥) for the 1-Wasserstein distance. Theorem 4.1 (Likelihood Stability Bound) Let two forward models yield likelihoods p(y|x) and p∗ (y|x), evidences p(y) and p∗ (y), and expected utilities UF and UF∗ , respectively. Then, the difference in their utility is bounded by the expected L1 error of the likelihoods, weighted by the growth profile: Z Z ∗ |UF − UF | ≤ ωF (x) + ω̄ p(y|x) − p∗ (y|x) µ(dx)dy. Y

Proof See Appendix A. The bound presented in Theorem 4.1 mathematically decouples the inherent topological properties of the metric from the approximation error of the surrogate model. Specifically, the utility perturbation is controlled linearly by the expected L1 error of the likelihoods, weighted by the growth profile ωF (x). This reveals a fundamental mechanism for stability where the penalty incurred by a surrogate’s error is strictly affected by how fast the metric’s chosen test functions are allowed to grow. Furthermore, our result clearly captures the distinct behaviors of bounded and unbounded metrics. For bounded IPMs, ωF (x) is globally capped, providing unconditional and uniform L1 control over utility errors. This guarantees that utility degrades gracefully alongside surrogate accuracy. Conversely, for unbounded IPMs, the growth of ωF (x) indicates that surrogate inaccuracies in the tails of the parameter space incur a higher penalty. This establishes a clear theoretical trade-off between the rich geometric transport of unbounded metrics and the universal stability guarantees of bounded ones. We provide a deeper investigation into the behavior of ωF (x) and its practical implications in Section 4.4. 9

Wu, Liang and Yang

Contrast with KL Divergence. To fully understand the robustness of IPMs, it is helpful to compare them with standard information-theoretic approaches. Consider the KL R p(x) divergence, namely, DKL (P ∥Q) := p(x) log q(x) dx, which implicitly utilizes the logratio as its test function. Because EIG relies on KL divergence, it scales with the logarithm of the likelihood ratio and reacts to multiplicative errors rather than additive L1 differences. Consequently, its stability cannot be neatly factored into a bounded growth profile. If a surrogate assigns near-zero likelihood where the true model assigns a non-negligible likelihood, the KL divergence blows up. This clearly demonstrates why EIG is highly unstable to even minor misspecifications in the tails of surrogate models. IPMs remove the ratio-based amplification mechanism of KL. For bounded IPMs this yields uniform control, while for unbounded IPMs rare observations may still matter, but their effect is governed by kernel geometry and moment conditions rather than logarithmic small-denominator instability. 4.3 Prior Stability While Section 4.2 establishes robustness against errors in the forward model, practical BOED faces another critical source of approximation error from the prior distribution. In most standard computational setups, the continuous prior µ is replaced by a discrete empirical measure µ̃ using Monte Carlo sampling. To bound the resulting utility error, we exploit the shift-invariance geometric property shared by standard IPMs. The following lemma demonstrates that we can always anchor our test functions to the approximate prior without changing the underlying distance. Lemma 4.1 (Prior-Anchored Supremum) Let F be a symmetric function class defining an IPM γF and µ̃ a perturbed probability measure. The prior-anchored subspace is defined as follows: Z Fprior =

g∈F

g(x)µ̃(dx) = 0 X

Then, the following statements hold: • The IPM between the two densities can be evaluated exactly via Z Z Z Z sup f (x)p(dx) − f (x)q(dx) = sup g(x)p(dx) − g(x)q(dx) . f ∈F

g∈Fprior

• Specifically, if we set X to be the parameter space X and p = µ̃, q = µ̃y , the IPM between the prior and posterior can be evaluated exactly over this restricted subspace Z y γF (µ̃, µ̃ ) = sup g(x)µ̃y (dx) ≜ D(y) g∈Fprior

For the convenience of the theoretical analysis, without loss of generality, we slightly x0 adjust the inducing function class for the 1-Wasserstein distance to FW := {f : f (x0 ) = 0, ∥f ∥L ≤ 1}. Because shifting a function by a constant does not change its Lipschitz conx0 stant, FW still perfectly defines the 1-Wasserstein distance. We denote the corresponding function classes for MMD and energy distance as Fk and Fkx0 , respectively. 10

Beyond EIG: Stable BOED with IPMs

To formalize how IPM utilities behave under prior perturbations, the observation noise model must be compatible with the chosen metric’s geometry. We assume the noise model is sufficiently smooth and bounded, ensuring the likelihood function interacts well with the function class F. Assumption 1 Let F be one of these three classes. For each y ∈ Y, define ey (·) :=

( p(y | ·),

F = Fk ,

p(y | ·) − p(y | x0 ),

x0 F = FW or F = FEx0 .

Assume there exist integrable functions CE (y) > 0 and CL (y) > 0 such that for almost every y ∈ Y and every f ∈ F, 1 ey ∈ F, CE (y)

1 f (·) p(y | ·) ∈ F. CL (y)

Assumption 1 imposes a pointwise analytic compatibility condition between the likelihood slice x 7→ p(y|x) and the function space defining the metric. It asserts that the evidence function belongs to the metric’s dual space, and that multiplication by the likelihood is a bounded operator on that space. Proposition 4.1 Consider the 1D parameter space X = [0, R] and the additive noise model Y = X + η, where η ∼ N (0, σ 2 ). The likelihood is given by p(y | x) = ϕσ (y − x), where ϕσ (z) = (2πσ 2 )−1/2 exp(−z 2 /2σ 2 ). We use H s (X ) to denote the standard fractional Sobolev space (Di Nezza et al., 2012) and As to be the universal constant of a Sobolev multiplier algebra when s > d/2 (Strichartz, 1967). Using x0 = 0 as the anchor point, Assumption 1 is satisfied for all three metrics with explicitly integrable local constants: • For Wasserstein-1 (W1 ): CE (y) = sup |ϕ′σ (y − x)|,

CL (y) = sup ϕσ (y − x) + R sup |ϕ′σ (y − x)|.

x∈[0,R]

• For Matérn-MMD (where Hk is norm-equivalent to H s ([0, R]) with s > 1/2, i.e., c1 ∥ · ∥H s ≤ ∥ · ∥Hk ≤ c2 ∥ · ∥H s ): CE (y) = c2 ∥ϕσ (y − ·)∥H s ([0,R]) ,

CL (y) =

c2 As ∥ϕσ (y − ·)∥H s ([0,R]) . c1

• For Energy Distance: CE (y) = ∥ϕ′σ (y − ·)∥L2 ([0,R]) ,

CL (y) = sup ϕσ (y − x) +

√

R∥ϕ′σ (y − ·)∥L2 ([0,R]) .

x∈[0,R]

Furthermore, for all three metrics, CE (y) and CL (y) exhibit Gaussian exponential tail decay as |y| → ∞, guaranteeing that they are integrable over y. 11

Wu, Liang and Yang

Proof See Appendix B. We emphasize that the example in Proposition 4.1 can be extended to the high-dimensional case. Consider a compact d-dimensional parameter space X ⊂ Rd , and the additive noise model Y = X + η, where η ∼ N (0, σ 2 Id ). The likelihood is given by p(y | x) = ϕσ (y − x), where ϕσ (z) = (2πσ 2 )−d/2 exp(−z 2 /2σ 2 ). Let R = supx∈X ∥x∥2 bound the domain radius. Under the same conditions as in Proposition 4.1, we have the following choices. For Wasserstein-1, we have CE (y) = supx∈X ∥∇x ϕσ (y − x)∥2 and CL (y) = supx∈X ϕσ (y − x) + R supx∈X ∥∇x ϕσ (y − x)∥2 . For Matérn-MMD, we have CE (y) = c2 ∥ϕσ (y − ·)∥H s (X ) and CL (y) = cc12 As ∥ϕσ (y − ·)∥H s (X ) . However, the corresponding derivation for the Euclidean Energy Distance is mathematically more delicate, requiring bounds on the posterior multiplier constant in a homogeneous fractional Sobolev space (Brasco et al., 2021) by verifying specialized Kato-Ponce inequalities (Grafakos and Oh, 2014). While conceptually analogous, we omit this derivation and focus on the Matérn-MMD case, which adequately captures the behavior of RKHS-based metrics in arbitrary dimensions via standard Sobolev multiplier algebras. We now state our main stability result for prior approximations. Theorem 4.2 (Prior Stability Bound) Let µ and µ̃ be the true and approximate prior distributions, yielding expected utilities UF and ŨF , respectively. Under Assumption 1, the difference in utility is linearly bounded by the IPM distance between the priors: |UF − ŨF | ≤ K̄γF (µ, µ̃), R where the stability constant is K̄ = 1 + Y CL (y) + 2CE (y)D(y) dy. Proof See Appendix B. Theorem 4.2 links the quality of the prior approximation directly to the resulting utility error. The utility perturbation is bounded linearly by the IPM discrepancy γF (µ, µ̃) between the approximate and true priors. The stability constant K̄ captures the sensitivity of this perturbation and depends primarily on the posterior discrepancy D(y). This reveals that the penalty incurred by approximating the prior depends strictly on how severely the observation y shifts the posterior distributions across the chosen test functions. In addition, this formulation highlights a key difference between bounded and unbounded metrics. For bounded IPMs , the test functions are capped, meaning the posterior discrepancy D(y) is universally bounded. As a result, the stability constant K̄ remains finite without additional restrictive assumptions on the true prior’s tail decay or the forward model’s growth rate. In contrast, for unbounded IPMs, D(y) can grow arbitrarily large if the forward model causes severe spatial warping. Ensuring a finite K̄ in these cases typically requires stronger conditions, such as assuming a sub-Gaussian prior and a globally Lipschitz forward model to restrict this growth. We provide a deeper investigation into the mechanics of D(y) and its practical implications in Section 4.4. 12

Beyond EIG: Stable BOED with IPMs

Contrast with KL Divergence: Standard EIG, which relies on KL divergence, is fundamentally incompatible with the discrete prior approximations commonly used in modern computational pipelines. Because KL divergence requires the distribution to be absolutely continuous with respect to one another, the divergence between a continuous prior and a discrete empirical measure is technically undefined. In contrast, because IPMs operate directly on probability measures via spatial transport or kernel embeddings rather than density ratios, they provide a much more natural theoretical justification for the empirical particle approximations used in practice. 4.4 Deeper Understanding of ωF (x) and D(y) Having established the main stability bounds, we now examine the two geometric quantities that govern them: the test-function growth profile ωF (x) and the posterior discrepancy D(y). In this section, we compare their behavior across a hierarchy of IPMs: bounded metrics (MMD, for which ωF (x) = O(1)), sublinear metrics (energy distance, for which ωF (x) = O(∥x∥1/2 )), and linear metrics (Wasserstein-1, for which ωF (x) = O(∥x∥)). We then show that controlling these quantities naturally leads to the assumptions used in related work, such as globally Lipschitz forward models and sub-Gaussian priors (Helin et al., 2025). A mathematically rigorous discussion is deferred to Appendix C. Likelihood Stability In the likelihood-stability bound, sensitivity to perturbations in the data misfit is governed by the test-function growth profile ωF (x). Assume a Gaussian observation model y = G(x) + η and a relative modeling perturbation satisfying ∥δ(x)∥ ≤ ϵ(1 + ∥G(x)∥). A first-order Gaussian shift estimate shows that R the likelihood perturbation enters the utility bound through the weighted prior integral X ωF (x) 1 + ∥G(x)∥ µ(dx). If G is globally Lipschitz, then ∥G(x)∥ ≲ 1 + ∥x∥, so the required prior integrability is determined by the growth of ωF (x). For bounded-kernel MMD, this yields a sufficient condition of a finite first moment. For the anchored energy-distance class, the corresponding sufficient condition is a finite 3/2-moment. For Wasserstein-1, the same argument gives a finite second moment. This makes the basic hierarchy precise: as the IPM becomes less bounded in space, the likelihood-stability bound requires stronger prior integrability. Prior Stability The prior-stability term behaves differently. The key quantity is D(y) = γF (µy , µ̃y ), which depends on the posterior reweighting induced by the likelihood. For bounded-kernel MMD, the test functions are uniformly bounded in the supremum norm, so D(y) is uniformly bounded and no tail-growth analysis is required. By contrast, unbounded IPMs involve transporting mass, and the local sensitivity is determined by the gradient of the effective integrand hy (x) = f (x) exp(−Φ(x; y)). Under a differentiable globally Lipschitz forward model, the gradient of the potential grows linearly, that is, as O(∥x∥). Consequently, the product rule introduces an additional factor of ∥∇x Φ(x; y)∥ into the effective test class. For Wasserstein-1, the envelope grows as O(∥x∥), so the resulting gradient grows quadratically, leading to a term of order O(∥x∥2 ). For the anchored energy-distance class, the envelope grows as O(∥x∥1/2 ), so the effective gradient grows superlinearly, namely as O(∥x∥1.5 ), depending on the precise formulation of the class. This leads to the same qualitative ordering as above: bounded-kernel MMD is the least sensitive to prior tails, the anchored energy-distance class is intermediate, and Wasserstein-1 is the most tail-sensitive. 13

Wu, Liang and Yang

In particular, strong tail assumptions, such as sub-Gaussian priors, arise naturally as convenient sufficient conditions in Wasserstein-based analyses.

5 Experiments In this section, we empirically evaluate our IPM-based BOED core framework and then its beyond-IPM plug-and-play extensions in high-dimensional settings. Building on benchmark tasks of Foster et al. (2019), we consider two settings of increasing complexity: (i) analytically tractable A/B testing problem, and (ii) preference learning model. In addition, to further illustrate the plug-and-play nature of our framework, we consider (iii) high-dimensional linear-Gaussian benchmark and (iv) high-dimensional sign-ambiguous benchmark. 5.1 Analytical Validation: A/B Testing We begin with an A/B testing scenario in which the expected utility can be computed exactly, allowing a rigorous validation of our formulation. Specifically, we consider two alternatives, denoted A and B. The design variable ξ = (nA , nB ) allocates a total budget of N = nA + nB samples between the two groups. The parameters xA and xB denote the mean responses of the respective groups. The quantity of interest is the treatment effect, δ := xB − xA , which measures the improvement of variant B over A and is invariant under global shifts. Our goal is therefore to identify a design ξ that maximally reduces uncertainty about δ. We assume a Gaussian prior x ∼ N (0, Σ0 ) and an observation model y|x, ξ ∼ N (Xξ x, I), where Xξ ∈ RN ×2 is the design matrix. The first nA rows of Xξ are (1, 0) and the remaining nB rows are (0, 1). By Bayes’ theorem, the posterior is given by x|y, ξ ∼ N (µ1 , Σ1 ), where −1 and µ = Σ X ⊤ y. ⊤ Σ1 = (Σ−1 1 1 ξ 0 + Xξ Xξ ) 2 , σ 2 ) with σ = 5 and For this experiment, we set the prior covariance to Σ0 = diag(σA A B 2 + σ 2 ). The σB = 1.8. Consequently, the prior over the treatment effect is δ ∼ N (0, σA B posterior parameters then simplify to:

Σ1 =

1 2 nA +1/σA

P

0 1 2 nB +1/σB

µ1 = 

i∈A yi 2 nP A +1/σA y i i∈B 2 nB +1/σB

 .

Denoting the marginal posterior distributions by xA | y, ξ ∼ N (mA , vA ) and xB | y, ξ ∼ N (mB , vB ), we obtain δ | y, ξ ∼ N (mB − mA , vA + vB ). The W1 distance between the prior and posterior distributions of δ is 2 2 W1 N (0, σA + σB ), N (mB − mA , vA + vB )

q √ 2 2 σ A + σ B − vA + v B Z , = E mA − m B +

where Z ∼ N (0, 1). Since mA − mB is a zero-mean Gaussian random variable, and the law of total variance gives Var(δ) = Var(E[δ | y]) + E[Var(δ | y)], we obtain 2 2 Var(mB − mA ) = σA + σB − (vA + vB ).

Beyond EIG: Stable BOED with IPMs

Therefore, the quantity inside the absolute value is itself a centered Gaussian random variable, which yields the closed-form expected utility r r q √ 2 2 + σ2 ) − 2 σ2 + σ2 v + v . U (ξ) = 2(σA A B A B B π 2 and σ 2 are constants, maximizing U (ξ) is equivalent to miniSince the prior variances σA B mizing the posterior variance 1 1 + vA + v B = 2 2 . nA + 1/σA nB + 1/σB

For N = 10, the continuous relaxation yields an optimal allocation of n∗A ≈ 5.13. The optimal discrete design is therefore nA = 5 and nB = 5. For comparison, we also evaluate the classical EIG based on KL divergence. By a similar derivation, the KL utility reduces to UKL (ξ) =

2 σ 2 + σB 1 log A . 2 vA + v B

Thus, maximizing the KL utility is also equivalent to minimizing vA + vB . Both criteria therefore identify the same optimal design, namely (nA , nB ) = (5, 5). This confirms that the Wasserstein utility preserves the correct global optimum in this standard setting while benefiting from the improved stability properties established in Section 4. We use this example mainly as a sanity check: it shows that replacing KL by an IPM does not distort the global optimum in a standard setting where the utility can be analyzed exactly. Differences in the shape of the resulting utility landscapes are explored empirically in the next subsection. 5.2 Computational Performance: Preference Learning To evaluate our framework in a setting without a closed-form posterior and with nonlinear observations, we next consider a preference learning model. In this more challenging setting, we compare IPM-based utilities with the classical KL-based criterion from the perspectives of both computational efficiency and optimization reliability. Model Formulation. We introduce a latent parameter x ∈ R representing a user’s intrinsic preference threshold, together with a continuous scalar design variable ξ ∈ [ξmin , ξmax ] representing a characteristic of the proposed item. We place a Gaussian prior x ∼ N (µx , σx2 ), with default values µx = −20 and σx = 20. Conditional on the parameter-design pair (x, ξ), we model the latent utility η as η | x, ξ ∼ N (ξ − x, ση2 (1 + |ξ|)2 ), where ση = 1. The noise variance therefore grows quadratically with the magnitude of the design, which naturally penalizes extreme design choices. The final observed response is obtained by passing the latent utility through a sigmoid link and then truncating the result for numerical stability: y = max(ϵ, min(1 − ϵ, sigmoid(η))), with ϵ = 0.01. This yields a bounded continuous observation y ∈ [ϵ, 1 − ϵ], which can be interpreted as a noisy preference probability. 15

Wu, Liang and Yang

Experimental Setup and Computational Speedup. We discretize the design space on a uniform grid from −80 to 80 with step size 2, yielding 81 candidate designs.

Figure 2: Computational runtime comparison between KL-based EIG and IPM-based utilities. The plots report the per-design evaluation time for both the A/B testing model (left) and the preference learning model (right), highlighting the substantial speedup achieved by density-free IPMs. As shown in Figure 2, the empirical results strongly support the computational advantage of the density-free IPM approach. Under the current implementation for this experiment, the KL-based baseline is substantially slower because it requires an additional density-estimation step before evaluating the log-ratio utility, together with bandwidth tuning that is both computationally expensive and numerically sensitive. By contrast, the IPM utilities considered here are evaluated directly from samples through pairwise distances or kernels, which avoids this intermediate density-estimation stage. As a result, in this experiment the W1-based utility requires only 0.3 seconds per design, whereas the KL baseline requires about 14 seconds per design. Figure 2 also includes the runtime for the A/B testing experiment in Section 5.1, showing that IPMs retain a consistent efficiency advantage in both analytically tractable and numerically approximated settings. Utility Landscape and Optimization Reliability. Beyond computational speed, a central practical observation is that IPMs also improve the geometry of the optimization landscape. Figure 3 displays the expected utility curves across the candidate designs and shows how the choice of utility metric changes the structure of the design objective.

Figure 3: Expected utility landscapes across candidate designs for preference learning. To formalize this optimization reliability, we define the high-utility design region Rt = {d : U (d) ≥ t · max U }. 16

Beyond EIG: Stable BOED with IPMs

Conceptually, this is analogous to a Highest Posterior Density (HPD) credible region in Bayesian inference (Lee, 1989), but applied to the experimental design space: it characterizes the basin of near-optimal and highly informative designs. The size of this region is measured by the number of discrete design points contained in Rt , while the concentration score SR quantifies how localized the utility mass is over the design space.

Metric KL Divergence Energy Distance MMD2 Wasserstein-1

SR 0.1213 0.3186 0.2158 0.2877

Region R0.80 ⊆ [−2, 0] ⊆ [−8, 6] ⊆ [−6, 0] ⊆ [−18, 2]

Details (2 pts ≥ 80% of max) (7 pts ≥ 80% of max) (4 pts ≥ 80% of max) (11 pts ≥ 80% of max)

Figure 4: Characterization of the high-utility design region R0.80 = {d : U (d) ≥ 0.8·max U }. At the threshold 0.8 max U , IPM-based utilities induce a substantially broader and more robust basin of near-optimal designs than the classical KL divergence.

Figure 5: Sensitivity analysis of the high-utility design region across threshold levels t. For a broad range of practically relevant thresholds, IPM-based utilities retain a larger number of acceptable design points (right) together with a smoother concentration profile (left), indicating greater optimization stability than KL divergence. Figures 4 and 5 show that KL- and IPM-based utilities differ not only in value, but also in landscape geometry. At the representative threshold t = 0.8, IPM-based utilities retain 17

Wu, Liang and Yang

a substantially larger superlevel set {d : U (d) ≥ t Umax }, indicating that highly informative designs occupy a much broader portion of the design space. By contrast, the KL-based criterion concentrates acceptable designs in a much narrower region around the maximizer. This distinction is important for downstream optimization. In practice, the selected design is affected by grid discretization, Monte Carlo error, approximate posterior computation, and imperfect search. A broader near-optimal basin therefore makes the design recommendation less sensitive to such perturbations, because many nearby designs still achieve utility close to the optimum. In this sense, IPM-based objectives transform the problem from finding a single sharp but unstable peak to identifying a much more stable set of highly competitive designs. We now analyze whether this geometric property translates into a practical advantage when the design optimization is approximate. To test whether a wider near-optimal region leads to more robust design selection under approximate optimization, we evaluate each metric under coarse-grid selection. We restrict the candidate set to every k-th point of Ξ for k ∈ {2, 4, 6, 8}. Denote the reference utility as Uref and the optimal design as ξ ∗ . For a given k there are k distinct coarse grids, obtained by starting at offset s ∈ {0, 1, . . . , k − 1}. For each offset we select the design ξˆ with the highest reference utility on that coarse set, and measure the normalized regret ˆ = r̄(ξ)

ˆ Uref (ξ ∗ ) − Uref (ξ) . ∗ Uref (ξ ) − minξ Uref (ξ)

We report the mean and maximum of r̄ over all k offsets. Because selection is based on the reference utility itself, this protocol isolates the effect of grid coarsening alone, without additional Monte Carlo noise. Averaging over offsets removes any dependence on a particular grid alignment and yields a more general measure of robustness. ˆ averaged over all k offsets. Lower Table 1: Mean and maximum normalised regret r̄(ξ) regret indicates better recovery under coarse-grid restriction. |R0.80 | stands for the number of points within the high-utility design region as shown in Figure 4. Mean norm. regret r̄ (max in parentheses) Metric KL divergence Energy distance MMD2 Wasserstein-1

|R0.80 |

k=2

k=4

k=6

k=8

2 7 4 11

0.018 (0.035) 0.068 (0.135) 0.080 (0.160) 0.001 (0.002)

0.121 (0.231) 0.093 (0.151) 0.118 (0.161) 0.023 (0.056)

0.180 (0.337) 0.117 (0.168) 0.150 (0.225) 0.038 (0.081)

0.234 (0.419) 0.148 (0.299) 0.195 (0.383) 0.054 (0.134)

Table 1 reports the resulting mean and maximum normalized regret across offsets. The results provide systematic, alignment-robust evidence in this preference-learning example that broader high-utility regions are associated with smaller regret under coarse-grid selection. Across all reported coarsening levels, Wasserstein-1 achieves the lowest mean and worst-case regret, while KL exhibits the strongest degradation as the grid becomes coarser. Energy Distance and MMD2 lie in between. Moreover, for each k, both the mean and maximum regret are ordered consistently with the size of the 0.80-superlevel set. This indi18

Beyond EIG: Stable BOED with IPMs

cates that the width of the near-optimal design region captures an operationally meaningful aspect of robustness to design discretization. The coarse-grid study above provides concrete one-step evidence that broader highutility regions can reduce the sensitivity of design selection to discretization. This observation may also be relevant for sequential OED, where design choices are repeatedly made under approximate posterior updates and imperfect numerical optimization. A broader near-optimal region may then help reduce the propagation of local optimization errors across stages. We emphasize, however, that this sequential interpretation remains suggestive in the present work, and a direct empirical study of multi-stage robustness is left for future work. 5.3 High-Dimensional Plug-and-Play Extensions Beyond IPMs: Scalability and Robustness The experiments in previous sections show that IPM-based utilities are numerically stable, computationally efficient, and more reliable than classical KL-based criteria. We now show that the same sample-based BOED pipeline can also accommodate geometry-aware discrepancies beyond the IPM class. First, we study whether the sample-based framework can be combined directly with scalable neural estimators in a high-dimensional problem where exact ground truth is available. Second, we test whether this approach remains reliable when the posterior geometry becomes strongly non-Gaussian. To this end, we consider two highdimensional benchmarks with complementary structure: a 64-dimensional linear-Gaussian model with a tractable unimodal posterior, and a 32-dimensional sign-ambiguous model with a bimodal posterior. Together, these experiments show that the proposed framework is not only scalable, but also robust to posterior features that are known to cause difficulty for KL-based approximations: density estimation fails under the curse of dimensionality (Li et al., 2024), and standard nested Monte Carlo (NMC) estimators can suffer from exponentially growing variance (Foster et al., 2019). Plug-and-Play Neural Estimators. To evaluate the design utility empirically (simulating a scenario where the closed-form equation is unknown), we utilize the OT-ICNN method in Makkuva et al. (2020). This approach parameterizes the Kantorovich dual potentials using Input Convex Neural Networks (ICNNs), which enforce the convexity structure required by optimal transport. Importantly, incorporating this neural estimator requires no modification of the BOED framework itself. For each design ξ and sampled observation y, we draw fresh minibatches at every optimization step. We then train two ICNN potentials under the semidual objective, alternating multiple updates of the transport network for each update of the critic, together with a convexity regularization penalty on the constrained weights. In these two experiments, we use three hidden layers of width 256 for benchmark I and 128 for benchmark II, a batch size of 128, the Adam optimizer with learning rate 10−4 , and a convexity penalty coefficient of 0.1. Final utility estimates are obtained by averaging the learned OT objective over several fresh evaluation minibatches. For comparison, we estimate the KL-based EIG using nested Monte Carlo with 64 outer samples and 512 inner prior samples per outer sample for benchmark I and Gaussian variational methods with 50000 samples in total for benchmark II. 19

Wu, Liang and Yang

Benchmark I: A High-Dimensional Linear-Gaussian Model. We first consider a linear-Gaussian design problem in dimension p = 64. The prior is x ∼ N (0, Ip ), and the design variable is a scalar gain ξ ∈ Ξ = {0.5, 2.0, 6.0}. Given a design ξ, the observation model is y | x, ξ ∼ N (ξx, σ 2 Ip ), with σ = 0.1. Because this model is conjugate (Raiffa and Schlaifer, 2000), the posterior is available in closed form: p(x | y, ξ) = N (µξ (y), Σξ ), where

ξ σ2 y, Σ = Ip . ξ ξ2 + σ2 ξ2 + σ2 This benchmark serves as a clean test of scalability. Since the posterior is unimodal and analytically tractable, it isolates the effect of dimension without introducing additional approximation error from complicated posterior geometry. It also allows us to compute the exact expected 2-Wasserstein utility in closed form: ! 2 σ U (ξ) = Ey|ξ W2 N (0, Ip ), N (µξ (y), Σξ ) = 2p 1 − p , ξ2 + σ2 µξ (y) =

which provides a ground-truth reference for evaluating the neural estimator. Results on Benchmark I. We compare the OT-ICNN Wasserstein estimator with a standard prior-sampling nested Monte Carlo estimator for the classical EIG. Table 2 reports the results for the candidate designs. Table 2: Comparison with exact utilities in the 64-dimensional linear-Gaussian BOED problem. ICNN is evaluated against the exact Wasserstein utility, while nested Monte Carlo is evaluated against the exact expected information gain (EIG). The relative error is defined by |Û − U ⋆ |/|U ⋆ |.

Design

Exact

W22 utility ICNN

ξ = 0.5 ξ=2 ξ=6

1.03 × 102 1.22 × 102 1.26 × 102

1.01 × 102 1.21 × 102 1.25 × 102

Rel. err.

Exact

EIG utility Nested MC

Rel. err.

2.18 × 10−2 2.70 × 10−3 3.10 × 10−3

1.04 × 102 1.92 × 102 2.62 × 102

9.49 × 102 1.54 × 104 1.34 × 105

8.10 × 100 7.91 × 101 5.12 × 102

The results show a clear contrast between the two approaches. The OT-ICNN estimator scales effectively to p = 64, closely matches the exact expected W22 utility, and recovers the correct ranking of the candidate designs, with relative errors of 0.0218, 0.0027, and 0.0031. By contrast, the naive nested Monte Carlo estimator for the KL-based EIG becomes highly inaccurate in this high-dimensional, low-noise regime, with relative errors ranging from 8.10 to 5.12 × 102 . Thus, even in a mild unimodal problem, standard EIG estimation can become unreliable in high dimension unless substantial problem-specific variance-reduction techniques are introduced. In contrast, the sample-based BOED formulation allows scalable neural estimators to be integrated directly into the framework, making high-dimensional utility estimation feasible without altering the underlying method. 20

Beyond EIG: Stable BOED with IPMs

Benchmark II: A High-Dimensional Sign-Ambiguous Bimodal Model. We now turn to a second high-dimensional benchmark that stresses a different issue. The challenge here is not conjugate high-dimensional scaling, but severe posterior non-Gaussianity. Specifically, we construct a 32-dimensional sign-ambiguous model in which the posterior is bimodal and therefore poorly matched by Gaussian variational approximations. Let x ∈ R32 . Each coordinate has an independent bimodal prior 1 1 p(xi ) = N (−µi , τ 2 ) + N (µi , τ 2 ), 2 2 with τ = 0.35. A design ξ ⊆ [32], with |ξ| = 8, observes the selected coordinates through yi = x2i + ϵi ,

ϵi ∼ N (0, σ 2 ),

σ = 0.20.

Because the map xi 7→ x2i removes the sign of xi , the likelihood is insensitive to the sign of xi . As a result, the posterior is sign-ambiguous, and the degree of non-Gaussianity increases with µi . The prior and likelihood both factorize across coordinates, so the exact design-level utilities decompose as sums of one-dimensional contributions. Table 3 reports the resulting per-coordinate reference values for four coordinate types. Table 3: Per-coordinate 1-D reference utilities for each µi type. Design-level utilities are sums over the eight observed coordinates. Type

Coord. indices

µi

EIG1d

2 W2,1d

Strong Medium Weak Null

1–8 9–16 17–24 25–32

4.0 2.0 1.0 0.1

2.637 1.941 1.251 0.265

0.243 0.200 0.167 0.031

We consider five candidate designs, each observing 8 coordinates, chosen to span a broad range of information profiles: A (8 strong), B (8 medium), C (8 weak), D (4 strong + 4 null), and E (8 null). Table 4: Utilities and design rankings on the sign-ambiguous benchmark (d = 32, |ξ| = 8, seed 42). Rankings in parentheses (1 = best). W22

EIG Design Composition A B C D E Ranking

8 strong 8 medium 8 weak 4 strong + 4 null 8 null

Exact Gaussian VB

Exact OT-ICNN

21.09 (1) 15.53 (2) 10.00 (4) 11.61 (3) 2.12 (5)

−13.91 (5) −8.27 (4) −2.53 (2) −6.34 (3) +1.23 (1)

1.943 (1) 1.596 (2) 1.337 (3) 1.097 (4) 0.251 (5)

1.608 (1) 1.585 (2) 1.249 (3) 0.793 (4) 0.161 (5)

ABDCE

ECDBA

ABCDE

Wu, Liang and Yang

Results on Benchmark II. To evaluate KL-based design under posterior misspecification, we use the Barber–Agakov lower bound (Barber and Agakov, 2004), EIG(ξ) ≥ H[p(x)] + Ep(x,y) [log qϕ (x | y)] , where qϕ is a Gaussian variational approximation. This lower bound is accurate only when qϕ approximates the true posterior well. In the present model, however, the posterior is bimodal, so any single Gaussian must place mass between the two modes rather than on them. A straightforward calculation can show that for large µi the resulting per-coordinate lower bound is negative. This explains the numerical results in Table 4 that the bound becomes most negative on the most informative coordinates, while remaining close to zero on the nearly uninformative null coordinates. This completely reverses the design ranking. By contrast, the OT-ICNN estimator within our Wasserstein BOED framework remains stable. Using 50 Monte Carlo samples, it produces the ranking A > B > C > D > E, which matches the true W22 ordering exactly. The relatively larger errors on designs A and D reflect the more difficult geometry of well-separated or mixed-type posteriors, but the separation between adjacent designs is large enough for all designs to be ranked correctly. This benchmark also highlights a conceptual difference between EIG and W22 . Exact EIG ranks D above C (11.61 versus 10.00) because the four null coordinates in design D still contribute 1.06. Exact W22 , however, ranks C above D because the same null coordinates contribute only 0.12, which is less than 3% of design C’s total utility. Thus, W22 penalizes the allocation of design budget to coordinates that induce only negligible posterior change more strongly than EIG.

6 Conclusion In this work, we have revisited Bayesian optimal experimental design through the lens of utility design. Our starting point is a structural limitation of the classical expected information gain (EIG) criterion: as a KL-based utility, it depends on a log-density ratio and can therefore be sensitive to rare events, tail underestimation, and model approximation errors. To address this issue, we have introduced a geometry-aware BOED framework based on Integral Probability Metrics (IPMs), with particular emphasis on the 1-Wasserstein distance, Maximum Mean Discrepancy, and Energy Distance. Our theoretical analysis has established systematic stability guarantees and has shown that IPM-based utilities naturally mitigate the tail sensitivity inherent in density-ratiobased metrics. In particular, we have shown that IPMs admit stability bounds controlled directly by the error in the surrogate likelihood or the prior distribution. Empirically, we have demonstrated that the proposed framework yields more stable and reliable behavior than classical EIG, measured via broader near-optimal design regions, smoother utility landscapes, and improved robustness. By incorporating a plug-and-play neural optimal transport estimator, we have also extended the same sample-based BOED template to high-dimensional design spaces where commonly used EIG estimators become unreliable. 22

Beyond EIG: Stable BOED with IPMs

Future Work The modular and plug-and-play structure of our IPM-based framework opens several directions for future research. A natural next step is to extend the present static BOED framework to sequential OED and active learning settings (He and Geng, 2008; Dror and Steinberg, 2008; MacKay, 1992). Our coarse-grid robustness results suggest that utility-landscape geometry may play an important role when design choices must be made repeatedly under approximate optimization, but a direct empirical study of sequential robustness remains to be carried out. Moreover, because our framework relies on samplebased discrepancy evaluation rather than density tracking, it is naturally compatible with modern generative modeling techniques. Incorporating advanced sampling methods, such as diffusion models (Ho et al., 2020) and flow matching (Lipman et al., 2022), into the utility estimators may further improve scalability and fidelity for highly complex and multimodal target distributions. More broadly, the modular structure of the proposed framework may also be relevant beyond classical BOED. For example, related ideas may be useful for data selection in large-scale model training and for active preference learning in human-in-the-loop systems (Sorscher et al., 2022; Xie et al., 2023; Rafailov et al., 2023; Muldrew et al., 2024). We hope that the present work provides a useful step toward more robust and scalable design criteria in resource-constrained and high-dimensional inference problems.

Acknowledgments and Disclosure of Funding HY was partially supported by the US National Science Foundation under awards IIS2520978, GEO/RISE-5239902, the Office of Naval Research Award N00014-23-1-2007, DOE (ASCR) Award DE-SC0026052, and the DARPA D24AP00325-00. DW thanks Chenguang Duan for the discussion of sampling. Approved for public release; distribution is unlimited.

Appendix A. Proof of Likelihood Stability Proof [Proof of Theorem 4.1] By definition, the expected utility can be written as: Z Z Z UF = sup f (x)p(y|x)µ(dx) − p(y) f (x)µ(dx) dy Y f ∈F

Define the unperturbed and perturbed inner operators for a given y and f : Z Z A(f, y) = f (x)p(y|x)µ(dx) − p(y) f (x)µ(dx), X X Z Z ∗ A (f, y) = f (x)p∗ (y|x)µ(dx) − p∗ (y) f (x)µ(dx). X

Then, the difference in utility can be bounded by Z ∗ |UF − UF | ≤ sup |A(f, y)| − sup |A∗ (f, y)| dy. Y

Using the reverse triangle inequality for suprema, we have: Z ∗ |UF − UF | ≤ sup |A(f, y) − A∗ (f, y)|dy. Y f ∈F

Wu, Liang and Yang

Now, we bound the difference inside the supremum: Z Z ∗ |f (x)|µ(dx). |f (x)||p(y|x) − p∗ (y|x)|µ(dx) + |p(y) − p∗ (y)| |A(f, y) − A (f, y)| ≤ X

We R bound both terms uniformly using the growth profile ωF (x). Note that X |f (x)|µ(dx) ≤ X ωF (x)µ(dx) = ω̄: Z ωF (x)|p(y|x) − p∗ (y|x)|µ(dx) + ω̄|p(y) − p∗ (y)|. |A(f, y) − A∗ (f, y)| ≤ X

Taking the supremum over f simply yields this same upper bound. Integrating over y: Z Z Z |p(y) − p∗ (y)|dy. ωF (x)|p(y|x) − p∗ (y|x)|µ(dx)dy + ω̄ |UF − UF∗ | ≤ Y

Notice that the difference in evidences is bounded by the difference in likelihoods: Z Z Z Z Z (p(y|x) − p∗ (y|x))µ(dx) dy ≤ |p(y|x) − p∗ (y|x)|µ(dx)dy. |p(y) − p∗ (y)|dy = Y

Substituting this into the previous inequality yields the desired bound: Z Z ∗ |UF − UF | ≤ (ωF (x) + ω̄)|p(y|x) − p∗ (y|x)|µ(dx)dy. Y

Appendix B. Proof of Prior Stability R R Proof [Proof of Lemma 4.1] Let H(f ) = X f (x)p(dx) − X f (x)q(dx) . We wish to prove supf ∈F H(f ) = supg∈Fprior H(g). Because Fprior ⊂ F, it trivially holds that supf ∈F H(f ) ≥ supg∈Fprior H(g). To prove the reverse inequality, consider any arbitrary function f ∈ F. R Let c = X f (x)p(dx).R By the shift-invariance of F, the function g(x) = f (x) − c belongs to F. By construction, X g(x)p(dx) = 0, meaning that g ∈ Fprior . We evaluate the posterior discrepancy score H(g) for this shifted function: Z Z (f (x) − c)p(dx) − (f (x) − c)q(dx) H(g) = X

y y By the linearity of the integral, we can separate the constant R y c. SinceR µ yand µ̃ are both valid probability measures, they integrate to 1 (i.e., dµ = 1 and dµ̃ = 1). The constants c perfectly cancel. Hence, it holds that Z Z f (x)p(dx) − f (x)q y (dx) = H(f ) H(g) = X

Because every function f ∈ F has a corresponding shifted function g ∈ Fprior that yields the exact same absolute difference, the supremum over the full class cannot exceed the 24

Beyond EIG: Stable BOED with IPMs

supremum over the restricted class. Thus, the supremums are strictly equal. Proof [Proof of Proposition 4.1] Let X = [0, R]. We verify the pointwise existence of CE (y) and CL (y) for each metric, followed by their global y-integrability. 0 = {f : f (0) = 0, ∥f ′ ∥ Part 1: Wasserstein-1 (W1 ). The dual class is FW ∞ ≤ 1}. The centered evidence function is ey (x) = ϕσ (y − x) − ϕσ (y − 0). By definition, ey (0) = 0. The 0 is the Lipschitz constant of e . By the minimal constant CE (y) such that CE1(y) ey ∈ FW y d ey (x) = −ϕ′σ (y − x), giving: chain rule, dx

d ey (x) = sup |ϕ′σ (y − x)|. dx x∈[0,R] x∈[0,R]

CE (y) = sup

0 . We must bound the Lipschitz constant of For the posterior constant CL (y), let f ∈ FW h(x) = f (x)ϕσ (y − x). By the product rule:

|h′ (x)| ≤ |f ′ (x)|ϕσ (y − x) + |f (x)||ϕ′σ (y − x)|. 0 , |f ′ (x)| ≤ 1, we have ∥f ∥ ′ Since f ∈ FW ∞ ≤ R∥f ∥∞ ≤ R. Taking the supremum over x ∈ [0, R] yields:

∥h′ ∥∞ ≤ sup ϕσ (y − x) + R sup |ϕ′σ (y − x)| = CL (y). x∈[0,R]

x∈[0,R]

Since both CE (y) and CL (y) are defined by the suprema of a Gaussian and its derivative over a compact interval, they exist pointwise for all y. Part 2: Matérn-MMD. The dual class is defined strictly by the native norm: Fk = {f ∈ Hk : ∥f ∥Hk ≤ 1}.The uncentered evidence function is ey (x) = ϕσ (y − x). Because the Gaussian is smooth, its restriction to [0, R] belongs to H s ([0, R]) for any s ≥ 0. To guarantee CE1(y) ey ∈ Fk , we require its Hk norm to be bounded by 1. Using the upper equivalence bound: ∥ey ∥Hk ≤ c2 ∥ey ∥H s ([0,R]) . Thus, defining CE (y) = c2 ∥ϕσ (y − ·)∥H s ([0,R]) satisfies the evidence condition. For the posterior constant CL (y), let f ∈ Fk . By the lower equivalence bound, we know ∥f ∥H s ≤ 1 1 s c1 ∥f ∥Hk ≤ c1 . Because s > 1/2, H ([0, R]) is a Banach algebra with multiplier constant As > 0. We bound the Hk norm of the product f (·)ϕσ (y − ·) by chaining the upper equivalence, the Banach algebra property, and the lower equivalence: ∥f · ey ∥Hk ≤ c2 ∥f · ey ∥H s ≤ c2 As ∥f ∥H s ∥ey ∥H s ≤

c2 As ∥ey ∥H s . c1

We define CL (y) = c2cA1 s ∥ϕσ (y − ·)∥H s ([0,R]) , which is strictly finite. Part 3: Energy Distance. The native space corresponding to the Euclidean energy distance in 1D, anchored at 0, yields the dual class FE0 = {f : f (0) = 0, f absolutely continuous, ∥f ′ ∥L2 ≤ 1}. The centered evidence function is ey (x) = ϕσ (y − x) − ϕσ (y − 0). We have Z R CE (y) = 0

d ey (x) dx

!1/2

2 dx

Z R =

(−ϕ′σ (y − x))2 dx

1/2

= ∥ϕ′σ (y − ·)∥L2 ([0,R]) .

Wu, Liang and Yang

For CL (y), let f ∈ FE0 . The multiplier h(x) = f (x)ϕσ (y − x) must satisfy h(0) = 0 and have a bounded L2 derivative. By the product rule and Minkowski’s inequality: ∥h′ ∥L2 = ∥f ′ ϕσ + f ϕ′σ ∥L2 ≤ ∥f ′ ϕσ ∥L2 + ∥f ϕ′σ ∥L2 . We bound the terms separately. First, ∥f ′ ϕσ ∥L2 ≤ ∥f ′ ∥L2 ∥ϕσR(y − ·)∥∞ ≤ 1 · supx ϕσ (y −√x). √ x Second, because f (0) = 0, Cauchy-Schwarz implies |f (x)| = | 0 f ′ (t)dt| ≤ x∥f ′ ∥L2 ≤ R. √ Thus, ∥f ϕ′σ ∥L2 ≤ R∥ϕ′σ (y − ·)∥L2 . Combining these gives: √ ∥h′ ∥L2 ≤ sup ϕσ (y − x) + R∥ϕ′σ (y − ·)∥L2 ([0,R]) = CL (y). x∈[0,R]

Part 4: Global Integrability in y. For all three metrics, CE (y) and CL (y) are linear combinations of the supremum, L2 , or H s norms of ϕσ (y − ·) and its derivatives over x ∈ [0, R]. As |y| → ∞, the distance between y and the compact set [0, R] grows as |y| − R. Consequently, the magnitudes of ϕσ (y − 2x) and its spatial derivatives on x ∈ [0, R] are uniformly (|y|−R) bounded above by K exp − 2σ2 for some polynomial K(y). Because the Gaussian tail decays faster than any polynomial, CE (y) and CL (y) decay exponentially as |y| → ∞. Thus they are integrable. Proof [Proof of Theorem 4.2] By expanding the definitions of expected utility and applying the triangle inequality to the integrals over Y, we can bound the difference by the supremum of the pointwise integrand errors: Z Z y y |p(y) − p̃(y)|γF (µ̃, µ̃y )dy |UF − ŨF | ≤ p(y) |γF (µ, µ ) − γF (µ̃, µ̃ )| dy + |Y {z } |Y {z } Term A

Term B

Bounding Term A: By the reverse triangle inequality for metrics, |γF (µ, µy ) − γF (µ̃, µ̃y )| ≤ γF (µ, µ̃) + γF (µy , µ̃y ). Multiplying by p(y) and integrating: Z Term A ≤ γF (µ, µ̃) +

p(y)γF (µy , µ̃y )dy.

Now we bound the expected difference between posteriors with IPM. Using the definition of the posterior and Lemma 4.1: Z Z p(y) y y p(y)γF (µ , µ̃ ) = sup g(x)p(y|x)µ(dx) − g(x)p(y|x)µ̃(dx) p̃(y) g∈Fprior Z f (x)p(y|x)(µ(dx) − µ̃(dx)) ≤ sup f ∈F

Z + sup

g(x)p(y|x)µ̃(dx) 1 −

g∈Fprior

p(y) p̃(y)

Beyond EIG: Stable BOED with IPMs

R R 1 g(x)p(y|x)µ̃(dx) = gdµ̃y . The supremum of this For the second part, note that p̃(y) within Fprior is exactly D(y). For the first part, applying Assumption 1, f (·)p(y|·) = CL (y)h(·) for some h ∈ F . Thus: Z Z sup f (x)p(y|x)(µ(dx) − µ̃(dx)) ≤ CL (y) sup h(x)(µ(dx) − µ̃(dx)) = CL (y)γF (µ, µ̃) f ∈F

h∈F

Integrating these bounds for Term A over y: Z Z D(y)|p̃(y) − p(y)|dy CL (y)γF (µ, µ̃)dy + Term A ≤ γF (µ, µ̃) + Y Y Z Z ≤ (1 + CL (y)dy)γF (µ, µ̃) + D(y)|p(y) − p̃(y)|dy Y

Bounding Term B: According to Lemma 4.1, we have: Z D(y)|p(y) − p̃(y)|dy Term B ≤ Y

Adding Term A and Term B gives Z |UF − ŨF | ≤ (1 + C̄L )γF (µ, µ̃) + 2

D(y)|p(y) − p̃(y)|dy Y

Again, we use Assumption 1 Z

|p(y) − p̃(y)| =

p(y|x)µ(dx) − p(y|x)µ̃(dx) Z Z ey (x)µ(dx) − ey (x)µ̃(dx) = Z Z ≤ CE (y) sup h(x)µ(dx) − h(x)µ̃(dx) h∈F

= CE (y)γF (µ, µ̃). Now, integrate this over all observations y: Z Z D(y)|p(y) − p̃(y)|dy ≤ CE (y)D(y)dy γF (µ, µ̃). Y

This gives the final bound.

Appendix C. Rigorous Discussion for Section 4.4 Likelihood Stability Consider the observation is generated by a forward model with Gaussian noise: y = G(x) + η, where η ∼ N (0, Σ). The true data misfit potential is: Φ(x; y) =

⊤ 1 y − G(x) Σ−1 y − G(x) 2 27

Wu, Liang and Yang

The true likelihood is p(y|x) = Z1 exp(−Φ(x; y)). Let the perturbed model be G∗ (x) = G(x) + δ(x). Assuming the modeling error scales with the model magnitude, the perturbation is bounded as ∥δ(x)∥ ≤ ϵ(1 + ∥G(x)∥). The difference in potentials is bounded by: |∆Φ| = |Φ∗ (x; y) − Φ(x; y)| ⊤ 1 = − y − G(x) Σ−1 δ(x) + δ(x)⊤ Σ−1 δ(x) 2 ≤ϵ(1 + ∥G(x)∥)∥Σ−1 ∥∥y − G(x)∥ + O(ϵ2 ). By the Mean Value Theorem, for a sufficiently small model perturbation ϵ → 0, the secondorder term vanishes, and we obtain the first-order bound on the likelihoods: h i |p∗ (y|x) − p(y|x)| ≤ p(y|x) · exp(|∆Φ|) · |∆Φ| ≲ p(y|x) · ϵ(1 + ∥G(x)∥)∥Σ−1 ∥∥y − G(x)∥ . Substituting this into the general expected utility bound yields:

|UF − UF∗ | ≲ ϵ∥Σ−1 ∥

Z ωF (x)(1 + ∥G(x)∥) X

∥y − G(x)∥p(y|x)dy µ(dx).

Notice that the inner integral over Y is simply the expected magnitude of the Gaussian noise E[∥η∥], which evaluates to a constant K. The bound simplifies to: Z |UF − UF∗ | ≲ ϵK∥Σ−1 ∥ ωF (x)(1 + ∥G(x)∥)µ(dx). X

If we assume G(x) is globally Lipschitz (L1 ), then ∥G(x)∥ ≤ L1 ∥x∥ + ∥G(0)∥. The polynomial order of the integrand dictates the prior moment requirements: • MMD (ωF (x) = O(1)): Integrand is O(∥x∥) =⇒ A finite 1st moment is sufficient. • Energy Distance (ωF (x) = O(∥x∥1/2 )): Integrand is O(∥x∥1.5 ) =⇒ A finite 1.5-th moment is sufficient. • Wasserstein-1 (ωF (x) = O(∥x∥)): Integrand is O(∥x∥2 ) =⇒ A finite 2nd moment is sufficient. Prior Stability

The posterior discrepancy is defined by the dual representation: Z Z y y y D(y) = γF (µ , µ̃ ) = sup f (x)µ (dx) − f (x)µ̃y (dx) . f ∈F

For bounded metrics (MMD), supf ∈F |f (x)| ≤ M . The difference is trivially bounded by 2M , eliminating the need to analyze spatial gradients or invoke transport bounds. For unbounded metrics (W1 , Energy), note that the posterior is µy (dx) =

1 exp(−Φ(x; y))µ(dx) Z 28

Beyond EIG: Stable BOED with IPMs

where Z is the marginal likelihood. Then Z Z f (x)µ̃y (dx) f (x)µy (dx) − X XZ Z 1 1 = f (x) exp(−Φ(x; y))µ̃(dx) f (x) exp(−Φ(x; y))µ(dx) − Z X Z̃ X For convenience, let I and I˜ represent those unnormalized integrals: Z

f (x)µ (dx) −

f (x)µ̃y (dx) =

I I˜ 1 I˜ − = I − I˜ + Z̃ − Z Z Z Z̃ Z Z̃

˜ as ∆I(y). The second term (Z̃ − Z) is the difference in We denote the first term (I − I) unnormalized integrals where the test function is a constant f (x) = 1 (Note that R the multi˜ plier |I| in the second term is trivially bounded by the prior’s first moment ωF (x)µ̃(dx), which is guaranteed to be finite under the higher-order moment assumptions required to bound ∆I(f )). Then the core mathematical challenge lies in bounding the difference in the unnormalized integrals: Z Z f (x) exp(−Φ(x; y))µ(dx) − f (x) exp(−Φ(x; y))µ̃(dx) . ∆I(y) = X

Let hy (x) = f (x) exp(−Φ(x; y)) denote the effective integrand. We can rewrite the difference using any joint coupling p(x, x∗ ) whose marginals are µ and µ̃: Z ∆I(y) =

hy (x) − hy (x∗ ) p(dx, dx∗ ) .

X ×X

By the Mean Value Theorem, for every pair (x, x∗ ), there exists a point z on the line segment between them such that: |hy (x) − hy (x∗ )| ≤ ∥∇hy (z)∥∥x − x∗ ∥. This yields the transport bound: Z ∆I(y) ≤

∥∇hy (z)∥∥x − x∗ ∥p(dx, dx∗ ).

X ×X

Applying the product rule, its gradient is: ∇hy (x) = ∇f (x) exp(−Φ(x; y)) − f (x) exp(−Φ(x; y))∇x Φ(x; y) Since exp(−Φ) ≤ 1, the norm is strictly bounded by: ∥∇hy (x)∥ ≤ ∥∇f (x)∥ + |f (x)|∥∇x Φ(x; y)∥ ≤ ∥∇f (x)∥ + ωF (x)∥∇x Φ(x; y)∥. For a Gaussian observation model y = G(x) + η with noise covariance Σ, the gradient of the data misfit potential is ∇x Φ(x; y) = −JG (x)⊤ Σ−1 y − G(x) . If we apply the Lipschitz 29

Wu, Liang and Yang

assumption for the forward model (L1 ), we have ∥JG (x)∥ ≤ L1 and ∥G(x)∥ ≤ L1 ∥x∥ + ∥G(0)∥. Substituting these into the gradient bound yields: ∥∇x Φ(x; y)∥ ≤ L1 ∥Σ−1 ∥ ∥y∥ + ∥G(x)∥

≤ L1 ∥Σ−1 ∥ L1 ∥x∥ + ∥y∥ + ∥G(0)∥ = O(L21 ∥x∥) + Cy . Because ∥∇f (x)∥ is locally integrable or strictly bounded depending on the metric, the tail behavior as ∥x∥ → ∞ is strictly dominated by the product ωF (x)∥∇x Φ(x; y)∥: • Energy Distance: With the energy distance, the test function envelope grows sublinearly, giving us ωE (x) = O(∥x∥1/2 ). Because of this, the gradient scales as: ∥∇hy (x)∥ ≈ O(∥x∥1/2 ) × O(L21 ∥x∥) = O(L21 ∥x∥1.5 ). R When we plug this into the transport bound ∥∇hy (z)∥∥x − x∗ ∥p(dx, dx∗ ), the resulting integrand scales as O(∥x∥2.5 ). In practice, this means the prior only needs a finite 2.5-th moment to keep the discrepancy finite, which gives us the flexibility to use heavy-tailed priors with polynomial decay. • Wasserstein-1: For Wasserstein-1, the test function envelope grows linearly, meaning ωW1 (x) = O(∥x∥). This pushes the gradient scaling to a quadratic rate: ∥∇hy (x)∥ ≈ O(∥x∥) × O(L21 ∥x∥) = O(L21 ∥x∥2 ). As a result, the transport integrand scales as O(∥x∥3 ). However, the real constraint comes from the overarching expected utility bound, which requires us to integrate D(y) over all possible observations y ∈ Y. Under extreme observation shifts, the polynomial warping inside D(y) scales with exp(O(L21 ∥x∥2 )). To keep this bound firmly under control, we are forced to use a prior with exponential tail decay. This naturally leads to the standard assumption of a strictly sub-Gaussian prior µ(x) ∝ exp(−L2 ∥x∥2 ), where the geometric condition L21 < CL2 ensures the prior decays fast enough to overpower the forward model’s warping.

References Alen Alexanderian. Optimal experimental design for infinite-dimensional Bayesian inverse problems governed by PDEs: A review. Inverse Problems, 37(4):043001, 2021. Alen Alexanderian, Noemi Petra, Georg Stadler, and Omar Ghattas. A fast and scalable method for A-optimal design of experiments for infinite-dimensional Bayesian nonlinear inverse problems. SIAM Journal on Scientific Computing, 38(1):A243–A272, 2016. Brandon Amos, Lei Xu, and J Zico Kolter. Input convex neural networks. In International Conference on Machine Learning, pages 146–155. PMLR, 2017. 30

Beyond EIG: Stable BOED with IPMs

Martin Arjovsky, Soumith Chintala, and Léon Bottou. Wasserstein generative adversarial networks. In International Conference on Machine Learning, pages 214–223. PMLR, 2017. David Barber and Felix Agakov. The IM algorithm: A variational approach to information maximization. Advances in Neural Information Processing Systems, 16(320):201, 2004. Joakim Beck, Ben Mansour Dia, Luis FR Espath, Quan Long, and Raul Tempone. Fast Bayesian experimental design: Laplace-based importance sampling for the expected information gain. Computer Methods in Applied Mechanics and Engineering, 334:523–553, 2018. Alain Berlinet and Christine Thomas-Agnan. Reproducing kernel Hilbert spaces in probability and statistics. Springer Science & Business Media, 2011. José M Bernardo. Expected information as expected utility. The Annals of Statistics, pages 686–690, 1979. Gantavya Bhatt, Yifang Chen, Arnav Das, Jifan Zhang, Sang Truong, Stephen Mussmann, Yinglun Zhu, Jeff Bilmes, Simon Du, Kevin Jamieson, et al. An experimental design framework for label-efficient supervised finetuning of large language models. In Findings of the Association for Computational Linguistics: ACL 2024, pages 6549–6560, 2024. Jeremiah Birrell, Paul Dupuis, Markos A Katsoulakis, Yannis Pantazis, and Luc Rey-Bellet. (f, Γ)-divergences: Interpolating between f-divergences and integral probability metrics. Journal of Machine Learning Research, 23(39):1–70, 2022. Lorenzo Brasco, David Gómez-Castro, and Juan Luis Vázquez. Characterisation of homogeneous fractional Sobolev spaces. Calculus of Variations and Partial Differential Equations, 60(2):60, 2021. Romain Camilleri, Kevin Jamieson, and Julian Katz-Samuels. High-dimensional experimental design and kernel bandits. In International Conference on Machine Learning, pages 1227–1237. PMLR, 2021. Zixuan Cang, Qing Nie, and Yanxiang Zhao. Supervised optimal transport. SIAM Journal on Applied Mathematics, 82(5):1851–1877, 2022. Kathryn Chaloner and Isabella Verdinelli. Bayesian experimental design: A review. Statistical Science, pages 273–304, 1995. Jinyuan Chang, Chenguang Duan, Yuling Jiao, Ruoxuan Li, Jerry Zhijian Yang, and Cheng Yuan. Provable diffusion posterior sampling for Bayesian inversion. arXiv preprint arXiv:2512.08022, 2025. Ke Chen, Haizhao Yang, and Chugang Yi. Data completion for electrical impedance tomography by conditional diffusion models. arXiv preprint arXiv:2602.07813, 2026. Peng Chen, Michael R Haberman, and Omar Ghattas. Optimal design of acoustic metamaterial cloaks under uncertainty. Journal of Computational Physics, 431:110114, 2021. 31

Wu, Liang and Yang

Shizhe Chen, Daniela M Witten, and Ali Shojaie. Selection and estimation for mixed graphical models. Biometrika, 102(1):47–64, 2015. Eleonora Di Nezza, Giampiero Palatucci, and Enrico Valdinoci. Hitchhiker’s guide to the fractional Sobolev spaces. Bulletin Des Sciences Mathématiques, 136(5):521–573, 2012. Hovav A Dror and David M Steinberg. Sequential experimental designs for generalized linear models. Journal of the American Statistical Association, 103(481):288–298, 2008. Richard M Dudley. Real analysis and probability. Chapman and Hall/CRC, 2018. Oliver RA Dunbar, Michael F Howland, Tapio Schneider, and Andrew M Stuart. Ensemblebased experimental design for targeting data acquisition to inform climate models. Journal of Advances in Modeling Earth Systems, 14(9):e2022MS002997, 2022. Tanner Fiez, Lalit Jain, Kevin G Jamieson, and Lillian Ratliff. Sequential experimental design for transductive linear bandits. Advances in Neural Information Processing Systems, 32, 2019. Adam Foster, Martin Jankowiak, Elias Bingham, Paul Horsfall, Yee Whye Teh, Thomas Rainforth, and Noah Goodman. Variational Bayesian optimal experimental design. Advances in Neural Information Processing Systems, 32, 2019. Adam Foster, Desi R Ivanova, Ilyas Malik, and Tom Rainforth. Deep adaptive design: Amortizing sequential bayesian experimental design. In International Conference on Machine Learning, pages 3384–3395. PMLR, 2021. Erdun Gao, Liang Zhang, Jake Fawkes, Aoqi Zuo, Wenqin Liu, Haoxuan Li, Mingming Gong, and Dino Sejdinovic. Observationally informed adaptive causal experimental design. arXiv preprint arXiv:2603.03785, 2026. Loukas Grafakos and Seungly Oh. The kato-ponce inequality. Communications in Partial Differential Equations, 39(6):1128–1157, 2014. Arthur Gretton, Kenji Fukumizu, Zaid Harchaoui, and Bharath K Sriperumbudur. A fast, consistent kernel two-sample test. Advances in Neural Information Processing Systems, 22, 2009. Arthur Gretton, Karsten M Borgwardt, Malte J Rasch, Bernhard Schölkopf, and Alexander Smola. A kernel two-sample test. The Journal of Machine Learning Research, 13(1):723– 773, 2012. Eldad Haber, Lior Horesh, and Luis Tenorio. Numerical methods for experimental design of large-scale linear ill-posed inverse problems. Inverse Problems, 24(5):055012, 2008. Leonid G Hanin. Kantorovich-Rubinstein norm and its application in the theory of Lipschitz spaces. Proceedings of the American Mathematical Society, 115(2):345–352, 1992. Yang-Bo He and Zhi Geng. Active learning of causal networks with intervention experiments and optimal designs. Journal of Machine Learning Research, 9(11), 2008. 32

Beyond EIG: Stable BOED with IPMs

Tapio Helin, Nuutti Hyvonen, and Juha-Pekka Puska. Edge-promoting adaptive Bayesian experimental design for X-ray imaging. SIAM Journal on Scientific Computing, 44(3): B506–B530, 2022. Tapio Helin, Youssef Marzouk, and Jose Rodrigo Rojo-Garcia. Bayesian optimal experimental design with Wasserstein information criteria. arXiv preprint arXiv:2504.10092, 2025. Kathrin Hellmuth, Ruhui Jin, Qin Li, and Stephen J Wright. Data selection: at the interface of PDE-based inverse problem and randomized linear algebra. arXiv preprint arXiv:2510.01567, 2025. Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. Advances in Neural Information Processing Systems, 33:6840–6851, 2020. Zixuan Hu, Li Shen, Zhenyi Wang, Yongxian Wei, and Dacheng Tao. Adaptive defense against harmful fine-tuning for large language models via Bayesian data scheduler. In The Thirty-ninth Annual Conference on Neural Information Processing Systems, 2025. Xun Huan and Youssef M Marzouk. Simulation-based optimal Bayesian experimental design for nonlinear systems. Journal of Computational Physics, 232(1):288–317, 2013. Xun Huan, Jayanth Jagalur, and Youssef Marzouk. Optimal experimental design: Formulations and computations. Acta Numerica, 33:715–840, 2024. Ruhui Jin, Qin Li, Stephen O Mussmann, and Stephen J Wright. Continuous nonlinear adaptive experimental design with gradient flow. arXiv preprint arXiv:2411.14332, 2024. Gavin Kerrigan, Christian A Naesseth, and Tom Rainforth. A geometric approach to optimal experimental design. arXiv preprint arXiv:2510.14848, 2025. Steven Kleinegesse and Michael U Gutmann. Efficient bayesian experimental design for implicit models. In The 22nd International Conference on Artificial Intelligence and Statistics, pages 476–485. PMLR, 2019. Alexander Korotin, Vage Egiazarian, Arip Asadulaev, Alexander Safin, and Evgeny Burnaev. Wasserstein-2 generative networks. arXiv preprint arXiv:1909.13082, 2019. Peter M Lee. Bayesian statistics. Oxford University Press London, 1989. Fengyi Li, Ricardo Baptista, and Youssef Marzouk. Expected information gain estimation via density approximations: Sample allocation and dimension reduction. arXiv preprint arXiv:2411.08390, 2024. Ling Liang and Haizhao Yang. PNOD: An efficient projected Newton framework for exact optimal experimental designs. https://arxiv.org/abs/2409.18392, 2024. Dennis V Lindley. On a measure of the information provided by an experiment. The Annals of Mathematical Statistics, 27(4):986–1005, 1956. 33

Wu, Liang and Yang

Yaron Lipman, Ricky TQ Chen, Heli Ben-Hamu, Maximilian Nickel, and Matt Le. Flow matching for generative modeling. arXiv preprint arXiv:2210.02747, 2022. Pangpang Liu and Yichuan Zhao. A review of recent advances in empirical likelihood. Wiley Interdisciplinary Reviews: Computational Statistics, 15(3):e1599, 2023. Shiao Liu, Xingyu Zhou, Yuling Jiao, and Jian Huang. Wasserstein generative learning of conditional distribution. arXiv preprint arXiv:2112.10039, 2021. Quan Long, Marco Scavino, Raúl Tempone, and Suojin Wang. Fast estimation of expected information gains for Bayesian experimental designs based on Laplace approximations. Computer Methods in Applied Mechanics and Engineering, 259:24–39, 2013. Turab Lookman, Prasanna V Balachandran, Dezhen Xue, and Ruihao Yuan. Active learning in materials science with emphasis on adaptive sampling using uncertainties for targeted design. npj Computational Materials, 5(1):21, 2019. David JC MacKay. Information-based objective functions for active data selection. Neural computation, 4(4):590–604, 1992. Ashok Makkuva, Amirhossein Taghvaei, Sewoong Oh, and Jason Lee. Optimal transport mapping via input convex neural networks. In International Conference on Machine Learning, pages 6672–6681. PMLR, 2020. William Muldrew, Peter Hayes, Mingtian Zhang, and David Barber. Active preference learning for large language models. In International Conference on Machine Learning, pages 36577–36590. PMLR, 2024. Alfred Müller. Integral probability metrics and their generating classes of functions. Advances in Applied Probability, 29(2):429–443, 1997. Jay I Myung, Daniel R Cavagnaro, and Mark A Pitt. A tutorial on adaptive design optimization. Journal of Mathematical Psychology, 57(3-4):53–67, 2013. Ilsang Ohn and Lizhen Lin. Adaptive variational Bayes: Optimality, computation and applications. The Annals of Statistics, 52(1):335–363, 2024. Antony M Overstall. Properties of Fisher information gain for Bayesian design of experiments. Journal of Statistical Planning and Inference, 218:138–146, 2022. Gabriel Peyré and Marco Cuturi. Computational optimal transport: With applications to data science. Now Foundations and Trends, 2019. Dennis Prangle, Sophie Harbisher, and Colin S Gillespie. Bayesian experimental design without posterior calculations: An adversarial approach. Bayesian Analysis, 18(1):133– 163, 2023. Rafael Rafailov, Archit Sharma, Eric Mitchell, Christopher D Manning, Stefano Ermon, and Chelsea Finn. Direct preference optimization: Your language model is secretly a reward model. Advances in Neural Information Processing Systems, 36:53728–53741, 2023. 34

Beyond EIG: Stable BOED with IPMs

Howard Raiffa and Robert Schlaifer. Applied statistical decision theory. John Wiley & Sons, 2000. Tom Rainforth, Rob Cornish, Hongseok Yang, Andrew Warrington, and Frank Wood. On nesting Monte Carlo estimators. In International Conference on Machine Learning, pages 4267–4276. PMLR, 2018. Tom Rainforth, Adam Foster, Desi R Ivanova, and Freddie Bickford Smith. Bayesian experimental design. Statistical Science, 39(1):100–114, 2024.

Modern

Lars Ruthotto, Julianne Chung, and Matthias Chung. Optimal experimental design for inverse problems with state constraints. SIAM Journal on Scientific Computing, 40(4): B1080–B1100, 2018. Elizabeth G Ryan, Christopher C Drovandi, James M McGree, and Anthony N Pettitt. A review of modern computational algorithms for Bayesian optimal design. International Statistical Review, 84(1):128–154, 2016. Jerome Sacks, William J Welch, Toby J Mitchell, and Henry P Wynn. Design and analysis of computer experiments. Statistical Science, 4(4):409–423, 1989. Paola Sebastiani and Henry P Wynn. Maximum entropy sampling and optimal Bayesian experimental design. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 62(1):145–157, 2000. Dino Sejdinovic, Bharath Sriperumbudur, Arthur Gretton, and Kenji Fukumizu. Equivalence of distance-based and RKHS-based statistics in hypothesis testing. The Annals of Statistics, pages 2263–2291, 2013. Claude Elwood Shannon. A mathematical theory of communication. The Bell System Technical Journal, 27(3):379–423, 1948. Michael C Shewry and Henry P Wynn. Maximum entropy sampling. Journal of Applied Statistics, 14(2):165–170, 1987. Alex Smola, Arthur Gretton, Le Song, and Bernhard Schölkopf. A Hilbert space embedding for distributions. In International Conference on Algorithmic Learning Theory, pages 13– 31. Springer, 2007. Shanshan Song, Tong Wang, Guohao Shen, Yuanyuan Lin, and Jian Huang. Wasserstein generative regression. Journal of the Royal Statistical Society Series B: Statistical Methodology, 88(1):330–351, 2026. Ben Sorscher, Robert Geirhos, Shashank Shekhar, Surya Ganguli, and Ari Morcos. Beyond neural scaling laws: Beating power law scaling via data pruning. Advances in Neural Information Processing Systems, 35:19523–19536, 2022. David M Steinberg and William G Hunter. Experimental design: Review and comment. Technometrics, 26(2):71–97, 1984. 35

Wu, Liang and Yang

Robert S Strichartz. Multipliers on fractional Sobolev spaces. Journal of Mathematics and Mechanics, 16(9):1031–1060, 1967. Ke Sun, Linglong Kong, Hongtu Zhu, and Chengchun Shi. Arma-design: Optimal treatment allocation strategies for a/b testing in partially observable time series experiments. arXiv preprint arXiv:2408.05342, 2024. Gábor J Székely and Maria L Rizzo. A new test for multivariate normality. Journal of Multivariate Analysis, 93(1):58–80, 2005. Gábor J Székely, Maria L Rizzo, et al. Testing for equal distributions in high dimension. InterStat, 5(16.10):1249–1272, 2004. Ilya Tolstikhin, Bharath K Sriperumbudur, and Krikamol Muandet. Minimax estimation of kernel mean embeddings. Journal of Machine Learning Research, 18(86):1–47, 2017. Cédric Villani et al. Optimal transport: Old and new, volume 338. Springer, 2009. Andrew Wagenmaker and Kevin G Jamieson. Instance-dependent near-optimal policy identification in linear MDPs via online experiment design. Advances in Neural Information Processing Systems, 35:5968–5981, 2022. Stephen G Walker. Bayesian information in an experiment and the Fisher information distance. Statistics & Probability Letters, 112:5–9, 2016. Di Wu, Ling Liang, and Haizhao Yang. PINS: Proximal iterations with sparse Newton and Sinkhorn for optimal transport. arXiv preprint arXiv:2502.03749, 2025. Keyi Wu, Peng Chen, and Omar Ghattas. A fast and scalable computational framework for large-scale high-dimensional Bayesian optimal experimental design. SIAM/ASA Journal on Uncertainty Quantification, 11(1):235–261, 2023. Sang Michael Xie, Shibani Santurkar, Tengyu Ma, and Percy S Liang. Data selection for language models via importance resampling. Advances in Neural Information Processing Systems, 36:34201–34227, 2023. Yun Yang, Martin J Wainwright, and Michael I Jordan. On the computational complexity of high-dimensional Bayesian variable selection. The Annals of Statistics, 44(6):2497–2532, 2016. Jin Zhu, Jingyi Li, Hongyi Zhou, Yinan Lin, Zhenhua Lin, and Chengchun Shi. Balancing interference and correlation in spatial experimental designs: A causal graph cut approach. arXiv preprint arXiv:2505.20130, 2025.

Record · ID 126504 · SHA-256 494f62dc980c9ea8
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.