# Adaptive Spectral Inversion: A new perspective to enhance the robustness in solving inverse problems

Nasrin Nikbakht, Marie Graff, Melissa Tacy

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

## To cite this version:

Nasrin Nikbakht, Marie Graff, Melissa Tacy. Adaptive Spectral Inversion: A new perspective to enhance the robustness in solving inverse problems. 2025. hal-05412427

HAL Id: hal-05412427

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

Preprint submitted on 12 Dec 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Adaptive Spectral Inversion: A new perspective to
enhance the robustness in solving inverse problems

Nasrin Nikbakht¹, Marie Graff¹,²,* and Melissa Tacy¹

¹ Department of Mathematics, The University of Auckland, Auckland 1010, New Zealand

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

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

12th December 2025

Abstract

In this paper, we propose to investigate in depth the robustness of the Adaptive Eigen-
space Inversion (AEI) method, developed over the last decade. Numerical results showed
very promising aspects of the method, but the lack of understanding and volatility of the
algorithm made it difficult to ensure the best result. Here, we study two regularisation
parameters involved in the method: (i) the number of eigenfunctions is determined by their
sensitivities, that is, how essential they are to build a descent direction in the optimisa-
tion algorithm, and (ii) the cancellation parameter in the diffusivity turns out to be more
essential and an optimal value of this parameter has been identified. In addition, we use
the Discrepancy Principle of Morozov to determine the number of adaptations required.
Backed up by numerical experiments, this new version of AEI, now called adaptive-ε ASI,
proves to be more efficient and more robust than its predecessor, and than standard tech-
niques, such as Tikhonov regularisation.

# 1 Introduction

Inverse problems are concerned with the recovery of the spatially varying properties of a source or scatterer, denoted by *p*, inside a bounded domain *Ω* ⊂ R<sup>*d*</sup> (*d* ≥ 1). Numerical methods for the solution of inverse problems essentially fall into two classes: qualitative and quantitative methods. Qualitative methods can estimate the location and the shape of parameter *p*, while quantitative methods are able to recover also some physical properties of it. Qualitative methods [3] generally require *a priori* knowledge about the type of scatterers *p* (sound-soft or sound-hard typically) and allow to estimate its location and shape efficiently. Examples are the MUSIC (Multiple Signal Classification) algorithm [5, 27], the DORT (Decomposition of the Time Reversal Operator) method [35, 4], the Linear Sampling method [6, 23], the probe method [34], and the factorization method [28, 29]. The Linear Sampling method, for instance,
provides a criterion for deciding whether a point is located inside or outside the scatterer $p$ from the knowledge of far-field measurements, while the DORT algorithm focusses on the decomposition of the time-dependent data matrix to extract the information. Hence, qualitative methods effectively determine the support of a scatterer or source and thus permit to detect material defects or obstacles quickly, and are dependent of the type of problems investigated (e.g., time domain vs frequency domain), but are not suitable for non-reflective scatterers or sources with varying wave propagation speed. When it comes to recover quantitative information, such as the wave propagation speed, the inverse problem is commonly reformulated as an optimization problem constrained by partial differential equations (PDE) [18, 21], in particular in applications such as medical imaging [31, 33, 30] or seismic imaging [22, 32, 37, 36, 40, 44]. In this setting, the unknown parameters are estimated iteratively by minimizing an objective functional that quantifies the misfit between simulated and observed data [41].

In this paper, we will focus on quantitative methods described as PDE-constrained optimisation problems. The minimisation process starts from an initial guess of the parameter $p$ and updates it during the optimisation as close as possible to the exact parameter. Regularisation techniques are essential for solving inverse problems because these problems are often ill-posed, meaning they lack stability, uniqueness, or even the existence of a solution. In practical applications, observational data are often noisy, incomplete, or inconsistent, making direct solutions to inverse problems highly unstable. Regularisation addresses these issues by introducing constraints or modifications that transform ill-posed problems into well-posed ones, ensuring stable and reliable solutions. For example, Tikhonov regularisation penalises large or unstable solutions by adding a smoothness constraint [14, 10, 43], while sparsity-promoting techniques leverage prior knowledge about the solution's structure, such as its sparsity or boundedness [11, 24, 25]. Additionally, regularisation mitigates numerical ill-conditioning, reducing the impact of small data errors on the solution [11, 24, 25, 42]. By balancing fidelity to the observed data with the imposition of desirable solution properties, regularisation ensures solutions are robust, physically meaningful, and computationally feasible.

The Adaptive Eigenspace Inversion (AEI) method was first introduced by de Buhan and Osses [9] in 2010. This approach restricts the search space to a low-dimensional eigenspace formed by the leading eigenfunctions of an elliptic differential operator $L[p]$. The operator $L$ is problem-independent and depends on $p$, the parameter we seek to reconstruct, making the search space inherently adaptive. This iterative adaptation allows the method to refine the search space at each step, aligning it with the evolving solution. To mitigate issues related to vanishing gradients, de Buhan and Kray [8] proposed a modified operator $L_ε[p]$, where a small parameter $ε > 0$ ensures the connectivity of the operator. This adjustment has been instrumental in improving numerical stability and robustness, making the method effective for ill-posed problems.

The AEI method has demonstrated its versatility across various applications. In inverse scattering problems, for example, it has been used to reconstruct media properties with high accuracy in acoustic [8, 15], electromagnetic [7], and seismic [12] settings. For seismic applica-
tions, the method leverages eigenfunctions to recover structures that are piecewise smooth or discontinuous, as demonstrated by Faucher et al. [12]. In optimal control, Sanders et al. [39] used the method to reduce design dimensions, enhancing computational efficiency for topology optimisation tasks. Beyond its practical success, the adaptive eigenspace method has been enriched by theoretical developments. Grote, Kray, and Nahum [19] extended the method by dynamically adjusting not only the eigenfunctions but also the dimension of the search space during iterations. This advancement ensures the search space is both computationally efficient and tailored to the problem's complexity at each step. Additionally, the connection between $L_ε$ and Total Variation (TV) regularisation [20] has opened avenues for its use in reconstructing piecewise constant and piecewise smooth solutions by tailoring the operator to the expected smoothness of the target medium.

The method's theoretical underpinnings have also been strengthened in recent years. Baffet et al. [1] established rigorous $L_2$-error estimates for adaptive spectral decomposition<fnref n="1" />, particularly in approximating piecewise constant functions. These results provide a solid mathematical foundation for the method's remarkable accuracy in such applications. In [13], a rigorous adaptive spectral inversion methodology is developed for inverse medium problems, in which the misfit is systematically projected onto a subspace of adaptively updated eigenfunctions of a carefully chosen elliptic operator. This framework inherently integrates regularisation and guarantees convergence under a precisely defined angle condition, providing a robust theoretical basis for accurate and stable reconstruction of spatially varying medium properties.

Our main contribution in this paper is twofold: (i) We implement the ASI version of the algorithm as proposed in Gleichmann and Grote [13] for an automatic update of the number of eigenfunctions in the truncation for the regularisation effect, but omitting the reshuffle and orthonormalisation of the basis to keep the inherent structure of the eigenvalue problem; (ii) We investigate the role of the second regularisation parameter, $\epsilon$, which has been chosen fixed and small in previous work. However, the fixed $\epsilon$ may not be optimal across spatial locations or iterative steps, potentially undermining both the stability and accuracy of the reconstruction. We propose a systematically updated, spatially and iteratively adaptive parameter $\epsilon$, a strategy that, to the best of our knowledge, has not been previously investigated. We also proved that the optimal value of $\epsilon$ depends on the norm of the gradient and the Laplacian of parameter $p$. The numerical results demonstrate that the proposed adaptive choice of $\epsilon$ is more efficient than the fixed strategy. In particular, the adaptive approach yields improved accuracy in terms of the relative $L^2$ error and provides a more accurate recovery of $p$. Moreover, when compared with the standard Tikhonov regularisation, the ASI with adaptive $\epsilon$ achieves superior reconstruction quality and enhanced stability.

The remaining part of this paper is structured as follows. In Section 2, we first reintroduce the AEI technique for solving inverse problems. Then we describe, ASI regularisation technique, extension of AEI for an automatic update of the number of eigenfunctions, and apply

<footnote n="1">¹The Adaptive Spectral Decomposition method is an extension to the Adaptive Eigenspace Inversion method.</footnote>
it to reconstruct $p$, with the goal of evaluating its effectiveness in addressing the eigenfunction selection challenge. It is important to note that, although we are inspired by [13], our approach to determining the basis $\Phi = \{\phi_m\}_{m \ge 1}$ for the eigenfunction representation does not involve the reshuffle of the basis as presented in [13]. In the remaining of this section, we present our new methodology in determining and adapting parameter $\varepsilon$, later denoted by $\varepsilon_a$, in operator $L_{\varepsilon}$. By allowing $\varepsilon$ to vary, we aim to achieve a more flexible and potentially improve regularisation strategy. Finally, in Section 3, we present several numerical experiments which illustrate the successful performance of the adaptive $\varepsilon_a$ in ASI method. In particular, we apply the ASI method with adaptive $\varepsilon_a$ to a time-harmonic inverse source and medium problems which demonstrates a significant improvement over previous versions from [2, 19, 13].

## 2 Methodology

In this section, we consider a time-harmonic acoustic inverse problem, where the wavefield $u(x)$ satisfies the Helmholtz equation inside a bounded domain $\Omega \subset \mathbb{R}^d$ ($d \ge 1$), together with a Sommerfeld-type radiation condition at the (artificial) boundary $\partial\Omega$.

