# Adaptive Spectral Inversion for Inverse Medium Problems

Yannik G. Gleichmann 111Department of Mathematics and Computer Science,
University of Basel, Basel, Switzerland<EMAIL_ADDRESS>Marcus J.
Grote 222Department of Mathematics and Computer Science, University of Basel,
Basel, Switzerland; Corresponding author<EMAIL_ADDRESS>

###### Abstract

A nonlinear optimization method is proposed for the solution of inverse medium
problems with spatially varying properties. To avoid the prohibitively large
number of unknown control variables resulting from standard grid-based
representations, the misfit is instead minimized in a small subspace spanned
by the first few eigenfunctions of a judicious elliptic operator, which itself
depends on the previous iteration. By repeatedly adapting both the dimension
and the basis of the search space, regularization is inherently incorporated
at each iteration without the need for extra Tikhonov penalization.
Convergence is proved under an angle condition, which is included into the
resulting _Adaptive Spectral Inversion_ (ASI) algorithm. The ASI approach
compares favorably to standard grid-based inversion using $L^{2}$-Tikhonov
regularization when applied to an elliptic inverse problem. The improved
accuracy resulting from the newly included angle condition is further
demonstrated via numerical experiments from time-dependent inverse scattering
problems.

Keywords: Inverse problem, inverse scattering problem, adaptive eigenspace
inversion, nonlinear optimization

##  1 Introduction

The solution of inverse medium problems entails the reconstruction of a
medium’s spatially varying properties, $u(x)$, inside a bounded region of
space $\Omega\subset\mathbb{R}^{d}$ from partially available, often noisy
observations of a state variable $y$. This is inverse to the direct or forward
problem of determining $y$ for a given medium $u(x)$, where typically $y$
satisfies a governing partial differential equation, be it stationary or time-
dependent, which involves $u(x)$. More abstractly, if we denote the solution
to the forward problem by a forward operator $F:H_{1}\to H_{2}$ acting between
two Hilbert spaces $H_{1},H_{2}$, the associated inverse problem then consists
in solving

$\displaystyle F(u)=y,$ (1.1)

for a given right-hand side $y^{\dagger}\in H_{2}$, that is, to determine a
solution $u^{\dagger}\in H_{1}$ such that $y^{\dagger}=F(u^{\dagger})$.

Due to measurement errors, only approximate (noisy) data $y^{\delta}$ is
available in practice,

$\displaystyle\|y^{\dagger}-y^{\delta}\|_{H_{2}}\leqslant\delta,$ (1.2)

for a noise level $\delta\geqslant 0$. Then, the inverse medium problem
consists in finding $u^{\dagger,\delta}\in H_{1}$ that satisfies

$\displaystyle F(u)=y^{\delta},$ (1.3)

which we reformulate as a least-squares minimization problem for the misfit
$\operatorname{\mathcal{J}}^{\delta}:H_{1}\to\mathbb{R}$ as

$\displaystyle u^{\dagger,\delta}=\operatorname*{argmin}_{u\in
H_{1}}\operatorname{\mathcal{J}}^{\delta}(u):=\operatorname*{argmin}_{u\in
H_{1}}\frac{1}{2}\|F(u)-y^{\delta}\|^{2}_{H_{2}}.$ (1.4)

In general, (1.3) or (1.4) is ill-posed and cannot be solved directly without
further regularization. Numerical methods for obtaining stable solutions to
(1.3) or (1.4) essentially fall into three categories. First, starting from an
initial guess $u^{(0)}$, _iterative regularization methods_ , such as the
Landweber iteration [21, 28], a special case of gradient descent, improve upon
$u^{(m)}$ at iteration $m$ by using the derivative of
$\operatorname{\mathcal{J}}^{\delta}$ [3, 25]. Due to semi-convergence, the
iteration needs to be stopped judiciously as the regularization parameter
corresponds to the reciprocal of $m$. Second, _Tikhonov regularization_ adds a
regularization term $R(u)$ to $\operatorname{\mathcal{J}}^{\delta}(u)$ [12,
27]. Hence, instead of solving (1.4), one now seeks a minimizer of
$\operatorname{\mathcal{J}}^{\delta}(u)+\alpha R(u)$, where $\alpha$ is the
regularization parameter which must be chosen carefully. Typical choices for
the regularization term $R(u)$ include the $L^{2}$-norm, the
$\operatorname{TV}$-norm, or the $H^{1}$-semi-norm, depending on the expected
smoothness in $u^{\dagger}$. Third, _regularization by projection_ _(or
discretization)_ in preimage-space replaces in (1.4) or (1.3) the infinite
dimensional space $H_{1}$ by a finite dimensional subspace $\Psi^{(m)}\subset
H_{1}$ and then solves for $u^{(m)}\in\Psi^{(m)}$ [16, 26]. Here, $\Psi^{(m)}$
typically corresponds to the underlying subspace of a finite difference or
finite element discretization whose mesh-size, $h_{m}>0$, then acts as the
regularization parameter. Alternative projection methods use finite-
dimensional subspaces in image-space [23, 24], or even projection into both
image- and preimage-space simultaneously [20, 33].

Unless an effective parametrization of $u(x)$ is known a priori, inverse
medium problems usually lead to prohibitively high-dimensional search spaces
where the number of unknown parameters (or control variables) is determined by
the degrees of freedom of the underlying finite difference or finite element
discretization. Alternatively, sparsity promoting strategies attempt to remain
sufficiently general while keeping the dimension of the search space small by
applying $\ell_{1}$-norm soft-thresholding, say, to promote a sparse
representation in a fixed wavelet, curvelet, etc. basis or frame [8, 22, 29,
31].

An adaptive inversion method was first introduced in [11, Section 5.4], where
the search space consists of the first few eigenfunctions of a judicious
elliptic differential operator $L[u^{(m)}]$. As $L[u^{(m)}]$ itself depends on
the previous iterate, $u^{(m)}$, so do its eigenfunctions; hence, the current
search space is repeatedly adapted at every iteration. To avoid division by
zero in the presence of vanishing gradients, De Buhan and Kray later
incorporated a small fixed parameter $\varepsilon$ which led to the coercive
elliptic operator $L_{\varepsilon}[u^{(m)}]$. Restricting the search space to
the first few eigenfunctions of $L_{\varepsilon}[u^{(m)}]$ has proved highly
effective in a number of acoustic [10, 15], electro-magnetic [9] and seismic
[13] inverse scattering problems, but also in optimal control [34].

To take advantage of the regularizing effect of a low-dimensional search
space, both the eigenfunctions and the dimension of the search space were
adapted in [17]. Moreover, the connection between $L_{\varepsilon}$ and TV-
regularization motivated the use of different elliptic operators depending on
the expected smoothness in the target medium [18]. Recently, a first step was
taken in developing a mathematical theory underpinning the remarkable accuracy
of that decomposition for the approximation of piecewise constant functions
[2], which led to rigorous $L^{2}$-error estimates in [1].

The remaining part of this paper is structured as follows. In Section 2, we
introduce the general _Adaptive Inversion_ iteration for the solution of
inverse problems, which proceeds by solving (1.4) in a sequence of finite
dimensional subspaces $\Psi^{(m)}\subset H_{1}$, not necessarily known a
priori. Here we identify a key angle condition, which yields convergence of
the Adaptive Inversion iteration, and also prove that it is a genuine
regularization method. In Section 3, we present the _Adaptive Spectral (AS)
decomposition_ and recall its approximation properties [1, 2]. By combining
the Adaptive Inversion iteration together with the AS decomposition, we
propose in Section 4 the _Adaptive Spectral Inversion (ASI)_ Algorithm, which
incorporates the new angle condition from Section 2 by including the
sensitivities of the gradient into the construction of the search space at the
subsequent iteration. Finally in Section 5, we present several numerical
experiments which illustrate the performance of the ASI method and verify the
theory from Section 2. In particular, we apply the ASI method to a time-
dependent inverse scattering problems which demonstrates that it yields a
significant improvement over previous versions from [17, 2].

##  2 Adaptive Inversion

To determine a solution to the inverse problem (1.3) for given data
$y^{\delta}$, we shall minimize the misfit
$\operatorname{\mathcal{J}}^{\delta}$ in (1.4) successively in a sequence of
(closed) finite-dimensional subspaces $\Psi^{(m)}\subset H_{1}$, $m\geqslant
1$, until its gradient $D\operatorname{\mathcal{J}}^{\delta}$ vanishes. In
doing so, we do not assume the entire sequence of subspaces
$\\{\Psi^{(m)}\\}_{m\geqslant 1}$ known a priori. Hence, we consider the
following adaptive inversion algorithm:

Algorithm 1 Adaptive Inversion

0. initial guess $u^{(0),\delta}$, search space $\Psi^{(1)}\subset H_{1}$, $m=1$.

1. while $\|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\|\neq 0$ do

2.  _Solve_ $\displaystyle u^{(m),\delta}=\operatorname*{argmin}_{u\in\Psi^{(m)}}\operatorname{\mathcal{J}}^{\delta}(u).$ (2.1)

3.  _Determine_ new search space $\Psi^{(m+1)}$, such that $u^{(m),\delta}\in\Psi^{(m+1)}$.

4. $m\leftarrow m+1$

By ensuring that $u^{(m),\delta}$ also belongs to the new search space
$\Psi^{(m+1)}$, we guarantee in every iteration that the misfit does not
increase,

$\displaystyle\operatorname{\mathcal{J}}^{\delta}(u^{(m+1),\delta})\leqslant\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta}),\qquad\forall
m\geqslant 1.$ (2.2)

Therefore the sequence $\\{u^{(m),\delta}\\}_{m\geqslant 1}$ obtained by the
Adaptive Inversion Algorithm 1 is a minimizing sequence of
$\operatorname{\mathcal{J}}^{\delta}$. Without further assumptions, however,
we do not know yet whether this sequence converges.

###  2.1 Convergence

To prove that the sequence $\\{u^{(m),\delta}\\}_{m\geqslant 1}$ generated by
the Adaptive Inversion Algorithm indeed converges, we proceed as follows.
First, we identify a key angle condition, which ensures that for fixed
$\delta\geqslant 0$,
$\|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\|\to 0$ as
$m\to\infty$. Then, under suitable assumptions from convex optimization
theory, we conclude that the sequence indeed converges to
$u^{\dagger,\delta}$.

######  Theorem 2.1.

Let the misfit $\operatorname{\mathcal{J}}^{\delta}$ be Fréchet differentiable
with derivative $D\operatorname{\mathcal{J}}^{\delta}(u)$, $u\in H_{1}$ and
Lipschitz-continuous for every direction $d\in H_{1}$, i.e.

$\displaystyle|D\operatorname{\mathcal{J}}^{\delta}(u+v)d-D\operatorname{\mathcal{J}}^{\delta}(u)d|\leqslant
L\|v\|_{H_{1}}\|d\|_{H_{1}}\qquad\forall u,v\in H_{1},\quad L>0,$ (2.3)

