# Bayesian Adaptive Eigenspace Inversion

Marie Graff, Bamdad Hosseini, Alfred Kim

► To cite this version:

Marie Graff, Bamdad Hosseini, Alfred Kim. Bayesian Adaptive Eigenspace Inversion. 2025. hal-05130206

HAL Id: hal-05130206

https://hal.science/hal-05130206v1

Preprint submitted on 26 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Bayesian Adaptive Eigenspace Inversion

Marie Graff¹,²,*, Bamdad Hosseini³ & Alfred Kim²

¹ Te Whai Ao — The Dodd-Walls Centre for Photonic and Quantum Technologies, New Zealand.

² Department of Mathematics, The University of Auckland, Auckland 1010, New Zealand.

³ Department of Applied Mathematics, University of Washington, Seattle WA, USA.

* Corresponding author: marie.graff@auckland.ac.nz

June 26, 2025

## Abstract

We develop a probabilistic analogue of the Adaptive Eigenspace Inversion (AEI) method, an algorithm for the adapting regularisation terms towards the recovery of sharp interfaces, from the perspective of Bayesian inverse problems to enable the quantification of uncertainties in the reconstructions. We show that the adaptations of the regularisation terms in AEI is analogous to adapting the prior distributions in a sequence of Bayesian inverse problems. We present a theoretical result, showing that appropriate convergence of the AEI loop will lead to convergence of the resulting Bayesian posteriors in a Wasserstein distance. Furthermore, extensive numerical experiments demonstrate the performance and wide applicability of the Bayesian extension of AEI in both linear and nonlinear inverse problems. In particular, we observe in our experiments that the posterior means recover sharp edges in the unknown fields, and the posterior variances are largely concentrated around sharp interfaces, a feature that is also observed in other models for recovery of piecewise constant fields.

## 1 Introduction

In recent years, uncertainty quantification (UQ) has been an integral part of scientific computing and engineering workflows in particular as it pertains to inverse problems; see [1] and references therein. UQ allows us to provide more than just an estimate of unknown parameters, but also add statistics regarding the reconstructed parameters that encode our confidence in the reconstructions. A workhorse of UQ for inverse problems is the Bayesian framework [2, 3, 4], which combines our prior knowledge of the unknown parameter with a likelihood term involving measurements to obtain a posterior distribution by combining the two. Bayesian Inverse Problems (BIPs) are now widely used in natural sciences and engineering, including image processing and computer vision [5], medical imaging [6, 7, 8], and geophysical prospecting [9, 10, 11], among others.

In most BIPs, the likelihood is assumed to be known due to well-known partial differential equation (PDE) models for physical processes of interest along with accepted models for the measured data. However, the choice of the prior distribution is less trivial and often involves a balancing of expert physical knowledge with mathematical and computational convenience. Gaussian priors, which are analogous to Tikhonov-type regularisation terms in deterministic inverse problems, are by far the most widely used class of priors [4, 2]. However, these priors are widely known to be inadequate for sparse or edge preserving recovery [12, 13, 14]. To
resolve this, various sparsity promoting priors have been proposed [12, 15, 16, 17, 18, 19] that can, for example, be used to impose sparsity in a wavelet basis towards the recovery of sharp images.

As is often the case for non-Gaussian priors, optimisation methods are used to compute a maximum a posteriori (MAP) estimator [9, 10, 20, 21, 22] as a pointwise estimator of the unknown parameter. To extract further statistics, a more complete computation of the posterior distribution is necessary. This can be achieved either by approximation, e.g., with the Laplace's approximation of the covariance [23, 24], or by using Markov chain Monte Carlo (MCMC) [25, 26, 27, 28]. However, the design of MCMC algorithms in the function space is largely limited to Gaussian priors [29, 30, 31] with consistent algorithms for non-Gaussian priors remaining extremely scarce [17, 32].

In this paper, we introduce an alternative approach for the construction of edge preserving priors motivated by the Adaptive Eigenspace Inversion (AEI) method [33, 34, 35, 36, 37, 38] and other ideas such as iterative thresholding [39] and nonlinear diffusion [40, 41]. Recently, AEI has reached growing interest from the theoretical perspective [42, 43, 44], as well as applications [45, 46, 47, 48, 49, 50, 51, 52].

Our original contribution in this paper is to propose an iterative algorithm using a Gaussian prior where, at each iteration, an estimator such as the posterior mean or MAP is used to update the prior covariance operator through an eigenvalue problem, as is done in the AEI algorithm. More specifically, our contributions are fourfold: (1) We introduce the Bayesian AEI method, a probabilistic formulation of the AEI procedure where an adaptive Gaussian prior is employed; (2) We prove a convergence result for the sequence of AEI posteriors, showing that if the eigenfunctions computed in the AEI loop converge appropriately then the resulting posteriors will also converge in the Wasserstein-1 distance; (3) We use the Bayesian formulation to suggest a modification of the AEI procedure using an adaptive penalty term; And (4) we present a number of numerical experiments for linear and nonlinear inverse problems in scattering where the Bayesian AEI procedure is applied and tested using variational inference techniques as well as MCMC. Most notably, our numerical results show that the posterior means obtained by AEI preserve edges, and posterior variances tend to concentrate around discontinuities, a behavior that does not happen under standard smoothness Gaussian priors.

The structure of the paper is as follows: In Section 2, we first recall the AEI method in its original version and then present the Bayesian extension by identifying an adaptive Gaussian prior distribution. In Section 3, we prove our convergence result for the posteriors. Finally, we summarise our numerical experiments in Section 4.

## 2 Methodology

Let $\Omega$ be an open subset of $\mathbb{R}^n$, $n = 2, 3$, and $\Gamma$ a subset of $\Omega$. Let us now define $\mathcal{M}$ a Hilbert space of functions on $\Omega$, referred to as the *model space*, and $\mathcal{D}$ a Hilbert space of functions on $\Gamma$, corresponding to the *data space*. Let $G$ denote the compact continuous forward operator mapping (linear or not) from $\mathcal{M}$ to $\mathcal{D}$. We aim to consider an imaging problem for a forward model $G$ that consists in the reconstruction of a parameter $\mathbf{m}_{\text{exact}} \in \mathcal{M}$ from given measurements $\mathbf{d} \in \mathcal{D}$. The resulting problem can be rewritten as the minimisation of the distance between the acquired data $\mathbf{d}$ and the predictions $G(\hat{\mathbf{m}})$, that is,

*Given the noisy data $\mathbf{d} \in \mathcal{D}$, find a parameter $\hat{\mathbf{m}} \in \mathcal{M}$ such that*

$$ \hat{\mathbf{m}} = \underset{\mathbf{m} \in \mathcal{M}}{\operatorname{argmin}} \frac{1}{2} \|\mathbf{d} - G(\mathbf{m})\|_{\mathcal{D}}^{2}. \qquad (1) $$
Such a problem is usually ill-posed in the sense of Hadamard and extra precautions are required to ensure existence, uniqueness, and stability of the solution. We now assume that existence and uniqueness is solved by referring to the pseudoinverse theory [53] and its nonlinear extensions [54]. In the following section, we recall the AEI methodology towards regularisation of $\hat{\mathbf{m}}$ [33, 34, 35, 37].

## 2.1 The AEI Method

The principle of AEI is to represent the parameter $\mathbf{m}$ in a basis that is adapted in a recursive manner. In particular, at iteration $k \ge 1$ we opt for the eigen basis obtained by solving the elliptic eigenvalue problem