$$ \left\{ \begin{array}{ll} -(k^2 + \nabla^2) u = f & \text{in } \Omega, \\ \frac{\partial u}{\partial n} + iku = 0 & \text{on } \partial\Omega. \end{array} \right. \qquad (2.1) $$

Here, $k$ corresponds the wavenumber and $f$ denotes the source term. To solve the inverse problem, we frame it as a PDE-constrained optimisation problem. Let us now denote by $p$ the parameter to estimate, either the wavenumber $k$ or the source term $f$. Then we seek for the minimiser $p$ of the standard $L^2$ data misfit functional,

$$ \underset{p \in L^2(\Omega)}{\operatorname{argmin}} \quad \mathcal{J}(p) = \underset{p \in L^2(\Omega)}{\operatorname{argmin}} \frac{1}{2} \|u_{\text{obs}} - u[p]\|_{L^2(\Gamma)}^2, \qquad (2.2) $$

where $\Gamma \subset \Omega$ denotes the region in which the receivers are positioned. The goal is to reconstruct $p$ such that the simulated state $u[p]$ best matches the observed data $u_{\text{obs}}$. The observed data $u_{\text{obs}}$ are modelled as

$$ u_{\text{obs}} = u[p_{\text{true}}]|_{\Gamma} + \eta_{\delta}. \qquad (2.3) $$

Here, $\eta_{\delta} \sim \mathcal{N}(0, \delta^2)$ represents additive white noise with variance $\delta^2$ and $\delta$ is the noise level.

As such inverse problems are usually ill-posed, regularisation techniques like AEI are addressed to solve (2.2). In AEI, the function $p: \Omega \subset \mathbb{R}^d \to \mathbb{R}$ ($d \ge 1$) is taken to be an arbitrary bounded function on a bounded Lipschitz domain. Within this framework, $p$ is represented using a basis of global eigenfunctions $\Phi = \{\phi_m\}_{m \ge 1}$, avoiding reliance on standard nodal basis functions typically used in finite difference or finite element methods. Specifically, $p$ is expressed as:

$$ p(x) = p_0 + \sum_{m=1}^{\infty} \beta_m \phi_m(x) \quad \forall x \in \Omega. \qquad (2.4) $$
Here, $p_0$ is defined as the background term. In this paper, we will only consider constant background $p_0$, typically for homogeneous media, but it may also be seen as a varying background, see [19]. For inverse scattering problems, $p_0$ coincides with the wavenumber, whereas for source problems $p_0 = 0$. In (2.4), the basis of eigenfunctions $\Phi = \{\phi_m\}_{m \ge 1}$ are obtained as solutions to the elliptic problem:

$$ \mathcal{L}_{\varepsilon}[p]\phi_m = -\nabla \cdot (\mu_{\varepsilon}[p]\nabla\phi_m) = \lambda_m \phi_m \quad \text{in } \Omega, \quad \phi_m = 0 \quad \text{on } \partial\Omega. \qquad (2.5) $$

By nature, this elliptic eigenvalue problem has a discrete but infinite spectrum with corresponding eigenvalues $0 < \lambda_m \le \lambda_{m+1}$, $m \ge 1$. The diffusivity $\mu$ is defined by

$$ \mu_{\varepsilon}[p] = \frac{1}{\sqrt{|\nabla p|^{2} + \varepsilon^{2}}}, \quad \text{in } \Omega. \qquad (2.6) $$

Here, the parameter $\varepsilon > 0$ ensures that the denominator remains well-defined, typically set to a small positive value such as $\varepsilon = 10^{-6}$. Since $p$ is precisely the unknown quantity we seek, in $\mu$ we use the value from the previous adaptation step $j$, hence the adaptive process.

$$ \mu_{\varepsilon}[p]^{(j)} = \frac{1}{\sqrt{|\nabla p^{(j-1)}|^{2} + \varepsilon^{2}}} \quad j \ge 2 \quad \text{in } \Omega. \qquad (2.7) $$

In the initial step, $j = 1$, where no prior information about $p$ within $\Omega$ is yet available, we set $\mu_{\varepsilon}[p] = 1$, which leads to harmonic prolongation of the boundary data and Laplacian-based eigenfunctions.

In (2.4), the coefficients $\beta = \{\beta_m\}_{m \ge 1}$ are computed by solving the following optimisation problem:

$$ \min_{\beta} \frac{1}{2} \|u_{\text{obs}} - u[p_0 + \Phi\beta]\|_{L^2(\Gamma)}^2, \quad \Gamma \subset \Omega. \qquad (2.8) $$

The computation of $\beta$ depends on the nature of the problem. In inverse source problems, it can be obtained directly using the linear least squares approach. For inverse scattering problems the Gauss-Newton method in conjunction with Armijo line search is employed to determine $\beta$, for instance. Actually, any optimisation method derived from the Newton's method can be chosen.

In practice, we truncate the infinite sum in Equation (2.4) at a finite number $M \ge 1$:

$$ p(x) = p_0 + \sum_{m=1}^{M} \beta_m \phi_m(x), \quad \forall x \in \Omega. \qquad (2.9) $$

We summarise the AEI algorithm, highlighting the adaptive process, in Algorithm 1.
**Algorithm 1:** AEI Algorithm

**Input:** Background $p_0$, observations $u_{\text{obs}}$

**Output:** $p^*$

1. Choose $M \ge 1$, and compute $\Phi_M^{(1)} = \{\phi_m^{(1)}\}_{m=1}^M$ from (2.5) with $\mu_{\varepsilon}[p] = 1$.

2. Solve the optimisation problem (2.8) to obtain $\beta^{(1)}$.

3. Find $p^{(1)}(x) = p_0 + \sum_{m=1}^M \beta_m^{(1)} \phi_m^{(1)}(x)$.

4. STOP: If misfit $r^{(j)} = \|u_{\text{obs}} - u[p^{(j)}]\| \le \delta$:

* Update $\mu_{\varepsilon}[p^{(j)}]$ from (2.7).

* Solve the eigenvalue problem (2.5) to get the basis $\Phi_M^{(j)}$.

* Solve the optimisation problem (2.8) to obtain $\beta^{(j)}$.

* Find $p^{(j)}(x) = p_0 + \sum_{m=1}^M \beta_m^{(j)} \phi_m^{(j)}(x)$.

* Set $j \leftarrow j + 1$.

5. Set $p^* = p^{(j)}$.

**Remark 1.** The number of adaptations to be performed in Algorithm 1 is determined according to the Discrepancy Principle of Morozov [26]. The main idea behind this principle is that the misfit cannot be better than the noise on the data, denoted by $\delta$. Therefore, the AEI algorithm is building a sequence of approximated solutions and stops when the lattest solution gives a misfit of about $\delta$. This pressuppose that we have some knowledge on the level of noise on the data, or at least some estimation of it. △

**Remark 2.** The elliptic operator specified by Equation (2.5) essentially coincides with the gradient of the penalised TV regularisation term [38, 45]:

$$R_{\text{TV}}(p) = \frac{1}{2} \int_{\Omega} \sqrt{|\nabla p|^2 + \varepsilon^2} dx, \qquad (2.10)$$

whose gradient is given by:

$$\nabla R_{\text{TV}}(p) = -\nabla \cdot \left( \frac{1}{\sqrt{|\nabla p|^2 + \varepsilon^2}} \nabla p \right) \approx -\nabla \cdot (\mu_{\varepsilon}[p] \nabla p). \qquad (2.11)$$

In the application of the AEI method, one of the significant challenges lies in selecting the appropriate number of eigenfunctions, denoted by $M$, to reconstruct the source or scatterer. To address this uncertainty, Gleichmann and Grote [13] introduced the Adaptive Spectral Inversion (ASI), a method that ranks eigenfunctions based on their sensitivities to improve the selection process.
## 2.1 ASI: Ranking Eigenfunctions

The sensitivity of each eigenfunction is a crucial metric that determines its influence on minimising the cost function $\mathcal{J}[p^{(j)}]$. Where $p^{(j)}$ is the current recovery of source or scatterer at adaptation $j$ with basis $\Phi_M = \{\phi_m\}_{m \ge 1}^M$. The sensitivity of an eigenfunction $\phi_m^{(j+1)}$ at the next iteration $j + 1$ is computed as:

$$ \sigma_m^{(j+1)} = \nabla \mathcal{J}[p^{(j)}]^\top \phi_m^{(j+1)}, \text{ where } 1 \le m \le M^+. \quad (2.12) $$

Here, $M^+ \ge M$ and $\nabla \mathcal{J}[p^{(j)}]$ denotes the gradient of the cost function with respect to the current recovery. This sensitivity measures the alignment between the gradient, which provides the direction of steepest descent to reduce the data misfit, and the eigenfunction $\phi_m^{(j+1)}$, derived from the spectral decomposition of the elliptic operator $\mathcal{L}_\epsilon[p^{(j+1)}]$ at the next iteration. A higher sensitivity magnitude indicates that the eigenfunction has a stronger impact on improving the reconstruction, while a lower magnitude suggests limited relevance.

The new eigenspace basis $\Phi_{M^+}^{(j+1)} = \{\phi_m^{(j+1)}\}_{m \ge 1}^{M^+}$ focuses on the most informative components, enhancing computational efficiency, stability, and accuracy. We then truncate the eigenbasis $\Phi_{M^+}^{(j+1)}$ and just use the first $N_\theta^{(j+1)} \le M^+$ for the next recovery. Therefore, $p^{(j+1)}$ is expressed as:

$$ p^{(j+1)}(x) = p_0 + \sum_{m=1}^{N_\theta^{(j+1)}} \beta_m^{(j+1)} \phi_m^{(j+1)}(x). \quad (2.13) $$

Here,

$$ N_{\theta}^{(j+1)} = \max\{N_{0}^{(j+1)}, N_{2}^{(j+1)}, N_{\infty}^{(j+1)}\}. \quad (2.14) $$

The finite numbers $N_0^{(j+1)}$, $N_2^{(j+1)}$, and $N_\infty^{(j+1)}$ are defined as outlined in the following.

**Index** $N_0^{(j+1)}$: Let $\gamma_m$ represent the orthogonal projection of $p^{(j)}$ onto the basis of eigenfunctions at the next iteration $\Phi_{M^+}^{(j+1)}$, i.e,

$$ \gamma_m^{(j+1)} = \langle p^{(j)}, \phi_m^{(j+1)} \rangle_{L^2(\Omega)}, \quad \text{where } 1 \le m \le M^+. \quad (2.15) $$

Then for a fixed $0 < \varepsilon_\Phi < 1$, $N_0^{(j+1)} \le M$ represents the smallest index such that

$$ \sum_{m=N_0}^{M} |\gamma_m^{(j+1)}|^2 \le \varepsilon_\Phi^2 \|\gamma^{(j+1)}\|_{\ell^2}^2. \quad (2.16) $$

Indices $N_2^{(j+1)}$ and $N_\infty^{(j+1)}$ are defined by the following theorem and lemmas. The corresponding proofs can be found in [13].

**Theorem 1.** (Theorem 2.1 in [13]). Let the misfit $\mathcal{J}$ be Fréchet differentiable with derivative $\nabla \mathcal{J}(p)$, $p \in H$. Here $H$ denotes a Hilbert space continuously embedded in $L^2(\Omega)$. Moreover, $\mathcal{J}$ is assumed to be Lipschitz-continuous for every direction $d \in H$, i.e.

$$ |\nabla \mathcal{J}(p+q)^T d - \nabla \mathcal{J}(p)^T d| \le L \|q\|_H \|d\|_H, \quad \forall p, q \in H, L > 0, \quad (2.17) $$
and further assume that $\mathcal{J}(p) \ge C > -\infty$ for all $p \in H$. Then if the corrections $d^{(j)} = p^{(j+1)} - p^{(j)}$ satisfies the angle condition

$$|\nabla \mathcal{J}(p^{(j)})^T d^{(j)}| \ge \varepsilon_{\theta} \|d^{(j)}\|_H \|\nabla \mathcal{J}(p^{(j)})\|_{L^2}, \quad 0 < \varepsilon_{\theta} < 1 \qquad (2.18)$$

uniformly in $j$, we have

$$\|\nabla \mathcal{J}(p^{(j)})\|_{L^2} \to 0, \quad j \to \infty. \qquad (2.19)$$

**Lemma 1.** (Lemma 4.1 in [13]). For a fixed $0 < \varepsilon_{\theta} < 1$, let $N_{\infty}^{(j+1)} \le M^+$ be the largest index such that

$$|\sigma_m^{(j+1)}| \ge \varepsilon_{\theta} \|\nabla \mathcal{J}(p^{(j)})\|_{L^2}, \quad m = 1, \dots, N_{\infty}^{(j+1)}, \qquad (2.20)$$

holds. Then, $d = \sum_{m=1}^{N_{\infty}} \sigma_m^{(j+1)} \phi_m^{(j+1)}$ satisfies the angle condition (2.18).

**Lemma 2.** (Lemma 4.2 in [13]). Let $N_2^{(j+1)} \le M^+$ be the smallest index such that

$$\sqrt{\sum_{m=1}^{N_2^{(j+1)}} |\sigma_m^{(j+1)}|^2} \ge \varepsilon_{\theta} \|\nabla \mathcal{J}(p^{(j)})\|_{L^2}, \qquad (2.21)$$

holds. Then, $d = \sum_{m=1}^{N_2} \sigma_m^{(j+1)} \phi_m^{(j+1)}$ satisfies the angle condition (2.18).

**Remark 3.** In [13], the new search space $\Psi^{(j+1)}$ is obtained by merging the current search space $\Psi^{(j)} \subset L^2(\Omega)$ with the newly computed eigenbasis $\Phi^{(j+1)}$. This space is then orthonormalised and ordered by their sensitivities. Finally, it is truncated to construct a reduced basis suitable for the adaptation $j+1$. In comparison, our approach to constructing the new eigenspace does not involve merging with the preceding basis and ordering eigenfunctions. Instead, we directly employ the new basis $\Phi_{M^+}^{(j+1)}$ and construct the $p^{(j+1)}$ using the basis $\Phi_{N_{\theta}}^{(j+1)} \subset \Phi_{M^+}^{(j+1)}$. The rational for keeping the eigenspace untouched is a potential Bayesian extension of the method, which requires to build a prior probability distribution based on the eigenvalues associated to the eigenfunctions in the basis [16, 17]. Therefore, we want to keep the eigenvalue structure and avoid reshuffling. As a result, our approach is actually a variation of [13]. △

Since we only compute a finite number of eigenfunctions $\Phi_{M^+} = \{\phi_m\}_{m=1}^{M^+}$, arranged by their eigenvalues, it is possible that none of them satisfy condition (2.20). In such cases, we set $N_{\infty} = 0$, which leads to $N_{\theta} = \max\{N_0, N_2\}$. If no $N_2 \le M^+$ satisfies condition (2.21), the value of $M^+$ must be increased until (2.21) is met. This ensures the existence of $d \in \Phi_{N_{\theta}}$ that satisfies condition (2.18). This yields the full ASI Algorithm 2.
**Algorithm 2:** ASI Algorithm

**Input:** Background $p_0$, observations $u_{\text{obs}}$

**Output:** $p^*$

1. Choose $M \ge 1$, and compute $\Phi_M^{(1)} = \{\phi_m^{(1)}\}_{m=1}^M$ from (2.5) with $\mu_{\varepsilon}[p] = 1$.

2. Solve the optimisation problem (2.8) to obtain $\beta^{(1)}$.

3. Find $p^{(1)}(x) = p_0 + \sum_{m=1}^{M} \beta_m^{(1)} \phi_m^{(1)}(x)$.

4. Choose $M^+ \ge M$.

5. STOP: If misfit $r^{(j)} = \|u_{\text{obs}} - u[p^{(j)}]\| \le \delta$:

• Update $\mu_{\varepsilon}[p^{(j)}]$ from (2.7).

• Solve the eigenvalue problem (2.5) to get the basis $\Phi_{M^+}^{(j)}$.

• Compute the sensitivities to determine $N_{\theta}^{(j)}$ from (2.14).

• Solve the optimisation problem (2.8) in the basis $\Phi_{N_{\theta}^{(j)}}^{(j)} \subset \Phi_{M^+}^{(j)}$ to obtain $\beta^{(j)}$.

• Find $p^{(j)}(x) = p_0 + \sum_{m=1}^{N_{\theta}^{(j)}} \beta_m^{(j)} \phi_m^{(j)}(x)$.

• Set $j \leftarrow j + 1$.

6. Set $p^* = p^{(j)}$.

## 2.2 ASI with Adaptive $\varepsilon_a$

Within the frameworks of AEI and ASI, the regularisation parameter $\varepsilon$ in the definition of the diffusivity function $\mu_{\varepsilon}[p]$ in (2.6) is typically assigned a fixed value. This parameter serves to guarantee numerical stability, specifically by avoiding division by zero. However, a fixed choice of $\varepsilon$ may not be optimal across different spatial regions or iterations, potentially affecting the accuracy and stability of the solution. To address this issue, we analytically determine the interval corresponding to $\varepsilon$, denoted by $\varepsilon_a$. By allowing $\varepsilon$ to vary, we aim to achieve a more flexible and potentially improved regularisation strategy.

In (2.7), we replace the fixed $\varepsilon$ with an adaptive regularisation parameter $\varepsilon_a > 0$ to get the following.

$$ \mu_{\varepsilon_a}[p] = \frac{1}{\sqrt{|\nabla p|^2 + \varepsilon_a^2}}, \quad \text{in } \Omega. \qquad (2.22) $$

Our goal is to determine an appropriate form for $\varepsilon_a$. To avoid affecting the behaviour of the diffusivity function $\mu_{\varepsilon_a}[p]$ negatively and to ensure that the adaptive method ASI remains effective, it is important to control the size of $\varepsilon_a$ appropriately. Specifically, $\varepsilon_a$ should not
exceed the maximum magnitude of the local gradient. In other words, we impose the condition

$$ \varepsilon_a \leq \max_{x \in \Omega} |\nabla p(x)|. \qquad (2.23) $$

This requirement ensures that $\varepsilon_a$ does not dominate the gradient term in the denominator of $\mu_{\varepsilon_a}[p]$, which could otherwise lead to excessive smoothing and a loss of important features in the reconstruction. We now apply a linear operator to $p$ to analyse its transformation under the given functional framework

$$ |\mathcal{L}_{\varepsilon_a}[p(x)]| = \left| p_0 + \sum_{m=1}^{N_{\theta}} \beta_m \mathcal{L} \phi_m(x) \right|, \quad \forall x \in \Omega. \qquad (2.24) $$

Based on the eigenvalue problem defined in Equation (2.5), the action of the elliptic operator $\mathcal{L}$ can be expressed as

$$ \mathcal{L} \phi_m = \lambda_m \phi_m, \quad 1 \le m \le N_{\theta}, \qquad (2.25) $$

where $\lambda_m$ are the corresponding eigenvalues. We obtain

$$ |\mathcal{L}_{\varepsilon_a}[p(x)]| = \left| p_0 + \sum_{m=1}^{N_{\theta}} \beta_m \lambda_m \phi_m(x) \right|, \quad \forall x \in \Omega. \qquad (2.26) $$

By defining $\lambda_{\max} = \max_{1 \le m \le N_{\theta}} \lambda_m$ and $\beta_{\max} = \max_{1 \le m \le N_{\theta}} |\beta_m|$, it follows that

$$ |\mathcal{L}_{\varepsilon_a}[p(x)]| \le N_{\theta} \cdot \beta_{\max} \cdot \lambda_{\max}. \qquad (2.27) $$

In what follows, we set

$$ \tau = N_{\theta} \cdot \beta_{\max} \cdot \lambda_{\max}. \qquad (2.28) $$

To motivate the construction of the adaptive $\varepsilon_a$ and understand its behaviour better, we begin by analysing the problem in one spatial dimension. This simplified setting offers a tractable framework for examining how the adaptive regularisation affects the solution. In developing $\varepsilon_a$, particular care is taken to ensure that the regularisation process does not interfere with the established upper bound on $|\mathcal{L}_{\varepsilon_a}[p(x)]|$ in Equation (2.27). Preserving this bound is essential for maintaining the theoretical consistency of the formulation. The insights gained in the one-dimensional analysis are subsequently extended to the higher-dimensional case, where the same guiding principles are applied.

For the sake of simplicity, the following analysis will be performed in one-space dimension. Justification for an extension to higher dimensions is given at the end of the proof. We now define a function $h: \Omega \subset \mathbb{R} \to \mathbb{R}$ by

$$ h = \frac{d}{dx} \left( \frac{1}{\sqrt{p_x^2 + \varepsilon_a^2}} p_x \right). \qquad (2.29) $$

Here $p_x$ denotes the first partial derivative of $p$ with respect to $x$, and we consider $\varepsilon_a$, to be independent of $x$. We aim at

$$ |h(x)| \le |\mathcal{L}_{\varepsilon_a}[p(x)]|, \quad \forall x \in \Omega. \qquad (2.30) $$
Applying the product and chain rules, we compute $h(x)$ as follows:

$$h = \left(\frac{d}{dx}p_x\right) \cdot (p_x^2 + \varepsilon_a^2)^{-\frac{1}{2}} + p_x \cdot \left(\frac{d}{dx}(p_x^2 + \varepsilon_a^2)^{-\frac{1}{2}}\right) = \frac{p_{xx} \cdot \varepsilon_a^2}{(p_x^2 + \varepsilon_a^2)^{\frac{3}{2}}}. \quad (2.31)$$

Then, taking absolute values, we have

$$|h(x)| = \frac{|p_{xx}(x)| \cdot \varepsilon_a^2}{(|p_x(x)|^2 + \varepsilon_a^2)^{\frac{3}{2}}}, \quad \forall x \in \Omega. \quad (2.32)$$

Following Equation (2.32), it is natural to partition the domain $\Omega$ into two subdomains according to the relation between $\varepsilon_a$ and $|p_x|$. This yields the following decomposition.

**Definition 1** (Domain decomposition based on $\varepsilon_a$). Following the derivation of admissible relations between $\varepsilon_a$ and $|p_x|$ in Equation (2.32), we partition the domain $\Omega$ into two disjoint subdomains, $\Omega_1$ and $\Omega_2$, such that $\Omega = \Omega_1 \cup \Omega_2$, and $\Omega_1 \cap \Omega_2 = \emptyset$.

$$\Omega_1 = \left\{ x \in \Omega \mid 0 \le \frac{|p_x|}{\varepsilon_a} < 1 \right\}, \quad \Omega_2 = \left\{ x \in \Omega \mid \frac{|p_x|}{\varepsilon_a} \ge 1 \right\}. \quad (2.33)$$

Based on the identified structure of the subdomains $\Omega_1$ and $\Omega_2$, the next step is to derive possible bounds on $\varepsilon_a$.

**Theorem 2** (Admissible bounds for $\varepsilon_a$). Let regularisation parameter $0 < \varepsilon_a < \infty$. If

$$\frac{\|p_{xx}\|_{L^{\infty}(\Omega)}}{\tau} \le \varepsilon_a \le \|p_x\|_{L^{\infty}(\Omega)}. \quad (2.34)$$

Then

$$|h(x)| \le \tau, \quad \forall x \in \Omega. \quad (2.35)$$

*Proof.* Let $\Omega_1$ and $\Omega_2$ be the subdomains defined in Definition 1. First, consider $x \in \Omega_1$, then we have

$$0 \le |p_x(x)|_{(\Omega_1)} < \varepsilon_a \Rightarrow 0 \le |p_x(x)|^2_{(\Omega_1)} < \varepsilon_a^2. \quad (2.36)$$

Adding $\varepsilon_a^2$ to each side, it follows that

$$\varepsilon_a^2 \le |p_x(x)|^2_{(\Omega_1)} + \varepsilon_a^2 < 2\varepsilon_a^2. \quad (2.37)$$

Subsequently taking the reciprocal of each side raised to the power $\frac{3}{2}$, we get

$$\frac{1}{[\sqrt{2}\varepsilon_a]^3} < \frac{1}{[|p_x(x)|^2_{(\Omega_1)} + \varepsilon_a^2]^{3/2}} \le \frac{1}{\varepsilon_a^3}, \quad (2.38)$$

and after multiplying by $|p_{xx}(x)|_{(\Omega_1)}\varepsilon_a^2$

$$\frac{|p_{xx}(x)|_{(\Omega_1)}\varepsilon_a^2}{[\sqrt{2}\varepsilon_a]^3} < \frac{|p_{xx}(x)|_{(\Omega_1)}\varepsilon_a^2}{[|p_x(x)|^2_{(\Omega_1)} + \varepsilon_a^2]^{3/2}} \le \frac{|p_{xx}(x)|_{(\Omega_1)}\varepsilon_a^2}{\varepsilon_a^3}. \quad (2.39)$$

By Equation (2.32), it follows that

$$|h(x)| \le \frac{|p_{xx}(x)|\varepsilon_a^2}{\varepsilon_a^3}, \quad \forall x \in \Omega_1. \quad (2.40)$$
Taking the maximum and noting that here $\Omega_1 \subset \Omega$

$$ \forall x \in \Omega_1 : |p_{xx}(x)| \leq \max |p_{xx}(x)| = \|p_{xx}\|_{L^{\infty}(\Omega_1)} \leq \|p_{xx}\|_{L^{\infty}(\Omega)}, \quad (2.41) $$

it follows that

$$ |h(x)| \leq \frac{\|p_{xx}\|_{L^{\infty}(\Omega)}}{\varepsilon_a}, \quad \forall x \in \Omega. \quad (2.42) $$

Finally, by selecting $\varepsilon_a$ to satisfy

$$ \varepsilon_a \geq \frac{\|p_{xx}\|_{L^{\infty}(\Omega)}}{\tau}, \quad (2.43) $$

we ensure that (2.42) holds and that the condition in (2.35) is fulfilled.

We next examine the case $x \in \Omega_2$, i.e.,

$$ \text{for } x \in \Omega_2 : \quad \varepsilon_a \leq |p_x(x)|. \quad (2.44) $$

It follows that

$$ 0 < \varepsilon_a \leq |p_x(x)|_{(\Omega_2)} \Rightarrow 0 < \varepsilon_a^2 \leq |p_x(x)|_{(\Omega_2)}^2. \quad (2.45) $$

Adding $|p_x(x)|_{\Omega_2}^2$ to each side yields

$$ |p_x(x)|_{(\Omega_2)}^2 < \varepsilon_a^2 + |p_x(x)|_{(\Omega_2)}^2 \leq 2|p_x(x)|_{(\Omega_2)}^2, \quad (2.46) $$

then, taking the reciprocal of each side raised to the power $\frac{3}{2}$, we get

$$ \frac{1}{[\sqrt{2}|p_x(x)|_{(\Omega_2)}]^3} \leq \frac{1}{[|p_x(x)|_{(\Omega_2)}^2 + \varepsilon_a^2]^{3/2}} < \frac{1}{|p_x(x)|_{(\Omega_2)}^3}. \quad (2.47) $$

After multiplying by $|p_{xx}(x)|_{(\Omega_2)}\varepsilon_a^2$, we have

$$ \frac{|p_{xx}(x)|_{(\Omega_2)}\varepsilon_a^2}{[\sqrt{2}|p_x(x)|_{(\Omega_2)}]^3} \leq \frac{|p_{xx}(x)|_{(\Omega_2)}\varepsilon_a^2}{[|p_x(x)|_{(\Omega_2)}^2 + \varepsilon_a^2]^{3/2}} < \frac{|p_{xx}(x)|_{(\Omega_2)}\varepsilon_a^2}{|p_x(x)|_{(\Omega_2)}^3}, \quad (2.48) $$

then, by Equation (2.32), it follows that

$$ |h(x)| < \frac{|p_{xx}(x)|\varepsilon_a^2}{|p_x(x)|^3}, \quad \text{where } x \in \Omega_2. \quad (2.49) $$

Taking the maximum and noting that here $\Omega_2 \subset \Omega$

$$ \forall x \in \Omega_2 : |p_{xx}(x)| \leq \max |p_{xx}(x)| = \|p_{xx}\|_{L^{\infty}(\Omega_2)} \leq \|p_{xx}\|_{L^{\infty}(\Omega)}, \quad (2.50) $$

it follows that

$$ |h(x)| < \frac{\|p_{xx}\|_{L^{\infty}(\Omega)}\varepsilon_a^2}{|p_x(x)|^3}, \quad \text{where } x \in \Omega_2. \quad (2.51) $$

Then, from (2.43), it follows that

$$ |h(x)| < \frac{(\tau \varepsilon_a) \varepsilon_a^2}{|p_x(x)|^3} = \frac{\tau \varepsilon_a^3}{|p_x(x)|^3}, \quad \text{where } x \in \Omega_2. \quad (2.52) $$

This implies that the condition in (2.35) is fulfilled. $\square$

<figure><img src="image_3.png" /></figure>
Building on the insights gained from Theorem 2, we now identify the following formula for adaptive $\epsilon_a$ based on geometric average of lower and upper bounds of $\epsilon_a$.

$$ \epsilon_a = \sqrt{\frac{\|p_{xx}\|_{L^{\infty}(\Omega)} \|p_x\|_{L^{\infty}(\Omega)}}{\tau}} \quad (2.53) $$

In higher spatial dimensions the effect of truncation by $\epsilon$ is, in general, more complicated. However, the structure of our sources/scatterers allow for a dimensional reduction. We are studying sources/scatterers that are piecewise constant. Therefore, locally near the boundary of a jump in value of $p$ the derivative in the direction normal to the boundary greatly dominates that in the tangential direction (indeed if the jump is across a straight line the tangential derivative is zero). Therefore, the two-dimensional problem reduces to a one dimensional one (in the normal direction). As a result, away from corners we can expect the one dimensional calculation to guide the behaviour in higher dimensions.

**Definition 2** (Higher-Dimensional Extension of Adaptive $\epsilon_a$). Let $p: \Omega \subset \mathbb{R}^d \to \mathbb{R}$ ($d \ge 2$), be a function defined on a bounded domain $\Omega$. Let $\|\nabla p\|_{L^{\infty}(\Omega)}$ and $\|\Delta p\|_{L^{\infty}(\Omega)}$ denote the maximum (infinity) norms of the gradient and Laplacian of $p$, respectively. Then, the adaptive regularisation parameter $\epsilon_a$ is given by

$$ \epsilon_a = \sqrt{\frac{\|\Delta p\|_{L^{\infty}(\Omega)} \|\nabla p\|_{L^{\infty}(\Omega)}}{\tau}} \quad (2.54) $$

Here, $\tau$ is given by Equation (2.28).

Since the function $p$ represents the unknown we seek, it cannot be used directly in the evaluation of $\epsilon_a$. Therefore, in each adaptation step $j \ge 2$, we compute $\epsilon_a$ using the approximation of $p$ obtained from the previous adaptation.

$$ \epsilon_a^{(j)} = \sqrt{\frac{\|\Delta p^{(j-1)}\|_{L^{\infty}(\Omega)} \|\nabla p^{(j-1)}\|_{L^{\infty}(\Omega)}}{\tau^{(j)}}} \quad (2.55) $$

Here

$$ \tau^{(j)} = N_{\theta}^{(j-1)} \cdot \beta_{\max}^{(j-1)} \cdot \lambda_{\max}^{(j-1)} \quad (2.56) $$

In the initial step $j = 1$, where no prior information about $p$ within $\Omega$ is yet available, we set $\epsilon_a^{(1)} = 1$. The final algorithm is proposed in Algorithm 3.
**Algorithm 3: ASI-ε<sub>a</sub> Algorithm**

**Input:** Background $p_0$, observations $u_{\text{obs}}$

**Output:** $p^*$

1. Choose $M \ge 1$, compute $\Phi_M^{(1)}$ and determine $\lambda_{\max}^{(1)}$ from (2.5) with $\mu_{\varepsilon_a}[p] = 1$, and $\varepsilon_a = 1$.

2. Solve the optimisation problem (2.8) to obtain $\beta^{(1)}$. Determine $\beta_{\max}^{(1)}$.

3. Find $p^{(1)}(x) = p_0 + \sum_{m=1}^{M} \beta_m^{(1)} \phi_m^{(1)}(x)$.

4. Choose $M^+ \ge M$, and set $N_{\theta}^{(1)} = M$.

5. STOP: If misfit $r^{(j)} = \|u_{\text{obs}} - u[p^{(j)}]\|_{L^2(\Gamma)} \le \delta$:

* Compute $\tau^{(j)}$ from (2.56).

* Adaptive $\varepsilon_a$. Update adaptive $\varepsilon_a^{(j)}$ from $\Delta p^{(j-1)}$ and $\nabla p^{(j-1)}$ via (2.55).

* Update $\mu_{\varepsilon_a}[p^{(j)}]$ from $p^{(j-1)}$ and $\varepsilon_a^{(j)}$ via (2.22).

* Compute $\Phi_{M^+}^{(j)}$ from (2.5), and find $\lambda_{\max}^{(j)}$.

* Compute the sensitivities to determine $N_{\theta}^{(j)}$ from (2.14).

* Solve the optimisation problem (2.8) in the basis $\Phi_{N_{\theta}^{(j)}}^{(j)} \subset \Phi_{M^+}^{(j)}$ to obtain $\beta^{(j)}$.

* Find $p^{(j)}(x) = p_0 + \sum_{m=1}^{N_{\theta}^{(j)}} \beta_m^{(j)} \phi_m^{(j)}(x)$.

* Set $j \leftarrow j + 1$.

6. Set $p^* = p^{(j)}$.

**Remark 4.** At each adaptation $j$ of ASI with $\varepsilon_a^{(j)}$, we enforce the condition

$$ \|\Delta p^{(j)}\|_{L^{\infty}(\Omega)} \le \tau^{(j)} \|\nabla p^{(j)}\|_{L^{\infty}(\Omega)}. $$

Whenever this requirement fails, the maximal number of eigenfunctions $M^+$, given before truncation using the sensitivities, is increased until the inequality is satisfied. The rationale for this procedure is to ensure that $\varepsilon_a^{(j)}$ fulfils condition (2.34) in Theorem 2. $\Delta$

# 3 Numerical experiments

In this section, we evaluate the impact of the adaptive $\varepsilon_a$ on the performance of ASI. We also use our approach of ASI as described in Section 2, thus avoiding the reshuffle of the basis proposed in [13]. For the sake of transparence, we compare our method to the fixed-$\varepsilon$ approach. Specifically, we first consider a fixed value $\varepsilon = 10^{-6}$. However, recognising that the adaptive regularisation parameter $\varepsilon_a$ may become significantly larger than this fixed value during the
adaptation process, we introduce three additional fixed ε values for comparison: ε_min, defined as the minimum of ε_a over all adaptations; ε_max, the maximum value of ε_a; and ε_ave, the average of ε_a across all adaptation steps. We perform numerical experiments using these three fixed parameters, as well as ε = 10⁻⁶ and the adaptive ε_a, in order to assess the effectiveness of the adaptive strategy. This allows us to compare the performance of the adaptive approach not only against a small fixed ε, but also against fixed values that are closely related to the adaptive ε_a itself. In all experiments, the ASI algorithm 2 is applied for fixed ε (ε, ε_min, ε_ave, ε_max), whereas the adaptive case ε_a employs the ASI-ε_a algorithm 3. To prevent any inverse crimes, the exact data u[p_true] is computed on a separate finer mesh, which does not contain the computational mesh used for the inversion, and white noise is added.

We compare the ASI further with adaptive ε_a (Algorithm 3) to a standard grid-based inversion method employing Tikhonov regularisation [14, 10, 43]. For the standard Tikhonov L²-regularisation approach, we solve the minimisation problem:

$$ \min_{p} \mathcal{J}(p), \quad \text{with} \quad \mathcal{J}(p) = \frac{1}{2} \|u_{\text{obs}} - u[p]\|_{L^2(\Gamma)}^2 + \frac{\alpha}{2} \|p\|_{L^2(\Gamma)}^2, \quad \Gamma \subset \Omega, \quad (3.1) $$

where α > 0 denotes the regularisation parameter. The optimal choice of α is determined using Discrepancy Principle of Morozov [26] to ensure consistency with the noise level in the observed data u obs.

To evaluate the accuracy of the ASI method, we monitor the following quantities:

* the dimension N (j) θ during the recovery process,

* the misfit

$$ r^{(j)} = \|u_{\text{obs}} - u[p^{(j)}]\|_{L^2(\Gamma)}, \quad (3.2) $$

* the relative error

$$ e_{L^2} = \frac{\|p_{\text{true}} - p^*\|^2}{\|p_{\text{true}}\|^2}, \quad (3.3) $$

* and the maximum norm error

$$ e_{L^\infty} = \max_{x \in \Gamma} |p_{\text{true}}(x) - p^*(x)|, \quad (3.4) $$

where p true denotes the true profile of the parameter p, while p* represents its final reconstruction obtained using the ASI technique.

## 3.1 Application to the inverse source problem

For this first example, we focus on inverse source problem for the Helmholtz equation. The forward problem (2.1) is solved in the bounded domain Ω = [0, 1]². We employ second-order staggered finite differences on a 80 × 80 Cartesian mesh for the discretisation of Equation (2.1). We assume that the wavenumber k is constant and known in (2.1). The goal is to reconstruct the source term p = f from measurements of the wavefield at specific sensor positions Γ ⊂ Ω, that is, we selected randomly 1,000 nodes in the computational domain. These measurements
are used to infer the unknown sources by ensuring that the modelled wavefield matches the observed data, on which 10% of white noise has been added. Note that these settings are not meant to be realistic but are designed to illustrate the performance of our method, in particular, in comparison with standard methods like Tikhonov regularisation.

We shall to recover the two distinct sources independently, $p = f_1$ and $p = f_2$, illustrated in Figure 3.1. The first exhibits a face toy-like profile, whereas $f_2$ comprises two inclusions, one a crescent moon-shaped and the other a heart-like. In all cases, we apply the ASI algorithm 2 and ASI-$\epsilon_a$ algorithm 3 with the following fixed parameter settings:

$$k = 50, p_0 = 0, M^+ = 400, \epsilon_\theta = 0.1, \epsilon_\Phi = 0.05. \qquad (3.5)$$

<figure><img src="image_4.png" /><figcaption>Figure 3.1: Inverse Source Problem – True sources $f_1$ (left) and $f_2$ (right), used to test the performance of the ASI algorithm. The forward model is the Helmholtz equation (2.1) with fixed wavenumber $k = 50$.</figcaption></figure>

### 3.1.1 Source Example 1 - Toy face

The ASI method is initialised with the background $p_0 = 0$. The initial eigenfunction basis $\Phi^{(1)}$ is chosen as the first $M = 6$ Laplace eigenfunctions. We recall that $M$ is the actual number of eigenfunctions, whereas $M^+ = 400$ is the maximum number of eigenfunctions before truncation of the basis according to the sensitivities. Since the inverse source problem is linear, the optimisation part of the algorithm consists simply in using the linear Least-Squares method. The stopping criterion for the adaptations is based on the Discrepancy Principle of Morozov, for which the level of noise is needed. In this example, we assume that we have 10% of white noise, so that the absolute noise is $\delta \approx 6.087 \times 10^{-3}$. The intensity of the source to be recovered is 1.5 in the ears, 1.2 in the face and 0.5 in the eyes and mouth, as pictured in Figure 3.1(a).

Figure 3.2 provides a comparative visualisation of source reconstructions corresponding to $f_1$. All results are obtained using the ASI method with either fixed or adaptive $\epsilon$ strategies. The original source is displayed in panel (a) in the top left. The reconstruction using an adaptive $\epsilon_a$ is shown in panel (b) of Figure 3.2, while the one obtained with the fixed strategy $\epsilon$, with $\epsilon = 10^{-6}$, is presented in panel (c). Panels (d), (e), and (f) illustrate the reconstructions based on $\epsilon_{\min}$ (the minimum value of $\epsilon_a$), $\epsilon_{\text{ave}}$ (the average value of $\epsilon_a$) and $\epsilon_{\max}$ (the maximum
value of $\epsilon_a$) during the adaptation process. The number of adaptations is restricted to a maximum of 50. Consequently, cases that reach this upper limit are regarded as unsuccessful in recovering the prescribed sources, as illustrated in panels (c), (d) in Figure 3.2.

<figure><img src="image_5.png" /><figcaption>Figure 3.2: Source Example 1 – Final reconstructions of source $p = f_1$ using ASI with adaptive and fixed $\epsilon$ alongside the exact source. The corresponding dimensions, values of $\epsilon$, final misfits, and errors for each panel are provided in Table 3.1. The exponent in brackets gives the number of adaptations needed to convergence, 50 meaning that the process did not converge. The stopping criterion is based on the level of noise, here, $\delta \approx 6.087 \times 10^{-3}[10\%]$.</figcaption></figure>

<table><caption>Table 3.1: Source Example 1 – Recovery of $p = f_1$ with ASI. Comparison of reconstruction performance for different parameter selection strategies in terms of final dimension, adaptation count, misfit, and errors, for each value of $\epsilon$ tested. The noise level is set to $\delta \approx 6.087 \times 10^{-3}[10\%]$.</caption><thead><tr><th></th><th>Value</th><th>Final N<sub>θ</sub></th><th>Final j</th><th>Misfit r</th><th>e<sub>L2</sub></th><th>e<sub>L∞</sub></th></tr></thead><tbody><tr><td>ε<sub>a</sub></td><td>-</td><td>119</td><td>9</td><td>5.23 × 10<sup>-3</sup></td><td>26.22%</td><td>1.288</td></tr><tr><td>ε</td><td>10<sup>-6</sup></td><td>46</td><td>50</td><td>1.410 × 10<sup>-2</sup></td><td>48.17%</td><td>2.145</td></tr><tr><td>ε<sub>min</sub></td><td>6.207 × 10<sup>-2</sup></td><td>74</td><td>50</td><td>1.018 × 10<sup>-2</sup></td><td>34.03%</td><td>1.155</td></tr><tr><td>ε<sub>ave</sub></td><td>2.337 × 10<sup>-1</sup></td><td>149</td><td>25</td><td>5.653 × 10<sup>-3</sup></td><td>27.51%</td><td>1.553</td></tr><tr><td>ε<sub>max</sub></td><td>2.987 × 10<sup>-1</sup></td><td>170</td><td>17</td><td>6.010 × 10<sup>-3</sup></td><td>28.1%</td><td>1.201</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>(a) $f_1$: Evolution of adaptive $\epsilon_a$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) $f_1$: Dimension $N_\theta$ vs. iteration $j$</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(c) $f_1$: Misfit $r$ vs. iteration $j$</figcaption></figure>