and further assume that $\operatorname{\mathcal{J}}^{\delta}(u)\geqslant
C>-\infty$ for all $u\in H_{1}$. Then, if the corrections
$d^{(m),\delta}=u^{(m+1),\delta}-u^{(m),\delta}$ satisfy the _angle condition_

$\displaystyle|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})d^{(m),\delta}|\geqslant{\varepsilon_{\theta}}\|d^{(m),\delta}\|_{H_{1}}\|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\|,\qquad
0<{\varepsilon_{\theta}}<1$ (2.4)

uniformly in $m$, we have

$\displaystyle\|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\|\to
0,\qquad m\to\infty.$ (2.5)

###### Proof.

First, we claim that

$\displaystyle\operatorname{\mathcal{J}}^{\delta}(u^{(m+1),\delta})-\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\leqslant-\frac{\alpha_{m}}{2}\mu_{m}D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})d^{(m),\delta},$
(2.6)

where

$\displaystyle\alpha_{m}=\frac{\mu_{m}}{L}\frac{D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})d^{(m),\delta}}{\|d^{(m),\delta}\|^{2}}>0,\qquad\mu_{m}=\operatorname{sign}(D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})d^{(m),\delta}).$
(2.7)

As both $u^{(m+1),\delta},u^{(m),\delta}\in\Psi^{(m+1)}$, we also have
$d^{(m),\delta}\in\Psi^{(m+1)}$. Since
$\operatorname{\mathcal{J}}^{\delta}({u}^{(m+1),\delta})\leqslant\operatorname{\mathcal{J}}^{\delta}(u)$
for all $u\in\Psi^{(m+1)}$, we may choose
$u=u^{(m),\delta}-\alpha_{m}\mu_{m}d^{(m),\delta}$. The linearity of the
derivative thus yields

$\displaystyle\operatorname{\mathcal{J}}^{\delta}(u^{(m+1),\delta})-\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})+\alpha_{m}\mu_{m}D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})d^{(m),\delta}$
(2.8) $\displaystyle\hskip
20.00003pt\leqslant\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta}-\alpha_{m}\mu_{m}d^{(m),\delta})-\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})+\alpha_{m}\mu_{m}D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})d^{(m),\delta}$
(2.9) $\displaystyle\hskip
20.00003pt=\alpha_{m}\mu_{m}\int_{0}^{1}\left(D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})-D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta}-\tau\alpha_{m}\mu_{m}d^{(m),\delta})\right)d^{(m),\delta}d\tau$
(2.10) $\displaystyle\hskip
20.00003pt\leqslant\frac{\alpha_{m}L}{2}\|d^{(m),\delta}\|^{2}\alpha_{m}=\frac{\alpha_{m}}{2}\mu_{m}D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})d^{(m),\delta},$
(2.11)

which proves (2.6).

Since $\operatorname{\mathcal{J}}^{\delta}$ is bounded from below and
decreases in every iteration, there exists a constant $C_{0}\geqslant 0$ such
that by (2.6)

$\displaystyle-C_{0}$ $\displaystyle\leqslant
C-\operatorname{\mathcal{J}}^{\delta}(u^{(1),\delta})\leqslant\operatorname{\mathcal{J}}^{\delta}(u^{(M+1),\delta})-\operatorname{\mathcal{J}}^{\delta}(u^{(1),\delta})$
(2.12)
$\displaystyle=\sum_{m=1}^{M}\left(\operatorname{\mathcal{J}}^{\delta}(u^{(m+1),\delta})-\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\right)\leqslant-\sum_{m=1}^{M}\frac{\alpha_{m}}{2}\mu_{m}D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})d^{(m),\delta}$
(2.13)
$\displaystyle=-\frac{1}{2L}\sum_{m=1}^{M}\frac{(D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta}))^{2}}{\|d^{(m),\delta}\|^{2}}\leqslant
0.$ (2.14)

Next, we define the angle $\theta_{m}$ via

$\displaystyle\cos(\theta_{m})=\frac{\mu_{m}D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})d^{(m),\delta}}{\|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\|\|d^{(m),\delta}\|}>0$
(2.15)

and rewrite (2.12) – (2.14) using (2.7) as

$\displaystyle-
C_{0}\leqslant-\frac{1}{2L}\sum_{m=1}^{M}\cos^{2}(\theta_{m})\|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\|^{2}\leqslant
0.$ (2.16)

Taking the limit $M\rightarrow\infty$ then yields the well known _Zoutendijk
condition_ [36]

$\displaystyle\sum_{m=1}^{\infty}\cos^{2}(\theta_{m})\|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\|^{2}<\infty,$
(2.17)

which immediately implies

$\displaystyle\cos^{2}(\theta_{m})\|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\|^{2}\to
0,\qquad m\to\infty.$ (2.18)

Since we optimize for $u^{(m),\delta}$ such that the angle condition (2.4)
holds uniformly in $m$, i.e. $\cos(\theta_{m})\geqslant{\varepsilon_{\theta}}$
for all $m\geqslant 1$, with $0<{\varepsilon_{\theta}}<1$, we thus conclude
that

$\displaystyle\|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\|\to
0,\qquad m\to\infty.$ (2.19)

∎

Theorem 2.1 implies that the Adaptive Inversion Algorithm yields a minimizing
sequence $\\{u^{(m),\delta}\\}_{m\geqslant 1}$ such that the gradient of the
misfit tends to zero; hence, the algorithm is well-defined and converges.
Without further assumptions, however, this does not imply that the sequence
$\\{u^{(m),\delta}\\}_{m\geqslant 1}$ converges (weakly) to a (local)
minimizer or accumulation point. Under further standard assumptions from
convex optimization theory, one can even prove convergence to a minimizer of
(1.4), as in [4, Corollary 11.30]. Those assumptions, however, seldom hold in
practice for inverse medium problems.

###  2.2 Regularization

In the presence of perturbed noisy data, it makes little sense to improve the
approximate solution $u^{(m),\delta}$ at iteration $m$ beyond the error
$\delta$ in the observations. Instead for $\tau>1$, we stop the Adaptive
Inversion Algorithm at iteration $m_{\ast}(\delta)$ when the _discrepancy
principle_ is satisfied:

$\displaystyle
m_{\ast}(\delta)=m_{\ast}=\min\\{m\in\mathbb{N}:\;\|F(u^{(m),\delta})-y^{\delta}\|_{H_{2}}\leqslant\tau\delta\\}.$
(2.20)

If we assume that the search spaces $\Psi^{(m)}$ satisfy

$\displaystyle\|(I-\Pi_{\Psi^{(m)}})u^{\dagger}\|_{H_{1}}\to 0,\qquad
m\to\infty,$ (2.21)

where $\Pi_{\Psi^{(m)}}$ denotes the projection into $\Psi^{(m)}$ and
$u^{\dagger}$ the exact (noise-free) solution, the following lemma guarantees
that the Adaptive Inversion Algorithm will always satisfy (2.20) after a
finite number of steps.

######  Lemma 2.2.

Suppose that the forward operator $F$ is continuous and that the spaces
$\\{\Psi^{(m)}\\}_{m\geqslant 1}$ satisfy (2.21). Then, for $\tau>1$, there
exists for every $\delta\geqslant 0$ an index $m_{\ast}(\delta)$ such that the
discrepancy principle (2.20) is satisfied.

###### Proof.

By minimality of $u^{(m),\delta}$, the residual is bounded by

$\displaystyle\|F(u^{(m),\delta})-y^{\delta}\|_{H_{2}}$
$\displaystyle\leqslant\|F(\Pi_{\Psi^{(m)}}u^{\dagger})-y^{\delta}\|_{H_{2}}$
(2.22)
$\displaystyle\leqslant\|F(\Pi_{\Psi^{(m)}}u^{\dagger})-F(u^{\dagger})\|_{H_{2}}+\|F(u^{\dagger})-y^{\delta}\|_{H_{2}}$
(2.23)
$\displaystyle\leqslant\|F(\Pi_{\Psi^{(m)}}u^{\dagger})-F(u^{\dagger})\|_{H_{2}}+\delta.$
(2.24)

From the continuity of $F$ and (2.21) we obtain

$\displaystyle\|F(u^{(m),\delta})-y^{\delta}\|_{H_{2}}\leqslant\tau\delta,\qquad\tau>1,$
(2.25)

for $m\geqslant M$ sufficiently large. Thus $m_{\ast}(\delta):=M$ is always
well defined. ∎

Next, we consider the sequence $\\{u^{(m_{\ast}),\delta}\\}_{\delta\geqslant
0}$ obtained from the Adaptive Inversion Algorithm when stopped via the
discrepancy principle (2.20). To show that it converges to the exact (noise-
free) solution $u^{\dagger}$ as $\delta\to 0$, that is that

$\displaystyle u^{(m_{\ast}),\delta}\to u^{\dagger},\qquad\delta\to 0,$ (2.26)

we assume that the forward operator $F$ is continuous with Fréchet derivative
$DF(u)$ and that the standard _tangential cone condition_ (aka _Scherzer
condition_),

$\displaystyle\|F(u)-F(v)-DF(u)(u-v)\|_{H_{2}}\leqslant\eta\|F(u)-F(v)\|_{H_{2}},\qquad\forall
u,v\in H_{1},$ (2.27)

holds for some $\eta\in(0,1)$. Moreover, we assume that $F$ is weakly
sequentially closed.

Under all the above assumptions, we conclude from [26, Theorem 3.4] that the
Adaptive Inversion Algorithm is a genuine regularization method in the sense
of [12, Definition 3.1]:

######  Theorem 2.3.

Let $\Psi^{(m)}$ satisfy (2.21), $F$ be weakly sequentially closed, Fréchet
differentiable and also satisfy (2.27). Then, the sequence
$\\{u^{(m_{\ast}),\delta}\\}_{\delta\geqslant 0}$ obtained from the Adaptive
Inversion Algorithm 1, stopped at iteration $m_{\ast}=m_{\ast}(\delta)$
according to (2.20), admits a subsequence converging to a minimizer of
$\operatorname{\mathcal{J}}^{\delta}$. Moreover, if the minimizer is unique,
then the sequence $\\{u^{(m_{\ast}),\delta}\\}_{\delta\geqslant 0}$ converges
to the exact (noise-free) solution $u^{\dagger}$ of the inverse problem (1.3),
that is

$\displaystyle u^{(m_{\ast}),\delta}\to u^{\dagger},\qquad\delta\to 0.$ (2.28)

######  Remark 1.

If the nullspace of the Fréchet derivative $DF(u)$ is trivial, the Scherzer
condition (2.27) immediately implies the uniqueness of $u^{\dagger}$. Indeed,
from (2.27) and the triangle inequality, we first deduce that

$\displaystyle\|DF(u)(u-v)\|\leqslant(1+\eta)\|F(u)-F(v)\|\qquad\forall u,v\in
H_{1}.$ (2.29)

Given two distinct solutions $u^{\dagger,\delta}\neq v^{\dagger,\delta}$ to
(1.3), we then infer that

$\displaystyle
0=\|F(u^{\dagger,\delta})-F(v^{\dagger,\delta})\|\geqslant\frac{1}{\eta+1}\|DF(u^{\dagger,\delta})(u^{\dagger,\delta}-v^{\dagger,\delta})\|,$
(2.30)

