Towards optimal sensor placement for inverse problems in spaces of measures

Despite the popularity of sparse inverse problems, most of the existing work, to the best of our knowledge, focuses on deterministic noise ε. Central objects in this context, are the (noiseless) minimum norm problem

$$\begin{align} \min_{\mu \in \mathcal{M}\left(\Omega_{s}\right)} \vert\vert{\mu}\vert\vert_{\mathcal{M}\left(\Omega_{s} \right)} \quad \text{subject to } K \mu = K \mu^{\dagger} \end{align} \tag{ 𝒫0 }$$

as well as the question whether $\mu^\dagger$ is identifiable, i.e. its unique solution. A sufficient condition for the latter, is, e.g. the injectivity of the restricted operator $K_{\vert \mathrm{supp}\, \mu^\dagger}$ as well as the existence of a so-called dual certificate $\eta^{\dagger} \in \mathcal{C}^2(\Omega_s)$, [11], i.e. a subgradient $\eta^\dagger \in \partial \|\mu^\dagger\|_{\mathcal{M}(\Omega_s)}$, which is in some sense minimal, satisfying a strengthened source condition

$$\begin{align*} |\eta^\dagger\left(y\right)|\unicode{x2A7D} 1 \quad \text{for all}~y \in \Omega_s,~\eta^\dagger\left(y^\dagger_n\right) = \operatorname{sign}\left(q^\dagger_n\right), \quad |\eta^\dagger\left(y\right)|\lt 1 \quad \text{for all}~y \in \Omega_s \setminus \left\{y^\dagger\right\}^{N_s}_{n = 1}. \end{align*}$$

For example, in particular settings, the groundbreaking paper [6] shows that $\mu^\dagger$ is identifiable if the source locations $y^\dagger_n$ are sufficiently well separated. In this context, several manuscripts, see e.g. [1, 11, 13, 34] for a non-exhaustive list, study the approximation of an identifiable $\mu^\dagger$ by solutions to the Tikhonov-regularized problem

$$\begin{align} \bar{\mu}\left(\varepsilon\right) \in \mathfrak{M}\left(\varepsilon\right) : = \mathop{\mathrm{arg}\,\mathrm{min}}_{\mu\in\mathcal{M}\left(\Omega_s\right)} \left \lbrack \dfrac{1}{2} \vert\vert{K \mu - z^d\left(\varepsilon\right)}\vert\vert^2_{\Sigma_0^{-1}} + \beta\vert\vert{\mu}\vert\vert_{\mathcal{M}\left(\Omega_s\right)} \right\rbrack, \end{align} \tag{ 𝒫βϵ }$$

where Σ₀ is a positive definite diagonal matrix and the regularization parameter $\beta = \beta(\|\varepsilon\|)\gt0$ is adapted to the strength of the noise. This represents a challenging nonsmooth minimization problem over the infinite-dimensional and non-reflexive space of Radon measures. Moreover, due to its lack of strict convexity, its solutions are typically not unique. Under mild conditions on the choice of β, arbitrary solutions $\bar{\mu}(\varepsilon)$ approximate $\mu^\dagger$ in the weak*-sense as ε goes to zero. Moreover, it was shown in [11] that if the minimal dual certificate $\eta^\dagger$ associated to problem (₀ ) satisfies the strengthened source condition and its curvature does not degenerate around $y^\dagger_n$, $\bar{\mu}(\varepsilon)$ is unique and of the form

$$\begin{align*} \bar{\mu}\left(\varepsilon\right) = \sum^{N^\dagger_s}_{n = 1} \bar{q}_n \left(\varepsilon\right) \delta_{\bar{y}_n \left(\varepsilon\right)} \quad \text{with} \quad |\bar{q}_n \left(\varepsilon\right)-q^\dagger_n|+ \vert\vert{\bar{y}_n \left(\varepsilon\right)-y^\dagger_n}\vert\vert = \mathcal{O}\left(\|\varepsilon\|\right) \end{align*}$$

provided that $\|\varepsilon\|$ and β are small enough.

From a practical perspective, assuming knowledge on the norm of the error is very restrictive or even unrealistic and a statistical model for the measurement error is more appropriate. While the literature on deterministic sparse inversion is very rich, there are only few works dealing with randomness in the problem. We point out, e.g. [5] in which the authors consider additive i.i.d. noise stemming from a low-pass filtering of the signal. A reconstruction $\bar{\mu}(\varepsilon)$ is obtained by solving a constrained version of (_β, ) and the authors show that, with high probability, there holds $Q_{\text{hi}} (\bar{\mu}(\varepsilon)) \approx Q_{\text{hi}} (\mu^\dagger)$ where Q_hi is a convolution with a high-resolution kernel. Moreover, in [34] the authors consider deterministic noise but allow for randomness in the forward operator K. Their main result provides an estimate on an optimal transport energy between two positive measures derived from source and reconstruction. These again hold with high probability. Finally, we also mention [12] in which the authors propose a first step towards Bayesian inversion for sparse problems, i.e. both measurement noise as well as the unknown $\mu^\dagger$ are considered to be random variables. A suitable prior is constructed and well-posedness of the associated Bayesian inverse problem is shown.

In this paper, similar to [5], we adopt a frequentist viewpoint on sparse inverse problems and assume that the measurement errors follow a known probability distribution. In contrast, the unknown signal $\mu^\dagger$ is treated as a deterministic object. More in detail, we assume unbiased independent Gaussian noise with diagonal covariance matrix $\Sigma = \mathrm{diag}(\sigma_j)$, corresponding to independent measurements with variable quality sensors at different locations. We consider the Tikhonov-type estimator (_β, ) with

$$\begin{align*} \Sigma^{-1}_0 = \Sigma^{-1} / p, \quad\text{where } p = \,\mathrm{tr}\left(\Sigma^{-1}\right) \end{align*}$$

and investigate its error to the ground truth, where we have to account for the randomness of the noise. In statistical terms, Σ⁻¹ is the precision matrix of the sensor array, and p can be interpreted as an overall precision of the combined measurement, roughly representing an analogue to $1/\vert\vert{\varepsilon}\vert\vert$ in the stochastic setting. First and foremost, the uncertainty of the noise propagates to the estimator and thus $\bar{\mu}$ has to be interpreted as a random variable. Second, unlike the deterministic setting of [11], the asymptotic analysis cannot exclusively rely on smallness assumptions on the Euclidean norm of the noise: some realizations of ε might be very large, albeit with small probability. Consequently, reconstructions can exhibit undesirable features such as clustering phenomena around $y^\dagger_n$ or spurious sources far away from the true support. In particular, the reconstructed signal may comprise more or less than $N^\dagger_s$ sources. Thus, we require a suitable distance between signed measures that is compatible with weak* convergence on bounded subsets of $\mathcal{M}(\Omega_s)$ . We find a suitable candidate in generalizations of optimal transport energies [24]; see also [9, 34].

Despite its various difficulties, stochastic noise also provides new opportunities. For example, unlike the deterministic case, we are given a whole distribution of the measurement data and not only one particular realization. Clearly, the uncertainty in the estimate critically depends on the appropriate choice of measurement locations $\boldsymbol{x} = (x_j)_{j = 1,\ldots,N_o}$, the overall precision p, and relative precision of each sensor $\Sigma_0^{-1}$. Formalizing this connection enables the mathematical program of optimal sensor placement or optimal design, i.e. an optimization of the measurement setup to mitigate the influence of the noise before any data is collected in a real experiment. This requires a cheap-to-evaluate design criterion which allows to compare the quality of different sensor setups. For linear inverse problems in Hilbert spaces, a popular performance indicator is the mean-squared error of the associated least-squares estimator, which admits a closed form representation through its decomposition into variance and bias; see, e.g. [18]. For nonlinear problems, locally optimal sensor placement approaches rely on a linearization of the forward model around a best guess for the unknown parameters; see, e.g. [36]. To the best of our knowledge, optimal sensor placement for nonsmooth estimation problems and for infinite dimensional parameter spaces beyond the Hilbert space setting is uncharted territory.

Taking the mentioned difficulties in the stochastic setting into consideration, we are led to the analysis of the worst-case mean-squared-error of the estimator

$$\begin{align} \operatorname{MSE}\left[\bar{\mu}\right] : = \mathbb{E} \left \lbrack \sup_{\bar{\mu} \in \mathfrak{M}} d_{\mathrm{HK}}\left(\bar{\mu},\mu^\dagger \right)^2\right \rbrack = \int_{\mathbb{R}^{N_o}} \sup_{\bar{\mu} \in \mathfrak{M}\left(\varepsilon\right)} d_{\mathrm{HK}}\left(\bar{\mu},\mu^\dagger \right)^2\mathrm{d}\gamma_p\left(\varepsilon\right), \end{align} \tag{ 1.3 }$$

where d_HK denotes an extension of the Hellinger–Kantorovich distance introduced in [24] to signed measures (see section 4) and γ_p is the noise distribution $\mathcal{N}(0,\Sigma)$. We point out that, in comparison to linear inverse problems in Hilbert space, $\operatorname{MSE}[\bar{\mu}]$ does not admit a closed form expression and its computation requires both, a sampling of the expected value, as well as an efficient way to calculate the Hellinger–Kantorovich distance. This prevents its direct use in the context of optimal sensor placement for sparse inverse problems.

To enable efficient sensor design, we first need to select an appropriate regularization parameter, depending on the noise level. Here, we focus on the a priori choice rule of $\beta(p) = \beta_0 / \sqrt{p}$ for some tunable $\beta_0 \gt 0$, that only takes into account the overall precision of the sensor. For this choice, we provide the following upper bound:

$$\begin{align} \mathrm{MSE}\left[\bar{\mu}\right] \unicode{x2A7D} \frac{8}{p} \, \psi_{\beta_0}\left(\boldsymbol{x},\Sigma_0\right) + \bar{c} \exp \left(-\bar{\lambda}\beta_0^2\right), \end{align} \tag{ 1.4 }$$

where the constant $\psi_{\beta_0}(\boldsymbol{x},\Sigma_0)$ (further detailed below) explicitly depends on the locations and relative precisions while the constants $\bar{c}$ and $\bar{\lambda}$ depend on the problem setup (the kernel and domain, some basic bounds on the ground truth), a non-degeneracy parameter of the dual certificate $\eta^\dagger$ (further detailed below), and quantities that can be bounded by $\psi_{\beta_0}(\boldsymbol{x},\Sigma_0)$, but do not depend on p or β₀ for $p \unicode{x2A7E} \bar{p}\gt0$ and $\beta_0 \unicode{x2A7E} \bar{\beta}_0\gt0$; see theorem 6.1. Thus, under these basic assumptions and by choosing β₀ large enough, the second term in (1.4) becomes negligible and the first term dominates and closely predicts the mean-squared error. This behavior is confirmed by numerical examples; see section 7.

To further illustrate the meaning of the constant $\psi_{\beta_0}(\boldsymbol{x},\Sigma_0)$, let us denote by $\boldsymbol{q} = (q_1;\ldots; q_{N_s})$ and $\boldsymbol{y} = (y_1;\ldots; y_{N_s})$ the vectors of coefficients and positions of sources, respectively. Additionally, we collect all the parameters of a given finite source µ in the vector $\boldsymbol{m} = (\boldsymbol{q}; \boldsymbol{y})$, and introduce the parameter-to-observation map $G(\boldsymbol{m}) = K\mu$, as well as its Jacobian $G^{^{\prime}}(\boldsymbol{m}^\dagger)$ evaluated at the parameters of the ground truth. Associated to this, we denote the Fisher information matrix $\mathcal{I}_0$ by

$$\begin{align} \mathcal{I}_0 : = G^{^{\prime}}\left(\boldsymbol{m}^\dagger\right)^\top \Sigma_0^{-1} G^{^{\prime}}\left(\boldsymbol{m}^\dagger\right) . \end{align} \tag{ 1.5 }$$

Then the constant in the estimate above is computed as

$$\begin{align*} \psi_{\beta_0}\left(\boldsymbol{x},\Sigma_0\right) = \,\mathrm{tr}\left(W_{\dagger}\mathcal{I}_0^{-1}\right) + \beta_0^2 \left\| \mathcal{I}_0^{-1} \left(\boldsymbol{\rho};\boldsymbol{0}\right) \right\|^2_{W_{\dagger}}, \end{align*}$$