Figure 3.3: Source Example 1 – (a) Evolution of $\epsilon_a$ at each adaptation in ASI method. (b) Different selection of eigenfunctions between adaptive and fixed $\epsilon$ settings. (c) Convergence patterns of recovery of source $f_1$ with ASI using adaptive and fixed $\epsilon$. The noise level set to $\delta \approx 6.087 \times 10^{-3}[10\%]$. The cases fixed $\epsilon$ and $\epsilon_{min}$ are omitted as they did not converge.
The results in Figure 3.2 highlight a clear performance advantage when employing an adaptive $\epsilon_a$. The ASI method with a fixed $\epsilon$ and $\epsilon_{min}$ fail to produce an accurate reconstruction. These reconstructions do not reflect a converged solution and underscore the limitations of the fixed strategy in these cases. Moreover, the recovery using the fixed strategies of $\epsilon_{ave}$ and $\epsilon_{max}$ are not as efficient as adaptive strategy $\epsilon_a$.

Table 3.1 presents a comparative evaluation of the different parameter choices in the performance of the fixed and adaptive strategies of $\epsilon$. The comparison is based on the final dimensions $N_\theta$, final number of adaptations $j$, misfit values, $e_{L2}$ and $e_{L\infty}$ corresponding to the reconstruction of the $f_1$ using ASI. ASI with fixed $\epsilon = 10^{-6}$ and $\epsilon_{min} = 6.207 \times 10^{-2}$ fail to provide satisfactory recovery, since their misfit does not converge to the noise level $\delta \approx 6.087 \times 10^{-3}[10\%]$ within the maximum number of 50 adaptations. By contrast, the adaptive $\epsilon_a$ yield a successful recovery while requiring fewer adaptations compared to the fixed $\epsilon_{ave} = 2.337 \times 10^{-1}$ and $\epsilon_{max} = 2.987 \times 10^{-1}$. Although the $e_{L\infty}$ associated with $\epsilon_a$ is slightly larger than those obtained with $\epsilon_{min}$ and $\epsilon_{max}$, it is important to note that, in terms of the $e_{L2}$, $\epsilon_a$ performs significantly better than the other parameter choices. Since the $e_{L2}$ reflects the global distribution of errors across the entire domain, this result indicates that $\epsilon_a$ provides a more accurate overall recovery, even if the pointwise maximum error is marginally higher.

