Batch-based Bayesian Optimal Experimental Design in Linear Inverse Problems
February 12, 2026
Abstract
Experimental design is central to science and engineering. A ubiquitous challenge is how to maximize the value of information obtained from expensive or constrained experimental settings. Bayesian optimal experimental design (OED) provides a principled framework for addressing such questions. In this paper, we study experimental design problems such as the optimization of sensor locations over a continuous domain in the context of linear Bayesian inverse problems. We focus in particular on batch design, that is, the simultaneous optimization of multiple design variables, which leads to a notoriously difficult non-convex optimization problem. We tackle this challenge using a promising strategy recently proposed in the frequentist setting, which relaxes A-optimal design to the space of finite positive measures. Our main contribution is the rigorous identification of the Bayesian inference problem corresponding to this relaxed A-optimal OED formulation. Moreover, building on recent work, we develop a Wasserstein gradient-flow–based optimization algorithm for the expected utility and introduce novel regularization schemes that guarantee convergence to an empirical measure. These theoretical results are supported by numerical experiments demonstrating both convergence and the effectiveness of the proposed regularization strategy.
1 Introduction↩︎
The design of experiments is a central theme across science and engineering. Data collection is often restricted by practical constraints such as physical restrictions (consider radiation exposure in medical imaging) or high costs (consider seismic imaging), to name a few. The allocation of measurement effort, the positioning of sensors, and the timing of observations are all decisions that shape the quality of inference one can draw from an experiment. This ubiquity has motivated the development of a wide range of approaches, from heuristic rules of thumb to formal optimization-based strategies, with the goal of extracting maximal value from every measurement.
A particularly principled framework for addressing such questions is Bayesian optimal experimental design (OED). In this approach, prior knowledge about the system under study is combined with a probabilistic model for future observations. Moreover, a utility function quantifies the expected benefit of making a particular set of measurements. Thus, Bayesian design naturally accounts for uncertainties in both parameters and data, allowing one to reason explicitly about information gain and decision quality. Standard choices for the optimality condition include the A- and D-optimality criteria, of which we consider the A-optimality [1]. While exact computation of these quantities is often infeasible, advances in computation and approximation methods have enabled Bayesian design to be applied to increasingly complex systems.
In this work, we consider experimental design problems defined over a continuous domain. Our focus is on batch-based optimization, where multiple design variables, such as sensor locations, are optimized simultaneously. Such batch design problems are notoriously nonconvex and exhibit multiple local minima, which makes them challenging to solve reliably. A promising recent strategy to address these difficulties is to convexify the optimization problem by relaxing it to the space of probability measures. This relaxation has been successfully employed in related contexts, both in Bayesian optimization [2] and in frequentist settings [3]–[5].
Our work is particularly closely aligned with the objectives studied in [3]. Let us introduce some notation to make this connection precise. Suppose function \(g\) in a reproducing kernel Hilbert space on the domain \({\mathcal{X}}\) is observed at a point \(x \in {\mathcal{X}}\), connecting it to an unknown function of interest \(f\) through the measurement model \[\label{eq:informal95obs95model} y(x) = (Af)(x) + \epsilon,\tag{1}\] where \(A\) is a linear operator and \(\epsilon\) is additive Gaussian noise with variance \(\sigma^2\), independent of both \(f\) and \(x\). In this setting, [3] investigates, in particular, the frequentist A-optimal measurement strategy, which under suitable assumptions on the measurement model amounts to solving \[\label{eq:informal95Aopt95freq} x_{Aopt} \in \operatorname*{argmin}_{x \in X} \operatorname{Tr}\left[(A_x^*A_x)^{-1}\right],\tag{2}\] where \(A_x\) denotes the operator \(A\) followed by evaluation at \(x\). The main idea in [3] is to relax the optimization task in 2 to the space of probability measures \({\mathcal{P}}({\mathcal{X}})\) by informally setting \[\label{eq:informal95Aopt95freq95meas} \mu_{Aopt} \in \operatorname*{argmin}_{\mu \in {\mathcal{P}}({\mathcal{X}})} \operatorname{Tr}\left[\left(\int_X A_z^*A_z \mu(dz)\right)^{-1}\right],\tag{3}\] where the minimization functional can be shown to be convex with respect to the probability measure \(\mu\). The relaxation can be understood by noting that the formulation 3 coincides with 2 when the search for \(\mu\) is restricted to Dirac measures \(\{\delta_x \in {\mathcal{P}}({\mathcal{X}}) \; | \; x \in {\mathcal{X}}\}\). The key contribution in [3] is to formulate a Wasserstein gradient flow scheme for optimizing \(\mu\) by particle approximations.
The relaxation extends naturally to a batch-based setting. If one seeks to optimize a batch of sensor locations \({x_1,\dots,x_B}\subset {\mathcal{X}}\) in 2 , then the corresponding object in 3 is the measure \(\sum_{j=1}^B \delta_{x_j}\). This suggests extending the relaxation to all positive measures satisfying \(\mu({\mathcal{X}})=B\). The resulting shift, from a highly non-convex search over \(\otimes_{j=1}^B {\mathcal{X}}\) to a convex optimization problem over measures on \({\mathcal{X}}\), is precisely what makes this approach appealing for large-scale multi-parameter design problems.
To develop this approach further, in this paper, we study whether the experimental design setting of 2 and 3 can be extended to a Bayesian framework. When the prior distribution of \(f\) is Gaussian, the posterior induced by 1 remains Gaussian. In this context, the A-optimality principle amounts to minimizing the posterior variance, however, a coherent relaxation of the problem and its interpretation in the Bayesian paradigm is not trivial. This observation motivates the central questions of our work: what observation does the relaxation correspond to in the Bayesian setting and is there a corresponding relaxation of the design problem aligned with 3 ?
1.1 Our Contribution↩︎
In this paper, we rigorously identify the Bayesian inference task associated to the frequentist batch-based A-optimal design problem proposed in [3], when relaxed to finite positive measures with fixed mass and interpreted through the Bayesian paradigm. This characterization is developed in Sections 2 and 3, where
we propose a non-parametric Bayesian inference task based on observing a Gaussian random process indexed over a Hilbert space associated to the relaxed design \(\mu\), and,
we derive the corresponding posterior covariance, which leads to an infinite-dimensional expected utility functional in Section 3.1.
Importantly, in Proposition 2 we prove that this expected utility is concave with respect to the design measure \(\mu\).
Following previous work [3], we introduce a Wasserstein gradient flow for the particle-based optimization of the expected utility. To this end, we derive the first variation and the gradient of the expected utility in Propositions 3 and 4.
While the optimization method is developed for positive measures with fixed mass, another key contribution of this work is the proposal of an alternative approach to accelerate convergence. Specifically, each of the \(B\) design variables is represented by a probability measure \(\mu_j\), \(j=1,\dots,B\), and the optimization problem is lifted back to \({\mathcal{X}}^B\), that is, the optimization is carried out over the tensor product \(\boldsymbol{\mu}= \bigotimes_{j=1}^B \mu_j\). This formulation enables us to regularize the expected utility through interplay of the design measures \(\mu_j\).
We introduce a two-fold regularization strategy, whose relative influence can be tuned by reweighting:
First, we penalize the expected utility by the empirical variance of each ensemble in order to ensure concentration toward a single-point measure.
Second, we penalize proximity between distinct ensembles via the negative maximum mean discrepancy, thereby discouraging multiple ensembles from collapsing to the same point measure.
While convergence to multiple point measures, or collapse to a single joint point measure, may be justified by the inference task and the resulting information gain, these regularization mechanisms provide practical tools for enforcing additional constraints arising in experimental settings, such as prescribed minimum distances between sensor locations. These effects are demonstrated through numerical examples in Section 6.
1.2 Literature overview↩︎
Bayesian experimental design is a mature research area supported by a substantial and evolving literature. Recent survey-style contributions include [6]–[8]. This work focuses on A-optimality criterion which is a classical utility construction discussed e.g. in [1]. We also highlight the recent work to promote Wasserstein distance based information criterion in [9], motivated by the general notion of a valid information measure proposed by Ginebra in [10].
Inverse problems form a class of high-dimensional inference tasks in which unknown parameters are linked to observations through complex forward models, often governed by partial differential equations. The requirement that experimental design criteria remain well defined across refinement levels has motivated extensions of classical Bayesian optimal experimental design to infinite-dimensional settings; for example, expected information gain and A-optimality criterion admit rigorous nonparametric formulations as developed in [11]. See also the review [12].
The joint optimization of multiple sensor locations prior to any data acquisition has been widely studied in the literature. Owing to the highly nonconvex nature of the resulting design problems, many existing approaches rely on greedy selection procedures, which consequently provide limited guarantees of global optimality. For recent contributions in the setting of general PDE-based inverse problems, see [13]–[15]. In particular, Bayesian sensor placement for electrical impedance tomography has been investigated in [16], [17]. This line of work should be distinguished from a different sequential paradigm in which data are collected after each newly selected sensor location, see e.g. [18], [19].
Gradient-flow–based approaches to optimal experimental design in the frequentist setting have been investigated in [3]–[5], where the optimization problem is formulated over the space of probability measures. In particular, [5] studies Wasserstein gradient flows for E-optimal design in regression models, while [3] and [4] examine A- and D-optimal continuous designs, focusing respectively on linear and nonlinear
settings. Separately, [2] explores the use of Wasserstein and Stein gradient flows for batch Bayesian optimization aimed at selecting
multiple evaluation points in parallel. For more broad discussion on Wasserstein gradient flows, we refer the reader to [20],
[21]. Optimization problems over probability measures via gradient flows have also been considered for general machine learning tasks, see [22].
Section 2 introduces the Bayesian inference task giving rise to the relaxed design problem over finite measures. The corresponding Bayesian OED task and the A-optimal expected utility are studied in Section 3. Section 4 introduces the Wasserstein gradient flow for the optimization to the expected utility, while Section 5 is devoted to
the practical regularization strategies in batch-based optimization. Finally, Section 6 illustrates the method for two inverse problems with one-dimensional Poisson and two-dimensional time-harmonic
Schrödinger equation, with a sensitivity study regarding the regularization in Section 6.3.
2 Relaxation of the Bayesian inference problem↩︎
Let \(H_k\) be a separable reproducing kernel Hilbert space on a bounded smooth domain \({\mathcal{X}}\subset \mathbb{R}^n\), generated by a bounded kernel \(k\) satisfying the required regularity conditions. Since the kernel remains fixed throughout the paper, we write \(H = H_k\) for brevity. We consider optimizing the sensor location \(x\in {\mathcal{X}}\) for recovering \(f\in U\) in a separable Hilbert space \(U\) from \[g(x) = \langle K_x, A f\rangle_{H} + \epsilon\in \mathbb{R}\] where \(A: U \to H\) is linear and bounded, \(K_x\in H\) is the reproducing element and \(\epsilon\) is normally distributed noise independent of the sensor location.
To relax the design problem, let \(\mu\) be a positive finite measure on \({\mathcal{X}}\), and for brevity write \(L^2(\mu) := L^2({\mathcal{X}}; \mu)\). Also, let us introduce the embedding operator \(\iota_\mu: H_k \to L^2(\mu)\). In what follows, we study an inference problem, where the observation \(g_\mu\) is a Gaussian random process indexed by \(L^2(\mu)\) satisfying \[\label{eq:relaxed95inference} g_\mu(h) = \langle A_\mu f, h\rangle_{L^2(\mu)} + \epsilon_\mu(h),\tag{4}\] where we write \(A_\mu = \iota_\mu A : U \to L^2(\mu)\) and \(\epsilon_\mu\) is a isonormal Gaussian process over \(L^2(\mu)\) [23].
The non-parametric Bayesian inference of the signal \(f \in U\) is well-understood [24]. In particular, for a Gaussian prior distribution with a linear Hilbert-space valued observational model, the posterior distribution was established in [25] and has been revisited in many subsequent work (see e.g. [26], [27]). While the random process formulation such as in 4 with signal contaminated with white noise has received less attention, it is well-known that there exists a regular version of the conditional distribution of \(f\) given \(g_\mu\) and that this conditional distribution is Gaussian [28].
For the reader’s convenience we recall a standard construction of the posterior covariance in Theorem 1 that simplifies the setting of [28]. Starting from the Hilbert space \(L^2(\mu)\), one introduces a Hilbert scale \(Z^\mu \subset L^2(\mu) \subset (Z^\mu)'\), where \(Z^\mu\) is a separable Hilbert space that embeds continuously and densely into \(L^2(\mu)\) with the property that realizations of \(\epsilon_\mu\) are ‘contained’ in \((Z^\mu)'\). This procedure is a Hilbert-space version of the abstract Wiener space construction that is standard in modern Gaussian measure theory [29].
More concretely, let \(\{e_j\}_{j\in\mathbb{J}}\) be the orthonormal basis of \(L^2(\mu)\), where \(\mathbb{J}\) is either finite or countable. We define \[Z_s^\mu = \left\{z = \sum_{j\in\mathbb{J}} z_j e_j \; \bigg| \; \sum_{j\in\mathbb{J}} j^{2s} z_j^2 < \infty \right\}, \quad s\in \mathbb{R},\] with the interpretation \(Z_0^\mu = L^2(\mu)\) and \((Z_s^\mu)' = Z_{-s}^\mu\). For the moment, suppose \(\mathbb{J}=\mathbb{N}\). We immediately observe that the random series \(\sum_{j=1}^J \xi_j e_j\), \(\xi_j\in{\mathcal{N}}(0,1)\) i.i.d., converges to a well-defined random variable \(\tilde{\epsilon}_\mu\) in \(L^2(\Omega; Z_{-s}^\mu)\) such that \[\epsilon_\mu(z) = \langle \tilde{\epsilon}_\mu, z\rangle_{Z_{-s}^\mu\times Z_s^\mu},\] in distribution for any \(z\in Z_s^\mu\), when \(s\) is large enough, i.e., \(\sum_{j\in\mathbb{J}} j^{-2s}<\infty.\) Clearly, for a finite index the claim is trivial.
We can now consider the joint distribution \((\tilde{g}_\mu, f)\) in the product space \(Z_{-s}^\mu \times U\), where \(\tilde{g}_\mu\) conditioned on \(f\) has the distribution of \(A_\mu f + \tilde{\epsilon}_\mu\) and \(s\) is fixed large enough. With the abuse of notation, we will identify \(g_\mu\) and \(\tilde{g}_\mu\) (as well as, \(\epsilon_\mu\) and \(\tilde{\epsilon}_\mu\)) in what follows. Moreover, below, we write \(A_\mu^*\) for the adjoint of \(A_\mu : U \to L^2(\mu)\), while \(A_\mu'\) denotes the adjoint of \(A_\mu : U \to Z_{-s}^\mu\), where the latter adjoint is obtained the canonical embedding of the range \(L^2(\mu)\) into \(Z_{-s}^\mu\).
Theorem 1. Suppose the random variable \((g_\mu, f) \in (Z_{-s}^\mu, U)\) has the joint distribution induced by the marginal \(f\in {\mathcal{N}}(0, C_f)\) on \(U\) and the likelihood induced by 4 . Then the conditional distribution of \(f\) given \(g_\mu\) has a regular Gaussian version with a trace-class covariance operator \[\label{eq:posterior95covariance95form1} C_{post}^\mu = C_f - C_f A_\mu^* \left(A_\mu C_f A_\mu^* + I\right)^{-1} A_\mu C_f : U\to U.\tag{5}\]
Proof. The regular conditional distribution of \(f\) given \(g_\mu\) is Gaussian with the covariance \(C_{post}^\mu : U\to U\) given by \[\label{eq:posterior95covariance95general} C_{post}^\mu = C_f - C_{f, g_\mu} C_{g_\mu}^{-1}C_{g_\mu, f},\tag{6}\] where \(C_{f,f'} = \mathbb{E}\left[f \otimes f'\right]\) stands for the cross-covariance of random variables \(f\) and \(f'\), and \(C_g\) is the covariance operator of \(g_\mu\) [25]. While \(C_g\) is not invertible in general, the product \(C_{f, g_\mu} C_{g_\mu}^{-\frac{1}{2}}\) is well-defined as a Hilbert–Schmidt operator, see e.g. [27].
It is straightforward to observe that \[C_{f, g_\mu} = C_f A_\mu' : Z_{-s}^\mu \to U \quad \text{and} \quad C_{g_\mu, f} = A_\mu C_f : U \to Z_{-s}^\mu.\] Furthermore, due to the Hilbert scale construction, there exists a linear bounded bijection \(R : Z_s^\mu \to Z_{-s}^\mu\) induced in general by the Riesz representation such that \[\langle a, b\rangle_{Z_{-s}^\mu \times Z_s^\mu} = \langle a, b \rangle_{Z_0^\mu} = \langle a, R b\rangle_{Z_{-s}^\mu},\] for any \(a\in Z_0^\mu\) and \(b\in Z_s^\mu\). It directly follows that \(A_\mu' = A_\mu^* R^{-1} : Z_{-s}^\mu \to U\).
Next, since \(f\) and \(\epsilon_\mu\) are independent, we have for any \(a,b\in Z_{-s}^\mu\) that \[\begin{align} \langle C_g a, b\rangle_{Z_{-s}^\mu} & = & \mathbb{E}\langle A_\mu u, a\rangle_{Z_{-s}^\mu} \langle A_\mu u, b\rangle_{Z_{-s}^\mu} + \mathbb{E}\langle \epsilon_\mu, a\rangle_{Z_{-s}^\mu} \langle \epsilon_\mu, b\rangle_{Z_{-s}^\mu} \\ & = & \langle C_f A_\mu' a, A_\mu' b\rangle_{U} + \mathbb{E}\langle \epsilon_\mu, R^{-1} a\rangle_{Z_{-s}^\mu\times Z_s^\mu} \langle \epsilon_\mu, R^{-1} b\rangle_{Z_{-s}^\mu \times Z_s^\mu} \\ & = & \langle A_\mu C_f A_\mu' a, b\rangle_{Z_{-s}^\mu} + \langle R^{-1} a, R^{-1} b\rangle_{Z_0^\mu} \\ & = & \langle A_\mu C_f A_\mu' a, b\rangle_{Z_{-s}^\mu} + \langle R^{-1} a, b\rangle_{Z_{-s}^\mu}. \end{align}\] In consequence, we have \[C_g = A_\mu C_f A_\mu' + R^{-1} : Z_{-s}^\mu \to Z_{-s}^\mu.\] Let us next show that the operator 6 has the same action on \(U\) as the operator 5 . Clearly, in the latter operator, the inverse \(\left(A_\mu C_f A_\mu^* + I\right)^{-1}\) is well-defined linear and bounded operator on \(Z_0^\mu\). To show that the actions coincide, let \(u\in {\mathcal{R}}(A_\mu C_f) \subset Z_0^\mu\). It remains to show that \[R^{-1} \left[\left(A_\mu C_f A_\mu^* + I\right)R^{-1}\right]^{-1} u = \left(A_\mu C_f A_\mu^* + I\right)^{-1} u.\] Notice that due to our Hilbert scale construction, \(R\) is a diagonal operator in the orthonormal basis \(\{e_j\}_{j\in\mathbb{J}}\), which extends to a bijective mapping between \(Z_t^\mu\) and \(Z_{t-2s}^\mu\) for any \(t\in\mathbb{R}\). It follows that there exists a unique vector \(w \in Z_{2s}^\mu\) such that \[(A_\mu C_f A_\mu^* + I)R^{-1} w = u \in Z_0^\mu.\] In consequence, we have \[R^{-1}w = (A_\mu C_f A_\mu^* + I)^{-1} u\] which proves the claim. ◻
To better characterize the self-adjoint operator \(A_\mu^*A_\mu\), observe that the adjoint \(\iota_\mu^*:L^2({\mathcal{X}};\mu)\to H\) satisfies \[\langle \iota_\mu^* f, g\rangle_{H} = \langle f, \iota_\mu g\rangle_{L^2(\mu)} = \int_{{\mathcal{X}}} f(x)\,\langle k_x,g\rangle_H\,\mu(dx) = \left\langle \int_{{\mathcal{X}}} f(x)\,k_x\,\mu(dx),\, g\right\rangle_H ,\] and therefore \[\label{eq:iotamuadj} \iota_\mu^* f = \int_{{\mathcal{X}}} k_x f(x)\,\mu(dx).\tag{7}\] The properties of the integral operator \(\iota_\mu^*\) are well known; see, for instance, [30]. In particular, \(\iota_\mu^*\) is bounded Hilbert–Schmidt operator with \(\left\|\iota_\mu^*\right\|_{HS} = \left\|k\right\|_{L^2(\mu)}\) [30]. To this vein, we also introduce the operator \[B_\mu := \int_{{\mathcal{X}}} k_x \otimes k_x \,\mu(dx) : H \to H,\] where the tensor product is understood in \(H\). The operator \(B_\mu\) is often called the covariance (or second-moment) operator in statistical learning theory. It is compact, positive, self-adjoint, and trace class , and satisfies \(B_\mu=\iota_\mu^*\iota_\mu\) with \(\operatorname{Tr}_H\, B_\mu = \left\|k\right\|_{L^2(\mu)}\) [30]. Consequently, \[\label{eq:normal95operator} A_\mu^* A_\mu = A^* B_\mu A = \int_{{\mathcal{X}}} (A^*k_x)\otimes(A^*k_x)\,\mu(dx) : U \to U\tag{8}\] is a bounded self-adjoint operator inheriting the properties of \(B_\mu\).
Corollary 1. The posterior covariance in 5 satisfies \[\label{eq:posterior95covariance95form2} C_{post}^\mu = C_f^{\frac{1}{2}} (I + C_f^{\frac{1}{2}} A_\mu^* A_\mu C_f^{\frac{1}{2}})^{-1}C_f^{\frac{1}{2}} : U \to U.\tag{9}\]
Proof. The operator defined by 9 is a linear bounded, self-adjoint and trace-class operator in \(U\), inheriting these properties from the prior covariance \(C_f\) and boundedness of \(A_\mu\). To see that the two operators coincide, abbreviate \(T_\mu = A_\mu C_f^{\frac{1}{2}} : U\to L^2(\mu)\) and observe \[(I + T_\mu^* T_\mu)^{-1} = I - T_\mu^* (T_\mu T_\mu^* + I)^{-1} T_\mu.\] This proves the claim. ◻
Example 1 (An empirical measure). Consider the inference task with design induced by the empirical measure \(\mu_B = \frac{1}{B} \sum_{j=1}^B \delta_{x_j} \in {\mathcal{P}}({\mathcal{X}})\). It follows that \[C_{post}^{\mu_B} = C_f^{\frac{1}{2}}\left(\frac{1}{B} \sum_{j=1}^B (C_f^{\frac{1}{2}}A^* k_{x_j}) \otimes (C_f^{\frac{1}{2}}A^* k_{x_j}) + I\right)^{-1} C_f^{\frac{1}{2}}.\] This expression coincides with the posterior covariance obtained from the observation model \[g_j = (Af)(x_j) + \sqrt{B} \epsilon_j, \quad, j=1,\ldots, B,\] where \(\epsilon_j\sim {\mathcal{N}}(0,1)\) i.i.d. Thus, the empirical measure \(\mu_B\) corresponds to a design in which independent observations are collected at the sensor locations \(x_j\), but each measurement is corrupted by noise whose standard deviation grows as \(\sqrt{B}\). In effect, spreading the sampling budget across \(B\) locations leads to a higher noise-to-signal ratio at each individual observation. Similarly, starting with \(\tilde{\mu}_B = \sum_{j=1}^B \delta_{x_j}\) leads to independent observations at same locations but with fixed noise-to-signal ratio.
3 Bayesian OED in the relaxed setting↩︎
3.1 A-optimality criterion↩︎
Bayesian optimal experimental design seeks to maximize a utility functional over the design space. In this work, the design space is the set of positive measures \(\mu\) satisfying \(\mu({\mathcal{X}})=B\), which we henceforth denote by \({\mathcal{P}}^B({\mathcal{X}})\) with the convention that \({\mathcal{P}}({\mathcal{X}}) = {\mathcal{P}}^1({\mathcal{X}})\). The optimal design is formalized as the maximizer \[\begin{align} \mu^* = \operatorname*{argmax}_{\mu \in {\mathcal{P}}^B({\mathcal{X}})}\, U(\mu) = \operatorname*{argmax}_{\mu \in {\mathcal{P}}^B({\mathcal{X}})}\, \mathbb{E}\, u(f,g^\mu; \mu), \end{align}\] where the expectation is taken over the joint Bayesian distribution of \(f\) and \(g^\mu\) constructed from the Gaussian prior on \(f\) and the likelihood in 4 . In this paper, we focus on the A-optimality criterion that emerges from the choice \[u(f,g^\mu; \mu) = - \left\|f - \hat{f}(g^\mu; \mu)\right\|^2,\] where \(\hat{f}(g^\mu; \mu)\) stands for the mean of the posterior distribution. This leads to the well-known formula for the expected utility, namely, \(U\) satisfies \[\begin{align} \label{eq:Umu} U(\mu) = - \operatorname{Tr}_U\, \left[C_{post}(\mu)\right] & = & - \operatorname{Tr}_U\, \left[C_f^{\frac{1}{2}} (I + C_f^{\frac{1}{2}} A_\mu^* A_\mu C_f^{\frac{1}{2}})^{-1}C_f^{\frac{1}{2}}\right] \nonumber\\ & = & - \operatorname{Tr}_{C_f^{1/2}(U)}\, \left(\int_{\mathcal{X}}(C_f^{\frac{1}{2}}A^* k_x) \otimes (C_f^{\frac{1}{2}} A^* k_x) \mu(dx) + I\right)^{-1} \end{align}\tag{10}\] by identity 8 . In what follows, we write \(U_B(\mu)=U(\mu)\) to emphasize the dependence on the domain \({\mathcal{P}}^B({\mathcal{X}})\). Motivated by the expression in 10 , let us abbreviate \[\label{eq:Smu} S_\mu = \int_{\mathcal{X}}(C_f^{\frac{1}{2}}A^* k_x) \otimes (C_f^{\frac{1}{2}} A^* k_x) \mu(dx).\tag{11}\] The operator \(S_\mu\) can be interpreted as the covariance \(B_\mu\) preconditioned by \(AC_f^{\frac{1}{2}}\).
Proposition 2. The mapping \(U_B : {\mathcal{P}}^B({\mathcal{X}})\to \mathbb{R}\) defined by 10 is concave.
Proof. Fix \(\mu_0,\mu_1\in{\mathcal{P}}^B({\mathcal{X}})\) and \(t\in[0,1]\), and set \[\mu_t:=(1-t)\mu_0+t\mu_1,\qquad M_t:=S_{\mu_t}+I.\] By linearity it follows that \(M_t=(1-t)M_0+tM_1\). Each \(M_i\) is bounded, self-adjoint, and strictly positive, hence invertible, and the same holds for \(M_t\).
We now use the operator convexity (in terms of the Loewner order) of the inverse on the cone of strictly positive bounded self-adjoint operators, i.e., it holds that \[\label{eq:operator-convex-inverse} \bigl((1-t)M_0+tM_1\bigr)^{-1}\;\leq (1-t)M_0^{-1}+tM_1^{-1}\tag{12}\] for \(0\le t\le 1\). Since the trace is monotone with respect to the Loewner order on positive trace-class operators [31], applying \(\operatorname{Tr}_{C_f^{1/2}(U)}(\cdot)\) to 12 yields \[\operatorname{Tr}_{C_f^{1/2}(U)}(M_t^{-1}) \;\le\;(1-t)\operatorname{Tr}_{C_f^{1/2}(U)}(M_0^{-1})+t\operatorname{Tr}_{C_f^{1/2}(U)}(M_1^{-1}).\] We note that the weight \(C_f^{1/2}\) on the trace does not affect the claim as it preserves the Loewner order. This yields the claim after applying the argument to the negative trace. ◻
3.2 Differentiability of the expected utility↩︎
Let us now consider the differentiability of the expected utility \(U_B(\mu)\) with respect to the argument \(\mu\). Notice that \(U_B\) is defined on \({\mathcal{P}}^B({\mathcal{X}})\), which differs slightly from much of the existing literature—particularly on gradient flows—where functionals are typically posed on the space of probability measures \({\mathcal{P}}({\mathcal{X}})\). For consistency, we introduce a functional \(V_B : {\mathcal{P}}({\mathcal{X}}) \to \mathbb{R}\) defined by \[V_B(\mu) = U_B(B \mu) = - \operatorname{Tr}_{C^{1/2}(U)}\left(B \cdot S_\mu + I\right)^{-1}.\] Clearly, optimizing \(U_B\) and \(V_B\) are equivalent up to the multiplicative constant \(\mu\mapsto B\mu\).
Recall that given a functional \(F: {\mathcal{P}}({\mathcal{X}}) \to \mathbb{R}\), we denote by \(\frac{\delta F(\mu)}{\delta \mu} : {\mathcal{X}}\to \mathbb{R}\), if it exists, the unique (up to additive constants) function such that \[\label{eq:first95variation} \frac{d}{dt} F(\mu + t \delta \mu)|_{t=0} = \int_{\mathcal{X}}\frac{\delta F(\mu)}{\delta \mu} (x)\delta \mu(dx)\tag{13}\] for every signed measure \(\delta \mu\) such that, for some \(t_0>0\) and \(t\in [0,t_0]\), the measure \(\mu + t \delta \mu \in {\mathcal{P}}({\mathcal{X}})\). The function \(\frac{\delta F}{\delta \mu} (\mu)\) is called first variation of the functional \(F\) at \(\mu\) [21]. In order to guarantee differentiability below, we make the following regularity assumption regarding the kernel \(k\).
Assumption 1. The feature map \(x\mapsto k_x\) belongs to \(C^2(\overline{\mathcal{X}}, H)\) and \[\sup_{x\in{\mathcal{X}}} \left\|k_x\right\|_H + \sup_{x\in{\mathcal{X}}} \left\|\partial_i k_x\right\|_H + \sup_{x\in{\mathcal{X}}} \left\|\partial_i \partial_j k_x\right\|_H < \infty\] for any \(i,j = 1,...,n\).
Proposition 3. Suppose Assumption 1 holds. The first variation of \(V_B\) at \(\mu\) satisfies \[\frac{\delta V_B}{\delta\mu}(\mu)(x) = B \left\|C_f^{\frac{1}{2}} (B S_\mu+I)^{-1} C_f^{\frac{1}{2}}A^* k_x\right\|_U^2\] for all \(x\in {\mathcal{X}}\).
Proof. Let \(\mu\) be fixed and \(\delta\mu\) a signed perturbation with \(t_0>0\) such that \(\mu+t\delta\mu \in {\mathcal{P}}({\mathcal{X}})\) for \(t\in [0,t_0]\). We first observe that \(S_{\delta\mu}\) is well-defined by 11 as a linear, self-adjoint and trace-class operator since due to linearity and since \(\mu, \mu+t\delta\mu \in {\mathcal{P}}({\mathcal{X}})\) one can set \[S_{\delta \mu} = \frac{S_{\mu + t_0\delta \mu} - S_{\mu}}{t_0}.\] Now it holds that \[\begin{align} V_B(\mu+t\delta\mu)-V_B(\mu) &= -\operatorname{Tr}_{C_f^{1/2}(U)}\,\left[(B(S_\mu+t S_{\delta\mu})+I)^{-1} - (BS_\mu + I)^{-1}\right] \\ &= tB\,\operatorname{Tr}_{C_f^{1/2}(U)}\!\left[(BS_\mu+I)^{-1} S_{\delta\mu} (BS_\mu+I)^{-1}\right] + \mathcal{E}(t), \end{align}\] using the standard identity \[(B(S_\mu+t S_{\delta\mu})+I)^{-1} = (BS_\mu+I)^{-1}-t B (BS_\mu+I)^{-1} S_{\delta\mu} (BS_\mu+I)^{-1}+\mathcal{E}(t),\] where \[\mathcal{E}(t) = t^2 B^2 (BS_\mu+I)^{-1} S_{\delta\mu} (B(S_\mu+t S_{\delta\mu})+I)^{-1} S_{\delta\mu} (BS_\mu+I)^{-1}.\] and, consequently, \(\left\|\mathcal{E}(t)\right\| = {\mathcal{O}}(t^2)\). Using the cyclic property of the trace and the identity \(\operatorname{Tr}(a\otimes b) = \langle a, b\rangle\) yields \[\begin{align} \frac{d}{dt} V_B(\mu + t \delta \mu)\big|_{t=0} &= B \operatorname{Tr}_{C_f^{1/2}(U)}\!\left[(BS_\mu+I)^{-1} S_{\delta\mu} (BS_\mu+I)^{-1}\right] \\ &= B \int_{\mathcal{X}} \operatorname{Tr}_{U}\!\left[(BS_\mu+I)^{-1} C_f (BS_\mu+I)^{-1} (C_f^{\frac{1}{2}}A^* k_x) \otimes (C_f^{\frac{1}{2}} A^* k_x) \right] \,\delta\mu(dx) \\ &= B \int_{\mathcal{X}}\langle (B S_\mu+I)^{-1} C_f (B S_\mu+I)^{-1} C_f^{\frac{1}{2}}A^* k_x ,\, C_f^{\frac{1}{2}} A^* k_x\rangle \,\delta\mu(dx) \\ &= B \int_{\mathcal{X}}\left\|C_f^{\frac{1}{2}} (B S_\mu+I)^{-1} C_f^{\frac{1}{2}}A^* k_x\right\|_U^2 \,\delta\mu(dx) \end{align}\] which proves the claim. ◻
Proposition 4. Suppose Assumption 1 holds. Then the first variation \(\frac{\delta V_B}{\delta \mu}\) is continuously differentiable and its gradient is given componentwise by \[\left[\partial_i \frac{\delta V_B}{\delta \mu}(\mu)\right](x) = 2 B \left\langle C_f^{\frac{1}{2}} (B S_\mu+I)^{-1} C_f^{\frac{1}{2}}A^* k_x, C_f^{\frac{1}{2}} (B S_\mu+I)^{-1} C_f^{\frac{1}{2}}A^* (\partial_i k_x)\right\rangle_U.\]
Proof. Since \(C_f^{\frac{1}{2}} (B S_\mu+I)^{-1} C_f^{\frac{1}{2}}A^* : H \to U\) is bounded, our assumption implies that the composition \[x\mapsto C_f^{\frac{1}{2}} (B S_\mu+I)^{-1} C_f^{\frac{1}{2}}A^* k_x : {\mathcal{X}}\to U\] is continuously as a map into \(U\) and \[\partial_i \left[C_f^{\frac{1}{2}} (B S_\mu+I)^{-1} C_f^{\frac{1}{2}}A^* k_x\right] = C_f^{\frac{1}{2}}(B S_\mu+I)^{-1} C_f^{\frac{1}{2}}A^* (\partial_i k_x).\] Now the result is given by the chain rule. ◻
4 Optimization by Wasserstein gradient flow↩︎
4.1 Preliminaries↩︎
We briefly recall the Wasserstein–\(2\) geometry on probability measures and the associated notion of gradient flow. Detailed expositions can be found in [20], [21], [32]. Throughout, let \({\mathcal{X}}\subset \mathbb{R}^d\) be a closed domain endowed with the Euclidean distance. We denote by \(\mathcal{P}_2({\mathcal{X}}) \subset \mathcal{P}({\mathcal{X}})\) the subset of probability measures with finite second moment, i.e., \[\mathcal{P}_2({\mathcal{X}}) = \Big\{ \mu\in\mathcal{P}({\mathcal{X}}): \int_{{\mathcal{X}}} |x|^2\,\mu(dx) < \infty \Big\}.\]
Definition 1 (Wasserstein \(2\)–distance). For \(\mu,\nu\in\mathcal{P}_2({\mathcal{X}})\), the Wasserstein–\(2\) distance is defined as \[W_2(\mu,\nu) = \Bigg( \inf_{\gamma\in\Gamma(\mu,\nu)} \int_{{\mathcal{X}}\times{\mathcal{X}}} |x-y|^2\,\gamma(dx,dy) \Bigg)^{1/2},\] where \(\Gamma(\mu,\nu)\) denotes the set of transport plans between \(\mu\) and \(\nu\), i.e. \[\begin{align} \Gamma(\mu,\nu) = \Big\{ \gamma\in\mathcal{P}({\mathcal{X}}\times{\mathcal{X}}): &\int_{{\mathcal{X}}\times{\mathcal{X}}}\phi(x)\,\gamma(dx,dy) = \int_{{\mathcal{X}}}\phi(x)\,\mu(dx),\\ &\int_{{\mathcal{X}}\times{\mathcal{X}}}\varphi(y)\,\gamma(dx,dy) = \int_{{\mathcal{X}}}\varphi(y)\,\nu(dy),\\ &\text{for all }\phi,\varphi\in C_b({\mathcal{X}}) \Big\}. \end{align}\]
The metric space \((\mathcal{P}_2({\mathcal{X}}),W_2)\) possesses a formal Riemannian structure in which tangent vectors at \(\mu\) are identified with velocity fields \(v\in L^2(\mu;\mathbb{R}^d)\) through the continuity equation \[\partial_t \mu_t + \nabla_x\cdot(\mu_t v_t)=0 .\] Now let \(F : \mathcal{P}_2({\mathcal{X}}) \to \mathbb{R}\) be a functional and recall the definition of the first variation in 13 .
Definition 2. (Wasserstein gradient) Assume that the first variation \(\frac{\delta F(\mu)}{\delta \mu} :{\mathcal{X}}\to \mathbb{R}\) is continuously differentiable. The Wasserstein gradient of \(F\) at \(\mu\) is \[\nabla_{W_2} F(\mu) = \nabla_x\frac{\delta F(\mu)}{\delta \mu}.\]
The associated Wasserstein gradient ascending flow is characterized by the evolution equation \[\label{eq:ascent} \partial_t \mu_t = -\nabla\cdot\left(\mu_t \nabla_{W_2}(F)(\mu_t)\right) = -\nabla_x\cdot\left(\mu_t\nabla_x\frac{\delta F(\mu_t)}{\delta \mu_t}\right).\tag{14}\] This evolution equation can be interpreted as transporting mass along the velocity field \[v_t(x) = \nabla_x\frac{\delta F(\mu_t)}{\delta \mu_t}(x).\] The corresponding characteristic curves, i.e., the trajectories followed by individual particles under this transport, satisfy the ordinary differential equation \[\label{grad95flow95eq} \frac{d x(t)}{dt} = v_t(x(t)) = \nabla_{x}\frac{\delta F(\mu_t)}{\delta \mu_t}(x(t)).\tag{15}\] For numerical computations, the probability measure \(\mu_t\) is unknown and is approximated by an empirical distribution of finitely many interacting particles \(\mu_t^N = \frac{1}{N}\sum_{i=1}^N \delta_{x_i(t)}\), where the particle positions \(x_i(t)\) evolve according to the characteristic dynamics above; see e.g. [3]. In particular, for the optimization problem \[\mu^* = \arg\max_{\mu\in\mathcal{P}_2({\mathcal{X}})} F(\mu),\] the measure \(\mu_t\) is evolved along the above ascent flow, which transports mass in the direction of increasing first variation \(\frac{\delta F}{\delta\mu}\) and converges, under suitable concavity and regularity assumptions, to maximizers of \(F\).
The following theorem establishes that the Wasserstein derivative of \(V_B\) exists.
Theorem 5 (Sufficient conditions for Wasserstein differentiability). Let \({\mathcal{X}}\subset\mathbb{R}^d\) be a bounded, connected domain with \(C^1\) boundary. Assume the setting of Section 2 and, in particular, that \(V_B:\mathcal{P}_2({\mathcal{X}})\to\mathbb{R}\) admits the first variation \(\phi_\mu=\frac{\delta V_B(\mu)}{\delta\mu}\) given in Proposition 3. Moreover, suppose Assumption 1 holds. Let \((\mu_t)_{t\in[0,T]}\) be any \(W_2\)-absolutely continuous curve with velocity field \(v_t\in L^2(\mu_t)\) solving the continuity equation with no-flux boundary condition, \[\partial_t\mu_t+\nabla\cdot(\mu_t v_t)=0 \quad \text{in }{\mathcal{X}}, \qquad (\mu_t v_t)\cdot n=0 \quad \text{on }\partial {\mathcal{X}}.\] Then \(t\mapsto V_B(\mu_t)\) is absolutely continuous and, for a.e.\(t\in[0,T]\), \[\frac{d}{dt}V_B(\mu_t)=\int_{\mathcal{X}}\nabla\phi_{\mu_t}(x)\cdot v_t(x)\,\mu_t(dx).\]
Proof. For convenience, we abbreviate \(g(x):=C_f^{1/2}A^*k_x\in U\) and recall that \(M_g := \sup_{x\in{\mathcal{X}}} \left\|g(x)\right\|_U < \infty\). We write the first variation from Proposition 3 in the form \[\phi_\mu(x)=B\|C_f^{1/2}T_\mu\, g(x)\|_U^2, \qquad T_\mu := (BS_\mu+I)^{-1},\] where \(S_\mu\) is given by 11 and, with notation above, satisfies \(S_\mu = \int_X g(x)\otimes g(x)\,\mu(dx)\).
Since the mapping \(x\mapsto k_x\) is \(C^1\) with bounded derivative and the domain \({\mathcal{X}}\) is bounded, the function \(x\mapsto g(x)\) is Lipschitz continuous; we denote its Lipschitz constant by \(L_g\). Moreover, Proposition 4 gives a formula for \(\nabla\phi_\mu\) in terms of \(g\) and \(\nabla g(x):=C_f^{1/2}A^*(\nabla k_x)\); by Assumption 1 we also have \(\sup_x\|\nabla g(x)\|<\infty\) and \(\nabla g\) is Lipschitz.
Let \(\mu,\nu\in\mathcal{P}_2(X)\) and let \(\pi\in\Pi(\mu,\nu)\) be any coupling. Then \[S_\mu-S_\nu = \int_{{\mathcal{X}}\times{\mathcal{X}}}\big(g(x)\otimes g(x)-g(y)\otimes g(y)\big)\,\pi(dx,\,dy).\] Using \(\|u\otimes u-v\otimes v\|\le (\|u\|+\|v\|)\,\|u-v\|\) for rank-one operators, we obtain \[\|S_\mu-S_\nu\| \le \int_{{\mathcal{X}}\times{\mathcal{X}}} (\|g(x)\|+\|g(y)\|)\,\|g(x)-g(y)\|\,\pi(dx,\,dy) \le 2M_g\,L_g \int_{{\mathcal{X}}\times{\mathcal{X}}} \|x-y\|\,\pi(dx,\,dy).\] By Cauchy–Schwarz, \[\int_{{\mathcal{X}}\times{\mathcal{X}}} \|x-y\|\,\pi(dx,\,dy) \le \left(\int_{{\mathcal{X}}\times{\mathcal{X}}} \|x-y\|^2\,\pi(dx,\,dy)\right)^{1/2}.\] Taking the infimum over couplings yields \[\label{eq:S95lip} \|S_\mu-S_\nu\|\le 2M_gL_g\, W_2(\mu,\nu).\tag{16}\] Since \(S_\mu\ge 0\) we have \(\|T_\mu\|\le 1\) for all \(\mu\). Using the resolvent identity \[T_\mu-T_\nu = (BS_\mu+I)^{-1}B(S_\nu-S_\mu)(BS_\nu+I)^{-1},\] we obtain \[\|T_\mu-T_\nu\|\le \|T_\mu\|\,B\|S_\mu-S_\nu\|\,\|T_\nu\|\le B\|S_\mu-S_\nu\| \le 2BM_g L_g\, W_2(\mu,\nu).\] Hence \(\mu\mapsto T_\mu\) is \(W_2\)-Lipschitz in operator norm.
From Proposition 4, \(\nabla\phi_\mu(x)\) is given by \[\nabla \phi_{\mu}(x) = 2B\langle C_f^{1/2}T_\mu g(x), C_f^{1/2}T_\mu \nabla g(x) \rangle_U\] and hence depends on \((T_\mu, g(x), \nabla g(x))\) polynomially of degree at most two. Since \(C_f^{1/2}\) is a bounded linear operator, there exists \(M_{C_f} = \lVert C_f^{1/2}\rVert < \infty\). Using the bounds \(\|T_\mu\|\le 1\), \(\sup_x\|g(x)\|\le M_g\), \(\sup_x\|\nabla g(x)\|\le M_{\nabla g}\) and the Lipschitz estimates for \(T_\mu\) and for \(g,\nabla g\), one obtains the following: there exists \(C<\infty\) such that for all \(\mu,\nu\in\mathcal{P}_2({\mathcal{X}})\) and all \(x,y\in {\mathcal{X}}\), \[\label{eq:grad95phi95pointwise} \|\nabla\phi_\mu(x)-\nabla\phi_\nu(y)\| \le C\big(\|x-y\| + W_2(\mu,\nu)\big).\tag{17}\] Let \(\pi\in\Pi_{\mathrm{opt}}(\mu,\nu)\) be any optimal coupling, then \(\int \|x-y\|^2\,d\pi = W_2(\mu,\nu)^2\), hence \[\label{eq:plan95lip} \int_{X\times X}\bigl\|\nabla \phi_\mu(x)-\nabla \phi_\nu(y)\bigr\|^2\,\pi(dx,dy) \le 4C^2\,W_2(\mu,\nu)^2.\tag{18}\] Therefore the hypotheses of [20] apply to \(V_B\). Let \((\mu_t)_{t\in[0,T]}\) be a \(W_2\)–absolutely continuous curve with velocity field \(v_t\in L^2(\mu_t)\) solving the continuity equation \(\partial_t\mu_t + \nabla\cdot(\mu_t v_t)=0\) in the sense of distributions (with the no-flux boundary condition imposed). Then \(t\mapsto V_B(\mu_t)\) is absolutely continuous and for a.e.\(t\in[0,T]\), \[\frac{d}{dt}V_B(\mu_t)=\int_X \nabla \phi_{\mu_t}(x)\cdot v_t(x)\,\mu_t(dx).\] ◻
4.2 First Order Optimality Condition↩︎
Having established conditions for Wasserstein differentiability, we now derive first order optimality conditions for (local) maximisers of \(V_B\) in \(\mathcal{P}_2\), following a similar approach to [3] and [33]. The following lemma establishes necessary and sufficient conditions for stationarity.
Lemma 1. Let \({\mathcal{X}}\subset\mathbb{R}^n\) be a bounded, connected domain with \(C^1\) boundary. Assume the hypotheses of Theorem 5. Let \((\mu_t)_{t\in[0,T]}\) be a (distributional) solution of the Wasserstein gradient ascent flow 14 . Then, the map \(t\mapsto V_B(\mu_t)\) is non-decreasing and, for a.e.\(t\in[0,T]\). A measure \(\mu^*\in\mathcal{P}_2({\mathcal{X}})\) is a stationary solution of 14 if and only if \[\label{eq:stationary95weak} \int_X \nabla\psi(x)\cdot \nabla\phi_{\mu^*}(x)\,\mu^*(dx)=0 \qquad \forall\,\psi\in C^1(\overline{\mathcal{X}}).\tag{19}\] In particular, if \(\mu^*(dx)=\rho^*(x)\,dx\) with \(\rho^*\in L^1({\mathcal{X}})\) and \(\rho^*(x)\ge c>0\) a.e.on \({\mathcal{X}}\), then stationarity implies \(\nabla\phi_{\mu^*}(x)=0\) for a.e.\(x\in {\mathcal{X}}\) (and hence also \(\mu^*\)-a.e.).
Proof. Write \(\phi_t:=\phi_{\mu_t}\). Formally differentiating \(V_B(\mu_t)\) along the flow and using 14 gives \[\frac{d}{dt}V_B(\mu_t) = -\int_{\mathcal{X}}\phi_t(x)\,\nabla\cdot\bigl(\mu_t \nabla \phi_t\bigr)(dx).\] Integrating by parts and using the boundary condition we obtain \[\frac{d}{dt}V_B(\mu_t) = \int_{\mathcal{X}}\nabla \phi_t(x)\cdot \nabla \phi_t(x)\,\mu_t(dx) = \int_{\mathcal{X}}\|\nabla \phi_t(x)\|^2\,\mu_t(dx)\ge 0,\] which establishes the monotonicity of \(t\mapsto V_B(\mu_t)\).
Stationarity of 14 means \(\partial_t\mu_t=0\), i.e. \(\nabla\cdot(\mu^*\nabla\phi_{\mu^*})=0\) with no-flux. In distributional form this is equivalent to 19 . For the final claim, assume \(\mu^*(dx)=\rho^*(x)\,dx\) with \(\rho^*\ge c>0\) a.e.and \(\phi_{\mu^*}\in H^1({\mathcal{X}})\). Then 19 implies that \[\int_{\mathcal{X}}\rho^*(x)\,\nabla\phi_{\mu^*}(x)\cdot \nabla\psi(x)\,dx=0 \qquad \forall \psi\in H^1({\mathcal{X}}).\] Taking \(\psi=\phi_{\mu^*}\) yields \[\int_{\mathcal{X}}\rho^*(x)\,\|\nabla\phi_{\mu^*}(x)\|^2\,dx = 0,\] and since \(\rho^*(x)\ge c>0\) a.e.it follows that \(\nabla\phi_{\mu^*}(x)=0\) for a.e.\(x\in {\mathcal{X}}\). ◻
Remark 6. Note that stationary points within the Wasserstein gradient flow do not necessarily correspond to maximisers of \(V_B\). While the functional \(V_B\) is concave in the space of signed measures, it is not geodesically / displacement-convex on \(\mathcal{P}_2({\mathcal{X}})\), and so we cannot rely on the typical arguments that are common for gradient flows of geodesically convex energies.
We now state the first order optimality condition for local maximisers or \(V_B\) in \(\mathcal{P}_2({\mathcal{X}})\).
Proposition 7. Assume the setting and hypotheses of Propositions 3.1–3.3. Then \(\mu^*\in\mathcal{P}({\mathcal{X}})\) is a maximiser of \(V_B\) if and only if there exists a constant \(c\in\mathbb{R}\) such that the shifted first variation \[g(x):=\phi_{\mu^*}(x)-c\] satisfies \(g(x)\le 0\) for all \(x\in {\mathcal{X}}\) and \(g(x)=0\) \(\mu^*\)-a.e.
Proof. Let \(\mu^*\) be a maximiser of \(V_B\) and fix any \(\nu\in\mathcal{P}({\mathcal{X}})\). Define \(\mu_t=(1-t)\mu^*+t\nu\). By concavity of \(V_B\), the map \(t\mapsto V_B(\mu_t)\) is concave on \([0,1]\), hence its right
derivative at \(t=0\) is non-positive: \[0\ge \frac{d}{dt}V_B(\mu_t)\Big|_{t=0^+}.\] By the first variation formula this derivative equals \[\frac{d}{dt}V_B(\mu_t)\Big|_{t=0^+}
=\int_{\mathcal{X}}\phi_{\mu^*}(x)\,(\nu-\mu^*)(dx)
=\int_{\mathcal{X}}\phi_{\mu^*}(x)\,\nu(dx)-\int_{\mathcal{X}}\phi_{\mu^*}(x)\,\mu^*(dx).\] Let \(c:=\int \phi_{\mu^*}d\mu^*\). Choosing \(\nu=\delta_x\) shows \(\phi_{\mu^*}(x)\le \int \phi_{\mu^*}\,d\mu^*\) for all \(x\), hence in particular \(\phi_{\mu^*}(x)\le c\) for all \(x\in{\mathcal{X}}\). Set \(g(x)=\phi_{\mu^*}(x)-c\le 0\) on \({\mathcal{X}}\). Since \(\mu^*\) is supported on \({\mathcal{X}}\), \[\int_{\mathcal{X}}g(x)\,\mu^*(dx)=\int_{\mathcal{X}}\phi_{\mu^*}(x)\,\mu^*(dx)-c \le 0.\] On the other hand, by definition of \(c\) and \(g\le 0\), we also have \[0\le \int_{\mathcal{X}}(c-\phi_{\mu^*}(x))\,\mu^*(dx)= -\int_{\mathcal{X}}g(x)\,\mu^*(dx),\] so \(\int_{\mathcal{X}}g\,d\mu^*=0\).
Since \(g\le 0\) everywhere, this forces \(g=0\) \(\mu^*\)- a.e.
Conversely, suppose there exists \(c\) such that \(g=\phi_{\mu^*}-c\le 0\) on \({\mathcal{X}}\) and \(g=0\) on \(\mathrm{supp}(\mu^*)\). Then for any \(\nu\in\mathcal{P}({\mathcal{X}})\) \[\int_{\mathcal{X}}\phi_{\mu^*}\,(\nu-\mu^*)(dx)
=\int_{\mathcal{X}}g(x)\,\nu(dx)-\int_{\mathcal{X}}g(x)\,\mu^*(dx)
\le 0,\] since \(g\le 0\) and \(g=0\) on \(\mathrm{supp}(\mu^*)\). By concavity of \(V_B\), \[V_B(\nu)\le V_B(\mu^*) + \frac{d}{dt}V_B((1-t)\mu^*+t\nu)\Big|_{t=0^+}
= V_B(\mu^*) + \int_{\mathcal{X}}\phi_{\mu^*}\,(\nu-\mu^*)(dx)
\le V_B(\mu^*).\] Thus \(\mu^*\) is a maximiser of \(V_B\). ◻
Generally, there are strong differences between stationarity and global optimality. However, we do know that the global maximiser of \(U\) is itself a stationary point of \(U\) in \(\mathcal{P}_2\).
Proposition 8. Let \({\mathcal{X}}\subset\mathbb{R}^n\) be a bounded, connected domain with \(C^1\) boundary and assume the hypotheses of Theorem 5. Let \(\mu^*\in\mathcal{P}_2({\mathcal{X}})\) be a global maximiser* of \(V_B\) over \(\mathcal{P}_2({\mathcal{X}})\). Then \(\mu^*\) is a stationary solution of 14 in the distributional sense.*
Proof. Fix \(\psi\in C^1(\overline{\mathcal{X}})\) with \(\partial_n\psi=0\) on \(\partial {\mathcal{X}}\) and consider the (smooth) velocity field \[v := \nabla\psi, \qquad\text{so that}\qquad v\cdot n = \partial_n\psi = 0\;\text{ on }\partial {\mathcal{X}}.\] Let \((\mu_t)_{t\in[0,T]}\) be any \(W_2\)-absolutely continuous curve solving the continuity equation \[\partial_t\mu_t+\nabla\cdot(\mu_t v)=0 \quad\text{in }{\mathcal{X}}, \qquad (\mu_t v)\cdot n=0 \quad\text{on }\partial {\mathcal{X}},\] with initial condition \(\mu_0=\mu^*\); such a curve can be realised by the flow map of \(v\) for small times. By global optimality of \(\mu^*\) we have \(V_B(\mu_t)\le V_B(\mu^*)\) for all sufficiently small \(t\ge 0\), hence the right derivative at \(t=0\) satisfies \[\frac{d}{dt}V_B(\mu_t)\Big|_{t=0^+}\le 0.\] Applying the chain rule from Theorem 5(ii) at \(t=0\) gives \[\frac{d}{dt}V_B(\mu_t)\Big|_{t=0^+} = \int_{\mathcal{X}}\nabla\phi_{\mu^*}(x)\cdot v(x)\,\mu^*(dx) = \int_{\mathcal{X}}\nabla\phi_{\mu^*}(x)\cdot \nabla\psi(x)\,\mu^*(dx),\] and therefore \[\label{eq:ineq95plus} \int_{\mathcal{X}}\nabla\phi_{\mu^*}(x)\cdot \nabla\psi(x)\,\mu^*(dx)\le 0 \qquad \forall\,\psi\in C^1(\overline{\mathcal{X}})\text{ with }\partial_n\psi=0.\tag{20}\] Applying the same argument to \(-\psi\) yields the reverse inequality, hence equality in 20 : \[\int_{\mathcal{X}}\nabla\phi_{\mu^*}(x)\cdot \nabla\psi(x)\,\mu^*(dx)=0 \qquad \forall\,\psi\in C^1(\overline{\mathcal{X}})\text{ with }\partial_n\psi=0.\] This is the distributional stationarity condition \(\nabla\cdot(\mu^*\nabla\phi_{\mu^*})=0\) in \({\mathcal{X}}\) with no-flux boundary condition, as required. ◻
Proposition 9. All local maximisers of \(V_B\) are global maximisers in \(\mathcal{P}_2({\mathcal{X}})\).
Proof. Suppose that \(\mu_0\) is a local maximum of \(V_B\), so that, for some \(\delta > 0\), for all \(\mu \in
\mathcal{P}(\mathcal{X})\) such that \(W_2(\mu, \mu_0) \leq \delta\) we have \(V_B(\mu_0) > V_B(\mu)\). Suppose also that \(V_B(\mu_0) <
V_B(\mu^*)\).
Consider the mixture \(\mu_t = (1-t)\mu_0 + t \mu^*\), for \(t\in [0,1]\), then by concavity \[V_B(\mu_t) \geq (1-t)V_B(\mu_0) + t V_B(\mu^*) >
V_B(\mu_0).\] Noting that \[W_2(\mu_t, \mu_0) \leq \sqrt{t} W_2(\mu^*, \mu_0),\] then choosing \(t\) sufficiently small yields a contradiction. It follows that \(V_B(\mu_0) = V_B(\mu^*)\). ◻
While concavity implies that all local maximisers of \(V_B\) are themselves global optimisers, this is not initself sufficient to establish a unique global maximiser. To see this, we consider the following toy example. For simplicity we work over the unit torus, to ignore effects of the boundary, however we note that similar counterexamples can be established in the more general setting.
Example 2. Consider \({\mathcal{X}}:= \mathbb{T} = \mathbb{R} / \mathbb{Z}\) and let \(k(x,y) := \cos\big(2\pi(x-y)\big)\) with \(x,y\in\mathbb{T}\). This kernel is positive semi-definite, translation invariant and admits the finite–dimensional feature representation \[k(x,y) = \phi(x)^\top \phi(y), \qquad \phi(x) := \begin{pmatrix} \cos(2\pi x)\\ \sin(2\pi x) \end{pmatrix} \in \mathbb{R}^2 .\] The associated RKHS \(H\) is two–dimensional and consists of functions of the form \[f(x)=\theta^\top \phi(x), \qquad \theta\in\mathbb{R}^2 .\] We consider the inverse problem, where \(A(x)f = f(x)\), i.e. \[y = A(x)f + \varepsilon = f(x) + \varepsilon = \theta^\top \phi(x) + \varepsilon, \qquad \varepsilon\sim\mathcal{N}(0,1),\] with Gaussian prior \[\theta \sim \mathcal{N}(0,\sigma^2 I_2).\] For a given probability measure \(\mu\in {\mathcal{P}}(\mathbb{T})\), the covariance operator \(B_\mu\) has the matrix representation \[\label{eq:example95Bmu} B_\mu = \int_{\mathbb{T}} \phi(x)\phi(x)^\top \,\mu(dx) = \frac{1}{2} \int_{\mathbb{T}} \begin{pmatrix} 1 + \cos(4\pi x) & \sin (4\pi x) \\ \sin(4 \pi x) & 1 - \cos(4 \pi x) \end{pmatrix} \,\mu(dx) \in \mathbb{R}^{2\times 2}.\tag{21}\] Since \(\|\phi(x)\|^2=1\) for all \(x\), every feasible \(B_\mu\) satisfies \(Tr(B_\mu)=1\), and with the posterior covariance of \(\theta\) given by \(C_{\mathrm{post}}(\mu) = \big(\sigma^{-2} I_2 + B_\mu\big)^{-1}\), the expected utility with \(B=1\) is \[U_1(\mu) = -\operatorname{Tr}\big(C_{\mathrm{post}}(\mu)\big) = -\operatorname{Tr}\big(\sigma^{-2} I_2 + B_\mu\big)^{-1}.\] Let \(\lambda_1,\lambda_2\) be the eigenvalues of \(B_\mu\). Since \(\lambda_1+\lambda_2=1\), convexity of \(t\mapsto (\sigma^{-2}+t)^{-1}\) implies \[U(\mu) \le -2(\sigma^{-2}+1/2)^{-1},\] with equality if and only if \[\label{eq:example95condition} B_\mu = \tfrac12 I_2.\tag{22}\] We note that with repeated eigenvalue, the symmetry of \(B_\mu\) implies the latter identity. Thus the covariance operator is unique, but observe that the maximising measure need not be. In fact, there exist infinitely many maximising measures \(\mu\), since the condition 22 is equivalent to requiring that \[\int_{\mathbb{T}} \cos(4 \pi x) \mu(dx) = 0\quad \text{and} \quad \int_{\mathbb{T}} \sin(4 \pi x) \mu(dx) = 0.\] Consequently, characterizing when the maximising measure is unique depends intrinsically on the observational model and lies beyond the scope of this work.
5 Particle Approximation and Regularisation↩︎
In previous sections, we have established a Wasserstein gradient flow approach to optimize \(U_B\) or, more precisely, the normalized utility \(V_B\). The method for optimizing \(U_B\) is summarized in Algorithm 1, where we utilize the forward Euler method similar to [3].
While Algorithm 1 has theoretical guarantees of convergence, there can be practical constraints in the sensor placement task that need to be considered. For a batch size \(B>1\), such constraints can involve requiring exactly \(B\) locations in the output. In the vanilla implementation of the method, this is currently not guaranteed and the support of the maximizing measure can vary. Secondly, the sensors may have a physical constraints related to their relative locations, e.g., we cannot place two sensors too close to each other. Such considerations require us to separate the mass associated to each design variable in \(\mu \in {\mathcal{P}}^B({\mathcal{X}})\) and enforce preferable solutions through regularization.
In order to enforce these ideas, we lift the problem back to the tensor space \({\mathcal{X}}^B = {\mathcal{X}}\times \ldots \times {\mathcal{X}}\) and consider optimization of a product measure \(\boldsymbol{\mu}= \otimes_{j=1}^B \mu_j\) over \(\otimes_{j=1}^B {\mathcal{P}}({\mathcal{X}})\). We formulate a utility function \({\boldsymbol{U}}_B\) for \(\boldsymbol{\mu}\) by setting \[\label{eq:UBbmu} {\boldsymbol{U}}_B(\boldsymbol{\mu}) := U_B\left(\sum_{j=1}^B \mu_j\right) = V_B\left(\frac{1}{B}\sum_{j=1}^B \mu_j\right).\tag{23}\] We immediately observe that \({\boldsymbol{U}}_B\) preserves the concavity and regularity properties of \(V_B\). In the same vein, we formally have \[\begin{align} \frac{d}{dt} {\boldsymbol{U}}_B(\boldsymbol{\mu}+ t \delta \boldsymbol{\mu}) \big|_{t=0} & = & \frac{d}{dt} V_B\left(\sum_{j=1}^B \mu_j + t \sum_{j=1}^B \delta \mu_j\right) \\ & = & \sum_{j=1}^B \int_{\mathcal{X}}B \cdot \left\|C_f^{\frac{1}{2}} (B S_\mu + I)^{-1} C_f^{\frac{1}{2}} A^* k_x\right\|_U^2 \delta \mu_j(dx) \\ & = & \int_{{\mathcal{X}}^B} B\sum_{j=1}^B \left\|C_f^{\frac{1}{2}} (B S_\mu + I)^{-1} C_f^{\frac{1}{2}} A^* k_{x_j}\right\|_U^2 \delta \boldsymbol{\mu}(d\boldsymbol{x}), \end{align}\] where we utilize notation \(\boldsymbol{x}= (x_1, ..., x_B)\). It follows that one has \[\frac{\delta {\boldsymbol{U}}_B(\boldsymbol{\mu})}{\delta \boldsymbol{\mu}}(\boldsymbol{x}) = B\sum_{j=1}^B \left\|C_f^{\frac{1}{2}} (B S_\mu + I)^{-1} C_f^{\frac{1}{2}} A^* k_{x_j}\right\|_U^2\] and, consequently, \[\partial_{j,i} \frac{\delta {\boldsymbol{U}}_B(\boldsymbol{\mu})}{\delta \boldsymbol{\mu}}(\boldsymbol{x}) = \left[\partial_i \frac{\delta V_B}{\delta \mu}(\mu)\right](x_j)\] where \(\partial_{j,i}\) stands for the derivative applied to the component \(i\) of vector \(x_j\).
Next, we introduce two regularization functionals acting on such product measures and derive their role in Wasserstein gradient flows. The first functional promotes concentration within each batch, while the second introduces repulsion between different batches in state space. The first regularization functional we utilize is the variance \[\mathcal{R}_v(\boldsymbol{\mu}) = \sum_{j=1}^B \mathbb{E}^{\mu_j} |x-\mathbb{E}^{\mu_j} x|^2,\] which naturally promotes the concentration of each individual ensemble in the gradient flow. The second regularizer is a repulsive interaction energy based on kernel embeddings, which encourages different measures to separate in state space. To define it, let \(q:{\mathcal{X}}\times{\mathcal{X}}\to\mathbb{R}\) be a symmetric, continuously differentiable positive definite kernel. Now we set \[\mathcal{R}_r(\boldsymbol{\mu}) = \sum_{j\neq j'} \iint_{{\mathcal{X}}\times {\mathcal{X}}} q(x,x')\, \mu_j(dx)\, \mu_{j'}(dx').\] Up to additive self–interaction terms \(\int q(x,x')\,\mu_j(dx)\mu_j(dx')\), this functional coincides with the negative squared maximum mean discrepancy between pairs of measures in the reproducing kernel Hilbert space associated with \(q\). As such, minimizing \(-\mathcal{R}_r\) encourages separation of the distributions \(\mu_j\) in feature space and induces mutual repulsion in the corresponding Wasserstein dynamics.
Adding the regularization terms to the utility function, we define \[H_{\alpha,\beta}(\boldsymbol{\mu}) = U_B(\boldsymbol{\mu}) - \alpha \mathcal{R}_r(\boldsymbol{\mu}) - \beta \mathcal{R}_v(\boldsymbol{\mu}),\] where \(\alpha,\beta \geq 0\) stand for the regularization parameters, and study the gradient flow \[\partial_t \boldsymbol{\mu}= -\nabla_{W_2}(\boldsymbol{\mu}) = \nabla_{x} \left( \boldsymbol{\mu}\nabla_{x} \frac{\delta H_{\alpha, \beta}(\boldsymbol{\mu})}{\delta \boldsymbol{\mu}}\right).\]
Let us next establish the first variation and gradients corresponding to the regularization terms.
Proposition 10. For \(\mu\in\mathcal{P}_2({\mathcal{X}})\), we have \[\frac{\delta \mathcal{R}_v(\boldsymbol{\mu})}{\delta\boldsymbol{\mu}}(\boldsymbol{x}) = \sum_{j=1}^B |x_j-\mathbb{E}^{\mu_j} x|^2 \quad \text{and} \quad \nabla_{x_i} \frac{\delta\mathcal{R}_v(\boldsymbol{\mu})}{\delta\boldsymbol{\mu}}(\boldsymbol{x}) = 2 (x_i-\mathbb{E}^{\mu_i} x).\]
Proof. Let \(\delta\boldsymbol{\mu}= \otimes_{j=1}^B \delta \mu_j\) be a signed measure such that there exists \(t_0>0\) such that \(\boldsymbol{\mu}+ t\delta \boldsymbol{\mu}\in \otimes_{j=1}^B\mathcal{P}_2({\mathcal{X}})\). Write \(m(\mu) = \int_{\mathcal{X}}x \mu(dx)\) for shorthand. We have \[\begin{align} \mathcal{R}_v(\boldsymbol{\mu}+t\delta \boldsymbol{\mu}) & = & \sum_{j=1}^B \int_{\mathcal{X}}| x - m(\mu_j) - t m(\delta \mu_j)|^2 (\mu_j(dx) + t \delta \mu_j(dx)) \\ & = & \mathcal{R}_v(\boldsymbol{\mu}) + \sum_{j=1}^B \left[t^2 |m(\delta \mu_j)|^2 + t \int_{\mathcal{X}}|x-m(\mu_j)|^2 \delta \mu_j(dx) + t^3 |m(\delta \mu_j)|^2\right] \end{align}\] and, therefore, \[\frac{d}{dt} \mathcal{R}_v(\boldsymbol{\mu}+ t \delta \boldsymbol{\mu})\Big|_{t=0} = \sum_{j=1}^B \left[\int_{\mathcal{X}}|x-m(\mu_j)|^2\,\mu_j(dx)\right] = \int_{{\mathcal{X}}^B} \left[\sum_{j=1}^B|x_j-m(\mu_j)|^2\right] \boldsymbol{\mu}(d\boldsymbol{x}).\] Now the formula for the gradient follows directly. ◻
Proposition 11. For \(\mu\in\mathcal{P}_2({\mathcal{X}})\), we have \[\frac{\delta {\mathcal{R}}_r(\boldsymbol{\mu})}{\delta \boldsymbol{\mu}}(\boldsymbol{x}) = 2 \sum_{j\neq j'} \iint_{{\mathcal{X}}\times {\mathcal{X}}} q(x_j,x)\, \mu_{j'}(dx).\] and \[\nabla_{x_i} \frac{\delta {\mathcal{R}}_r(\boldsymbol{\mu})}{\delta \boldsymbol{\mu}}(\boldsymbol{x}) = 2 \sum_{\substack{j=1 \\ j\neq i}}^B\int_{\mathcal{X}}\nabla_{x_i}q(x_i,x) \mu_{j}(dx).\]
Proof. Let \(\delta\boldsymbol{\mu}= \otimes_{j=1}^B \delta \mu_j\) be a signed measure such that there exists \(t_0>0\) such that \(\boldsymbol{\mu}+ t\delta \boldsymbol{\mu}\in \otimes_{j=1}^B\mathcal{P}_2({\mathcal{X}})\). We first observe that \[\begin{gather} {\mathcal{R}}_r(\boldsymbol{\mu}+ t \delta \boldsymbol{\mu}) - {\mathcal{R}}_r(\boldsymbol{\mu}) \\ = \sum_{j\neq j'} \iint_{{\mathcal{X}}\times {\mathcal{X}}} q(x,x') \left(t \mu_j(dx) \delta \mu_{j'}(dx') + t \delta \mu_j(dx) \mu_{j'}(dx) + t^2 \delta \mu_j(dx) \delta \mu_{j'}(dx') \right) \end{gather}\] and, consequently, \[\begin{align} \frac{d}{dt} {\mathcal{R}}_r(\boldsymbol{\mu}+ t \delta \boldsymbol{\mu})\big|_{t=0} & = & 2 \sum_{j\neq j'} \iint_{{\mathcal{X}}\times {\mathcal{X}}} q(x,x') \delta \mu_j(dx) \mu_j(dx') \\ & = & 2 \sum_{j=1}^B \int_{\mathcal{X}}\left[\sum_{\substack{j'=1 \\ j'\neq j}}^B\int_{\mathcal{X}}q(x,x') \mu_{j'}(dx') \right] \delta \mu_j(dx) \\ & = & \int_{\mathcal{X}}\left[2 \sum_{j=1}^B\sum_{\substack{j'=1 \\ j'\neq j}}^B\int_{\mathcal{X}}q(x_j,x') \mu_{j'}(dx') \right] \delta \boldsymbol{\mu}(d\boldsymbol{x}) \\ & = & \int_{\mathcal{X}}\left[2 \sum_{j\neq j'}\int_{\mathcal{X}}q(x_j,x') \mu_{j'}(dx') \right] \delta \boldsymbol{\mu}(d\boldsymbol{x}). \end{align}\] This yields the claim. ◻
Algorithm 2 summarizes the regularized Wasserstein gradient flow.
6 Numerical Simulations↩︎
In this section, we demonstrate the use of the algorithms proposed in previous section with numerical examples in the setting of inverse source term problem for the Poisson’s equation and the time-harmonic Schrödinger equation. Furthermore, we test the sensitivity of the \(\mathcal{R}_v\) and \(\mathcal{R}_r\) regularizers to show how the selection of the regularization parameters affects.
6.1 Inverse source problem for one-dimensional Poisson’s equation↩︎
Consider the one-dimensional Poisson’s equation with homogeneous Dirichlet boundary conditions \[\label{eq:poisson951D} \left\{ \begin{align} -\dfrac{\partial^2 u}{\partial x^2}(x) &= f(x), && x \in {\mathcal{X}}=[0,1],\\ u(x) &= 0, && x \in \partial{\mathcal{X}}. \end{align} \right.\tag{24}\] Here, the inverse problem is to identify the source term \(f\) given noisy point observations of the solution \(u\) in \({\mathcal{X}}\). For a source \(f\in H^1(0,1)\), it is well-known that the problem has a unique solution \(u\in H\), where \[H := H^3(0,1)\cap H_0^1(0,1), \qquad \langle u,v\rangle_H := \int_{\mathcal{X}}\frac{\partial^3 u}{\partial x^3}(x)\,\frac{\partial^3 v}{\partial x^3}(x)\,dx .\] By Sobolev embedding theorem, the embedding \(H\hookrightarrow C^2({\mathcal{X}})\) is continuous. Consequently, point evaluation and its first two derivatives are continuous linear functionals on \(H\), so \(H\) is a reproducing kernel Hilbert space. In particular, the mapping \(x\mapsto k_x\in H\) is twice continuously Fréchet differentiable. We note that \(C^2({\mathcal{X}})\) stands here for functions on \(C^2(0,1)\) that are continuously extended to the closed interval.
Our objective is to identify optimal configuration of \(B\) sensor locations for measuring \(u\) simultaneously. In the context of previous analysis, the linear operator \(A\) describes the mapping \(f \mapsto u : H^1(0,1) \to H\). Each independent point observation is contaminated with additive zero-mean Gaussian noise with variance 0.01. For the numerical implementation, the operator \(A\) is discretized on a grid of \(N=100\) equidistant points. Since the Green’s function for 24 is classical, the output function \(f\) can be evaluated continuously on \({\mathcal{X}}\). In this example, we impose a non-stationary zero-mean Gaussian prior distribution for \(f\) based on covariance function \[c(x,z) = \left(1 + 50\left(x-\frac{1}{2}\right)^2\right)\left(1 + 50\left(z-\frac{1}{2}\right)^2\right) \exp\left(-\frac{|x-z|^2}{2\sigma_0^2}\right),\] where \(\sigma_0 = 0.01\). We simulate Algorithm 1 using 120 particles initialized from a uniform distribution over \({\mathcal{X}}\). For simplicity, we take a batch size \(B=2\) and regularize the empirical variance with weight \(\alpha=0.008\). Here, no repulsive regularizer was applied, i.e. \(\beta=0\). Using 500 time steps of size \(4e-3\), we illustrate the results in Figures 3 and 4.
First, we display the expected utility for the two design variables together with the two true optimal configurations marked with red dots \((x_1,x_2)\approx(0.1616,0.8384)\) and \((x_1,x_2)\approx(0.8384, 0.1616)\) (note the symmetry). The optimal configuration found via Algorithm 1 is marked with the cross. Second, Figure 4 shows the evolution of the particle trajectories as well as the resulting empirical distribution of the measure \(\mu\). We observe that the empirical distribution converges to the optimal configuration.
6.2 Time-harmonic Schrödinger equation↩︎
Consider the two-dimensional boundary value problem \[\label{eq:timeharmonicschrodinger} \begin{cases} (-\Delta + V - \omega^2)u = f & \text{in } \Omega=(0,1)^2,\\ u = 0 & \text{on } \partial\Omega, \end{cases}\tag{25}\] where \(V\in L^\infty(\Omega)\) is real-valued and satisfies \(V\ge 0\) a.e. in \(\Omega\). By the Rayleigh–Ritz (min–max) principle, \(\lambda_1(-\Delta+V) \geq \lambda_1(-\Delta)=2\pi^2\). Hence, if \(\omega^2 < 2\pi^2\), then \(\omega^2\) cannot be a Dirichlet eigenvalue of \(-\Delta+V\), and 25 admits a unique solution. The inverse problem is, given \(V\) and noisy point observations of \(u\), to identify \(f\).
The mapping properties of 25 are classical. Following the template of Section 6.1, our aim is to ensure that the solution belongs to a reproducing kernel Hilbert space that is continuously embedded into \(C^1({\mathcal{X}})\). To this end, we recall from [34] that for any \(V \in C^\infty({\mathcal{X}})\) and \(f\in H^s({\mathcal{X}})\), \(s\geq 0\), the unique solution satisfies \(u\in H^{s+2}({\mathcal{X}}) \cap H^1_0({\mathcal{X}})\). By the Sobolev embedding theorem, any \(s>1\) guarantees that \(H = H^{s+2}({\mathcal{X}}) \cap H^1_0({\mathcal{X}}) \hookrightarrow C^2({\mathcal{X}})\) continuously.
We define the smooth potential term by setting \[\begin{align} V = h_\epsilon \ast W|_{{\mathcal{X}}}, \quad \text{where} \quad W(z,y) = 200\cdot \mathbf{1}\left(\left|z-\frac{1}{2}\right|\leq 0.08\right) \cdot \mathbf{1}\left(\left|y- \frac{1}{2}\right| \leq 0.08\right) \in L^2(\mathbb{R}^2), \end{align}\] where \(h_\epsilon\) is Gaussian kernel with a small variance \(\epsilon\) convolving \(W\), which is an indicator function of a cross-shaped field.
Similarly to the Poisson example, our objective is to identify the optimal configuration of sensor locations, here, for the batch size \(B=4\). The independent point observations are contaminated with zero-mean Gaussian noise with variance 0.01. In the numerical implementation of the forward mapping \(A\), both the input and output domains are discretized on a regular grid with \(N=60^2\) points. For continuous evaluation of the output function, the values are interpolated based on the grid values. We impose a zero-mean squared-exponential Gaussian prior with a covariance function \[c(x,z) = \exp\left(-\frac{|x-z|^2}{2\sigma_0^2}\right),\] where \(\sigma_0 = 0.5\). We simulate Algorithm 1 using 100 particles initialized from a uniform distribution over \({\mathcal{X}}\). We regularize the empirical variance with weight \(\alpha=5e-4\) and apply no repulsive regularizer, i.e. \(\beta=0\). Using 120 time steps of size \(8e-3\), we illustrate the results in Figure 5. First, the expected utility for a single observation is mapped on the left-hand side of the image with the reconstructed optimal sensor locations. The right-hand side plot visualizes the evolution of the ensembles towards these optimal locations.
6.3 Sensitivity of regularization↩︎
The regularization of the particle gradient flow relies on two distinct terms: the regularizer \(\mathcal{R}_v\), which attracts particles within the same ensemble, and the regularizer \(\mathcal{R}_r\), which enforces separation between two different ensembles. The best choice of the corresponding regularization parameters is influeced by various factors including the underlying model, the discretization, the number of design variables, and the time step \(dt\), and, below, our aim is to illustrate these effects numerically.
6.3.1 Regularization with \(\mathcal{R}_v\)↩︎
Consider the Poisson source identification problem in Section 6.1 for batch-size \(B=3\). We use the same parameters as in Section 6.1 with the exception that we use 300 time steps of size \(1e-5\) and initialize 900 particles from the uniform distribution such that \(\mu_N^1\sim \mathcal{U}(0,\frac{1}{3})\), \(\mu_N^2\sim \mathcal{U}(\frac{1}{3},\frac{2}{3})\) and \(\mu_N^3\sim \mathcal{U}(\frac{2}{3},1)\). Figure 6 (left) displays the magnitude of the regularization for several parameter choices. For sufficiently large values, the measure converges to the desired number of locations before the regularization becomes dominant. In the empirical distribution plot, we observe that for smaller values of parameter \(\alpha\) the particle gradient flow has not yet converged to the desired number of optimal measurement locations. To isolate the effect of \(\alpha\), the parameters for the time step and its size are chosen so that the particle flow has not yet fully converged to the optimum. This is done since for this problem, when \(B>2\) the gradient flow would converge to the optimal measurement locations of which some would overlap, without the regularizer \(\mathcal{R}_r\), thus making the effect of only \(\alpha\) not apparent.
6.3.2 Regularization with \(\mathcal{R}_v\) and \(\mathcal{R}_r\)↩︎
As the number of design parameters increases, the ensembles corresponding to different design variables may begin to collapse into one another. While such configurations may be optimal from an inference or information–gain perspective, practical constraints, such as the need to enforce distinct sensor locations, can make them undesirable or infeasible. To enforce separations of design variables, a repulsive regularization term can be employed in addition to \(\mathcal{R}_v\). The choice of the corresponding regularization parameter \(\beta\) depends on the kernel used in \(\mathcal{R}_r\) and its parameters; in practice, however, \(\beta\) is typically small, in comparison to the parameter \(\alpha\). For instance \(\beta<10^{-4}\), for the Poisson problem in Section 6.1.
In this example, we take \(\mathcal{R}_r\) to be the MMD distance associated with a Gaussian kernel of standard deviation \(\sigma=0.009\). We examine the effect of \(\beta\) using the same one–dimensional Poisson problem as before. The batch size is set to \(B=8\), and the regularization parameter for \(\mathcal{R}_v\) is fixed to \(\alpha=0.05\). We use 1000 time steps of size \(4e-3\) and initialize 600 particles similarly as in the example 6.3.1. The left panel of Figure 7 shows the pairwise distances \(d = (\overline{x}_j-\overline{x}_l)^2, j,l = 1,...,B,\) \(j\neq l\) between the particle ensembles when the repulsive term \(\mathcal{R}_r\) is applied. The right panel compares the resulting empirical distributions, from which it is evident that when \(\beta\) is chosen too small the particle clouds collapse to only six optimal locations.
7 Conclusion↩︎
In this work, we established a rigorous connection between the frequentist batch-based A-optimal design problem formulated over measures and a corresponding Bayesian inference framework. By interpreting the relaxation to positive measures of fixed mass as arising from a non-parametric Gaussian observation model, we derived an infinite-dimensional expected utility functional and proved its concavity with respect to the design measure. This result provides a principled Bayesian justification for the convex relaxation previously proposed in the frequentist setting and clarifies the statistical meaning of optimizing over measures rather than discrete sensor configurations. The derivation of first variations and Wasserstein gradients further enabled the formulation of a particle-based optimization scheme grounded in gradient flow dynamics.
Beyond the theoretical characterization, we proposed a regularized tensorized formulation that lifts the optimization back to the product space while retaining the advantages of the relaxed representation. The introduced variance and maximum mean discrepancy penalties offer practical mechanisms to promote concentration of individual design measures and enforce separation between distinct sensor ensembles. Together, these elements yield a flexible framework for batch-based Bayesian experimental design over continuous domains. Future work may investigate convergence properties in greater generality, extensions to other optimality criteria, and applications to high-dimensional inverse problems arising in large-scale scientific and engineering systems.
Acknowledgement↩︎
This work was supported by the Research Council of Finland (decisions 353094, 348504, 359183). SM was supported by Emil Aaltonen foundation.