with the sign vector $\boldsymbol{\rho} = \mathrm{sign}\ \boldsymbol{q}^\dagger$ and a weighted Euclidean norm $\|\cdot\|_{W_\dagger}$ which is induced by a positive definite matrix $W_\dagger$ connected to the ground truth $\boldsymbol{m}^\dagger$. This clarifies how the multiplicative constant in the estimate explicitly depends on the measurement setup and we note that it closely resembles the 'classical' A-optimal design criterion; see [18]. Together with the estimate (1.4), and the smallness of the second term, this suggests that $\psi_{\beta_0}(\boldsymbol{x},\Sigma_0)$ is a suitable criterion to quantify the quality of a given design in terms of the MSE (1.3).

Concerning the smallness of the second term, we note that the constant $\bar{\lambda}$ also depends on a non-degeneracy constant θ > 0, which is a further tightening of the assumption on the dual certificate. This non-degenerate source condition on $\mu^\dagger$ requires the associated minimal norm dual certificate $\eta^\dagger$ to fulfill

$$\begin{align} \vert \eta^\dagger\left(y\right)\vert \unicode{x2A7D} 1 - \theta \min\left\{\theta,\, \min_{n = 1,\ldots,N_s} \left\| \sqrt{|q_n^{\dagger}|}\left(y-y_n^{\dagger}\right) \right\|_{2}^2 \right\} \quad\text{for all } y\in \Omega_s \end{align} \tag{ 1.6 }$$

for some θ > 0. This condition has been employed in many previous works, and is known to uniformly hold for several settings under general assumptions on the measurement and a separation of the condition of the sources; see, e.g. [33] and the references therein.

The proof of the main result relies on a splitting of the set of measurement errors $\mathbb{R}^{N_o}$ into a set of 'nice' events $\mathcal{A}_{\text{nice}}$ as well as an estimate of the probability of its complement $\mathbb{R}^{N_o} \setminus \mathcal{A}_{\text{nice}}$, related to the second term in (1.4). On $\mathcal{A}_{\text{nice}}$, there is a unique optimal parameter $\bar{\boldsymbol{m}} = (\bar{\boldsymbol{q}},\bar{\boldsymbol{y}})$ with the correct number of sources that parametrizes $\bar{\mu}$. Then, the distance between the reconstruction and the ground truth in the Hellinger–Kantorovich distance can be estimated by a weighted Euclidean distance of the parameters. Those can be further estimated with a linearization of G, which leads to (1.4) after explicitly computing the expectation. This estimate is specific to the choice of d_HK and relies on its interpretation as an unbalanced Wasserstein-2 distance. While similar estimates can be derived for other popular metrics such as the Kantorovich-Rubinstein distance (related to the Wasserstein-1 distance; see appendix C) this would introduce additional constants in the first term of (1.4) stemming from an inverse inequality of discrete $\ell_1$ and weighted $\ell_2$ norms. Thus, the first term in the modified estimate would overestimate the true error by a potentially substantial factor. In contrast, the first term in (1.4) is sharp in the sense that the convenient factor of 8 can, mutatis mutandis, be replaced by any c > 1, at the cost of increasing the constant in the second term.

1.4.1. Sparse minimization problems beyond inverse problems.

1.4.2. Inverse problems with random noise.

The paper is organized as follows: In section 3, we recall some properties of the minimum norm problem (_β, ) and the Tikhonov regularized problem (_β, ) as well as its solutions. In section 4, we define the Hellinger–Kantorovich distance and investigate its properties. Section 5 is devoted to study the linearized estimate $\delta \widehat{\boldsymbol{m}}$. Using these results, we then investigate sparse inverse problems with random noise in section 6 and provide a sharp upper bound for $\mathrm{MSE}[\bar{\mu}]$ in section 6.2. Finally, in section 7 we present some numerical examples to verify our theory.

Throughout the paper, $c_i, C_i$, $i = 1,2,\ldots$ denote generic constants that may vary from line to line. By $C = C(a,b,\ldots)$, we indicate that C depends on $a,b,\ldots$. We denote by $\Omega_s \subset \mathbb{R}^{d}$ and $\Omega_o \subset \mathbb{R}^{d_o}$ the compact location and observation set, where $d_o, d \unicode{x2A7E} 1$ and $\Omega_s$ has a nonempty interior. A vector in X^m for a set X and m > 1, will be written in bold face, for instance $\boldsymbol{y} = (y_1;\ldots;y_{N_s}) \in \Omega_s^{N_s}$, $\boldsymbol{q} = (q_1;\ldots;q_{N_s}) \in \mathbb{R}^{N_s}$ and $\boldsymbol{x} = (x_1;\ldots; x_{N_o}) \in \Omega_o^{N_o}$ are vectors of coefficients, positions of sources and positions of observations, respectively, where the formal definitions are introduced in the sequel. We write $(\boldsymbol{a}_1, \ldots, \boldsymbol{a}_n)$ and $(\boldsymbol{a}_1;\ldots; \boldsymbol{a}_n)$ to stack vectors $\boldsymbol{a}_1, \ldots, \boldsymbol{a}_n$ horizontally and vertically, respectively. We write $\left\| \cdot \right\|_p$ for the usual $\ell^p$-norm on $\mathbb{R}^m$. For a vector $x \in \mathbb{R}^m$ and a positively defined matrix $W \in \mathbb{R}^{m \times m}$, we define the weighted W-norm of x as $\left\| x \right\|_{W}: = \left\| W^{1/2}x \right\|_2$. The closed ball in this weighted norm is denoted by $B_W(x,r) : = \{\, x^{^{\prime}} \in \mathbb{R}^m\,:\, \left\| x^{^{\prime}} - x \right\|_W \unicode{x2A7D} r \,\}$. For a linear map $A: X \to Y$, the operator norm of A is given by $\left\| A \right\|_{X \to Y} = \sup_{\left\| x \right\|_X \unicode{x2A7D} 1} \left\| Ax \right\|_Y$. Similarly, any bilinear map $A: X_1 \times X_2 \to Y$ has a natural operator norm $\left\| A \right\|_{X_1 \times X_2 \to Y}: = \sup_{\left\| x_1 \right\|_{X_1} \unicode{x2A7D} 1, \left\| x_2 \right\|_{X_2} \unicode{x2A7D} 1} \left\| A(x_1,x_2) \right\|_{Y}.$

Furthermore, let $k : \Omega_o \times \Omega_s \to \mathbb{R}$ be a real-valued kernel. We introduce the following notations which turn k into vector-valued kernels: $k[\boldsymbol{x}](y) = k[\boldsymbol{x},y]$ is a column vector with

$$\begin{align} k\left[\boldsymbol{x},y\right] : = \left(k\left(x_1,y\right);\ldots; k\left(x_{N_o}, y\right)\right), \quad \boldsymbol{x} = \left(x_1;\ldots; x_{N_o}\right) \in \Omega_o^{N_o}, \quad y \in \Omega_s, \end{align} \tag{ 2.1 }$$

while $k[x,\boldsymbol{y}]$ is a row vector with

$$\begin{align} k\left[x,\boldsymbol{y}\right] : = \left(k\left(x,y_1\right), \ldots, k\left(x,y_{N_s}\right)\right), \quad x \in \Omega_o, \quad \boldsymbol{y} = \left(y_1;\ldots; y_{N_s}\right) \in \Omega_s^{N_s}. \end{align} \tag{ 2.2 }$$

Similarly, we also have the matrix $k[\boldsymbol{x}, \boldsymbol{y}]$ defined as

$$\begin{align} k\left[\boldsymbol{x}, \boldsymbol{y}\right] : = \left(k\left(x_1, \boldsymbol{y}\right); \ldots; k\left(x_{N_o}, \boldsymbol{y}\right)\right). \end{align} \tag{ 2.3 }$$

When $k = k(x,\cdot)$ is a smooth function in variable y, we consider the r^th-derivative of k the tensor of partial derivatives is y by $\nabla^r_{y\cdots y} k(x,y)$. In particular, $\nabla_y k(x,y)$ and $\nabla^2_{yy} k(x,y)$ are the gradient and Hessian of k (with respect to variable y,) respectively. We note that $\nabla_y k \colon \Omega_o \times\Omega_s \to \mathbb{R}^{N_s}$ is a vector valued kernel and thus we define $\nabla_y^\top k[x,\boldsymbol{y}]$ as a matrix defined by

$$\begin{align} \nabla_y^\top k\left[x,\boldsymbol{y}\right] & = \left(\nabla_{y}k\left(x,y_1\right)^\top, \nabla_{y}k\left(x,y_2\right)^\top,\ldots, \nabla_{y}k\left(x,y_{N_s}\right)^\top\right). \end{align} \tag{ 2.4 }$$

Similarly, $\nabla_y^\top k[\boldsymbol{x},\boldsymbol{y}]$ is a block matrix defined by

$$\begin{align} \nabla_y^\top k\left[\boldsymbol{x},\boldsymbol{y}\right] & = \left(\nabla_y^\top k\left[x_1,\boldsymbol{y}\right], \ldots, \nabla_y^\top k\left[x_{N_o},\boldsymbol{y}\right]\right). \end{align} \tag{ 2.5 }$$

Throughout the paper, by a slight abuse of notation, we denote by ε a variable deterministic noise, a random variable, or its realization, which will be clear from the context. By γ_p we denote the density of a multivariate Gaussian random variable with expectation zero and covariance Σ. Further notation, specific to the present manuscript, will be introduced at first appearance. For quicker reference, a notation table can be found in appendix D.

We also recall some basic facts and assumptions for inverse source location.

2.2.1. Integral kernels.

Throughout the paper, we assume that the kernel is sufficiently regular:

• (A1)  
  The kernel $k \in \mathcal{C}(\Omega_o \times \Omega_s)$ is three-times differentiable in the variable y. For abbreviation, we further set

$$\begin{align*} C_{k} &: = \sup_{x\in\Omega_o,y\in\Omega_s} \vert k\left(x,y\right)\vert, & C^{^{\prime}}_{k} &: = \sup_{x\in\Omega_o,y\in\Omega_s} \left\| \nabla_y k\left(x,y\right) \right\|_2, \\ C^{^{\prime\prime}}_{k} &: = \sup_{x\in\Omega_o,y\in\Omega_s} \left\| \nabla^2_{yy} k\left(x,y\right) \right\|_{2 \to 2}, & C^{^{\prime\prime\prime}}_{k} &: = \sup_{x\in\Omega_o,y\in\Omega_s} \left\| \nabla^3_{yyy} k\left(x,y\right) \right\|_{2 \times 2 \to 2}. \end{align*}$$

By means of the kernel k, we introduce the weak* continuous source-to-measurements operator $K \colon \mathcal{M}(\Omega_s) \to \mathbb{R}^{N_o}$ with

$$\begin{align} K\mu = \left(\int_{\Omega_s}k\left(x_1,y\right)\mathrm{d} \mu\left(y\right);\ldots;\int_{\Omega_s}k\left(x_{N_o},y\right)\mathrm{d} \mu\left(y\right)\right). \end{align} \tag{ 2.6 }$$

Moreover, consider the operator $K^* \colon \mathbb{R}^{N_o} \to \mathcal{C}^2(\Omega_s)$ given by

$$\begin{align} \lbrack K^*z \rbrack \left(y\right) = \sum^{N_o}_{j = 1} z_j k\left(x_j,y\right) \quad \text{for all} \quad z \in \mathbb{R}^{N_o}. \end{align} \tag{ 2.7 }$$

Then $K^*$ is linear and continuous and there holds

$$\begin{align*} \int_\Omega \lbrack K^*z \rbrack \left(y\right)~\mathrm{d}\mu\left(y\right) = z^\top \lbrack K\mu \rbrack \quad \text{for all}\quad \mu \in \mathcal{M}\left(\Omega_s\right),~z \in \mathbb{R}^{N_o}. \end{align*}$$

2.2.2. Space of Radon measures.

We recall some properties of Radon measures. Let $\Omega \subset \mathbb{R}^d$, $d \unicode{x2A7E} 1$ be a compact set. We define the space of Radon measures $\mathcal{M}(\Omega)$ as the topological dual of the space $\mathcal{C}(\Omega)$ of continuous functions on Ω endowed with the supremum norm. It is then a Banach space equipped with the dual norm