Figure 3.3 provides a detailed examination of the behaviour of key parameters during the ASI procedure, with particular emphasis on $\epsilon_a$, $\epsilon_{ave}$, and $\epsilon_{max}$ which successfully recover $f_1$. Specifically, panel (a) presents the evolution of $\epsilon_a$ throughout the adaptation process, beginning from iteration $j = 2$. Panel (b) compares the selection of eigenfunctions under two approaches—adaptive versus fixed $\epsilon$—emphasising the differences in basis selection strategies. Lastly, panel (c) illustrates the convergence patterns observed in the recovery of source $f_1$, for both adaptive and fixed $\epsilon$ settings. As shown in panel (c) in Figure 3.3, $\epsilon_a$ reaches the noise threshold $\delta \approx 6.087 \times 10^{-3}[10\%]$ at an earlier stage of the adaptations compared with $\epsilon_{ave}$ and $\epsilon_{max}$. This behaviour demonstrates that the adaptive strategy $\epsilon_a$ enhances the convergence of the ASI procedure and thereby improves its overall effectiveness in achieving accurate recovery.

### 3.1.2 Source Example 2 - Crescent and heart shapes

As for the previous example, we initialise with the background $p_0 = 0$. The initial eigenfunction basis $\Phi^{(1)}$ is now chosen as the first $M = 4$ Laplace eigenfunctions. Since the problem is still linear, we use linear least-squares for the optimisation and the stopping criterion for the adaptation is based on the level noise, here, $\delta \approx 2.776 \times 10^{-3}$, corresponding to 10%.