$$ \left\{ \begin{array}{ll} -\nabla \cdot (\mu[\mathbf{m}^{k-1}]\nabla) \varphi_j^{k-1} = \lambda_j^{k-1} \varphi_j^{k-1}, & \text{in } \Omega, \\ \varphi_j^{k-1} = 0, & \text{on } \partial\Omega, \end{array} \right. \qquad (2) $$

where $\varphi_j^{k-1}$ are the eigenfunctions and $\lambda_j^{k-1}$ are the eigenvalues, and $\mu$ is given by

$$ \mu[\mathbf{m}^{k-1}] = \frac{1}{\sqrt{|\nabla \mathbf{m}^{k-1}|^2 + \varepsilon^2}}, \quad \varepsilon > 0. \qquad (3) $$

Here $\mathbf{m}^{k-1}$ is our estimate of the unknown parameter from the previous iteration. When $k = 1$ we simply take $\mu \equiv 1$, which gives the Dirichlet eigenvalue problem. The small parameter $\varepsilon$ is used for stability and ensures that the denominator never vanishes. We note that (2) and (3) are reminiscent of the celebrated model of Perona and Malik [40, 41, 55].

The choice of Dirichlet boundary conditions in (2) is innocuous. In fact, we can consider the problem

$$ \left\{ \begin{array}{ll} -\nabla \cdot (\mu[\mathbf{m}^{k-1}]\nabla) \mathbf{m}_0^{k-1} = 0 & \text{in } \Omega, \\ \mathbf{m}_0^{k-1} = \mathbf{m}_{\text{exact}} & \text{on } \partial\Omega, \end{array} \right. \qquad (4) $$

where we assume that the recovered parameter is known on the boundary $\partial\Omega$ and $\mathbf{m}_0^{k-1}$ is called the background. In this representation, our aim is to compute the next estimator $\mathbf{m}^k$ in the form

$$ \mathbf{m}^k := \mathbf{m}_0^{k-1} + \sum_{j=1}^{K} \beta_j \varphi_j^{k-1}, \qquad (5) $$

with $K$ the number of eigenfunctions taken in the basis and $\beta_j$ the coefficients of $\mathbf{m}^k$ in the eigen basis of (2). Note that in practice, the basis will be truncated for two reasons: (i) it is too expensive to compute all the eigenfunctions as there are as many of them as mesh gridpoints, (ii) the regularisation is induced by truncation, hence parameter $K$ acts as a regularisation parameter. In the following, $\Phi_k$ will represent the concatenation of the first $K$ eigenfunctions [$\varphi_1^k, \dots, \varphi_K^k$] and $\beta \in \mathbb{C}^K$ will denote the vector containing the coefficients $\beta_j$ so that we can introduce the shorthand notation

$$ \mathbf{m}^k = \mathbf{m}_0^{k-1} + \Phi_{k-1} \beta. $$

With these notations, the AEI method now reads:

Given the noisy data $\mathbf{d} \in \mathcal{D}$ and an index $k_{\max} \ge 1$, find a parameter $\hat{\mathbf{m}} \in \mathcal{M}$ such that

$$ \hat{\mathbf{m}} = \mathbf{m}_0^{k_{\max}-1} + \Phi_{k_{\max}-1} \hat{\beta}, \quad \text{and} \quad \hat{\beta} = \underset{\beta \in \mathbb{C}^K}{\operatorname{argmin}} \frac{1}{2} \|\mathbf{d} - G(\Phi_{k_{\max}-1} \beta)\|_{\mathcal{D}}^2. \qquad (6) $$
This procedure can be summarised in Algorithm 1 where for a known level of noise $\delta$ on the data $\mathbf{d}$, the number of adaptations $k_{\max}$ to perform can be determined using the Morozov discrepancy principle, that is, the adaptations stop when the misfit reaches the noise level to avoid overfitting.

<caption>**Algorithm 1 Adaptive Eigenspace Inversion**</caption>

1: **Input:** $\mathbf{d}, G, K, \delta$
2: **Output:** $\hat{\mathbf{m}}$
3: Set $\mu[\hat{\mathbf{m}}] = 1$
4: Solve the BVP for $\mathbf{m}_0^0$
5: Solve the eigenvalue problem for $\Phi_0 = [\varphi_1^0, \dots, \varphi_K^0]$
6: Set $k = 1$ and $\hat{\mathbf{m}} \leftarrow \mathbf{m}_0$
7: **while** $\|\mathbf{d} - G(\hat{\mathbf{m}})\|_D > \delta$ **do**
8: &nbsp;&nbsp;&nbsp;&nbsp;Solve the least-squares problem to get $\hat{\beta}$
9: &nbsp;&nbsp;&nbsp;&nbsp;Update $\mathbf{m}^k \leftarrow \mathbf{m}_0^{k-1} + \Phi_{k-1} \hat{\beta}$
10: Solve the PDE for $\mathbf{m}_0^k$, with $\mu[\mathbf{m}^k]$
11: Solve the eigenvalue problem for $\Phi_k = [\varphi_1^k, \dots, \varphi_K^k]$, with $\mu[\mathbf{m}^k]$
12: **end while**
13: Return $\hat{\mathbf{m}} = \mathbf{m}^k$

In Algorithm 1, Line 8 remains purposely vague as any minimisation algorithm can be implemented depending on the problem at hand. For example, simple least-squares solvers for linear inverse problems, or Gauss-Newton or similar quasi-Newton methods if the inverse problem is nonlinear. The AEI method acts as an implicit regularisation scheme through the choice of $K$, the total number of basis functions, as well as the adaptation procedure for modifying the matrix $\Phi^k$ at each iteration. Finally, we note that various discretisation methods including finite differences or finite elements can be used to solve the subsequent eigenvalue problems (2), see [35] for example.

## 2.2 The Bayesian AEI Method

We now outline the Bayesian analogue of the AEI algorithm recalled above. We recall that Bayes' rule identifies the posterior distribution of the unknown parameter $\mathbf{m}$ as the product of a likelihood distribution and a prior [2]:

$$ \nu(\mathbf{m}) \propto \pi(\mathbf{d} | \mathbf{m})\eta(\mathbf{m}), \qquad (7) $$

where $\nu \in \mathbb{P}(M)$ denotes the posterior distribution, $\pi(\mathbf{d} | \mathbf{m})$ is the likelihood function, and $\eta \in \mathbb{P}(M)$ is the prior. Assuming $D \subseteq \mathbb{R}^m$ and an additive Gaussian noise model for $\mathbf{d}$, i.e.,

$$ \mathbf{d} := G(\mathbf{m}) + \epsilon, \quad \epsilon \sim \mathcal{N}(0, \delta^2 I), \qquad (8) $$

we obtain the familiar quadratic likelihood

$$ \pi(\mathbf{d}|\mathbf{m}) \propto \exp \left( -\frac{1}{2\delta^2} \|\mathbf{d} - G(\mathbf{m})\|_D^2 \right), $$

which corresponds to the least squares problem outlined in (1) and (6). Henceforth, we consider this quadratic likelihood model in mind but note that the rest of exposition naturally includes other types of likelihoods as well.

Our proposed Bayesian AEI method considers a sequence of posteriors defined via (7), in particular

$$ \nu^k(\mathbf{m}) \propto \pi(\mathbf{d} | \mathbf{m})\eta^k(\mathbf{m}), \qquad (9) $$
where the priors $\eta^k = \mathcal{N}(\mathbf{m}_0^{k-1}, C_k)$ are Gaussian distributions with covariance operators of the form

$$C_k := \Phi_{k-1} \Lambda_{k-1}^{-2s} \Phi_{k-1}^{\top},$$

where we emphasize that $s \ge 0$ is an exponent and not an index. Here $\Phi_{k-1}$ denotes the unitary matrices associated to the eigenfunctions $\varphi_j^{k-1}$ and $\Lambda_{k-1}$ is the diagonal matrix of the eigenvalues $\lambda_j^{k-1}$ as in (2).

Observing that the matrices $\Phi_{k-1}$ are computed using a previous pointwise estimator $\mathbf{m}^{k-1}$ in the original AEI framework we propose to follow a similar idea in the Bayesian approach where we take $\mathbf{m}^k$ to be either the MAP estimator of $\nu^k$,

$$\mathbf{m}_{\text{MAP}}^{k} = \mathbf{m}_{0}^{k-1} + \Phi_{k-1} \beta_{\text{MAP}}^{k}, \quad \beta_{\text{MAP}}^{k} = \underset{\beta \in \mathbb{C}^{K}}{\operatorname{argmin}} \frac{1}{2\sigma^2} \|\mathbf{d} - G(\mathbf{m}_{0}^{k-1} + \Phi_{k-1}\beta)\|_{D}^{2} + \frac{1}{2} \|\Lambda_{k-1}^{s}\beta\|_{2}^{2}. \quad (10)$$

or its conditional mean

$$\mathbf{m}_{\text{COM}}^{k} = \int \mathbf{m} \, d\nu^{k}(\mathbf{m}).$$

The MAP estimator is analogous to the standard AEI algorithm and can be used in settings where an optimization algorithm for solving (10) is available. The conditional mean estimator is desirable when sampling algorithms such as MCMC are used for interrogation of $\nu^k$. In our numerical experiments in Section 4 we will use both of these estimators.

To this end, the goal of Bayesian AEI can be summarised as follows:

*Given the noisy data $\mathbf{d} \in \mathcal{D}$, provide statistics on the posterior*

$$\nu^{k_{\text{max}}}(\mathbf{m}) \propto \pi(\mathbf{d} | \mathbf{m})\eta^{k_{\text{max}}}(\mathbf{m}), \quad \eta^{k_{\text{max}}} = \mathcal{N}(\mathbf{m}_0^{K_{\text{max}}-1}, \Phi_{k_{\text{max}}-1} \Lambda_{k_{\text{max}}-1}^{-2s} \Phi_{k_{\text{max}}-1}^{\top}).$$

# 3 Theory

In this section, we outline some theoretical results in support of the Bayesian AEI method. Our main result states that the sequence of AEI posterior distributions will converge, with respect to a Wasserstein distance, so long as the matrices $\Lambda_k, \Phi_k$ converge in an appropriate sense. For simplicity, here we take $\mathbf{m}_0^k \equiv 0$ in our analysis so that we can simply take our priors $\eta^k = \mathcal{N}(0, C_k)$, where we recall $C_k = \Phi_k \Lambda_k^{-2s} \Phi_k^\top$. Note that the likelihood distribution is fixed in the iterates (9). Then our goal is to show that, for an appropriate probability metric $d$, it holds that $d(\nu^k, \nu^\infty) \to 0$ with $\nu^\infty$ denoting the limiting posterior assuming that the pairs $(\Lambda_k, \Phi_k) \to (\Lambda_\infty, \Phi_\infty)$ in a sense that will be made clear shortly.

Following the perturbation theory of [56, 57], we take $d$ to be a Wasserstein metric as it enables the quantification of posterior perturbations with respect to prior perturbations. For an overview of Wasserstein distances we refer the reader to [58]. Consider a Hilbert space $H$ and let $\mathbb{P}(H)$ denote the space of Borel probability measures on $H$. Further define the space

$$\mathbb{P}^p(H) := \{\rho \in \mathbb{P}(H) \mid \int \|\mathbf{m}\|_H^p \rho(d\mathbf{m}) < +\infty\}.$$

Then the Wasserstein-$p$ metric for $p \in [1, +\infty]$, denoted as $W_p$, on $\mathbb{P}^p(H)$ is defined as follows:

$$W_p(\rho_1, \rho_2) = \left( \inf_{\gamma \in \Pi(\rho_1, \rho_2)} \int_{H \times H} \|\mathbf{m} - \mathbf{m}'\|_H^p \gamma(d\mathbf{m}, d\mathbf{m}') \right)^{1/p}, \quad \rho_1, \rho_2 \in \mathbb{P}^p(H). \quad (11)$$
Here $\Pi(\mu_1, \mu_2) := \{\gamma \in \mathbb{P}(H \times H) \mid \int_H \gamma(\cdot, dm') = \rho_1$ and $\int_H \gamma(d\mathbf{m}, \cdot) = \rho_2\}$ denotes the space of all couplings between $\rho_1, \rho_2$, i.e., probability measures on $H \times H$ where the first marginal coincides with $\rho_1$ while the second marginal coincides with $\rho_2$.

For the analysis of Bayesian AEI we will take $H = L^2(\Omega)^1$ for a Lipschitz domain $\Omega \subset \mathbb{R}^n$ and consider likelihoods of the form

$$ \pi(\mathbf{m} | \mathbf{d}) = \frac{1}{2\delta^2} \|G(\mathbf{m}) - \mathbf{d}\|_2^2, $$

where $G: L^2(\Omega) \to \mathbb{R}^m$ is a bounded and linear operator, i.e., $G(\mathbf{m}) \le C_G \|\mathbf{m}\|_{L^2(\Omega)}$ for some constant $C_G \ge 0$ $^2$. In this setting we have the following theorem that allows us to control the $W_1$ distance between the AEI posteriors $\nu^k$ and $\nu^\infty$ with the $W_2$ distance between the priors $\eta^k$ and $\eta^\infty$.

**Theorem 1.** Consider the above setting for Bayesian AEI with a bounded and linear forward map $G$ and a sequence of Gaussian priors $\eta^k$ and a presumed limiting Gaussian prior $\eta^\infty$. Then there exists a constant $C = C(\mathbf{d}, \delta, C_G) > 0$ so that, whenever $W_2^2(\eta^k, \eta^\infty) < \min\{1, \int \exp(-\frac{1}{\delta}C_G \|\mathbf{m}\|_{L^2(\Omega)}^2) \eta^\infty(d\mathbf{m})\}$, then

$$ W_1(\nu^k, \nu^\infty) \le C \frac{1 + \int \|\mathbf{m}\|_{L^2(\Omega)}^2 \eta^\infty(d\mathbf{m})}{\left(\int \exp(-\frac{1}{\delta}C_G \|\mathbf{m}\|_{L^2(\Omega)}^2) \eta^\infty(d\mathbf{m})\right)^2} W_2(\eta^k, \eta^\infty). $$

*Proof.* Following [56, Thm. 3.8] as well as the calculations in Sec. 4.1 and 4.2 of that paper, we have the bound,

$$ W_1(\nu^k, \nu^\infty) \le C \frac{\left(1 + \int \|\mathbf{m}\|_{L^2(\Omega)}^2 \eta^k(d\mathbf{m}) + \int \|\mathbf{m}\|_{L^2(\Omega)}^2 \eta^\infty(d\mathbf{m})\right)^{1/2}}{\left(\int \exp\left(-\frac{1}{\delta}C_G \|\mathbf{m}\|_{L^2(\Omega)}^2\right) \eta^k(d\mathbf{m}) \cdot \left(\int \exp\left(-\frac{1}{\delta}C_G \|\mathbf{m}\|_{L^2(\Omega)}^2\right) \eta^\infty(d\mathbf{m})\right)\right)^2} W_2(\eta^k, \eta^\infty) \quad (12) $$

where $C(\mathbf{d}, \delta, C_G) > 0$ is a constant. Now let $\gamma^* \in \Pi(\eta^k, \eta^\infty)$ be the optimal coupling achieving $W_2(\eta^k, \eta^\infty)$. Then we can write

$$ \begin{aligned} \int \|\mathbf{m}\|_{L^2(\Omega)}^2 \eta^k(d\mathbf{m}) &= \int \|\mathbf{m}\|_{L^2(\Omega)}^2 \gamma^*(d\mathbf{m}, d\mathbf{m}') \\ &\le 2 \int \|\mathbf{m} - \mathbf{m}'\|_{L^2(\Omega)}^2 \gamma^*(d\mathbf{m}, d\mathbf{m}') + 2 \int \|\mathbf{m}'\|_{L^2(\Omega)}^2 \gamma(d\mathbf{m}') \\ &= 2 \int \|\mathbf{m} - \mathbf{m}'\|_{L^2(\Omega)}^2 \gamma^*(d\mathbf{m}, d\mathbf{m}') + 2 \int \|\mathbf{m}\|_{L^2(\Omega)}^2 \eta^\infty(d\mathbf{m}) \\ &= 2W_2^2(\eta^k, \eta^\infty) + 2 \int \|\mathbf{m}\|_{L^2(\Omega)}^2 \eta^\infty(d\mathbf{m}). \end{aligned} $$

Moreover, using the strong duality for $W_2$ and the fact that $\exp(-\|\mathbf{m}\|_{L^2(\Omega)}^2)$ is bounded and globally Lipschitz, we infer, using a similar calculation to the above, that

$$ \int \exp \left( -\frac{1}{\delta} C_G \|\mathbf{m}\|_{L^2(\Omega)}^2 \right) \eta^k(d\mathbf{m}) \le C \left( W_2^2(\eta^k, \eta^\infty) + \int \exp \left( -\frac{1}{\delta} C_G \|\mathbf{m}\|_{L^2(\Omega)}^2 \right) \eta^\infty(d\mathbf{m}) \right) $$

<footnote n="1">¹One can also choose a Sobolev space here to impose further regularity on the parameter space as needed.</footnote>

<footnote n="2">²Our analysis can be generalized to nonlinear $G$ under sufficient local Lipschitz continuity assumptions using the framework of [56]. However, we work with the linear forward map for simplicity as it lends itself to the use of common Wasserstein distances rather than those defined with more complicated cost functions.</footnote>
for some constant $C > 0$ that is independent of $\eta^k$ and $\eta^\infty$. Substituting these bounds into (12) we obtain the inequality

$$W_1(\nu^k, \nu^\infty) \le C \frac{\left(1 + W_2^2(\eta^k, \eta^\infty) + \int \|\mathbf{m}\|_{L^2(\Omega)}^2 \eta^\infty(d\mathbf{m})\right)^{1/2}}{\left(W_2^2(\eta^k, \eta^\infty) + \int \exp\left(-\frac{1}{\delta}C_G \|\mathbf{m}\|_{L^2(\Omega)}^2\right) \eta^\infty(d\mathbf{m})\right) \cdot \left(\int \exp\left(-\frac{1}{\delta}C_G \|\mathbf{m}\|_{L^2(\Omega)}^2\right) \eta^\infty(d\mathbf{m})\right)} \cdot W_2(\eta^k, \eta^\infty).$$

Provided that $W_2^2(\eta^k, \eta^\infty) \le \min\{1, \int \exp\left(-\frac{1}{\delta}C_G \|\mathbf{m}\|_{L^2(\Omega)}^2\right) \eta^\infty(d\mathbf{m})\}$ the bound then simplifies to

$$W_1(\nu^k, \nu^\infty) \le C \frac{\left(1 + \int \|\mathbf{m}\|_{L^2(\Omega)}^2 \eta^\infty(d\mathbf{m})\right)^{1/2}}{\left(\int \exp\left(-\frac{2}{\delta}C_G \|\mathbf{m}\|_{L^2(\Omega)}^2\right) \eta^\infty(d\mathbf{m})\right)^2} W_2(\eta^k, \eta^\infty),$$

for a constant $C > 0$ that is independent of $\eta^k$ and $\eta^\infty$ as desired. $\square$

This bound suggests that the convergence of the posterior measures $\nu^k$ to $\nu^\infty$ is controlled by the rate at which the priors $\eta^k$ converge to $\eta^\infty$. The following theorem gives a quantitative rate for this error.

**Proposition 1.** Let $\eta^k = \mathcal{N}(0, C_k)$ be the sequence of AEI priors defined as above with $\eta^\infty = \mathcal{N}(0, C_\infty)$ denoting the limiting prior. Then it holds that

$$W_2^2(\eta^k, \eta^\infty) = \mathrm{Tr} C_k + \mathrm{Tr} C_\infty - 2\mathrm{Tr}(C_k^{1/2}C_\infty C_k^{1/2})^{1/2} =: d_{\mathrm{WB}}(C_k, C_\infty), \quad (13)$$

where Tr C denotes the trace of a Hermitian operator. We further have the explicit upper bound

$$W_2^2(\eta^k, \eta^\infty) \le 2 \left( \sum_{j=1}^{\infty} (\lambda_j^\infty)^{-s} \|\varphi_j^k - \varphi_j^\infty\|_{L^2(\Omega)}^2 + \sum_{j=1}^{\infty} |(\lambda_j^\infty)^{-s} - (\lambda_j^k)^{-s}|^2 \right). \quad (14)$$

*Proof.* The first display (13) is a well-known expression for the $W_2$-distance between Gaussian measures on Hilbert spaces [59, 60]. The upper bound (14) can be shown directly using a coupling argument as in [56, Lem. 4.1]. $\square$

The quantity $d_{\mathrm{WB}}$ defined above is known as the Wasserstein-Bures (WB) distance between the Hermitian operators. The above theorem presents an avenue for analysis of Bayesian AEI by connecting it to the deterministic analysis of the AEI method. For example, if the operators $C_k$, produced by Bayesian AEI, converge to a limit $C_\infty$ in the WB metric, then we readily obtain convergence of the $\eta^k$ to $\eta^\infty$ in $W_2$ which in turn yields the convergence of the posterior measures.

**Corollary 1.** Consider the above setting and suppose that Bayesian AEI produces a sequence of covariance operators $C_k$ converging to $C_\infty$ with respect to $d_{\mathrm{WB}}$. Then $W_1(\nu^k, \nu^\infty) \to 0$.

# 4 Numerical experiments

In this section, we illustrate the performance and UQ capabilities of the Bayesian AEI method. We start with a linear problem to observe essentially the AEI reconstruction, using linear least squares and a closed form for the resulting posteriors [2]. Then, we extend to a nonlinear problem and sample the posteriors using MCMC [9]. These examples show the computational performance and interesting properties of Bayesian AEI for recovering unknown piecewise constant fields with sharp edges.
## 4.1 Example 1: A linear inverse problem

For this first example, we consider the Helmholtz equation with a constant wavenumber $k$, a source term $f$ and homogeneous Neumann boundary conditions in the open set $\Omega = (0, 1) \times (0, 1)$ of $\mathbb{R}^2$

$$ \left\{ \begin{array}{ll} -k^2 u - \Delta u = f, & \text{in } \Omega, \\ \partial_n u = 0, & \text{on } \partial\Omega, \end{array} \right. \qquad (15) $$

where $\partial_n$ is the outer normal derivative operator on $\partial\Omega$. Using finite differences or finite element methods, $\Omega$ is discretised as a mesh with $N$ gridpoints and Problem (15) can be reformulated as

$$ H[k]\mathbf{u} = \mathbf{f}, $$

where $H[k]$ represents the discretisation of the Helmholtz operator by an $N \times N$ matrix, for a given $k$ including Neumann conditions.

For the inverse problem, we assume that the wavenumber $k$ is known, but we are seeking for the source term $f$: here the source problem is linear. We also opt for the discretise-then-optimise strategy and decide to select our data as random points in the discretised domain. A projection matrix $P: \mathbb{R}^N \to \mathbb{R}^d$ restricts the data vector to the values selected randomly.

### 4.1.1 A Toy example

<figure><img src="image_2.png" /><figcaption>Figure 1: Example 1 – Exact source (a), solution of the forward problem solved for $k = 100$ (b) and measurements with 1000 randomly chosen data point with 5% of noise (c).</figcaption></figure>

Figure 1 displays three panels. On the left, we show the source distribution $\mathbf{f}$ composed by three shapes of different size and amplitude used for this toy example. In the middle, we display the wavefield $\mathbf{u}$, solution to (15) with $k = 100$ and the source term on the left. The meshgrid has $N = 16641$ points. On the right, we highlight the location and value of the 1000 measurements $\mathbf{d}$, randomly selected, to which 5% of white noise has been added. From the data displayed on the right panel, we aim to reconstruct the source profile of the left panel using the AEI method, and even quantify the uncertainty on this reconstruction.

For the sake of comparison, we solve this inverse source problem using standard Tikhonov regularisations as well:

$$ \hat{\mathbf{f}}_{\text{Tikh},n} = \underset{\mathbf{f} \in \mathbb{R}^N}{\operatorname{argmin}} \frac{1}{2} \|\mathbf{d} - P H[k]^{-1} \mathbf{f}\|_2^2 + \frac{\alpha_{\text{Tikh}}}{2} \|L^n \mathbf{f}\|_2^2, \quad n = 0, 1, $$
where matrix $L$ is the discretisation of the gradient operator to correspond to the $H^1$-semi-norm, so that $L^0 = I$ the identity matrix for the $L^2$-norm equivalent. The choice of the regularisation parameter $\alpha_{\text{Tikh}}$ is determined by the Discrepancy Principle of Morozov [61]. Recall that, similarly, the number of adaptations performed by the AEI method is also determined by the same criterion. As a result, the misfits of all three reconstructions are about the level of noise.

<figure><img src="image_3.png" /><figcaption>Figure 2: Example 1 – Reconstruction of the source term (a) using three different regularization techniques: (b) Bayesian AEI, (c) $H^1$-Tikhonov and (d) $L^2$-Tikhonov. Same colour scale for all four images.</figcaption></figure>

In Figure 2, we show the source reconstructions obtained by the three regularisation techniques, compared to the true source term. With 5% white noise, the criterion value $\delta$ is equal to $5.8127 \times 10^{-4}$. For the reconstructions involving Tikhonov regularisation, we chose regularisation parameters $\alpha_{\text{Tikh}}$ such that the discrepancy principle of Morozov is satisfied: $\alpha_{\text{Tikh}} = 3.4514 \times 10^{-5}$ for the $L^2$ regularisation, leading to a misfit of $5.8281 \times 10^{-4}$, while $\alpha_{\text{Tikh}} = 1.1940 \times 10^{-8}$ for the $H^1$ regularisation and a misfit of $5.8162 \times 10^{-4}$. By comparison, our AEI algorithm converged after 12 adaptations and gave a misfit of $5.8212 \times 10^{-4}$, for $\alpha = 10^{-7}$ fixed at each iteration, $K = 150$ eigenfunctions for the first adaptation, then $K = 160$ fixed for the eleven following adaptations before convergence. For the sake of clarity, we do not display here the AEI solution at each adaptation: a panel is available in Appendix, Figure 14. At comparable misfits, all three constructions show all three features of the source term, but the AEI reconstruction outweighs the other two.

<figure><img src="image_4.png" /><figcaption>Figure 3: Example 1 – Uncertainty quantification of the AEI reconstruction: (a) exact source term; (b) mean or MAP estimate; (c) pointwise variance of the reconstruction.</figcaption></figure>

Finally, we compute the covariance matrix in the basis $\Phi_{11}$, which is then transformed back to a pointwise covariance in the mesh. We extract and display the resulting pointwise variance in Figure 3 (c), where (a) is the exact source and (b) is the MAP estimate. It is interesting to see that the variance is more informative than its counterpart for the Tikhonov regularisation. Clearly, the errors of reconstruction is located on the boundaries of the different features, while the reconstruction is certain on the flat areas.
4.1.2 Shepp-Logan phantom

To challenge our method, we used the modified Shepp-Logan phantom [62], scaled times
two, so that the background is zero and the skull ellipse is 2. For this example, we opted for
a increasing number of eigenfunctions at each adaptation, except of the very first adaptation
with Laplacian eigenfunctions. In Figure 4, we present the resulting reconstruction using
Bayesian AEI together with the L²- and H¹-Tikhonov reconstructions. All three methods
satisfy the discrepancy principle of Morozov, which means that for the Tikhonov methods,
we selected the parameter αTikh (2.5410 × 10⁻⁵ and 5.3827 × 10⁻⁹, respectively) so that the
misfits (5.1859 × 10⁻⁴ and 5.1908 × 10⁻⁴, respectively) are about the level of noise. Clearly,
the Bayesian AEI method is outperforming the other two methods in the reconstruction of
the source term.

<figure><img src="image_5.png" /><figcaption>Figure 4: Shepp-Logan phantom 1 – Reconstruction comparisons with the same display scale: (a) exact source term, (b) using Bayesian AEI, (c) using L²-Tikhonov regularisation with α = 2.541 × 10⁻⁵, and (d) using H¹-Tikhonov regularisation with α = 5.3827 × 10⁻⁹. Colour scale applies to all four images.</figcaption></figure>

To show the value of uncertainty quantification, we also computed the pointwise variance
matrix that highlights the regions in the source term where the uncertainty is high; see
Figure 5. As in the previous example, the borders of the different regions are highlighted
mostly, but also the regions that are actually not featured in the exact source term. This is
a good indication that these parts may not be correct in the reconstruction.

<figure><img src="image_6.png" /><figcaption>Figure 5: Shepp-Logan phantom 1 – Uncertainty quantification of the Bayesian AEI reconstruction: (a) exact source term for comparison, (b) MAP estimate, and (c) pointwise variance matrix.</figcaption></figure>

Finally, we show the evolution of the reconstruction through Bayesian AEI in Appendix,
Figure 15. We started with 150 Laplacian eigenfunctions, then restarting with only 50
eigenfunctions for the next adaptation, increasing this number by 50 at each subsequent
adaptation. In the 21st and final adaptation, the number of eigenfunctions is 1000. Note
that the panel only displays adaptations 1 to 20, as the final one is already available in
Figures 4 and 5.
## 4.2 Example 2: An inverse scattering problem

Now we consider the time-dependent wave equation with wave speed c, a source term f, and absorbing boundary conditions in the same set $\Omega = (0, 1) \times (0, 1)$ of $\mathbb{R}^2$:

$$ egin{cases} 
\partial_t^2 u - \nabla \cdot (c^2 \nabla u) = f, & \text{in } (0, T) \times \Omega, \\
u(x, 0) = 0, & \text{in } \Omega, t = 0, \\
\partial_t u(x, 0) = 0, & \text{in } \Omega, t = 0, \\
\partial_t u + c \partial_n u = 0, & \text{on } (0, T) \times \partial \Omega,
\end{cases} 
\quad (16) $$

with $c^2 = c_0^2 + p$, where $c_0$ is constant. To find the solution $u = G(\mathbf{c})$ of the forward problem, we discretise the variational formulation in space with finite elements, and in time with finite differences. In particular, we use the Newmark scheme in time.

The inverse scattering problem is nonlinear; we assume that the source term f is known and look for the wave speed c. The contrast p is represented on its AEI basis with coefficients $\beta$. We can now express the inverse scattering problem as

Given the noisy data $\mathbf{d} \in \partial\Omega$, find a parameter $\hat{\mathbf{c}} \in \mathbb{R}^N$ such that

$$ \hat{\mathbf{c}}^2 = \mathbf{c}_0^2 + \Phi\beta, \quad \text{and} \quad \beta = \underset{\beta \in \mathbb{R}^N}{\operatorname{argmin}} \frac{1}{2} \|\mathbf{d} - G(\mathbf{c}_0^2 + \Phi\beta)\|_2^2 + \frac{\alpha}{2} \|\Lambda^s \beta\|_2^2, \quad (17) $$

where we assume $c_0 = 1$ and choose $s = 1/2$. We opt for the optimise-then-discretise strategy, using the Lagrangian formulation for optimisation and the finite element method for discretisation. Following the reduced-space approach, the minimisation problem simplifies to finding the zero gradient of the Lagrangian with respect to c [63]. To do so, we must also have the Lagrange multiplier v governed by the adjoint equation:

$$ egin{cases} 
\partial_t^2 v - \nabla \cdot (c^2 \nabla v) = f, & \text{in } (0, T) \times \Omega, \\
v(x, T) = 0, & \text{in } \Omega, t = T, \\
\partial_t v(x, T) = 0, & \text{in } \Omega, t = T, \\
-\partial_t v + c \partial_n v = d - u_{|_{\partial\Omega}}, & \text{on } (0, T) \times \partial\Omega,
\end{cases} 
\quad (18) $$

in which the solution is available because the wave equation is time reversible. As a result, the gradient of the cost function to minimise can be well expressed in terms of the state u and the adjoint v. In Algorithm 1, we choose the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method [64] to solve our reduced minimisation problem with an approximated Hessian.

### 4.2.1 Another Toy example

Figure 6 shows the toy scatterer we would like to recover and the total wave, that is, the incident field plus the scattered response. We set a single source term $f(x, t) = g(x)r(t)$ centred at $x_l = (0.8, 0.2)$ of $\Omega$, where g is a Gaussian point source in space and r a Ricker wavelet in time:

$$ g(x) = \kappa \exp^{-\frac{(x-x_l)^2}{s}}, \quad x \in \Omega, \\ r(t) = (1 - 2\pi^2(\nu t - 1)^2) \exp^{-\pi^2(\nu t - 1)^2}, \quad t \in [0, T], $$

with central frequency $\nu = 10$, $s = 10^{-2}$ and $\kappa = 200$. We decide to have the data $\mathbf{d}$ at all points of the discretised boundary for each discrete time value and add 5% white noise. The meshgrid is identical to Example 1 with $N = 16641$ points.
<figure><img src="image_7.png" /><figcaption>Figure 6: Toy example 2 – Exact wave speed (a) with a maximum value of 7.8 and solution of the forward problem at progressing values of time: (b), (c) and (d) show three different snapshots as time increases.</figcaption></figure>

Again for the sake of comparison, we solve this inverse scattering problem with other regularisation methods as well, this time with Tikhonov and total variation (TV) regularisation:

$$ \hat{\mathbf{c}}_{\text{Tikh}} = \underset{\mathbf{c} \in \mathbb{R}^N}{\operatorname{argmin}} \frac{1}{2} \| \mathbf{d} - G(\mathbf{c})_{|\partial\Omega} \|_2^2 + \frac{\alpha_{\text{Tikh}}}{2} \|Lp\|_2^2, $$

$$ \hat{\mathbf{c}}_{\text{TV}} = \underset{\mathbf{c} \in \mathbb{R}^N}{\operatorname{argmin}} \frac{1}{2} \| \mathbf{d} - G(\mathbf{c})_{|\partial\Omega} \|_2^2 + \frac{\alpha_{\text{TV}}}{2} \| \sqrt{|\nabla p|^2 + \epsilon^2} \|_1, $$

where matrix L is the discretisation of the gradient operator to correspond to the H¹ semi-norm. We choose regularisation parameters α_Tikh and α-TV that satisfy the Discrepancy Principle of Morozov [61].

<figure><img src="image_8.png" /><figcaption>Figure 7: Toy example 2 – Reconstruction comparisons with the same display scale: (a) exact scatterer for comparison, (b) using Bayesian AEI, (c) using H¹-Tikhonov regularisation with α = 10⁻⁹, and (d) using TV regularisation with α = 10⁻¹².</figcaption></figure>

In Figure 7, we show the scatterer reconstructions obtained by the three regularisation techniques, compared to the true scatterer term. With 5% white noise, the criterion value δ
is equal to 5.2946. In terms of both contrast and finer details, Bayesian AEI performs significantly better than both Tikhonov and TV regularisation. The selected evolution of the reconstructions is shown in Figure 8. The AEI algorithm converged in 8 adaptations.

<figure><img src="image_9.png" /><figcaption>Figure 8: Toy example 2 – The second (a), fourth (b), sixth (c) and eighth (d) adaptation to recover the scatterer. The misfit decays from 5.6738 to 5.2946, reaching the level of noise. The number of eigenfunctions is incremented by 20 each adaptation, and $\alpha$ is gradually decreased in the latter half of adaptations from $10^{-4}$ to $10^{-6}$.</figcaption></figure>

Now we would like to quantify uncertainties in a similar fashion to Example 1. Unlike the linear case, we cannot find the covariance matrix of non-linear Bayesian inverse problems analytically. Hence, we implement Markov-Chain-Monte-Carlo (MCMC) methods; Specifically, we use the Metropolis-Adjusted-Langevin algorithm (MALA) [9]. MALA has an advantage over standard MCMC methods: The proposed samples are driven to areas of high probability in the posterior density using the density gradient. We gather a summary of our final results in Figure 9. The step size for MALA is fixed at $\tau = 5.2184 \times 10^{-13}$ with an acceptance rate of 69.1%, and we take 50,000 samples. Note that the misfits for the BFGS MAP, MCMC MAP and the conditional mean differ only in the sixth significative digit. All three results look very similar to each other, indicating that, for Bayesian AEI, a Gaussian approximation to the posterior may be sufficient. However, the skewness of the posterior suggests that such an approximation may have limited meaning: We observe high skewness for almost the entire field and extremely high skewness at the edges of the scatterers, see Figure 9 (g).

<figure><img src="image_10.png" /><figcaption>Figure 9: Toy example 2 – Uncertainty quantification of the AEI reconstruction: (a) exact scatterer term; (b) BFGS MAP estimate; (c) MCMC MAP estimate; (d) mean; (e) sample covariance matrix for the first 70 eigenfunctions; (f) MCMC pointwise variance of the reconstruction; (g) skewness of the reconstruction.</figcaption></figure>
For a closer look at the statistics, we display in Figure 10 the pairwise marginal posterior distributions for three AEI coefficients, $\beta_1$, $\beta_3$ and $\beta_4$, and compare them with the MAP estimate obtained from BFGS. We observe that $\beta_1$ and $\beta_4$ have left-skewed marginal posterior distributions, and their joint posterior is non-Gaussian. On the other hand, $\beta_3$ has a bimodal posterior and its joint posterior distributions with $\beta_1$ and $\beta_4$ are also non-Gaussian, consequently to the non-linearity of the problem.

<figure><img src="image_11.png" /><figcaption>Figure 10: Toy example 2 – MCMC joint posterior plots of $\beta_1$, $\beta_3$ and $\beta_4$. Top row: Eigenfunctions $\phi_1$, $\phi_3$ and $\phi_4$ at the final AEI adaptation. The red stars indicate the location of the BFGS solution for comparison purposes.</figcaption></figure>

Finally, there are numerous different dependencies that can be seen in the covariance matrix, as a result of the non-linear problem. The pointwise variance is high at the edges of the three scatterers. In Figures 16 and 17 of the appendix, we show the diagnostics for the convergence of example posterior distributions obtained from MCMC; the trace and ACF plots indicate poor mixing.
4.2.2 Shepp-Logan phantom revisited

We tested the algorithm with the modified Shepp-Logan phantom, scaled so that the maxi-
mum contrast is 2 and the background is 1. The data is at all points of the boundary and
the noise level has increased to 10% compared to 5% in the toy example. We now set four
sources surrounding the skull at points (0.1, 0.5), (0.9, 0.5), (0.5, 0.975) and (0.5, 0.025). In
the Bayesian AEI algorithm, we keep α fixed at 10⁻⁴ but increase the number of eigenfunc-
tions by 50 for each adaptation, starting from 50.

In Figure 11, we present the resulting reconstruction using Bayesian AEI together with $H^1$-Tikhonov and TV reconstructions, with regularisation parameters satisfying the discrepancy principle of Morozov. With Bayesian AEI, we find a smoother reconstruction than both Tikhonov and TV. We also have some idea of the uncertainty in our reconstruction, see Figure 12. We obtain the variance through MALA with step size $1.0858 \times 10^{-14}$, acceptance rate 66.3%, and 50000 samples. For more details about the trace and ACF plots of this sampling, we refer to Figures 18 and 19 in Appendix C.

<figure><img src="image_12.png" /><figcaption>Figure 11: Shepp-Logan phantom 2 – Reconstruction comparisons with the same display scale: (a) exact scatterer for comparison, (b) using Bayesian AEI, (c) using $H^1$-Tikhonov regularisation with $\alpha = 10^{-6}$, and (d) using TV regularisation with $\alpha = 5 \times 10^{-9}$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: Shepp-Logan phantom 2 – Uncertainty quantification of the AEI reconstruction: (a) exact scatterer term; (b) BFGS optimal solution; (c) MCMC MAP estimate; (d) conditional mean; (e) sample covariance matrix for the first 70 eigenfunctions; (f) MCMC pointwise variance of the reconstruction; (g) skewness of the reconstruction.</figcaption></figure>
For the sake of illustration, the selected evolution of the reconstructions is shown in Figure 13. The Bayesian AEI converged in 8 adaptations.

<figure><img src="image_14.png" /><figcaption>Figure 13: Shepp-Logan phantom 2 – The second (a), fourth (b), sixth (c) and eighth (d) adaptation to recover the scatterer. Each adaptation has a clarifying effect.</figcaption></figure>

As in the toy example, the algorithm is more uncertain about the edges of the scatterer. Finer details such as the skull are unclear in all reconstructions. However, choosing eigenfunctions based on their sensitivities is being further investigated and has the potential to recover these small-scale features [44]. Moreover, there is potential to improve performance using a new or existing variant of the BFGS algorithm alongside AEI, such as the limited-memory BFGS algorithm (L-BFGS) [65].

Similarly to the toy example, the skewness is overall high despite near identical MAP estimates from BFGS and MCMC, shown in Figure 11. We show examples of the Shepp-Logan phantom posterior distributions in Figure 20 of Appendix C. Upon closer look at the posterior distributions obtained from MCMC, the BFGS and MCMC MAP estimates generally do not match. The posterior distributions are non-Gaussian for both the toy and Shepp-Logan phantom examples. However, we have small-scale posterior distributions and BFGS MAP estimates within a probable range of these distributions. If we assume that the MCMC algorithm has sufficiently converged, then the histogram plots of the posterior distribution act as a close-up of the local minimum of the functional. In other words, the posterior histogram plots represent small imperfections in the BFGS reconstruction.

We find that in both the toy and Shepp-Logan phantom examples, the acceptance rate is restrictive to very small step sizes. MCMC methods generally have high computational costs; the calculation of the gradient in MALA is especially expensive. With enough computational budget, additional samples can be taken to obtain a more accurate posterior, but the qualitative difference in the MAP may still be insignificant. Currently, it is computationally infeasible to increase the number of samples. More work is needed to compare the performance of different sampling methods with Bayesian AEI. Furthermore, adapting Bayesian AEI to machine learning methods such as auto-differentiation [66] may accelerate the calculation of the gradient, which would offer faster sampling with MALA.

# 5 Conclusion

In this paper, we proposed a Bayesian extension to the Adaptive Eigenspace Inversion method. We first identified a suitable sequence of priors, then we proved that, under appropriate assumptions, the sequences of AEI posterior distributions converge. We also illustrated through various numerical examples the performance of Bayesian AEI.

In future research, we propose to explore several research directions. Some promising improvements of the AEI have already been proposed in [42, 44], namely in selecting an optimal number *K* of eigenfunctions at each adaptation, and could be implemented in the
Bayesian framework of AEI after some adjustment. Furthermore, ongoing studies about the
choice of parameter α in the penalty term and parameter ε in the diffusivity µ reveal that
the persistence of these parameters in the algorithm hinders performance. We observe that
a systematic update of these parameters will enhance the performance and almost certainly
guarantee convergence. This is under investigation.

Finally, from the MCMC results shown in the last section, the skewness results indicated
that the Laplace approximation to compute the covariance matrix may not be satisfactory.
However, the calculation of the gradient in MALA being particularly expensive, we plan
to explore machine learning techniques to help enhance the performance. Machine learning
libraries like PyTorch and Tensorflow could be useful in our future research [67, 68].

Acknowledgement

MG and BH were supported by the Marsden Fund from Government funding, Project
Nr 3725886 "Adaptive Eigenspace Inversion: Analysis, improved efficiency, and extensions",
managed by Royal Society Te Apārangi. BH was also supported by the National Science
Foundation Grants NSF-DMS-2208535, "Machine Learning for Bayesian Inverse Problems"
and NSF-CAREER-37678, "Gaussian Processes for Scientific Machine Learning: Theoretical
Analysis and Computational Algorithms".

References

[1] Moloud Abdar, Farhad Pourpanah, Sadiq Hussain, Dana Rezazadegan, Li Liu, Mo-
hammad Ghavamzadeh, Paul Fieguth, Xiaochun Cao, Abbas Khosravi, U. Rajendra
Acharya, Vladimir Makarenkov, and Saeid Nahavandi. A review of uncertainty quan-
tification in deep learning: Techniques, applications and challenges. *Information Fusion*,
76:243–297, 2021.

[2] Andrew M. Stuart. Inverse problems: A Bayesian perspective. *Acta Numerica*, 19:451–559, 2010.

[3] Daniela Calvetti and Erkki Somersalo. Inverse problems: From regularization to Bayesian inference. *Wiley Interdisciplinary Reviews: Computational Statistics*, 10(3):e1427, 2018.

[4] Jari Kaipio and Erkki Somersalo. *Statistical and computational inverse problems*, volume 160 of *Applied Mathematical Sciences*. Springer-Verlag New York, 2006.

[5] Ali Harakeh, Michael Smart, and Steven L. Waslander. BayesOD: A Bayesian Approach for Uncertainty Estimation in Deep Object Detectors. In *2020 IEEE International Conference on Robotics and Automation (ICRA)*, pages 87–93, 2020.

[6] Qingping Zhou, Tengchao Yu, Xiaoqun Zhang, and Jinglai Li. Bayesian inference and uncertainty quantification for medical image reconstruction with Poisson data. *SIAM Journal on Imaging Sciences*, 13(1):29–52, 2020.

[7] Yao-Tien Chen and Din-Chang Tseng. Medical image segmentation based on the Bayesian level set method. In *Medical Imaging and Informatics*, pages 25–34. Springer Berlin Heidelberg, 2008.
[8] P.A. Bromiley, N.A. Thacker, M.L.J. Scott, M. Pokrić, A.J. Lacey, and T.F. Cootes. Bayesian and non-Bayesian probabilistic models for medical image analysis. *Image and Vision Computing*, 21(10):851-864, 2003.

[9] Muhammad Izzatullah, Tristan van Leeuwen, and Daniel Peter. Bayesian seismic inversion: a fast sampling Langevin dynamics Markov chain Monte Carlo method. *Geophys. J. Int.*, 227:1523-1553, 2021.

[10] Tan Bui-Thanh, Omar Ghattas, James Martin, and Georg Stadler. A computational framework for infinite-dimensional Bayesian inverse problems Part I: The linearized case, with application to global seismic inversion. *SIAM Journal on Scientific Computing*, 35(6):A2494-A2523, 2013.

[11] Sean Berti, Mattia Aleardi, and Eusebio Stucchi. A computationally efficient Bayesian approach to full-waveform inversion. *Geophysical Prospecting*, 72(2):580-603, 2024.

[12] Matti Lassas and Samuli Siltanen. Can one use total variation prior for edge-preserving Bayesian inversion? *Inverse problems*, 20(5):1537, 2004.

[13] David L Donoho. Compressed sensing. *IEEE Transactions on information theory*, 52(4):1289-1306, 2006.

[14] Stan Z Li. *Markov random field modeling in image analysis*. Springer Science & Business Media, 2009.

[15] Masoumeh Dashti, Stephen Harris, and Andrew Stuart. Besov priors for Bayesian inverse problems. *Inverse Problems and Imaging*, 6(2):183-200, 2012.

[16] Matthew M Dunlop, Marco A Iglesias, and Andrew M Stuart. Hierarchical Bayesian level set inversion. *Statistics and Computing*, 27:1555-1584, 2017.

[17] Bamdad Hosseini. Two Metropolis-Hastings algorithms for posterior measures with non-Gaussian priors in infinite dimensions. *SIAM/ASA Journal on Uncertainty Quantification*, 7(4):1185-1223, 2019.

[18] Daniela Calvetti, Erkki Somersalo, and A Strang. Hierarchical Bayesian models and sparsity: $\ell_2$-magic. *Inverse Problems*, 35(3):035003, 2019.

[19] Daniela Calvetti, Monica Pragliola, Erkki Somersalo, and Alexander Strang. Sparse reconstructions from few noisy data: Analysis of hierarchical Bayesian models with generalized gamma hyperpriors. *Inverse Problems*, 36(2):025010, 2020.

[20] Zhilong Fang, Curt Da Silva, Rachel Kuske, and Felix J Herrmann. Uncertainty quantification for inverse problems with weak partial-differential-equation constraints. *Geophysics*, 83(6):R629-R647, 2018.

[21] Wences P Gouveia and John A Scales. Bayesian seismic waveform inversion: Parameter estimation and uncertainty analysis. *Journal of Geophysical Research: Solid Earth*, 103(B2):2759-2779, 1998.

[22] Albert Tarantola. Theoretical background for the inversion of seismic waveforms including elasticity and attenuation. *Pure and Applied Geophysics*, 128:365-399, 1988.
[23] Dimitris G. Tzikas, Aristidis C. Likas, and Nikolaos P. Galatsanos. The variational approximation for Bayesian inference. *IEEE Signal Processing Magazine*, 25(6):131–146, 2008.

[24] Aki Pulkkinen, Ben T Cox, Simon R Arridge, Hwan Goh, Jari P Kaipio, and Tanja Tarvainen. Direct estimation of optical parameters from photoacoustic time series in quantitative photoacoustic tomography. *IEEE transactions on medical imaging*, 35(11):2497–2508, 2016.

[25] Don S. Lemons and Anthony Gythiel. Paul Langevin's 1908 paper "On the Theory of Brownian Motion" ["Sur la théorie du mouvement brownien," C. R. Acad. Sci. (Paris) 146, 530-533 (1908)]. *American Journal of Physics*, 65(11):1079-1081, 11 1997.

[26] S. Brooks, A. Gelman, G. Jones, and X.-L. Meng. *Handbook of Markov Chain Monte Carlo*. Chapman and Hall/CRC, 1 edition, 2011.

[27] Max Welling and Yee W Teh. Bayesian learning via stochastic gradient Langevin dynamics. In *Proceedings of the 28th international conference on machine learning (ICML-11)*, pages 681–688. Citeseer, 2011.

[28] Arnak S Dalalyan and Lionel Riou-Durand. On sampling from a log-concave density using kinetic Langevin diffusions. *Bernoulli*, 26(3):1956–1988, 2020.

[29] SL Cotter, GO Roberts, AM Stuart, and D White. MCMC methods for functions: Modifying old algorithms to make them faster. *Statistical Science*, 28(3):424–446, 2013.

[30] Alexandros Beskos, Frank J Pinski, Jesús María Sanz-Serna, and Andrew M Stuart. Hybrid Monte Carlo on Hilbert spaces. *Stochastic Processes and their Applications*, 121(10):2201–2230, 2011.

[31] Alexandros Beskos, Mark Girolami, Shiwei Lan, Patrick E Farrell, and Andrew M Stuart. Geometric MCMC for infinite-dimensional inverse problems. *Journal of Computational Physics*, 335:327–351, 2017.

[32] Bamdad Hosseini and James E Johndrow. Spectral gaps and error estimates for infinite-dimensional Metropolis-Hastings with non-Gaussian priors. *The Annals of Applied Probability*, 33(3):1827–1873, 2023.

[33] Maya de Buhan and Axel Osses. Logarithmic stability in determination of a 3D viscoelastic coefficient and a numerical example. *Inverse Problems*, 26(9):95006, 2010.

[34] Maya de Buhan and Marie Kray. A new approach to solve the inverse scattering problem for waves: combining the TRAC and the Adaptive Inversion methods. *Inverse Problems*, 29(8):085009, 2013.

[35] Marcus J. Grote, Marie Kray, and Uri Nahum. Adaptive eigenspace method for inverse scattering problems in the frequency domain. *Inverse Problems*, 33(2):025006, 2017.

[36] Uri Nahum. *Adaptive eigenspace for inverse problems in the frequency domain*. PhD thesis, University of Basel, Switzerland, 2016.

[37] Marie Graff, Marcus J. Grote, Frédéric Nataf, and Franck Assous. How to solve inverse scattering problems without knowing the source term: a three-step strategy. *Inverse Problems*, 35(10):104001, 2019.
[38] Marie Graff and Bamdad Hosseini. Recent advances on the Adaptive Eigenspace Inversion method. In *Proceedings of The 16th International Conference on Mathematical and Numerical Aspects of Wave Propagation*, pages 458–459, 2024.

[39] Thomas Blumensath and Mike E Davies. Iterative thresholding for sparse approximations. *Journal of Fourier analysis and Applications*, 14:629–654, 2008.

[40] Francine Catté, Pierre-Louis Lions, Jean-Michel Morel, and Tomeu Coll. Image selective smoothing and edge detection by nonlinear diffusion. *SIAM Journal on Numerical analysis*, 29(1):182–193, 1992.

[41] Luis Alvarez, Pierre-Louis Lions, and Jean-Michel Morel. Image selective smoothing and edge detection by nonlinear diffusion. II. *SIAM Journal on numerical analysis*, 29(3):845–866, 1992.

[42] Daniel H Baffet, Marcus J Grote, and Jet Hoe Tang. Adaptive spectral decompositions for inverse medium problems. *Inverse Problems*, 37(2):025006, jan 2021.

[43] Daniel H. Baffet, Yannik G. Gleichmann, and Marcus J. Grote. Error estimates for adaptive spectral decompositions. *Journal of Scientific Computing*, 93(3):68, 2022.

[44] Yannik G Gleichmann and Marcus J Grote. Adaptive spectral inversion for inverse medium problems. *Inverse problems*, 39(12):125007, 2023.

[45] Maya de Buhan and Marion Darbas. Numerical resolution of an electromagnetic inverse medium problem at fixed frequency. *Computers and Mathematics with Applications*, 74(12):3111–3128, 2017.

[46] Marcus J. Grote and Uri Nahum. Adaptive eigenspace for multi-parameter inverse scattering problems. *Computers & Mathematics with Applications*, 77(12):3264–3280, June 2019.

[47] Florian Faucher, Otmar Scherzer, and Hélène Barucq. Eigenvector models for solving the seismic inverse problem for the Helmholtz equation. *Geophysical Journal International*, 221(1):394–414, 2020.

[48] Clay Sanders, Marc Bonnet, and Wilkins Aquino. An Adaptive Eigenfunction Basis Strategy to Reduce Design Dimension in Topology Optimization. *International Journal for Numerical Methods in Engineering*, 122:7452–7481, 2021.

[49] Izar Azpiroz. *Contribution to the numerical reconstruction in inverse elasto-acoustic scattering*. Theses, UPPA (LMA-Pau), February 2018.

[50] Abdalkalag Hamad and Mohsen Tadi. Inverse scattering based on proper solution space. *Journal of Theoretical and Computational Acoustics*, 27(03):1850033, 2019.

[51] Uri Nahum, Azhar Zam, and Philippe C. Cattin. Bone reconstruction and depth control during laser ablation. In *Computational Methods and Clinical Applications in Musculoskeletal Imaging*, pages 126–135. Springer International Publishing, 2019.

[52] Carlo Seppi, Uri Nahum, Peter A. von Niederhäusern, Simon Pezold, Michael Rissi, Stephan K. Haerle, and Philippe C. Cattin. Compressed sensing on multi-pinhole collimator SPECT camera for sentinel lymph node biopsy. In *Medical Image Computing and Computer-Assisted Intervention – MICCAI 2017*, pages 415–423. Springer International Publishing, 2017.
[53] João Carlos Alves Barata and Mahir Saleh Hussein. The Moore-Penrose pseudoinverse: A tutorial review of the theory. *Braz J Phys*, 42:146—165, 2012.

[54] Laurence Richard Lines and Sven Treitel. Tutorial: A Review of Least-Squares Inversion and its Application to Geophysical Problems. *Geophysical Prospecting*, 32:159–186, 1984.

[55] Pietro Perona and Jitendra Malik. Scale-space and edge detection using anisotropic diffusion. *IEEE Transactions on pattern analysis and machine intelligence*, 12(7):629–639, 2002.

[56] Alfredo Garbuno-Inigo, Tapio Helin, Franca Hoffmann, and Bamdad Hosseini. Bayesian posterior perturbation analysis with integral probability metrics. *arXiv preprint arXiv:2303.01512*, 2023.

[57] Björn Sprungk. On the local Lipschitz stability of Bayesian inverse problems. *Inverse Problems*, 36(5):055015, 2020.

[58] Cédric Villani et al. *Optimal transport: old and new*, volume 338. Springer, 2008.

[59] Anton Mallasto and Aasa Feragen. Learning from uncertain curves: The 2-Wasserstein metric for Gaussian processes. *Advances in Neural Information Processing Systems*, 30, 2017.

[60] Matthias Gelbrich. On a formula for the L2 Wasserstein metric between measures on Euclidean and Hilbert spaces. *Mathematische Nachrichten*, 147(1):185–203, 1990.

[61] Andreas Kirsch. *An Introduction to the Mathematical Theory of Inverse Problems*, volume 120 of *Applied Mathematical Sciences*. Springer-Verlag New York, 2nd edition, 2011.

[62] L. A. Shepp and B. F. Logan. The Fourier reconstruction of a head section. *IEEE Transactions on Nuclear Science*, 21(3):21–43, 1974.

[63] R-E Plessix. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications. *Geophysical Journal International*, 167(2):495–503, 2006.

[64] Jorge Nocedal and Stephen J Wright. *Numerical optimization*. Springer NY, 2nd edition, 2006.

[65] Dong C Liu and Jorge Nocedal. On the limited memory BFGS method for large scale optimization. *Mathematical programming*, 45(1):503–528, 1989.

[66] Weiqiang Zhu, Kailai Xu, Eric Darve, and Gregory C Beroza. A general approach to seismic inversion with automatic differentiation. *Computers & Geosciences*, 151:104751, 2021.

[67] Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer. Automatic differentiation in PyTorch. In *NIPS-W*, 2017.

[68] Martín Abadi, Paul Barham, Jianmin Chen, Zhifeng Chen, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Geoffrey Irving, Michael Isard, et al. {TensorFlow}: a system for {Large-Scale} machine learning. In *12th USENIX symposium on operating systems design and implementation (OSDI 16)*, pages 265–283, 2016.
# Appendices

## A Example 1 – AEI evolution

<figure><img src="image_15.png" /><figcaption>Figure 14: Example 1 – All 12 adaptations to recover the source term. The misfit is decaying from $6.2 \times 10^{-3}$ to $5.8 \times 10^{-4}$, reaching the level of noise. We started with 150 eigenfunctions for the first computation, then we kept the number of eigenfunctions fixed and equal to 160 for all subsequent adaptations. Colour scale fixed for all images.</figcaption></figure>
B Shepp-Logan Phantom – AEI evolution

<figure><img src="image_16.png" /><figcaption>Figure 15: Shepp-Logan phantom 1 – From top left to bottom right, all 20 adaptations for evolution visualization of the reconstruction.</figcaption></figure>
C Example 2 – MCMC diagnostics

<figure><img src="image_17.png" /><figcaption>Figure 16: Toy example 2 – Trace plots for $β_1, β_2, β_4, β_{38}$ and $β_{39}$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 17: Toy example 2 – ACF plot for $β_1, β_2, β_4, β_{38}$ and $β_{39}$.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 18: Shepp-Logan phantom 2 – Trace plots for $β_1, β_9, β_{10}, β_{41}$ and $β_{42}$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 19: Shepp-Logan phantom 2 – ACF plot for $β_1, β_9, β_{10}, β_{41}$ and $β_{42}$.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 20: Shepp-Logan phantom 2 – MCMC joint posterior plots of $β_1, β_9$ and $β_{10}$. Comparison with the BFGS MAP estimate indicated by red stars. Top row are the respective eigenfunctions for each AEI coefficient.</figcaption></figure>