$$\begin{align*} \left\| \mu \right\|_{\mathcal{M}\left(\Omega\right)} : = \sup \left\{\int_{\Omega} f \mathrm{d} \mu: f \in \mathcal{C}\left(\Omega\right), \left\| f \right\|_{\mathcal{C}\left(\Omega\right)} \unicode{x2A7D} 1 \right\}. \end{align*}$$

Weak* convergence of a sequence in $\mathcal{M}(\Omega)$ will be denoted by '$\rightharpoonup^*$'. More specifically, we have

$$\begin{align*} \mu_n \rightharpoonup^* \mu \quad\text{if and only if} \quad \int_{\Omega} f \mathrm{d} \mu_n \to \int_{\Omega} f \mathrm{d}\mu \quad \text{for all} \quad f \in \mathcal{C}\left(\Omega\right). \end{align*}$$

Next, by the definition of the total variation norm, its subdifferential is defined by

$$\begin{align*} \partial \left\| \mu \right\|_{\mathcal{M}\left(\Omega_s\right)} : = \left\{\eta \in \mathcal{C}\left(\Omega_s\right): |\eta\left(y\right)| \unicode{x2A7D} 1,\forall y \in \Omega_s \text{ and } \int_{\Omega_s} \eta \mathrm{d}\mu = \left\| \mu \right\|_{\mathcal{M}\left(\Omega_s\right)} \right\}, \end{align*}$$

see for instance [11]. In particular, for a discrete measure $\mu = \sum_{n = 1}^{N} q_n\delta_{y_n}$ one has

$$\begin{align*} \partial \left\| \mu \right\|_{\mathcal{M}\left(\Omega_s\right)} = \left\{\eta \in \mathcal{C}\left(\Omega_s\right): |\eta\left(y\right)| \unicode{x2A7D} 1,\forall y \in \Omega_s \,\textrm{ and } \eta\left(y_n\right) = \mathrm{sign}\left(q_n\right),\forall n = 1,\ldots, N \right\}. \end{align*}$$

Finally, by $\mathcal{M}^+(\Omega)$ we refer to the set of positive Radon measures on Ω.

Sparse inverse problems appear in a variety of interesting applications. In the following, we give some examples which fit into our setting.

Example 3.1. Consider the advection-diffusion equation

$$\begin{align} \partial_t u - \nabla\left(\boldsymbol{D}\cdot \nabla u\right) + \nabla \cdot \left(\kappa u\right) = 0\, \textrm{ in}\,\left(0,T\right) \times \mathbb{R}^d, \end{align} \tag{ 3.2 }$$

together with the initial value $u(0,\cdot) = \mu$. The boundary condition is given by u → 0 as $x \to \infty$. This equation describes the rate of change of the concentration of the contaminant $u(t,x)$. For simplicity, we consider a two-dimensional medium, and both $\kappa = (\kappa_1,\kappa_2)$ and $\boldsymbol{D} = \mathrm{diag}(D_1,D_2)$ are independent of x. Here the solution to (3.2) is given by

$$\begin{align*} u\left(t,x\right) = \int_{\mathbb{R}^2} G\left(x - y, t\right) \mathrm{d}\mu\left(y\right) \end{align*}$$

where $G(x,t)$ is the Green's function of the advection-diffusion equation, which is given by

$$\begin{align*} G\left(x,t\right) = \frac{1}{4\pi \sqrt{D_{1}D_{2}t}} \exp\left(-\left\| x - \kappa t \right\|_{\boldsymbol{D}^{-1}}^2/\left(4t\right)\right). \end{align*}$$

Here, if one seeks to identify the initial value µ from finite number of measurements at time $T_o \gt 0$ in the observation set $\Omega_o \subset \mathbb{R}^2$, the kernel is given by $k(x,y) = G(x-y, T_o)$.

Example 3.2. Consider the advection-diffusion equation on a bounded smooth domain Ω, together with the Dirichlet boundary conditions $u\rvert_{(0,T) \times \partial \Omega} = 0$, then there exists a kernel $G(x,y,t)$ such that

$$\begin{align*} u\left(t,x\right) = \int_{\Omega} G\left(x, y, t\right) \mathrm{d}\mu\left(y\right), \end{align*}$$

see, e.g. [15]. In this case, for observations at time T_o we choose $k = G(\cdot, \cdot, T_o)$. For $\Omega_o \subset \Omega$ (i.e. no observation near the boundary), the regularity requirements on $\partial\Omega$ are not necessary since one can employ interior regularity arguments; see, e.g. [17].

In this section, we briefly summarize some preliminary results concerning the regularized problem (_β, ) as well as its solution set. We start by discussing its well-posedness.

Proposition 3.3. Problem (_β, ) admits a solution $\bar{\mu}$. Furthermore, any solution $\bar{\mu}$ to (_β, ) satisfies $\left\| \bar{\mu} \right\|_{\mathcal{M}(\Omega_s)} \unicode{x2A7D} \left\| \varepsilon \right\|_{\Sigma_0^{-1}}^2/(2\beta) + \left\| \mu^\dagger \right\|_{\mathcal{M}(\Omega_s)}$ and the solution set

$$\begin{align*} \mathfrak{M}\left(\varepsilon\right) = \mathrm{arg}\,\mathrm{min}\, \left(\mathcal{P}_{\beta,\varepsilon}\right) \end{align*}$$

is weak* compact.

Proof. Existence of a minimizer of (_β, ) is guaranteed by [4, proposition 3.1] noticing that the forward operator $K: \mathcal{M}(\Omega_s) \to \mathbb{R}^{N_o}$ of (_β, ) is weak*-to-strong continuous. For the upper bound we use the optimality of $\bar{\mu}$ compared to $\mu^\dagger$ as well as the definition of $z^d(\varepsilon)$ to get

$$\begin{align*} \beta\left\| \bar{\mu} \right\|_{\mathcal{M}\left(\Omega_s\right)}\unicode{x2A7D} \frac{1}{2}\left\| K \bar{\mu} - z^d \right\|_{\Sigma_0^{-1}}^2 + \beta\left\| \bar{\mu} \right\|_{\mathcal{M}\left(\Omega_s\right)} \unicode{x2A7D} \frac{1}{2}\left\| \varepsilon \right\|_{\Sigma_0^{-1}}^2 + \beta\left\| \mu^\dagger \right\|_{\mathcal{M}\left(\Omega_s\right)}. \end{align*}$$

Moreover, $\mathfrak{M}(\varepsilon)$ is weak* closed since the objective functional in (_β, ) is weak* lower semicontinuous. Combining both observations, we conclude the weak* compactness of $\mathfrak{M}(\varepsilon)$. □

In particular, note that $\mathfrak{M}(\varepsilon)$ is, in general, not a singleton due to the lack of strict convexity in (_β, ). Moreover, we recall that the inverse problem was introduced as a lifting of the nonconvex and combinatorial integral equation (1.1). From the same perspective, (_β, ) can be interpreted as a convex relaxation of the parametrized problem

$$\begin{align} \inf_{\substack{\boldsymbol{y}\in \Omega^N_s,~\boldsymbol{q} \in \mathbb{R}^N, \\ N \in \mathbb{N}}} \left \lbrack \dfrac{1}{2} \left\| k\left[\boldsymbol{x},\boldsymbol{y}\right]\boldsymbol{q} - z^d \right\|_{\Sigma_0^{-1}}^2 + \beta \left\| \boldsymbol{q} \right\|_{1} \right\rbrack, \end{align} \tag{ 3.3 }$$

In the following proposition, we show that this relaxation is exact, i.e. there exists at least one solution to (3.3) and its minimizers parametrize sparse solutions to (_β, ).

Proposition 3.4. There holds $\mathrm{min}$ (_β, ) $= \,\mathrm{inf}$ (3.3). For a triple $(\bar{N},\bar{\boldsymbol{y}},\bar{{\boldsymbol{q}}})$ with $\bar{y}_i \neq \bar{y}_j$, i ≠ j, the following statements are equivalent:

• The triple $(\bar{N},\bar{\boldsymbol{y}},\bar{{\boldsymbol{q}}})$ is a solution of (3.3).
• The parametrized measure $\bar{\mu} = \sum^{\bar{N}}_{n = 1} \bar{q}_n \delta_{\bar{y}_n}$ is a solution of (_β, ).

Moreover, (_β, ) admits at least one solution of this form with $\bar{N}\unicode{x2A7D} N_o$.

Proof. Given $({N},{\boldsymbol{y}},{{\boldsymbol{q}}})$ with ${y}_i \neq {y}_j$, i ≠ j, note that the sparse measure

$$\begin{align*} \mu\left({\boldsymbol{y}},{{\boldsymbol{q}}}\right) = \sum^N_{n = 1} q_n \delta_{y_n} \quad \text{satisfies} \quad K\mu\left({\boldsymbol{y}},{{\boldsymbol{q}}}\right) = k\left[\boldsymbol{x},\boldsymbol{y}\right]\boldsymbol{q},~ \|\mu\left({\boldsymbol{y}},{{\boldsymbol{q}}}\right)\|_{\mathcal{M}\left(\Omega_s\right)} = \left\| \boldsymbol{q} \right\|_{1}. \end{align*}$$

Hence, one readily verifies min (_β, ) $= \, \mathrm{inf}$ (3.3) as well as the claimed equivalence due to the weak* density of the set of sparse measures in $\mathcal{M}(\Omega_s)$ and since the objective functional in (_β, ) is weakly* lower semicontinuous. The existence of a sparse solution to (_β, ) follows similarly to [30, theorem 3.7]. □

The equivalence between both of these problems will play a significant role in our subsequent analysis. Additional insight on the structure of solutions to (_β, ) can be gained through the study of its first order necessary and sufficient optimality conditions. Since our interest lies in sparse solutions, we restrict the following proposition to this particular case.

Proposition 3.5. A measure $\bar{\mu} = \sum_{n = 1}^{\bar N} \bar{q}_n \delta_{\bar{y}_n}$ is a solution of (_β, ) if and only if

$$\begin{align*} |\bar{\eta}\left(y\right)| \unicode{x2A7D} 1 \;\text{for all } y \in \Omega_s, \quad \bar{\eta}\left(\bar{y}_n\right) = \operatorname{sign}\left(\bar{q}_n\right), \quad \forall n = 1, \ldots, \bar{N}, \end{align*}$$

where

$$\begin{align} \bar{\eta} = -K^*\Sigma_0^{-1}\left(K \bar{\mu} - z^d\right)/\beta = K^*\Sigma_0^{-1} \left(z^d - k\left[\boldsymbol{x},\bar{\boldsymbol{y}}\right]\bar{\boldsymbol{q}}\right)/\beta. \end{align} \tag{ 3.4 }$$

Note that $\bar{\eta}$ is independent of the particular choice of the solution to (_β, ). We will refer to it as the dual certificate associated to (_β, ) in the following. Finally, we give a connection between (_β, ) and the minimum norm problem (₀ ) in the vanishing noise limit. The following general convergence property follows directly from [19].

Proposition 3.6. Assume that $\beta = \beta(\varepsilon)$ is chosen such that

$$\begin{align*} \beta \to 0 \,\textrm{ and } \dfrac{\left\| \varepsilon \right\|_{\Sigma_0^{-1}}^2}{\beta} \to 0\, \textrm{ as }\, \left\| \varepsilon \right\|_{\Sigma_0^{-1}} \to 0. \end{align*}$$

Then solutions to (_β, ) subsequentially converge weakly-* towards solutions of (₀ ).

Following proposition 3.6, guaranteed recovery of the ground truth measure requires that $\mu^\dagger$ is identifiable, i.e. the unique solution of (₀ ). In this section, we briefly summarize some key concepts regarding (₀ ) and state sufficient assumptions for the latter. For this purpose, introduce the associated Fenchel dual problem

$$\begin{align} \min_{\zeta \in \mathbb{R}^{N_o}} \left \lbrack -\langle \mu^\dagger, K^* \Sigma_0^{-1} \zeta \rangle + \mathbb{I}_{\|K^* \Sigma_0^{-1} \zeta\|_{C\left(\Omega_s\right)}\unicode{x2A7D} 1} \right\rbrack. \end{align} \tag{ 3.5 }$$

as well as the minimal-norm dual certificate