As previously, Figure 3.4 provides a comparative visualisation of source reconstructions corresponding to $f_2$. All results are obtained using the ASI method with either fixed or adaptive $\epsilon$ strategies. The original source is displayed in panel (a) in the top left. The reconstruction using an adaptive $\epsilon_a$ is shown in panel (b), while the one obtained with the fixed strategy $\epsilon$ is presented in panel (c). Panels (d), (e), and (f) illustrate the reconstructions based on $\epsilon_{min}$, $\epsilon_{ave}$ and $\epsilon_{max}$, respectively. Again, the number of adaptations is restricted to a maximum of 50. As a result, the unsuccessful cases in recovering the prescribed sources are now for $\epsilon_{min}$ and $\epsilon_{ave}$,
as illustrated in panels (d) and (e) in Figure 3.4.

<figure><img src="image_9.png" /><figcaption>Figure 3.4: Source Example 2 – Final reconstructions of source $p = f_2$ using ASI with adaptive and fixed $\varepsilon$ alongside the exact source. The corresponding dimensions, values of $\varepsilon$, final misfits, and errors for each panel are provided in Table 3.2. The exponent in brackets gives the number of adaptations needed to convergence, 50 meaning that the process did not converge. The stopping criterion is based on the level of noise, here, $\delta \approx 2.776 \times 10^{-3}$[10%].</figcaption></figure>

<table><caption>Table 3.2: Source Example 2 – Recovery of source $p = f_2$ with ASI. Comparison of reconstruction performance for different parameter selection strategies in terms of final dimension, adaptation count, misfit and errors, for each value of $\varepsilon$ tested. The noise level is set to $\delta \approx 2.776 \times 10^{-3}$[10%].</caption><thead><tr><th></th><th>Value</th><th>Final N<sub>θ</sub></th><th>Final j</th><th>Misfit r</th><th>e<sub>L2</sub></th><th>e<sub>L∞</sub></th></tr></thead><tbody><tr><td>ε<sub>a</sub></td><td>-</td><td>128</td><td>8</td><td>2.66 × 10<sup>-3</sup></td><td>25.6%</td><td>1.961</td></tr><tr><td>ε</td><td>10<sup>-6</sup></td><td>142</td><td>10</td><td>2.447 × 10<sup>-3</sup></td><td>28.13%</td><td>1.554</td></tr><tr><td>ε<sub>min</sub></td><td>7.631 × 10<sup>-2</sup></td><td>23</td><td>50</td><td>8.291 × 10<sup>-3</sup></td><td>40.04%</td><td>1.721</td></tr><tr><td>ε<sub>ave</sub></td><td>5.020 × 10<sup>-1</sup></td><td>72</td><td>50</td><td>7.590 × 10<sup>-3</sup></td><td>36.61%</td><td>1.523</td></tr><tr><td>ε<sub>max</sub></td><td>1.065</td><td>196</td><td>9</td><td>1.883 × 10<sup>-3</sup></td><td>29.40%</td><td>3.302</td></tr></tbody></table>