which yields the sought contradiction and hence the uniqueness of
$u^{\dagger}$. Thus, if the nullspace of the Fréchet derivative $DF(u)$ is
trivial $\forall u\in H_{1}$, then Theorem 2.3 immediately implies the
convergence of $\\{u^{(m_{\ast}),\delta}\\}_{\delta\geqslant 0}$ to the exact
solution. Clearly, those assumptions may not hold in practice, for instance,
in a situation of limited indirect observations.

In summary, the Adaptive Inversion Algorithm generates a minimizing sequence
$\\{u^{(m),\delta}\\}_{m\geqslant 1}$ for
$\operatorname{\mathcal{J}}^{\delta}$ and, under standard assumptions, Theorem
2.1 implies that the gradient of $\operatorname{\mathcal{J}}^{\delta}$ tends
to zero; hence, the algorithm is well defined and terminates. In fact, under
standard assumptions from convex optimization theory, there exists for every
$\delta$ a minimizer $u^{\dagger,\delta}$ with $u^{(m),\delta}\to
u^{\dagger,\delta}$ as $m\to\infty$. Moreover, under standard assumptions for
inverse problems, we deduce that the Adaptive Inversion Algorithm yields a
genuine regularization method; thus, $u^{(m_{\ast}),\delta}$ converges to the
exact (noise-free) solution $u^{\dagger}$ of (1.3) as $\delta\to 0$.

##  3 Adaptive Spectral Decomposition

To apply the Adaptive Inversion Algorithm from Section 2, we must specify the
search spaces $\Psi^{(m)}$ in practice. Here we shall consider low-dimensional
search spaces spanned by the first few eigenfunctions of a judicious elliptic
operator, which itself depends on the previous iterate. By combining those
search spaces with the Adaptive Inversion Algorithm, we shall devise in
Section 4 the Adaptive Spectral Inversion Algorithm for the solution of
inverse medium problems.

Consider an arbitrary bounded function $u:\Omega\to\mathbb{R}$ on a bounded
Lipschitz domain $\Omega\subset\mathbb{R}^{d}$, $d\geqslant 2$, which vanishes
on the boundary. Next, let $V_{h}\subset H_{0}^{1}(\Omega)$ be the finite
element space of continuous, piecewise polynomials of degree $r\geqslant 1$,
on a regular and quasi-uniform mesh with mesh-size $h>0$. Denoting by
$u_{h}\in V_{h}$ the standard $\mathbb{P}^{r}$-FE interpolant of $u$, we
introduce the differential operator

$\displaystyle
L_{\varepsilon}[u_{h}]v=-\nabla\cdot\left(\mu_{\varepsilon}[u_{h}]\nabla
v\right),\quad\mu_{\varepsilon}[u_{h}]=\frac{1}{\sqrt{|\nabla
u_{h}|^{2}+\varepsilon^{2}}},\qquad\varepsilon>0.$ (3.1)

In contrast to $u$, which may be discontinuous, $u_{h}$ is continuous and
piecewise polynomial with $\nabla u_{h}\in L^{\infty}(\Omega)$. For
$\varepsilon$ small, $L_{\varepsilon}[u_{h}]$ in (3.1) is well defined and for
$h$ small even uniformly elliptic in $\Omega$ [1]; in practice, we always set
$\varepsilon=10^{-8}$. Hence, there exists a non-decreasing sequence of
strictly positive eigenvalues $\\{\lambda_{k}\\}_{k\geqslant 1}$ with
corresponding eigenfunctions $\\{\varphi_{k}\\}_{k\geqslant 1}$ satisfying

$\displaystyle
L_{\varepsilon}[u_{h}]\varphi_{k}=\lambda_{k}\varphi_{k}\quad\text{in}\;\Omega,\qquad\varphi_{k}=0\quad\text{on}\;\partial\Omega.$
(3.2)

Moreover, the eigenfunctions of $L_{\varepsilon}[u_{h}]$ form an
$L^{2}$-orthonormal basis of $L^{2}(\Omega)$.

(a) $u$ (or $u_{h}$)

(b) $\Pi_{K}[u_{h}]u$ for $K=3$

Figure 1: AS decomposition: A piecewise constant function
$u:\Omega\to\mathbb{R}$ and its projection $\Pi_{K}[u_{h}]u$ onto the first
$K=3$ eigenfunctions from (3.2).

(a) $\varphi_{1}$, $\lambda_{1}\approx 1.5$

(b) $\varphi_{2}$, $\lambda_{2}\approx 4.7$

(c) $\varphi_{3}$, $\lambda_{3}\approx 5.3$

(d) $\varphi_{4}$, $\lambda_{4}\approx 1.0\cdot 10^{8}$

(e) $\varphi_{5}$, $\lambda_{5}\approx 1.1\cdot 10^{8}$

(f) $\varphi_{6}$, $\lambda_{6}\approx 1.3\cdot 10^{8}$

Figure 2: AS decomposition: The first six eigenfunctions $\varphi_{k}$,
$k=1,\ldots,6$ of $L_{\varepsilon}[u_{h}]$, with its corresponding eigenvalues
$\lambda_{k}$ satisfying (3.2) for $u_{h}$ as in Figure 1(a).

For $K\geqslant 1$, we let
$\Phi_{K}=\operatorname*{span}\\{\varphi_{1},\ldots,\varphi_{K}\\}$ and denote
by $\Pi_{K}[u_{h}]$ the $L^{2}$-projection into $\Phi_{K}$:

$\displaystyle\Pi_{K}[u_{h}]:L^{2}(\Omega)\to\mathbb{R},\qquad(v-\Pi_{K}[u_{h}]v,\varphi)_{L^{2}(\Omega)}=0\qquad\forall\varphi\in\Phi_{K},\forall
v\in L^{2}(\Omega).$ (3.3)

We call

$\displaystyle\Pi_{K}[u_{h}]u=\sum_{k=1}^{K}\beta_{k}\varphi_{k},\qquad\beta_{k}\in\mathbb{R}$
(3.4)

the (truncated) _Adaptive Spectral (AS) decomposition_ of $u$. Note that the
operator $L_{\varepsilon}[u_{h}]$ depends on $u_{h}$ (and $\varepsilon$), and
hence so do its eigenfunctions as well as $\Pi_{K}[u_{h}]$.

To illustrate the AS decomposition, consider the piecewise constant function
$u:\Omega\subset\mathbb{R}^{2}\to\mathbb{R}$ shown in Figure 1(a). Next, we
(numerically) solve the eigenvalue problem (3.2) with $\varepsilon=10^{-8}$
using a $\mathbb{P}^{1}$-FE discretization with mesh-size $h=0.01$. As shown
in Figure 2, the first three eigenfunctions $\varphi_{k}$, $k=1,2,3$ capture
the inclusions rather well, whereas the subsequent eigenfunctions resemble
eigenfunctions of the Laplacian with eigenvalues that scale as
$1/\varepsilon$. Next, we project $u$ into
$\Phi_{3}=\operatorname*{span}\\{\varphi_{1},\varphi_{2},\varphi_{3}\\}$. As
shown in Figure 1(b), the projection $\Pi_{3}[u_{h}]u$ is hardly
distinguishable from $u$.

In [1, 2], rigorous $L^{2}$-error bounds were derived for the AS
decomposition. More precisely, consider a piecewise constant function with
zero boundary values of the form

$\displaystyle
u=\sum_{k=1}^{K}\alpha_{k}\chi_{A_{k}},\qquad\alpha_{k}\in\mathbb{R}\smallsetminus\\{0\\},$
(3.5)

where $\chi_{A_{k}}$ is the characteristic function of Lipschitz domains
$A_{k}\subset\subset\Omega$ with mutually disjoint and connected boundaries
$\partial A_{k}$.

Next we partition $\Omega=\overline{M_{h}\cup D_{h}}$ into the two open sets

$\displaystyle
M_{h}=\bigcup_{k=1}^{K}\left\\{x\in\Omega:\operatorname{dist}(x,\partial
A_{k})<h\right\\},\qquad D_{h}=\Omega\smallsetminus\overline{M_{h}},$ (3.6)

where $M_{h}$ is the open tube with diameter $h>0$ around the jump
discontinuities of $u$. In Figure 3, for instance, the sets $A_{k}$, $M_{h}$,
and $D_{h}$ are shown for $u$ as in Figure 1(a).

Figure 3: AS decomposition: The sets $A_{k}$, $k=1,2,3$, $M_{h}$, and $D_{h}$
from (3.6) for $u$ as in Figure 1(b).

Clearly, the FE-interpolant $u_{h}\in V_{h}$ of $u$ satisfies $u_{h}=u$ and
$\nabla u_{h}=0$ in $D_{h}$. Moreover, for a sequence of regular and quasi-
uniform meshes $(\mathcal{T}_{h})_{h>0}$, we have

$\displaystyle\lim_{h\to 0}\|u_{h}-u\|_{L^{2}(\Omega)}=0,\qquad\nabla u_{h}\in
L^{\infty},\qquad\operatorname{supp}(\nabla u_{h})\subset\overline{M_{h}},$
(3.7)

and $h\|u_{h}\|_{L^{\infty}}\leqslant C$ with $C>0$ independent of $h$, [1,
Proposition 2.2]. In [2], those properties led to the following upper bound
for the first $K$ eigenfunctions.

######  Theorem 3.1.

Let $u$ be given by (3.5), $u_{h}\in V_{h}$ be its FE interpolant, and
$\varphi_{k}$, $k\geqslant 1$ the eigenfunctions of the AS operator
$L_{\varepsilon}[u_{h}]$. Then, for every $\varepsilon>0$ and $h>0$
sufficiently small, there exists a constant $C>0$, independent of
$\varepsilon$ and $h$, such that

$\displaystyle\|\nabla\varphi_{k}\|_{L^{2}(D_{h})}\leqslant
C\sqrt{\varepsilon},\qquad k=1,\ldots,K.$ (3.8)

From Theorem 3.1 we conclude that the first $K$ eigenfunctions of
$L_{\varepsilon}[u_{h}]$ are almost piecewise constant in $D_{h}$, that is
away from any discontinuities, when $u$ is piecewise constant with $K$
inclusions. This indicates that the AS decomposition (3.3) is able to
approximate $u$ well throughout $\Omega$, which was proved in [1]:

######  Theorem 3.2.

Let $u$ be given by (3.5), $u_{h}\in V_{h}$ be its FE interpolant,
$\varphi_{k}$, $k\geqslant 1$ the eigenfunctions of the operator
$L_{\varepsilon}[u_{h}]$ and $\Pi_{K}[u_{h}]u$ be its projection. Then, for
every $\varepsilon>0$ and $h>0$ sufficiently small, there exists a constant
$C>0$, independent of $\varepsilon$ and $h$, but possibly depending on $u$,
such that

$\displaystyle\|u-\Pi_{K}[u_{h}]u\|_{L^{2}(\Omega)}\leqslant
C\sqrt{h+\varepsilon}.$ (3.9)