$$\begin{align} \eta^{\dagger} : = K^*\Sigma^{-1}_0 \zeta^\dagger \in \mathcal{C}^2\left(\Omega_s\right) \quad \text{where} \quad \zeta^\dagger = \mathop{\mathrm{arg}\,\mathrm{min}}_{\zeta \in \mathbb{R}^{N_o}} \left\{\,\|\zeta\|_2\;:\; \zeta \in \mathrm{arg}\,\mathrm{min}\, \left(3.5\right)\,\right\}. \end{align} \tag{ 3.6 }$$

Note that the existence of $\zeta^{\dagger}$, and therefore the minimum-norm dual certificate $\eta^{\dagger}$, is guaranteed in this setting following [30, proposition A.2] as well as due to $K^* \colon \mathbb{R}^{N_o} \to \mathcal{C}^2(\Omega_s)$. Moreover, by standard results from convex analysis, a given $\mu \in \mathcal{M}(\Omega_s)$ is a solution to (₀ ) if and only if $\eta^\dagger \in \partial \|\mu\|_{\mathcal{M}(\Omega_s)}$. The following assumptions on $\mu^\dagger$ and $\eta^\dagger$ are made throughout the paper:

• (A2)  
  Structure of $\mu^\dagger$ : We assume that there holds

  $$\begin{align} \mu^\dagger = \sum^{N_s^{\dagger}}_{n = 1} q^\dagger_n \delta_{y^\dagger_n} \quad \text{where} \quad q^\dagger_n \neq 0,~y^\dagger_n \in \operatorname{int}\left(\Omega_s\right) \quad \text{for all} \quad n = 1, \ldots, N_s^{\dagger}. \end{align}$$
• (A3)  
  Source condition: We assume that the minimum-norm dual certificate $\eta^{\dagger}$ satisfies

  $$\begin{align*} |\eta^\dagger\left(y\right)| \unicode{x2A7D} 1 \quad \text{for all} \quad y \in \Omega_s \quad \text{and} \quad \eta^\dagger\left(y^\dagger_n\right) = \mathrm{sign}\left(q^\dagger_n\right) \quad \text{for all} \quad n = 1,\dots,N_s. \end{align*}$$
• (A4)  
  Strengthened source condition: We assume that

  $$\begin{align*} |\eta^\dagger\left(y\right)| < 1 \quad \text{for all} \quad y \in \Omega_s \setminus \left\{y^\dagger_n\right\}^{N^\dagger_s}_{n = 1} \end{align*}$$

  and the operator $K_{|\mathrm{supp}\, \mu^\dagger} : = k[\boldsymbol{x},\boldsymbol{y}^\dagger]$ is injective.

Here, assumption A3 is equivalent to $\eta^\dagger \in \partial \left\| \mu^\dagger \right\|_{\mathcal{M}(\Omega_s)}$, i.e. $\mu^\dagger$ is indeed a solution to (₀ ), whereas assumptions A2 and A3 imply its uniqueness. While assumption A4 seems very strong at first glance, it can be explicitly verified in some settings (see, e.g. [6]) and is often numerically observed in practice. According to [11, proposition 5] we have the following:

Proposition 3.7. Let assumptions A2–A4 hold. Then $\mu^\dagger$ is the unique solution of (₀ ).

As a consequence, proposition 3.6 implies $\bar{\mu} \rightharpoonup^* \mu^\dagger$. Moreover, according to [11, Proposition 1], the dual certificates $\bar{\eta}$ associated to (_β, ) approximate the minimal norm dual certificate $\eta^\dagger$ in a suitable sense. Taking into account assumption A3 as well as proposition 3.5, we thus conclude that the reconstruction of $\mu^\dagger$ from (3.3) is governed by the convergence of the global extrema of $\bar{\eta}$ towards those of $\eta^\dagger$. However, in order to capitalize on this observation in our analysis, we need to compute a closed form expression for $\eta^\dagger$. In general, this is intractable due to the global constraint $|\eta^{\dagger}(z)| \unicode{x2A7D} 1$, $z \in \Omega_s$. As a remedy, the authors of [11] introduce a simpler proxy replacing this constraint by finitely many linear ones noting that

$$\begin{align*} \nabla \eta^\dagger\left(y^\dagger_n\right) = 0, \quad \eta^\dagger\left(y^\dagger_n\right) = \mathrm{sign}\left(q^\dagger_n\right) \quad \text{for all } \quad n = 1,\ldots, N_s^{\dagger}. \end{align*}$$

The computation of the associated vanishing derivative pre-certificate $\eta_{\mathrm{PC}} : = K^* \Sigma_0^{-1} \zeta_{\mathrm{PC}} \in \mathcal{C}^2(\Omega_s)$ where

$$\begin{align} \zeta_{\mathrm{PC}} = \mathop{\mathrm{arg}\,\mathrm{min}}_{\zeta \in \mathbb{R}^{N_o}} \left\{\left\| \zeta \right\|_2: \nabla \eta_{\mathrm{PC}}\left(y^\dagger_i\right) = 0, \quad \eta_{\mathrm{PC}}\left(y^\dagger_n\right) = \operatorname{sign}\left(q^\dagger_n\right) \quad \text{for all } \quad n = 1,\ldots, N_s^{\dagger} \right\} \end{align} \tag{ 3.7 }$$

only requires the solution of a linear systems of equations and coincides with $\eta^\dagger$ under appropriate conditions, see [11, proposition 7]. Finally, in order to derive quantitative statements on the reconstruction error between $\bar{\mu}$ and $\mu^\dagger$, we require the non-degeneracy of the minimal norm dual certificate of $\mu^\dagger$ in the sense of [11]. Since we aim to use (1.4) in the context of optimal sensor placement, that is, we need to track the dependence of the involved constants on the measurement setting, we utilize the following quantitative definition; see [33].

Definition 3.8. We say that $\eta \in \mathcal{C}^2(\Omega_s)$ is $\theta-$non-degenerate or $\theta-$admissible for the sparse measure $\mu = \sum_{n = 1}^{N_s} q_n \delta_{y_n}$ and $\theta \in (0,1]$ if there holds

$$\begin{align} \vert \eta\left(y\right)\vert \unicode{x2A7D} 1 - \theta \min\left\{\theta,\, \min_{n = 1,\ldots,N_s} \left\| w^\dagger_n\left(y-y_n\right) \right\|_{2}^2 \right\}, \quad \eta\left(y_n\right) = \mathrm{sign}\left(q_n\right) \quad \quad\text{for all } \quad y\in \Omega_s \end{align} \tag{ 3.8 }$$

and weights $w^\dagger_n = \sqrt{|q_n^{\dagger}|}$.

Due to the regularity of η one readily verifies that (3.8) is equivalent to

$$\begin{align} -\mathrm{sign}{\eta}\left(y_n\right) \nabla^2 \eta\left(y_n\right) \unicode{x2A7E} 2\theta |w_n^{\dagger}|^2\operatorname{Id} \quad \text{for every} \quad n = 1,2,\ldots, N_s, \end{align} \tag{ 3.9 }$$

as well as

$$\begin{align} \vert \eta\left(y\right)\vert \unicode{x2A7D} 1-\theta^2, \quad \text{for all} \quad y \in \Omega_s \setminus \bigcup_{n = 1,\ldots,N_s} B_{w_n^{\dagger}}\left(y_n, \sqrt{\theta}\right). \end{align} \tag{ 3.10 }$$

In this section, we show that for sufficiently small noise ε, problem (5.1) admits a unique stationary point $\widehat{\boldsymbol{m}}(\varepsilon)$ in the vicinity of $\boldsymbol{m}^\dagger$. Moreover, loosely speaking, $\boldsymbol{m}^\dagger$ and $\boldsymbol{m}^\dagger+\delta \widehat{\boldsymbol{m}}(\varepsilon)$ provide Taylor expansions of zeroth and first order, respectively, for $\widehat{\boldsymbol{m}}(\varepsilon)$.

Proposition 5.2. Suppose that $G^{^{\prime}}(\boldsymbol{m}^\dagger)$ has full column rank. Then, for some constant $C_1 = C_1(k,\mu^{\dagger},\left\| \mathcal{I}_0^{-1} \right\|_{W_{\dagger}^{-1} \to W_{\dagger}})$ and radius $\hat{r}\gt0$ and all ε with $C_1(\left\| \varepsilon \right\|_{\Sigma^{-1}_0} + \beta) \unicode{x2A7D} 1$, the stationarity condition (5.1) admits a unique solution $\widehat{\boldsymbol{m}} = \widehat{\boldsymbol{m}}(\varepsilon)$ on ${B_{W^{\dagger}}}(\boldsymbol{m}^{\dagger},(3/2)\hat{r})$. Moreover, the stationary point satisfies $\widehat{\boldsymbol{m}} \in {B_{W^{\dagger}}}(\boldsymbol{m}^{\dagger},\hat{r})$ as well as

$$\begin{align*} \left\| \widehat{\boldsymbol{m}} - \boldsymbol{m}^\dagger \right\|_{W_\dagger} &\unicode{x2A7D} 2 \left\| \delta \widehat{\boldsymbol{m}} \right\|_{W_\dagger} \unicode{x2A7D} C_1 \left(\left\| \varepsilon \right\|_{\Sigma_0^{-1}} + \beta\right), \\ \left\| \widehat{\boldsymbol{m}} - \boldsymbol{m}^\dagger - \delta\widehat{\boldsymbol{m}} \right\|_{W_\dagger} &\unicode{x2A7D} C_1^2 \left(\left\| \varepsilon \right\|_{\Sigma_0^{-1}} + \beta\right)^2. \end{align*}$$

For the sake of brevity, we omit by now the proof of proposition 5.2, which is then presented in appendix B.

Remark 5.3. We note that C₁ depends monotonically on the norm of the inverse Fisher information matrix; see Remark B.3. Moreover, the dependency on the ground truth $\mu^\dagger$ is only in terms of the norm $\left\| \boldsymbol{q}^\dagger \right\|_{1}$, and distances of $y_n^\dagger$ to the boundary and $q_n^\dagger$ to zero.

As mentioned in the preceding section, solving the stationarity equation (5.4) for $\widehat{\boldsymbol{m}} = (\widehat{\boldsymbol{y}},\widehat{\boldsymbol{q}})$ is not feasible in practice since it presupposes knowledge of $N^\dagger_s$. Moreover, recalling that $\widehat{\boldsymbol{m}}$ is merely a stationary point, the parametrized measure

$$\begin{align} \widehat{\mu} = \sum^{N^\dagger_s}_{n = 1} \widehat{q}_n \delta_{\widehat{y}_n} \end{align} \tag{ 5.7 }$$

is not necessarily a minimizer of (_β, ). In this section, our primary goal is to show that $\widehat{\boldsymbol{m}}$ indeed parametrizes the unique solution of problem (_β, ) if the minimum norm dual certificate $\eta^\dagger$ associated to (₀ ) is θ-admissible and if the set of admissible noises ε is further restricted. A fully-explicit estimate for the reconstruction error between $\widehat{\mu}$ and the ground truth $\mu^\dagger$ in the Hellinger–Kantorovich distance then follows immediately. For this purpose, recall from [11, proposition 7] that the non-degeneracy of $\eta^\dagger$ implies

$$\begin{align} \eta^\dagger = \eta_{\text{PC}} = K^* \Sigma^{-1/2}_0 \left(G^{^{\prime}}\left(\boldsymbol{m}^\dagger\right) \Sigma^{-1/2}_0\right)^+ \left(\boldsymbol{\rho};\boldsymbol{0}\right) = K^* \Sigma_0^{-1} G^{^{\prime}}\left(\boldsymbol{m}^\dagger\right)\mathcal{I}_0^{-1}\left(\boldsymbol{\rho};\boldsymbol{0}\right). \end{align} \tag{ 5.8 }$$

where η_PC denotes the vanishing derivative pre-certificate from section 3.3.

We first prove that

$$\begin{align*} \widehat{\eta} = \beta^{-1}K^*\Sigma_0^{-1}\left(z^d\left(\varepsilon\right)-K\widehat{\mu}\right) = \beta^{-1}K^*\Sigma_0^{-1}\left(G\left(\boldsymbol{m}^\dagger\right)+\varepsilon-G\left(\widehat{\boldsymbol{m}}\right)\right) \end{align*}$$

is $\theta/2$-admissible for certain ε and β.