Figure 3.4, in conjunction with Table 3.2, presents a comparative visualisation of the reconstructions of $f_2$ obtained using the ASI method with fixed $\varepsilon$ and adaptive $\varepsilon_a$ strategies. In this configuration, the noise level is set to $\delta \approx 2.776 \times 10^{-3}$[10%], that is, 10% of white noise on the data. The result obtained with the ASI with $\varepsilon_a$ (panel (b) in Figure 3.4), demonstrates
<figure><img src="image_10.png" /><figcaption>(a) $f_2$: Evolution of adaptive $\epsilon_a$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) $f_2$: Dimension $N_\theta$ vs. iteration $j$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(c) $p_2$: Misfit $r$ vs. iteration $j$</figcaption></figure>

Figure 3.5: Source Example 2 – (a) Evolution of $\epsilon_a$ at each adaptation in ASI method. (b) Different selection of eigenfunctions between adaptive and fixed $\epsilon$ settings. (c) Convergence patterns of recovery of source $f_2$ with ASI using adaptive and fixed $\epsilon$. The noise level set to $\delta \approx 2.776 \times 10^{-3}$[10%]. The cases fixed $\epsilon_{\min}$ and $\epsilon_{\text{ave}}$ are omitted as they did not converge.
successful recovery. In contrast, the reconstructions associated with $\epsilon_{\min}$ and $\epsilon_{\text{ave}}$, displayed in panels (d) and (e), respectively, do not achieve a successful recovery of $f_2$ within the prescribed maximum of 50 adaptations. Among the strategies that succeed, the ASI method with $\epsilon_a$ yields the most accurate reconstruction, characterised by a smooth representation along the boundary and, in particular, a well-defined recovery of the heart-shaped structure.

The results in Table 3.2 show that for both $e_{L^2}$ and $e_{L^\infty}$, the ASI method using the adaptive $\epsilon_a$ for the recovery of $f_2$ exhibits the same behaviour as in the recovery of $f_1$. The $e_{L^\infty}$ error corresponding to $\epsilon_a$ is marginally larger than those obtained with $\epsilon = 10^{-6}$, $\epsilon_{\min} = 7.631 \times 10^{-2}$ and $\epsilon_{\text{ave}} = 5.020 \times 10^{-1}$. In contrast, the $e_{L^2}$ of $\epsilon_a$ is substantially smaller. This indicates that $\epsilon_a$ achieves a more accurate recovery overall, despite the marginally higher pointwise maximum error.

A detailed analysis of the behaviour of the key parameters in the ASI procedure is provided in Figure 3.5. Particular emphasis is placed on $\epsilon_a$, $\epsilon$, and $\epsilon_{\max}$, all of which achieve successful recovery of $f_2$. As shown in panel (c) of Figure 3.5, the misfit in the case $\epsilon_a$ attains the noise threshold $\delta \approx 2.776 \times 10^{-3}[10\%]$ at an earlier stage of the adaptations compared with $\epsilon$, and $\epsilon_{\max}$. This observation indicates that the adaptive strategy governed by $\epsilon_a$ accelerates the convergence of the ASI procedure and, thereby, enhances its overall effectiveness in producing an accurate recovery.

### 3.1.3 Comparison with Tikhonov regularisation

Finally, we propose to compare our ASI approach with standard methods, namely, the Tikhonov regularisation as in (3.1). For a fair comparison, the regularisation parameter $\alpha$ in the penalty term is determined using the Discrepancy Principle of Morozov, that is, we select the parameter $\alpha$ that minimises the cost function and for which the misfit reaches the level of noise $\delta$.

Figure 3.6 compares the recovered profiles for $p = f_1$ and $p = f_2$, as depicted in Figure 3.1, obtained with ASI using adaptive $\epsilon_a$ and the standard grid-based inversion method using Tikhonov regularisation. In all cases, the results obtained using the Tikhonov method (panels (c) and (f)) demonstrate inferior performance relatively to those achieved with ASI using adaptive $\epsilon_a$ (panels (b) and (e)). Clearly, with the $L^2$-Tikhonov regularisation, the background appears noisy and the coefficients inside the sources are not accurately reconstructed, more precisely, even though the shape is good, the value within each feature is not correct. Note for instance how the eyes and mouth of the toy-face example keep the background value while they should not.

Remarkably, the $e_{L^2}$ obtained by the ASI with adaptive $\epsilon_a$ method is smaller than that achieved by standard $L^2$-Tikhonov regularisation, as presented in Table 3.3. Both methods successfully provide a reconstruction such that the misfit has reached the level of noise, that is, satisfies the Discrepancy Principle of Morozov. From the $e_{L^\infty}$, we can conclude that there must exist some nodes where the error is higher for ASI, but the $e_{L^2}$ relative error shows that globally the ASI reconstruction matches better the exact source, which was confirmed in Figure 3.6.
We can conclude that ASI outperformed the Tikhonov regularisation.

<figure><img src="image_13.png" /><figcaption>Figure 3.6: ASI vs Tikhonov regularisation – Final reconstruction of the selected sources using ASI with adaptive $\epsilon_a$ and Tikhonov regularisation, alongside the exact sources $p = f_1$ and $p = f_2$. The analysis is performed in the presence of $\delta \approx 6.087 \times 10^{-3}$[10%] for recovery of $p = f_1$, and $\delta \approx 2.776 \times 10^{-3}$[10%] for recovery of $p = f_2$, see Table 3.3 for the respective misfits and errors.</figcaption></figure>

<table><caption>Table 3.3: ASI vs Tikhonov regularisation – Comparison of reconstruction performance for different parameter selection strategies in terms of misfit $r$, regularisation parameter $\alpha$ (Tikhonov regularisation only), and errors. The noise level is set to $\delta \approx 6.087 \times 10^{-3}$[10%] for recovery of $p = f_1$, and $\delta \approx 2.776 \times 10^{-3}$[10%] for recovery of $p = f_2$.</caption><thead><tr><th rowspan="2"></th><th colspan="2">p = f<sub>1</sub></th><th colspan="2">p = f<sub>2</sub></th></tr><tr><th>ASI-ε<sub>a</sub></th><th>Tikhonov</th><th>ASI-ε<sub>a</sub></th><th>Tikhonov</th></tr></thead><tbody><tr><td>Misfit r</td><td>5.23 × 10<sup>-3</sup></td><td>5.943 × 10<sup>-3</sup></td><td>2.66 × 10<sup>-3</sup></td><td>2.728 × 10<sup>-3</sup></td></tr><tr><td>Parameter α</td><td>-</td><td>4.962 × 10<sup>-8</sup></td><td>-</td><td>2.662 × 10<sup>-8</sup></td></tr><tr><td>e<sub>L<sup>2</sup></sub></td><td>26.22%</td><td>43.60%</td><td>25.6%</td><td>35.08%</td></tr><tr><td>e<sub>L<sup>∞</sup></sub></td><td>1.288</td><td>1.277</td><td>1.961</td><td>1.335</td></tr></tbody></table>

## 3.2 Application to inverse scattering problem

To illustrate the accuracy and usefulness of the ASI-$\epsilon_a$ Algorithm 3, we shall now apply it to inverse scattering (medium) problems. Similarly to the source-recovery settings, the forward
problem is posed on the bounded computational domain $\Omega = [0, 1]^2$. Following the same approach discretisation of Equation (2.1) is carried out using second-order staggered finite differences on a uniform Cartesian grid of size 80×80. We use observations from eight equispaced Gaussian sources located along the boundary $\partial\Omega$. The goal is to reconstruct the medium term $p = k$ from measurements of the wave amplitude and phase (frequency domain) at specific sensor positions $\Gamma = \{(x, y) \in \Omega : \text{dist}((x, y), \partial\Omega) \le 0.1\} \subset \Omega$, that is, we selected randomly 320 nodes in the $\Gamma$.

For given the eight sources $f$ and noisy data $u_{\text{obs}}$ as in (2.3), the recovery is carried out independently for the two distinct unknown media $p = k_1$ and $p = k_2$ illustrated in Figure 3.7. The first medium exhibits two distinct inclusions, one a crescent moon-shaped and the other a heart-like, whereas $p = k_2$ consists of five distinct inclusions of various sizes, shapes and values. Both profiles $k_1$ and $k_2$ are set to a uniform background $p_0 = k$ where $k$ is a central wavenumber corresponding to the main frequency. In all cases, we apply the ASI algorithm 2 and ASI-$\epsilon_a$ algorithm 3 with the following fixed parameter settings:

$$p_0 = 50, M^+ = 400, \epsilon_{\theta} = 0.1, \epsilon_{\Phi} = 0.05. \qquad (3.6)$$

<figure><img src="image_14.png" /><figcaption>Figure 3.7: Inverse Medium Problem – True media $k_1$ (left) and $k_2$ (right), used to test the performance of the ASI algorithm. The forward model is the Helmholtz equation (2.1) with eight fixed equispaced Gaussian sources located along the boundary $\partial\Omega$.</figcaption></figure>

### 3.2.1 Medium Example 1 - Crescent and heart shapes

The ASI method is initialised with the background $p_0 = 50$. The initial eigenfunction basis $\Phi^{(1)}$ is chosen as the first $M = 5$ Laplace eigenfunctions. Since the inverse medium problem is nonlinear, within the $j^{th}$ step of the ASI Algorithm, the optimisation problem (2.8) is solved using the Gauss Newton algorithm in conjunction with Armijo line search to determine the coefficients $\beta$, where we recall that $p = \Phi\beta$ in our representation. As for the source recovery, the stopping criterion for the adaptations is based on the Discrepancy Principle of Morozov, for which the level of noise is needed. In this example, we assume that we have 10% of white noise, so that the absolute noise is $\delta \approx 2.561 \times 10^{-4}$. The wavenumber of the medium to be recovered is 62 in the crescent and 56 in the heart shape, as pictured in Figure 3.7(a).
<figure><img src="image_15.png" /><figcaption>Figure 3.8: Medium Example 1 – Final reconstruction of medium $p = k_1$ using ASI with adaptive and fixed $\varepsilon$ alongside the exact medium $k_1$. The corresponding dimensions, values of $\varepsilon$, final misfits, and errors for each panel are provided in Table 3.4. The exponent in brackets gives the number of adaptations needed to convergence, 50 meaning that the process did not converge. The stopping criterion is based on the level of noise, here, $\delta \approx 2.561 \times 10^{-4}[10\%]$.</figcaption></figure>

<table><caption>Table 3.4: Medium Example 1 – Recovery of medium $p = k_1$ with ASI. Comparison of reconstruction performance for different parameter selection strategies in terms of final dimension, adaptation count, misfit, and errors, for each value of $\varepsilon$ tested. The noise level set to $\delta \approx 2.561 \times 10^{-4}[10\%]$.</caption><thead><tr><th></th><th>Value</th><th>Final N<sub>θ</sub></th><th>Final j</th><th>Misfit r</th><th>e<sub>L</sub><sup>2</sup></th><th>e<sub>L</sub><sup>∞</sup></th></tr></thead><tbody><tr><td>ε<sub>a</sub></td><td>-</td><td>26</td><td>6</td><td>2.486 × 10<sup>-4</sup></td><td>2.81%</td><td>11.362</td></tr><tr><td>ε</td><td>10<sup>-6</sup></td><td>23</td><td>4</td><td>2.518 × 10<sup>-4</sup></td><td>3.02%</td><td>11.058</td></tr><tr><td>ε<sub>min</sub></td><td>4.417 × 10<sup>-1</sup></td><td>29</td><td>10</td><td>2.467 × 10<sup>-4</sup></td><td>3.07%</td><td>14.983</td></tr><tr><td>ε<sub>ave</sub></td><td>1.348</td><td>24</td><td>15</td><td>2.556 × 10<sup>-4</sup></td><td>3.41%</td><td>12.463</td></tr><tr><td>ε<sub>max</sub></td><td>2.753</td><td>36</td><td>11</td><td>2.522 × 10<sup>-4</sup></td><td>3.45%</td><td>19.093</td></tr></tbody></table>
<figure><img src="image_16.png" /><figcaption>(a) $k_1$: Evolution of adaptive $\epsilon_a$</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(b) $k_1$: Dimension $N_\theta$ vs. iteration $j$</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(c) $k_1$: Misfit $r$ vs. iteration $j$</figcaption></figure>