In summary, the AS decomposition $\Pi_{K}[u_{h}]u$ of $u$ approximates $u$
arbitrarily well, when $u$ is piecewise constant and consists of $K$
inclusions. Note that more general versions of Theorem 3.1 and Theorem 3.2
were proved in [2, Theorem 5] and [1, Theorem 3.6], which also apply to
piecewise constant functions with nonzero boundary values.

Finally, we remark that the operator $L_{\varepsilon}[u]$ corresponds to the
Fréchet derivative of the regularized TV-functional,
$\operatorname{TV}_{\varepsilon}(u)=\int_{\Omega}\sqrt{|\nabla
u|^{2}+\varepsilon^{2}}$, which reduces to the standard TV-functional
$\operatorname{TV}(u)=\int_{\Omega}|\nabla u|$ for $\varepsilon=0$ – see [17,
Remark 1]. In fact, for $v\in H_{0}^{1}$, $v\approx u_{h}$, the “energy-norm”
associated to the operator $L_{\varepsilon}[u_{h}]$ essentially corresponds to
the TV-“energy”,

$\displaystyle(L_{\varepsilon}[u_{h}]v,v)_{L^{2}(\Omega)}=\int_{\Omega}\mu_{\varepsilon}[u_{h}]|\nabla
v|^{2}=\int_{\Omega}\frac{|\nabla v|^{2}}{\sqrt{|\nabla
u_{h}|^{2}+\varepsilon^{2}}}\approx\int_{\Omega}|\nabla
v|=\operatorname{TV}(v).$ (3.10)

The AS decomposition also bears a remarkable resemblance to the spectral
decomposition of the nonlinear $\operatorname{TV}$-functional [6, 14].

##  4 Adaptive Spectral Inversion

We shall now combine the simple Adaptive Inversion iteration from Section 2
with the AS decomposition (3.3) from Section 3. Hence, at each iteration, we
shall determine the new search space $\Psi^{(m+1)}$ from the current minimizer
$u^{(m),\delta}\in\Psi^{(m)}$ using the eigenfunctions $\varphi_{k}$ of the
elliptic operator $L_{\varepsilon}[u^{(m),\delta}]$ in (3.2). Here, we recall
from Theorems 3.1 and 3.2 in Section 3 that a piecewise constant medium $u$
with $K$ inclusions is approximated with high accuracy by the first $K$
eigenfunctions of $L_{\varepsilon}[u_{h}]$ – see also [1, 2] for further
details.

First, we merge the current search space $\Psi^{(m)}$ with those
eigenfunctions and reduce its dimension, while ensuring that $u^{(m),\delta}$
is still well represented in the merged space. This first step corresponds to
the algorithm previously used in [2] – See Remark 2 below for further details.

Next, to incorporate the new angle condition (2.4) required for Theorem 2.1,
we shall include in $\Psi^{(m+1)}$ the most _sensitive_ eigenfunctions
$\varphi_{k}$ reordered according to their _sensitivities_ :

$\displaystyle|\sigma_{1}|\geqslant|\sigma_{2}|\geqslant\ldots\geqslant|\sigma_{k}|\geqslant\ldots\
,\qquad\sigma_{k}=D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\varphi_{k}.$
(4.1)

We shall now construct a subspace
$\Phi_{N_{\theta}}=\operatorname*{span}\\{\varphi_{1},\ldots,\varphi_{N_{\theta}}\\}\subset
L^{2}(\Omega)$, which contains at least one $d\in\Phi_{N_{\theta}}$ that
satisfies the angle condition (2.4), for a fixed $0<{\varepsilon_{\theta}}<1$.

######  Lemma 4.1.

Let $N_{\infty}$ be the largest index such that

$\displaystyle|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\varphi_{k}|\geqslant{\varepsilon_{\theta}}\|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\|,\qquad
k=1,\ldots,N_{\infty},$ (4.2)

holds. Then, $d=\sum_{k=1}^{N_{\infty}}\sigma_{k}\varphi_{k}$ satisfies the
angle condition (2.4).

###### Proof.

Since the eigenfunctions $\varphi_{k}$ are orthonormal, we have
$\|d\|_{L^{2}(\Omega)}=\|\sigma\|_{\ell^{2}}$ which yields

$\displaystyle|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})d|$
$\displaystyle=\left|\sum_{k=1}^{N_{\infty}}\sigma_{k}D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\varphi_{k}\right|=\sum_{k=1}^{N_{\infty}}|D\operatorname{\mathcal{J}}^{\delta}(u^{(m)})\varphi_{k}|^{2}$
(4.3)
$\displaystyle\geqslant{\varepsilon_{\theta}}\|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\|\sum_{k=1}^{N_{\infty}}|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\varphi_{k}|$
(4.4)
$\displaystyle={\varepsilon_{\theta}}\|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\|\|\sigma\|_{\ell^{1}}\geqslant{\varepsilon_{\theta}}\|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\|\|\sigma\|_{\ell^{2}}$
(4.5)
$\displaystyle={\varepsilon_{\theta}}\|d\|_{L^{2}(\Omega)}\|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\|.$
(4.6)

Thus, the angle condition (2.4) is satisfied. ∎

Note that for ${\varepsilon_{\theta}}<1$ there always exists
$N_{\infty}\geqslant 1$ sufficiently large such that (4.2) is satisfied.

######  Lemma 4.2.

Let $N_{2}$, be the smallest index such that

$\displaystyle\sqrt{\sum_{k=1}^{N_{2}}|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\varphi_{k}|^{2}}\geqslant{\varepsilon_{\theta}}\|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\|$
(4.7)

holds. Then, $d=\sum_{k=1}^{N_{2}}\sigma_{k}\varphi_{k}$ satisfies the angle
condition (2.4).

###### Proof.

Similarly, since the eigenfunctions $\varphi_{k}$ are orthonormal, we have

$\displaystyle|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})d|$
$\displaystyle=\sum_{k=1}^{N_{2}}|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\varphi_{k}|^{2}$
(4.8)
$\displaystyle=\|d\|_{L_{2}(\Omega)}\sqrt{\sum_{k=1}^{N_{2}}|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\varphi_{k}|^{2}}$
(4.9)
$\displaystyle\geqslant{\varepsilon_{\theta}}\|d\|_{L^{2}(\Omega)}\|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\|,$
(4.10)

because of (4.7). Hence, the angle condition (2.4) is satisfied. ∎

We refer to (4.2) as the $\ell^{\infty}$-criterion and to (4.7) as the
$\ell^{2}$-criterion. If we let

$\displaystyle N_{\theta}=\max\\{N_{\infty},N_{2}\\},$ (4.11)

we ensure that there exists at least one element in the subspace
$\Phi_{N_{\theta}}$ that satisfies the newly introduced angle condition (2.4).
Clearly, (4.2) is more stringent than (4.7) and thus $N_{2}\leqslant
N_{\infty}$. As we only compute in practice a finite number of eigenfunctions
$\varphi_{k}$, $1\leqslant k\leqslant N$, ordered according to their
eigenvalues, none might in fact satisfy (4.2). Then we set $N_{\infty}=0$ and
thereby $N_{\theta}=N_{2}$ in (4.11). In the unlikely case that no
$N_{2}\leqslant N$ satisfies (4.7), one needs to increase $N$ until (4.7)
holds to ensure the existence of $d\in\Phi_{N_{\theta}}$ which satisfies
(2.4). This yields the full Adaptive Spectral Inversion (ASI) Algorithm below.

Algorithm 2 Adaptive Spectral Inversion

0. initial guess $u^{(0),\delta}$, search space $\Psi^{(1)}$ set $m=1$, $\tau_{0}\geqslant 1$

0. reconstruction $u^{(m_{\ast}),\delta}$

1. while $\|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\|\neq 0$ do

2.  _Minimize_ $\operatorname{\mathcal{J}}^{\delta}$ in the current search space $\Psi^{(m)}$, $\dim\Psi^{(m)}=K_{m}$:  $\displaystyle u^{(m),\delta}=\operatorname*{argmin}_{u\in\Psi^{(m)}}\operatorname{\mathcal{J}}^{\delta}(u).$ (4.12)

3.  _Compute_ $\tau_{m}=\frac{1}{\delta}\|F(u^{(m),\delta})-y^{\delta}\|$.

4. if $\tau_{m}\leqslant\tau_{0}$ then

5. Discrepancy principle (2.20) is satisfied, return $u^{(m_{\ast}),\delta}=u^{(m-1),\delta}$.

6. Determine AS basis $\Phi^{(m+1)}$ from $L_{\varepsilon}[u^{(m),\delta}]$ via (3.2).

7. Merge AS basis and current search space: $\hat{\Psi}^{(m+1)}=\Phi^{(m+1)}\cup\Psi^{(m)}$.

8. Truncate $\hat{\Psi}^{(m+1)}$ while maintaining the accuracy of $u^{(m),\delta}$. This yields $\tilde{\Psi}^{(m+1)}$.

9. Add sensitivities. Include AS basis functions $\varphi_{k}\in\tilde{\Phi}^{(m+1)}\subset\Phi^{(m+1)}$ with maximal $|\sigma_{k}|$ in (4.1). This yields the new search space $\Psi^{(m+1)}=\tilde{\Psi}^{(m+1)}\cup\tilde{\Phi}^{(m+1)}$ with $\dim\Psi^{(m+1)}=K_{m+1}$.

10. $m\leftarrow m+1$

By combining the previous search space $\Psi^{(m)}$ with promising new search
directions determined from the current iterate $u^{(m),\delta}$, Steps 6 – 8
construct a low-dimensional subspace $\tilde{\Psi}^{(m+1)}$ able to represent
$u^{(m),\delta}$ up to a small error. In Step 9, we then add the basis
functions which will probably contribute most to the reduction of the misfit.
Below, we discuss in detail Steps 6 – 9 of the ASI Algorithm 2.

#### Step 6: Determine $\Phi^{(m+1)}$.

Once the minimizer $u^{(m),\delta}$ of $\operatorname{\mathcal{J}}^{\delta}$
in the current search space $\Psi^{(m)}$ of dimension $\dim\Psi^{(m)}=K_{m}$
has been found, we compute the first few eigenfunctions of the elliptic
operator $L_{\varepsilon}[u^{(m),\delta}]$, that is, we numerically solve the
(linear, symmetric and positive definite) eigenvalue problem (3.2), which
yields the new AS basis
$\Phi^{(m+1)}=\operatorname*{span}\left\\{\varphi_{1},\ldots,\varphi_{K_{m}}\right\\}$.
Here we arbitrarily set the dimension of $\Phi^{(m+1)}$ to that of
$\Psi^{(m)}$, though any other sufficiently large number of eigenfunctions
could be chosen.

#### Step 7: Merge.