Proposition 5.4. Let the assumptions in proposition 5.2 be satisfied and $\eta^\dagger$ be $\theta-$admissible for $\mu^{\dagger}$, $\theta \in (0,1]$. Then there exists a constant $C_2 = C_2(k,\mu^{\dagger},\left\| \mathcal{I}_0^{-1} \right\|_{W_{\dagger}^{-1} \to W_{\dagger}})$ such that if

$$\begin{align} &C_1\left(\left\| \varepsilon \right\|_{\Sigma^{-1}_0} + \beta\right) \unicode{x2A7D} \sqrt{\theta/32}, \end{align} \tag{ 5.9 }$$

$$\begin{align} &C_2 \beta^{-1}\left(\left(\left\| \varepsilon \right\|_{\Sigma^{-1}_0} + \beta\right)^2 + \left\| \varepsilon \right\|_{\Sigma_0^{-1}}\right) \unicode{x2A7D} \theta^2/32, \end{align} \tag{ 5.10 }$$

then the function $\widehat{\eta}$ is $\theta/2-$admissible for $\widehat{\mu}$.

The proof of proposition 5.4 is then provided in appendix B.

Remark 5.5. We note that C₂ depends monotonically on the norm of the inverse Fisher information matrix; see remark B.5. Moreover, the dependency on the ground truth $\mu^\dagger$ is only in terms of the norm $\left\| \boldsymbol{q}^\dagger \right\|_{1}$, and distances of $y_n^\dagger$ to the boundary and $q_n^\dagger$ to zero.

As a consequence, we conclude that the solution to (_β, ) is unique and parametrized by $\widehat{\boldsymbol{m}}$. Moreover, its H-K distance to $\mu^\dagger$ can be bounded in terms of the linearization $\delta \widehat{\boldsymbol{m}}$.

Theorem 5.6. Let the assumptions of proposition 5.4 hold. Then the solution of (_β, ) is unique and given by $\widehat{\mu}$ from (5.7). There holds

$$\begin{align} d_{\mathrm{HK}}\left(\widehat \mu,\mu^\dagger\right)^2 \unicode{x2A7D} 8 \left\| \delta \widehat{\boldsymbol{m}} \right\|^2_{W_\dagger}. \end{align} \tag{ 5.11 }$$

Proof. From proposition 5.4, we conclude that $\widehat{\eta}$ is $\theta/2$-admissible for $\widehat{\mu}$. Consequently, we have $\widehat{\eta} \in \partial \|\widehat{\mu}\|_{\mathcal{M}(\Omega_s)}$, i.e. $\widehat{\mu}$ is a solution of (_β, ). It remains to show its uniqueness. For this purpose, it suffices to argue that

$$\begin{align*} K_{| \mathrm{supp}\, \widehat{\mu} } = k\lbrack \boldsymbol{x}, \widehat{\boldsymbol{y}} \rbrack \in \mathbb{R}^{N_o \times N^\dagger_s} \end{align*}$$

is injective, see, e.g. the proof of [31, proposition 3.6]. Assume that this is not the case. Then, following [30, theorem B.4], there is $\boldsymbol{v}\neq 0$ with $k\lbrack \boldsymbol{x}, \boldsymbol{y} \rbrack \boldsymbol{v} = 0$ and τ ≠ 0 such that the measure $\tilde{\mu}$ parametrized by $\tilde{\boldsymbol{m}} = (\tilde{\boldsymbol{q}};\widehat{\boldsymbol{y}})$ with $\tilde{\boldsymbol{q}} = \widehat{\boldsymbol{q}} + \tau\boldsymbol{v}$ is also a solution of (_β, ) (choose the sign of τ to not increase the $\ell_1$-regularization, and the magnitude small not to change the sign of $\tilde{\boldsymbol{q}}$) and $\tilde{\boldsymbol{q}} \neq \widehat{\boldsymbol{q}}$. For $s \in (0,1)$, set $\boldsymbol{q}_s = (1-s) \widehat{\boldsymbol{q}} + s \tilde{\boldsymbol{q}}$. By convexity of (_β, ), the measure parametrized by $\boldsymbol{m}_s = ( \boldsymbol{q}_s;\widehat{\boldsymbol{y}})$ is also a minimizer of (_β, ). Consequently, m _s is a solution of (5.1) and thus also a stationary point. Finally, noting that $\boldsymbol{m}_s \neq \widehat{\boldsymbol{m}}$, $s \in (0,1)$, and $\lim_{s\rightarrow 0} \boldsymbol{m}_s = \widehat{\boldsymbol{m}}$, we arrive at a contradiction to the uniqueness of stationary points in the vicinity of $\boldsymbol{m}^\dagger$. The estimate in (5.11) immediately follows from

$$\begin{align*} d_{\mathrm{HK}}\left(\widehat{\mu},\mu^\dagger\right)^2 \unicode{x2A7D} R\left(\widehat{\boldsymbol{q}},\boldsymbol{q}^\dagger \right) \left\| \widehat{\boldsymbol{m}}-\boldsymbol{m}^\dagger \right\|^2_{W_\dagger} \unicode{x2A7D} 2 \left\| \widehat{\boldsymbol{m}}-\boldsymbol{m}^\dagger \right\|^2_{W_\dagger} \unicode{x2A7D} 8 \left\| \delta \widehat{\boldsymbol{m}} \right\|^2_{W_\dagger} \end{align*}$$

 □

Before stating the main result of the manuscript, let us briefly motivate the particular choice of the misfit term in (_β, ) as well as the employed parameter choice rule from the perspective of the stochastic noise model. Since we consider independent measurements, their covariance matrix $\Sigma = p^{-1}\Sigma_0$ is diagonal with $\Sigma_{jj} = \sigma_j^2$ for variances $\sigma^2_j \gt 0$, $j = 1,\ldots,N_o$. This corresponds to performing the individual measurements with independent sensors of variable precision $p_j = 1/\sigma_j^2$. We call

$$\begin{align*} p = \sum_{j = 1}^{N_o} p_j = \sum_{n = 1}^{N_o} \sigma_j^{-2} = \,\mathrm{tr}\left(\Sigma^{-1}\right). \end{align*}$$

the total precision of the sensor array. It can be seen that its reciprocal $\sigma^2_{\text{tot}} = 1/p$ corresponds to the harmonic average of the variances divided by the number of sensors N_o . Therefore, the misfit in (_β, ) satisfies

$$\begin{align*} \left\| K\mu - z^d\left(\varepsilon\right) \right\|^2_{\Sigma_0^{-1}} = \frac{1}{p}\sum_{j = 1}^{N_o} \sigma^{-2}_j \vert \lbrack K \mu \rbrack_j - z^d\left(\varepsilon\right)_j\vert^2. \end{align*}$$

For identical sensors and measurements $\varepsilon \sim \mathcal{N}(0, \operatorname{Id}_{N_o} )$ this simply leads to the scaled Euclidean norm $(1/N_o)\left\| K\mu - z^d(\varepsilon) \right\|^2_2$. In general, by increasing the total precision of the sensor setup p, we improve the measurements by proportionally decreasing the variances by $\sigma^2_{\text{tot}}$. While this will decrease the expected level of noise through its distribution, it will not affect the misfit functional, which is just influenced by Σ₀, or the normalized variances $\sigma^2_{0,j} = \sigma^2_j/\sigma^2_{\text{tot}}$.

Moreover, since $\varepsilon \sim \mathcal{N}(0,\Sigma)$, we have $\Sigma^{-1/2}\varepsilon \sim \mathcal{N}(0, _{N_o})$ and by direct calculations, the following estimate holds

$$\begin{align*} \frac{N_o}{\sqrt{N_o + 1}} \unicode{x2A7D} \mathbb{E}_{\gamma_p}\left[\left\| \varepsilon \right\|_{\Sigma^{-1}}\right] \unicode{x2A7D} \sqrt{N_o}. \end{align*}$$

Hence, with high probability, realizations of the error fulfill the estimate

$$\begin{align*} \sqrt{\sum_{j = 1}^{N_o} \varepsilon_j^2/\sigma^2_j} = \left\| \varepsilon \right\|_{\Sigma^{-1}} = \sqrt{p} \left\| \varepsilon \right\|_{\Sigma^{-1}_0} \unicode{x2A7D} C\sqrt{N_o} \end{align*}$$

and thus $\left\| \varepsilon \right\|_{\Sigma^{-1}_0} \lesssim 1/\sqrt{p}$. Thus, we consider the expected noise $\sigma_{\text{tot}} = 1/\sqrt{p}$ as an (expected) upper bound for the noise. This motivates the parameter choice rule

$$\begin{align} \beta\left(p\right) = \beta_0 / \sqrt{p} = \beta_0 \,\mathrm{tr} \left(\Sigma^{-1}\right)^{-1/2} \end{align}$$

for some $\beta_0\gt0$ large enough.

We are now prepared to prove a quantitative estimate on the worst-case mean-squared error by lifting the deterministic result of theorem 5.6 to the stochastic setting.

Theorem 6.1. Assume that $\eta^\dagger$ is θ-admissible for $\theta \in (0,1)$ and set $\beta(p) = \beta_0/\sqrt{p}$. Then there exists

$$\begin{align*} \overline{p} = \beta_0^2 c_{p}\left( \theta, k, \mu^{\dagger}, \left\| \mathcal{I}_0^{-1} \right\|_{W_{\dagger}^{-1} \to W_{\dagger}}\right) \end{align*}$$

such that for $p \unicode{x2A7E} \overline{p}$, there holds

$$\begin{align} \mathbb{E}_{\gamma_p}\left[\sup_{{\mu} \in \mathfrak{M}\left(\cdot\right)} d_{\mathrm{HK}}\left({\mu},\mu^{\dagger}\right)^2\right] \unicode{x2A7D} 8 \mathbb{E}_{\gamma_p}\left[\left\| \delta\widehat{\boldsymbol{m}} \right\|^2_{W_\dagger}\right] + C_3\exp \left[-\left(\dfrac{ \theta^2 \beta_0}{64 C_4}\right)^2/\left(2N_o\right)\right], \end{align} \tag{ 6.1 }$$

where $C_3 = 2\| \mu^\dagger \|_{\mathcal{M}(\Omega_s)} + \sqrt{2N_o}/(2\beta_0 \sqrt{p})$ and $C_4 = \max \{C_1,C_2 \}$. In addition, the expectation $\mathbb{E}_{\gamma_p}[\left\| \delta\widehat{\boldsymbol{m}} \right\|^2_{W_\dagger}]$ has the closed form

$$\begin{align} \mathbb{E}_{\gamma_p}\left[\left\| \delta\widehat{\boldsymbol{m}} \right\|^2_{W_\dagger}\right] = \frac{1}{p} \left(\mathrm{tr}\left(W_{\dagger}\mathcal{I}_0^{-1}\right) + \beta_0^2 \left\| \mathcal{I}_0^{-1} \left(\boldsymbol{\rho};\boldsymbol{0}\right) \right\|^2_{W_{\dagger}}\right). \end{align} \tag{ 6.2 }$$

Proof. Define the sets

$$\begin{align*} A_1 = \left\{\varepsilon : C_4 \beta\left(p\right)^{-1}\left\| \varepsilon \right\|_{\Sigma_0^{-1}} \unicode{x2A7D} \frac{\theta^2}{64}\right\}, \quad A_2 = \left\{\varepsilon : C_4 \beta\left(p\right)^{-1}\left(\left\| \varepsilon \right\|_{\Sigma^{-1}_0} + \beta\left(p\right)\right)^2 \unicode{x2A7D} \frac{\theta^2}{64}\right\}. \end{align*}$$

By a case distinction, we readily verify

$$\begin{align*} \mathbb{R}^{N_o} \setminus \left(A_1 \cap A_2\right) \subset \left(\mathbb{R}^{N_o}\setminus A_1\right) \cup \left(\left(\mathbb{R}^{N_o}\setminus A_2\right) \cap A_1\right) \end{align*}$$

and thus