Figure 3.9: Medium Example 1 - (a) Evolution of $\epsilon_a$ at each adaptation in ASI method. (b) Different selection of eigenfunctions between adaptive and fixed $\epsilon$ settings. (c) Convergence patterns of recovery of Medium $k_1$ with ASI using adaptive and fixed $\epsilon$. The noise level set to $\delta \approx 2.561 \times 10^{-4}$[10%].
Figure 3.8 depicts the final reconstructions corresponding to the medium $k_1$. The noise level set to $\delta \approx 2.561 \times 10^{-4}[10\%]$. As previously, the number of adaptations is limited to 50. Under these settings, all the strategies involving the fixed regularisation parameter $\varepsilon$ and the adaptive strategy $\varepsilon_a$ within ASI achieve a successful recovery of $k_1$. However, the reconstructions obtained with adaptive $\varepsilon_a$ (panel (b)) and fixed $\varepsilon = 10^{-6}$ (panel (c)) exhibit higher accuracy than the remaining cases. In particular, ASI with fixed $\varepsilon$ reconstruction smooths out sharp corner features, producing rounded edges in the crescent moon-shaped, whereas the adaptive $\varepsilon_a$ strategy preserves these corners more effectively (panel (b)).

A comparison of the reconstruction performance for the medium $k_1$ under ASI is provided in Table 3.4. Five regularisation parameters, $\varepsilon_a$, $\varepsilon$, $\varepsilon_{\min}$, $\varepsilon_{\text{ave}}$, and $\varepsilon_{\max}$, are evaluated across five key metrics: final dimensionality $N_\theta$, final adaptation count $j$, final misfit $r$, $e_{L^2}$, and $e_{L^\infty}$ errors. The reconstruction using the ASI method with $\varepsilon_a$ achieves the lowest $e_{L^2}$ of 2.81% at $N_\theta = 26$ after $j = 6$ adaptations, with a final misfit of $r = 2.486 \times 10^{-4}$. The reconstruction using ASI with fixed $\varepsilon = 10^{-6}$ follows closely, with the $e_{L^2}$ of 3.02% and the lowest $e_{L^\infty}$ of 11.058, converging at $j = 4$ adaptations. Among other fixed parameters, $\varepsilon_{\max} = 2.753$ produces the highest errors ($e_{L^2} = 3.45\%$ and $e_{L^\infty} = 19.093$). Taken together, these results indicate that ASI with $\varepsilon_a$ achieves a more accurate overall recovery, despite the marginally higher pointwise maximum error.

In Figure 3.9, Panel (a) depicts the evolution of the regularisation parameter $\varepsilon_a$ over the adaptation process. Panel (b) provides a comparative analysis of eigenfunction selection under the adaptive and fixed $\varepsilon$ strategies, thereby highlighting the differences in their respective basis selection mechanisms. Finally, Panel (c) displays the convergence trajectories associated with the recovery of medium $k_1$ in both adaptive and fixed $\varepsilon$ configurations.

### 3.2.2 Medium Example 2 - Five inclusions

Following the approach outlined in the previous example, we initialise the background parameter as $p_0 = 50$. The initial eigenfunction basis $\Phi^{(1)}$ is selected as the first $M = 3$ Laplace eigenfunctions. Given that the underlying problem remains non-linear, the optimisation is performed using the Gauss Newton algorithm in conjunction with Armijo line search. The adaptation process employs a stopping criterion based on the noise level, denoted here as $\delta \approx 1.129 \times 10^{-4}$, which corresponds to 2%.

As previously, Figure 3.10 shows the final recovery of Medium $p = k_2$ using ASI as a regularisation technique with adaptive $\varepsilon_a$, fixed $\varepsilon = 10^{-6}$, $\varepsilon_{\min} = 2.654 \times 10^{-1}$, $\varepsilon_{\text{ave}} = 6.205$, and $\varepsilon_{\max} = 9.889$. Panel (b) in Figure 3.10 shows that ASI using adaptive $\varepsilon_a$ captures the corners of the inclusions more effectively compared with using a fixed $\varepsilon$, whatever values we choose. In particular, the original profile showed in Panel (a) contains two triangles, and the results demonstrate that the ASI method with adaptive $\varepsilon_a$ is able to recover this feature with markedly higher accuracy. This improvement is especially pronounced when contrasted with the fixed $\varepsilon = 10^{-6}$, $\varepsilon_{\min} = 2.654 \times 10^{-1}$, and $\varepsilon_{\text{ave}} = 6.205 \times 10^{-1}$ as illustrated in Panels (c), (d), and (e), respectively.
The quantitative results presented in Table 3.5 are consistent with the visual evidence provided in Figure 3.10. This table provides a rigorous comparison of reconstruction performance for the medium $k_2$ using ASI with both the fixed and adaptive regularisation parameter $\varepsilon$ and a prescribed noise level of $\delta \approx 1.129 \times 10^{-4}$[2%]. The analysis evaluates again five regularisation parameters, $\varepsilon_a$, $\varepsilon = 10^{-6}$, $\varepsilon_{\min} = 2.654 \times 10^{-1}$, $\varepsilon_{\text{ave}} = 6.205$, and $\varepsilon_{\max} = 9.889$, across five diagnostic metrics: final dimension $N_{\theta}$, final adaptation count $j$, misfit residual $r$, and reconstruction error norms ($e_{L^2}$ and $e_{L^{\infty}}$). The results concerning the number of adaptation steps $j$, demonstrate that ASI using both adaptive and fixed strategies of $\varepsilon$ successfully converged within the $j \le 50$ limit. Although ASI with adaptive $\varepsilon_a$ exhibits a slightly higher $e_{L^{\infty}}$ (15.643) compared to ASI using fixed $\varepsilon_{\max} = 9.889$ (15.138) and $\varepsilon_{\text{ave}} = 6.205$ (15.244). It provides the best average accuracy, achieving the lowest $e_{L^2}$ (2.83%) among all other strategies. This indicates that ASI with adaptive $\varepsilon_a$ yields the reconstruction that is globally closest to the true medium $p = k_2$ overall, where the global accuracy and contrast are prioritised over controlling localised peak deviations or sharp edges.

<figure><img src="image_19.png" /><figcaption>Figure 3.10: Medium Example 2 – Final reconstruction of medium $p = k_2$ using ASI with adaptive and fixed $\varepsilon$ alongside the exact medium $k_2$. The corresponding dimensions, values of $\varepsilon$, final misfits, and errors for each panel are provided in Table 3.5. The exponent in brackets gives the number of adaptations needed to convergence. The stopping criterion is based on the level of noise, here, $\delta \approx 1.129 \times 10^{-4}$[2%].</figcaption></figure>

Figure 3.11 illustrates the performance of the ASI for the recovery of the medium $p = k_2$ with a noise level of $\delta \approx 1.129 \times 10^{-4}$[2%]. Panel (a) depicts the evolution of the adaptive $\varepsilon_a$ across
<table><caption>Table 3.5: Medium Example 2 – Recovery of medium $p = k_2$ with ASI. Comparison of reconstruction performance for different parameter selection strategies in terms of final dimension, adaptation count, misfit, and $\delta \approx 1.129 \times 10^{-4}$ [2%].</caption><thead><tr><th></th><th>Value</th><th>Final N<sub>θ</sub></th><th>Final j</th><th>Misfit r</th><th>e<sub>L2</sub></th><th>e<sub>L∞</sub></th></tr></thead><tbody><tr><td>ε<sub>a</sub></td><td>-</td><td>52</td><td>10</td><td>1.114 × 10<sup>-4</sup></td><td>2.83%</td><td>15.643</td></tr><tr><td>ε</td><td>10<sup>-6</sup></td><td>28</td><td>6</td><td>1.123 × 10<sup>-4</sup></td><td>3.20%</td><td>15.673</td></tr><tr><td>ε<sub>min</sub></td><td>2.654 × 10<sup>-1</sup></td><td>29</td><td>14</td><td>1.117 × 10<sup>-4</sup></td><td>3.41%</td><td>31.123</td></tr><tr><td>ε<sub>ave</sub></td><td>6.205</td><td>32</td><td>5</td><td>1.010 × 10<sup>-4</sup></td><td>3.06%</td><td>15.244</td></tr><tr><td>ε<sub>max</sub></td><td>9.889</td><td>63</td><td>6</td><td>1.111 × 10<sup>-4</sup></td><td>3.10%</td><td>15.138</td></tr></tbody></table>

adaptation steps. Panel (b) compares the selected eigenfunctions $N_θ$ between the adaptive $ε_a$ and the fixed approaches $ε_{ave}$, $ε_{max}$ and $ε$, highlighting differences in subspace enrichment and dimensionality. Panel (c) presents the convergence trajectories of the reconstruction misfit $r$ for $k_2$ with both adaptive and fixed $ε$ configurations.

### 3.2.3 Comparison with Tikhonov regularisation

As for the source recovery examples, we compare our ASI approach using adaptive $ε_a$ with the standard Tikhonov regularisation, as in (3.1).

Figure 3.12 presents a visual comparison of the final reconstructions obtained for the medium $p = k_1$ (top row) and the medium $p = k_2$ (bottom row) using two regularisation strategies: ASI with adaptive $ε_a$, and classical Tikhonov regularisation. Panel (a) displays the exact profile $p = k_1$, characterized by a crescent and heart-shaped feature, while Panel (d) shows the true profile of $p = k_2$ including five inclusions. Final reconstructions show that ASI with $ε_a$ in both cases accurately recovers shape, location, and intensity of all geometric components with high spatial resolution – see Panels (b) and (e). In contrast, Panels (c) and (f) exhibit noisy background and the coefficients inside the media are not accurately reconstructed. More precisely, in the five-inclusion case, even though the shape looks good, the value within each feature is not accurate – see, e.g., the triangles and small circle.

Results in Table 3.6 are fully consistent with the visual evidence provided in Figure 3.12. For both $p = k_1$ and $p = k_2$, the ASI method with $ε_a$ outperforms the standard Tikhonov regularisation in both global and peak error metrics. Specifically, ASI with $ε_a$ achieves the $e_{L2}$ of 2.81% for $k_1$ and of 2.83% for $k_2$, compared to 5.06% and 5.47% under Tikhonov regularisation, respectively, demonstrating superior average reconstruction fidelity. The same trend holds for the $e_{L∞}$. Therefore, ASI with $ε_a$ demonstrates a clear advantage over Tikhonov regularisation, achieving a more accurate recovery overall.

Finally, in Figure 3.13 we evaluate the impact of increasing the noise level on the performance of our proposed ASI approach with adaptive $ε_a$ and standard Tikhonov regularisation. For ASI with adaptive $ε_a$ (top row), the recovered medium $p = k_2$ remains visually stable across noise levels 2%, 5% and 10%. At 2%, the ASI method produces a reconstruction that
<figure><img src="image_20.png" /><figcaption>(a) $k_2$: Evolution of adaptive $\epsilon_a$</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(b) $k_2$: Dimension $N_\theta$ vs. iteration $j$</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(c) $k_2$: Misfit $r$ vs. iteration $j$</figcaption></figure>

Figure 3.11: Medium Example 2 - (a) Evolution of $\epsilon_a$ at each adaptation in ASI method. (b) Different selection of eigenfunctions between adaptive and fixed $\epsilon$ settings. (c) Convergence patterns of recovery of Medium $k_2$ with ASI using adaptive and fixed $\epsilon$. The noise level set to $\delta \approx 1.129 \times 10^{-4}$[2%].
<figure><img src="image_23.png" /><figcaption>Figure 3.12: ASI vs Tikhonov regularisation – Final reconstruction of the selected media using ASI with adaptive $\epsilon_a$ and Tikhonov regularisation, alongside the exact media $p = k_1$ and $p = k_2$. The noise level set to $\delta \approx 2.561 \times 10^{-4}$[10%] for recovery of $p = k_1$, and $\delta \approx 1.129 \times 10^{-4}$[2%] for recovery of $p = k_2$, see Table 3.6 for the respective misfits and errors.</figcaption></figure>