Next we merge the previous search space $\Psi^{(m)}$ with $\Phi^{(m+1)}$ as
$\hat{\Psi}^{(m+1)}=\Psi^{(m)}\cup\Phi^{(m+1)}$ and compute an
$L^{2}$-orthonormal basis $\hat{\psi}_{k}$, $k=1,\ldots,\hat{K}_{m+1}$, of
$\hat{\Psi}^{(m+1)}$ via modified Gram-Schmidt. Since the dimension
$\hat{K}_{m+1}$ of the merged subspace $\hat{\Psi}^{(m+1)}$ may be as large as
$2K_{m}$, we now need to truncate $\hat{\Psi}^{(m+1)}$ to keep the number of
control variables small.

#### Step 8: Truncate.

To truncate $\hat{\Psi}^{(m+1)}$ and retain only those basis functions
$\hat{\psi}_{k}$ essential for representing $u^{(m),\delta}$, we proceed in
two steps. First, we compute an indicator $v$ close to $u^{(m),\delta}$, but
with minimal $\operatorname{TV}$-“energy”, to remove noise and preserve edges.
In doing so, we keep the computational costs low by minimizing the linearized
$\operatorname{TV}$-functional (3.10). Thus, the indicator $v$ satisfies

$\displaystyle\begin{split}\min_{v\in\hat{\Psi}^{(m+1)}}\quad&\int_{\Omega}\mu_{\varepsilon}[u^{(m),\delta}]|\nabla
v|^{2}\\\
\text{s.t.}\quad&\|v-u^{(m),\delta}\|_{L^{2}(\Omega)}\leqslant{\varepsilon_{\Psi}}\|u^{(m),\delta}\|_{L^{2}(\Omega)},\end{split}$
(4.13)

for a prescribed truncation tolerance ${\varepsilon_{\Psi}}>0$; typically,
${\varepsilon_{\Psi}}=5\%$. Since (4.13) is a quadratic optimization problem
with quadratic inequality constraints, computing $v$ is cheap.

Second, we reduce the dimension of $\hat{\Psi}^{(m+1)}$ by discarding all
basis functions that do not contribute much to the indicator $v$. Let
$\gamma_{k}$ denote the Fourier coefficients $\gamma_{k}$ of $v$,

$\displaystyle\gamma_{k}=(v,\hat{\psi}_{k})_{L^{2}(\Omega)},\qquad
k=1,\ldots,\hat{K}_{m+1},$ (4.14)

sorted in decreasing order and sort the $L^{2}$-orthonormal basis functions
$\hat{\psi}^{(m+1)}_{k}$ accordingly. Next, we determine the index

$\displaystyle N_{0}=\min\left\\{1\leqslant
K\leqslant\hat{K}_{m+1}:\;\sum_{k=K+1}^{\hat{K}_{m+1}}\gamma_{k}^{2}\leqslant{\varepsilon_{\Psi}}^{2}\|\gamma\|^{2}_{\ell^{2}}\right\\},$
(4.15)

such that the relative $L^{2}$-error in the Fourier expansion truncated at
$N_{0}$ is below ${\varepsilon_{\Psi}}$. To avoid drastic changes in the
dimension of the search space, we now calculate $\rho=N_{0}/K_{m}$. For
prescribed $\rho_{0},\rho_{1}$ such that $0<\rho_{0}\leqslant
1\leqslant\rho_{1}$, typically $\rho_{0}=0.8$ and $\rho_{1}=1.2$, we choose
the dimension $\tilde{K}_{m+1}$ of the truncated space $\tilde{\Psi}^{(m+1)}$
as follows:

  1. (i)

If $\rho\in[\rho_{0},\rho_{1}]$, set $\tilde{K}_{m+1}=N_{0}$.

  2. (ii)

If $\rho<\rho_{0}$, the dimension decreases too fast. Set
$\tilde{K}_{m+1}=\left\lceil\rho_{0}K_{m}\right\rceil$ and halve
${\varepsilon_{\Psi}}$ to avoid a rapid decrease in the dimension at the next
iteration.

  3. (iii)

If $\rho>\rho_{1}$, the dimension increases too fast. Set
$\tilde{K}_{m+1}=\left\lceil\rho_{1}K_{m}\right\rceil$ and double
${\varepsilon_{\Psi}}$ to avoid a rapid increase in the dimension at the next
iteration.

This yields the truncated space as
$\tilde{\Psi}^{(m+1)}=\operatorname*{span}\left\\{\hat{\psi}_{1},\ldots,\hat{\psi}_{\tilde{K}_{m+1}}\right\\}$.

#### Step 9: Add sensitivities.

Finally, we determine the most sensitive AS basis functions
$\tilde{\Phi}^{(m+1)}=\operatorname*{span}\left\\{\varphi_{1},\ldots,\varphi_{N_{\theta}}\right\\}$
from $\Phi^{(m+1)}$, reorderd according to their sensitivities as in (4.1). In
doing so, we use the $\ell^{\infty}$\- and $\ell^{2}$-criteria (4.2) and (4.7)
from Lemmas 4.1 and 4.2, respectively, to determine $N_{\theta}$. In the
unlikely event that the AS space $\Phi^{(m+1)}$ from Step 6 contains no
eigenfunction that satisfies (4.2) or (4.7), we can either increase the
dimension of $\Phi^{(m+1)}$, or set $\tilde{\Phi}^{(m+1)}=\emptyset$ (and
simply proceed). By combining $\tilde{\Psi}^{(m+1)}$ and
$\tilde{\Phi}^{(m+1)}$, we eventually obtain the subsequent search space

$\displaystyle\Psi^{(m+1)}=\tilde{\Psi}^{(m+1)}\cup\tilde{\Phi}^{(m+1)}=\operatorname*{span}\left\\{\psi_{1},\ldots,\psi_{K_{m+1}}\right\\}.$
(4.16)

Note that the sensitivity based selection procedure in Step 9 only ensures
that the angle condition (2.4) is satisfied by at least one element $d$ in
$\tilde{\Phi}^{(m+1)}\subset\Psi^{(m+1)}$, but not necessarily by the defect
$d^{(m),\delta}=u^{(m),\delta}-u^{(m-1),\delta}$ itself. The latter, stronger
condition would require verifying (2.4) at every iteration, and possibly
rejecting the new minimizer of (4.12) when not satisfied, while further
increasing the search space. Instead, Lemmas 4.1 and 4.2 guarantee that at
least one element in the search space satisfies the angle condition, thereby
making it rather likely that it will also be satisfied by the correction
$d^{(m),\delta}$.

######  Remark 2.

In the above Adaptive Spectral Inversion (ASI) Algorithm, Steps 1 – 8
correspond to the ASI Algorithm previously introduced in [2], yet with the
added growth control on $K_{m}$ in Step 8 (iii). Step 9, however, where
further eigenfunctions are included based on their sensitivites (4.1), is new
and will prove crucial for detecting small-scale features in the medium.
Henceforth we denote by $\text{ASI}_{0}$ the above ASI Algorithm without Step
9, similar to that from [2], to distinguish it from the present ASI Algorithm.

##  5 Numerical Results

To illustrate the accuracy and usefulness of the ASI Algorithm from Section 4,
we shall now apply it to two inverse medium problems of the form: Find $u\in
H_{1},y\in H_{2}$ that satisfy

$\displaystyle\min_{u,y}\quad$
$\displaystyle\frac{1}{2}\|y-y^{\delta}\|_{H_{2}}^{2}$ (5.1) s.t.
$\displaystyle A[u]y-f=0,$ (5.2)

for a given source $f$ and noisy data $y^{\delta}\in H_{2}$ as in (1.2). Each
inverse problem is governed by a distinct _forward problem_ (5.2) whose
solution, for any given medium333We refer to a function
$u:\mathbb{R}^{2}\to\mathbb{R}$ as a medium. $u$, is $y=A[u]^{-1}f$. Thus, we
can eliminate the constraint (5.2), which leads to the equivalent
unconstrained minimization problem: Find $u^{\dagger,\delta}\in H_{1}$ such
that

$\displaystyle u^{\dagger,\delta}=\operatorname*{argmin}_{u\in
H_{1}}\frac{1}{2}\|y[u]-y^{\delta}\|_{H_{2}}^{2}.$ (5.3)

For each inverse problem, we shall attempt to recover separately the two
different (unknown) media shown in Figure 4. The first consists of six discs,
each with a different value and radius. The second consists of three
inclusions: an open wedge with a sharp $90^{\circ}$ interior angle, a convex
drop-like inclusion with a sharp tip, and a kite-shaped inclusion, which is
non-convex with a smooth boundary. Since the unknown medium $u^{\dagger}$ must
be strictly bounded away from zero, we write $u^{\dagger}=1+u^{\dagger}_{0}$
for $u^{\dagger}_{0}$ with zero boundary. Then, we apply the ASI Algorithm 2
to minimize the misfit

$\displaystyle\min_{u}\operatorname{\mathcal{J}}^{\delta}(1+u)$ (5.4)

with $u$ equal to zero at the boundary.

Figure 4: Inverse Problem: The two media $u^{\dagger}$ used in the numerical
experiments. Left: six discs; right: three inclusions.

In all cases, we apply the ASI Algorithm 2 from Section 4 with the following
fixed parameter settings:

$\displaystyle\rho_{0}=0.8,\quad\rho_{1}=1.2,\quad{\varepsilon_{\theta}}=10^{-4},\quad\varepsilon_{\Psi}=0.05,\quad\tau_{0}=1.$
(5.5)

As initial guess, we always choose $u^{(0),\delta}=1$ constant and set the
initial search space $\Psi^{(1)}$ equal to the first $K_{1}=50$ or $K_{1}=100$
$L^{2}$-orthonormal eigenfunctions of the Laplacian sorted in non-decreasing
order w.r.t. their eigenvalues.

To assess the accuracy of the ASI method, we shall monitor the following
quantities: The dimension $K_{m}$ of the search space $\Psi^{(m)}$, the
relative error

$\displaystyle
e_{m}=\frac{\|u^{(m),\delta}-u^{\dagger}\|_{H_{1}}}{\|u^{\dagger}\|_{H_{1}}},$
(5.6)

the ratio $\tau_{m}$ from (2.20),

$\displaystyle\tau_{m}=\frac{\|y[u^{(m),\delta}]-y^{\delta}\|_{H_{2}}}{\delta},$
(5.7)

and the total number of iterations $m_{\ast}$, where $m=m_{\ast}+1$ is the
first index such that $\tau_{m}\leqslant\tau_{0}$; hence, the discrepancy
principle (2.20) is then satisfied with $\tau=\tau_{m_{\ast}}$.

###  5.1 Elliptic Inverse Problem

First, we consider as forward problem the elliptic differential equation

$\displaystyle\begin{aligned} -\nabla\cdot(u\nabla
y)&=f\qquad&&\text{in}\;\Omega,\\\
y&=0\qquad&&\text{on}\;\partial\Omega,\end{aligned}$ (5.8)

in the unit square $\Omega=(0,1)^{2}$ with constant right-hand side $f=100$;
hence, the observations $y^{\delta}$ are available throughout $\Omega$. In
[24, 25, 35], the one-dimensional version of (5.8) was considered and shown to
satisfy the Scherzer condition (2.27).