$$\begin{align} \begin{aligned} \mathbb{E}_{\gamma_p}\left\lbrack \sup_{\mu \in \mathfrak{M}\left(\cdot\right)} d_{\mathrm{HK}}\left(\mu,\mu^{\dagger}\right)^2 \right \rbrack \unicode{x2A7D} & \int_{A_1 \cap A_2} \sup_{\mu \in \mathfrak{M}\left(\varepsilon\right)} d_{\mathrm{HK}}\left(\mu,\mu^{\dagger}\right)^2 \mathrm{d}\gamma_p\left(\varepsilon\right)\\ &+ \underbrace{\int_{\mathbb{R}^{N_o}\backslash A_1} \sup_{\mu \in \mathfrak{M}\left(\varepsilon\right)} d_{\mathrm{HK}}\left({\mu},\mu^{\dagger}\right)^2 \mathrm{d}\gamma_p\left(\varepsilon\right)}_{I_1} \\ &+ \underbrace{\int_{\left(\mathbb{R}^{N_o}\backslash A_2\right) \cap A_1} \sup_{\mu \in \mathfrak{M}\left(\varepsilon\right)} d_{\mathrm{HK}}\left(\mu,\mu^{\dagger}\right)^2 \mathrm{d}\gamma_p\left(\varepsilon\right)}_{I_2}. \end{aligned} \end{align} \tag{ 6.3 }$$

For $\varepsilon \in A_1 \cap A_2$, we have

$$\begin{align} C_2 \beta\left(p\right)^{-1} \left(\left(\left\| \varepsilon \right\|_{\Sigma^{-1}_0} + \beta\left(p\right)\right)^2 + \left\| \varepsilon \right\|_{\Sigma_0^{-1}}\right) &\unicode{x2A7D} C_4 \beta^{-1}\left(p\right)\left(\left\| \varepsilon \right\|_{\Sigma^{-1}_0} + \beta\left(p\right)\right)^2 + C_4 \beta^{-1}\left(p\right)\left\| \varepsilon \right\|_{\Sigma_0^{-1}} \nonumber\\ & \unicode{x2A7D} \frac{\theta^2}{64}+\frac{\theta^2}{64} = \frac{\theta^2}{32}, \end{align}$$

i.e. ε satisfies (5.10). Moreover, expanding the square in the definition of A₂, we conclude that (5.9) also holds due to

$$\begin{align*} 2C_4 \left(\left\| \varepsilon \right\|_{\Sigma^{-1}_0} + \beta\left(p\right)\right) \unicode{x2A7D} \frac{\theta^2}{32} \unicode{x2A7D} \frac{\sqrt{\theta}}{2 \sqrt{32}}. \end{align*}$$

Hence, for $\varepsilon \in A_1 \cap A_2$, there holds $\mathfrak{M}(\varepsilon) = \{\widehat{\mu}\}$ and

$$\begin{align} \sup_{\mu \in \mathfrak{M}\left(\varepsilon\right)} d_{\mathrm{HK}}\left(\mu,\mu^{\dagger}\right)^2 = d_{\mathrm{HK}}\left(\widehat{\mu},\mu^{\dagger}\right)^2 \unicode{x2A7D} 8 \left\| \delta{\widehat{\boldsymbol{m}}} \right\|^2_{W_\dagger} \end{align} \tag{ 6.4 }$$

by proposition 5.6. Next, we estimate I₁ by

$$\begin{align*} d_{\mathrm{HK}}\left(\widehat{\mu}, \mu^{\dagger}\right)^2 \unicode{x2A7D} \left\| \mu^\dagger \right\|_{\mathcal{M}\left(\Omega_s\right)} + \left\| \widehat{\mu} \right\|_{\mathcal{M}\left(\Omega_s\right)} \unicode{x2A7D} 2\left\| \mu^\dagger \right\|_{\mathcal{M}\left(\Omega_s\right)} + \left\| \varepsilon \right\|_{\Sigma_0^{-1}}^2/\left(2\beta_0/\sqrt{p}\right) \end{align*}$$

applying proposition 3.3 and [24, proposition 7.8]. Together with lemma A.1 this yields

$$\begin{align} \begin{aligned} I_1 = \int_{\mathbb{R}^{N_o}\backslash A_1}d_{\mathrm{HK}}\left( \widehat{\mu}, \mu^{\dagger}\right)^2 \mathrm{d}\gamma_p\left(\varepsilon\right) &\unicode{x2A7D} \int_{\mathbb{R}^{N_o}\backslash A_1} \left(2\left\| \mu^\dagger \right\|_{\mathcal{M}\left(\Omega_s\right)} + \left\| \varepsilon \right\|_{\Sigma^{-1}}^2/\left(2\beta_0\sqrt{p}\right)\right) \mathrm{d}\gamma_p\left(\varepsilon\right) \\ &\unicode{x2A7D} \left(2\left\| \mu^\dagger \right\|_{\mathcal{M}\left(\Omega_s\right)} + \dfrac{\sqrt{2N_o}}{2\beta_0 \sqrt{p}} \right) \exp \left[-\left(\dfrac{ \theta^2 \beta_0}{64 C_4 }\right)^2\Bigg/\left(2N_o\right) \right]. \end{aligned} \end{align} \tag{ 6.5 }$$

Finally, for $\varepsilon \in (\mathbb{R}^{N_o}\backslash A_2) \cap A_1$, one has

$$\begin{align*} p^{1/4} \left(\frac{\theta^2 \beta_0}{64C_4}\right)^{1/2} - \beta_0 \lt \left\| \varepsilon \right\|_{\Sigma^{-1}} = p^{-1/2} \left\| \varepsilon \right\|_{\Sigma^{-1}_0} \unicode{x2A7D} \dfrac{\theta^2 \beta_0}{64C_4} \end{align*}$$

where the first inequality follows from $\varepsilon \not \in A_2$ and the second follows from $\varepsilon \in A_1$. Hence, if we choose

$$\begin{align*} p \unicode{x2A7E} \beta_0^2\left(\frac{\theta^2}{64 C_4} + 1\right)^4 \Bigg/ \left(\frac{\theta^2 }{64 C_4}\right)^2 : = \beta_0^2 \overline c_{p}: = \overline{p} \end{align*}$$

then $(\mathbb{R}^{N_o}\backslash A_2) \cap A_1$ is empty and I₂ = 0. Together with (6.3)–(6.5), we obtain (6.1) for every $p \unicode{x2A7E} \overline{p}$. The equality in (6.2) follows immediately from the closed form expression (5.6) for $\delta \widehat{\boldsymbol{m}}$ and $\varepsilon \sim \mathcal{N}(0, p^{-1}\Sigma_0)$. □

Let us interpret this result: By choosing β₀ large enough, the second term on the right hand side of (6.6) becomes negligible, i.e.

$$\begin{align} \mathbb{E}_{\gamma_p}\left[\sup_{{\mu} \in \mathfrak{M}\left(\cdot\right)} d_{\mathrm{HK}}\left({\mu},\mu^{\dagger}\right)^2\right] \unicode{x2A7D} 8 \mathbb{E}_{\gamma_p}\left[\left\| \delta\widehat{\boldsymbol{m}} \right\|^2_{W_\dagger}\right]+ \delta \end{align} \tag{ 6.6 }$$

for some $0 \lt\delta \ll 1$. As a consequence, due to its closed form representation (6.2), $\mathbb{E}_{\gamma_p}[\left\| \delta\widehat{\boldsymbol{m}} \right\|^2_{W_\dagger}]$ provides a computationally inexpensive, approximate upper surrogate for the worst-case mean-squared error which vanishes as $p \rightarrow \infty$. Moreover, due to its explicit dependence on the measurement setup, it represents a suitable candidate for an optimal design criterion in the context of optimal sensor placement for the class of sparse inverse problems under consideration. This potential will be further investigated in a follow-up paper.

Remark 6.2. It is worth mentioning that the constant 8 appearing on the right hand side of (6.6) is not optimal and is primarily a result of the proof technique. In fact, by appropriately selecting constants in propositions B.2 and 5.2, it is possible to replace 8 by $1 + \delta$, where $0 \lt \delta \ll 1$ at the cost of increasing $\bar{p}$. We will illustrate the sharpness of the estimate of the worst-case mean-squared error by $\mathbb{E}_{\gamma_p}[\left\| \delta\widehat{\boldsymbol{m}} \right\|^2_{W_\dagger}]$ in the subsequent numerical results.

Remark 6.3. Relying on similar arguments as in the proof of theorem 6.1, we are also able to derive pointwise estimates on the Hellinger–Kantorovich distance which hold with high probability. Indeed, noticing that (6.4) holds in the set $A_1 \cap A_2$, we derive a lower probability bound for $\mathbb{P}( \varepsilon \in A_1 \cap A_2)$ by noticing that

$$\begin{align*} \mathbb{P}\left( \varepsilon \in A_1 \cap A_2\right) &\unicode{x2A7E} \mathbb{P}\left( \varepsilon \in A_1\right) + \mathbb{P}\left( \varepsilon \in A_2\right) - 1 \\ &\unicode{x2A7E} 1 - \mathbb{P}\left(\varepsilon \in \mathbb{R}^{N_o}\backslash A_1\right) - \mathbb{P}\left(\varepsilon \in \mathbb{R}^{N_o}\backslash A_2\right). \end{align*}$$

By invoking lemma A.1, one has

$$\begin{align*} &\mathbb{P}\left(\varepsilon \in \mathbb{R}^{N_o}\backslash A_1\right) = \mathbb{P} \left(\left\| \varepsilon \right\|_{\Sigma^{-1}} > \dfrac{\theta^2 \beta_0}{64 C_4}\right) \unicode{x2A7D} 2\exp\left[- \left(\dfrac{\theta^2 \beta_0}{64 C_4}\right)^2\Bigg/\left(2 N_o\right) \right],\\ &\mathbb{P}\left(\varepsilon \in \mathbb{R}^{N_o}\backslash A_2\right) = \mathbb{P} \left( \left\| \varepsilon \right\|_{\Sigma^{-1}} > \dfrac{p^{1/2} \theta \beta_0^{1/2} }{8 C_4^{1/2}} - \beta_0 \right) \unicode{x2A7D} 2 \exp\left[\!- \left(\dfrac{p^{1/2} \theta \beta_0^{1/2} }{8 C_4^{1/2}} - \beta_0\right)^2\Bigg/\left(2 N_o\right) \!\right]. \end{align*}$$

Hence, since $\exp(-x^2) \to 0$ as $x \to \infty$, we can see that for every $\delta \in (0,1)$, one can choose β₀ and p large enough such that

$$\begin{align*} \exp\left[- \left(\dfrac{\theta^2 \beta_0}{64 C_4}\right)^2\Bigg/\left(2 N_o\right) \right] \lt \delta/4,\quad \exp\left[- \left(\dfrac{p^{1/2} \theta \beta_0^{1/2} }{8 C_4^{1/2}} - \beta_0\right)^2\Bigg/\left(2 N_o\right) \right] \lt \delta/ 4, \end{align*}$$

which implies $\mathbb{P}(\varepsilon \in A_1 \cap A_2) \unicode{x2A7E} 1 - \delta$. Therefore, with probability at least $1-\delta$, we have

$$\begin{align} \sup_{{\mu} \in \mathfrak{M}\left(\varepsilon\right)} d_{\mathrm{HK}}\left({\mu},\mu^{\dagger}\right)^2 \unicode{x2A7D} 8\left\| \delta \widehat{\boldsymbol{m}} \right\|^2_{W_{\dagger}} \end{align}$$

for realization ε of the noise. Furthermore, employing lemma A.1 again, we know that with probability at least $1- \delta$, and independently from p, one has $\left\| \varepsilon \right\|_{\Sigma^{-1}} \unicode{x2A7D} \sqrt{-2N_o\ln(\delta/2)}$. Hence, by proposition 5.2 together with $\varepsilon \in A_1 \cap A_2$, we have

$$\begin{align*} \sup_{{\mu} \in \mathfrak{M}\left(\varepsilon\right)} d_{\mathrm{HK}}\left({\mu},\mu^{\dagger}\right)^2 &\unicode{x2A7D} 8C_1 p^{-1/2}\left(\left\| \varepsilon \right\|_{\Sigma^{-1}} + \beta_0\right) \\ &\unicode{x2A7D} 8C_1 p^{-1/2} \left(\sqrt{-2N_o\ln\left(\delta/2\right)} + \beta_0 \right) \end{align*}$$

with probability at least $1-2 \delta$.