<table><caption>Table 3.6: ASI vs Tikhonov regularisation – Comparison of reconstruction performance for different parameter selection strategies in terms of misfit $r$, regularisation parameter $\alpha$ (Tikhonov regularisation only), and errors. The noise level is set to $\delta \approx 2.561 \times 10^{-4}$[10%] for recovery of $p = k_1$, and $\delta \approx 1.129 \times 10^{-4}$[2%] for recovery of $p = k_2$.</caption><thead><tr><th rowspan="2"></th><th colspan="2">p = f<sub>1</sub></th><th colspan="2">p = f<sub>2</sub></th></tr><tr><th>ASI-ε<sub>a</sub></th><th>Tikhonov</th><th>ASI-ε<sub>a</sub></th><th>Tikhonov</th></tr></thead><tbody><tr><th>Misfit r</th><td>2.486 × 10<sup>-4</sup></td><td>2.493 × 10<sup>-4</sup></td><td>1.114 × 10<sup>-4</sup></td><td>1.125 × 10<sup>-4</sup></td></tr><tr><th>Parameter α</th><td>-</td><td>4.302 × 10<sup>-9</sup></td><td>-</td><td>2.493 × 10<sup>-8</sup></td></tr><tr><th>e<sub>L</sub><sup>2</sup></th><td>2.81%</td><td>5.06%</td><td>2.83%</td><td>5.47%</td></tr><tr><th>e<sub>L</sub><sup>∞</sup></th><td>11.362</td><td>15.117</td><td>15.643</td><td>15.747</td></tr></tbody></table>

closely matches the true medium, with consistent contrast across all scatterers. When the noise level increases to 5%, the overall structure of the reconstruction remains stable: the shapes, positions, and intensities of the scatterers are preserved, although slight smoothing and boundary blurring appear at this noise level. Even under the highest tested noise level of 10%, the ASI method with adaptive $\epsilon_a$ continues to capture the fundamental geometric characteristics of all five inclusions. Despite some boundary degradation, the essential geometry, contrast,
and location of all inclusions are preserved, indicating that ASI with adaptive $\epsilon_a$ effectively mitigates the effect of noise. In contrast, the Tikhonov reconstructions (bottom row) exhibit a rapid degradation in the quality of recovered medium $p = k_2$ as noise increases. Even at 2%, Tikhonov reconstructions contain noticeable artefacts. As the noise increases to 5%, these effects intensify: the reconstructed medium becomes increasingly dominated by noise-induced patterns, and the boundaries of the inclusions begin to merge with the background, significantly diminishing the contrast. At the highest noise level 10%, the method almost entirely fails to recover any structural detail. At this noise level, the scatterers become almost indistinguishable from the background. These results collectively highlight the superior noise resilience of our ASI approach with adaptive $\epsilon_a$, underscoring its capability to recover meaningful structural information in scenarios where standard Tikhonov regularisation fails.

<figure><img src="image_24.png" /><figcaption>Figure 3.13: ASI vs Tikhonov regularisation for different noise – Reconstruction of medium $p = k_2$ using the ASI with adaptive $\epsilon_a$ method (top) and standard $L^2$-Tikhonov regularisation (bottom) for different noise levels. The noise level is set to $\delta \approx 1.129 \times 10^{-4}$.</figcaption></figure>

# 4 Conclusion and discussion

We studied more closely different parameters involved in the Adaptive Eigenspace Inversion method, that is, the number of adaptations $j$ performed, the number of eigenfunctions $M = N_\theta$ needed at each adaptation for an optimal reconstruction, and the cancellation parameter $\epsilon_a$ for a better control in the diffusivity. Even though the first two parameters were inspired by the
results of Gleichmann and Grote [13], we first confirm that their approach using sensitivities to build an optimal descent direction in the optimisation algorithm is valid, but we also reduce the computation struggle of merging and orthonormalising the resulting basis by keeping our eigenvalue structure.

For the last parameter, that is, $\epsilon_a$, we conducted a thorough analysis to determine an optimal value. It is obvious that $\epsilon_a$ cannot be larger than the gradient of $p$ to avoid excessive smoothing. However, the actual proportion of the norm of the gradient was unclear. Through this analysis, we found an upper and lower bound for $\epsilon_a$, so that the diffusivity contains the right level of information on $p$ for the next adaptation. An optimal value, as the geometric average of the lower and upper bound, has been identified and successfully tested through multiple numerical experiments.

Current investigations are exploring further the ASI method, in particular, for its stochastic extension using Bayesian inference [16, 17]. One of the consequences is that we rely on the eigenvalue structure, hence our choice to not reshuffle the eigenfunction basis. Adding this probabilistic dimension to the method also involves another parameter, denoted by $\alpha$, in front of a penalty term in the spirit of Tikhonov regularisation. Finding an optimal value for $\alpha$ is our next step to enhance the robustness of ASI while providing uncertainty quantification.

## Acknowledgments

MG and NN were supported by the Marsden Fund from Government funding, Project UOA 2225 "Adaptive Eigenspace Inversion: Analysis, improved efficiency, and extensions", managed by Royal Society Te Apārangi.

## References

[1] D. H. Baffet, Y. G. Gleichmann, and M. J. Grote. Error estimates for adaptive spectral decompositions. *Journal of Scientific Computing*, 93(3):68, 2022.

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

[3] F. Cakoni and D. L. Colton. *Qualitative methods in inverse scattering theory: an introduction*. Interaction of mechanics and mathematics series. Springer, 2006.

[4] M. Cassier and C. Hazard. Multiple scattering of acoustic waves by small sound-soft obstacles in two dimensions: mathematical justification of the foldy-lax model. *Wave Motion*, 50(1):18-28, 2013.

[5] M. Cheney. The linear sampling method and the MUSIC algorithm. *Inverse Problems*, 17(4):591-595, 2001.
[6] D. Colton and A. Kirsch. A simple method for solving inverse scattering problems in the resonance region. *Inverse Problems*, 12(4):383–393, 1996.

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

[8] M. de Buhan and M. 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.

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

[10] S. C. Eisenstat and H. F. Walker. Choosing the forcing terms in an inexact newton method. *SIAM Journal on Scientific Computing*, 17(1):16–32, 1996.

[11] H. Engl, M. Hanke, and A. Neubauer. *Regularization of Inverse Problems*. Mathematics and Its Applications. Springer Netherlands, 1996.

[12] F. Faucher, O. Scherzer, and H. Barucq. Eigenvector models for solving the seismic inverse problem for the Helmholtz equation. *Geophysical Journal International*, 221(1):394–414, 2020.

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

[14] G. H. Golub, P. C. Hansen, and D. P. O'Leary. Tikhonov regularization and total least squares. *SIAM Journal on Matrix Analysis and Applications*, 21(1):185–194, 1999.

[15] M. Graff, M. J. Grote, F. Nataf, and F. Assous. How to solve inverse scattering problems without knowing the source term: a three-step strategy. *Inverse Problems*, 35(10):104001, 2019.

[16] M. Graff and B. 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.

[17] M. Graff, B. Hosseini, and A. Kim. Bayesian Adaptive Eigenspace Inversion. working paper or preprint, June 2025.

[18] M. J. Grote, J. Huber, D. Kourounis, and O. Schenk. Inexact interior-point method for pde-constrained nonlinear optimization. *SIAM Journal on Scientific Computing*, 36(3):A1251–A1276, 2014.

[19] M. J. Grote, M. Kray, and U. Nahum. Adaptive eigenspace method for inverse scattering problems in the frequency domain. *Inverse Problems*, 33(2):025006, 2017.
[20] M. J. Grote and U. Nahum. Adaptive eigenspace for multi-parameter inverse scattering problems. *Computers & Mathematics with Applications*, 77(12):3264–3280, June 2019.

[21] E. Haber, U. M. Ascher, and D. Oldenburg. On optimization techniques for solving non-linear inverse problems. *Inverse Problems*, 16:1263, 2000.

[22] E. Haber and M. Chung. Simultaneous source for non-uniform data variance and missing data. available on arXiv:1404.5254, submitted in 2014.

[23] H. Haddar, A. Lechleiter, and S. Marmorat. An improved time domain linear sampling method for Robin and Neumann obstacles. *Applicable Analysis*, 0:1–22, 2013.

[24] P. C. Hansen. *Discrete Inverse Problems*. Society for Industrial and Applied Mathematics, 2010.

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

[26] A. Kirsch. *An Introduction to the Mathematical Theory of Inverse Problems*. Applied Mathematical Sciences. Springer New York, 1996.

[27] A. Kirsch. The MUSIC-algorithm and the factorization method in inverse scattering theory for inhomogeneous media. *Inverse Problems*, 18(4):1025–1040, 2002.

[28] A. Kirsch and N. Grinberg. *The Factorization Method for Inverse Problems*. Oxford Lecture Series in Mathematics and Its Applications. OUP Oxford, 2008.

[29] A. Lechleiter. The factorization method is independent of transmission eigenvalues. *Inverse Problems and Imaging*, 3(1):123–138, 2009.

[30] F. Lucka, N. Huynh, M. Betcke, E. Zhang, P. Beard, B. Cox, and S. Arridge. Enhancing compressed sensing 4D photoacoustic tomography by simultaneous motion estimation. *SIAM J. Imaging Sciences*, 11(4):2224–2253, 2018.

[31] A. Mang, A. Gholami, C. Davatzikos, and G. Biros. PDE-constrained optimization in medical image analysis. *Optimization and Engineering*, 19(3):765–812, 2018.

[32] L. Métivier, R. Brossier, J. Virieux, and S. Operto. Full Waveform Inversion and the Truncated Newton Method. *SIAM J. Sci. Comput.*, 35(2):B401–B437, 2013.

[33] P. Mojabi and J. LoVetri. Microwave Biomedical Imaging Using the Multiplicative Regularized Gauss-Newton Inversion. *Antennas and Wireless Propagation Letters*, IEEE, 8:645–648, 2009.

[34] R. Potthast. A survey on sampling and probe methods for inverse problems. *Inverse Problems*, 22(2):R1, 2006.
[35] C. Prada, S. Manneville, D. Spoliansky, and M. Fink. Decomposition of the Time Reversal Operator: Application to detection and selective focusing on two scatterers. *J. Acoust. Soc. Am.*, 99(4):2067–2076, 1996.

[36] R. G. Pratt. Seismic waveform inversion in the frequency domain; part 1, theory and verification in a physical scale model. *Geophysics*, 64(3):888–901, 1999.

[37] R. G. Pratt, C. Shin, and G. J. Hick. Gauss–newton and full newton methods in frequency–space seismic waveform inversion. *Geophysical Journal International*, 133(2):341–362, 1998.

[38] L. I. Rudin, S. Osher, and E. Fatemi. Nonlinear total variation based noise removal algorithms. *Physica D: Nonlinear Phenomena*, 60(1):259–268, 1992.

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

[40] A. Tarantola. Inversion of seismic reflection data in the acoustic approximation. *Geophysics*, 49(8):1259–1266, 1984.

[41] A. Tarantola. *Inverse Problem Theory And Methods For Model Parameter Estimation*. Society for Industrial and Applied Mathematics, 2005.

[42] A. Tikhonov and V. Arsenin. *Solutions of Ill-posed Problems*. Halsted Press book. Winston, 1977.

[43] A. N. Tikhonov. On the stability of inverse problems. *Dokl. Akad. Nauk SSSR*, 39(5):195–198, 1943.

[44] J. Virieux and S. Operto. An overview of full-waveform inversion in exploration geophysics. *Geophysics*, 74(6):WCC1–WCC26, 2009.

[45] C. R. Vogel and M. E. Oman. Iterative methods for total variation denoising. *SIAM Journal on Scientific Computing*, 17(1):227–238, 1996.