Here we shall compare the ASI Algorithm 2 from Section 4 with a standard grid-
based inversion method using Tikhonov regularization. We initialize the ASI
method with the constant initial guess $u^{(0),\delta}=1$ and let $\Psi^{(1)}$
equal the first $K_{1}=100$ Laplace eigenfunctions. The forward problem (5.8)
is solved in $H_{1}=H_{2}=V_{h}(\Omega)$, where $V_{h}$ corresponds to the
subspace of continuous, piecewise linear $\mathbb{P}^{1}$ finite elements (FE)
with mesh size $h>0$. Here, we use for both $u$ and $y$ the same fixed
triangular mesh with vertices located on a $400\times 400$ equidistant
Cartesian grid. Despite the large number ($160^{\prime}000$) of dof’s in the
FE representation, we recall that the ASI Algorithm determines the $m$-th
iterate only in the much smaller subspace $\Psi^{(m)}$ of dimension $K_{m}$.
In the $m$-th step of the ASI Algorithm the minimizer $u^{(m),\delta}$ to
(4.12) is determined using standard BFGS together with Armijo [36] instead of
Wolfe-Powell line search to keep the number of (expensive) gradient evaluation
small.

For the standard Tikhonov $L^{2}$-regularization approach, we minimize

$\displaystyle\min_{u\in
L^{2}(\Omega)}\frac{1}{2}\|y[u]-y^{\delta}\|_{L^{2}(\Omega)}^{2}+\frac{\alpha}{2}\|u\|_{L^{2}(\Omega)}^{2},$
(5.9)

where $\alpha>0$ denotes the regularization parameter. Again, we discretize
$u$ and $y$ with $\mathbb{P}^{1}$-FE on the same triangular mesh as for the
ASI method, yet due to the resulting large number of unknowns, we now solve
(5.9) using standard _limited memory BFGS_ (L-BFGS) [36] together with Armijo
line search. Following [12, Section $11.2$], we set the regularization
parameter $\alpha_{n}=2^{-n}$ at the $n$-th L-BFGS iteration.

To avoid any potential inverse crime, the exact data
$y^{\dagger}=y[u^{\dagger}]$ was computed from a $20\%$ finer mesh and
perturbed at each grid point $x_{j}$ as

$\displaystyle
y^{\delta}(x_{j})=y^{\dagger}(x_{j})+\hat{\delta}\,\eta_{j},\qquad
x_{j}\in\Omega,$ (5.10)

for a _noise level_ $\hat{\delta}\geqslant 0$. Here, $\eta_{j}$ corresponds to
Gaussian noise normalized such that the data misfit with
$y^{\dagger},y^{\delta}\in V_{h}(\Omega)$ satisfies exactly

$\displaystyle\|y^{\dagger}-y^{\delta}\|_{L^{2}(\Omega)}=\hat{\delta}=\delta.$
(5.11)

Table 1: Elliptic inverse problem, six discs: the relative error $e_{m_{\ast}}$, the total number of iterations $m_{\ast}$, and the dimension $K_{m_{\ast}}$ of the search space are shown for the ASI method and $L^{2}$-Tikhonov regularization. Method | ASI |  $L^{2}$-regularization
---|---|---
$\delta=\hat{\delta}$ | $1\%$ | $2\%$ | $5\%$ | $10\%$ | $1\%$ | $2\%$ | $5\%$ | $10\%$
$e_{m_{\ast}}$ | $2.7\%$ | $2.8\%$ | $3.9\%$ | $5.3\%$ | $5.3\%$ | $6.1\%$ | $8.0\%$ | $8.7\%$
$m_{\ast}$ | $50$ | $19$ | $15$ | $10$ | $5.3$ | $39$ | $15$ | $11$
$K_{m_{\ast}}$ | $67$ | $95$ | $76$ | $51$ | $\mathrm{nDof}=160^{\prime}000$
$\tau_{m_{\ast}}$ | $1.005$ | $1.0001$ | $1.000$ | $1.0003$ | $1.000$ | $1.0008$ | $1.005$ | $1.001$

In Table 1 we list the relative error $e_{m_{\ast}}$ at the final iteration
$m_{\ast}$ for the ASI method and $L^{2}$-Tikhonov regularization together
with the dimension $K_{m_{\ast}}$ of the search space $\Psi^{(m_{\ast})}$. As
expected, both methods require fewer iterations $m_{\ast}$ to achieve (2.20)
as the noise level $\delta$ increases, while the relative error obtained by
the ASI method is always smaller than that achieved by standard
$L^{2}$-Tikhonov regularization. Even with $\delta=10\%$ the ASI method yields
a smaller relative error than standard Tikhonov regularization using smaller
$\delta$. In Figure 5, we compare the reconstructed media
$u^{(m_{\ast}),\delta}$ for the two methods. Clearly, the reconstruction
obtained by the ASI method displays sharper contrasts and crisper edges while
the coefficients inside each disc are more accurate. Moreover, for
$L^{2}$-Tikhonov regularization, the background appears noisy and the
coefficients inside the inclusions are not accurately reconstructed.
Remarkably, the better reconstruction obtained by the ASI method is achieved
with as few as $K_{m_{\ast}}=100$ control variables only, in contrast to grid-
based Tikhonov regularization with more than $160^{\prime}000$ unknowns.

In Figure 6, we observe that the relative error $e_{m}$ and
$\|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\|$ decrease throughout
all iterations, as expected from Theorem 2.1 in Section 2, while $\tau_{m}$
tends to $1$. Note that the dimension $K_{m}$, or equivalently the number of
control variables, remains small for the ASI method regardless of $\delta$,
which keeps the computational effort low.

(a) ASI: $u^{(m_{\ast}),\delta}$ for $\hat{\delta}=1\%$

(b) ASI: $u^{(m_{\ast}),\delta}$ for $\hat{\delta}=2\%$

(c) ASI: $u^{(m_{\ast}),\delta}$ for $\hat{\delta}=5\%$

(d) Tikhonov: $u^{(m_{\ast}),\delta}$ for $\hat{\delta}=1\%$

(e) Tikhonov: $u^{(m_{\ast}),\delta}$ for $\hat{\delta}=2\%$

(f) Tikhonov: $u^{(m_{\ast}),\delta}$ for $\hat{\delta}=5\%$

Figure 5: Elliptic inverse problem, six discs: Reconstructed media using the
ASI method (top) or standard $L^{2}$-Tikhonov regularization (bottom) for
different noise levels $\hat{\delta}$.

(a) relative $L^{2}$ error $e_{m}$

(b) ratio $\tau_{m}$ from the discrepancy principle

(c) $\|D\operatorname{\mathcal{J}}(u^{(m),\delta})\|$

(d) dimension $K_{m}$

Figure 6: Elliptic inverse problem, six discs: The relative error (5.6), the
ratio $\tau_{m}$ from the discrepancy principle (5.7), the norm of the
gradient, and the dimension of the search space at iteration $m$ for different
noise levels $\delta$.

Next, we consider the (unknown) medium $u^{\dagger}$ shown in the right frame
of Figure 4, which consists of three distinct inclusions. Again, the ASI
Algorithm from Section 4 is able to recover the medium at the various noise
levels $\delta$, as shown in Figure 7: all three inclusions are clearly
visible with good contrast and sharp edges, except for the reentrant corner of
the open wedge with $5\%$ noise. $L^{2}$-Tikhonov regularization, however,
results in more noisy and blurred reconstructions, where the inclusions become
hardly visible beyond $5\%$ noise.

From Table 2, we again infer that the relative error $e_{m_{\ast}}$ of the ASI
method alway remains below that obtained with $L^{2}$-Tikhonov regularization:
Even with $\hat{\delta}=10\%$ noise, the ASI method remains more accurate than
$L^{2}$-Tikhonov regularization with as little as $1\%$ noise. Moreover, the
number of control variables $K_{m_{\ast}}$ used in the ASI method never
exceeds $160$, in comparison to approximately $160^{\prime}000$ control
variables used in the nodal FE representation for $L^{2}$-Tikhonov
regularization. As a consequence, the computational effort of the ASI method
always remains well below that with standard Tikhonov regularization.

In Figure 8, we observe that the relative error $e_{m}$ and
$\|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\|$ decrease with each
iteration, while the ratio $\tau_{m}$ tends to $1$ and the number of basis
functions $K_{m}$ of the search space $\Psi^{(m)}$, i.e. the number of control
variables, remains small, which again translates to a low computational
effort. Note that $\tau_{m}\approx 1$ implies that $u^{(m),\delta}$ (nearly)
yields an optimal data misfit because
$\|y^{\dagger}-y^{\delta}\|_{L_{2}(\Omega)}\leqslant\delta$.

Table 2: Elliptic inverse problem, three inclusions: the relative error $e_{m_{\ast}}$, the total number of iterations $m_{\ast}$, and the dimension $K_{m_{\ast}}$ of the search space are shown for the ASI method and $L^{2}$-Tikhonov regularization. Method | ASI |  $L^{2}$-regularization
---|---|---
$\delta=\hat{\delta}$ | $1\%$ | $2\%$ | $5\%$ | $10\%$ | $1\%$ | $2\%$ | $5\%$ | $10\%$
$e_{m_{\ast}}$ | $3.3\%$ | $3.7\%$ | $5.0\%$ | $6.4\%$ | $7\%$ | $8.3\%$ | $12.9\%$ | $13.0\%$
$m_{\ast}$ | $50$ | $25$ | $12$ | $8$ | $69$ | $48$ | $15$ | $15$
$K_{m_{\ast}}$ | $153$ | $120$ | $110$ | $72$ | $\mathrm{nDof}=160^{\prime}000$
$\tau_{m_{\ast}}$ | $1.001$ | $1.0001$ | $1.0001$ | $1.000$ | $1.0004$ | $1.0001$ | $1.029$ | $1.0014$

(a) ASI: $u^{(m_{\ast}),\delta}$ for $\hat{\delta}=1\%$

(b) ASI: $u^{(m_{\ast}),\delta}$ for $\hat{\delta}=2\%$

(c) ASI: $u^{(m_{\ast}),\delta}$ for $\hat{\delta}=5\%$

(d) Tikhonov: $u^{(m_{\ast}),\delta}$ for $\hat{\delta}=1\%$

(e) Tikhonov: $u^{(m_{\ast}),\delta}$ for $\hat{\delta}=2\%$

(f) Tikhonov: $u^{(m_{\ast}),\delta}$ for $\hat{\delta}=5\%$

Figure 7: Elliptic inverse problem, three inclusions: reconstructed medium
using the ASI method (top) or standard $L^{2}$-Tikhonov regularization
(bottom) for different noise levels $\hat{\delta}$.

(a) relative $L^{2}$ error $e_{m}$

(b) ratio $\tau_{m}$ from the discrepancy principle

(c) $\|D\operatorname{\mathcal{J}}(u^{(m),\delta})\|$

(d) dimension $K_{m}$

Figure 8: Elliptic inverse problem, three inclusions: The relative error
(5.6), the ratio $\tau_{m}$ from the discrepancy principle (5.7), the norm of
the gradient, and dimension of the search space, for every iteration $m$ for
different noise levels $\hat{\delta}>0$.

###  5.2 Time Dependent Inverse Scattering Problem