In the first example, we illustrate the reconstruction capabilities of the proposed ansatz for different measurement setups and with and without noise in the observations. To this end, we attempt to recover the reference measure $\mu^{\dagger}$ using a variable number N_o of uniformly distributed sensors. For noisy data, the regularization parameter is selected as $\beta = \beta_0 / \sqrt{p}$ where β₀ = 2 and $p = 10^4$. We first consider the exact measurement data with $N_o \in \{6,9,11\}$ and try to obtain $\mu^\dagger$ by solving (₀ ). The results are shown in figure 1. We observe that with 6 sensors, the pre-certificate $\eta_{\mathrm{PC}}$ is not admissible. Recalling [11, proposition 7], this implies that $\mu^\dagger$ is not a minimum norm solution. In contrast, the experiments with 9 and 11 uniform sensors provide admissible pre-certificates. In these situations, the pre-certificates coincide with the minimum norm dual certificates and the ground truth $\mu^\dagger$ is indeed an identifiable minimum norm solution.

Zoom In Zoom Out Reset image size

Next, we consider noisy data and solve (_β, ) for the aforementioned choice of $\beta(p)$. Following the observation in the first example, we only evaluate the reconstruction results obtained by 9 and 11 uniform sensors. In the absence of the measurement data obtained from experiments, we generate synthetic noisy measurements where the noise vector ε is a realization of the Gaussian random noise $\varepsilon \sim \mathcal{N}(0,\Sigma)$. The results are shown in figure 2. Since $\mu^\dagger$ is identifiable in these cases, $\widehat{\mu}$ and $\bar{\mu}$ coincide and closely approximate $\mu^\dagger$ with high probability for an appropriate choice of β₀ and p large enough. Both properties can be clearly observed in the plots, where β₀ = 2.

Zoom In Zoom Out Reset image size

In the second example we study the influence of the parameter choice rule on the reconstruction result. To this end, we fix the measurement setup to 9 uniformly distributed sensors. We recall that the a priori parameter choice rule is given by $\beta(p) = \beta_0/\sqrt{p}$. According to section 6.2, selecting a sufficiently large value for β₀ is recommended to achieve a high quality reconstruction. To determine a useful range of regularization parameters, we solve problem (_β, ) for a sequence of regularization parameters using PDAP. Here, we choose $\beta_0 \in \{0.5, 1, 2\}$ and $p \in \{10^4, 10^5, 10^6\}$.

In figure 3, different reconstruction results are shown for the same realization of noise, $\beta_0 \in \{0.5, 1, 2\}$ and $p = 10^{4}$. As one can see, for this particular realization of the noise, the number of spikes is recovered exactly in the case β₀ = 2 and we again observe that $\widehat{\mu} = \bar{\mu}$. In contrast, for smaller β₀, the noisy pre-certificate is not admissible. Hence, while $\widehat{\mu}$ still provides a good approximation of $\mu^\dagger$, $\bar{\mu}$ admits two additional spikes away from the support of $\mu^\dagger$. These observations can be explained by looking at theorem 6.1 the second term on the right hand side of the inequality becomes negligible for increasing β₀ and large enough p. Thus, roughly speaking, the parameter β₀ controls the probability of the 'good events' in which $\widehat{\mu}$ is the unique solution of (_β, ).

Zoom In Zoom Out Reset image size

Finally, we address the reconstruction error from a quantitative perspective. For this purpose, we simplify the evaluation of the maximum mean-squared error (MSE) by inserting the solution $\bar{\mu}$ computed algorithmically. We note that this could only lead to an under-estimation of the maximum error in the case of non-unique solutions of (_β, ); a degenerate case that is unlikely to occur in practice. Moreover, the expectation is approximated using 10³ Monte-Carlo samples. Additionally, we use the closed form expression (6.2) for evaluating the linearized estimate $\mathbb{E}_{\gamma_p}[\left\| \delta \widehat{\boldsymbol{m}} \right\|_{W_{\dagger}}^2]$ exactly. Here, the expectations are computed for $\beta_0 \in\{2, 0.5\}$. The results are collected in table 1. We make several observations: Clearly, the MSE decreases for increasing p, i.e. lower noise level. For increased β₀, the behavior differs: For the theoretical quantities $\widehat{\boldsymbol{m}}$ and $\delta\widehat{\boldsymbol{m}}$ increased β₀ only introduces additional bias and thus increases error. For the estimator $\bar{\mu}$, the increased regularization however leads to generally improved results, since the probability of $\widehat{\mu}\neq \bar{\mu}$ is decreased. We highlight in bold the estimator which performed best for each β₀. Here, the results conform to theorem 6.1: For larger β₀, the second term on the right-hand side of (6.1) is negligible and the linearized estimate provides an excellent bound on the MSE for both $\widehat\mu$ and $\bar\mu$. We also note that the estimate is closer to the MSE in the limiting case for larger p. In contrast, for β = 0.5, the linearized estimate and the MSE of $\widehat\mu$ are much smaller than the MSE of the estimator $\bar{\mu}$. This underlines the observation that theorem 5.6 requires further restrictions on the admissible noises in comparison to proposition 5.2.

Table 1. Reconstruction results with β₀ = 2 and $\beta_0 = 0.5$.

		$p = 10^4$	$p = 10^5$	$p = 10^6$	$p = 10^7$
β0 = 2	$\mathbb{E}_{\gamma_p}[\left\| \delta \widehat{\boldsymbol{m}} \right\|^2_{W_{\dagger}}]$	$7.97 \cdot 10^{-3}$	$7.97 \cdot 10^{-4}$	$7.97 \cdot 10^{-5}$	$7.97 \cdot 10^{-6}$
	$\mathbb{E}_{\gamma_p}[d_{\mathrm{HK}}(\mu^{\dagger},\widehat{\mu})^2]$	$5.43 \cdot 10^{-3}$	$6.49 \cdot 10^{-4}$	$7.35 \cdot 10^{-5}$	$7.76 \cdot 10^{-6}$
	$\mathbb{E}_{\gamma_p}[d_{\mathrm{HK}}(\mu^{\dagger},\bar{\mu})^2]$	$5.44 \cdot 10^{-3}$	$\boldsymbol{6.51 \cdot 10^{-4}}$	$\boldsymbol{7.42 \cdot 10^{-5}}$	$\boldsymbol{7.99 \cdot 10^{-6}}$
$\beta_0 = 0.5$	$\mathbb{E}_{\gamma_p}[\left\| \delta \widehat{\boldsymbol{m}} \right\|^2_{W_{\dagger}}]$	$2.07 \cdot 10^{-3}$	$2.07 \cdot 10^{-4}$	$2.07 \cdot 10^{-5}$	$2.07 \cdot 10^{-6}$
	$\mathbb{E}_{\gamma_p}[d_{\mathrm{HK}}(\mu^{\dagger},\widehat{\mu})^2]$	$1.71 \cdot 10^{-3}$	$1.94 \cdot 10^{-4}$	$2.03 \cdot 10^{-5}$	$2.06 \cdot 10^{-6}$
	$\mathbb{E}_{\gamma_p}[d_{\mathrm{HK}}(\mu^{\dagger},\bar{\mu})^2]$	$\boldsymbol{4.12 \cdot 10^{-3}}$	$9.83 \cdot 10^{-4}$	$2.65 \cdot 10^{-4}$	$7.94 \cdot 10^{-5}$

The final example is devoted to compare the reconstruction results obtained by uniform designs and an improved design chosen by heuristics. To this end, we consider three measurement setups: uniformly distributed setups with 6 and 11 sensors, respectively, and one with 6 sensors selected on purpose. More precisely, in the later case, we place the sensors at $\Omega_o = \{-0.8, -0.6, -0.4, -0.1, 0.1, 0.4 \}$. The different error measures are computed as in the previous example and the results are gathered in table 2. We observe that the measurement setup with 6 selected sensors performs better than the uniform ones. Moreover, the linearized estimate again provides a sharp upper bound on the error for both ten uniform and six selected sensors but yields numerically singular Fisher information matrices for six uniform sensors (denoted as Inf in the table), i.e. $\mu^\dagger$ is not stably identifiable in this case. Note that the estimator $\bar{\mu}$ still yields somewhat useful results, which are however affected by a constant error due to the difference in minimum norm solution and exact source as depicted in figure 1 and do not improve with lower noise level. These results suggest that the reconstruction quality does not only rely on the amount of measurements taken but also on their specific setup. In this case, we point out that the selected sensors are chosen to be adapted to the sources as every two sensors are placed on the two sides of every source. Thus the obtained results imply that if we have some reasonable prior information on the source positions and amplitudes, one may obtain a better sensor placement setup by incorporating it in the design of the measurement setup. This leads to the concept of optimal sensor placement problems for sparse inversion which we will consider in a future work.

Table 2. Reconstruction results with different sensor setups.

		11 sensors	'selected' 6 sensors	6 sensors
$\mathbb{E}_{\gamma_p}[d_{\mathrm{HK}}(\mu^{\dagger},\bar{\mu})^2]$	$p = 10^4$	$5.03 \cdot 10^{-3}$	$\boldsymbol{4.25 \cdot 10^{-3}}$	$1.54 \cdot 10^{-2}$
	$p = 10^5$	$5.61 \cdot 10^{-4}$	$\boldsymbol{4.58 \cdot 10^{-4}}$	$1.19 \cdot 10^{-2}$
	$p = 10^6$	$6.31 \cdot 10^{-5}$	$\boldsymbol{4.65 \cdot 10^{-5}}$	$1.18 \cdot 10^{-2}$
$\mathbb{E}_{\gamma_p}[\left\| \delta \widehat{\boldsymbol{m}} \right\|^2_{W_{\dagger}}]$	$p = 10^4$	$6.09 \cdot 10^{-3}$	$4.77 \cdot 10^{-3}$	Inf
	$p = 10^5$	$6.09 \cdot 10^{-4}$	$4.77 \cdot 10^{-4}$
	$p = 10^6$	$6.09 \cdot 10^{-5}$	$4.77 \cdot 10^{-5}$

The following lemma gives an estimate on the tail probabilities of a Gaussian random variables as well as on its moments.

Lemma A.1. Assume that $\varepsilon_0 \sim \gamma_E = \mathcal{N}(0, \operatorname{Id}_{N_o})$. Then for every α > 0 there holds

$$\begin{align*} \mathbb{P}\left(\left\| \varepsilon_0 \right\|_2 > \alpha\right) \unicode{x2A7D} 2 \exp \left( -\frac{\alpha^2}{2N_o}\right). \end{align*}$$

Moreover for $l\unicode{x2A7E} 1$, with $C_l = (2l-1)!! = (2l-1)(2l-3)\cdots 1$, we have

$$\begin{align*} \int_{\left\| \varepsilon_0 \right\|_2 > \alpha} \left\| \varepsilon_0 \right\|^l_2 \mathrm{d}\gamma_E\left(\varepsilon_0\right) \unicode{x2A7D} \sqrt{2 N_o C_l} \exp \left( -\frac{\alpha^2}{4N_o}\right) = \exp \left( -\frac{\alpha^2}{4N_o} + \log\left(2 N_o C_l\right)/2\right) . \end{align*}$$

Proof. According to remark 4 in [32], we get for any λ > 0 that

$$\begin{align*} \mathbb{P}\left(\left\| \varepsilon_0 \right\|_2 > \alpha\right) = \int_{\left\| \varepsilon_0 \right\| > \alpha} \mathrm{d}\gamma_E\left(\varepsilon_0\right) \unicode{x2A7D} 2 \exp\left(-\lambda \alpha\right) \exp\left(\lambda^2 N_o/2 \right). \end{align*}$$

Minimizing the right-hand side with respect to λ yields $\lambda = \alpha/N_o$ and the first estimate. The second inequality is due to

$$\begin{align*} \int_{\left\| \varepsilon_0 \right\| > \alpha} \left\| \varepsilon_0 \right\|^l_2 \mathrm{d}\gamma_E\left(\varepsilon_0\right) \unicode{x2A7D} \sqrt{\int \left\| \varepsilon_0 \right\|^{2l}_2 \mathrm{d}\gamma_E\left(\varepsilon_0\right)} \sqrt{2 \exp \left( -\frac{\alpha^2}{2N_o}\right)} \end{align*}$$

with Cauchy-Schwarz. The proof is finished noting that $\mathbb{E}[ \left\| \varepsilon_0 \right\|^{2l}_2 ] = N_o C_l$ where $C_l = (2l-1)!!$ denotes the 2l-the moment of the univariate standard normal distribution. □

In this section we address the measurability of the worst-case distance function