Next, we consider wave scattering from an unknown spatially distributed medium
illuminated by surrounding point sources. Hence, the forward problem (5.2) now
corresponds to the time-dependent wave equation in $\Omega=(0,1)^{2}$,

$\displaystyle\begin{aligned} \frac{\partial^{2}}{\partial
t^{2}}y_{\ell}(x,t)-\nabla\cdot\left(u(x)\nabla
y_{\ell}(x,t)\right)&=f_{\ell}(x,t),&&\quad x\in\Omega,\;t\in(0,T),\\\
y_{\ell}(x,0)=\tfrac{\partial}{\partial t}y_{\ell}(x,0)&=0,&&\quad
x\in\Omega,\\\ \tfrac{\partial}{\partial
t}y_{\ell}(x,t)+\sqrt{u(x)}\tfrac{\partial}{\partial
n}y_{\ell}(x,t)&=0,&&\quad x\in\partial\Omega,\;t\in(0,T),\end{aligned}$
(5.12)

with homogeneous initial conditions and first-order absorbing boundary
conditions. Here, $u(x)$ denotes the squared wave speed whereas the sources
$f_{\ell}(x,t)=g_{\ell}(x)r(t)$, $\ell=1,\ldots,N_{s}$, correspond to smoothed
Gaussian point sources

$\displaystyle g_{\ell}(x)=\kappa e^{\frac{(x-x_{\ell})^{2}}{s}},\qquad
x\in\Omega,$ (5.13)

in space, centered about distinct locations $x_{\ell}\in\Omega$, with
$s=10^{-2}$ and $\kappa=200$, and a Ricker wavelet [10, 30] in time,

$\displaystyle r(t)=(1-2\pi^{2}(\nu-t)^{2})e^{-\pi^{2}(\nu t-1)^{2}},\qquad
t\in[0,T],$ (5.14)

with central frequency $\nu=10$. In Figure 9, snapshots of the solution to the
forward problem (5.12) are shown at different times for the source located at
the top left corner.

(a) $t=0.4$

(b) $t=0.7$

(c) $t=1.0$

Figure 9: Snapshots of the solution to the wave equation (5.12).

The forward problem (5.12) is solved by a standard Galerkin-FE method, where
we again discretize $u(x)\in V_{h}$ using piecewise linear $\mathbb{P}^{1}$-FE
on a triangular mesh with vertices located on a $400\times 400$ equidistant
Cartesian grid. For $y(\cdot,t)\in\widetilde{V}_{h}$ in (5.12), however, we
use quadratic $\mathbb{P}^{2}$-FE with mass-lumping [7, 32] on a separate
triangular mesh with about $10$ elements per wavelength, resulting in
approximately $60^{\prime}000$ nodes. For the time integration, we use the
standard (fully explicit) leapfrog method with time-step $\Delta t\approx
4.5\cdot 10^{-4}$.

To generate the (synthetic) observations, we now place $N_{s}=32$ sources
located at $x_{\ell}$ equidistributed near the boundary and illuminate the
medium, one source at a time. In contrast to the previous example from Section
5.1, here the data $y^{\delta}_{\ell}$ is only available at the boundary
nodes, yet for all discrete time-steps $t_{n}=n\Delta t$,
$n=0,1,\ldots,N_{T}$, until the final time $T=2$ when the incident wave has
essentially left $\Omega$. Hence, we set $H_{1}=V_{h}(\Omega)$ and
$H_{2}=\widetilde{V}_{h}(\partial\Omega)\times\\{t_{n}\\}_{n=0}^{N_{T}}$ in
(5.3), where the misfit

$\displaystyle\operatorname{\mathcal{J}}^{\delta}(u)=\frac{1}{2}\sum_{\ell=1}^{N_{s}}\|y_{\ell}[u]-y_{\ell}^{\delta}\|_{L^{2}(\partial\Omega\times(0,T))}^{2}$
(5.15)

now accounts for the data $y_{\ell}^{\delta}$ from multiple sources
$\ell=1,\ldots,N_{s}$.

To avoid any potential inverse crime, the exact data
$y^{\dagger}_{\ell}=y_{\ell}[u^{\dagger}]$ was computed from a $20\%$ finer
mesh. The perturbed (noisy) data $y_{\ell}^{\delta}\in H_{2}$ was then
obtained at each boundary node $x_{j}$ and time-step $t_{n}$ as

$\displaystyle
y_{\ell}^{\delta}(x_{j},t_{n})=y_{\ell}^{\dagger}(x_{j},t_{n})\cdot(1+\hat{\delta}\cdot\left(\eta_{\ell})_{j,n}\right),\qquad(x_{j},t_{n})\in\partial\Omega\times(0,T),$
(5.16)

for a _noise level_ $\hat{\delta}\geqslant 0$. Here, $(\eta_{\ell})_{j,n}$
corresponds to normally distributed Gaussian noise with

$\displaystyle\sum_{\ell=1}^{N_{s}}\|y_{\ell}^{\dagger}-y_{\ell}^{\delta}\|_{L^{2}(\partial\Omega\times(0,T))}\leqslant\delta.$
(5.17)

To reduce computational cost during the inverse iteration, we do not minimize
(5.15) directly, but instead use a standard _sample average approximation_
(SAA) [19]: at each iteration, we combine all sources $f_{\ell}$ into a single
“super-shot” $f^{(m)}$,

$\displaystyle f^{(m)}=\sum_{\ell=1}^{N_{s}}\xi_{\ell}^{(m)}f_{\ell},$ (5.18)

where $\xi_{\ell}^{(m)}=\pm 1$ follow a Rademacher distribution with zero
mean. Thus, at each iteration $m$, we solve the minimization problem

$\displaystyle\min_{u\in\Psi^{(m)}}\frac{1}{2}\|y[u]-y^{(m),\delta}\|_{L^{2}(\partial\Omega\times(0,T))},$
(5.19)

with corresponding boundary observations

$\displaystyle
y^{(m),\delta}=\sum_{\ell=1}^{N_{s}}\xi_{\ell}^{(m)}y_{\ell}^{\delta}.$ (5.20)

First, we compare the ASI Algorithm from Section 4 with the former
$\text{ASI}_{0}$ Algorithm, where we omit Step 9 and thus ignore the most
sensitive AS basis functions required for the angle condition, see Remark 2.
To do so, we consider the (unknown) medium $u^{\dagger}$ consisting of six
discs shown in Figure 4. In Table 3, we observe that the ASI and
$\text{ASI}_{0}$ Algorithms perform similarly in terms of the relative error
$e_{m_{\ast}}$ and dimension of the search space $K_{m_{\ast}}$. However, when
we compare the reconstructed media $u^{(m_{\ast}),\delta}$ from Figure 10, we
observe that the $\text{ASI}_{0}$ Algorithm (bottom row) fails to reconstruct
the smallest disc with increasing noise, whereas the ASI Algorithm (top row)
always recovers even that smallest disc.

In Figure 11, the relative error $e_{m}$ decreases throughout all iterations,
while $K_{m}$, and hence the number of control variables, remains small, which
keeps the overall computational cost low. As expected from Theorem 2.1, the
norm of the gradient
$\|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\|$ decreases and the
ratio $\tau_{m}$ from the discrepancy principle (5.7) tends to 1.

Table 3: Inverse scattering problem, six discs: the relative error $e_{m_{\ast}}$, the total number of iterations $m_{\ast}$, and the dimension $K_{m_{\ast}}$ of the search space are shown for the ASI and $\text{ASI}_{0}$ method. Method | ASI | $\text{ASI}_{0}$
---|---|---
$\hat{\delta}$ | $1\%$ | $2\%$ | $5\%$ | $10\%$ | $1\%$ | $2\%$ | $5\%$ | $10\%$
$e_{m_{\ast}}$ | $2.6\%$ | $2.4\%$ | $2.0\%$ | $1.8\%$ | $2.4\%$ | $3.4\%$ | $4.4\%$ | $3.2\%$
$m_{\ast}$ | $50$ | $48$ | $50$ | $32$ | $50$ | $50$ | $31$ | $49$
$K_{m_{\ast}}$ | $81$ | $133$ | $85$ | $69$ | $83$ | $28$ | $57$ | $83$
$\tau_{m_{\ast}}$ | $1.076$ | $1.006$ | $1.0003$ | $1.0001$ | $1.047$ | $1.26$ | $1.061$ | $1.017$
$\delta$ | $0.18\%$ | $0.35\%$ | $0.88\%$ | $1.76\%$ | $0.18\%$ | $0.35\%$ | $0.88\%$ | $1.76\%$

(a) ASI: $u^{(m_{\ast}),\delta}$ for $\hat{\delta}=1\%$

(b) ASI: $u^{(m_{\ast}),\delta}$ for $\hat{\delta}=2\%$

(c) ASI: $u^{(m_{\ast}),\delta}$ for $\hat{\delta}=5\%$

(d) $\text{ASI}_{0}$: $u^{(m_{\ast}),\delta}$ for $\hat{\delta}=1\%$

(e) $\text{ASI}_{0}$: $u^{(m_{\ast}),\delta}$ for $\hat{\delta}=2\%$

(f) $\text{ASI}_{0}$: $u^{(m_{\ast}),\delta}$ for $\hat{\delta}=5\%$

Figure 10: Inverse scattering problem, six discs: Comparison of ASI (top) and
$\text{ASI}_{0}$ (bottom) for different noise levels $\hat{\delta}$.

(a) relative $L^{2}$ error $e_{m}$

(b) ratio $\tau_{m}$ from the discrepancy principle

(c) $\|D\operatorname{\mathcal{J}}(u^{(m),\delta})\|$

(d) dimension $K_{m}$

Figure 11: Inverse scattering problem, six discs: The relative error (5.6),
the ratio $\tau_{m}$ from the discrepancy principle (5.7), the norm of the
gradient, and the dimension of the search space are shown at every iteration
$m$ for the ASI method and the different noise levels $\hat{\delta}$.

Finally, we consider the medium $u^{\dagger}$ shown in Figure 4 with three
geometric inclusions. As shown in Figure 12, the ASI Algorithm recovers the
shape and height of all three inclusions with high fidelity and regardless of
the noise level $\hat{\delta}$. In Table 4 and Figure 13, we observe that the
relative error $e_{m}$ remains low for all $\hat{\delta}$ and decreases
throughout all iterations. Again, the ratio $\tau_{m}$ from (5.7) tends to $1$
while $\|D\operatorname{\mathcal{J}}^{\delta}(u^{(m),\delta})\|$ decreases.
For all noise levels $\hat{\delta}$, the number of basis functions $K_{m}$
never exceeds $140$, which greatly reduces the computational effort compared
to a standard nodal based optimization approach with $160^{\prime}000$ control
variables. Since the ASI and $\text{ASI}_{0}$ Algorithms performed similarly,
the results from the latter are omitted here.