$$\begin{align} \varepsilon \mapsto \sup_{\mu \in\mathfrak{M}\left(\varepsilon\right)} d_{\mathrm{HK}}\left(\mu ,\mu^\dagger\right) \end{align} \tag{ A.1 }$$

as well as the boundedness of its second moment. For this purpose, recall the definition of the solution set

$$\begin{align*} \mathfrak{M}\left(\varepsilon\right) = \mathop{\mathrm{arg}\,\mathrm{min}}_{\mu\in\mathcal{M}\left(\Omega_s\right)} \left \lbrack \dfrac{1}{2} \left\| K \mu - z^d\left(\varepsilon\right) \right\|_{\Sigma_0^{-1}}^2 + \beta\left\| \mu \right\|_{\mathcal{M}\left(\Omega_s\right)} \right\rbrack \end{align*}$$

and note that $\mathfrak{M}(\varepsilon)$ is weak* compact. We first show that the supremum in (A.1) is attained and give a useful upper bound on it.

Lemma A.2. For every $\varepsilon \in \mathbb{R}^{N_o}$, there is $\bar{\mu}(\varepsilon)\in \mathfrak{M}(\varepsilon)$ with

$$\begin{align*} d_{\mathrm{HK}}\left(\bar{\mu}\left(\varepsilon\right),\mu^\dagger\right) = \sup_{\mu \in \mathfrak{M}\left(\varepsilon\right)} d_{\mathrm{HK}}\left(\mu,\mu^\dagger\right). \end{align*}$$

Moreover, we have

$$\begin{align} \sup_{\mu \in \mathfrak{M}\left(\varepsilon\right)} d_{\mathrm{HK}}\left(\mu,\mu^\dagger\right)^2 \unicode{x2A7D} 2\left\| \mu^\dagger \right\|_{\mathcal{M}\left(\Omega_s\right)} + \frac{1}{2\beta} \left\| \varepsilon \right\|_{\Sigma_0^{-1}}^2. \end{align} \tag{ A.2 }$$

Proof. The inequality (A.2) follows from proposition 3.3 and [24, proposition 7.18]. Hence, there exists a supremizing sequence $\{\mu_k\}_k \subset \mathfrak{M}(\varepsilon)$, i.e,

$$\begin{align*} \lim_{k\rightarrow \infty} d_{\mathrm{HK}}\left(\mu_k,\mu^\dagger\right) = \sup_{\mu \in \mathfrak{M}\left(\varepsilon\right)} d_{\mathrm{HK}}\left(\mu,\mu^\dagger\right). \end{align*}$$

Due to weak* compactness of $\mathfrak{M}(\varepsilon)$, it admits a subsequence, denoted by the same subscript, as well as $\bar{\mu}(\varepsilon) \in \mathfrak{M}(\varepsilon)$ with $\mu_k \rightharpoonup^* \bar{\mu}(\varepsilon)$. Since $d_{\mathrm{HK}}$ metrizes weak* convergence, the inverse triangle inequality yields

$$\begin{align*} |d_{\mathrm{HK}}\left(\mu_k,\mu^\dagger\right)-d_{\mathrm{HK}}\left(\bar{\mu}\left(\varepsilon\right),\mu^\dagger\right)| \unicode{x2A7D} d_{\mathrm{HK}}\left(\bar{\mu}\left(\varepsilon\right),\mu_k\right) \rightarrow 0 \end{align*}$$

and thus

$$\begin{align*} d_{\mathrm{HK}}\left(\bar{\mu}\left(\varepsilon\right),\mu^\dagger\right) = \lim_{k\rightarrow \infty} d_{\mathrm{HK}}\left(\mu_k,\mu^\dagger\right) = \sup_{\mu \in \mathfrak{M}\left(\varepsilon\right)} d_{\mathrm{HK}}\left(\mu,\mu^\dagger\right). \end{align*}$$

 □

Using lemma A.2, we conclude the measurability of the worst-case distance.

Proposition A.3. The function defined in (A.1) is γ_p -measurable. Moreover, there holds

$$\begin{align*} \mathbb{E}_p \left \lbrack \sup_{\mu \in \mathfrak{M}\left(\cdot\right)} d_{\mathrm{HK}}\left(\mu,\mu^\dagger \right)^2\right \rbrack \unicode{x2A7D} 2\left\| \mu^\dagger \right\|_{\mathcal{M}\left(\Omega_s\right)} + \frac{N_0}{2\beta p} < \infty \end{align*}$$

Proof. For abbreviation, define

$$\begin{align*} \operatorname{Err}\left[\bar{\mu}\right]\left(\varepsilon\right) = \sup_{\mu \in \mathfrak{M}\left(\varepsilon\right)} d_{\mathrm{HK}}\left(\mu,\mu^\dagger\right) \end{align*}$$

and let $\{\varepsilon_k\}_k$ denote a convergent sequence with limit ɛ. Note that the set $\bigcup_k \mathfrak{M}(\varepsilon_k)$ as well as the sequence $\operatorname{Err}[\bar{\mu}](\varepsilon_k)$ are bounded. Now, choose a subsequence $\{\varepsilon_{k,i}\}_i$ such that

$$\begin{align*} \lim_{i \rightarrow \infty} \operatorname{Err}\left[\bar{\mu}\right]\left(\varepsilon_{k,i}\right) = \limsup_{k \rightarrow \infty} \operatorname{Err}\left[\bar{\mu}\right]\left(\varepsilon_{k}\right) \end{align*}$$

and let $\{\bar{\mu}(\varepsilon_{k,i})\}_i$ be a corresponding sequence of maximizers from lemma A.2. By possibly extracting another subsequence, there is $\tilde{\mu}$ with $\bar{\mu}(\varepsilon_{k,i}) \rightharpoonup^* \tilde{\mu}$. By standard arguments, we verify that $\tilde{\mu} \in \mathfrak{M}(\varepsilon)$. Consequently, we have

$$\begin{align*} \limsup_{k \rightarrow \infty} \operatorname{Err}\left[\bar{\mu}\right]\left(\varepsilon_{k}\right) = \lim_{i \rightarrow \infty} \operatorname{Err}\left[\bar{\mu}\right]\left(\varepsilon_{k,i}\right) = \lim_{i \rightarrow \infty} d_{\mathrm{HK}}\left(\mu\left(\varepsilon_{k,i}\right), \mu^\dagger\right)^2 = d_{\mathrm{HK}}\left(\tilde{\mu}, \mu^\dagger\right)^2 \unicode{x2A7D} \operatorname{Err}\left[\bar{\mu}\right]\left(\varepsilon\right). \end{align*}$$

Hence $\operatorname{Err}[\bar{\mu}]$ is upper semicontinuous and thus measurable w.r.t γ_p . Finally, we apply (A.2) to conclude

$$\begin{align*} \mathbb{E}_p \left \lbrack \sup_{\mu \in \mathfrak{M}\left(\cdot\right)} d_{\text{HK}}\left(\mu,\mu^\dagger \right)^2\right \rbrack \unicode{x2A7D} 2\left\| \mu^\dagger \right\|_{\mathcal{M}\left(\Omega_s\right)} + \frac{1}{2\beta} \mathbb{E}_p \left \lbrack \left\| \varepsilon \right\|_{\Sigma_0^{-1}}^2 \right \rbrack \unicode{x2A7D} 2\left\| \mu^\dagger \right\|_{\mathcal{M}\left(\Omega_s\right)} + \frac{N_0}{2\beta p}. \end{align*}$$

 □

The distance induced by the Kantorovich–Rubinstein (KR) norm (equivalent to the 'Bounded-Lipschitz' norm) is also referred to as the 'flat' metric, and metricises weak* convergence on bounded sets in $\mathcal{M}(\Omega_s)$. The norm can be defined by

$$\begin{align*} \left\| \mu \right\|_{\mathrm{KR}} : = \sup \left\{\int_{\Omega_s} f \mathrm{d}\mu: f \in C^{0,1}\left(\Omega_s\right), \left\| f \right\|_{C\left(\Omega_s\right)} \unicode{x2A7D} 1 \text{ and } \operatorname{Lip} \left(f\right) \unicode{x2A7D} 1 \right\}. \end{align*}$$

Here, $C^{0,1}(\Omega_s)$ is the space of Lipschitz functions on $\Omega_s$ and

$$\begin{align*} \operatorname{Lip}\left(f\right) : = \sup_{x \neq y} \dfrac{\left\| f\left(x\right) - f\left(y\right) \right\|}{\left\| x-y \right\|}, \quad x, y \in \Omega_s. \end{align*}$$

The KR distance is then set to

$$\begin{align*} d_{\mathrm{KR}}\left(\mu_1,\mu_2\right) = \left\| \mu_1 - \mu_2 \right\|_{\mathrm{KR}}. \end{align*}$$

Although its evaluation for general sparse measures requires the solution of a minimization problem (see [21]) it has many useful properties. For positive measures $\mu_1,\mu_2 \unicode{x2A7E} 0$ it is equivalent to a generalized Wasserstein-1 distance for measures not necessarily of the same TV-norm as in [28];

$$\begin{align} d_{\mathrm{KR}}\left(\mu_1,\mu_2\right) = \inf_{\left\{\,\tilde{\mu}_1,\tilde{\mu}_2 \colon \; \left\| \tilde{\mu}_1 \right\|_{\mathcal{M}} = \left\| \tilde{\mu}_2 \right\|_{\mathcal{M}} \,\right\}} \left[ \left\| \tilde{\mu}_1 - \mu_1 \right\|_{\mathcal{M}} + \left\| \tilde{\mu}_2 - \mu_2 \right\|_{\mathcal{M}} + W_1\left(\tilde{\mu}_1,\tilde{\mu}_2\right) \right]. \end{align} \tag{ C.1 }$$

For signed measures, we use the Jordan decomposition $\mu_i = \mu^+_i - \mu^-_i$ and observe that

$$\begin{align*} d_{\mathrm{KR}}\left(\mu_1, \mu_2\right) = \left\| \left(\mu_1^+ + \mu_2^-\right) - \left(\mu_2^+ + \mu_1^-\right) \right\|_{\mathrm{KR}} = d_{\mathrm{KR}} \left(\mu_1^+ + \mu_2^-, \mu_2^+ + \mu_1^-\right), \end{align*}$$

which then allows to apply the characterization (C.1). This representation, together with the help of [28, proposition 2] allows to characterize the KR distance of single point sources with weight of equal sign:

$$\begin{align} d_{\mathrm{KR}}\left(q_1\delta_{y_1},q_2\delta_{y_2}\right) & = \min_{0 \unicode{x2A7D} \theta \mathrm{sign}{q_1} \unicode{x2A7D} \min \left\{|q_1|, |q_2| \right\}}\left[ |\theta - q_1| + |\theta - q_2| + \theta \left\| y_1 - y_2 \right\|_2 \right] \nonumber\\ & = |q_1 - q_2| + \min \left\{|q_1|, |q_2| \right\} \min\left\{\left\| y_1 - y_2 \right\|_2, 2\right\}. \end{align} \tag{ C.2 }$$

For the case of $\mathrm{sign} q_1 \neq \mathrm{sign} q_2$, we instead have $d_{\mathrm{KR}}(q_1\delta_{y_1},q_2\delta_{y_2}) = |q_1| + |q_2|$. The above formula can be used for all finitely supported measures with the same number of support points by applying the triangle inequality. Motivated by property (C.2), the Kantorovich-Rubinstein distance $d_{\mathrm{KR}}$ is also appropriate to quantify the distance between two discrete measures, and a similar upper bound as in proposition 4.2 can be obtained, i.e.

$$\begin{align*} d_{\mathrm{KR}}\left(\mu,\mu^\dagger\right) \unicode{x2A7D} \sum_{n = 1}^{N}\left( \vert q_n - q_n^\dagger\vert + \vert q_n^\dagger\vert \, \left\| y_n - y_n^\dagger \right\|_2 \right) \unicode{x2A7D} 2\sqrt{2 \left\| \mu^\dagger \right\|} \left\| \boldsymbol{m} - \boldsymbol{m}^\dagger \right\|_{W_\dagger}, \end{align*}$$

However, this bound either uses an $\ell_1$ like sum over the weighted errors, which is not suitable for the rest of our analysis, or the equivalence between a weighted $\ell_1$ and $\ell_2$ norm (by Hölder's inequality), and is not an asymptotically sharp bound.