Table 4: Inverse scattering problem, three inclusions: the relative error $e_{m_{\ast}}$, the total number of iterations $m_{\ast}$, and the dimension $K_{m_{\ast}}$ of the search space are shown for the ASI method. Method | ASI
---|---
$\hat{\delta}$ | $1\%$ | $2\%$ | $5\%$ | $10\%$
$e_{m_{\ast}}$ | $4.0\%$ | $3.4\%$ | $3.2\%$ | $3.4\%$
$m_{\ast}$ | $50$ | $49$ | $50$ | $50$
$K_{m_{\ast}}$ | $114$ | $83$ | $61$ | $45$
$\tau_{m_{\ast}}$ | $1.075$ | $1.015$ | $1.005$ | $1.002$
$\delta$ | $0.18\%$ | $0.35\%$ | $0.88\%$ | $1.46\%$

(a) ASI: $u^{(m_{\ast}),\delta}$ for $\hat{\delta}=1\%$

(b) ASI: $u^{(m_{\ast}),\delta}$ for $\hat{\delta}=2\%$

(c) ASI: $u^{(m_{\ast}),\delta}$ for $\hat{\delta}=5\%$

Figure 12: Inverse scattering problem, three inclusions: reconstructed medium
using the ASI method for different noise levels $\hat{\delta}$.

(a) relative $L^{2}$ error $e_{m}$

(b) ratio $\tau_{m}$ from the discrepancy principle

(c) $\|D\operatorname{\mathcal{J}}(u^{(m),\delta})\|$

(d) dimension $K_{m}$

Figure 13: Inverse scattering problem, three inclusions: The relative error
(5.6), the ratio $\tau_{m}$ from the discrepancy principle (5.7), the norm of
the gradient, and the dimension of the search space, at iteration $m$ for the
ASI method and different noise levels $\hat{\delta}=1\%,\;2\%\;5\%,\;10\%$.

##  6 Concluding Remarks

The Adaptive Spectral Inversion (ASI) method has proved remarkably effective
for the solution of PDE-constrained inverse medium problems. At the $m$-th
iteration, the ASI method minimizes the data misfit with added noise $\delta$
in a small subspace $\Psi^{(m)}$, which inherits key information from the
previous step while adding promising new search directions from the first few
eigenfunctions of the elliptic operator $L_{\varepsilon}[u^{(m),\delta}]$ in
(3.2). Since the operator $L_{\varepsilon}[u^{(m),\delta}]$ itself depends on
the current iterate, $u^{(m),\delta}$, so do its eigenfunctions and thus also
the new search space $\Psi^{(m+1)}$. The convergence of the ASI method hinges
upon the angle condition (2.4) – see Theorem 2.1 – which has been newly
included as a final step into each iteration. The full ASI Algorithm 2 is
listed in Section 4.

Under suitable assumptions, the ASI iteration stops after (finitely many)
$m_{*}$ iterations when the discrepancy principle (2.20) is satisfied. Hence,
the ASI Algorithm yields a genuine regularization method whose solution
$u^{(m_{*}),\delta}$ converges to the exact (noise-free) solution as
$\delta\rightarrow 0$, without the need for extra Tikhonov regularization –
see Theorem 2.3. By adapting the search space at each iteration while keeping
its dimension low, the ASI method achieves more accurate reconstructions than
standard grid-based Tikhonov $L^{2}$-regularization together with a
thousandfold decrease in the number of unknowns. Thanks to the newly
incorporated angle condition (2.4), the ASI Algorithm is able to detect even
the smallest inclusions in the medium, which previous versions of the
algorithm [18] at times failed to identify with increasing noise.

The added cost from the numerical solution of the eigenvalue problem (3.2) is
rather small. On the one hand, the Galerkin FE discretization leads to a
generalized eigenvalue problem which is sparse, symmetric, and positive
definite. On the other hand, we only require a few eigenfunctions which can be
efficiently computed by a standard Lanczos iteration. A further reduction in
the computational cost can easily be achieved by using an adaptively refined
FE mesh when solving (3.2), as in [17], which need not coincide with the mesh
used for discretizing the governing PDE. As the higher eigenfunctions become
increasingly localized, they can easily be “sparsified” simply by setting to
zero the smallest coefficients in their discrete FE representation [17].

Although we have concentrated here on two-dimensional inverse medium problems,
the AS decomposition in fact applies to arbitrary many space dimensions [1].
When the target medium is not piecewise constant but instead smoothly varying,
adaptive spectral bases resulting from different elliptic operators may be
more effective [18]. Clearly, the ASI approach would probably prove useful for
inverse source problems, for instance, or could be be combined with
alternative (globally convergent) inversion methods [5]. It could also be
applied to other inverse problems, unrelated to wave scattering, where the
state variable is governed by a different partial differential, or possibly
even integral, equation.

#### Acknowledgments.

We thank Daniel Baffet for useful comments and suggestions.

## References

  * [1] D. H. Baffet, Y. G. Gleichmann, and M. J. Grote, Error estimates for aptive spectral decompositions, Journal of Scientific Computing, 93 (2022).
  * [2] D. H. Baffet, M. J. Grote, and J. H. Tang, Adaptive spectral decompositions for inverse medium problems, Inverse Problems, 37 (2021), p. 025006.
  * [3] A. B. Bakushinsky and M. Y. Kokurin, Iterative Methods for Approximate Solution of Inverse Problems, Springer Netherlands, 2004\.
  * [4] H. H. Bauschke and P. L. Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces, Springer International Publishing, 2017.
  * [5] L. Beilina and M. V. Klibanov, Approximate Global Convergence and Adaptivity for Coefficient Inverse Problems, Springer US, 2012\.
  * [6] M. Burger, G. Gilboa, M. Moeller, L. Eckardt, and D. Cremers, Spectral decompositions using one-homogeneous functionals, SIAM Journal on Imaging Sciences, 9 (2016), pp. 1374–1408.
  * [7] G. Cohen, P. Joly, J. E. Roberts, and N. Tordjman, Higher order triangular finite elements with mass lumping for the wave equation, SIAM Journal on Numerical Analysis, 38 (2001), pp. 2047–2078.
  * [8] I. Daubechies, M. Defrise, and C. De Mol, An iterative thresholding algorithm for linear inverse problems with a sparsity constraint, Communications on Pure and Applied Mathematics, 57 (2004), pp. 1413–1457.
  * [9] M. de Buhan and M. Darbas, Numerical resolution of an electromagnetic inverse medium problem at fixed frequency, Computers & Mathematics with Applications, 74 (2017), pp. 3111–3128.
  * [10] 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 (2013), p. 085009.
  * [11] M. de Buhan and A. Osses, Logarithmic stability in determination of a 3d viscoelastic coefficient and a numerical example, Inverse Problems, 26 (2010), p. 095006.
  * [12] H. W. Engl, M. Hanke, and A. Neubauer, Regularization of Inverse Problems, vol. 375, Springer Science & Business Media, 1996.
  * [13] F. Faucher, O. Scherzer, and H. Barucq, Eigenvector models for solving the seismic inverse problem for the Helmholtz equation, Geophysical Journal International, (2020).
  * [14] G. Gilboa, M. Moeller, and M. Burger, Nonlinear spectral analysis via one-homogeneous functionals: Overview and future prospects, Journal of Mathematical Imaging and Vision, 56 (2016), pp. 300–319.
  * [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 (2019), p. 104001.
  * [16] C. W. Groetsch and A. Neubauer, Convergence of a general projection method for an operator equation of the first kind, Houston J. Math., 14 (1988), pp. 201–208.
  * [17] M. J. Grote, M. Kray, and U. Nahum, Adaptive eigenspace method for inverse scattering problems in the frequency domain, Inverse Problems, 33 (2017), p. 025006.
  * [18] M. J. Grote and U. Nahum, Adaptive eigenspace for multi-parameter inverse scattering problems, Computers & Mathematics with Applications, 77 (2019), pp. 3264–3280.
  * [19] E. Haber, M. Chung, and F. Herrmann, An effective method for parameter estimation with PDE constraints with multiple right-hand sides, SIAM Journal on Optimization, 22 (2012), pp. 739–757.
  * [20] U. Hämarik, E. Avi, and A. Ganina, On the solution of ill-posed problems by projection methods with a posteriori choice of the discretization level, Math. Model. Anal., 7 (2002), pp. 241–252.
  * [21] M. Hanke, A. Neubauer, and O. Scherzer, A convergence analysis of the Landweber iteration for nonlinear ill-posed problems, Numerische Mathematik, 72 (1995), pp. 21–37.
  * [22] F. J. Herrmann and G. Hennenfent, Non-parametric seismic data recovery with curvelet frames, Geophys. J. Int., (2008), pp. 233–248.
  * [23] B. Hofmann, P. Mathé, and S. V. Pereverzev, Regularization by projection: Approximation theoretic aspects and distance functions, Journal of Inverse and Ill-posed Problems, 15 (2007).
  * [24] B. Kaltenbacher, Regularization by projection with a posteriori discretization level choice for linear and nonlinear ill-posed problems, Inverse Problems, 16 (2000), pp. 1523–1539.
  * [25] B. Kaltenbacher, A. Neubauer, and O. Scherzer, Iterative Regularization Methods for Nonlinear Ill-Posed Problems, De Gruyter, 2008.
  * [26] B. Kaltenbacher and J. Offtermatt, A convergence analysis of regularization by discretization in preimage space, Math. Comput., 81 (2012), pp. 2049–2069.
  * [27] A. Kirsch, An Introduction to the Mathematical Theory of Inverse Problems, vol. 120 of Applied Mathematical Sciences, Springer-Verlag, New York, 1996.
  * [28] L. Landweber, An iteration formula for Fredholm integral equations of the first kind, American Journal of Mathematics, 73 (1951), pp. 615–624.
  * [29] Y. Lin, A. Abubakar, and T. M. Habashy, Seismic full-waveform inversion using truncated wavelet representations, (2012), pp. 1–6.  SEG Annual meeting 2012, Las Vegas.
  * [30] C. D. Lines and S. N. Chandler-Wilde, A time domain point source method for inverse scattering by rough surfaces, Computing, 75 (2005), pp. 157–180.
  * [31] I. Loris, H. Douma, G. Nolet, I. Daubechies, and C. Regone, Nonlinear regularization techniques for seismic tomography, Journal of Computational Physics, (2010), pp. 890–905.
  * [32] W. A. Mulder, Higher-order mass-lumped finite elements for the wave equation, Journal of Computational Acoustics, 09 (2001), pp. 671–680.
  * [33] F. Natterer, Regularisierung schlecht gestellter Probleme durch Projektionsverfahren, Numerische Mathematik, 28 (1977), pp. 329–341.
  * [34] C. Sanders, M. Bonnet, and W. Aquino, An adaptive eigenfunction basis strategy to reduce design dimension in topology optimization, Internat. J. Numer. Methods Engrg., 122 (2021), pp. 7452–7481.
  * [35] O. Scherzer, H. W. Engl, and K. Kunisch, Optimal a posteriori parameter choice for Tikhonov regularization for solving nonlinear ill-posed problems, SIAM Journal on Numerical Analysis, 30 (1993), pp. 1796–1838.
  * [36] S. Wright and J. Nocedal, Numerical Optimization, Springer New York, 2006.